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 clptrg.f90 28!> 29!> \brief Boundary conditions for rough walls (icodcl = 6). 30!> 31!> The wall functions may change the value of the diffusive flux. 32!> 33!> The values at a boundary face \f$ \fib \f$ stored in the face center 34!> \f$ \centf \f$ of the variable \f$ P \f$ and its diffusive flux \f$ Q \f$ 35!> are written as: 36!> \f[ 37!> P_\centf = A_P^g + B_P^g P_\centi 38!> \f] 39!> and 40!> \f[ 41!> Q_\centf = A_P^f + B_P^f P_\centi 42!> \f] 43!> where \f$ P_\centi \f$ is the value of the variable \f$ P \f$ at the 44!> neighboring cell. 45!> 46!> Warning: 47!> 48!> - for a vector field such as the velocity \f$ \vect{u} \f$ the boundary 49!> conditions may read: 50!> \f[ 51!> \vect{u}_\centf = \vect{A}_u^g + \tens{B}_u^g \vect{u}_\centi 52!> \f] 53!> and 54!> \f[ 55!> \vect{Q}_\centf = \vect{A}_u^f + \tens{B}_u^f \vect{u}_\centi 56!> \f] 57!> where \f$ \tens{B}_u^g \f$ and \f$ \tens{B}_u^f \f$ are 3x3 tensor matrix 58!> which coupled veclocity components next to a boundary. 59!> 60!> Please refer to the <a href="../../theory.pdf#cpltrg"><b>clptrg</b></a> section 61!> of the theory guide for more informations. 62!------------------------------------------------------------------------------- 63 64!------------------------------------------------------------------------------- 65! Arguments 66!______________________________________________________________________________. 67! mode name role ! 68!______________________________________________________________________________! 69!> \param[in] nscal total number of scalars 70!> \param[in] isvhb indicator to save exchange coeffient 71!> \param[in,out] icodcl face boundary condition code: 72!> - 1 Dirichlet 73!> - 3 Neumann 74!> - 4 sliding and 75!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 76!> - 5 smooth wall and 77!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 78!> - 6 rough wall and 79!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 80!> - 9 free inlet/outlet 81!> (input mass flux blocked to 0) 82!> \param[in,out] rcodcl boundary condition values: 83!> - rcodcl(1) value of the dirichlet 84!> - rcodcl(2) value of the exterior exchange 85!> coefficient (infinite if no exchange) 86!> - rcodcl(3) value flux density 87!> (negative if gain) in w/m2 or roughness 88!> in m if icodcl=6 89!> -# for the velocity \f$ (\mu+\mu_T) 90!> \gradv \vect{u} \cdot \vect{n} \f$ 91!> -# for the pressure \f$ \Delta t 92!> \grad P \cdot \vect{n} \f$ 93!> -# for a scalar \f$ cp \left( K + 94!> \dfrac{K_T}{\sigma_T} \right) 95!> \grad T \cdot \vect{n} \f$ 96!> \param[in] velipb value of the velocity at \f$ \centip \f$ 97!> of boundary cells 98!> \param[in] rijipb value of \f$ R_{ij} \f$ at \f$ \centip \f$ 99!> of boundary cells 100!> \param[out] visvdr viscosite dynamique ds les cellules 101!> de bord apres amortisst de v driest 102!> \param[out] hbord coefficients d'echange aux bords 103!> 104!> \param[in] theipb boundary temperature in \f$ \centip \f$ 105!> (more exaclty the energetic variable) 106!_______________________________________________________________________________ 107 108subroutine clptrg & 109 ( nscal , isvhb , icodcl , & 110 rcodcl , & 111 velipb , rijipb , visvdr , & 112 hbord , theipb ) 113 114!=============================================================================== 115 116!=============================================================================== 117! Module files 118!=============================================================================== 119 120use paramx 121use numvar 122use optcal 123use cstphy 124use cstnum 125use pointe 126use entsor 127use albase 128use parall 129use ppppar 130use ppthch 131use ppincl 132use radiat 133use cplsat 134use mesh 135use field 136use lagran 137use turbomachinery 138use cs_c_bindings 139use atincl 140 141!=============================================================================== 142 143implicit none 144 145! Arguments 146 147integer nscal, isvhb 148 149integer, pointer, dimension(:,:) :: icodcl 150 151double precision, pointer, dimension(:,:,:) :: rcodcl 152double precision, dimension(:,:) :: velipb 153double precision, pointer, dimension(:,:) :: rijipb 154double precision, pointer, dimension(:) :: visvdr, hbord, theipb 155 156! Local variables 157 158integer ifac, iel, isou, ii, jj, kk 159integer iscal, clsyme 160integer modntl 161integer iuntur, f_id, iustar 162integer nlogla, nsubla, iuiptn 163integer kdflim 164integer f_id_rough 165 166double precision rnx, rny, rnz 167double precision tx, ty, tz, txn, txn0, t2x, t2y, t2z 168double precision utau, upx, upy, upz, usn 169double precision uiptn, uiptmn, uiptmx 170double precision uetmax, uetmin, ukmax, ukmin, yplumx, yplumn 171double precision tetmax, tetmin, tplumx, tplumn 172double precision dlmomax, dlmomin 173double precision uk, uet, yplus, uplus, phit 174double precision gredu, temp 175double precision cfnns, cfnnk, cfnne 176double precision sqrcmu, ek 177double precision xmutlm 178double precision rcprod, rcflux 179double precision hflui, hint, pimp, qimp 180double precision eloglo(3,3), alpha(6,6) 181double precision rcodcx, rcodcy, rcodcz, rcodcn 182double precision visclc, visctc, romc , distbf, srfbnf 183double precision cofimp 184double precision distb0, rough_d , ydep 185double precision duplus 186double precision dtplus, rough_t, yplus_t 187double precision dsa0 188double precision rinfiv(3) 189double precision visci(3,3), fikis, viscis, distfi 190double precision fcoefa(6), fcoefb(6), fcofaf(6) 191double precision fcofbf(6), fcofad(6), fcofbd(6) 192 193double precision rxx, rxy, rxz, ryy, ryz, rzz, rnnb 194double precision rttb, alpha_rnn, liqwt, totwt 195double precision cpp 196double precision sigmak, sigmae 197double precision coef_mom,coef_momm 198double precision one_minus_ri 199double precision dlmo,dt,theta0,flux 200 201double precision, dimension(:), pointer :: crom 202double precision, dimension(:), pointer :: viscl, visct, cpro_cp, yplbr, ustar 203double precision, dimension(:), allocatable :: byplus, buk 204double precision, dimension(:), allocatable, target :: buet, bcfnns_loc 205double precision, dimension(:), allocatable :: bdlmo 206 207double precision, dimension(:), pointer :: cvar_k, bcfnns 208double precision, dimension(:,:), pointer :: cvar_rij 209double precision, dimension(:), pointer :: cvara_nusa 210 211double precision, dimension(:), pointer :: cvar_totwt, cvar_t, cpro_liqwt 212double precision, dimension(:), pointer :: bpro_rough_d 213double precision, dimension(:), pointer :: bpro_rough_t 214double precision, dimension(:), pointer :: cpro_diff_lim_k 215double precision, dimension(:), pointer :: cpro_diff_lim_eps 216double precision, dimension(:), pointer :: cpro_diff_lim_rij 217 218double precision, dimension(:,:), pointer :: coefau, cofafu, visten 219double precision, dimension(:,:,:), pointer :: coefbu, cofbfu 220double precision, dimension(:), pointer :: coefa_k, coefb_k, coefaf_k, coefbf_k 221double precision, dimension(:), pointer :: coefa_ep, coefaf_ep 222double precision, dimension(:), pointer :: coefb_ep, coefbf_ep 223double precision, dimension(:,:), pointer :: coefa_rij, coefaf_rij, coefad_rij 224double precision, dimension(:,:,:), pointer :: coefb_rij, coefbf_rij, coefbd_rij 225double precision, dimension(:), pointer :: coefa_omg, coefaf_omg 226double precision, dimension(:), pointer :: coefb_omg, coefbf_omg 227double precision, dimension(:), pointer :: coefa_al, coefaf_al 228double precision, dimension(:), pointer :: coefb_al, coefbf_al 229double precision, dimension(:), pointer :: coefa_phi, coefaf_phi 230double precision, dimension(:), pointer :: coefb_phi, coefbf_phi 231double precision, dimension(:), pointer :: coefa_fb, coefaf_fb 232double precision, dimension(:), pointer :: coefb_fb, coefbf_fb 233double precision, dimension(:), pointer :: coefa_nu, coefaf_nu 234double precision, dimension(:), pointer :: coefb_nu, coefbf_nu 235 236integer ntlast , iaff 237data ntlast , iaff /-1 , 0/ 238save ntlast , iaff 239 240type(var_cal_opt) :: vcopt 241type(var_cal_opt) :: vcopt_rij, vcopt_ep 242 243!=============================================================================== 244! Interfaces 245!=============================================================================== 246 247interface 248 249 subroutine clptrg_scalar(iscal, isvhb, icodcl, rcodcl, & 250 byplus , buk, buet , bcfnns, bdlmo, & 251 hbord , theipb , & 252 tetmax , tetmin , tplumx , tplumn ) 253 254 implicit none 255 integer iscal, isvhb 256 integer, pointer, dimension(:,:) :: icodcl 257 double precision, pointer, dimension(:,:,:) :: rcodcl 258 double precision, dimension(:) :: byplus, buk, buet, bcfnns 259 double precision, pointer, dimension(:) :: hbord, theipb 260 double precision tetmax, tetmin, tplumx, tplumn 261 double precision, dimension(:) :: bdlmo 262 263 end subroutine clptrg_scalar 264 265 end interface 266 267!=============================================================================== 268! 1. Initializations 269!=============================================================================== 270 271! Initialize variables to avoid compiler warnings 272 273ek = 0.d0 274phit = 0.d0 275uiptn = 0.d0 276 277rinfiv(1) = rinfin 278rinfiv(2) = rinfin 279rinfiv(3) = rinfin 280 281uet = 1.d0 282utau = 1.d0 283 284! --- Constants 285sqrcmu = sqrt(cmu) 286 287! --- Correction factors for stratification (used in atmospheric models) 288cfnns = 1.d0 289cfnnk = 1.d0 290cfnne = 1.d0 291dlmo = 0.d0 292 293! Alpha constant for a realisable BC for R12 with the SSG model 294alpha_rnn = 0.47d0 295 296! pointers to y+ if saved 297 298yplbr => null() 299 300if (iyplbr.ge.0) then 301 call field_get_val_s(iyplbr, yplbr) 302endif 303if (itytur.eq.3 .and. idirsm.eq.1) call field_get_val_v(ivsten, visten) 304 305! Diffusion limiter 306call field_get_key_id("diffusion_limiter_id", kdflim) 307 308! --- Save wall friction velocity 309 310call field_get_id_try('ustar', iustar) 311if (iustar.ge.0) then !TODO remove, this information is in cofaf cofbf 312 call field_get_val_s(iustar, ustar) 313else 314 allocate(buet(nfabor)) 315 ustar => buet 316endif 317 318! --- Gradient and flux boundary conditions 319 320call field_get_coefa_v(ivarfl(iu), coefau) 321call field_get_coefb_v(ivarfl(iu), coefbu) 322call field_get_coefaf_v(ivarfl(iu), cofafu) 323call field_get_coefbf_v(ivarfl(iu), cofbfu) 324 325if (ik.gt.0) then 326 call field_get_coefa_s(ivarfl(ik), coefa_k) 327 call field_get_coefb_s(ivarfl(ik), coefb_k) 328 call field_get_coefaf_s(ivarfl(ik), coefaf_k) 329 call field_get_coefbf_s(ivarfl(ik), coefbf_k) 330 call field_get_key_double(ivarfl(ik), ksigmas, sigmak) 331 332 ! Diffusion limiter 333 call field_get_key_int(ivarfl(ik), kdflim, f_id) 334 if (f_id.ge.0) then 335 call field_get_val_s(f_id, cpro_diff_lim_k) 336 else 337 cpro_diff_lim_k => null() 338 endif 339else 340 coefa_k => null() 341 coefb_k => null() 342 coefaf_k => null() 343 coefbf_k => null() 344endif 345 346if (iep.gt.0) then 347 call field_get_coefa_s(ivarfl(iep), coefa_ep) 348 call field_get_coefb_s(ivarfl(iep), coefb_ep) 349 call field_get_coefaf_s(ivarfl(iep), coefaf_ep) 350 call field_get_coefbf_s(ivarfl(iep), coefbf_ep) 351 call field_get_key_double(ivarfl(iep), ksigmas, sigmae) 352 ! Diffusion limiter 353 call field_get_key_int(ivarfl(iep), kdflim, f_id) 354 if (f_id.ge.0) then 355 call field_get_val_s(f_id, cpro_diff_lim_eps) 356 else 357 cpro_diff_lim_eps => null() 358 endif 359else 360 coefa_ep => null() 361 coefb_ep => null() 362 coefaf_ep => null() 363 coefbf_ep => null() 364endif 365 366if (itytur.eq.3) then! Also have boundary conditions for the momentum equation 367 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt_rij) 368 call field_get_key_struct_var_cal_opt(ivarfl(iep), vcopt_ep) 369 370 call field_get_coefa_v(ivarfl(irij), coefa_rij) 371 call field_get_coefb_v(ivarfl(irij), coefb_rij) 372 call field_get_coefaf_v(ivarfl(irij), coefaf_rij) 373 call field_get_coefbf_v(ivarfl(irij), coefbf_rij) 374 call field_get_coefad_v(ivarfl(irij), coefad_rij) 375 call field_get_coefbd_v(ivarfl(irij), coefbd_rij) 376 377 ! Diffusion limiter 378 call field_get_key_int(ivarfl(irij), kdflim, f_id) 379 if (f_id.ge.0) then 380 call field_get_val_s(f_id, cpro_diff_lim_rij) 381 else 382 cpro_diff_lim_rij => null() 383 endif 384endif 385 386if (ial.gt.0) then 387 call field_get_coefa_s(ivarfl(ial), coefa_al) 388 call field_get_coefb_s(ivarfl(ial), coefb_al) 389 call field_get_coefaf_s(ivarfl(ial), coefaf_al) 390 call field_get_coefbf_s(ivarfl(ial), coefbf_al) 391else 392 coefa_al => null() 393 coefb_al => null() 394 coefaf_al => null() 395 coefbf_al => null() 396endif 397 398if (iomg.gt.0) then 399 call field_get_coefa_s(ivarfl(iomg), coefa_omg) 400 call field_get_coefb_s(ivarfl(iomg), coefb_omg) 401 call field_get_coefaf_s(ivarfl(iomg), coefaf_omg) 402 call field_get_coefbf_s(ivarfl(iomg), coefbf_omg) 403else 404 coefa_omg => null() 405 coefb_omg => null() 406 coefaf_omg => null() 407 coefbf_omg => null() 408endif 409 410if (iphi.gt.0) then 411 call field_get_coefa_s(ivarfl(iphi), coefa_phi) 412 call field_get_coefb_s(ivarfl(iphi), coefb_phi) 413 call field_get_coefaf_s(ivarfl(iphi), coefaf_phi) 414 call field_get_coefbf_s(ivarfl(iphi), coefbf_phi) 415else 416 coefa_phi => null() 417 coefb_phi => null() 418 coefaf_phi => null() 419 coefbf_phi => null() 420endif 421 422if (ifb.gt.0) then 423 call field_get_coefa_s(ivarfl(ifb), coefa_fb) 424 call field_get_coefb_s(ivarfl(ifb), coefb_fb) 425 call field_get_coefaf_s(ivarfl(ifb), coefaf_fb) 426 call field_get_coefbf_s(ivarfl(ifb), coefbf_fb) 427else 428 coefa_fb => null() 429 coefb_fb => null() 430 coefaf_fb => null() 431 coefbf_fb => null() 432endif 433 434if (inusa.gt.0) then 435 call field_get_val_prev_s(ivarfl(inusa), cvara_nusa) 436 call field_get_coefa_s(ivarfl(inusa), coefa_nu) 437 call field_get_coefaf_s(ivarfl(inusa), coefaf_nu) 438 call field_get_coefb_s(ivarfl(inusa), coefb_nu) 439 call field_get_coefbf_s(ivarfl(inusa), coefbf_nu) 440 call field_get_key_struct_var_cal_opt(ivarfl(inusa), vcopt) 441else 442 cvara_nusa => null() 443 coefa_nu => null() 444 coefb_nu => null() 445 coefaf_nu => null() 446 coefbf_nu => null() 447endif 448 449! --- Physical quantities 450call field_get_val_s(icrom, crom) 451 452if (itytur.eq.2 .or. itytur.eq.5 & 453 .or. iturb.eq.60 .or. iturb.eq.50 .or. iturb.eq.51) then 454 call field_get_val_s(ivarfl(ik), cvar_k) 455endif 456 457if (itytur.eq.3) then 458 call field_get_val_v(ivarfl(irij), cvar_rij) 459endif 460 461call field_get_val_s(iviscl, viscl) 462call field_get_val_s(ivisct, visct) 463if (icp.ge.0) then 464 call field_get_val_s(icp, cpro_cp) 465endif 466 467! min. and max. of wall tangential velocity 468uiptmx = -grand 469uiptmn = grand 470 471! min. and max. of wall friction velocity 472uetmax = -grand 473uetmin = grand 474ukmax = -grand 475ukmin = grand 476 477! min. and max. of y+ 478yplumx = -grand 479yplumn = grand 480 481! min. and max. of wall friction of the thermal scalar 482tetmax = -grand 483tetmin = grand 484 485! min. and max. of inverse of MO length 486dlmomax = -grand 487dlmomin = grand 488 489! min. and max. of T+ 490tplumx = -grand 491tplumn = grand 492 493! Counters (turbulent, laminar, reversal, scale correction) 494nlogla = 0 495nsubla = 0 496iuiptn = 0 497 498 499! With v2f type model, (phi-fbar et BL-v2/k) u=0 is set directly, so 500! uiptmx and uiptmn are necessarily 0 501if (itytur.eq.5) then 502 uiptmx = 0.d0 503 uiptmn = 0.d0 504endif 505 506! Pointers to specific fields 507allocate(byplus(nfabor)) 508allocate(buk(nfabor)) 509allocate(bdlmo(nfabor)) 510 511call field_get_id_try("non_neutral_scalar_correction", f_id) 512if (f_id.ge.0) then 513 call field_get_val_s(f_id, bcfnns) 514else 515 allocate(bcfnns_loc(nfabor)) 516 bcfnns => bcfnns_loc 517endif 518 519cvar_t => null() 520cvar_totwt => null() 521cpro_liqwt => null() 522bpro_rough_d => null() 523bpro_rough_t => null() 524 525call field_get_id_try("boundary_roughness", f_id_rough) 526if (f_id_rough.ge.0) then 527 call field_get_val_s(f_id_rough, bpro_rough_d) 528 529 ! same thermal roughness if not specified 530 call field_get_val_s(f_id_rough, bpro_rough_t) 531endif 532 533call field_get_id_try("boundary_thermal_roughness", f_id_rough) 534if (f_id_rough.ge.0) then 535 call field_get_val_s(f_id_rough, bpro_rough_t) 536endif 537 538if (ippmod(iatmos).ge.1) then 539 theta0 = t0 * (ps / p0)**(rair/cp0) 540 call field_get_val_s(ivarfl(isca(iscalt)), cvar_t) 541 if (ippmod(iatmos).eq.2) then 542 call field_get_val_s(ivarfl(isca(iymw)), cvar_totwt) 543 call field_get_val_s(iliqwt, cpro_liqwt) 544 endif 545endif 546 547! --- Loop on boundary faces 548do ifac = 1, nfabor 549 550 ! Test on the presence of a rough wall 551 if (icodcl(ifac,iu).eq.6) then 552 553 iel = ifabor(ifac) 554 555 ! Physical properties 556 visclc = viscl(iel) 557 visctc = visct(iel) 558 romc = crom(iel) 559 560 ! Geometric quantities 561 distbf = distb(ifac) 562 srfbnf = surfbn(ifac) 563 564 !=========================================================================== 565 ! 1. Local framework 566 !=========================================================================== 567 568 ! Unit normal 569 570 rnx = surfbo(1,ifac)/srfbnf 571 rny = surfbo(2,ifac)/srfbnf 572 rnz = surfbo(3,ifac)/srfbnf 573 574 ! Handle displacement velocity 575 576 rcodcx = rcodcl(ifac,iu,1) 577 rcodcy = rcodcl(ifac,iv,1) 578 rcodcz = rcodcl(ifac,iw,1) 579 580 ! If we are not using ALE, force the displacement velocity for the face 581 ! to be tangential (and update rcodcl for possible use) 582 ! In frozen rotor (iturbo = 1), the velocity is neither tangential to the 583 ! wall (absolute velocity solved in a relative frame of reference) 584 if (iale.eq.0.and.iturbo.eq.0) then 585 rcodcn = rcodcx*rnx+rcodcy*rny+rcodcz*rnz 586 rcodcx = rcodcx -rcodcn*rnx 587 rcodcy = rcodcy -rcodcn*rny 588 rcodcz = rcodcz -rcodcn*rnz 589 rcodcl(ifac,iu,1) = rcodcx 590 rcodcl(ifac,iv,1) = rcodcy 591 rcodcl(ifac,iw,1) = rcodcz 592 endif 593 594 ! Relative tangential velocity 595 596 upx = velipb(ifac,1) - rcodcx 597 upy = velipb(ifac,2) - rcodcy 598 upz = velipb(ifac,3) - rcodcz 599 600 usn = upx*rnx+upy*rny+upz*rnz 601 tx = upx -usn*rnx 602 ty = upy -usn*rny 603 tz = upz -usn*rnz 604 txn = sqrt(tx**2 +ty**2 +tz**2) 605 utau= txn 606 607 ! Unit tangent 608 609 if (txn.ge.epzero) then 610 611 txn0 = 1.d0 612 613 tx = tx/txn 614 ty = ty/txn 615 tz = tz/txn 616 617 else 618 619 ! If the velocity is zero, 620 ! Tx, Ty, Tz is not used (we cancel the velocity), so we assign any 621 ! value (zero for example) 622 623 txn0 = 0.d0 624 625 tx = 0.d0 626 ty = 0.d0 627 tz = 0.d0 628 629 endif 630 631 ! Complete if necessary for Rij-Epsilon 632 633 if (itytur.eq.3) then 634 635 ! --> T2 = RN X T (where X is the cross product) 636 637 t2x = rny*tz - rnz*ty 638 t2y = rnz*tx - rnx*tz 639 t2z = rnx*ty - rny*tx 640 641 ! --> Orthogonal matrix for change of reference frame ELOGLOij 642 ! (from local to global reference frame) 643 644 ! |TX -RNX T2X| 645 ! ELOGLO = |TY -RNY T2Y| 646 ! |TZ -RNZ T2Z| 647 648 ! Its transpose ELOGLOt is its inverse 649 650 eloglo(1,1) = tx 651 eloglo(1,2) = -rnx 652 eloglo(1,3) = t2x 653 eloglo(2,1) = ty 654 eloglo(2,2) = -rny 655 eloglo(2,3) = t2y 656 eloglo(3,1) = tz 657 eloglo(3,2) = -rnz 658 eloglo(3,3) = t2z 659 660 ! Compute Reynolds stress transformation matrix 661 662 clsyme = 0 663 call turbulence_bc_rij_transform(clsyme, eloglo, alpha) 664 665 endif 666 667 !=========================================================================== 668 ! 2. Friction velocities 669 !=========================================================================== 670 671 ! ---> Compute Uet depending if we are in the log zone or not 672 ! in 1 or 2 velocity scales 673 ! and uk based on ek 674 675 676 if (abs(utau).le.epzero) utau = epzero 677 678 ! rough_d: roughness length scale for dynamics 679 rough_d = bpro_rough_d(ifac) 680 681 ! NB: for rough walls, yplus is computed from the roughness and not uk. 682 yplus = distbf/rough_d 683 684 ! Compute turbulent velocity scale 685 if (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) then 686 ek = cvar_k(iel) 687 else if(itytur.eq.3) then 688 ek = 0.5d0*(cvar_rij(1,iel)+cvar_rij(2,iel)+cvar_rij(3,iel)) 689 rxx = cvar_rij(1,iel) 690 rxy = cvar_rij(4,iel) 691 rxz = cvar_rij(6,iel) 692 ryy = cvar_rij(2,iel) 693 ryz = cvar_rij(5,iel) 694 rzz = cvar_rij(3,iel) 695 rnnb = rnx * (rxx * rnx + rxy * rny + rxz * rnz) & 696 + rny * (rxy * rnx + ryy * rny + ryz * rnz) & 697 + rnz * (rxz * rnx + ryz * rny + rzz * rnz) 698 699 rttb = tx * (rxx * tx + rxy * ty + rxz * tz) & 700 + ty * (rxy * tx + ryy * ty + ryz * tz) & 701 + tz * (rxz * tx + ryz * ty + rzz * tz) 702 endif 703 704 ! Neutral value, might be overwritten after 705 uk = cmu025*sqrt(ek) 706 707 ! Pseudo shift of wall by rough_d ((distbf+rough_d)/rough_d) 708 if (iwalfs.ne.3) then 709 ! ustar for neutral, may be modified after 710 uet = utau/log(yplus+1.d0)*xkappa 711 ! Dimensionless velocity, neutral wall function, may be modified after 712 uplus = log(yplus+1.d0)/xkappa 713 714 ! Atmospheric Louis wall functions 715 if (ippmod(iatmos).ge.1) then 716 717 ! Compute reduced gravity for non horizontal walls : 718 gredu = gx*rnx + gy*rny + gz*rnz 719 720 temp = cvar_t(iel) 721 totwt = 0.d0 722 liqwt = 0.d0 723 724 if (ippmod(iatmos).eq.2) then 725 totwt = cvar_totwt(iel) 726 liqwt = cpro_liqwt(iel) 727 endif 728 729 ! 1/U+ for neutral 730 duplus = 1.d0 / uplus 731 732 rough_t = bpro_rough_t(ifac) 733 yplus_t = distbf/rough_t 734 ! 1/T+ 735 dtplus = xkappa/log((distbf+rough_t)/rough_t) 736 737 call atmcls & 738 !========== 739 ( ifac , & 740 utau , rough_d, duplus , dtplus , & 741 yplus_t, & 742 uet , & 743 gredu , & 744 cfnns , cfnnk , cfnne , & 745 dlmo , & 746 temp , totwt , liqwt , & 747 icodcl , rcodcl ) 748 749 endif 750 751 ! Monin Obukhov wall function 752 else 753 754 ! Compute local LMO 755 if (ippmod(iatmos).ge.1) then 756 gredu = gx*rnx + gy*rny + gz*rnz 757 758 if (icodcl(ifac,isca(iscalt)).eq.6) then 759 760 dt = theipb(ifac)-rcodcl(ifac,isca(iscalt),1) 761 call mo_compute_from_thermal_diff(distbf,rough_d,utau,dt, & 762 theta0, gredu, & 763 dlmo, uet) 764 765 elseif (icodcl(ifac,isca(iscalt)).eq.3) then 766 if (icp.ge.0) then 767 cpp = cpro_cp(iel) 768 else 769 cpp = cp0 770 endif 771 772 flux = rcodcl(ifac, isca(iscalt),3)/romc/cpp 773 call mo_compute_from_thermal_flux(distbf,rough_d,utau,flux, & 774 theta0, gredu, & 775 dlmo, uet) 776 777 endif 778 779 else 780 781 ! No temperature delta: neutral 782 call mo_compute_from_thermal_diff(distbf,rough_d,utau,0.d0,0.d0,0.d0, & 783 dlmo,uet) 784 785 endif 786 787 ! Take stability into account for the turbulent velocity scale 788 coef_mom = cs_mo_phim(distbf+rough_d,dlmo) 789 ! Ri = z/L / Phim 790 one_minus_ri = 1.d0-(distbf+rough_d) * dlmo/coef_mom 791 if (one_minus_ri.gt.0) then 792 uk = uk / one_minus_ri**0.25d0 793 794 ! Epsilon should be modified as well to get P+G = P(1-Ri) = epsilon 795 ! P = -R_tn dU/dn = uk^2 uet Phi_m / (kappa z) 796 cfnne = one_minus_ri * coef_mom 797 ! Nothing done for the moment for really high stability 798 else 799 cfnne = 1.d0 800 endif 801 802 endif ! End Monin Obukhov 803 ! Dimensionless velocity, recomputed and therefore may take stability 804 ! into account 805 uplus = utau / uet 806 807 808 ! One velocity scale: set uk to uet 809 if (iwallf.le.2) then 810 uk = uet 811 endif 812 813 uetmax = max(uet,uetmax) 814 uetmin = min(uet,uetmin) 815 ukmax = max(uk,ukmax) 816 ukmin = min(uk,ukmin) 817 yplumx = max(yplus,yplumx) 818 yplumn = min(yplus,yplumn) 819 dlmomin= min(dlmo, dlmomin) 820 dlmomax= max(dlmo, dlmomax) 821 822 ! save turbulent subgrid viscosity after van Driest damping in LES 823 ! care is taken to not dampen it twice at boundary cells having more 824 ! one boundary face 825 if (itytur.eq.4.and.idries.eq.1) then 826 if (visvdr(iel).lt.-900.d0) then 827 ! FIXME amortissement de van Driest a revoir en rugueux : 828 ! visct(iel) = visct(iel)*(1.d0-exp(-yplus/cdries))**2 829 visvdr(iel) = visct(iel) 830 visctc = visct(iel) 831 endif 832 endif 833 834 ! Save yplus if post-processed 835 if (iyplbr.ge.0) then 836 yplbr(ifac) = yplus 837 endif 838 839 !=========================================================================== 840 ! 3. Velocity boundary conditions 841 !=========================================================================== 842 843 ! uiptn respecte la production de k 844 ! de facon conditionnelle --> Coef RCPROD 845 846 ! --> All turbulence models (except v2f and EBRSM) 847 !------------------------------------------------- 848 if (itytur.eq.2 .or. iturb.eq.60 .or. & 849 iturb.eq.0 .or. iturb.eq.10 .or. & 850 iturb.eq.30.or. iturb.eq.31 .or. & 851 itytur.eq.4 .or. & 852 iturb.eq.70 ) then 853 854 if (visctc.gt.epzero) then 855 856 ! Pseudo shift of wall by rough_d ((distbf+rough_d)/rough_d) 857 distb0=distbf+rough_d 858 ! FIXME uk not modified for Louis yet.... 859 xmutlm = xkappa*uk*distb0*romc 860 861 if (iwalfs.ne.3) then 862 rcprod = distbf/distb0*max(1.d0, & 863 2.d0*sqrt(xmutlm/visctc) - distb0/distbf/(2.d0+rough_d/distb0)) 864 865 ! Ground apparent velocity (for log only) 866 uiptn = max(utau - uet/xkappa*rcprod,0.d0) 867 iuntur = 1 868 nlogla = nlogla + 1 869 870 ! Coupled solving of the velocity components 871 ! The boundary term for velocity gradient is implicit 872 ! modified for non neutral boundary layer (in uplus) 873 cofimp = max(1.d0 - 1.d0/(xkappa*uplus)*rcprod, 0.d0) 874 ! The term (rho*uet*uk) is implicit 875 rcflux = max(xmutlm,visctc)/distb0 ! TODO merge with MO without this max 876 hflui = rcflux/(xkappa*uplus) 877 878 !Monin Obukhov 879 else 880 ! Boundary condition on the velocity to have approximately the good 881 ! turbulence production 882 coef_mom = cs_mo_phim(distbf+rough_d,dlmo) 883 coef_momm = cs_mo_phim(2.d0*distbf+rough_d,dlmo) 884 rcprod = 2.d0*distbf*sqrt(xkappa*uk*romc*coef_mom/visctc/distb0) & 885 - coef_momm/(2.d0+rough_d/distbf) 886 887 ! Ground apparent velocity (for log only) 888 uiptn = max(utau - uet/xkappa*rcprod,0.d0) 889 iuntur = 1 890 nlogla = nlogla + 1 891 892 ! Coupled solving of the velocity components 893 ! The boundary term for velocity gradient is implicit 894 ! modified for non neutral boundary layer (in uplus) 895 cofimp = min(max(1.d0 - 1.d0/(xkappa*uplus)*rcprod,0.d0),1.d0) 896 ! The term (rho*uet*uk) is implicit 897 hflui = romc * uk / uplus 898 endif 899 900 ! In the viscous sub-layer 901 else 902 uiptn = 0.d0 903 iuntur = 0 904 nsubla = nsubla + 1 905 906 ! Coupled solving of the velocity components 907 cofimp = 0.d0 908 hflui = visclc / distbf 909 910 endif 911 912 ! Clipping : 913 ! On borne U_f,grad entre 0 et Utau (il y a surement mieux...) 914 ! - 0 : on interdit le retournement en face de bord, qui est en 915 ! contradiction avec l'hypoth\E8se de loi log. 916 ! - Utau : la production turbulente ne peut etre nulle 917 ! On empeche U_f,flux d'etre negatif 918 919 ! --> v2f and EBRSM !FIXME EBRSM 920 !------------------ 921 elseif (itytur.eq.5) then 922 923 ! Avec ces conditions, pas besoin de calculer uiptmx, uiptmn 924 ! et iuiptn qui sont nuls (valeur d'initialisation) 925 iuntur = 0 926 uiptn = 0.d0 927 928 ! Coupled solving of the velocity components 929 hflui = (visclc + visctc) / distbf 930 cofimp = 0.d0 931 932 endif 933 934 ! Min and Max and counter of reversal layer 935 uiptmn = min(uiptn*iuntur,uiptmn) 936 uiptmx = max(uiptn*iuntur,uiptmx) 937 if (uiptn*iuntur.lt.-epzero) iuiptn = iuiptn + 1 938 939 if (itytur.eq.3) then 940 hint = visclc /distbf 941 else 942 hint = (visclc + visctc)/distbf 943 endif 944 945 ! Coupled solving of the velocity components 946 947 ! Gradient boundary conditions 948 !----------------------------- 949 rcodcn = rcodcx*rnx+rcodcy*rny+rcodcz*rnz 950 951 coefau(1,ifac) = (1.d0-cofimp)*(rcodcx - rcodcn*rnx) + rcodcn*rnx 952 coefau(2,ifac) = (1.d0-cofimp)*(rcodcy - rcodcn*rny) + rcodcn*rny 953 coefau(3,ifac) = (1.d0-cofimp)*(rcodcz - rcodcn*rnz) + rcodcn*rnz 954 955 ! Projection in order to have the velocity parallel to the wall 956 ! B = cofimp * ( IDENTITY - n x n ) 957 958 coefbu(1,1,ifac) = cofimp*(1.d0-rnx**2) 959 coefbu(2,2,ifac) = cofimp*(1.d0-rny**2) 960 coefbu(3,3,ifac) = cofimp*(1.d0-rnz**2) 961 coefbu(1,2,ifac) = -cofimp*rnx*rny 962 coefbu(1,3,ifac) = -cofimp*rnx*rnz 963 coefbu(2,1,ifac) = -cofimp*rny*rnx 964 coefbu(2,3,ifac) = -cofimp*rny*rnz 965 coefbu(3,1,ifac) = -cofimp*rnz*rnx 966 coefbu(3,2,ifac) = -cofimp*rnz*rny 967 968 ! Flux boundary conditions 969 !------------------------- 970 971 cofafu(1,ifac) = -hflui*(rcodcx - rcodcn*rnx) - hint*rcodcn*rnx 972 cofafu(2,ifac) = -hflui*(rcodcy - rcodcn*rny) - hint*rcodcn*rny 973 cofafu(3,ifac) = -hflui*(rcodcz - rcodcn*rnz) - hint*rcodcn*rnz 974 975 ! Projection in order to have the shear stress parallel to the wall 976 ! B = hflui*( IDENTITY - n x n ) 977 978 cofbfu(1,1,ifac) = hflui*(1.d0-rnx**2) + hint*rnx**2 979 cofbfu(2,2,ifac) = hflui*(1.d0-rny**2) + hint*rny**2 980 cofbfu(3,3,ifac) = hflui*(1.d0-rnz**2) + hint*rnz**2 981 982 cofbfu(1,2,ifac) = (hint - hflui)*rnx*rny 983 cofbfu(1,3,ifac) = (hint - hflui)*rnx*rnz 984 cofbfu(2,1,ifac) = (hint - hflui)*rny*rnx 985 cofbfu(2,3,ifac) = (hint - hflui)*rny*rnz 986 cofbfu(3,1,ifac) = (hint - hflui)*rnz*rnx 987 cofbfu(3,2,ifac) = (hint - hflui)*rnz*rny 988 989 ! In case of transient turbomachinery computations, save the coefficents 990 ! associated to rough wall velocity BC, in order to update the wall velocity 991 ! after the geometry update (between prediction and correction step) 992 if (iturbo.eq.2) then 993 if (irotce(iel).ne.0) then 994 coftur(ifac) = cofimp 995 hfltur(ifac) = hflui 996 endif 997 endif 998 999 !=========================================================================== 1000 ! 4. Boundary conditions on k and epsilon 1001 !=========================================================================== 1002 1003 ydep = distbf*0.5d0+rough_d 1004 1005 if (itytur.eq.2) then 1006 1007 ! Dirichlet Boundary Condition on k 1008 !---------------------------------- 1009 1010 pimp = uk**2/sqrcmu*cfnnk 1011 hint = (visclc+visctc/sigmak)/distbf 1012 1013 call set_dirichlet_scalar & 1014 !==================== 1015 ( coefa_k(ifac), coefaf_k(ifac), & 1016 coefb_k(ifac), coefbf_k(ifac), & 1017 pimp , hint , rinfin ) 1018 1019 1020 ! No diffusion reconstruction when using wall functions 1021 if (associated(cpro_diff_lim_k)) cpro_diff_lim_k(ifac) = 0.d0 1022 1023 ! Neumann Boundary Condition on epsilon 1024 !-------------------------------------- 1025 1026 hint = (visclc+visctc/sigmae)/distbf 1027 1028 pimp = uk**3/(xkappa*ydep**2)*distbf*cfnne 1029 qimp = -pimp*hint !TODO transform it to use d eps / d y directly 1030 1031 call set_neumann_scalar & 1032 !================== 1033 ( coefa_ep(ifac), coefaf_ep(ifac), & 1034 coefb_ep(ifac), coefbf_ep(ifac), & 1035 qimp , hint ) 1036 1037 ! No diffusion reconstruction when using wall functions 1038 if (associated(cpro_diff_lim_eps)) cpro_diff_lim_eps(ifac) = 0.d0 1039 1040 !=========================================================================== 1041 ! 5. Boundary conditions on Rij-epsilon 1042 !=========================================================================== 1043 1044 elseif (itytur.eq.3) then 1045 1046 ! Exchange coefficient 1047 1048 ! Symmetric tensor diffusivity (Daly Harlow -- GGDH) 1049 if (iand(vcopt_rij%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 1050 1051 visci(1,1) = visclc + visten(1,iel) 1052 visci(2,2) = visclc + visten(2,iel) 1053 visci(3,3) = visclc + visten(3,iel) 1054 visci(1,2) = visten(4,iel) 1055 visci(2,1) = visten(4,iel) 1056 visci(2,3) = visten(5,iel) 1057 visci(3,2) = visten(5,iel) 1058 visci(1,3) = visten(6,iel) 1059 visci(3,1) = visten(6,iel) 1060 1061 ! ||Ki.S||^2 1062 viscis = ( visci(1,1)*surfbo(1,ifac) & 1063 + visci(1,2)*surfbo(2,ifac) & 1064 + visci(1,3)*surfbo(3,ifac))**2 & 1065 + ( visci(2,1)*surfbo(1,ifac) & 1066 + visci(2,2)*surfbo(2,ifac) & 1067 + visci(2,3)*surfbo(3,ifac))**2 & 1068 + ( visci(3,1)*surfbo(1,ifac) & 1069 + visci(3,2)*surfbo(2,ifac) & 1070 + visci(3,3)*surfbo(3,ifac))**2 1071 1072 ! IF.Ki.S 1073 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1074 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1075 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1076 )*surfbo(1,ifac) & 1077 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1078 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1079 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1080 )*surfbo(2,ifac) & 1081 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1082 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1083 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1084 )*surfbo(3,ifac) 1085 1086 distfi = distb(ifac) 1087 1088 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 1089 ! NB: eps =1.d-1 must be consistent with vitens.f90 1090 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1091 1092 hint = viscis/surfbn(ifac)/fikis 1093 1094 ! Scalar diffusivity 1095 else 1096 hint = (visclc+visctc*csrij/cmu)/distbf 1097 endif 1098 1099 ! ---> Tensor Rij (Partially implicited) 1100 1101 do isou = 1, 6 1102 fcoefa(isou) = 0.0d0 1103 fcoefb(isou) = 0.0d0 1104 fcofad(isou) = 0.0d0 1105 fcofbd(isou) = 0.0d0 1106 enddo 1107 1108 do isou = 1, 6 1109 1110 if (isou.eq.1) then 1111 jj = 1 1112 kk = 1 1113 else if (isou.eq.2) then 1114 jj = 2 1115 kk = 2 1116 else if (isou.eq.3) then 1117 jj = 3 1118 kk = 3 1119 else if (isou.eq.4) then 1120 jj = 1 1121 kk = 2 1122 else if (isou.eq.5) then 1123 jj = 2 1124 kk = 3 1125 else if (isou.eq.6) then 1126 jj = 1 1127 kk = 3 1128 endif 1129 1130 ! Coupled version: we implicit as much as possible 1131 if (irijco.eq.1) then 1132 coefa_rij(isou, ifac) = - ( eloglo(jj,1)*eloglo(kk,2) & 1133 + eloglo(jj,2)*eloglo(kk,1)) & 1134 * alpha_rnn * sqrt(rnnb * rttb) * cfnnk 1135 coefaf_rij(isou, ifac) = -hint * coefa_rij(isou, ifac) 1136 coefad_rij(isou, ifac) = 0.d0 1137 do ii = 1, 6 1138 coefb_rij(isou,ii, ifac) = alpha(ii,isou) 1139 if (ii.eq.isou) then 1140 coefbf_rij(isou,ii, ifac) = hint & 1141 * (1.d0 - coefb_rij(isou,ii, ifac)) 1142 else 1143 coefbf_rij(isou,ii, ifac) = - hint & 1144 * coefb_rij(isou,ii, ifac) 1145 endif 1146 coefbd_rij(isou,ii, ifac) = coefb_rij(isou,ii, ifac) 1147 enddo 1148 1149 else if (iclptr.eq.1) then 1150 do ii = 1, 6 1151 if (ii.ne.isou) then 1152 fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii) 1153 endif 1154 enddo 1155 fcoefb(isou) = alpha(isou,isou) 1156 else 1157 do ii = 1, 6 1158 fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii) 1159 enddo 1160 fcoefb(isou) = 0.d0 1161 endif 1162 1163 ! Boundary conditions for the momentum equation 1164 fcofad(isou) = fcoefa(isou) 1165 fcofbd(isou) = fcoefb(isou) 1166 1167 fcoefa(isou) = fcoefa(isou) & 1168 - ( eloglo(jj,1)*eloglo(kk,2) & 1169 + eloglo(jj,2)*eloglo(kk,1))*uet*uk*cfnnk 1170 1171 ! Translate into Diffusive flux BCs 1172 fcofaf(isou) = -hint*fcoefa(isou) 1173 fcofbf(isou) = hint*(1.d0-fcoefb(isou)) 1174 enddo 1175 1176 if (irijco.ne.1) then 1177 do isou = 1, 6 1178 coefa_rij(isou,ifac) = fcoefa(isou) 1179 coefaf_rij(isou,ifac) = fcofaf(isou) 1180 coefad_rij(isou,ifac) = fcofad(isou) 1181 do ii = 1,6 1182 coefb_rij(isou,ii,ifac) = 0 1183 coefbf_rij(isou,ii,ifac) = 0 1184 coefbd_rij(isou,ii,ifac) = 0 1185 enddo 1186 coefb_rij(isou,isou,ifac) = fcoefb(isou) 1187 coefbf_rij(isou,isou,ifac) = fcofbf(isou) 1188 coefbd_rij(isou,isou,ifac) = fcofbd(isou) 1189 enddo 1190 endif 1191 1192 ! No diffusion reconstruction when using wall functions 1193 if (associated(cpro_diff_lim_rij)) cpro_diff_lim_rij(ifac) = 0.d0 1194 1195 ! ---> Epsilon 1196 1197 ! Symmetric tensor diffusivity (Daly Harlow -- GGDH) 1198 if (iand(vcopt_ep%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1199 1200 visci(1,1) = visclc + visten(1,iel)/sigmae 1201 visci(2,2) = visclc + visten(2,iel)/sigmae 1202 visci(3,3) = visclc + visten(3,iel)/sigmae 1203 visci(1,2) = visten(4,iel)/sigmae 1204 visci(2,1) = visten(4,iel)/sigmae 1205 visci(2,3) = visten(5,iel)/sigmae 1206 visci(3,2) = visten(5,iel)/sigmae 1207 visci(1,3) = visten(6,iel)/sigmae 1208 visci(3,1) = visten(6,iel)/sigmae 1209 1210 ! ||Ki.S||^2 1211 viscis = ( visci(1,1)*surfbo(1,ifac) & 1212 + visci(1,2)*surfbo(2,ifac) & 1213 + visci(1,3)*surfbo(3,ifac))**2 & 1214 + ( visci(2,1)*surfbo(1,ifac) & 1215 + visci(2,2)*surfbo(2,ifac) & 1216 + visci(2,3)*surfbo(3,ifac))**2 & 1217 + ( visci(3,1)*surfbo(1,ifac) & 1218 + visci(3,2)*surfbo(2,ifac) & 1219 + visci(3,3)*surfbo(3,ifac))**2 1220 1221 ! IF.Ki.S 1222 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1223 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1224 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1225 )*surfbo(1,ifac) & 1226 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1227 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1228 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1229 )*surfbo(2,ifac) & 1230 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1231 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1232 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1233 )*surfbo(3,ifac) 1234 1235 distfi = distb(ifac) 1236 1237 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 1238 ! NB: eps =1.d-1 must be consistent with vitens.f90 1239 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1240 1241 hint = viscis/surfbn(ifac)/fikis 1242 1243 ! Scalar diffusivity 1244 else 1245 hint = (visclc+visctc/sigmae)/distbf 1246 endif 1247 1248 ! Neumann Boundary Condition on epsilon 1249 !-------------------------------------- 1250 1251 pimp = uk**3/(xkappa*ydep**2)*distbf*cfnne 1252 qimp = -pimp*hint !TODO transform it to use d eps / d y directly 1253 1254 call set_neumann_scalar & 1255 !================== 1256 ( coefa_ep(ifac), coefaf_ep(ifac), & 1257 coefb_ep(ifac), coefbf_ep(ifac), & 1258 qimp , hint ) 1259 1260 ! No diffusion reconstruction when using wall functions 1261 if (associated(cpro_diff_lim_eps)) cpro_diff_lim_eps(ifac) = 0.d0 1262 1263 !=========================================================================== 1264 ! 6a.Boundary conditions on k, epsilon, f_bar and phi in the phi_Fbar model 1265 !=========================================================================== 1266 1267 elseif (iturb.eq.50) then 1268 1269 ! Dirichlet Boundary Condition on k 1270 !---------------------------------- 1271 1272 pimp = 0.d0 1273 hint = (visclc+visctc/sigmak)/distbf 1274 1275 call set_dirichlet_scalar & 1276 !==================== 1277 ( coefa_k(ifac), coefaf_k(ifac), & 1278 coefb_k(ifac), coefbf_k(ifac), & 1279 pimp , hint , rinfin ) 1280 1281 ! Dirichlet Boundary Condition on epsilon 1282 !---------------------------------------- 1283 1284 pimp = 2.0d0*visclc/romc*cvar_k(iel)/distbf**2 1285 hint = (visclc+visctc/sigmae)/distbf 1286 1287 call set_dirichlet_scalar & 1288 !==================== 1289 ( coefa_ep(ifac), coefaf_ep(ifac), & 1290 coefb_ep(ifac), coefbf_ep(ifac), & 1291 pimp , hint , rinfin ) 1292 1293 ! Dirichlet Boundary Condition on Phi 1294 !------------------------------------ 1295 1296 pimp = 0.d0 1297 hint = (visclc+visctc/sigmak)/distbf 1298 1299 call set_dirichlet_scalar & 1300 !==================== 1301 ( coefa_phi(ifac), coefaf_phi(ifac), & 1302 coefb_phi(ifac), coefbf_phi(ifac), & 1303 pimp , hint , rinfin ) 1304 1305 ! Dirichlet Boundary Condition on Fb 1306 !----------------------------------- 1307 1308 pimp = 0.d0 1309 hint = 1.d0/distbf 1310 1311 call set_dirichlet_scalar & 1312 !==================== 1313 ( coefa_fb(ifac), coefaf_fb(ifac), & 1314 coefb_fb(ifac), coefbf_fb(ifac), & 1315 pimp , hint , rinfin ) 1316 1317 !=========================================================================== 1318 ! 6b.Boundary conditions on k, epsilon, phi and alpha in the Bl-v2/k model 1319 !=========================================================================== 1320 1321 elseif (iturb.eq.51) then 1322 1323 ! Dirichlet Boundary Condition on k 1324 !---------------------------------- 1325 1326 pimp = 0.d0 1327 hint = (visclc+visctc/sigmak)/distbf 1328 1329 call set_dirichlet_scalar & 1330 !==================== 1331 ( coefa_k(ifac), coefaf_k(ifac), & 1332 coefb_k(ifac), coefbf_k(ifac), & 1333 pimp , hint , rinfin ) 1334 1335 ! Dirichlet Boundary Condition on epsilon 1336 !---------------------------------------- 1337 1338 pimp = visclc/romc*cvar_k(iel)/distbf**2 1339 hint = (visclc+visctc/sigmae)/distbf 1340 1341 call set_dirichlet_scalar & 1342 !==================== 1343 ( coefa_ep(ifac), coefaf_ep(ifac), & 1344 coefb_ep(ifac), coefbf_ep(ifac), & 1345 pimp , hint , rinfin ) 1346 1347 ! Dirichlet Boundary Condition on Phi 1348 !------------------------------------ 1349 1350 pimp = 0.d0 1351 hint = (visclc+visctc/sigmak)/distbf 1352 1353 call set_dirichlet_scalar & 1354 !==================== 1355 ( coefa_phi(ifac), coefaf_phi(ifac), & 1356 coefb_phi(ifac), coefbf_phi(ifac), & 1357 pimp , hint , rinfin ) 1358 1359 ! Dirichlet Boundary Condition on alpha 1360 !-------------------------------------- 1361 1362 pimp = 0.d0 1363 hint = 1.d0/distbf 1364 1365 call set_dirichlet_scalar & 1366 !==================== 1367 ( coefa_al(ifac), coefaf_al(ifac), & 1368 coefb_al(ifac), coefbf_al(ifac), & 1369 pimp , hint , rinfin ) 1370 1371 1372 !=========================================================================== 1373 ! 7. Boundary conditions on k and omega 1374 !=========================================================================== 1375 1376 elseif (iturb.eq.60) then 1377 1378 ! Always out of the viscous sub-layer 1379 1380 ! Dirichlet Boundary Condition on k 1381 !---------------------------------- 1382 1383 pimp = uk**2/sqrcmu 1384 1385 !FIXME it is wrong because sigma is computed within the model 1386 ! see turbkw.f90 1387 hint = (visclc+visctc/ckwsk2)/distbf 1388 1389 call set_dirichlet_scalar & 1390 !==================== 1391 ( coefa_k(ifac), coefaf_k(ifac), & 1392 coefb_k(ifac), coefbf_k(ifac), & 1393 pimp , hint , rinfin ) 1394 1395 ! Neumann Boundary Condition on omega 1396 !------------------------------------ 1397 1398 !FIXME it is wrong because sigma is computed within the model 1399 ! see turbkw.f90 (So the flux is not the one we impose!) 1400 hint = (visclc+visctc/ckwsw2)/distbf 1401 1402 pimp = distbf*4.d0*uk**3*romc**2/ & 1403 (sqrcmu*xkappa*visclc**2*yplus**2) 1404 qimp = -pimp*hint !TODO transform it to use d eps / d y directly 1405 1406 call set_neumann_scalar & 1407 !================== 1408 ( coefa_omg(ifac), coefaf_omg(ifac), & 1409 coefb_omg(ifac), coefbf_omg(ifac), & 1410 qimp , hint ) 1411 1412 !=========================================================================== 1413 ! 7.1 Boundary conditions on the Spalart Allmaras turbulence model 1414 !=========================================================================== 1415 1416 elseif (iturb.eq.70) then 1417 1418 dsa0 = rough_d ! FIXME is it the sand grain roughness or the length scale as here? 1419 hint = (visclc + vcopt%idifft*cvara_nusa(iel)*romc*dsa0/(distbf+dsa0) ) & 1420 / distbf / csasig 1421 1422 ! If we have a rough wall then: 1423 ! nusa_wall*(1- I'F/d0)=nusa_I' 1424 ! which is a Robin type BC 1425 coefa_nu(ifac) = 0.d0 1426 coefb_nu(ifac) = dsa0/(dsa0+distbf) 1427 1428 1429 coefaf_nu(ifac) = 0.d0 1430 coefbf_nu(ifac) = hint*distbf/(dsa0+distbf) 1431 endif 1432 1433 byplus(ifac) = yplus 1434 buk(ifac) = uk 1435 ustar(ifac) = uet 1436 bcfnns(ifac) = cfnns 1437 bdlmo(ifac) = dlmo 1438 1439 endif 1440 ! Test on the presence of a rough wall (End) 1441 1442enddo 1443! --- End of loop over faces 1444 1445!=========================================================================== 1446! 8. Boundary conditions on the other scalars 1447! (Specific treatment for the variances of the scalars next to walls: 1448! see condli) 1449!=========================================================================== 1450 1451do iscal = 1, nscal 1452 1453 if (iscavr(iscal).le.0) then 1454 1455 call clptrg_scalar(iscal, isvhb, icodcl, rcodcl, & 1456 byplus, buk, ustar, bcfnns, bdlmo, & 1457 hbord, theipb, & 1458 tetmax, tetmin, tplumx, tplumn) 1459 endif 1460 1461enddo 1462 1463if (irangp.ge.0) then 1464 call parmin (uiptmn) 1465 call parmax (uiptmx) 1466 call parmin (uetmin) 1467 call parmax (uetmax) 1468 call parmin (ukmin) 1469 call parmax (ukmax) 1470 call parmin (yplumn) 1471 call parmax (yplumx) 1472 call parcpt (nlogla) 1473 call parcpt (nsubla) 1474 call parcpt (iuiptn) 1475 if (iscalt.gt.0) then 1476 call parmin (tetmin) 1477 call parmax (tetmax) 1478 call parmin (tplumn) 1479 call parmax (tplumx) 1480 endif 1481endif 1482 1483deallocate(byplus) 1484deallocate(buk) 1485if (allocated(buet)) deallocate(buet) 1486if (allocated(bcfnns_loc)) deallocate(bcfnns_loc) 1487deallocate(bdlmo) 1488 1489!=============================================================================== 1490! 9. Writings 1491!=============================================================================== 1492 1493! Remarque : afin de ne pas surcharger les logs dans le cas ou 1494! quelques yplus ne sont pas corrects, on ne produit le message 1495! qu'aux deux premiers pas de temps ou le message apparait et 1496! aux deux derniers pas de temps du calcul, ou si IWARNI est >= 2 1497! On indique aussi le numero du dernier pas de temps auquel on 1498! a rencontre des yplus hors bornes admissibles 1499 1500call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 1501 1502if (vcopt%iwarni.ge.0) then 1503 if (ntlist.gt.0) then 1504 modntl = mod(ntcabs,ntlist) 1505 elseif (ntlist.eq.-1.and.ntcabs.eq.ntmabs) then 1506 modntl = 0 1507 else 1508 modntl = 1 1509 endif 1510 1511 if ((modntl.eq.0 .or. vcopt%iwarni.ge.2).and.iscalt.gt.0 & 1512 .and.ippmod(iatmos).ge.1) then 1513 write(nfecra,2012) & 1514 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1515 tetmin, tetmax, tplumn, tplumx, dlmomin, dlmomax, & 1516 iuiptn,nsubla,nsubla+nlogla 1517 else if ((modntl.eq.0 .or. vcopt%iwarni.ge.2).and.iscalt.gt.0) then 1518 write(nfecra,2011) & 1519 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1520 tetmin, tetmax, tplumn, tplumx, iuiptn,nsubla,nsubla+nlogla 1521 elseif (modntl.eq.0 .or. vcopt%iwarni.ge.2) then 1522 write(nfecra,2010) & 1523 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1524 iuiptn, nsubla,nsubla+nlogla 1525 endif 1526 1527endif 1528 1529!=============================================================================== 1530! 10. Formats 1531!=============================================================================== 1532 1533 2010 format(/, & 1534 3X,'** BOUNDARY CONDITIONS FOR ROUGH WALLS',/, & 1535 ' --------------------------------------',/, & 1536 '------------------------------------------------------------',/,& 1537 ' Minimum Maximum',/,& 1538 '------------------------------------------------------------',/,& 1539 ' Rel velocity at the wall uiptn : ',2E12.5 ,/,& 1540 ' Friction velocity uet : ',2E12.5 ,/,& 1541 ' Friction velocity uk : ',2E12.5 ,/,& 1542 ' Rough dimensionless dist yplus : ',2E12.5 ,/,& 1543 ' ------------------------------------------------------' ,/,& 1544 ' Nb of reversal of the velocity at the wall : ',I10 ,/,& 1545 ' Nb of faces within the viscous sub-layer : ',I10 ,/,& 1546 ' Total number of wall faces : ',I10 ,/,& 1547 '------------------------------------------------------------', & 1548 /,/) 1549 1550 2011 format(/, & 1551 3X,'** BOUNDARY CONDITIONS FOR ROUGH WALLS',/, & 1552 ' --------------------------------------',/, & 1553 '------------------------------------------------------------',/,& 1554 ' Minimum Maximum',/,& 1555 '------------------------------------------------------------',/,& 1556 ' Rel velocity at the wall uiptn : ',2E12.5 ,/,& 1557 ' Friction velocity uet : ',2E12.5 ,/,& 1558 ' Friction velocity uk : ',2E12.5 ,/,& 1559 ' Rough dimensionless dist yplus : ',2E12.5 ,/,& 1560 ' Friction thermal sca. tstar : ',2E12.5 ,/,& 1561 ' Rough dim-less th. sca. tplus : ',2E12.5 ,/,& 1562 ' ------------------------------------------------------' ,/,& 1563 ' Nb of reversal of the velocity at the wall : ',I10 ,/,& 1564 ' Nb of faces within the viscous sub-layer : ',I10 ,/,& 1565 ' Total number of wall faces : ',I10 ,/,& 1566 '------------------------------------------------------------', & 1567 /,/) 1568 1569 2012 format(/, & 1570 3X,'** BOUNDARY CONDITIONS FOR ROUGH WALLS',/, & 1571 ' --------------------------------------',/, & 1572 '------------------------------------------------------------',/,& 1573 ' Minimum Maximum',/,& 1574 '------------------------------------------------------------',/,& 1575 ' Rel velocity at the wall uiptn : ',2E12.5 ,/,& 1576 ' Friction velocity uet : ',2E12.5 ,/,& 1577 ' Friction velocity uk : ',2E12.5 ,/,& 1578 ' Rough dimensionless dist yplus : ',2E12.5 ,/,& 1579 ' Friction thermal sca. tstar : ',2E12.5 ,/,& 1580 ' Rough dim-less th. sca. tplus : ',2E12.5 ,/,& 1581 ' Inverse Monin-Ob. length dlmo : ',2E12.5 ,/,& 1582 ' ------------------------------------------------------' ,/,& 1583 ' Nb of reversal of the velocity at the wall : ',I10 ,/,& 1584 ' Nb of faces within the viscous sub-layer : ',I10 ,/,& 1585 ' Total number of wall faces : ',I10 ,/,& 1586 '------------------------------------------------------------', & 1587 /,/) 1588 1589!---- 1590! End 1591!---- 1592 1593return 1594end subroutine 1595 1596!=============================================================================== 1597! Local functions 1598!=============================================================================== 1599 1600!------------------------------------------------------------------------------- 1601! Arguments 1602!______________________________________________________________________________. 1603! mode name role ! 1604!______________________________________________________________________________! 1605!> \param[in] iscal scalar id 1606!> \param[in] isvhb indicator to save exchange coeffient 1607!> \param[in,out] icodcl face boundary condition code: 1608!> - 1 Dirichlet 1609!> - 3 Neumann 1610!> - 4 sliding and 1611!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 1612!> - 5 smooth wall and 1613!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 1614!> - 6 rough wall and 1615!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 1616!> - 9 free inlet/outlet 1617!> (input mass flux blocked to 0) 1618!> \param[in,out] rcodcl boundary condition values: 1619!> - rcodcl(1) value of the dirichlet 1620!> - rcodcl(2) value of the exterior exchange 1621!> coefficient (infinite if no exchange) 1622!> - rcodcl(3) value flux density 1623!> (negative if gain) in w/m2 1624!> -# for the velocity \f$ (\mu+\mu_T) 1625!> \gradv \vect{u} \cdot \vect{n} \f$ 1626!> -# for the pressure \f$ \Delta t 1627!> \grad P \cdot \vect{n} \f$ 1628!> -# for a scalar \f$ cp \left( K + 1629!> \dfrac{K_T}{\sigma_T} \right) 1630!> \grad T \cdot \vect{n} \f$ 1631!> \param[in] byplus dimensionless distance to the wall 1632!> \param[in] buk dimensionless velocity 1633!> \param[in] buet boundary ustar value 1634!> \param[in] bcfnns boundary correction factor 1635!> \param[in] bdlmo boundary Monin Obukhov length inverse 1636!> \param[in,out] hbord exchange coefficient at boundary 1637!> \param[in] theipb boundary temperature in \f$ \centip \f$ 1638!> (more exaclty the energetic variable) 1639!> \param[out] tetmax maximum local ustar value 1640!> \param[out] tetmin minimum local ustar value 1641!> \param[out] tplumx maximum local tplus value 1642!> \param[out] tplumn minimum local tplus value 1643!_______________________________________________________________________________ 1644 1645subroutine clptrg_scalar & 1646 ( iscal , isvhb , icodcl , & 1647 rcodcl , & 1648 byplus , buk , buet , bcfnns , & 1649 bdlmo , & 1650 hbord , theipb , & 1651 tetmax , tetmin , tplumx , tplumn) 1652 1653!=============================================================================== 1654! Module files 1655!=============================================================================== 1656 1657use paramx 1658use numvar 1659use optcal 1660use cstphy 1661use cstnum 1662use pointe 1663use entsor 1664use albase 1665use parall 1666use ppppar 1667use ppthch 1668use ppincl 1669use radiat 1670use cplsat 1671use mesh 1672use field 1673use lagran 1674use turbomachinery 1675use cs_c_bindings 1676use field_operator 1677use atincl 1678 1679!=============================================================================== 1680 1681implicit none 1682 1683! Arguments 1684 1685integer iscal, isvhb 1686integer, pointer, dimension(:,:) :: icodcl 1687double precision, pointer, dimension(:,:,:) :: rcodcl 1688double precision, dimension(:) :: byplus, buk, buet, bcfnns 1689double precision, pointer, dimension(:) :: hbord, theipb 1690double precision tetmax, tetmin, tplumx, tplumn 1691double precision, dimension(:) :: bdlmo 1692 1693! Local variables 1694 1695integer ivar, f_id, b_f_id, isvhbl 1696integer f_id_ut 1697integer ifac, iel, isou, jsou 1698integer iscacp, ifcvsl, itplus, itstar 1699integer f_id_rough 1700integer kturt, turb_flux_model, turb_flux_model_type 1701 1702double precision cpp, rkl, prdtl, visclc, romc, tplus, cpscv 1703double precision distfi, distbf, fikis, hint, heq, hflui, hext 1704double precision yplus, phit, pimp, temp, tet, uk 1705double precision viscis, visctc, cofimp 1706double precision dtplus, rough_t, visls_0 1707double precision rinfiv(3), pimpv(3) 1708double precision visci(3,3), hintt(6) 1709double precision turb_schmidt, exchange_coef 1710double precision rcprod 1711double precision coef_mom,coef_moh,coef_mohh,dlmo 1712 1713character(len=80) :: fname 1714 1715double precision, dimension(:), pointer :: val_s, bval_s, crom, viscls 1716double precision, dimension(:), pointer :: viscl, visct, cpro_cp, cpro_cv 1717double precision, dimension(:), pointer :: bpro_rough_t 1718 1719double precision, dimension(:), pointer :: bfconv, bhconv 1720double precision, dimension(:), pointer :: tplusp, tstarp 1721double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp 1722double precision, dimension(:,:), pointer :: coefaut, cofafut, cofarut, visten 1723double precision, dimension(:,:,:), pointer :: coefbut, cofbfut, cofbrut 1724 1725integer, save :: kbfid = -1 1726 1727type(var_cal_opt) :: vcopt 1728 1729!=============================================================================== 1730 1731ivar = isca(iscal) 1732f_id = ivarfl(ivar) 1733 1734call field_get_val_s(ivarfl(ivar), val_s) 1735 1736call field_get_val_s(iviscl, viscl) 1737call field_get_val_s(ivisct, visct) 1738 1739call field_get_key_int (f_id, kivisl, ifcvsl) 1740 1741if (ifcvsl .ge. 0) then 1742 call field_get_val_s(ifcvsl, viscls) 1743endif 1744 1745call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1746 1747! If we have no diffusion, no boundary face should have a wall BC type 1748! (this is ensured in typecl) 1749 1750if (vcopt%idiff .eq. 0) then 1751 tetmax = 0.d0 1752 tetmin = 0.d0 1753 tplumx = 0.d0 1754 tplumn = 0.d0 1755 return 1756endif 1757 1758! Get the turbulent flux model for the scalar 1759call field_get_key_id('turbulent_flux_model', kturt) 1760call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 1761turb_flux_model_type = turb_flux_model / 10 1762 1763if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0.or.turb_flux_model_type.eq.3) then 1764 if (iturb.ne.32.or.turb_flux_model_type.eq.3) then 1765 call field_get_val_v(ivsten, visten) 1766 else ! EBRSM and (GGDH or AFM) 1767 call field_get_val_v(ivstes, visten) 1768 endif 1769endif 1770 1771call field_get_coefa_s(f_id, coefap) 1772call field_get_coefb_s(f_id, coefbp) 1773call field_get_coefaf_s(f_id, cofafp) 1774call field_get_coefbf_s(f_id, cofbfp) 1775 1776call field_get_val_s(icrom, crom) 1777if (icp.ge.0) then 1778 call field_get_val_s(icp, cpro_cp) 1779endif 1780 1781if (ippmod(icompf) .ge. 0) then 1782 if (icv.ge.0) then 1783 call field_get_val_s(icv, cpro_cv) 1784 endif 1785endif 1786 1787isvhbl = 0 1788if (iscal.eq.isvhb) then 1789 isvhbl = isvhb 1790endif 1791 1792if (iscal.eq.iscalt) then 1793 ! min. and max. of wall friction of the thermal scalar 1794 tetmax = -grand 1795 tetmin = grand 1796 ! min. and max. of T+ 1797 tplumx = -grand 1798 tplumn = grand 1799endif 1800 1801rinfiv(1) = rinfin 1802rinfiv(2) = rinfin 1803rinfiv(3) = rinfin 1804 1805if (turb_flux_model_type.eq.3) then 1806 1807 ! Name of the scalar ivar 1808 call field_get_name(ivarfl(ivar), fname) 1809 1810 ! Index of the corresponding turbulent flux 1811 call field_get_id(trim(fname)//'_turbulent_flux', f_id_ut) 1812 1813 call field_get_coefa_v(f_id_ut,coefaut) 1814 call field_get_coefb_v(f_id_ut,coefbut) 1815 call field_get_coefaf_v(f_id_ut,cofafut) 1816 call field_get_coefbf_v(f_id_ut,cofbfut) 1817 call field_get_coefad_v(f_id_ut,cofarut) 1818 call field_get_coefbd_v(f_id_ut,cofbrut) 1819 1820endif 1821 1822! pointers to T+ and T* if saved 1823 1824itplus = -1 1825itstar = -1 1826 1827tplusp => null() 1828tstarp => null() 1829 1830if (iscal.eq.iscalt) then 1831 call field_get_id_try('tplus', itplus) 1832 if (itplus.ge.0) then 1833 call field_get_val_s (itplus, tplusp) 1834 endif 1835 call field_get_id_try('tstar', itstar) 1836 if (itstar.ge.0) then 1837 call field_get_val_s (itstar, tstarp) 1838 endif 1839endif 1840 1841bpro_rough_t => null() 1842 1843call field_get_id_try("boundary_roughness", f_id_rough) 1844if (f_id_rough.ge.0) then 1845 ! same thermal roughness if not specified 1846 call field_get_val_s(f_id_rough, bpro_rough_t) 1847endif 1848 1849call field_get_id_try("boundary_thermal_roughness", f_id_rough) 1850if (f_id_rough.ge.0) then 1851 call field_get_val_s(f_id_rough, bpro_rough_t) 1852endif 1853 1854 1855! Pointers to specific fields 1856 1857if (iirayo.ge.1 .and. iscal.eq.iscalt) then 1858 call field_get_val_s_by_name("rad_convective_flux", bfconv) 1859 call field_get_val_s_by_name("rad_exchange_coefficient", bhconv) 1860endif 1861 1862if (kbfid.lt.0) call field_get_key_id("boundary_value_id", kbfid) 1863 1864call field_get_key_int(f_id, kbfid, b_f_id) 1865 1866! If thermal variable has no boundary but temperature does, use it 1867if (b_f_id .lt. 0 .and. iscal.eq.iscalt .and. itherm.eq.2) then 1868 b_f_id = itempb 1869endif 1870 1871if (b_f_id .ge. 0) then 1872 call field_get_val_s(b_f_id, bval_s) 1873else 1874 bval_s => null() 1875endif 1876 1877! Does the scalar behave as a temperature ? 1878call field_get_key_int(f_id, kscacp, iscacp) 1879 1880! Retrieve turbulent Schmidt value for current scalar 1881call field_get_key_double(f_id, ksigmas, turb_schmidt) 1882 1883! Reference diffusivity 1884call field_get_key_double(f_id, kvisl0, visls_0) 1885 1886! --- Loop on boundary faces 1887do ifac = 1, nfabor 1888 1889 yplus = byplus(ifac) 1890 uk = buk(ifac) 1891 dlmo = bdlmo(ifac) 1892 1893 ! Test on the presence of a rough wall condition (start) 1894 if (icodcl(ifac,iu).eq.6) then 1895 1896 iel = ifabor(ifac) 1897 1898 ! Physical quantities 1899 1900 visclc = viscl(iel) 1901 visctc = visct(iel) 1902 romc = crom(iel) 1903 1904 ! Geometric quantities 1905 distbf = distb(ifac) 1906 1907 cpp = 1.d0 1908 if (iscacp.eq.1) then 1909 if (icp.ge.0) then 1910 cpp = cpro_cp(iel) 1911 else 1912 cpp = cp0 1913 endif 1914 endif 1915 1916 if (ifcvsl.lt.0) then 1917 rkl = visls_0 1918 prdtl = cpp*visclc/rkl 1919 else 1920 rkl = viscls(iel) 1921 prdtl = cpp*visclc/rkl 1922 endif 1923 1924 ! --> Compressible module: 1925 ! On suppose que le nombre de Pr doit etre 1926 ! defini de la meme facon que l'on resolve 1927 ! en enthalpie ou en energie, soit Mu*Cp/Lambda. 1928 ! Si l'on resout en energie, on a calcule ci-dessus 1929 ! Mu*Cv/Lambda. 1930 1931 if (iscal.eq.iscalt .and. itherm.eq.3) then 1932 if (icp.ge.0) then 1933 prdtl = prdtl*cpro_cp(iel) 1934 else 1935 prdtl = prdtl*cp0 1936 endif 1937 if (icv.ge.0) then 1938 prdtl = prdtl/cpro_cv(iel) 1939 else 1940 prdtl = prdtl/cv0 1941 endif 1942 endif 1943 1944 ! Scalar diffusivity 1945 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 1946 ! En compressible, pour l'energie LAMBDA/CV+CP/CV*(MUT/TURB_SCHMIDT) 1947 if (iscal.eq.iscalt .and. itherm.eq.3) then 1948 if (icp.ge.0) then 1949 cpscv = cpro_cp(iel) 1950 else 1951 cpscv = cp0 1952 endif 1953 if (icv.ge.0) then 1954 cpscv = cpscv/cpro_cv(iel) 1955 else 1956 cpscv = cpscv/cv0 1957 endif 1958 hint = (rkl+vcopt%idifft*cpscv*visctc/turb_schmidt)/distbf 1959 else 1960 hint = (rkl+vcopt%idifft*cpp*visctc/turb_schmidt)/distbf 1961 endif 1962 1963 ! Symmetric tensor diffusivity (GGDH or AFM) 1964 else if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1965 ! En compressible, pour l'energie LAMBDA/CV+CP/CV*(MUT/SIGMAS) 1966 if (iscal.eq.iscalt .and. itherm.eq.3) then 1967 if (icp.ge.0) then 1968 cpscv = cpro_cp(iel) 1969 else 1970 cpscv = cp0 1971 endif 1972 if (icv.ge.0) then 1973 cpscv = cpscv/cpro_cv(iel) 1974 else 1975 cpscv = cpscv/cv0 1976 endif 1977 temp = vcopt%idifft*cpscv*ctheta(iscal)/csrij 1978 else 1979 temp = vcopt%idifft*cpp*ctheta(iscal)/csrij 1980 endif 1981 visci(1,1) = temp*visten(1,iel) + rkl 1982 visci(2,2) = temp*visten(2,iel) + rkl 1983 visci(3,3) = temp*visten(3,iel) + rkl 1984 visci(1,2) = temp*visten(4,iel) 1985 visci(2,1) = temp*visten(4,iel) 1986 visci(2,3) = temp*visten(5,iel) 1987 visci(3,2) = temp*visten(5,iel) 1988 visci(1,3) = temp*visten(6,iel) 1989 visci(3,1) = temp*visten(6,iel) 1990 1991 ! ||Ki.S||^2 1992 viscis = ( visci(1,1)*surfbo(1,ifac) & 1993 + visci(1,2)*surfbo(2,ifac) & 1994 + visci(1,3)*surfbo(3,ifac))**2 & 1995 + ( visci(2,1)*surfbo(1,ifac) & 1996 + visci(2,2)*surfbo(2,ifac) & 1997 + visci(2,3)*surfbo(3,ifac))**2 & 1998 + ( visci(3,1)*surfbo(1,ifac) & 1999 + visci(3,2)*surfbo(2,ifac) & 2000 + visci(3,3)*surfbo(3,ifac))**2 2001 2002 ! IF.Ki.S 2003 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 2004 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 2005 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 2006 )*surfbo(1,ifac) & 2007 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 2008 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 2009 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 2010 )*surfbo(2,ifac) & 2011 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 2012 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 2013 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 2014 )*surfbo(3,ifac) 2015 2016 distfi = distb(ifac) 2017 2018 ! Take I" so that I"F= eps*||FI||*Ki.n when I" is not in cell i 2019 ! NB: eps =1.d-1 must be consistent with vitens.f90 2020 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 2021 2022 hint = viscis/surfbn(ifac)/fikis 2023 endif 2024 2025 ! Note: for Neumann, Tplus is chosen for post-processing 2026 rough_t = bpro_rough_t(ifac) 2027 2028 ! Modified wall function from Louis 2029 if (iwalfs.ne.3) then 2030 2031 ! T+ = (T_I - T_w) / Tet 2032 tplus = log((distbf+rough_t)/rough_t)/ (xkappa * bcfnns(ifac)) 2033 else 2034 ! Dry atmosphere, Monin Obukhov 2035 coef_moh = cs_mo_psih(distbf+rough_t,rough_t,dlmo) 2036 ! T+ 2037 tplus = coef_moh / xkappa 2038 endif 2039 2040 ! Dirichlet on the scalar, with wall function 2041 if (iturb.ne.0.and.icodcl(ifac,ivar).eq.6) then 2042 ! 1/T+ 2043 dtplus = 1.d0 / tplus 2044 !FIXME apparently buet should be buk 2045 hflui = romc*cpp*buet(ifac) * dtplus 2046 2047 ! Neumann on the scalar, with wall function (for post-processing) 2048 else 2049 hflui = hint 2050 endif 2051 2052 hext = rcodcl(ifac,ivar,2) 2053 pimp = rcodcl(ifac,ivar,1) 2054 2055 if (abs(hext).gt.rinfin*0.5d0) then 2056 heq = hflui 2057 else 2058 heq = hflui*hext/(hflui+hext) 2059 endif 2060 2061 ! Dirichlet BC with wall function 2062 ! Gradients and flux BCs are imposed 2063 if (icodcl(ifac,ivar).eq.6) then 2064 !FIXME this should also be done for Neumann, but overwritten in condli for now 2065 ! Same remark for smooth wall... 2066 !if ((icodcl(ifac,ivar).eq.6).or.(icodcl(ifac,ivar).eq.3)) then 2067 2068 ! Modified wall function from Louis 2069 if (iwalfs.ne.3) then 2070 cofimp = 1.d0 - heq/hint 2071 2072 ! Monin obukhov 2073 else 2074 2075 ! To approximately respect thermal turbulent production with 2 hypothesis 2076 coef_mom = cs_mo_phim(distbf+rough_t,dlmo) 2077 coef_mohh = cs_mo_phih (2.0*distbf+rough_t,dlmo) 2078 ! Gradient BCs 2079 coefap(ifac) = 0.d0 2080 2081 rcprod = 2.d0*romc/visctc*distbf*uk*tplus/coef_mom & 2082 - coef_mohh/(2.d0+rough_t/distbf) 2083 2084 cofimp = 1.d0 - rcprod / (xkappa*tplus) 2085 endif 2086 2087 ! To be coherent with a wall function, clip it to 0 2088 cofimp = max(cofimp, 0.d0) 2089 2090 ! Gradient BCs 2091 coefap(ifac) = (1.d0 - cofimp)*pimp 2092 coefbp(ifac) = cofimp 2093 2094 ! Flux BCs 2095 cofafp(ifac) = -heq*pimp 2096 cofbfp(ifac) = heq 2097 2098 ! Storage of the thermal exchange coefficient 2099 ! (conversion in case of energy or enthalpy) 2100 ! the exchange coefficient is in W/(m2 K) 2101 ! Useful for thermal coupling or radiative transfer 2102 if (iirayo.ge.1 .and. iscal.eq.iscalt.or.isvhbl.gt.0) then 2103 ! Enthalpy 2104 if (itherm.eq.2) then 2105 ! If Cp is variable 2106 if (icp.ge.0) then 2107 exchange_coef = hflui*cpro_cp(iel) 2108 else 2109 exchange_coef = hflui*cp0 2110 endif 2111 2112 ! Total energy (compressible module) 2113 elseif (itherm.eq.3) then 2114 ! If Cv is variable 2115 if (icv.ge.0) then 2116 exchange_coef = hflui*cpro_cv(iel) 2117 else 2118 exchange_coef = hflui*cv0 2119 endif 2120 2121 ! Temperature 2122 elseif (iscacp.eq.1) then 2123 exchange_coef = hflui 2124 endif 2125 endif 2126 2127 ! ---> Thermal coupling, store h = lambda/d 2128 if (isvhbl.gt.0) hbord(ifac) = exchange_coef 2129 2130 ! ---> Radiative transfer 2131 if (iirayo.ge.1 .and. iscal.eq.iscalt) then 2132 bhconv(ifac) = exchange_coef 2133 2134 ! The outgoing flux is stored (Q = h(Ti'-Tp): negative if 2135 ! gain for the fluid) in W/m2 2136 bfconv(ifac) = cofafp(ifac) + cofbfp(ifac)*theipb(ifac) 2137 endif 2138 2139 endif ! End if icodcl=6 2140 2141 !--> Turbulent heat flux 2142 if (turb_flux_model_type.eq.3) then 2143 2144 phit = (cofafp(ifac) + cofbfp(ifac)*val_s(iel)) 2145 2146 hintt(1) = 0.5d0*(visclc+rkl)/distbf & 2147 + visten(1,iel)*ctheta(iscal)/distbf/csrij 2148 hintt(2) = 0.5d0*(visclc+rkl)/distbf & 2149 + visten(2,iel)*ctheta(iscal)/distbf/csrij 2150 hintt(3) = 0.5d0*(visclc+rkl)/distbf & 2151 + visten(3,iel)*ctheta(iscal)/distbf/csrij 2152 hintt(4) = visten(4,iel)*ctheta(iscal)/distbf/csrij 2153 hintt(5) = visten(5,iel)*ctheta(iscal)/distbf/csrij 2154 hintt(6) = visten(6,iel)*ctheta(iscal)/distbf/csrij 2155 2156 ! Dirichlet Boundary Condition 2157 !----------------------------- 2158 2159 ! Add rho*uk*Tet to T'v' in High Reynolds 2160 do isou = 1, 3 2161 pimpv(isou) = surfbo(isou,ifac)*phit/(surfbn(ifac)*cpp*romc) 2162 enddo 2163 2164 call set_dirichlet_vector_aniso & 2165 !======================== 2166 ( coefaut(:,ifac) , cofafut(:,ifac) , & 2167 coefbut(:,:,ifac), cofbfut(:,:,ifac), & 2168 pimpv , hintt , rinfiv ) 2169 2170 ! Boundary conditions used in the temperature equation 2171 do isou = 1, 3 2172 cofarut(isou,ifac) = 0.d0 2173 do jsou = 1, 3 2174 cofbrut(isou,jsou,ifac) = 0.d0 2175 enddo 2176 enddo 2177 2178 endif 2179 2180 2181 2182 ! Save the values of T^star and T^+ for post-processing 2183 2184 if (b_f_id.ge.0 .or. iscal.eq.iscalt) then 2185 2186 ! Rough wall function 2187 if (icodcl(ifac,ivar).eq.6) then 2188 phit = cofafp(ifac)+cofbfp(ifac)*theipb(ifac) 2189 ! Imposed flux with wall function for post-processing 2190 elseif (icodcl(ifac,ivar).eq.3) then 2191 phit = rcodcl(ifac,ivar,3) 2192 else 2193 phit = 0.d0 2194 endif 2195 2196 tet = phit/(romc*cpp*max(buk(ifac)*bcfnns(ifac),epzero)) 2197 !FIXME Should be uk rather than ustar? 2198 tet = phit/(romc*cpp*max(buet(ifac),epzero)) 2199 2200 if (b_f_id .ge. 0) bval_s(ifac) = bval_s(ifac) - tplus*tet 2201 2202 if (itplus .ge. 0) tplusp(ifac) = tplus 2203 if (itstar .ge. 0) tstarp(ifac) = tet 2204 2205 if (iscal.eq.iscalt) then 2206 tetmax = max(tet, tetmax) 2207 tetmin = min(tet, tetmin) 2208 tplumx = max(tplus,tplumx) 2209 tplumn = min(tplus,tplumn) 2210 endif 2211 2212 endif 2213 2214 endif ! rough wall condition 2215 2216enddo 2217 2218!---- 2219! End 2220!---- 2221 2222return 2223end subroutine 2224