1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23!=============================================================================== 24! Function: 25! --------- 26 27!> \file resssg2.f90 28!> 29!> \brief This subroutine prepares the solving of the coupled Reynolds stress 30!> components in \f$ R_{ij} - \varepsilon \f$ RANS (SSG) turbulence model. 31!> 32!> \remark 33!> - cvar_var(1,*) for \f$ R_{11} \f$ 34!> - cvar_var(2,*) for \f$ R_{22} \f$ 35!> - cvar_var(3,*) for \f$ R_{33} \f$ 36!> - cvar_var(4,*) for \f$ R_{12} \f$ 37!> - cvar_var(5,*) for \f$ R_{23} \f$ 38!> - cvar_var(6,*) for \f$ R_{13} \f$ 39!------------------------------------------------------------------------------- 40 41!------------------------------------------------------------------------------- 42! Arguments 43!______________________________________________________________________________. 44! mode name role 45!______________________________________________________________________________! 46!> \param[in] ivar variable number 47!> \param[in] gradv work array for the velocity grad term 48!> only for iturb=31 49!> \param[in] produc work array for production 50!> \param[in] gradro work array for grad rom 51!> (without rho volume) only for iturb=30 52!> \param[in] viscf visc*surface/dist at internal faces 53!> \param[in] viscb visc*surface/dist at edge faces 54!> \param[out] viscce Daly Harlow diffusion term 55!> \param[in] smbr working array 56!> \param[in] rovsdt working array 57!> \param[out] weighf working array 58!> \param[out] weighb working array 59!_______________________________________________________________________________ 60 61subroutine resssg2 & 62 ( ivar , & 63 gradv , produc , gradro , & 64 viscf , viscb , viscce , & 65 smbr , rovsdt , & 66 weighf , weighb ) 67 68!=============================================================================== 69 70!=============================================================================== 71! Module files 72!=============================================================================== 73 74use, intrinsic :: iso_c_binding 75 76use paramx 77use numvar 78use entsor 79use optcal 80use cstphy 81use cstnum 82use parall 83use period 84use lagran 85use mesh 86use field 87use field_operator 88use cs_f_interfaces 89use rotation 90use turbomachinery 91use cs_c_bindings 92 93!=============================================================================== 94 95implicit none 96 97! Arguments 98 99integer ivar 100 101double precision gradv(3, 3, ncelet) 102double precision produc(6, ncelet) 103double precision gradro(3, ncelet) 104double precision viscf(nfac), viscb(nfabor), viscce(6, ncelet) 105double precision smbr(6, ncelet) 106double precision rovsdt(6, 6, ncelet) 107double precision weighf(2, nfac), weighb(nfabor) 108 109! Local variables 110 111integer iel, isou, jsou 112integer ii , jj , kk , iii , jjj 113integer iwarnp 114integer imvisp 115integer st_prv_id 116integer iprev , inc, iccocg, ll 117integer t2v(3,3) 118integer iv2t(6), jv2t(6) 119integer key_t_ext_id, f_id, rot_id 120integer iroext 121 122double precision trprod, trrij 123double precision deltij(6), dij(3,3) 124double precision thets , thetv 125double precision aiksjk, aikrjk, aii ,aklskl, aikakj 126double precision xaniso(3,3), xstrai(3,3), xrotac(3,3), xprod(3,3), matrot(3,3) 127double precision sym_strain(6) 128double precision xrij(3,3), xnal(3), xnoral, xnnd 129double precision d1s2, d1s3, d2s3 130double precision alpha3 131double precision pij, phiij1, phiij2, epsij 132double precision phiijw, epsijw, epsijw_imp 133double precision ccorio 134double precision rctse 135double precision eigen_max 136double precision eigen_vals(3) 137double precision turb_schmidt 138double precision gradchk, gradro_impl, const 139double precision kseps 140double precision matrn(6), oo_matrn(6) 141double precision ceps_impl, cphi3impl, cphiw_impl 142double precision impl_drsm(6,6) 143double precision implmat2add(3,3) 144double precision impl_lin_cst, impl_id_cst 145double precision gkks3 146double precision grav(3) 147double precision, dimension(3,3) :: cvara_r 148 149character(len=80) :: label 150double precision, allocatable, dimension(:,:) :: grad 151double precision, allocatable, dimension(:) :: w1, w2 152double precision, allocatable, dimension(:,:), target :: buoyancy 153double precision, dimension(:), pointer :: crom, cromo 154double precision, dimension(:,:), pointer :: visten 155double precision, dimension(:), pointer :: cvara_ep, cvar_al 156double precision, dimension(:,:), pointer :: cvar_var, cvara_var 157double precision, dimension(:), pointer :: viscl, visct 158double precision, dimension(:,:), pointer :: c_st_prv 159double precision, dimension(:,:), pointer :: cpro_buoyancy 160 161type(var_cal_opt) :: vcopt 162 163!=============================================================================== 164 165!=============================================================================== 166! Initialization 167!=============================================================================== 168 169! Time extrapolation? 170call field_get_key_id("time_extrapolated", key_t_ext_id) 171 172call field_get_key_int(icrom, key_t_ext_id, iroext) 173 174! Allocate work arrays 175allocate(w1(ncelet), w2(ncelet)) 176 177! Generating the tensor to vector (t2v) and vector to tensor (v2t) mask arrays 178 179! a) t2v 180t2v(1,1) = 1; t2v(1,2) = 4; t2v(1,3) = 6; 181t2v(2,1) = 4; t2v(2,2) = 2; t2v(2,3) = 5; 182t2v(3,1) = 6; t2v(3,2) = 5; t2v(3,3) = 3; 183 184! b) i index of v2t 185iv2t(1) = 1; iv2t(2) = 2; iv2t(3) = 3; 186iv2t(4) = 1; iv2t(5) = 2; iv2t(6) = 1; 187 188! c) j index of v2t 189jv2t(1) = 1; jv2t(2) = 2; jv2t(3) = 3; 190jv2t(4) = 2; jv2t(5) = 3; jv2t(6) = 3; 191 192! d) kronecker symbol 193dij(:,:) = 0.0d0; 194dij(1,1) = 1.0d0; dij(2,2) = 1.0d0; dij(3,3) = 1.0d0; 195 196call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 197 198if (vcopt%iwarni.ge.1) then 199 call field_get_label(ivarfl(ivar), label) 200 write(nfecra,1000) label 201endif 202 203call field_get_val_s(icrom, crom) 204call field_get_val_s(iviscl, viscl) 205call field_get_val_s(ivisct, visct) 206 207call field_get_val_prev_s(ivarfl(iep), cvara_ep) 208if (iturb.ne.31) call field_get_val_s(ivarfl(ial), cvar_al) 209 210call field_get_val_v(ivarfl(ivar), cvar_var) 211call field_get_val_prev_v(ivarfl(ivar), cvara_var) 212 213d1s2 = 1.d0/2.d0 214d1s3 = 1.d0/3.d0 215d2s3 = 2.d0/3.d0 216 217do isou = 1, 3 218 deltij(isou) = 1.0d0 219enddo 220do isou = 4, 6 221 deltij(isou) = 0.0d0 222enddo 223 224! S as Source, V as Variable 225thets = thetst 226thetv = vcopt%thetav 227 228call field_get_key_int(ivarfl(ivar), kstprv, st_prv_id) 229if (st_prv_id .ge. 0) then 230 call field_get_val_v(st_prv_id, c_st_prv) 231else 232 c_st_prv=> null() 233endif 234 235if (st_prv_id.ge.0.and.iroext.gt.0) then 236 call field_get_val_prev_s(icrom, cromo) 237else 238 call field_get_val_s(icrom, cromo) 239endif 240 241! Coefficient of the "Coriolis-type" term 242if (icorio.eq.1) then 243 ! Relative velocity formulation 244 ccorio = 2.d0 245elseif (iturbo.eq.1) then 246 ! Mixed relative/absolute velocity formulation 247 ccorio = 1.d0 248else 249 ccorio = 0.d0 250endif 251 252!=============================================================================== 253! Production, Pressure-Strain correlation, dissipation, Coriolis 254!=============================================================================== 255 256! Source term 257! -rho*epsilon*( Cs1*aij + Cs2*(aikajk -1/3*aijaij*deltaij)) 258! -Cr1*P*aij + Cr2*rho*k*sij - Cr3*rho*k*sij*sqrt(aijaij) 259! +Cr4*rho*k(aik*sjk+ajk*sik-2/3*akl*skl*deltaij) 260! +Cr5*rho*k*(aik*rjk + ajk*rik) 261! -2/3*epsilon*deltaij 262 263! EBRSM 264if (iturb.eq.32) then 265 266 allocate(grad(3,ncelet)) 267 268 ! Compute the gradient of Alpha 269 iprev = 1 270 inc = 1 271 iccocg = 1 272 273 call field_gradient_scalar(ivarfl(ial), iprev, 0, inc, iccocg, grad) 274 275endif 276 277do iel = 1, ncel 278 279 ! EBRSM 280 if (iturb.eq.32) then 281 ! Compute the magnitude of the Alpha gradient 282 xnoral = ( grad(1,iel)*grad(1,iel) & 283 + grad(2,iel)*grad(2,iel) & 284 + grad(3,iel)*grad(3,iel) ) 285 xnoral = sqrt(xnoral) 286 ! Compute the unitary vector of Alpha 287 if (xnoral.le.epzero) then 288 xnal(1) = 0.d0 289 xnal(2) = 0.d0 290 xnal(3) = 0.d0 291 else 292 xnal(1) = grad(1,iel)/xnoral 293 xnal(2) = grad(2,iel)/xnoral 294 xnal(3) = grad(3,iel)/xnoral 295 endif 296 endif 297 298 ! Pij 299 xprod(1,1) = produc(1, iel) 300 xprod(1,2) = produc(4, iel) 301 xprod(1,3) = produc(6, iel) 302 xprod(2,2) = produc(2, iel) 303 xprod(2,3) = produc(5, iel) 304 xprod(3,3) = produc(3, iel) 305 306 ! Rotating frame of reference => "Coriolis production" term 307 308 rot_id = icorio 309 if (iturbo.eq.1) rot_id = irotce(iel) 310 311 if (rot_id .ge. 1) then 312 313 call coriolis_t(rot_id, 1.d0, matrot) 314 cvara_r(1,1) = cvara_var(1,iel) 315 cvara_r(2,2) = cvara_var(2,iel) 316 cvara_r(3,3) = cvara_var(3,iel) 317 cvara_r(1,2) = cvara_var(4,iel) 318 cvara_r(2,3) = cvara_var(5,iel) 319 cvara_r(1,3) = cvara_var(6,iel) 320 cvara_r(2,1) = cvara_var(4,iel) 321 cvara_r(3,2) = cvara_var(5,iel) 322 cvara_r(3,1) = cvara_var(6,iel) 323 324 do ii = 1, 3 325 do jj = ii, 3 326 do kk = 1, 3 327 xprod(ii,jj) = xprod(ii,jj) & 328 - ccorio*( matrot(ii,kk)*cvara_r(jj,kk) & 329 + matrot(jj,kk)*cvara_r(ii,kk) ) 330 enddo 331 enddo 332 enddo 333 endif 334 335 xprod(2,1) = xprod(1,2) 336 xprod(3,1) = xprod(1,3) 337 xprod(3,2) = xprod(2,3) 338 339 trprod = d1s2 * (xprod(1,1) + xprod(2,2) + xprod(3,3) ) 340 trrij = d1s2 * (cvara_var(1,iel) + cvara_var(2,iel) + cvara_var(3,iel)) 341 342 !-----> aII = aijaij 343 aii = 0.d0 344 aklskl = 0.d0 345 aiksjk = 0.d0 346 aikrjk = 0.d0 347 aikakj = 0.d0 348 349 ! aij 350 xaniso(1,1) = cvara_var(1,iel)/trrij - d2s3 351 xaniso(2,2) = cvara_var(2,iel)/trrij - d2s3 352 xaniso(3,3) = cvara_var(3,iel)/trrij - d2s3 353 xaniso(1,2) = cvara_var(4,iel)/trrij 354 xaniso(1,3) = cvara_var(6,iel)/trrij 355 xaniso(2,3) = cvara_var(5,iel)/trrij 356 xaniso(2,1) = xaniso(1,2) 357 xaniso(3,1) = xaniso(1,3) 358 xaniso(3,2) = xaniso(2,3) 359 360 ! Sij 361 xstrai(1,1) = gradv(1, 1, iel) 362 xstrai(1,2) = d1s2*(gradv(2, 1, iel)+gradv(1, 2, iel)) 363 xstrai(1,3) = d1s2*(gradv(3, 1, iel)+gradv(1, 3, iel)) 364 xstrai(2,1) = xstrai(1,2) 365 xstrai(2,2) = gradv(2, 2, iel) 366 xstrai(2,3) = d1s2*(gradv(3, 2, iel)+gradv(2, 3, iel)) 367 xstrai(3,1) = xstrai(1,3) 368 xstrai(3,2) = xstrai(2,3) 369 xstrai(3,3) = gradv(3, 3, iel) 370 371 ! omegaij 372 xrotac(1,1) = 0.d0 373 xrotac(1,2) = d1s2*(gradv(2, 1, iel)-gradv(1, 2, iel)) 374 xrotac(1,3) = d1s2*(gradv(3, 1, iel)-gradv(1, 3, iel)) 375 xrotac(2,1) = -xrotac(1,2) 376 xrotac(2,2) = 0.d0 377 xrotac(2,3) = d1s2*(gradv(3, 2, iel)-gradv(2, 3, iel)) 378 xrotac(3,1) = -xrotac(1,3) 379 xrotac(3,2) = -xrotac(2,3) 380 xrotac(3,3) = 0.d0 381 382 do ii=1,3 383 do jj = 1,3 384 ! aii = aij.aij 385 aii = aii+xaniso(ii,jj)*xaniso(ii,jj) 386 ! aklskl = aij.Sij 387 aklskl = aklskl + xaniso(ii,jj)*xstrai(ii,jj) 388 enddo 389 enddo 390 391 if (irijco .ne. 0) then 392 393 ! Initalize implicit matrices at 0 394 do isou = 1, 6 395 do jsou = 1, 6 396 impl_drsm(isou, jsou) = 0.0d0 397 end do 398 end do 399 do isou = 1, 3 400 do jsou = 1, 3 401 implmat2add(isou, jsou) = 0.0d0 402 end do 403 end do 404 405 impl_lin_cst = 0.0d0 406 impl_id_cst = 0.0d0 407 408 ! Computation of implicit components 409 ! ----------------------------------- 410 411 sym_strain(1) = xstrai(1,1) 412 sym_strain(2) = xstrai(2,2) 413 sym_strain(3) = xstrai(3,3) 414 sym_strain(4) = xstrai(1,2) 415 sym_strain(5) = xstrai(2,3) 416 sym_strain(6) = xstrai(1,3) 417 418 ! Global variables needed for SSG and EBRSM 419 ! ------------- 420 ! Computing the inverse matrix of R^n 421 ! Scaling by tr(R) in order to dodge inversion errors 422 do isou = 1, 6 423 matrn(isou) = cvara_var(isou,iel)/trrij 424 oo_matrn(isou) = 0.0d0 425 end do 426 427 ! Inversing the matrix 428 call symmetric_matrix_inverse(matrn, oo_matrn) 429 do isou = 1, 6 430 oo_matrn(isou) = oo_matrn(isou)/trrij 431 end do 432 433 ! Computing the maximal eigenvalue (in terms of norm!) of S 434 call calc_symtens_eigvals(sym_strain, eigen_vals) 435 eigen_max = maxval(abs(eigen_vals)) 436 437 ! Constant for the dissipation 438 ceps_impl = d1s3 * cvara_ep(iel) 439 440 if (iturb .eq. 31 ) then ! SSG - Epsilon 441 442 ! Identity constant for phi3 443 cphi3impl = abs(cssgr2 - cssgr3*sqrt(aii)) 444 445 ! Identity constant 446 impl_id_cst = - d2s3*cssgr1*min(trprod,0.0d0) & ! Phi1 447 - d1s3*cssgs2*cvara_ep(iel)*aii & ! Phi2 448 + cphi3impl * trrij * eigen_max & ! Phi3 449 + 2.0d0*d2s3*cssgr4*trrij*eigen_max & ! Phi4 450 + d2s3*trrij*cssgr4*max(aklskl,0.0d0) ! Phi4 451 452 ! Linear constant 453 impl_lin_cst = eigen_max * ( & 454 1.0d0 & ! Production 455 + cssgr4 & ! Phi 4 linear part 456 + cssgr5 ) ! Phi 5 linear part 457 458 do jsou = 1, 3 459 do isou = 1 ,3 460 iii = t2v(isou,jsou) 461 implmat2add(isou,jsou) = xrotac(isou,jsou) & 462 + impl_lin_cst*deltij(iii) & 463 + impl_id_cst*d1s2*oo_matrn(iii) & 464 + ceps_impl*oo_matrn(iii) 465 end do 466 end do 467 468 impl_drsm(:,:) = 0.0d0 469 call reduce_symprod33_to_66(implmat2add, impl_drsm) 470 471 else ! EBRSM 472 473 alpha3 = cvar_al(iel)**3 474 475 ! Phi3 constant 476 cphi3impl = abs(cebmr2 - cebmr3*sqrt(aii)) 477 478 ! PhiWall + epsilon_wall constants for EBRSM 479 cphiw_impl = 6.0d0*(1.0d0-alpha3)*cvara_ep(iel)/trrij 480 481 ! ------------- 482 ! The implicit components of Phi (pressure-velocity fluctuations) 483 ! are split into the linear part (A*R) and Id part (A*Id). 484 ! ------------- 485 486 ! Identity constant 487 impl_id_cst = alpha3 * ( & 488 - d2s3*cebmr1*min(trprod,0.0d0) & ! Phi1 489 + cphi3impl * trrij * eigen_max & ! Phi3 490 + 2.0d0*d2s3*cebmr4*trrij*eigen_max & ! Phi4 491 + d2s3*trrij*cebmr4*max(aklskl,0.d0) ) ! Phi4 492 493 ! Linear constant 494 impl_lin_cst = eigen_max * ( & 495 1.0d0 & ! Production 496 + cebmr4 * alpha3 & ! Phi4 Linear part 497 + cebmr5 * alpha3 ) & ! Phi5 Linear part 498 + cphiw_impl ! Epsilon + Phi wall 499 500 do jsou = 1, 3 501 do isou = 1, 3 502 iii = t2v(isou,jsou) 503 implmat2add(isou,jsou) = xrotac(isou,jsou) & 504 + impl_lin_cst*deltij(iii) & 505 + impl_id_cst*d1s2*oo_matrn(iii) & 506 + alpha3*ceps_impl*oo_matrn(iii) 507 end do 508 end do 509 510 impl_drsm(:,:) = 0.0d0 511 call reduce_symprod33_to_66(implmat2add, impl_drsm) 512 513 end if ! EBRSM 514 515 endif ! (irijco .ne. 0) 516 517 ! Rotating frame of reference => "absolute" vorticity 518 if (icorio.eq.1) then 519 do ii = 1, 3 520 do jj = 1, 3 521 xrotac(ii,jj) = xrotac(ii,jj) + matrot(ii,jj) 522 enddo 523 enddo 524 endif 525 526 do isou = 1, 6 527 iii = iv2t(isou) 528 jjj = jv2t(isou) 529 aiksjk = 0 530 aikrjk = 0 531 aikakj = 0 532 do kk = 1,3 533 ! aiksjk = aik.Sjk+ajk.Sik 534 aiksjk = aiksjk + xaniso(iii,kk)*xstrai(jjj,kk) & 535 + xaniso(jjj,kk)*xstrai(iii,kk) 536 ! aikrjk = aik.Omega_jk + ajk.omega_ik 537 aikrjk = aikrjk + xaniso(iii,kk)*xrotac(jjj,kk) & 538 + xaniso(jjj,kk)*xrotac(iii,kk) 539 ! aikakj = aik*akj 540 aikakj = aikakj + xaniso(iii,kk)*xaniso(kk,jjj) 541 enddo 542 543 ! If we extrapolate the source terms (rarely), 544 ! we put all in the previous ST. 545 ! We do not implicit the term with Cs1*aij neither the term with Cr1*P*aij. 546 ! Otherwise, we put all in smbr and we can implicit Cs1*aij 547 ! and Cr1*P*aij. Here we store the second member and the implicit term 548 ! in W1 and W2, to avoid the test(ST_PRV_ID.GE.0) 549 ! in the ncel loop 550 ! In the term with W1, which is dedicated to be extrapolated, we use 551 ! cromo. 552 ! The implicitation of the two terms can also be done in the case of 553 ! extrapolation, by isolating those two terms and by putting it in 554 ! the RHS but not in the prev. ST and by using ipcrom .... to be modified 555 ! if needed. 556 557 if (iturb.eq.31) then 558 559 ! Explicit terms 560 pij = xprod(iii,jjj) 561 phiij1 = -cvara_ep(iel)* & 562 (cssgs1*xaniso(iii,jjj)+cssgs2*(aikakj-d1s3*deltij(isou)*aii)) 563 phiij2 = - cssgr1*trprod*xaniso(iii,jjj) & 564 + trrij*xstrai(iii,jjj)*(cssgr2-cssgr3*sqrt(aii)) & 565 + cssgr4*trrij*(aiksjk-d2s3*deltij(isou)*aklskl) & 566 + cssgr5*trrij* aikrjk 567 epsij = -d2s3*cvara_ep(iel)*deltij(isou) 568 569 w1(iel) = cromo(iel)*cell_f_vol(iel)*(pij+phiij1+phiij2+epsij) 570 571 ! Implicit terms 572 w2(iel) = cell_f_vol(iel)/trrij*crom(iel) & 573 * (cssgs1*cvara_ep(iel) + cssgr1*max(trprod,0.d0)) 574 575 ! EBRSM 576 else 577 578 xrij(1,1) = cvara_var(1,iel) 579 xrij(2,2) = cvara_var(2,iel) 580 xrij(3,3) = cvara_var(3,iel) 581 xrij(1,2) = cvara_var(4,iel) 582 xrij(2,3) = cvara_var(5,iel) 583 xrij(1,3) = cvara_var(6,iel) 584 xrij(2,1) = xrij(1,2) 585 xrij(3,1) = xrij(1,3) 586 xrij(3,2) = xrij(2,3) 587 588 ! Compute the explicit term 589 590 ! Calculation of the terms near the walls and et almost homogeneous 591 ! of phi and epsilon 592 593 ! Calculation of the term near the wall \f$ \Phi_{ij}^w \f$ --> W3 594 ! 595 ! Phiw = -5.0 * (eps/k) * [ R*Xn + Xn^T*R - 0.5*tr(Xn*R)*(Xn + Id) ] 596 ! 597 phiijw = 0.d0 598 xnnd = d1s2*( xnal(iii)*xnal(jjj) + deltij(isou) ) 599 do kk = 1, 3 600 phiijw = phiijw + xrij(iii,kk)*xnal(jjj)*xnal(kk) 601 phiijw = phiijw + xrij(jjj,kk)*xnal(iii)*xnal(kk) 602 do ll = 1, 3 603 phiijw = phiijw - xrij(kk,ll)*xnal(kk)*xnal(ll)*xnnd 604 enddo 605 enddo 606 phiijw = -5.d0*cvara_ep(iel)/trrij * phiijw 607 608 ! Calculation of the almost homogeneous term \f$ \phi_{ij}^h \f$ --> W4 609 phiij1 = -cvara_ep(iel)*cebms1*xaniso(iii,jjj) 610 phiij2 = -cebmr1*trprod*xaniso(iii,jjj) & 611 +trrij*xstrai(iii,jjj)*(cebmr2-cebmr3*sqrt(aii)) & 612 +cebmr4*trrij *(aiksjk-d2s3*deltij(isou)*aklskl) & 613 +cebmr5*trrij * aikrjk 614 615 ! Calculation of \f $\e_{ij}^w \f$ --> W5 (Rotta model) 616 ! Rij/k*epsilon 617 epsijw = xrij(iii,jjj)/trrij *cvara_ep(iel) 618 619 ! Calcul de \e_{ij}^h --> W6 620 epsij = d2s3*cvara_ep(iel)*deltij(isou) 621 622 ! Calcul du terme source explicite de l'equation des Rij 623 ! \f[ P_{ij} + (1-\alpha^3)\Phi_{ij}^w + \alpha^3\Phi_{ij}^h 624 ! - (1-\alpha^3)\e_{ij}^w - \alpha^3\e_{ij}^h ]\f$ --> W1 625 alpha3 = cvar_al(iel)**3 626 627 w1(iel) = cell_f_vol(iel)*crom(iel)*( & 628 xprod(iii,jjj) & 629 + (1.d0-alpha3)*phiijw + alpha3*(phiij1+phiij2) & 630 - (1.d0-alpha3)*epsijw - alpha3*epsij) 631 632 ! Implicit term 633 634 if (irijco.eq.0) then 635 ! Implicitation of epsijw 636 ! (the factor 5 appears when we calculate \f$ Phi_{ij}^w - epsijw\f$) 637 ! epsijw_imp = 5.d0 * (1.d0-alpha3)*cvara_ep(iel)/trrij & 638 ! + (1.d0-alpha3)*cvara_ep(iel)/trrij 639 epsijw_imp = 6.d0 * (1.d0-alpha3)*cvara_ep(iel)/trrij 640 else 641 ! FIXME 642 epsijw_imp = 0.d0 643 endif 644 645 ! The term below corresponds to the implicit part of SSG 646 ! in the context of elliptical weighting, it is multiplied by 647 ! \f$ \alpha^3 \f$ 648 w2(iel) = cell_f_vol(iel)*crom(iel) & 649 * ( cebms1*cvara_ep(iel)/trrij*alpha3 & 650 + cebmr1*max(trprod/trrij,0.d0)*alpha3 & 651 + epsijw_imp) 652 653 endif ! EBRSM 654 655 if (st_prv_id.ge.0) then 656 c_st_prv(isou,iel) = c_st_prv(isou,iel) + w1(iel) 657 else 658 smbr(isou,iel) = smbr(isou,iel) + w1(iel) 659 rovsdt(isou,isou,iel) = rovsdt(isou,isou,iel) + w2(iel) 660 661 if (irijco.ne.0) then 662 ! Careful ! Inversion of the order of the coefficients since 663 ! rovsdt matrix is then used by a c function for the linear solving 664 do jsou = 1, 6 665 rovsdt(jsou,isou,iel) = rovsdt(jsou,isou,iel) + cell_f_vol(iel) & 666 *crom(iel) * impl_drsm(isou,jsou) 667 end do 668 endif 669 endif 670 enddo 671enddo 672 673if (iturb.eq.32) then 674 deallocate(grad) 675endif 676 677!=============================================================================== 678! Buoyancy source term 679!=============================================================================== 680 681if (igrari.eq.1) then 682 683 call field_get_id_try("rij_buoyancy", f_id) 684 if (f_id.ge.0) then 685 call field_get_val_v(f_id, cpro_buoyancy) 686 else 687 ! Allocate a work array 688 allocate(buoyancy(6,ncelet)) 689 cpro_buoyancy => buoyancy 690 endif 691 692 call rijthe2(gradro, cpro_buoyancy) 693 694 ! If we extrapolate the source terms: previous ST 695 if (st_prv_id.ge.0) then 696 do iel = 1, ncel 697 do isou = 1, 6 698 c_st_prv(isou,iel) = c_st_prv(isou,iel) + cpro_buoyancy(isou,iel) * cell_f_vol(iel) 699 enddo 700 enddo 701 ! Otherwise smbr 702 else 703 do iel = 1, ncel 704 do isou = 1, 6 705 smbr(isou,iel) = smbr(isou,iel) + cpro_buoyancy(isou,iel) * cell_f_vol(iel) 706 enddo 707 enddo 708 endif 709 710 ! Free memory 711 if (allocated(buoyancy)) deallocate(buoyancy) 712 713 if (irijco.ne.0) then 714 715 grav(1) = gx 716 grav(2) = gy 717 grav(3) = gz 718 719 ! Implicit buoyancy term 720 if (iscalt.gt.0) then 721 call field_get_key_double(ivarfl(isca(iscalt)), ksigmas, turb_schmidt) 722 const = -1.5d0 * cmu / turb_schmidt 723 else 724 const = -1.5d0 * cmu 725 end if 726 727 do iel = 1, ncel 728 729 do jsou = 1, 6 730 do isou = 1, 6 731 impl_drsm(isou,jsou) = 0.0d0 732 end do 733 end do 734 implmat2add(:,:) = 0.0d0 735 736 kseps = (cvara_var(1,iel) + cvara_var(2,iel) + cvara_var(3,iel)) & 737 / (2.0d0*cvara_ep(iel)) 738 739 xrij(1,1) = cvara_var(1,iel) 740 xrij(2,2) = cvara_var(2,iel) 741 xrij(3,3) = cvara_var(3,iel) 742 xrij(1,2) = cvara_var(4,iel) 743 xrij(2,3) = cvara_var(5,iel) 744 xrij(1,3) = cvara_var(6,iel) 745 xrij(2,1) = xrij(1,2) 746 xrij(3,1) = xrij(1,3) 747 xrij(3,2) = xrij(2,3) 748 749 trrij = 0.5d0 * (xrij(1,1) + xrij(2,2) + xrij(3,3)) 750 751 gkks3 = 0.d0 752 do jsou = 1, 3 753 do isou = 1, 3 754 gkks3 = gkks3 + grav(isou) * gradro(jsou,iel) * xrij(isou, jsou) 755 enddo 756 enddo 757 gkks3 = const * kseps * gkks3 * crij3 * 2.d0 / 3.d0 758 759 if (gkks3.le.0.d0) then 760 ! Term "C3 tr(G) Id" 761 ! Computing the inverse matrix of R^n 762 ! Scaling by tr(R) in order to dodge inversion errors 763 do isou = 1, 6 764 matrn(isou) = cvara_var(isou,iel)/trrij 765 oo_matrn(isou) = 0.0d0 766 end do 767 768 ! Inversing the matrix 769 call symmetric_matrix_inverse(matrn, oo_matrn) 770 do isou = 1, 6 771 oo_matrn(isou) = oo_matrn(isou)/trrij 772 end do 773 774 do jsou = 1, 3 775 do isou = 1, 3 776 iii = t2v(isou,jsou) 777 implmat2add(isou,jsou) = - 0.5d0 * gkks3 * oo_matrn(iii) 778 end do 779 end do 780 781 endif 782 783 gradchk = grav(1)*gradro(1,iel) + grav(2)*gradro(2,iel) + grav(3)*gradro(3,iel) 784 if (gradchk .gt. 0.0d0) then 785 786 ! Implicit term written as: 787 ! Po . R^n+1 + R^n+1 . Po^t 788 ! with Po proportional to "g (x) Grad rho" 789 gradro_impl = const * (1.0d0-crij3) * kseps 790 do jsou = 1, 3 791 do isou = 1, 3 792 implmat2add(isou,jsou) = implmat2add(isou,jsou) & 793 - gradro_impl * grav(isou) * gradro(jsou,iel) 794 end do 795 end do 796 end if 797 798 ! Compute the 6x6 matrix A which verifies 799 ! A . R = M . R + R . M^t 800 call reduce_symprod33_to_66(implmat2add, impl_drsm) 801 802 do isou = 1, 6 803 do jsou = 1, 6 804 ! Careful ! Inversion of the order of the coefficients since 805 ! rovsdt matrix is then used by a c function for the linear solving 806 rovsdt(jsou,isou,iel) = rovsdt(jsou,isou,iel) & 807 + cell_f_vol(iel) * impl_drsm(isou,jsou) 808 end do 809 end do 810 811 end do 812 813 endif ! (irijco .ne. 0) 814 815endif 816 817!=============================================================================== 818! Diffusion term (Daly Harlow: generalized gradient hypothesis method) 819!=============================================================================== 820 821! Symmetric tensor diffusivity (GGDH) 822if (iand(vcopt%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 823 824 call field_get_val_v(ivsten, visten) 825 826 do iel = 1, ncel 827 viscce(1,iel) = visten(1,iel) + viscl(iel) 828 viscce(2,iel) = visten(2,iel) + viscl(iel) 829 viscce(3,iel) = visten(3,iel) + viscl(iel) 830 viscce(4,iel) = visten(4,iel) 831 viscce(5,iel) = visten(5,iel) 832 viscce(6,iel) = visten(6,iel) 833 enddo 834 835 iwarnp = vcopt%iwarni 836 837 call vitens(viscce, iwarnp, weighf, weighb, viscf, viscb) 838 839! Scalar diffusivity 840else 841 842 do iel = 1, ncel 843 rctse = csrij * visct(iel) / cmu 844 w1(iel) = viscl(iel) + vcopt%idifft*rctse 845 enddo 846 847 imvisp = vcopt%imvisf 848 849 call viscfa(imvisp, w1, viscf, viscb) 850 851endif 852 853! Free memory 854 855deallocate(w1, w2) 856 857!-------- 858! Formats 859!-------- 860 861 1000 format(/,' Solving variable ', a8 ,/) 862 863!---- 864! End 865!---- 866 867return 868 869end subroutine 870