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 resrij2.f90 28!> 29!> \brief This subroutine prepares the solving of the coupled Reynolds stress 30!> components in \f$ R_{ij} - \varepsilon \f$ RANS (LRR) 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 resrij2 & 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 cs_f_interfaces 88use rotation 89use turbomachinery 90use cs_c_bindings 91 92!=============================================================================== 93 94implicit none 95 96! Arguments 97 98integer ivar 99 100double precision gradv(3, 3, ncelet) 101double precision produc(6, ncelet) 102double precision gradro(3, ncelet) 103double precision viscf(nfac), viscb(nfabor), viscce(6, ncelet) 104double precision smbr(6, ncelet) 105double precision rovsdt(6, 6, ncelet) 106double precision weighf(2, nfac), weighb(nfabor) 107 108! Local variables 109 110integer iel 111integer isou, jsou 112integer ii , jj , kk 113integer iflmas, iflmab 114integer iwarnp 115integer imvisp 116integer st_prv_id 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 matrot(3,3) 126double precision d1s2, d1s3, d2s3 127double precision ccorio 128double precision rctse 129double precision, dimension(3,3) :: cvara_r 130 131character(len=80) :: label 132double precision, allocatable, dimension(:,:), target :: buoyancy 133double precision, allocatable, dimension(:) :: w1 134double precision, allocatable, dimension(:,:) :: w2 135double precision, dimension(:), pointer :: imasfl, bmasfl 136double precision, dimension(:), pointer :: crom, cromo 137double precision, dimension(:,:), pointer :: visten 138double precision, dimension(:), pointer :: cvara_ep 139double precision, dimension(:,:), pointer :: cvar_var, cvara_var 140double precision, dimension(:), pointer :: viscl 141double precision, dimension(:,:), pointer :: c_st_prv 142double precision, dimension(:,:), pointer :: cpro_buoyancy 143 144integer iii, jjj 145 146double precision impl_drsm(6,6) 147double precision implmat2add(3,3) 148double precision impl_lin_cst, impl_id_cst 149double precision aiksjk, aikrjk, aii ,aklskl, aikakj 150double precision xaniso(3,3), xstrai(3,3), xrotac(3,3), xprod(3,3) 151double precision sym_strain(6) 152double precision matrn(6), oo_matrn(6) 153double precision eigen_vals(3) 154double precision ceps_impl 155double precision eigen_max 156double precision pij, phiij1, phiij2, epsij 157 158type(var_cal_opt) :: vcopt 159 160!=============================================================================== 161 162!=============================================================================== 163! Initialization 164!=============================================================================== 165 166! Time extrapolation? 167call field_get_key_id("time_extrapolated", key_t_ext_id) 168 169call field_get_key_int(icrom, key_t_ext_id, iroext) 170 171! Allocate work arrays 172allocate(w1(ncelet)) 173allocate(w2(6,ncelet)) 174 175! Generating the tensor to vector (t2v) and vector to tensor (v2t) mask arrays 176 177! a) t2v 178t2v(1,1) = 1; t2v(1,2) = 4; t2v(1,3) = 6; 179t2v(2,1) = 4; t2v(2,2) = 2; t2v(2,3) = 5; 180t2v(3,1) = 6; t2v(3,2) = 5; t2v(3,3) = 3; 181 182! b) i index of v2t 183iv2t(1) = 1; iv2t(2) = 2; iv2t(3) = 3; 184iv2t(4) = 1; iv2t(5) = 2; iv2t(6) = 1; 185 186! c) j index of v2t 187jv2t(1) = 1; jv2t(2) = 2; jv2t(3) = 3; 188jv2t(4) = 2; jv2t(5) = 3; jv2t(6) = 3; 189 190! d) kronecker symbol 191dij(:,:) = 0.0d0; 192dij(1,1) = 1.0d0; dij(2,2) = 1.0d0; dij(3,3) = 1.0d0; 193 194call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 195 196if (vcopt%iwarni.ge.1) then 197 call field_get_label(ivarfl(ivar), label) 198 write(nfecra,1000) label 199endif 200 201call field_get_val_s(icrom, crom) 202call field_get_val_s(iviscl, viscl) 203 204call field_get_val_prev_s(ivarfl(iep), cvara_ep) 205 206call field_get_val_v(ivarfl(ivar), cvar_var) 207call field_get_val_prev_v(ivarfl(ivar), cvara_var) 208 209call field_get_key_int(ivarfl(iu), kimasf, iflmas) 210call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 211call field_get_val_s(iflmas, imasfl) 212call field_get_val_s(iflmab, bmasfl) 213 214d1s2 = 1.d0/2.d0 215d1s3 = 1.d0/3.d0 216d2s3 = 2.d0/3.d0 217 218do isou = 1, 3 219 deltij(isou) = 1.0d0 220enddo 221do isou = 4, 6 222 deltij(isou) = 0.0d0 223enddo 224 225! S as Source, V as Variable 226thets = thetst 227thetv = vcopt%thetav 228 229call field_get_key_int(ivarfl(ivar), kstprv, st_prv_id) 230if (st_prv_id .ge. 0) then 231 call field_get_val_v(st_prv_id, c_st_prv) 232else 233 c_st_prv=> null() 234endif 235 236if (st_prv_id.ge.0.and.iroext.gt.0) then 237 call field_get_val_prev_s(icrom, cromo) 238else 239 call field_get_val_s(icrom, cromo) 240endif 241 242! Coefficient of the "Coriolis-type" term 243if (icorio.eq.1) then 244 ! Relative velocity formulation 245 ccorio = 2.d0 246elseif (iturbo.eq.1) then 247 ! Mixed relative/absolute velocity formulation 248 ccorio = 1.d0 249else 250 ccorio = 0.d0 251endif 252 253!=============================================================================== 254! Production, Pressure-Strain correlation, dissipation 255!=============================================================================== 256 257do iel = 1, ncel 258 259 ! Initalize implicit matrices at 0 260 do isou = 1, 6 261 do jsou = 1, 6 262 impl_drsm(isou, jsou) = 0.0d0 263 end do 264 end do 265 do isou = 1, 3 266 do jsou = 1, 3 267 implmat2add(isou, jsou) = 0.0d0 268 end do 269 end do 270 271 impl_lin_cst = 0.0d0 272 impl_id_cst = 0.0d0 273 274 ! Pij 275 xprod(1,1) = produc(1, iel) 276 xprod(1,2) = produc(4, iel) 277 xprod(1,3) = produc(6, iel) 278 xprod(2,2) = produc(2, iel) 279 xprod(2,3) = produc(5, iel) 280 xprod(3,3) = produc(3, iel) 281 282 xprod(2,1) = xprod(1,2) 283 xprod(3,1) = xprod(1,3) 284 xprod(3,2) = xprod(2,3) 285 286 trprod = d1s2 * (xprod(1,1) + xprod(2,2) + xprod(3,3) ) 287 trrij = d1s2 * (cvara_var(1,iel) + cvara_var(2,iel) + cvara_var(3,iel)) 288 289 !-----> aII = aijaij 290 aii = 0.d0 291 aklskl = 0.d0 292 aiksjk = 0.d0 293 aikrjk = 0.d0 294 aikakj = 0.d0 295 296 ! aij 297 xaniso(1,1) = cvara_var(1,iel)/trrij - d2s3 298 xaniso(2,2) = cvara_var(2,iel)/trrij - d2s3 299 xaniso(3,3) = cvara_var(3,iel)/trrij - d2s3 300 xaniso(1,2) = cvara_var(4,iel)/trrij 301 xaniso(1,3) = cvara_var(6,iel)/trrij 302 xaniso(2,3) = cvara_var(5,iel)/trrij 303 xaniso(2,1) = xaniso(1,2) 304 xaniso(3,1) = xaniso(1,3) 305 xaniso(3,2) = xaniso(2,3) 306 307 ! Sij 308 xstrai(1,1) = gradv(1, 1, iel) 309 xstrai(1,2) = d1s2*(gradv(2, 1, iel)+gradv(1, 2, iel)) 310 xstrai(1,3) = d1s2*(gradv(3, 1, iel)+gradv(1, 3, iel)) 311 xstrai(2,1) = xstrai(1,2) 312 xstrai(2,2) = gradv(2, 2, iel) 313 xstrai(2,3) = d1s2*(gradv(3, 2, iel)+gradv(2, 3, iel)) 314 xstrai(3,1) = xstrai(1,3) 315 xstrai(3,2) = xstrai(2,3) 316 xstrai(3,3) = gradv(3, 3, iel) 317 318 ! omegaij 319 xrotac(1,1) = 0.d0 320 xrotac(1,2) = d1s2*(gradv(2, 1, iel)-gradv(1, 2, iel)) 321 xrotac(1,3) = d1s2*(gradv(3, 1, iel)-gradv(1, 3, iel)) 322 xrotac(2,1) = -xrotac(1,2) 323 xrotac(2,2) = 0.d0 324 xrotac(2,3) = d1s2*(gradv(3, 2, iel)-gradv(2, 3, iel)) 325 xrotac(3,1) = -xrotac(1,3) 326 xrotac(3,2) = -xrotac(2,3) 327 xrotac(3,3) = 0.d0 328 329 do ii=1,3 330 do jj = 1,3 331 ! aii = aij.aij 332 aii = aii+xaniso(ii,jj)*xaniso(ii,jj) 333 ! aklskl = aij.Sij 334 aklskl = aklskl + xaniso(ii,jj)*xstrai(ii,jj) 335 enddo 336 enddo 337 338 ! Computation of implicit components 339 340 sym_strain(1) = xstrai(1,1) 341 sym_strain(2) = xstrai(2,2) 342 sym_strain(3) = xstrai(3,3) 343 sym_strain(4) = xstrai(1,2) 344 sym_strain(5) = xstrai(2,3) 345 sym_strain(6) = xstrai(1,3) 346 347 do isou = 1, 6 348 matrn(isou) = cvara_var(isou,iel)/trrij 349 oo_matrn(isou) = 0.0d0 350 end do 351 352 ! Inversing the matrix 353 call symmetric_matrix_inverse(matrn, oo_matrn) 354 do isou = 1, 6 355 oo_matrn(isou) = oo_matrn(isou)/trrij 356 end do 357 358 ! Computing the maximal eigenvalue (in terms of norm!) of S 359 call calc_symtens_eigvals(sym_strain, eigen_vals) 360 eigen_max = maxval(abs(eigen_vals)) 361 362 ! Constant for the dissipation 363 ceps_impl = d1s3 * cvara_ep(iel) 364 365 ! Identity constant 366 impl_id_cst = -d1s3*crij2*min(trprod,0.0d0) 367 368 ! Linear constant 369 impl_lin_cst = eigen_max * ( & 370 (1.d0-crij2) ) ! Production + Phi2 371 372 do jsou = 1, 3 373 do isou = 1 ,3 374 iii = t2v(isou,jsou) 375 implmat2add(isou,jsou) = (1.d0-crij2)*xrotac(isou,jsou) & 376 + impl_lin_cst*deltij(iii) & 377 + impl_id_cst*d1s2*oo_matrn(iii) & 378 + ceps_impl*oo_matrn(iii) 379 end do 380 end do 381 382 impl_drsm(:,:) = 0.0d0 383 call reduce_symprod33_to_66(implmat2add, impl_drsm) 384 385 ! Rotating frame of reference => "absolute" vorticity 386 if (icorio.eq.1) then 387 do ii = 1, 3 388 do jj = 1, 3 389 xrotac(ii,jj) = xrotac(ii,jj) + matrot(ii,jj) 390 enddo 391 enddo 392 endif 393 394 do isou = 1, 6 395 iii = iv2t(isou) 396 jjj = jv2t(isou) 397 aiksjk = 0 398 aikrjk = 0 399 aikakj = 0 400 do kk = 1,3 401 ! aiksjk = aik.Sjk+ajk.Sik 402 aiksjk = aiksjk + xaniso(iii,kk)*xstrai(jjj,kk) & 403 + xaniso(jjj,kk)*xstrai(iii,kk) 404 ! aikrjk = aik.Omega_jk + ajk.omega_ik 405 aikrjk = aikrjk + xaniso(iii,kk)*xrotac(jjj,kk) & 406 + xaniso(jjj,kk)*xrotac(iii,kk) 407 ! aikakj = aik*akj 408 aikakj = aikakj + xaniso(iii,kk)*xaniso(kk,jjj) 409 enddo 410 411 ! Explicit terms 412 pij = (1.d0 - crij2)*xprod(iii,jjj) 413 phiij1 = -cvara_ep(iel)*crij1*xaniso(iii,jjj) 414 phiij2 = d2s3*crij2*trprod*deltij(isou) 415 epsij = -d2s3*cvara_ep(iel)*deltij(isou) 416 417 if (st_prv_id.ge.0) then 418 c_st_prv(isou,iel) = c_st_prv(isou,iel) & 419 + cromo(iel)*cell_f_vol(iel)*(pij+phiij1+phiij2+epsij) 420 else 421 smbr(isou,iel) = smbr(isou,iel) & 422 + cromo(iel)*cell_f_vol(iel)*(pij+phiij1+phiij2+epsij) 423 ! Implicit terms 424 rovsdt(isou,isou,iel) = rovsdt(isou,isou,iel) & 425 + cell_f_vol(iel)/trrij*crom(iel)*(crij1*cvara_ep(iel)) 426 427 ! Careful ! Inversion of the order of the coefficients since 428 ! rovsdt matrix is then used by a c function for the linear solving 429 do jsou = 1, 6 430 rovsdt(jsou,isou,iel) = rovsdt(jsou,isou,iel) + cell_f_vol(iel) & 431 *crom(iel) * impl_drsm(isou,jsou) 432 end do 433 434 endif 435 436 enddo 437 438enddo 439 440!=============================================================================== 441! Coriolis terms in the Phi1 and production 442!=============================================================================== 443 444if (icorio.eq.1 .or. iturbo.eq.1) then 445 446 do iel = 1, ncel 447 448 rot_id = icorio 449 if (iturbo.eq.1) rot_id = irotce(iel) 450 451 if (rot_id .lt. 1) cycle 452 453 call coriolis_t(rot_id, 1.d0, matrot) 454 455 cvara_r(1,1) = cvara_var(1,iel) 456 cvara_r(2,2) = cvara_var(2,iel) 457 cvara_r(3,3) = cvara_var(3,iel) 458 cvara_r(1,2) = cvara_var(4,iel) 459 cvara_r(2,3) = cvara_var(5,iel) 460 cvara_r(1,3) = cvara_var(6,iel) 461 cvara_r(2,1) = cvara_var(4,iel) 462 cvara_r(3,2) = cvara_var(5,iel) 463 cvara_r(3,1) = cvara_var(6,iel) 464 465 ! Compute Gij: (i,j) component of the Coriolis production 466 do isou = 1, 6 467 ii = iv2t(isou) 468 jj = jv2t(isou) 469 470 w2(isou,iel) = 0.d0 471 do kk = 1, 3 472 w2(isou,iel) = w2(isou,iel) - ccorio*( matrot(ii,kk)*cvara_r(jj,kk) & 473 + matrot(jj,kk)*cvara_r(ii,kk) ) 474 enddo 475 enddo 476 enddo 477 478 ! Coriolis contribution in the Phi1 term: (1-C2/2)Gij 479 if (icorio.eq.1) then 480 do iel = 1, ncel 481 do isou = 1, 6 482 w2(isou,iel) = crom(iel)*cell_f_vol(iel)*(1.d0 - 0.5d0*crij2)*w2(isou,iel) 483 enddo 484 enddo 485 endif 486 487 ! If source terms are extrapolated 488 if (st_prv_id.ge.0) then 489 do iel = 1, ncel 490 do isou = 1, 6 491 c_st_prv(isou,iel) = c_st_prv(isou,iel) + w2(isou,iel) 492 enddo 493 enddo 494 ! Otherwise, directly in smbr 495 else 496 do iel = 1, ncel 497 do isou = 1, 6 498 smbr(isou,iel) = smbr(isou,iel) + w2(isou,iel) 499 enddo 500 enddo 501 endif 502 503endif 504 505!=============================================================================== 506! Wall echo terms 507!=============================================================================== 508 509if (irijec.eq.1) then !todo 510 511 do iel = 1, ncel 512 do isou = 1, 6 513 w2(isou,iel) = 0.d0 514 enddo 515 enddo 516 517 call rijech2(produc, w2) 518 519 ! If we extrapolate the source terms: c_st_prv 520 if (st_prv_id.ge.0) then 521 do iel = 1, ncel 522 do isou = 1, 6 523 c_st_prv(isou,iel) = c_st_prv(isou,iel) + w2(isou,iel) 524 enddo 525 enddo 526 ! Otherwise smbr 527 else 528 do iel = 1, ncel 529 do isou = 1, 6 530 smbr(isou,iel) = smbr(isou,iel) + w2(isou,iel) 531 enddo 532 enddo 533 endif 534 535endif 536 537!=============================================================================== 538! Buoyancy source term 539!=============================================================================== 540 541if (igrari.eq.1) then 542 543 call field_get_id_try("rij_buoyancy", f_id) 544 if (f_id.ge.0) then 545 call field_get_val_v(f_id, cpro_buoyancy) 546 else 547 ! Allocate a work array 548 allocate(buoyancy(6,ncelet)) 549 cpro_buoyancy => buoyancy 550 endif 551 552 call rijthe2(gradro, cpro_buoyancy) 553 554 ! If we extrapolate the source terms: previous ST 555 if (st_prv_id.ge.0) then 556 do iel = 1, ncel 557 do isou = 1, 6 558 c_st_prv(isou,iel) = c_st_prv(isou,iel) + cpro_buoyancy(isou,iel) * cell_f_vol(iel) 559 enddo 560 enddo 561 ! Otherwise smbr 562 else 563 do iel = 1, ncel 564 do isou = 1, 6 565 smbr(isou,iel) = smbr(isou,iel) + cpro_buoyancy(isou,iel) * cell_f_vol(iel) 566 enddo 567 enddo 568 endif 569 570 ! Free memory 571 if (allocated(buoyancy)) deallocate(buoyancy) 572 573endif 574 575!=============================================================================== 576! Diffusion term (Daly Harlow: generalized gradient hypothesis method) 577!=============================================================================== 578 579! Symmetric tensor diffusivity (GGDH) 580if (iand(vcopt%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 581 582 call field_get_val_v(ivsten, visten) 583 584 do iel = 1, ncel 585 viscce(1,iel) = visten(1,iel) + viscl(iel) 586 viscce(2,iel) = visten(2,iel) + viscl(iel) 587 viscce(3,iel) = visten(3,iel) + viscl(iel) 588 viscce(4,iel) = visten(4,iel) 589 viscce(5,iel) = visten(5,iel) 590 viscce(6,iel) = visten(6,iel) 591 enddo 592 593 iwarnp = vcopt%iwarni 594 595 call vitens(viscce, iwarnp, weighf, weighb, viscf, viscb) 596 597! Scalar diffusivity 598else 599 600 do iel = 1, ncel 601 trrij = 0.5d0 * (cvara_var(1,iel) + cvara_var(2,iel) + cvara_var(3,iel)) 602 rctse = crom(iel) * csrij * trrij**2 / cvara_ep(iel) 603 w1(iel) = viscl(iel) + vcopt%idifft*rctse 604 enddo 605 606 imvisp = vcopt%imvisf 607 608 call viscfa(imvisp, w1, viscf, viscb) 609 610endif 611 612! Free memory 613 614deallocate(w1, w2) 615 616!-------- 617! Formats 618!-------- 619 620 1000 format(/,' Solving variable ', a8 ,/) 621 622!---- 623! End 624!---- 625 626return 627 628end subroutine 629