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 resrit.f90 28!> 29!> \brief This subroutine perform the solving of the transport equation 30!> of the turbulent heat fluxes. 31!------------------------------------------------------------------------------- 32 33!------------------------------------------------------------------------------- 34! Arguments 35!______________________________________________________________________________. 36! mode name role 37!______________________________________________________________________________! 38!> \param[in] nscal total number of scalars 39!> \param[in] iscal number of the scalar used 40!> \param[in] xcpp \f$ C_p \f$ 41!> \param[in,out] xut, xuta calculated variables at cell centers 42!> (at current and previous time steps) 43!> \param[in] dt time step (per cell) 44!> \param[in] gradv mean velocity gradient 45!> \param[in] gradt mean scalar gradient 46!> \param[in] grad_al alpha scalar gradient 47!______________________________________________________________________________! 48 49subroutine resrit & 50 ( nscal , & 51 iscal , xcpp , xut , xuta , & 52 dt , & 53 gradv , gradt , grad_al) 54 55!=============================================================================== 56! Module files 57!=============================================================================== 58 59use, intrinsic :: iso_c_binding 60 61use paramx 62use numvar 63use entsor 64use optcal 65use cstnum 66use cstphy 67use parall 68use period 69use field 70use mesh 71use cs_f_interfaces 72use cs_c_bindings 73 74!=============================================================================== 75 76implicit none 77 78! Arguments 79 80integer nscal , iscal 81 82double precision dt(ncelet) 83double precision xcpp(ncelet), xut(3,ncelet), xuta(3,ncelet) 84double precision gradv(3,3,ncelet) 85double precision gradt(3,ncelet) 86double precision grad_al(3,ncelet) 87 88! Local variables 89 90integer iel 91integer ii, ivar 92integer iflmas, iflmab 93integer iwarnp 94integer imvisp 95integer iescap 96integer st_prv_id 97integer ivisep, ifcvsl 98integer isou, jsou 99integer itt 100integer icvflb 101integer f_id 102integer keyvar, iut 103integer init 104integer kturt, turb_flux_model 105 106integer ivoid(1) 107 108double precision trrij 109double precision thets , thetv , thetp1 110double precision xttke , prdtl 111double precision grav(3) 112double precision xrij(3,3),phiith(3), phiitw(3) 113double precision xnal(3), xnoral 114double precision alpha, xttdrbt, xttdrbw 115double precision pk, gk, xxc1, xxc2, xxc3, imp_term 116double precision rctse, visls_0 117 118double precision rvoid(1) 119 120character(len=80) :: fname, name 121 122double precision, allocatable, dimension(:,:) :: viscce 123double precision, allocatable, dimension(:) :: viscb 124double precision, allocatable, dimension(:) :: viscf 125double precision, allocatable, dimension(:,:) :: weighf 126double precision, allocatable, dimension(:) :: weighb 127double precision, allocatable, dimension(:,:) :: smbrut 128double precision, allocatable, dimension(:,:,:) :: fimp 129double precision, allocatable, dimension(:) :: w1 130 131double precision, dimension(:,:), pointer :: coefav, cofafv, visten 132double precision, dimension(:,:,:), pointer :: coefbv, cofbfv 133double precision, dimension(:), pointer :: imasfl, bmasfl 134double precision, dimension(:), pointer :: crom, cpro_beta 135double precision, dimension(:), pointer :: cvar_al 136double precision, dimension(:), pointer :: cvar_ep 137double precision, dimension(:,:), pointer :: cvar_rij 138double precision, dimension(:), pointer :: cvar_tt, cvara_tt 139double precision, dimension(:), pointer :: viscl, visct, viscls, c_st_prv 140 141type(var_cal_opt) :: vcopt, vcopt_ut 142type(var_cal_opt), target :: vcopt_loc 143type(var_cal_opt), pointer :: p_k_value 144type(c_ptr) :: c_k_value 145 146!=============================================================================== 147 148!=============================================================================== 149! 1. Initialization 150!=============================================================================== 151 152! Allocate work arrays 153allocate(w1(ncelet)) 154allocate(viscce(6,ncelet)) 155allocate(smbrut(3,ncelet)) 156allocate(fimp(3,3,ncelet)) 157allocate(viscf(nfac), viscb(nfabor)) 158allocate(weighf(2,nfac)) 159allocate(weighb(nfabor)) 160 161if ((itytur.eq.2).or.(itytur.eq.5).or.(iturb.eq.60)) then 162 write(nfecra,*)'Utiliser un modele Rij avec ces modeles de thermiques'!FIXME 163 call csexit(1) 164endif 165 166call field_get_val_s(icrom, crom) 167call field_get_val_s(iviscl, viscl) 168call field_get_val_s(ivisct, visct) 169 170call field_get_val_s(ivarfl(iep), cvar_ep) 171call field_get_val_v(ivarfl(irij), cvar_rij) 172 173call field_get_key_int(ivarfl(iu), kimasf, iflmas) 174call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 175call field_get_val_s(iflmas, imasfl) 176call field_get_val_s(iflmab, bmasfl) 177 178if (ibeta.ge.0) then 179 call field_get_val_s(ibeta, cpro_beta) 180endif 181 182call field_get_val_v(ivsten, visten) 183 184ivar = isca(iscal) 185 186call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 187 188call field_get_key_id("variable_id", keyvar) 189! Get name of the reference scalar 190call field_get_name(ivarfl(ivar), fname) 191! Index of the corresponding turbulent flux 192call field_get_id(trim(fname)//'_turbulent_flux', f_id) 193! Set pointer values of turbulent fluxes 194call field_get_key_int(f_id, keyvar, iut) 195 196call field_get_key_struct_var_cal_opt(ivarfl(iut), vcopt_ut) 197 198if (vcopt%iwarni.ge.1) then 199 call field_get_name(ivarfl(ivar), name) 200 write(nfecra,1000) trim(name)//'_turbulent_flux' 201endif 202 203! Get the turbulent flux model 204call field_get_key_id('turbulent_flux_model', kturt) 205call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 206 207! S pour Source, V pour Variable 208thets = thetst 209thetv = vcopt%thetav 210 211call field_get_key_int(ivarfl(ivar), kstprv, st_prv_id) 212if (st_prv_id.ge.0) then 213 call field_get_val_s(st_prv_id, c_st_prv) 214else 215 c_st_prv=> null() 216endif 217 218call field_get_key_int (ivarfl(ivar), kivisl, ifcvsl) 219if (ifcvsl .ge. 0) then 220 call field_get_val_s(ifcvsl, viscls) 221else 222 call field_get_key_double(ivarfl(ivar), kvisl0, visls_0) 223endif 224 225do iel = 1, ncelet 226 do isou = 1, 3 227 smbrut(isou,iel) = 0.d0 228 do jsou = 1, 3 229 fimp(isou,jsou,iel) = 0.d0 230 enddo 231 enddo 232enddo 233 234! Find the corresponding variance of the scalar iscal 235itt = -1 236if ((abs(gx)+abs(gy)+abs(gz)).gt.0) then 237 grav(1) = gx 238 grav(2) = gy 239 grav(3) = gz 240 do ii = 1, nscal 241 if (iscavr(ii).eq.iscal) itt = ii 242 enddo 243endif 244 245!=============================================================================== 246! 2. Mass source terms FIXME 247!=============================================================================== 248 249if (st_prv_id.ge.0) then 250 do iel = 1, ncel 251 do isou = 1,3 252 smbrut(isou,iel) = fimp(isou,isou,iel)*xuta(isou,iel) 253 fimp(isou,isou,iel) = - thetv*fimp(isou,isou,iel) 254 enddo 255 enddo 256! If we do not extrapolate the source terms 257else 258 do iel = 1, ncel 259 do isou = 1, 3 260 ! User source term 261 smbrut(isou,iel) = smbrut(isou,iel) + fimp(isou,isou,iel)*xuta(isou,iel) 262 ! Diagonal 263 fimp(isou,isou,iel) = max(-fimp(isou,isou,iel),zero) 264 enddo 265 enddo 266endif 267 268!=============================================================================== 269! 3. Unsteady term 270!=============================================================================== 271 272do iel = 1, ncel 273 do isou = 1, 3 274 fimp(isou,isou,iel) = fimp(isou,isou,iel) & 275 + vcopt%istat*(crom(iel)/dt(iel))*cell_f_vol(iel) 276 enddo 277enddo 278 279!=============================================================================== 280! 4. Right Hand Side of the thermal fluxes: 281! rho*(Pit + Git + Phi*_it - eps_it) 282!=============================================================================== 283 284if (itt.gt.0) then 285 call field_get_val_s(ivarfl(isca(itt)), cvar_tt) 286 call field_get_val_prev_s(ivarfl(isca(itt)), cvara_tt) 287endif 288if (turb_flux_model.eq.31) then 289 ! Name of the scalar 290 call field_get_name(ivarfl(isca(iscal)), fname) 291 292 ! Index of the corresponding turbulent flux 293 call field_get_id(trim(fname)//'_alpha', f_id) 294 295 call field_get_val_s(f_id, cvar_al) 296endif 297 298do iel = 1, ncel 299 300 xrij(1,1) = cvar_rij(1,iel) 301 xrij(2,2) = cvar_rij(2,iel) 302 xrij(3,3) = cvar_rij(3,iel) 303 xrij(1,2) = cvar_rij(4,iel) 304 xrij(2,3) = cvar_rij(5,iel) 305 xrij(1,3) = cvar_rij(6,iel) 306 307 xrij(2,1) = xrij(1,2) 308 xrij(3,1) = xrij(1,3) 309 xrij(3,2) = xrij(2,3) 310 311 if (ifcvsl.ge.0) then 312 prdtl = viscl(iel)*xcpp(iel)/viscls(iel) 313 else 314 prdtl = viscl(iel)*xcpp(iel)/visls_0 315 endif 316 317 trrij = 0.5d0*(xrij(1,1) + xrij(2,2) + xrij(3,3)) 318 ! --- Compute Durbin time scheme 319 xttke = trrij/cvar_ep(iel) 320 321 if (turb_flux_model.eq.31) then 322 alpha = cvar_al(iel) 323 !FIXME Warning / rhebdfm**0.5 compared to F Dehoux 324 xttdrbt = xttke * sqrt( (1.d0-alpha)*prdtl/ rhebdfm + alpha ) 325 xttdrbw = xttdrbt * sqrt(rhebdfm/prdtl) ! And so multiplied by (R/Prandt)^0.5 326 327 ! Compute the unite normal vector 328 xnoral = & 329 ( grad_al(1, iel)*grad_al(1, iel) & 330 + grad_al(2, iel)*grad_al(2, iel) & 331 + grad_al(3, iel)*grad_al(3, iel)) 332 xnoral = sqrt(xnoral) 333 334 if (xnoral.le.epzero/cell_f_vol(iel)**(1.d0/3.d0)) then 335 xnal(1) = 0.d0 336 xnal(2) = 0.d0 337 xnal(3) = 0.d0 338 else 339 xnal(1) = grad_al(1, iel)/xnoral 340 xnal(2) = grad_al(2, iel)/xnoral 341 xnal(3) = grad_al(3, iel)/xnoral 342 endif 343 344 ! Production and buoyancy for TKE 345 pk = -( xrij(1,1)*gradv(1,1,iel) +& 346 xrij(1,2)*gradv(1,2,iel) +& 347 xrij(1,3)*gradv(1,3,iel) +& 348 xrij(2,1)*gradv(2,1,iel) +& 349 xrij(2,2)*gradv(2,2,iel) +& 350 xrij(2,3)*gradv(2,3,iel) +& 351 xrij(3,1)*gradv(3,1,iel) +& 352 xrij(3,2)*gradv(3,2,iel) +& 353 xrij(3,3)*gradv(3,3,iel) ) 354 355 if (ibeta.ge.0) then 356 gk = cpro_beta(iel)*(xuta(1,iel)*gx + xuta(2,iel)*gy + xuta(3,iel)*gz)!FIXME make buoyant term coherent elsewhere 357 else 358 gk = 0.d0 359 endif 360 xxc1 = 1.d0+2.d0*(1.d0 - cvar_al(iel))*(pk+gk)/cvar_ep(iel) 361 xxc2 = 0.5d0*(1.d0+1.d0/prdtl)*(1.d0-0.3d0*(1.d0 - cvar_al(iel)) & 362 *(pk+gk)/cvar_ep(iel)) 363 xxc3 = xxc2 364 365 else 366 xttdrbt = xttke 367 xttdrbw = xttke 368 alpha = 1.d0 369 xxc1 = 0.d0 370 xxc2 = 0.d0 371 xxc3 = 0.d0 372 xnal(1) = 0.d0 373 xnal(2) = 0.d0 374 xnal(3) = 0.d0 375 endif 376 377 do isou = 1, 3 378 phiith(isou) = -c1trit / xttdrbt * xuta(isou,iel) & 379 + c2trit*(xuta(1,iel)*gradv(1,isou,iel) & 380 +xuta(2,iel)*gradv(2,isou,iel) & 381 +xuta(3,iel)*gradv(3,isou,iel)) & 382 + c4trit*(-xrij(isou,1)*gradt(1,iel) & 383 -xrij(isou,2)*gradt(2,iel) & 384 -xrij(isou,3)*gradt(3,iel)) 385 if (itt.gt.0.and.ibeta.ge.0) then 386 phiith(isou) = phiith(isou) & 387 + c3trit*(cpro_beta(iel)*grav(isou)*cvar_tt(iel)) 388 endif 389 390 phiitw(isou) = -1.d0/xttdrbw * xxc1 * & !FIXME full implicit 391 ( xuta(1,iel)*xnal(1)*xnal(isou) & 392 + xuta(2,iel)*xnal(2)*xnal(isou) & 393 + xuta(3,iel)*xnal(3)*xnal(isou)) 394 395 ! Pressure/thermal fluctuation correlation term 396 !---------------------------------------------- 397 smbrut(isou,iel) = smbrut(isou,iel) + & 398 cell_f_vol(iel)*crom(iel)*( alpha * phiith(isou) & 399 + (1.d0 - alpha) * phiitw(isou)) 400 401 imp_term = max(cell_f_vol(iel)*crom(iel)*( & 402 alpha * (c1trit/xttdrbt-c2trit*gradv(isou,isou,iel)) & 403 + (1.d0-alpha) * (xxc1*xnal(isou)*xnal(isou) / xttdrbw) & 404 ), 0.d0) 405 406 fimp(isou,isou,iel) = fimp(isou,isou,iel) + imp_term 407 408 ! Production terms 409 !----------------- 410 smbrut(isou,iel) = smbrut(isou,iel) & 411 + cell_f_vol(iel)*crom(iel) & 412 ! Production term due to the mean velcoity 413 *( -gradv(1,isou,iel)*xuta(1,iel) & 414 -gradv(2,isou,iel)*xuta(2,iel) & 415 -gradv(3,isou,iel)*xuta(3,iel) & 416 ! Production term due to the mean temperature 417 -xrij(1,isou)*gradt(1,iel) & 418 -xrij(2,isou)*gradt(2,iel) & 419 -xrij(3,isou)*gradt(3,iel) & 420 ) 421 422 ! Production term due to the gravity 423 if (itt.gt.0.and.ibeta.ge.0) then 424 smbrut(isou,iel) = smbrut(isou,iel) & 425 + cell_f_vol(iel)*crom(iel)*( & 426 -grav(isou)*cpro_beta(iel)*cvara_tt(iel)) 427 endif 428 429 ! Dissipation (Wall term only because "h" term is zero 430 smbrut(isou,iel) = smbrut(isou,iel) - & 431 cell_f_vol(iel)*crom(iel)*(1.d0 - alpha) / xttdrbw * & 432 ( xxc2 * xuta(isou, iel) & 433 + xxc3*( xuta(1,iel)*xnal(1)*xnal(isou) & 434 + xuta(2,iel)*xnal(2)*xnal(isou) & 435 + xuta(3,iel)*xnal(3)*xnal(isou))) 436 437 ! TODO we can implicite more terms 438 imp_term = max(cell_f_vol(iel)*crom(iel)*(1.d0 - alpha) / xttdrbw * & 439 ( xxc2 + xxc3 * xnal(isou)*xnal(isou)), 0.d0) 440 fimp(isou,isou,iel) = fimp(isou,isou,iel) + imp_term 441 442 enddo 443enddo 444 445!=============================================================================== 446! 5. Tensor diffusion 447!=============================================================================== 448! Symmetric tensor diffusivity (GGDH) 449if (iand(vcopt_ut%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 450 451 call field_get_key_double(ivarfl(isca(iscal)), kvisl0, visls_0) 452 453 do iel = 1, ncel 454 455 if (ifcvsl.ge.0) then 456 prdtl = viscl(iel)*xcpp(iel)/viscls(iel) 457 else 458 prdtl = viscl(iel)*xcpp(iel)/visls_0 459 endif 460 461 do isou = 1, 6 462 if (isou.le.3) then 463 viscce(isou,iel) = 0.5d0*(viscl(iel)*(1.d0+1.d0/prdtl)) & 464 + ctheta(iscal)*visten(isou,iel)/csrij 465 else 466 viscce(isou,iel) = ctheta(iscal)*visten(isou,iel)/csrij 467 endif 468 enddo 469 enddo 470 471 iwarnp = vcopt%iwarni 472 473 call vitens & 474 ( viscce , iwarnp , & 475 weighf , weighb , & 476 viscf , viscb ) 477 478! Scalar diffusivity 479else 480 481 do iel = 1, ncel 482 rctse = ctheta(iscal) * visct(iel) / cmu 483 w1(iel) = viscl(iel) + vcopt%idifft*rctse 484 enddo 485 486 imvisp = vcopt%imvisf 487 488 call viscfa & 489 ( imvisp , & 490 w1 , & 491 viscf , viscb ) 492 493endif 494 495!=============================================================================== 496! 6. Vectorial solving of the turbulent thermal fluxes 497!=============================================================================== 498 499if (st_prv_id.ge.0) then 500 thetp1 = 1.d0 + thets 501 do iel = 1, ncel 502 do isou = 1, 3 503 smbrut(isou,iel) = smbrut(isou,iel) + thetp1*c_st_prv(iel) !FIXME 504 enddo 505 enddo 506endif 507 508! Name of the scalar ivar 509call field_get_name(ivarfl(ivar), fname) 510 511! Index of the corresponding turbulent flux 512call field_get_id(trim(fname)//'_turbulent_flux', f_id) 513 514call field_get_coefa_v(f_id,coefav) 515call field_get_coefb_v(f_id,coefbv) 516call field_get_coefaf_v(f_id,cofafv) 517call field_get_coefbf_v(f_id,cofbfv) 518 519iescap = 0 520 521! We do not take into account transpose of grad 522ivisep = 0 523 524! all boundary convective flux with upwind 525icvflb = 0 526init = 1 527 528vcopt_loc = vcopt 529 530vcopt_loc%istat = -1 531vcopt_loc%idifft = -1 532vcopt_loc%iwgrec = 0 533vcopt_loc%thetav = thetv 534vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 535 536p_k_value => vcopt_loc 537c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 538 539call cs_equation_iterative_solve_vector & 540 ( idtvar , init , & 541 f_id , c_null_char , & 542 ivisep , iescap , c_k_value , & 543 xuta , xuta , & 544 coefav , coefbv , cofafv , cofbfv , & 545 imasfl , bmasfl , viscf , & 546 viscb , viscf , viscb , rvoid , & 547 rvoid , viscce , weighf , weighb , & 548 icvflb , ivoid , & 549 fimp , smbrut , xut , rvoid ) 550 551!=============================================================================== 552! 7. Writings 553!=============================================================================== 554 555! Free memory 556deallocate(viscce) 557deallocate(viscf, viscb) 558deallocate(smbrut) 559deallocate(fimp) 560deallocate(weighf, weighb) 561deallocate(w1) 562!-------- 563! Formats 564!-------- 565 566 1000 format(/,' Solving variable ',A23 ,/) 567 568!---- 569! End 570!---- 571 572return 573end subroutine 574