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 turrij.f90 28!> 29!> \brief Solving the \f$ R_{ij} - \epsilon \f$ for incompressible flows or 30!> slightly compressible flows for one time step. 31!> 32!> Please refer to the 33!> <a href="../../theory.pdf#rijeps"><b>\f$ R_{ij} - \epsilon \f$ model</b></a> 34!> section of the theory guide for more informations, as well as the 35!> <a href="../../theory.pdf#turrij"><b>turrij</b></a> section. 36!------------------------------------------------------------------------------- 37 38!------------------------------------------------------------------------------- 39! Arguments 40!______________________________________________________________________________. 41! mode name role 42!______________________________________________________________________________! 43!> \param[in] nvar total number of variables 44!> \param[in] nscal total number of scalars 45!> \param[in] ncepdp number of cells with head loss 46!> \param[in] ncesmp number of cells with mass source term 47!> \param[in] icepdc index of the ncepdp cells with head loss 48!> \param[in] icetsm index of cells with mass source term 49!> \param[in] itypsm mass source type for the variables 50!> (cf. cs_user_mass_source_terms) 51!> \param[in] dt time step (per cell) 52!> \param[in] tslagr coupling term of the lagangian module 53!> \param[in] ckupdc work array for the head loss 54!> \param[in] smacel values of the variables associated to the 55!> mass source 56!> (for ivar=ipr, smacel is the mass flux) 57!_______________________________________________________________________________ 58 59subroutine turrij & 60 ( nvar , nscal , ncepdp , ncesmp , & 61 icepdc , icetsm , itypsm , & 62 dt , & 63 tslagr , & 64 ckupdc , smacel ) 65 66!=============================================================================== 67 68!=============================================================================== 69! Module files 70!=============================================================================== 71 72use paramx 73use lagran 74use numvar 75use entsor 76use cstphy 77use cstnum 78use optcal 79use ppincl 80use mesh 81use field 82use field_operator 83use cs_c_bindings 84 85!=============================================================================== 86 87implicit none 88 89! Arguments 90 91integer nvar , nscal 92integer ncepdp , ncesmp 93 94integer icepdc(ncepdp) 95integer icetsm(ncesmp) 96 97integer, dimension(ncesmp,nvar), target :: itypsm 98 99double precision dt(ncelet) 100double precision ckupdc(6,ncepdp) 101 102double precision, dimension(ncesmp,nvar), target :: smacel 103double precision, dimension(ncelet,ntersl), target :: tslagr 104double precision, dimension(:), pointer :: bromo, cromo 105 106! Local variables 107 108integer ifac , iel , ivar , isou, jsou 109integer iflmas, iflmab 110integer inc , iccocg 111integer iwarnp, iclip 112integer imrgrp, nswrgp, imligp 113integer f_id0 , f_id, st_prv_id 114integer iprev 115integer key_t_ext_id 116integer iroext 117integer f_id_phij 118integer icvflb 119integer ivoid(1) 120double precision epsrgp, climgp 121double precision rhothe 122double precision utaurf, ut2, ypa, ya, tke, xunorm, limiter, nu0, alpha 123double precision xnoral, xnal(3) 124double precision k, p, thets, thetv, tuexpr, thetp1 125double precision d1s3, d2s3 126 127double precision, allocatable, dimension(:) :: viscf, viscb 128double precision, allocatable, dimension(:) :: smbr, rovsdt 129double precision, allocatable, dimension(:,:,:) :: gradv 130double precision, allocatable, dimension(:,:), target :: produc 131double precision, allocatable, dimension(:,:) :: gradro 132double precision, allocatable, dimension(:,:) :: grad 133double precision, allocatable, dimension(:,:) :: smbrts, gatinj 134double precision, allocatable, dimension(:,:,:) :: rovsdtts 135double precision, allocatable, dimension(:,:) :: viscce 136double precision, allocatable, dimension(:,:) :: weighf 137double precision, allocatable, dimension(:) :: weighb 138 139double precision, pointer, dimension(:) :: tslagi 140double precision, dimension(:,:), pointer :: coefau 141double precision, dimension(:,:,:), pointer :: coefbu 142double precision, dimension(:), pointer :: brom, crom 143double precision, dimension(:), pointer :: cvara_scalt 144double precision, dimension(:), pointer :: cvar_ep, cvar_al 145double precision, dimension(:,:), pointer :: cvara_rij, cvar_rij, vel 146double precision, dimension(:,:), pointer :: c_st_prv, lagr_st_rij 147double precision, dimension(:,:), pointer :: cpro_produc 148double precision, dimension(:,:), pointer :: cpro_press_correl 149double precision, dimension(:), pointer :: cpro_beta 150 151double precision, dimension(:), pointer :: imasfl, bmasfl 152double precision, dimension(:,:), pointer :: coefa_rij, cofaf_rij 153double precision, dimension(:,:,:), pointer :: coefb_rij, cofbf_rij 154 155type(var_cal_opt) :: vcopt 156type(var_cal_opt), target :: vcopt_loc 157type(var_cal_opt), pointer :: p_k_value 158type(c_ptr) :: c_k_value 159 160!=============================================================================== 161 162!=============================================================================== 163! 1. Initialization 164!=============================================================================== 165 166call field_get_coefa_v(ivarfl(iu), coefau) 167call field_get_coefb_v(ivarfl(iu), coefbu) 168 169! Allocate temporary arrays for the turbulence resolution 170allocate(viscf(nfac), viscb(nfabor)) 171allocate(smbr(ncelet), rovsdt(ncelet)) 172allocate(gradv(3, 3, ncelet)) 173allocate(smbrts(6,ncelet)) 174allocate(rovsdtts(6,6,ncelet)) 175 176! Allocate other arrays, depending on user options 177call field_get_id_try("rij_production", f_id) 178if (f_id.ge.0) then 179 call field_get_val_v(f_id, cpro_produc) 180else 181 allocate(produc(6,ncelet)) 182 cpro_produc => produc 183endif 184 185call field_get_id_try("rij_pressure_strain_correlation", f_id_phij) 186if (f_id_phij.ge.0) then 187 call field_get_val_v(f_id_phij, cpro_press_correl) 188endif 189call field_get_val_s(ivarfl(iep), cvar_ep) 190 191call field_get_val_s(icrom, crom) 192call field_get_val_s(ibrom, brom) 193 194call field_get_val_v(ivarfl(irij), cvar_rij) 195call field_get_val_prev_v(ivarfl(irij), cvara_rij) 196call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt) 197 198thets = thetst 199thetv = vcopt%thetav 200 201if (vcopt%iwarni.ge.1) then 202 if (iturb.eq.30) then 203 write(nfecra,1000) 204 elseif (iturb.eq.31) then 205 write(nfecra,1001) 206 else 207 write(nfecra,1002) 208 endif 209endif 210 211! Time extrapolation? 212call field_get_key_id("time_extrapolated", key_t_ext_id) 213 214call field_get_key_int(ivarfl(irij), kstprv, st_prv_id) 215if (st_prv_id .ge. 0) then 216 call field_get_val_v(st_prv_id, c_st_prv) 217else 218 c_st_prv=> null() 219endif 220 221! Mass flux 222call field_get_key_int(ivarfl(iu), kimasf, iflmas) 223call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 224call field_get_val_s(iflmas, imasfl) 225call field_get_val_s(iflmab, bmasfl) 226 227!=============================================================================== 228! 1.1 Advanced init for EBRSM 229!=============================================================================== 230 231! Automatic reinitialization at the end of the first iteration: 232! wall distance y^+ is computed with -C log(1-alpha), where C=CL*Ceta*L*kappa, 233! then y so we have an idea of the wall distance in complex geometries. 234! Then U is initialized with a Reichard layer, 235! Epsilon by 1/(kappa y), clipped next to the wall at its value for y^+=15 236! k is given by a blending between eps/(2 nu)*y^2 and utau/sqrt(Cmu). 237! The blending function is chosen so that the asymptotic behavior 238! and give the correct peak of k. 239 240!TODO FIXME Are the BC uncompatible? 241if (ntcabs.eq.1.and.reinit_turb.eq.1.and.iturb.eq.32) then 242 243 allocate(grad(3,ncelet)) 244 245 ! Compute the gradient of Alpha 246 iprev = 0 247 inc = 1 248 iccocg = 1 249 250 call field_gradient_scalar(ivarfl(ial), iprev, 0, inc, iccocg, grad) 251 252 call field_get_val_v(ivarfl(irij), cvar_rij) 253 utaurf=0.05d0*uref 254 255 call field_get_val_s(ivarfl(iep), cvar_ep) 256 call field_get_val_s(ivarfl(ial), cvar_al) 257 call field_get_val_v(ivarfl(iu), vel) 258 259 nu0 = viscl0/ro0 260 261 do iel = 1, ncel 262 ! Compute the velocity magnitude 263 xunorm = vel(1,iel)**2 + vel(2,iel)**2 + vel(3,iel)**2 264 xunorm = sqrt(xunorm) 265 266 ! y+ is bounded by 400, because in the Reichard profile, 267 ! it corresponds to saturation (u>uref) 268 cvar_al(iel) = max(min(cvar_al(iel),(1.d0-exp(-400.d0/50.d0))),0.d0) 269 ! Compute the magnitude of the Alpha gradient 270 xnoral = ( grad(1,iel)*grad(1,iel) & 271 + grad(2,iel)*grad(2,iel) & 272 + grad(3,iel)*grad(3,iel) ) 273 xnoral = sqrt(xnoral) 274 ! Compute the unitary vector of Alpha 275 if (xnoral.le.epzero/cell_f_vol(iel)**(1.d0/3.d0)) then 276 xnal(1) = 1.d0/sqrt(3.d0) 277 xnal(2) = 1.d0/sqrt(3.d0) 278 xnal(3) = 1.d0/sqrt(3.d0) 279 else 280 xnal(1) = grad(1,iel)/xnoral 281 xnal(2) = grad(2,iel)/xnoral 282 xnal(3) = grad(3,iel)/xnoral 283 endif 284 285 alpha = cvar_al(iel) 286 287 ! Compute YA, therefore alpha is given by 1-exp(-YA/(50 nu/utau)) 288 ! NB: y^+ = 50 give the best compromise 289 ya = -dlog(1.d0-cvar_al(iel))*50.d0*nu0/utaurf 290 ypa = ya/(nu0/utaurf) 291 ! Velocity magnitude is imposed (limitted only), the direction is 292 ! conserved 293 if (xunorm.le.1.d-12*uref) then 294 limiter = 1.d0 295 else 296 limiter = min( utaurf/xunorm*(2.5d0*dlog(1.d0+0.4d0*ypa) & 297 +7.8d0*(1.d0-dexp(-ypa/11.d0) & 298 -(ypa/11.d0)*dexp(-0.33d0*ypa))), & 299 1.d0) 300 endif 301 302 vel(1,iel) = limiter*vel(1,iel) 303 vel(2,iel) = limiter*vel(2,iel) 304 vel(3,iel) = limiter*vel(3,iel) 305 306 ut2 = 0.05d0*uref 307 cvar_ep(iel) = utaurf**3*min(1.d0/(xkappa*15.d0*nu0/utaurf), & 308 1.d0/(xkappa*ya)) 309 tke = cvar_ep(iel)/2.d0/nu0*ya**2 & 310 * exp(-ypa/25.d0)**2 & 311 + ut2**2/0.3d0*(1.d0-exp(-ypa/25.d0))**2 312 313 cvar_rij(1,iel) = alpha**3 *2.d0/3.d0 *tke & 314 + (1.d0-alpha**3)*(1.d0-xnal(1)**2)*tke 315 cvar_rij(2,iel) = alpha**3 *2.d0/3.d0 *tke & 316 + (1.d0-alpha**3)*(1.d0-xnal(2)**2)*tke 317 cvar_rij(3,iel) = alpha**3 *2.d0/3.d0 *tke & 318 + (1.d0-alpha**3)*(1.d0-xnal(3)**2)*tke 319 cvar_rij(4,iel) = -(1.d0-alpha**3)*(xnal(1)*xnal(2))*tke 320 cvar_rij(5,iel) = -(1.d0-alpha**3)*(xnal(2)*xnal(3))*tke 321 cvar_rij(6,iel) = -(1.d0-alpha**3)*(xnal(1)*xnal(3))*tke 322 enddo 323 324 call field_current_to_previous(ivarfl(iu)) 325 call field_current_to_previous(ivarfl(irij)) 326 327 deallocate(grad) 328endif 329 330!=============================================================================== 331! 2.1 Compute the velocity gradient 332!=============================================================================== 333 334inc = 1 335iprev = 1 336 337call field_gradient_vector(ivarfl(iu), iprev, 0, inc, gradv) 338 339!=============================================================================== 340! 2.2 Compute the production term for Rij 341!=============================================================================== 342 343do iel = 1, ncel 344 345 ! Pij = - (Rik dUk/dXj + dUk/dXi Rkj) 346 ! Pij is stored as (P11, P22, P33, P12, P23, P13) 347 cpro_produc(1,iel) = & 348 - 2.0d0*(cvara_rij(1,iel)*gradv(1, 1, iel) + & 349 cvara_rij(4,iel)*gradv(2, 1, iel) + & 350 cvara_rij(6,iel)*gradv(3, 1, iel) ) 351 352 cpro_produc(4,iel) = & 353 - (cvara_rij(4,iel)*gradv(1, 1, iel) + & 354 cvara_rij(2,iel)*gradv(2, 1, iel) + & 355 cvara_rij(5,iel)*gradv(3, 1, iel) ) & 356 - (cvara_rij(1,iel)*gradv(1, 2, iel) + & 357 cvara_rij(4,iel)*gradv(2, 2, iel) + & 358 cvara_rij(6,iel)*gradv(3, 2, iel) ) 359 360 cpro_produc(6,iel) = & 361 - (cvara_rij(6,iel)*gradv(1, 1, iel) + & 362 cvara_rij(5,iel)*gradv(2, 1, iel) + & 363 cvara_rij(3,iel)*gradv(3, 1, iel) ) & 364 - (cvara_rij(1,iel)*gradv(1, 3, iel) + & 365 cvara_rij(4,iel)*gradv(2, 3, iel) + & 366 cvara_rij(6,iel)*gradv(3, 3, iel) ) 367 368 cpro_produc(2,iel) = & 369 - 2.0d0*(cvara_rij(4,iel)*gradv(1, 2, iel) + & 370 cvara_rij(2,iel)*gradv(2, 2, iel) + & 371 cvara_rij(5,iel)*gradv(3, 2, iel) ) 372 373 cpro_produc(5,iel) = & 374 - (cvara_rij(6,iel)*gradv(1, 2, iel) + & 375 cvara_rij(5,iel)*gradv(2, 2, iel) + & 376 cvara_rij(3,iel)*gradv(3, 2, iel) ) & 377 - (cvara_rij(4,iel)*gradv(1, 3, iel) + & 378 cvara_rij(2,iel)*gradv(2, 3, iel) + & 379 cvara_rij(5,iel)*gradv(3, 3, iel) ) 380 381 cpro_produc(3,iel) = & 382 - 2.0d0*(cvara_rij(6,iel)*gradv(1, 3, iel) + & 383 cvara_rij(5,iel)*gradv(2, 3, iel) + & 384 cvara_rij(3,iel)*gradv(3, 3, iel) ) 385 386enddo 387 388!=============================================================================== 389! 2.3 Compute the pressure correlation term for Rij 390!=============================================================================== 391! Phi,ij = Phi1,ij+Phi2,ij 392! Phi,ij = -C1 k/eps (Rij-2/3k dij) - C2 (Pij-2/3P dij) 393! Phi,ij is stored as (Phi11, Phi22, Phi33, Phi12, Phi23, Phi13) 394 395! TODO : coherent with the model 396if (f_id_phij.ge.0) then 397 d1s3 = 1.0d0/3.0d0 398 d2s3 = 2.0d0/3.0d0 399 do iel = 1, ncel 400 k=0.5*(cvara_rij(1,iel)+cvara_rij(2,iel)+cvara_rij(3,iel)) 401 p=0.5*(cpro_produc(1,iel)+cpro_produc(2,iel)+cpro_produc(3,iel)) 402 do isou=1,3 403 cpro_press_correl(isou, iel) = -crij1*cvar_ep(iel)/k*(cvara_rij(isou,iel)-d2s3*k) & 404 -crij2*(cpro_produc(isou,iel)-d2s3*p) 405 enddo 406 do isou=4,6 407 cpro_press_correl(isou, iel) = -crij1*cvar_ep(iel)/k*(cvara_rij(isou,iel)) & 408 -crij2*(cpro_produc(isou,iel)) 409 enddo 410 enddo 411endif 412 413!=============================================================================== 414! 3. Compute the density gradient for buoyant terms 415!=============================================================================== 416 417! Note that the buoyant term is normally expressed in temr of 418! (u'T') or (u'rho') here modelled with a GGDH: 419! (u'rho') = C * k/eps * R_ij Grad_j(rho) 420 421! Buoyant term for the Atmospheric module 422! (function of the potential temperature) 423if (igrari.eq.1 .and. ippmod(iatmos).ge.1) then 424 ! Allocate a temporary array for the gradient calculation 425 ! Warning, grad(theta) here 426 allocate(gradro(3,ncelet)) 427 428 call field_get_val_s(icrom, cromo) 429 430 call field_get_val_prev_s(ivarfl(isca(iscalt)), cvara_scalt) 431 432 inc = 1 433 iccocg = 1 434 435 call field_gradient_scalar(ivarfl(isca(iscalt)), 1, 0, inc, iccocg, gradro) 436 437 ! gradro stores: - rho grad(theta)/theta 438 ! grad(rho) and grad(theta) have opposite signs 439 do iel = 1, ncel 440 rhothe = cromo(iel)/cvara_scalt(iel) 441 gradro(1, iel) = -rhothe*gradro(1, iel) 442 gradro(2, iel) = -rhothe*gradro(2, iel) 443 gradro(3, iel) = -rhothe*gradro(3, iel) 444 enddo 445 446else if (igrari.eq.1) then 447 ! Allocate a temporary array for the gradient calculation 448 allocate(gradro(3,ncelet)) 449 450 ! Boussinesq approximation, only for the thermal scalar for the moment 451 if (idilat.eq.0) then 452 453 iccocg = 1 454 inc = 1 455 456 ! Use the current value... 457 call field_gradient_scalar(ivarfl(isca(iscalt)), 0, 0, inc, iccocg, gradro) 458 459 !FIXME make it dependant on the scalar and use is_buoyant field 460 call field_get_val_s(ibeta, cpro_beta) 461 462 ! gradro stores: - beta grad(T) 463 ! grad(rho) and grad(T) have opposite signs 464 do iel = 1, ncel 465 gradro(1, iel) = - ro0 * cpro_beta(iel) * gradro(1, iel) 466 gradro(2, iel) = - ro0 * cpro_beta(iel) * gradro(2, iel) 467 gradro(3, iel) = - ro0 * cpro_beta(iel) * gradro(3, iel) 468 enddo 469 470 else 471 472 ! Boundary conditions: Dirichlet romb 473 ! We use viscb to store the relative coefficient of rom 474 ! We impose in Dirichlet (coefa) the value romb 475 476 do ifac = 1, nfabor 477 viscb(ifac) = 0.d0 478 enddo 479 480 ! The choice below has the advantage to be simple 481 482 imrgrp = vcopt%imrgra 483 nswrgp = vcopt%nswrgr 484 imligp = vcopt%imligr 485 iwarnp = vcopt%iwarni 486 epsrgp = vcopt%epsrgr 487 climgp = vcopt%climgr 488 489 f_id0 = -1 490 iccocg = 1 491 492 call field_get_key_int(icrom, key_t_ext_id, iroext) 493 ! If we extrapolate the source terms and rho, we use cpdt rho^n 494 if(isto2t.gt.0.and.iroext.gt.0) then 495 call field_get_val_prev_s(icrom, cromo) 496 call field_get_val_prev_s(ibrom, bromo) 497 else 498 call field_get_val_s(icrom, cromo) 499 call field_get_val_s(ibrom, bromo) 500 endif 501 502 call gradient_s(f_id0, imrgrp, inc, iccocg, nswrgp, imligp, & 503 iwarnp, epsrgp, climgp, cromo, bromo, viscb, & 504 gradro) 505 506 endif 507endif 508 509!=============================================================================== 510! 4.1 Prepare to solve Rij 511! We solve the equation in a routine similar to covofi.f90 512!=============================================================================== 513 514!=============================================================================== 515! 4.1.1 Source terms for Rij 516!=============================================================================== 517 518do iel = 1, ncel 519 do isou = 1 ,6 520 smbrts(isou,iel) = 0.d0 521 enddo 522enddo 523do iel = 1, ncel 524 do isou = 1, 6 525 do jsou = 1, 6 526 rovsdtts(isou,jsou,iel) = 0.d0 527 enddo 528 enddo 529enddo 530 531! User source terms 532!------------------ 533 534call cs_user_turbulence_source_terms2(nvar, nscal, ncepdp, ncesmp, & 535 ivarfl(irij), & 536 icepdc, icetsm, itypsm, & 537 ckupdc, smacel, & 538 smbrts, rovsdtts) 539 540! C version 541call user_source_terms(ivarfl(irij), smbrts, rovsdtts) 542 543! If we extrapolate the source terms 544if (st_prv_id.ge.0) then 545 ! S as Source, V as Variable 546 do iel = 1, ncel 547 do isou = 1, 6 548 ! Save for exchange 549 tuexpr = c_st_prv(isou,iel) 550 ! For continuation and the next time step 551 c_st_prv(isou,iel) = smbrts(isou,iel) 552 553 smbrts(isou,iel) = - thets*tuexpr 554 ! Right hand side of the previous time step 555 ! We suppose -rovsdt > 0: we implicit 556 ! the user source term (the rest) 557 do jsou = 1, 6 558 smbrts(isou,iel) = smbrts(isou,iel) & 559 + rovsdtts(jsou,isou,iel)*cvara_rij(jsou,iel) 560 ! Diagonal 561 rovsdtts(jsou,isou,iel) = - thetv*rovsdtts(jsou,isou,iel) 562 enddo 563 enddo 564 enddo 565else 566 do iel = 1, ncel 567 do isou = 1, 6 568 do jsou = 1, 6 569 smbrts(isou,iel) = smbrts(isou,iel) & 570 + rovsdtts(jsou,isou,iel)*cvara_rij(jsou,iel) 571 enddo 572 rovsdtts(isou,isou,iel) = max(-rovsdtts(isou,isou,iel), 0.d0) 573 enddo 574 enddo 575endif 576 577! Lagrangian source terms 578!------------------------ 579 580! 2nd order is not taken into account 581if (iilagr.eq.2 .and. ltsdyn.eq.1) then 582 583 call field_get_val_v_by_name('rij_st_lagr', lagr_st_rij) 584 tslagi => tslagr(1:ncelet,itsli) 585 586 do iel = 1,ncel 587 do isou = 1, 6 588 smbrts(isou, iel) = smbrts(isou, iel) + lagr_st_rij(isou,iel) 589 rovsdtts(isou,isou,iel) = rovsdtts(isou,isou, iel) + max(-tslagi(iel),zero) 590 enddo 591 enddo 592 593endif 594 595! Mass source term 596!----------------- 597 598if (ncesmp.gt.0) then 599 600 allocate(gatinj(6,ncelet)) 601 602 ! We increment smbr with -Gamma.var_prev and rovsdr with Gamma 603 call catsmt(ncesmp, 6, icetsm, itypsm(:,ivar), & 604 cell_f_vol, cvara_rij, smacel(:,ivar+isou-1), smacel(:,ipr), & 605 smbrts, rovsdtts, gatinj) 606 607 do isou = 1, 6 608 609 ! If we extrapolate the source terms we put Gamma Pinj in the previous st 610 if (st_prv_id.ge.0) then 611 do iel = 1, ncel 612 c_st_prv(isou,iel) = c_st_prv(isou,iel) + gatinj(isou,iel) 613 enddo 614 ! Otherwise we put it directly in the RHS 615 else 616 do iel = 1, ncel 617 smbrts(isou, iel) = smbrts(isou, iel) + gatinj(isou,iel) 618 enddo 619 endif 620 621 enddo 622 623 deallocate(gatinj) 624 625endif 626 627!=============================================================================== 628! 4.1.2 Unsteady term 629!=============================================================================== 630 631! ---> Added in the matrix diagonal 632 633if (vcopt%istat .eq. 1) then 634 do iel = 1, ncel 635 do isou = 1, 6 636 rovsdtts(isou,isou,iel) = rovsdtts(isou,isou,iel) & 637 + (crom(iel)/dt(iel))*cell_f_vol(iel) 638 enddo 639 enddo 640endif 641 642ivar = irij 643 644!=============================================================================== 645! 4.1.3 Rij-epsilon model-specific terms 646!=============================================================================== 647 648allocate(viscce(6,ncelet)) 649allocate(weighf(2,nfac)) 650allocate(weighb(nfabor)) 651 652! Rij-epsilon standard (LRR) 653if (iturb.eq.30) then 654 655 if (irijco.eq.1) then 656 call resrij2(ivar, & 657 gradv, cpro_produc, gradro, & 658 viscf, viscb, viscce, & 659 smbrts, rovsdtts, & 660 weighf, weighb) 661 else 662 call resrij(ivar, & 663 cpro_produc, gradro, & 664 viscf, viscb, viscce, & 665 smbrts, rovsdtts, & 666 weighf, weighb) 667 endif 668 669! Rij-epsilon SSG or EBRSM 670elseif (iturb.eq.31.or.iturb.eq.32) then 671 672 call resssg2(ivar, & 673 gradv, cpro_produc, gradro, & 674 viscf, viscb, viscce, & 675 smbrts, rovsdtts, & 676 weighf, weighb) 677 678endif 679 680!=============================================================================== 681! 4.2 Solve Rij 682!=============================================================================== 683 684if (st_prv_id.ge.0) then 685 thetp1 = 1.d0 + thets 686 do iel = 1, ncel 687 do isou = 1, 6 688 smbrts(isou,iel) = smbrts(isou,iel) + thetp1*c_st_prv(isou,iel) 689 enddo 690 enddo 691endif 692 693! all boundary convective flux with upwind 694icvflb = 0 695 696call field_get_coefa_v(ivarfl(ivar), coefa_rij) 697call field_get_coefb_v(ivarfl(ivar), coefb_rij) 698call field_get_coefaf_v(ivarfl(ivar), cofaf_rij) 699call field_get_coefbf_v(ivarfl(ivar), cofbf_rij) 700 701vcopt_loc = vcopt 702 703vcopt_loc%istat = -1 704vcopt_loc%idifft = -1 705vcopt_loc%iwgrec = 0 ! Warning, may be overwritten if a field 706vcopt_loc%thetav = thetv 707vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 708 709p_k_value => vcopt_loc 710c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 711 712call cs_equation_iterative_solve_tensor(idtvar, ivarfl(ivar), c_null_char, & 713 c_k_value, cvara_rij, cvara_rij, & 714 coefa_rij, coefb_rij, & 715 cofaf_rij, cofbf_rij, & 716 imasfl, bmasfl, viscf, & 717 viscb, viscf, viscb, viscce, & 718 weighf, weighb, icvflb, ivoid, & 719 rovsdtts, smbrts, cvar_rij) 720 721!=============================================================================== 722! 5. Solve Epsilon 723!=============================================================================== 724 725call reseps(nvar, ncesmp, icetsm, itypsm, & 726 dt, gradv, cpro_produc, gradro, & 727 smacel, viscf, viscb, & 728 smbr, rovsdt) 729 730!=============================================================================== 731! 6. Clipping 732!=============================================================================== 733 734if (iturb.eq.32) then 735 iclip = 1 736else 737 iclip = 2 738endif 739 740if (irijco.eq.1) then 741 call clprij2(ncel) 742else 743 call clprij(ncel, iclip) 744endif 745 746! Free memory 747deallocate(viscf, viscb) 748deallocate(smbr, rovsdt) 749if (allocated(gradro)) deallocate(gradro) 750if (allocated(produc)) deallocate(produc) 751deallocate(gradv) 752 753deallocate(smbrts) 754deallocate(rovsdtts) 755 756!-------- 757! Formats 758!-------- 759 760 1000 format(/, & 761' ** Solving Rij-EPSILON LRR' ,/,& 762' -----------------------' ,/) 763 1001 format(/, & 764' ** Solving Rij-EPSILON SSG' ,/,& 765' -----------------------' ,/) 766 1002 format(/, & 767' ** Solving Rij-EPSILON EBRSM' ,/,& 768' -------------------------' ,/) 769 770!---- 771! End 772!---- 773 774return 775 776end subroutine 777