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 clptur.f90 28!> 29!> \brief Boundary conditions for smooth walls (icodcl = 5). 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 veclocity \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 61!> <a href="../../theory.pdf#wallboundary"><b>wall boundary conditions</b></a> 62!> section of the theory guide for more informations, as well as the 63!> <a href="../../theory.pdf#clptur"><b>clptur</b></a> section. 64!------------------------------------------------------------------------------- 65 66!------------------------------------------------------------------------------- 67! Arguments 68!______________________________________________________________________________. 69! mode name role ! 70!______________________________________________________________________________! 71!> \param[in] nscal total number of scalars 72!> \param[in] isvhb indicator to save exchange coeffient 73!> \param[in,out] icodcl face boundary condition code: 74!> - 1 Dirichlet 75!> - 3 Neumann 76!> - 4 sliding and 77!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 78!> - 5 smooth wall and 79!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 80!> - 6 rough wall and 81!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 82!> - 9 free inlet/outlet 83!> (input mass flux blocked to 0) 84!> \param[in,out] rcodcl boundary condition values: 85!> - rcodcl(1) value of the dirichlet 86!> - rcodcl(2) value of the exterior exchange 87!> coefficient (infinite if no exchange) 88!> - rcodcl(3) value flux density 89!> (negative if gain) in w/m2 90!> -# for the velocity \f$ (\mu+\mu_T) 91!> \gradv \vect{u} \cdot \vect{n} \f$ 92!> -# for the pressure \f$ \Delta t 93!> \grad P \cdot \vect{n} \f$ 94!> -# for a scalar \f$ cp \left( K + 95!> \dfrac{K_T}{\sigma_T} \right) 96!> \grad T \cdot \vect{n} \f$ 97!> \param[in] velipb value of the velocity at \f$ \centip \f$ 98!> of boundary cells 99!> \param[in] rijipb value of \f$ R_{ij} \f$ at \f$ \centip \f$ 100!> of boundary cells 101!> \param[out] visvdr dynamic viscosity after V. Driest damping in 102!> boundary cells 103!> \param[out] hbord exchange coefficient at boundary 104!> \param[in] theipb value of thermal scalar at \f$ \centip \f$ 105!> of boundary cells 106!_______________________________________________________________________________ 107 108subroutine clptur & 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_dim 162integer nlogla, nsubla, iuiptn 163integer f_id_rough, f_id, iustar 164integer f_id_uet, f_id_uk 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 uk, uet, nusury, yplus, dplus, yk 173double precision gredu, temp 174double precision cfnns, cfnnk, cfnne 175double precision sqrcmu, ek 176double precision xnuii, xnuit, xmutlm, mut_lm_dmut 177double precision rcprod 178double precision hflui, hint, pimp, qimp 179double precision eloglo(3,3), alpha(6,6) 180double precision rcodcx, rcodcy, rcodcz, rcodcn 181double precision visclc, visctc, romc , distbf, srfbnf 182double precision cofimp, ypup 183double precision bldr12 184double precision xkip 185double precision rinfiv(3) 186double precision visci(3,3), fikis, viscis, distfi 187double precision fcoefa(6), fcoefb(6), fcofaf(6), fcofbf(6), fcofad(6), fcofbd(6) 188double precision rxx, rxy, rxz, ryy, ryz, rzz, rnnb 189double precision rttb, alpha_rnn 190double precision cpp 191double precision sigmak, sigmae 192double precision liqwt, totwt 193double precision rough_d, duplus, uplus 194double precision rough_t 195double precision dtplus, yplus_t 196double precision coef_mom, coef_momm 197double precision one_minus_ri 198double precision dlmo,dt,theta0,flux 199 200double precision, dimension(:), pointer :: crom 201double precision, dimension(:), pointer :: viscl, visct, cpro_cp, yplbr, ustar 202double precision, dimension(:), pointer :: bpro_rough_d 203double precision, dimension(:), pointer :: bpro_rough_t 204double precision, dimension(:), pointer :: f_uet, f_uk 205double precision, dimension(:), allocatable :: byplus, bdplus, buk 206double precision, dimension(:), allocatable, target :: buet, bcfnns_loc 207double precision, dimension(:), pointer :: cvar_k, cvar_ep, bcfnns 208double precision, dimension(:,:), pointer :: cvar_rij 209 210double precision, dimension(:), pointer :: cvar_totwt, cvar_t, cpro_liqwt 211double precision, dimension(:,:), pointer :: coefau, cofafu, visten 212double precision, dimension(:,:,:), pointer :: coefbu, cofbfu 213double precision, dimension(:), pointer :: coefa_k, coefb_k, coefaf_k, coefbf_k 214double precision, dimension(:), pointer :: coefa_ep, coefaf_ep 215double precision, dimension(:), pointer :: coefb_ep, coefbf_ep 216double precision, dimension(:), pointer :: coefa_omg, coefaf_omg 217double precision, dimension(:), pointer :: coefb_omg, coefbf_omg 218double precision, dimension(:), pointer :: coefa_al, coefaf_al 219double precision, dimension(:), pointer :: coefb_al, coefbf_al 220double precision, dimension(:), pointer :: coefa_phi, coefaf_phi 221double precision, dimension(:), pointer :: coefb_phi, coefbf_phi 222double precision, dimension(:), pointer :: coefa_fb, coefaf_fb 223double precision, dimension(:), pointer :: coefb_fb, coefbf_fb 224double precision, dimension(:), pointer :: coefa_nu, coefaf_nu 225double precision, dimension(:), pointer :: coefb_nu, coefbf_nu 226 227double precision, dimension(:,:), pointer :: coefa_rij, coefaf_rij, coefad_rij 228double precision, dimension(:,:,:), pointer :: coefb_rij, coefbf_rij, coefbd_rij 229 230double precision pimp_lam, pimp_turb, gammap, fep, dep, falpg, falpv, ypsd, fct_bl 231 232integer ntlast , iaff 233data ntlast , iaff /-1 , 0/ 234save ntlast , iaff 235 236type(var_cal_opt) :: vcopt_u, vcopt_rij, vcopt_ep 237 238!=============================================================================== 239! Interfaces 240!=============================================================================== 241 242interface 243 244 subroutine clptur_scalar(iscal , isvhb , icodcl , rcodcl , & 245 byplus , bdplus , buk , bcfnns , & 246 hbord , theipb , & 247 tetmax , tetmin , tplumx , tplumn ) 248 249 implicit none 250 integer iscal, isvhb 251 integer, pointer, dimension(:,:) :: icodcl 252 double precision, pointer, dimension(:,:,:) :: rcodcl 253 double precision, dimension(:) :: byplus, bdplus, buk, bcfnns 254 double precision, pointer, dimension(:) :: hbord, theipb 255 double precision tetmax, tetmin, tplumx, tplumn 256 257 end subroutine clptur_scalar 258 259 subroutine clptur_vector(iscal, isvhb, icodcl, rcodcl, & 260 byplus, bdplus, buk) 261 262 implicit none 263 integer iscal, isvhb 264 integer, pointer, dimension(:,:) :: icodcl 265 double precision, pointer, dimension(:,:,:) :: rcodcl 266 double precision, dimension(:) :: byplus, bdplus, buk 267 268 end subroutine clptur_vector 269 270 end interface 271 272!=============================================================================== 273! 1. Initializations 274!=============================================================================== 275 276! Initialize variables to avoid compiler warnings 277 278cofimp = 0.d0 279ek = 0.d0 280rcprod = 0.d0 281uiptn = 0.d0 282rnnb = 0.d0 283 284rinfiv(1) = rinfin 285rinfiv(2) = rinfin 286rinfiv(3) = rinfin 287 288uet = 1.d0 289utau = 1.d0 290 291! --- Constants 292sqrcmu = sqrt(cmu) 293 294! --- Correction factors for stratification (used in atmospheric models) 295cfnns = 1.d0 296cfnnk = 1.d0 297cfnne = 1.d0 298 299bpro_rough_d => null() 300bpro_rough_t => null() 301 302call field_get_id_try("boundary_roughness", f_id_rough) 303if (f_id_rough.ge.0) then 304 call field_get_val_s(f_id_rough, bpro_rough_d) 305 306 ! same thermal roughness if not specified 307 call field_get_val_s(f_id_rough, bpro_rough_t) 308endif 309 310call field_get_id_try("boundary_thermal_roughness", f_id_rough) 311if (f_id_rough.ge.0) then 312 call field_get_val_s(f_id_rough, bpro_rough_t) 313endif 314 315f_uet => null() 316f_uk => null() 317 318call field_get_id_try("boundary_ustar", f_id_uet) 319if (f_id_uet.ge.0) call field_get_val_s(f_id_uet, f_uet) 320call field_get_id_try("boundary_uk", f_id_uk) 321if (f_id_uk.ge.0) call field_get_val_s(f_id_uk, f_uk) 322 323 324! pointers to y+ if saved 325yplbr => null() 326if (iyplbr.ge.0) call field_get_val_s(iyplbr, yplbr) 327 328if (itytur.eq.3 .and. idirsm.eq.1) call field_get_val_v(ivsten, visten) 329 330ustar => null() 331 332! --- Save wall friction velocity 333 334call field_get_id_try('ustar', iustar) 335if (iustar.ge.0) then 336 call field_get_val_s(iustar, ustar) 337else 338 allocate(buet(nfabor)) 339 ustar => buet 340endif 341 342! --- Gradient and flux boundary conditions 343 344call field_get_coefa_v(ivarfl(iu), coefau) 345call field_get_coefb_v(ivarfl(iu), coefbu) 346call field_get_coefaf_v(ivarfl(iu), cofafu) 347call field_get_coefbf_v(ivarfl(iu), cofbfu) 348 349if (ik.gt.0) then 350 call field_get_coefa_s(ivarfl(ik), coefa_k) 351 call field_get_coefb_s(ivarfl(ik), coefb_k) 352 call field_get_coefaf_s(ivarfl(ik), coefaf_k) 353 call field_get_coefbf_s(ivarfl(ik), coefbf_k) 354 call field_get_key_double(ivarfl(ik), ksigmas, sigmak) 355else 356 coefa_k => null() 357 coefb_k => null() 358 coefaf_k => null() 359 coefbf_k => null() 360endif 361 362if (iep.gt.0) then 363 call field_get_coefa_s(ivarfl(iep), coefa_ep) 364 call field_get_coefb_s(ivarfl(iep), coefb_ep) 365 call field_get_coefaf_s(ivarfl(iep), coefaf_ep) 366 call field_get_coefbf_s(ivarfl(iep), coefbf_ep) 367 call field_get_key_double(ivarfl(iep), ksigmas, sigmae) 368else 369 coefa_ep => null() 370 coefb_ep => null() 371 coefaf_ep => null() 372 coefbf_ep => null() 373endif 374 375if (itytur.eq.3) then! Also have boundary conditions for the momentum equation 376 call field_get_coefa_v(ivarfl(irij), coefa_rij) 377 call field_get_coefb_v(ivarfl(irij), coefb_rij) 378 call field_get_coefaf_v(ivarfl(irij), coefaf_rij) 379 call field_get_coefbf_v(ivarfl(irij), coefbf_rij) 380 call field_get_coefad_v(ivarfl(irij), coefad_rij) 381 call field_get_coefbd_v(ivarfl(irij), coefbd_rij) 382endif 383 384if (ial.gt.0) then 385 call field_get_coefa_s(ivarfl(ial), coefa_al) 386 call field_get_coefb_s(ivarfl(ial), coefb_al) 387 call field_get_coefaf_s(ivarfl(ial), coefaf_al) 388 call field_get_coefbf_s(ivarfl(ial), coefbf_al) 389else 390 coefa_al => null() 391 coefb_al => null() 392 coefaf_al => null() 393 coefbf_al => null() 394endif 395 396if (iomg.gt.0) then 397 call field_get_coefa_s(ivarfl(iomg), coefa_omg) 398 call field_get_coefb_s(ivarfl(iomg), coefb_omg) 399 call field_get_coefaf_s(ivarfl(iomg), coefaf_omg) 400 call field_get_coefbf_s(ivarfl(iomg), coefbf_omg) 401else 402 coefa_omg => null() 403 coefb_omg => null() 404 coefaf_omg => null() 405 coefbf_omg => null() 406endif 407 408if (iphi.gt.0) then 409 call field_get_coefa_s(ivarfl(iphi), coefa_phi) 410 call field_get_coefb_s(ivarfl(iphi), coefb_phi) 411 call field_get_coefaf_s(ivarfl(iphi), coefaf_phi) 412 call field_get_coefbf_s(ivarfl(iphi), coefbf_phi) 413else 414 coefa_phi => null() 415 coefb_phi => null() 416 coefaf_phi => null() 417 coefbf_phi => null() 418endif 419 420if (ifb.gt.0) then 421 call field_get_coefa_s(ivarfl(ifb), coefa_fb) 422 call field_get_coefb_s(ivarfl(ifb), coefb_fb) 423 call field_get_coefaf_s(ivarfl(ifb), coefaf_fb) 424 call field_get_coefbf_s(ivarfl(ifb), coefbf_fb) 425else 426 coefa_fb => null() 427 coefb_fb => null() 428 coefaf_fb => null() 429 coefbf_fb => null() 430endif 431 432if (inusa.gt.0) then 433 call field_get_coefa_s(ivarfl(inusa), coefa_nu) 434 call field_get_coefaf_s(ivarfl(inusa), coefaf_nu) 435 call field_get_coefb_s(ivarfl(inusa), coefb_nu) 436 call field_get_coefbf_s(ivarfl(inusa), coefbf_nu) 437else 438 coefa_nu => null() 439 coefb_nu => null() 440 coefaf_nu => null() 441 coefbf_nu => null() 442endif 443 444! --- Physical quantities 445call field_get_val_s(icrom, crom) 446call field_get_val_s(iviscl, viscl) 447call field_get_val_s(ivisct, visct) 448if (icp.ge.0) then 449 call field_get_val_s(icp, cpro_cp) 450endif 451 452if (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) then 453 call field_get_val_s(ivarfl(ik), cvar_k) 454endif 455if ((iturb.eq.30).or.(iturb.eq.31)) then 456 call field_get_val_s(ivarfl(iep), cvar_ep) 457endif 458 459if (itytur.eq.3) then 460 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt_rij) 461 call field_get_key_struct_var_cal_opt(ivarfl(iep), vcopt_ep) 462 call field_get_val_v(ivarfl(irij), cvar_rij) 463endif 464 465! min. and max. of wall tangential velocity 466uiptmx = -grand 467uiptmn = grand 468 469! min. and max. of wall friction velocity 470uetmax = -grand 471uetmin = grand 472ukmax = -grand 473ukmin = grand 474 475! min. and max. of y+ 476yplumx = -grand 477yplumn = grand 478 479! min. and max. of wall friction of the thermal scalar 480tetmax = -grand 481tetmin = grand 482 483! min. and max. of T+ 484tplumx = -grand 485tplumn = grand 486 487! Counters (turbulent, laminar, reversal, scale correction) 488nlogla = 0 489nsubla = 0 490iuiptn = 0 491 492! Alpha constant for a realisable BC for R12 with the SSG model 493alpha_rnn = 0.47d0 494 495! With v2f type model, (phi-fbar et BL-v2/k) u=0 is set directly, so 496! uiptmx and uiptmn are necessarily 0 497if (itytur.eq.5) then 498 uiptmx = 0.d0 499 uiptmn = 0.d0 500endif 501 502! Pointers to specific fields 503allocate(byplus(nfabor)) 504allocate(bdplus(nfabor)) 505allocate(buk(nfabor)) 506 507! Correction for atmospheric wall functions 508call field_get_id_try("non_neutral_scalar_correction", f_id) 509if (f_id.ge.0) then 510 call field_get_val_s(f_id, bcfnns) 511else 512 allocate(bcfnns_loc(nfabor)) 513 bcfnns => bcfnns_loc 514endif 515 516cvar_t => null() 517cvar_totwt => null() 518cpro_liqwt => null() 519 520if (ippmod(iatmos).ge.1) then 521 theta0 = t0 * (ps / p0)**(rair/cp0) 522 call field_get_val_s(ivarfl(isca(iscalt)), cvar_t) 523 if (ippmod(iatmos).eq.2) then 524 call field_get_val_s(ivarfl(isca(iymw)), cvar_totwt) 525 call field_get_val_s(iliqwt, cpro_liqwt) 526 endif 527endif 528 529! --- Loop on boundary faces 530do ifac = 1, nfabor 531 532 ! Test on the presence of a smooth wall condition (start) 533 if (icodcl(ifac,iu).eq.5) then 534 535 iel = ifabor(ifac) 536 537 ! Physical properties 538 visclc = viscl(iel) 539 visctc = visct(iel) 540 romc = crom(iel) 541 542 ! Geometric quantities 543 distbf = distb(ifac) 544 srfbnf = surfbn(ifac) 545 546 !=========================================================================== 547 ! 1. Local framework 548 !=========================================================================== 549 550 ! Unit normal 551 552 rnx = surfbo(1,ifac)/srfbnf 553 rny = surfbo(2,ifac)/srfbnf 554 rnz = surfbo(3,ifac)/srfbnf 555 556 ! Handle displacement velocity 557 558 rcodcx = rcodcl(ifac,iu,1) 559 rcodcy = rcodcl(ifac,iv,1) 560 rcodcz = rcodcl(ifac,iw,1) 561 562 ! If we are not using ALE, force the displacement velocity for the face 563 ! to be tangential (and update rcodcl for possible use) 564 ! In frozen rotor (iturbo = 1), the velocity is neither tangential to the 565 ! wall (absolute velocity solved in a relative frame of reference) 566 if (iale.eq.0.and.iturbo.eq.0) then 567 rcodcn = rcodcx*rnx+rcodcy*rny+rcodcz*rnz 568 rcodcx = rcodcx -rcodcn*rnx 569 rcodcy = rcodcy -rcodcn*rny 570 rcodcz = rcodcz -rcodcn*rnz 571 rcodcl(ifac,iu,1) = rcodcx 572 rcodcl(ifac,iv,1) = rcodcy 573 rcodcl(ifac,iw,1) = rcodcz 574 endif 575 576 ! Relative tangential velocity 577 578 upx = velipb(ifac,1) - rcodcx 579 upy = velipb(ifac,2) - rcodcy 580 upz = velipb(ifac,3) - rcodcz 581 582 usn = upx*rnx+upy*rny+upz*rnz 583 tx = upx -usn*rnx 584 ty = upy -usn*rny 585 tz = upz -usn*rnz 586 txn = sqrt(tx**2 +ty**2 +tz**2) 587 utau= txn 588 589 ! Unit tangent 590 591 if (txn.ge.epzero) then 592 593 txn0 = 1.d0 594 595 tx = tx/txn 596 ty = ty/txn 597 tz = tz/txn 598 599 else 600 601 ! If the velocity is zero, 602 ! Tx, Ty, Tz is not used (we cancel the velocity), so we assign any 603 ! value (zero for example) 604 605 txn0 = 0.d0 606 607 tx = 0.d0 608 ty = 0.d0 609 tz = 0.d0 610 611 endif 612 613 ! Complete if necessary for Rij-Epsilon 614 615 if (itytur.eq.3) then 616 617 ! --> T2 = RN X T (where X is the cross product) 618 619 t2x = rny*tz - rnz*ty 620 t2y = rnz*tx - rnx*tz 621 t2z = rnx*ty - rny*tx 622 623 ! --> Orthogonal matrix for change of reference frame ELOGLOij 624 ! (from local to global reference frame) 625 626 ! |TX -RNX T2X| 627 ! ELOGLO = |TY -RNY T2Y| 628 ! |TZ -RNZ T2Z| 629 630 ! Its transpose ELOGLOt is its inverse 631 632 eloglo(1,1) = tx 633 eloglo(1,2) = -rnx 634 eloglo(1,3) = t2x 635 eloglo(2,1) = ty 636 eloglo(2,2) = -rny 637 eloglo(2,3) = t2y 638 eloglo(3,1) = tz 639 eloglo(3,2) = -rnz 640 eloglo(3,3) = t2z 641 642 ! Compute Reynolds stress transformation matrix 643 644 clsyme = 0 645 call turbulence_bc_rij_transform(clsyme, eloglo, alpha) 646 647 endif 648 649 !=========================================================================== 650 ! 2. Friction velocities 651 !=========================================================================== 652 653 ! ---> Compute Uet depending if we are in the log zone or not 654 ! in 1 or 2 velocity scales 655 ! and uk based on ek 656 657 658 if (abs(utau).le.epzero) utau = epzero 659 660 nusury = visclc/(distbf*romc) 661 xnuii = visclc/romc 662 xnuit = visctc/romc 663 664 if (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) then 665 ek = cvar_k(iel) 666 ! TODO: we could add 2*nu_T dv/dy to rnnb 667 rnnb = 2.d0 / 3.d0 * ek 668 else if (itytur.eq.3) then 669 ek = 0.5d0*(cvar_rij(1,iel)+cvar_rij(2,iel)+cvar_rij(3,iel)) 670 rxx = cvar_rij(1,iel) 671 rxy = cvar_rij(4,iel) 672 rxz = cvar_rij(6,iel) 673 ryy = cvar_rij(2,iel) 674 ryz = cvar_rij(5,iel) 675 rzz = cvar_rij(3,iel) 676 rnnb = rnx * (rxx * rnx + rxy * rny + rxz * rnz) & 677 + rny * (rxy * rnx + ryy * rny + ryz * rnz) & 678 + rnz * (rxz * rnx + ryz * rny + rzz * rnz) 679 680 rttb = tx * (rxx * tx + rxy * ty + rxz * tz) & 681 + ty * (rxy * tx + ryy * ty + ryz * tz) & 682 + tz * (rxz * tx + ryz * ty + rzz * tz) 683 endif 684 685 if (f_id_rough.ge.0) then 686 rough_d = bpro_rough_d(ifac) 687 else 688 rough_d = 0.d0 689 endif 690 691 call wallfunctions & 692 ( iwallf, ifac , & 693 xnuii , xnuit , utau , distbf, rough_d, rnnb, ek, & 694 iuntur, nsubla, nlogla, & 695 uet , uk , yplus , ypup , cofimp, dplus ) 696 697 ! Louis or Monin Obukhov wall function for atmospheric flows 698 if (ippmod(iatmos).ge.1.and.(iwalfs.eq.2.or.iwalfs.eq.3)) then 699 ! Louis 700 if (iwalfs.eq.2) then 701 702 ! Compute reduced gravity for non horizontal walls : 703 gredu = gx*rnx + gy*rny + gz*rnz 704 705 temp = cvar_t(iel) 706 totwt = 0.d0 707 liqwt = 0.d0 708 709 if (ippmod(iatmos).eq.2) then 710 totwt = cvar_totwt(iel) 711 liqwt = cpro_liqwt(iel) 712 endif 713 714 yk = distbf * uk / xnuii 715 ! 1/U+ for neutral 716 duplus = ypup / yk 717 718 rough_t = bpro_rough_t(ifac) 719 ! 1/T+ 720 ! "y+_t" tends to "y/rough_t" for rough regime and to "y+k" 721 ! times a shift for smooth regime 722 ! 723 ! Rough regime reads: 724 ! T+ = 1/kappa ln(y/rough_t) = 1/kappa ln(y/zeta) + 8.5 725 ! = 1/kappa ln[y/zeta * exp(8.5 kappa)] 726 ! 727 ! Note zeta_t = rough_t * exp(8.5 kappa) 728 ! 729 ! Smooth regime reads: 730 ! T+ = 1/kappa ln(y uk/nu) + 5.2 = 1/kappa ln[y uk*exp(5.2 kappa) / nu] 731 ! 732 ! Mixed regime reads: 733 ! T+ = 1/kappa ln[y uk*exp(5.2 kappa)/(nu + alpha uk zeta)] 734 ! = 1/kappa ln[y uk*exp(5.2 kappa)/(nu + alpha uk rough_t * exp(8.5 kappa))] 735 ! = 1/kappa ln[y uk*exp(5.2 kappa)/(nu + alpha uk rough_t * exp(8.5 kappa))] 736 ! with 737 ! alpha * exp(8.5 kappa) / exp(5.2 kappa) = 1 738 ! ie 739 ! alpha = exp(-(8.5-5.2) kappa) = 0.25 740 ! so 741 ! T+ = 1/kappa ln[y uk*exp(5.2 kappa)/(nu + uk rough_t * exp(5.2 kappa))] 742 ! = 1/kappa ln[y+k / (exp(-5.2 kappa) + uk rough_t/nu)] 743 ! 744 745 ! shifted y+ 746 yplus_t = yk / (exp(-xkappa * 5.2d0) + uk *rough_t / xnuii)! FIXME use log constant 747 ! 1/T+ for neutral 748 dtplus = xkappa / log(yplus_t) 749 750 call atmcls & 751 !========== 752 ( ifac , & 753 utau , rough_d, duplus , dtplus , & 754 yplus_t, & 755 uet , & 756 gredu , & 757 cfnns , cfnnk , cfnne , & 758 dlmo , & 759 temp , totwt , liqwt , & 760 icodcl , rcodcl ) 761 762 ! Monin Obukhov wall function 763 else if (iwalfs.eq.3) then 764 765 ! Compute local LMO 766 if (ippmod(iatmos).ge.1) then 767 gredu = gx*rnx + gy*rny + gz*rnz 768 769 if (icodcl(ifac,isca(iscalt)).eq.6 & 770 .or.icodcl(ifac,isca(iscalt)).eq.5) then 771 772 dt = theipb(ifac)-rcodcl(ifac,isca(iscalt),1) 773 call mo_compute_from_thermal_diff(distbf,rough_d,utau,dt, & 774 theta0, gredu, & 775 dlmo, uet) 776 777 elseif (icodcl(ifac,isca(iscalt)).eq.3) then 778 if (icp.ge.0) then 779 cpp = cpro_cp(iel) 780 else 781 cpp = cp0 782 endif 783 784 flux = rcodcl(ifac, isca(iscalt),3)/romc/cpp 785 call mo_compute_from_thermal_flux(distbf,rough_d,utau,flux, & 786 theta0, gredu, & 787 dlmo, uet) 788 789 endif 790 791 else 792 793 ! No temperature delta: neutral 794 call mo_compute_from_thermal_diff(distbf,rough_d,utau,0.d0,0.d0,0.d0, & 795 dlmo,uet) 796 797 endif 798 799 ! Take stability into account for the turbulent velocity scale 800 coef_mom = cs_mo_phim(distbf+rough_d,dlmo) 801 one_minus_ri = 1.d0-(distbf+rough_d) * dlmo/coef_mom 802 803 if (one_minus_ri.gt.0) then 804 ! Warning: overwritting uk, yplus should be recomputed 805 uk = uk / one_minus_ri**0.25d0 806 yplus = distbf * uk / xnuii 807 808 ! Epsilon should be modified as well to get P+G = P(1-Ri) = epsilon 809 ! P = -R_tn dU/dn = uk^2 uet Phi_m / (kappa z) 810 cfnne = one_minus_ri * coef_mom 811 ! Nothing done for the moment for really high stability 812 else 813 cfnne = 1.d0 814 endif 815 ! Boundary condition on the velocity to have approximately the good 816 ! turbulence production 817 coef_mom = cs_mo_phim(distbf+rough_d,dlmo) 818 coef_momm = cs_mo_phim(2.d0*distbf+rough_d,dlmo) 819 rcprod = 2.d0*distbf*sqrt(xkappa*uk*romc*coef_mom/visctc/(distbf+rough_d)) & 820 - coef_momm/(2.d0+rough_d/distbf) 821 822 iuntur = 1 823 824 uplus = utau / uet 825 ! Coupled solving of the velocity components 826 ! The boundary term for velocity gradient is implicit 827 ! modified for non neutral boundary layer (in uplus) 828 cofimp = min(max(1.d0 - 1.d0/(xkappa*uplus)*rcprod,0.d0),1.d0) 829 yk = distbf * uk / xnuii 830 endif ! End Monin Obukhov 831 ! Dimensionless velocity, recomputed and therefore may take stability 832 ! into account 833 uplus = utau / uet 834 ! y+/U+ for non neutral is recomputed 835 ypup = yk / max(uplus, epzero) 836 837 endif 838 ! Store u_star and uk if needed 839 if(f_id_uet.ge.0) then 840 f_uet(ifac) = uet 841 f_uk(ifac) = uk 842 endif 843 844 uetmax = max(uet,uetmax) 845 uetmin = min(uet,uetmin) 846 ukmax = max(uk,ukmax) 847 ukmin = min(uk,ukmin) 848 yplumx = max(yplus,yplumx) 849 yplumn = min(yplus,yplumn) 850 851 if (iustar.ge.0) then 852 ustar(ifac) = uet !TODO remove, this information is in cofaf cofbf 853 endif 854 855 ! save turbulent subgrid viscosity after van Driest damping in LES 856 ! care is taken to not dampen it twice at boundary cells having more 857 ! one boundary face 858 if (itytur.eq.4.and.idries.eq.1) then 859 if (visvdr(iel).lt.-900.d0) then 860 visct(iel) = visct(iel)*(1.d0-exp(-(yplus)/cdries))**2 861 visvdr(iel) = visct(iel) 862 visctc = visct(iel) 863 endif 864 endif 865 866 ! Save yplus if post-processed or condensation modelling 867 if (iyplbr.ge.0) then 868 yplbr(ifac) = yplus 869 endif 870 871 !=========================================================================== 872 ! 3. Velocity boundary conditions 873 !=========================================================================== 874 875 ! Deprecated power law (Werner & Wengle) 876 if (iwallf.eq.1) then 877 uiptn = utau + uet*apow*bpow*yplus**bpow*(2.d0**(bpow-1.d0)-2.d0) 878 879 ! Dependant on the turbulence Model 880 else 881 882 ! uiptn respecte la production de k 883 ! de facon conditionnelle --> Coef RCPROD 884 885 ! --> k-epsilon and k-omega 886 !-------------------------- 887 if (itytur.eq.2.or.iturb.eq.60) then 888 889 xmutlm = xkappa*visclc*(yplus + dplus)!FIXME should be efvisc... 890 if (cell_is_active(iel).eq.1) then 891 mut_lm_dmut = xmutlm/visctc 892 else 893 mut_lm_dmut = 0.d0 894 endif 895 896 ! If yplus=0, uiptn is set to 0 to avoid division by 0. 897 ! By the way, in this case: iuntur=0 898 if (yplus.gt.epzero) then !TODO use iuntur.eq.1 899 rcprod = min(xkappa , max(1.0d0,sqrt(mut_lm_dmut))/(yplus+dplus))!FIXME not valid for rough 900 901 uiptn = utau - distbf*uet*uk*romc/xkappa/visclc*( & 902 2.0d0*rcprod - 1.0d0/(2.0d0*yplus+dplus)) 903 else 904 uiptn = 0.d0 905 endif 906 907 ! --> No turbulence, mixing length or Rij-espilon 908 !------------------------------------------------ 909 elseif (iturb.eq.0.or.iturb.eq.10.or.itytur.eq.3) then 910 911 ! Dans le cadre de la ponderation elliptique, on ne doit pas 912 ! tenir compte des lois de paroi. On fait donc un test sur le modele 913 ! de turbulence : 914 ! si on est en LRR ou SSG on laisse les lois de paroi, si on est en 915 ! EBRSM, on impose l adherence. 916 if (iturb.eq.32.or.iturb.eq.0) then 917 uiptn = 0.d0 918 else 919 920 ! If yplus=0, uiptn is set to 0 to avoid division by 0. 921 ! By the way, in this case: iuntur=0 922 if (yplus.gt.epzero) then !FIXME use iuntur 923 uiptn = utau - distbf*romc*uet*uk/xkappa/visclc & 924 *(2.0d0/(yplus+dplus) - 1.0d0/(2.0d0*yplus+dplus)) 925 else 926 uiptn = 0.d0 927 endif 928 929 endif 930 931 ! --> LES and Spalart Allmaras 932 !----------------------------- 933 elseif (itytur.eq.4.or.iturb.eq.70) then 934 935 uiptn = utau - uet/xkappa*1.5d0 936 937 ! If (mu+mut) becomes zero (dynamic models), an arbitrary value is set 938 ! (zero flux) but without any problems because the flux 939 ! is really zero at this face. 940 if (visctc+visclc.le.0) then 941 hflui = 0.d0 942 endif 943 944 ! --> v2f 945 !-------- 946 elseif (itytur.eq.5) then 947 948 ! Avec ces conditions, pas besoin de calculer uiptmx, uiptmn 949 ! et iuiptn qui sont nuls (valeur d'initialisation) 950 uiptn = 0.d0 951 952 endif 953 endif 954 955 ! Min and Max and counter of reversal layer 956 uiptmn = min(uiptn*iuntur,uiptmn) 957 uiptmx = max(uiptn*iuntur,uiptmx) 958 if (uiptn*iuntur.lt.-epzero) iuiptn = iuiptn + 1 959 960 ! To be coherent with a wall function, clip it to 0 961 cofimp = max(cofimp, 0.d0) 962 963 ! On implicite le terme (rho*uet*uk) 964 hflui = visclc / distbf * ypup 965 966 if (itytur.eq.3) then 967 hint = visclc /distbf 968 else 969 hint = (visclc + visctc)/distbf 970 endif 971 972 ! Gradient boundary conditions 973 !----------------------------- 974 rcodcn = rcodcx*rnx+rcodcy*rny+rcodcz*rnz 975 976 coefau(1,ifac) = (1.d0-cofimp)*(rcodcx - rcodcn*rnx) + rcodcn*rnx 977 coefau(2,ifac) = (1.d0-cofimp)*(rcodcy - rcodcn*rny) + rcodcn*rny 978 coefau(3,ifac) = (1.d0-cofimp)*(rcodcz - rcodcn*rnz) + rcodcn*rnz 979 980 ! Projection in order to have the velocity parallel to the wall 981 ! B = cofimp * ( IDENTITY - n x n ) 982 983 coefbu(1,1,ifac) = cofimp*(1.d0-rnx**2) 984 coefbu(2,2,ifac) = cofimp*(1.d0-rny**2) 985 coefbu(3,3,ifac) = cofimp*(1.d0-rnz**2) 986 coefbu(1,2,ifac) = -cofimp*rnx*rny 987 coefbu(1,3,ifac) = -cofimp*rnx*rnz 988 coefbu(2,1,ifac) = -cofimp*rny*rnx 989 coefbu(2,3,ifac) = -cofimp*rny*rnz 990 coefbu(3,1,ifac) = -cofimp*rnz*rnx 991 coefbu(3,2,ifac) = -cofimp*rnz*rny 992 993 ! Flux boundary conditions 994 !------------------------- 995 996 cofafu(1,ifac) = -hflui*(rcodcx - rcodcn*rnx) - hint*rcodcn*rnx 997 cofafu(2,ifac) = -hflui*(rcodcy - rcodcn*rny) - hint*rcodcn*rny 998 cofafu(3,ifac) = -hflui*(rcodcz - rcodcn*rnz) - hint*rcodcn*rnz 999 1000 ! Projection in order to have the shear stress parallel to the wall 1001 ! B = hflui*( IDENTITY - n x n ) 1002 1003 cofbfu(1,1,ifac) = hflui*(1.d0-rnx**2) + hint*rnx**2 1004 cofbfu(2,2,ifac) = hflui*(1.d0-rny**2) + hint*rny**2 1005 cofbfu(3,3,ifac) = hflui*(1.d0-rnz**2) + hint*rnz**2 1006 1007 cofbfu(1,2,ifac) = (hint - hflui)*rnx*rny 1008 cofbfu(1,3,ifac) = (hint - hflui)*rnx*rnz 1009 cofbfu(2,1,ifac) = (hint - hflui)*rny*rnx 1010 cofbfu(2,3,ifac) = (hint - hflui)*rny*rnz 1011 cofbfu(3,1,ifac) = (hint - hflui)*rnz*rnx 1012 cofbfu(3,2,ifac) = (hint - hflui)*rnz*rny 1013 1014 ! In case of transient turbomachinery computations, save the coefficents 1015 ! associated to turbulent wall velocity BC, in order to update the wall 1016 ! velocity after the geometry update(between prediction and correction step) 1017 if (iturbo.eq.2) then 1018 if (irotce(iel).ne.0) then 1019 coftur(ifac) = cofimp 1020 hfltur(ifac) = hflui 1021 endif 1022 endif 1023 1024 !=========================================================================== 1025 ! 4. Boundary conditions on k and epsilon 1026 !=========================================================================== 1027 1028 if (itytur.eq.2) then 1029 1030 ! ================================== 1031 ! Launder Sharma boundary conditions 1032 ! ================================== 1033 if(iturb.eq.22) then 1034 1035 ! Dirichlet Boundary Condition on k 1036 !---------------------------------- 1037 if (iwallf.eq.0) then 1038 ! No wall functions forced by user 1039 pimp = 0.d0 1040 else 1041 ! Use of wall functions 1042 if (iuntur.eq.1) then 1043 pimp = uk**2/sqrcmu 1044 else 1045 pimp = 0.d0 1046 endif 1047 endif 1048 pimp = pimp * cfnnk 1049 hint = (visclc+visctc/sigmak)/distbf 1050 1051 call set_dirichlet_scalar & 1052 !==================== 1053 ( coefa_k(ifac), coefaf_k(ifac), & 1054 coefb_k(ifac), coefbf_k(ifac), & 1055 pimp , hint , rinfin ) 1056 1057 1058 ! Dirichlet Boundary Condition on epsilon tilda 1059 !--------------------------------------------- 1060 pimp_lam = 0.d0 1061 1062 if (iwallf.eq.0) then 1063 ! No wall functions forced by user 1064 pimp = pimp_lam 1065 else 1066 ! Use of wall functions 1067 if (yplus .gt. epzero) then 1068 pimp_turb = 5.d0*uk**4*romc/ & 1069 (xkappa*visclc*yplus) 1070 1071 ! Blending function, from JF Wald PhD (2016) 1072 fct_bl = exp(-0.674d-3*yplus**3.d0) 1073 pimp = pimp_lam*fct_bl + pimp_turb*(1.d0-fct_bl) 1074 fep = exp(-((yplus+dplus)/4.d0)**1.5d0) 1075 dep = 1.d0- exp(-((yplus+dplus)/9.d0)**2.1d0) 1076 pimp = fep*pimp_lam + (1.d0-fep)*dep*pimp_turb 1077 else 1078 pimp = pimp_lam 1079 endif 1080 endif 1081 1082 hint = (visclc+visctc/sigmae)/distbf 1083 pimp = pimp * cfnne 1084 call set_dirichlet_scalar & 1085 !==================== 1086 ( coefa_ep(ifac), coefaf_ep(ifac), & 1087 coefb_ep(ifac), coefbf_ep(ifac), & 1088 pimp , hint , rinfin ) 1089 1090 1091 ! =================================== 1092 ! Quadratic Baglietto k-epsilon model 1093 ! =================================== 1094 else if(iturb.eq.23) then 1095 1096 ! Dirichlet Boundary Condition on k 1097 !---------------------------------- 1098 if (iwallf.eq.0) then 1099 ! No wall functions forces by user 1100 pimp = 0.d0 1101 else 1102 ! Use of wall functions 1103 if (iuntur.eq.1) then 1104 pimp = uk**2/sqrcmu 1105 else 1106 pimp = 0.d0 1107 endif 1108 endif 1109 1110 hint = (visclc+visctc/sigmak)/distbf 1111 pimp = pimp * cfnnk 1112 call set_dirichlet_scalar & 1113 !==================== 1114 ( coefa_k(ifac), coefaf_k(ifac), & 1115 coefb_k(ifac), coefbf_k(ifac), & 1116 pimp , hint , rinfin ) 1117 1118 ! Dirichlet Boundary Condition on epsilon 1119 !--------------------------------------- 1120 if(iwallf.ne.0) then 1121 1122 pimp_lam = 2.0d0*visclc/romc*cvar_k(iel)/distbf**2 1123 1124 if (yplus.gt.epzero) then 1125 pimp_turb = 5.d0*uk**4*romc/ & 1126 (xkappa*visclc*yplus) 1127 1128 ! Blending between wall and homogeneous layer 1129 fep = exp(-((yplus+dplus)/4.d0)**1.5d0) 1130 dep = 1.d0- exp(-((yplus+dplus)/9.d0)**2.1d0) 1131 pimp = fep*pimp_lam + (1.d0-fep)*dep*pimp_turb 1132 else 1133 pimp = pimp_lam 1134 end if 1135 1136 else 1137 1138 pimp = 2.0d0*visclc/romc*cvar_k(iel)/distbf**2 1139 end if 1140 pimp = pimp * cfnne 1141 call set_dirichlet_scalar & 1142 !==================== 1143 ( coefa_ep(ifac), coefaf_ep(ifac), & 1144 coefb_ep(ifac), coefbf_ep(ifac), & 1145 pimp , hint , rinfin ) 1146 1147 ! ============================================== 1148 ! k-epsilon and k-epsilon LP boundary conditions 1149 ! ============================================== 1150 else 1151 1152 ! Dirichlet Boundary Condition on k 1153 !---------------------------------- 1154 1155 if (iuntur.eq.1) then 1156 pimp = uk**2/sqrcmu 1157 else 1158 pimp = 0.d0 1159 endif 1160 hint = (visclc+visctc/sigmak)/distbf 1161 pimp = pimp * cfnnk 1162 call set_dirichlet_scalar & 1163 !==================== 1164 ( coefa_k(ifac), coefaf_k(ifac), & 1165 coefb_k(ifac), coefbf_k(ifac), & 1166 pimp , hint , rinfin ) 1167 1168 1169 ! Neumann Boundary Condition on epsilon 1170 !-------------------------------------- 1171 1172 hint = (visclc+visctc/sigmae)/distbf 1173 1174 ! If yplus=0, uiptn is set to 0 to avoid division by 0. 1175 ! By the way, in this case: iuntur=0 1176 if (yplus.gt.epzero.and.iuntur.eq.1) then !FIXME use only iuntur 1177 xnuii = visclc/romc 1178 pimp = distbf*4.d0*uk**5/ & 1179 (xkappa*xnuii**2*(yplus+2.d0*dplus)**2) 1180 1181 qimp = -pimp*hint !TODO transform it, it is only to be fully equivalent 1182 else 1183 qimp = 0.d0 1184 endif 1185 pimp = pimp * cfnne 1186 call set_neumann_scalar & 1187 !================== 1188 ( coefa_ep(ifac), coefaf_ep(ifac), & 1189 coefb_ep(ifac), coefbf_ep(ifac), & 1190 qimp , hint ) 1191 1192 end if 1193 1194 !=========================================================================== 1195 ! 5. Boundary conditions on Rij-epsilon 1196 !=========================================================================== 1197 1198 elseif (itytur.eq.3) then 1199 1200 ! Exchange coefficient 1201 1202 ! Symmetric tensor diffusivity (Daly Harlow -- GGDH) 1203 if (iand(vcopt_rij%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 1204 1205 visci(1,1) = visclc + visten(1,iel) 1206 visci(2,2) = visclc + visten(2,iel) 1207 visci(3,3) = visclc + visten(3,iel) 1208 visci(1,2) = visten(4,iel) 1209 visci(2,1) = visten(4,iel) 1210 visci(2,3) = visten(5,iel) 1211 visci(3,2) = visten(5,iel) 1212 visci(1,3) = visten(6,iel) 1213 visci(3,1) = visten(6,iel) 1214 1215 ! ||Ki.S||^2 1216 viscis = ( visci(1,1)*surfbo(1,ifac) & 1217 + visci(1,2)*surfbo(2,ifac) & 1218 + visci(1,3)*surfbo(3,ifac))**2 & 1219 + ( visci(2,1)*surfbo(1,ifac) & 1220 + visci(2,2)*surfbo(2,ifac) & 1221 + visci(2,3)*surfbo(3,ifac))**2 & 1222 + ( visci(3,1)*surfbo(1,ifac) & 1223 + visci(3,2)*surfbo(2,ifac) & 1224 + visci(3,3)*surfbo(3,ifac))**2 1225 1226 ! IF.Ki.S 1227 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1228 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1229 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1230 )*surfbo(1,ifac) & 1231 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1232 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1233 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1234 )*surfbo(2,ifac) & 1235 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1236 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1237 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1238 )*surfbo(3,ifac) 1239 1240 distfi = distb(ifac) 1241 1242 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 1243 ! NB: eps =1.d-1 must be consistent with vitens.f90 1244 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1245 1246 hint = viscis/surfbn(ifac)/fikis 1247 1248 ! Scalar diffusivity 1249 else 1250 hint = (visclc+visctc*csrij/cmu)/distbf 1251 endif 1252 1253 ! ---> Tensor Rij (Partially or totally implicited) 1254 1255 do isou = 1, 6 1256 fcoefa(isou) = 0.0d0 1257 fcoefb(isou) = 0.0d0 1258 fcofad(isou) = 0.0d0 1259 fcofbd(isou) = 0.0d0 1260 enddo 1261 1262 ! blending factor so that the component R(n,tau) have only 1263 ! -mu_T/(mu+mu_T)*uet*uk 1264 bldr12 = visctc/(visclc + visctc) 1265 1266 do isou = 1, 6 1267 1268 if (isou.eq.1) then 1269 jj = 1 1270 kk = 1 1271 else if (isou.eq.2) then 1272 jj = 2 1273 kk = 2 1274 else if (isou.eq.3) then 1275 jj = 3 1276 kk = 3 1277 else if (isou.eq.4) then 1278 jj = 1 1279 kk = 2 1280 else if (isou.eq.5) then 1281 jj = 2 1282 kk = 3 1283 else if (isou.eq.6) then 1284 jj = 1 1285 kk = 3 1286 endif 1287 1288 ! LRR and the Standard SGG or EB-RSM + wall functions 1289 if (((iturb.eq.30).or.(iturb.eq.31).and.iuntur.eq.1).or. & 1290 (iturb.eq.32.and.iwallf.ne.0.and.yplus.gt.epzero)) then 1291 1292 if (irijco.eq.1) then 1293 coefa_rij(isou, ifac) = - ( eloglo(jj,1)*eloglo(kk,2) & 1294 + eloglo(jj,2)*eloglo(kk,1)) & 1295 * alpha_rnn * sqrt(rnnb * rttb)*cfnnk 1296 coefaf_rij(isou, ifac) = -hint * coefa_rij(isou, ifac) 1297 coefad_rij(isou, ifac) = 0.d0 1298 do ii = 1, 6 1299 coefb_rij(isou,ii, ifac) = alpha(ii, isou) 1300 if (ii.eq.isou) then 1301 coefbf_rij(isou,ii, ifac) = hint * (1.d0 - coefb_rij(isou,ii, ifac)) 1302 else 1303 coefbf_rij(isou,ii, ifac) = - hint * coefb_rij(isou,ii, ifac) 1304 endif 1305 coefbd_rij(isou,ii, ifac) = coefb_rij(isou,ii, ifac) 1306 enddo 1307 1308 else if ((iclptr.eq.1)) then 1309 1310 do ii = 1, 6 1311 if (ii.ne.isou) then 1312 fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii) 1313 endif 1314 enddo 1315 fcoefb(isou) = alpha(isou,isou) 1316 else 1317 do ii = 1, 6 1318 fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii) 1319 enddo 1320 fcoefb(isou) = 0.d0 1321 endif 1322 1323 ! Boundary conditions for the momentum equation 1324 fcofad(isou) = fcoefa(isou) 1325 fcofbd(isou) = fcoefb(isou) 1326 1327 fcoefa(isou) = fcoefa(isou) & 1328 - ( eloglo(jj,1)*eloglo(kk,2) & 1329 + eloglo(jj,2)*eloglo(kk,1))*bldr12*uet*uk*cfnnk 1330 1331 ! In the viscous sublayer or for EBRSM: zero Reynolds' stresses 1332 else 1333 if (irijco.eq.1) then 1334 coefa_rij(isou, ifac) = 0.d0 1335 coefaf_rij(isou, ifac) = 0.d0 1336 coefad_rij(isou, ifac) = 0.d0 1337 do ii = 1, 6 1338 coefb_rij(isou,ii, ifac) = 0.d0 1339 if (ii.eq.isou) then 1340 coefbf_rij(isou,ii, ifac) = hint 1341 else 1342 coefbf_rij(isou,ii, ifac) = 0.d0 1343 endif 1344 coefbd_rij(isou,ii, ifac) = 0.d0 1345 enddo 1346 else 1347 fcoefa(isou) = 0.d0 1348 fcofad(isou) = 0.d0 1349 fcoefb(isou) = 0.d0 1350 fcofbd(isou) = 0.d0 1351 endif 1352 1353 endif 1354 1355 ! Translate into Diffusive flux BCs 1356 fcofaf(isou) = -hint*fcoefa(isou) 1357 fcofbf(isou) = hint*(1.d0-fcoefb(isou)) 1358 enddo 1359 1360 if (irijco.ne.1) then 1361 do isou = 1, 6 1362 coefa_rij(isou,ifac) = fcoefa(isou) 1363 coefaf_rij(isou,ifac) = fcofaf(isou) 1364 coefad_rij(isou,ifac) = fcofad(isou) 1365 do ii = 1,6 1366 coefb_rij(isou,ii,ifac) = 0 1367 coefbf_rij(isou,ii,ifac) = 0 1368 coefbd_rij(isou,ii,ifac) = 0 1369 enddo 1370 coefb_rij(isou,isou,ifac) = fcoefb(isou) 1371 coefbf_rij(isou,isou,ifac) = fcofbf(isou) 1372 coefbd_rij(isou,isou,ifac) = fcofbd(isou) 1373 enddo 1374 endif 1375 1376 ! ---> Epsilon 1377 ! NB: no reconstruction, possibility of partial implicitation 1378 1379 ! Symmetric tensor diffusivity (Daly Harlow -- GGDH) 1380 if (iand(vcopt_ep%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1381 1382 visci(1,1) = visclc + visten(1,iel)/sigmae 1383 visci(2,2) = visclc + visten(2,iel)/sigmae 1384 visci(3,3) = visclc + visten(3,iel)/sigmae 1385 visci(1,2) = visten(4,iel)/sigmae 1386 visci(2,1) = visten(4,iel)/sigmae 1387 visci(2,3) = visten(5,iel)/sigmae 1388 visci(3,2) = visten(5,iel)/sigmae 1389 visci(1,3) = visten(6,iel)/sigmae 1390 visci(3,1) = visten(6,iel)/sigmae 1391 1392 ! ||Ki.S||^2 1393 viscis = ( visci(1,1)*surfbo(1,ifac) & 1394 + visci(1,2)*surfbo(2,ifac) & 1395 + visci(1,3)*surfbo(3,ifac))**2 & 1396 + ( visci(2,1)*surfbo(1,ifac) & 1397 + visci(2,2)*surfbo(2,ifac) & 1398 + visci(2,3)*surfbo(3,ifac))**2 & 1399 + ( visci(3,1)*surfbo(1,ifac) & 1400 + visci(3,2)*surfbo(2,ifac) & 1401 + visci(3,3)*surfbo(3,ifac))**2 1402 1403 ! IF.Ki.S 1404 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 1405 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 1406 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 1407 )*surfbo(1,ifac) & 1408 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 1409 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 1410 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 1411 )*surfbo(2,ifac) & 1412 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 1413 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 1414 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 1415 )*surfbo(3,ifac) 1416 1417 distfi = distb(ifac) 1418 1419 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 1420 ! NB: eps =1.d-1 must be consistent with vitens.f90 1421 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 1422 1423 hint = viscis/surfbn(ifac)/fikis 1424 1425 ! Scalar diffusivity 1426 else 1427 hint = (visclc+visctc/sigmae)/distbf 1428 endif 1429 1430 if ((iturb.eq.30).or.(iturb.eq.31)) then 1431 1432 ! Si yplus=0, on met coefa a 0 directement pour eviter une division 1433 ! par 0. 1434 if (yplus.gt.epzero.and.iuntur.eq.1) then 1435 xnuii = visclc/romc 1436 pimp = distbf*4.d0*uk**5/ & 1437 (xkappa*xnuii**2*(yplus+2.d0*dplus)**2) 1438 else 1439 pimp = 0.d0 1440 endif 1441 1442 ! Neumann Boundary Condition 1443 !--------------------------- 1444 1445 if (iclptr.eq.1) then !TODO not available for k-eps 1446 1447 qimp = -pimp*hint !TODO transform it, it is only to be fully equivalent 1448 pimp = pimp * cfnne 1449 call set_neumann_scalar & 1450 !================== 1451 ( coefa_ep(ifac), coefaf_ep(ifac), & 1452 coefb_ep(ifac), coefbf_ep(ifac), & 1453 qimp , hint ) 1454 1455 ! Dirichlet Boundary Condition 1456 !----------------------------- 1457 1458 else 1459 1460 pimp = pimp + cvar_ep(iel) 1461 pimp = pimp * cfnne 1462 call set_dirichlet_scalar & 1463 !==================== 1464 ( coefa_ep(ifac), coefaf_ep(ifac), & 1465 coefb_ep(ifac), coefbf_ep(ifac), & 1466 pimp , hint , rinfin ) 1467 1468 endif 1469 1470 elseif (iturb.eq.32) then 1471 1472 if(iwallf.ne.0) then 1473 ! Use k at I' 1474 xkip = 0.5d0*(rijipb(ifac,1)+rijipb(ifac,2)+rijipb(ifac,3)) 1475 1476 pimp_lam = 2.d0*visclc*xkip/(distbf**2*romc) 1477 1478 if (yplus.gt.epzero) then 1479 pimp_turb = 5.d0*uk**4*romc/ & 1480 (xkappa*visclc*(yplus+2.d0*dplus)) 1481 1482 ! Blending between wall and homogeneous layer 1483 ! from JF Wald PhD (2016) 1484 fep = exp(-((yplus+dplus)/4.d0)**1.5d0) 1485 dep = 1.d0- exp(-((yplus+dplus)/9.d0)**2.1d0) 1486 pimp = fep*pimp_lam + (1.d0-fep)*dep*pimp_turb 1487 else 1488 pimp = pimp_lam 1489 end if 1490 1491 else 1492 ! Use k at I' 1493 xkip = 0.5d0*(rijipb(ifac,1)+rijipb(ifac,2)+rijipb(ifac,3)) 1494 pimp = 2.d0*visclc*xkip/(distbf**2*romc) 1495 end if 1496 pimp = pimp * cfnne 1497 call set_dirichlet_scalar & 1498 !==================== 1499 ( coefa_ep(ifac), coefaf_ep(ifac), & 1500 coefb_ep(ifac), coefbf_ep(ifac), & 1501 pimp , hint , rinfin ) 1502 1503 ! ---> Alpha 1504 1505 ! Dirichlet Boundary Condition 1506 !----------------------------- 1507 if(iwallf.ne.0) then 1508 1509 if(yplus.gt.epzero) then 1510 ypsd = (yplus+dplus)/2.d0 1511 falpg = 16.d0 /(16.d0+4.d-2*ypsd)**2 * exp(- ypsd / (16.d0 + 4.d-2*ypsd) ) 1512 falpv = 1.d0 - exp(- (yplus+dplus) / (16.d0 + 4.d-2*(yplus+dplus)) ) 1513 pimp = falpv - (yplus+dplus)*falpg 1514 else 1515 pimp = 0.d0 1516 end if 1517 else 1518 pimp = 0.d0 1519 end if 1520 1521 hint = 1.d0/distbf 1522 pimp = pimp * cfnne 1523 call set_dirichlet_scalar & 1524 !==================== 1525 ( coefa_al(ifac), coefaf_al(ifac), & 1526 coefb_al(ifac), coefbf_al(ifac), & 1527 pimp , hint , rinfin ) 1528 1529 endif 1530 1531 !=========================================================================== 1532 ! 6a.Boundary conditions on k, epsilon, f_bar and phi in the phi_Fbar model 1533 !=========================================================================== 1534 1535 elseif (iturb.eq.50) then 1536 1537 ! Dirichlet Boundary Condition on k 1538 !---------------------------------- 1539 1540 pimp = 0.d0 1541 hint = (visclc+visctc/sigmak)/distbf 1542 1543 call set_dirichlet_scalar & 1544 !==================== 1545 ( coefa_k(ifac), coefaf_k(ifac), & 1546 coefb_k(ifac), coefbf_k(ifac), & 1547 pimp , hint , rinfin ) 1548 1549 ! Dirichlet Boundary Condition on epsilon 1550 !---------------------------------------- 1551 1552 pimp = 2.0d0*visclc/romc*cvar_k(iel)/distbf**2 1553 hint = (visclc+visctc/sigmae)/distbf 1554 1555 call set_dirichlet_scalar & 1556 !==================== 1557 ( coefa_ep(ifac), coefaf_ep(ifac), & 1558 coefb_ep(ifac), coefbf_ep(ifac), & 1559 pimp , hint , rinfin ) 1560 1561 ! Dirichlet Boundary Condition on Phi 1562 !------------------------------------ 1563 1564 pimp = 0.d0 1565 hint = (visclc+visctc/sigmak)/distbf 1566 1567 call set_dirichlet_scalar & 1568 !==================== 1569 ( coefa_phi(ifac), coefaf_phi(ifac), & 1570 coefb_phi(ifac), coefbf_phi(ifac), & 1571 pimp , hint , rinfin ) 1572 1573 ! Dirichlet Boundary Condition on Fb 1574 !----------------------------------- 1575 1576 pimp = 0.d0 1577 hint = 1.d0/distbf 1578 1579 call set_dirichlet_scalar & 1580 !==================== 1581 ( coefa_fb(ifac), coefaf_fb(ifac), & 1582 coefb_fb(ifac), coefbf_fb(ifac), & 1583 pimp , hint , rinfin ) 1584 1585 !=========================================================================== 1586 ! 6b.Boundary conditions on k, epsilon, phi and alpha in the Bl-v2/k model 1587 !=========================================================================== 1588 1589 elseif (iturb.eq.51) then 1590 1591 ! Dirichlet Boundary Condition on k 1592 !---------------------------------- 1593 1594 pimp = 0.d0 1595 hint = (visclc+visctc/sigmak)/distbf 1596 1597 call set_dirichlet_scalar & 1598 !==================== 1599 ( coefa_k(ifac), coefaf_k(ifac), & 1600 coefb_k(ifac), coefbf_k(ifac), & 1601 pimp , hint , rinfin ) 1602 1603 ! Dirichlet Boundary Condition on epsilon 1604 !---------------------------------------- 1605 1606 pimp = visclc/romc*cvar_k(iel)/distbf**2 1607 hint = (visclc+visctc/sigmae)/distbf 1608 1609 call set_dirichlet_scalar & 1610 !==================== 1611 ( coefa_ep(ifac), coefaf_ep(ifac), & 1612 coefb_ep(ifac), coefbf_ep(ifac), & 1613 pimp , hint , rinfin ) 1614 1615 ! Dirichlet Boundary Condition on Phi 1616 !------------------------------------ 1617 1618 pimp = 0.d0 1619 hint = (visclc+visctc/sigmak)/distbf 1620 1621 call set_dirichlet_scalar & 1622 !==================== 1623 ( coefa_phi(ifac), coefaf_phi(ifac), & 1624 coefb_phi(ifac), coefbf_phi(ifac), & 1625 pimp , hint , rinfin ) 1626 1627 ! Dirichlet Boundary Condition on alpha 1628 !-------------------------------------- 1629 1630 pimp = 0.d0 1631 hint = 1.d0/distbf 1632 1633 call set_dirichlet_scalar & 1634 !==================== 1635 ( coefa_al(ifac), coefaf_al(ifac), & 1636 coefb_al(ifac), coefbf_al(ifac), & 1637 pimp , hint , rinfin ) 1638 1639 1640 !=========================================================================== 1641 ! 7. Boundary conditions on k and omega 1642 !=========================================================================== 1643 1644 elseif (iturb.eq.60) then 1645 1646 ! Dirichlet Boundary Condition on k 1647 !---------------------------------- 1648 1649 ! Si on est hors de la sous-couche visqueuse (reellement ou via les 1650 ! scalable wall functions) 1651 if (iuntur.eq.1) then 1652 pimp = uk**2/sqrcmu 1653 1654 ! Si on est en sous-couche visqueuse 1655 else 1656 pimp = 0.d0 1657 endif 1658 1659 !FIXME it is wrong because sigma is computed within the model 1660 ! see turbkw.f90 1661 hint = (visclc+visctc/ckwsk2)/distbf 1662 1663 call set_dirichlet_scalar & 1664 !==================== 1665 ( coefa_k(ifac), coefaf_k(ifac), & 1666 coefb_k(ifac), coefbf_k(ifac), & 1667 pimp , hint , rinfin ) 1668 1669 ! Dirichlet Boundary Condition on omega 1670 !-------------------------------------- 1671 1672 !FIXME it is wrong because sigma is computed within the model 1673 ! see turbkw.f90 (So the flux is not the one we impose!) 1674 hint = (visclc+visctc/ckwsw2)/distbf 1675 1676 if (ikwcln.eq.1) then 1677 ! In viscous sub layer 1678 pimp_lam = 60.d0*visclc/(romc*ckwbt1*distbf**2) 1679 1680 ! If we are outside the viscous sub-layer (either naturally, or 1681 ! artificialy using scalable wall functions) 1682 1683 if (yplus.gt.epzero) then 1684 pimp_turb = 5.d0*uk**2*romc/ & 1685 (sqrcmu*xkappa*visclc*(yplus+2.d0*dplus)) 1686 1687 ! Use gamma function of Kader to weight 1688 !between high and low Reynolds meshes 1689 1690 gammap = -0.01d0*(yplus+2.d0*dplus)**4.d0/(1.d0+5.d0*(yplus+2.d0*dplus)) 1691 1692 pimp = pimp_lam*exp(gammap) + exp(1.d0/gammap)*pimp_turb 1693 else 1694 pimp = pimp_lam 1695 endif 1696 1697 call set_dirichlet_scalar & 1698 !==================== 1699 ( coefa_omg(ifac), coefaf_omg(ifac), & 1700 coefb_omg(ifac), coefbf_omg(ifac), & 1701 pimp , hint , rinfin ) 1702 1703 1704 ! If ikwcln is equal to 0, switch to deprecated Neumann 1705 ! condition on omega. 1706 else 1707 ! In viscous sub layer 1708 pimp_lam = 120.d0*8.d0*visclc/(romc*ckwbt1*distbf**2) 1709 1710 ! If we are outside the viscous sub-layer (either naturally, or 1711 ! artificialy using scalable wall functions) 1712 1713 if (yplus.gt.epzero) then 1714 pimp_turb = distbf*4.d0*uk**3*romc**2/ & 1715 (sqrcmu*xkappa*visclc**2*(yplus+2.d0*dplus)**2) 1716 1717 ! Use gamma function of Kader to weight 1718 !between high and low Reynolds meshes 1719 1720 gammap = -0.01d0*(yplus+2.d0*dplus)**4.d0 & 1721 / (1.d0+5.d0*(yplus+2.d0*dplus)) 1722 1723 pimp = pimp_lam*exp(gammap) + exp(1.d0/gammap)*pimp_turb 1724 else 1725 pimp = pimp_lam 1726 endif 1727 1728 qimp = -pimp*hint !TODO transform it, it is only to be fully equivalent 1729 1730 call set_neumann_scalar & 1731 !================== 1732 ( coefa_omg(ifac), coefaf_omg(ifac), & 1733 coefb_omg(ifac), coefbf_omg(ifac), & 1734 qimp , hint ) 1735 end if 1736 1737 !=========================================================================== 1738 ! 7.1 Boundary conditions on the Spalart Allmaras turbulence model 1739 !=========================================================================== 1740 1741 elseif (iturb.eq.70) then 1742 1743 ! Dirichlet Boundary Condition on nusa 1744 !------------------------------------- 1745 1746 pimp = 0.d0 1747 1748 hint = visclc / distbf / csasig ! Note: nusa is zero at the wall 1749 1750 call set_dirichlet_scalar & 1751 !==================== 1752 ( coefa_nu(ifac), coefaf_nu(ifac), & 1753 coefb_nu(ifac), coefbf_nu(ifac), & 1754 pimp , hint , rinfin ) 1755 1756 endif 1757 1758 byplus(ifac) = yplus 1759 bdplus(ifac) = dplus 1760 buk(ifac) = uk 1761 ustar(ifac) = uet 1762 bcfnns(ifac) = 1.d0 ! FIXME note taken into account yet in hturbp cfnns 1763 1764 endif 1765 ! Test on the presence of a smooth wall (End) 1766 1767enddo 1768! --- End of loop over faces 1769 1770 1771!=========================================================================== 1772! 8. Boundary conditions on the other scalars 1773! (Specific treatment for the variances of the scalars next to walls: 1774! see condli) 1775!=========================================================================== 1776 1777do iscal = 1, nscal 1778 1779 if (iscavr(iscal).le.0) then 1780 1781 f_id = ivarfl(isca(iscal)) 1782 1783 call field_get_dim(f_id, f_dim) 1784 1785 if (f_dim.eq.1) then 1786 1787 call clptur_scalar(iscal , isvhb , icodcl , rcodcl , & 1788 byplus , bdplus, buk , bcfnns , & 1789 hbord , theipb , & 1790 tetmax , tetmin, tplumx , tplumn ) 1791 1792 ! Vector field 1793 else 1794 1795 call clptur_vector(iscal, isvhb, icodcl, rcodcl, & 1796 byplus, bdplus, buk) 1797 1798 endif 1799 1800 endif 1801 1802enddo 1803 1804if (irangp.ge.0) then 1805 call parmin (uiptmn) 1806 call parmax (uiptmx) 1807 call parmin (uetmin) 1808 call parmax (uetmax) 1809 call parmin (ukmin) 1810 call parmax (ukmax) 1811 call parmin (yplumn) 1812 call parmax (yplumx) 1813 call parcpt (nlogla) 1814 call parcpt (nsubla) 1815 call parcpt (iuiptn) 1816 if (iscalt.gt.0) then 1817 call parmin (tetmin) 1818 call parmax (tetmax) 1819 call parmin (tplumn) 1820 call parmax (tplumx) 1821 endif 1822endif 1823 1824deallocate(byplus) 1825deallocate(buk) 1826if (allocated(buet)) deallocate(buet) 1827if (allocated(bcfnns_loc)) deallocate(bcfnns_loc) 1828 1829!=============================================================================== 1830! 9. Writings 1831!=============================================================================== 1832 1833! Remarque : afin de ne pas surcharger les logs dans le cas ou 1834! quelques yplus ne sont pas corrects, on ne produit le message 1835! qu'aux deux premiers pas de temps ou le message apparait et 1836! aux deux derniers pas de temps du calcul, ou si IWARNI est >= 2 1837! On indique aussi le numero du dernier pas de temps auquel on 1838! a rencontre des yplus hors bornes admissibles 1839 1840call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt_u) 1841 1842if (vcopt_u%iwarni.ge.0) then 1843 if (ntlist.gt.0) then 1844 modntl = mod(ntcabs,ntlist) 1845 elseif (ntlist.eq.-1.and.ntcabs.eq.ntmabs) then 1846 modntl = 0 1847 else 1848 modntl = 1 1849 endif 1850 1851 if ( (iturb.eq.0.and.nlogla.ne.0) .or. & 1852 (itytur.eq.5.and.nlogla.ne.0) .or. & 1853 ((itytur.eq.2.or.itytur.eq.3) .and. nsubla.gt.0) ) & 1854 ntlast = ntcabs 1855 1856 if ( (ntlast.eq.ntcabs.and.iaff.lt.2 ) .or. & 1857 (ntlast.ge.0 .and.ntcabs.ge.ntmabs-1) .or. & 1858 (ntlast.eq.ntcabs.and.vcopt_u%iwarni.ge.2) ) then 1859 iaff = iaff + 1 1860 1861 if (iscalt.gt.0) then 1862 write(nfecra,2011) & 1863 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1864 tetmin, tetmax, tplumn, tplumx, iuiptn,nsubla,nsubla+nlogla 1865 else 1866 write(nfecra,2010) & 1867 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1868 iuiptn,nsubla,nsubla+nlogla 1869 endif 1870 1871 if (iturb.eq. 0) write(nfecra,2020) ntlast,ypluli 1872 if (itytur.eq.5) write(nfecra,2030) ntlast,ypluli 1873 ! No warnings in EBRSM 1874 if ((itytur.eq.2.and.iturb.ne.22).or.iturb.eq.30.or.iturb.eq.31) & 1875 write(nfecra,2040) ntlast,ypluli 1876 if (vcopt_u%iwarni.lt.2.and.(iturb.ne.32.and.iturb.ne.22)) then 1877 write(nfecra,2050) 1878 elseif (iturb.ne.32.and.iturb.ne.22) then 1879 write(nfecra,2060) 1880 endif 1881 1882 else if (modntl.eq.0 .or. vcopt_u%iwarni.ge.2) then 1883 1884 if (iscalt.gt.0) then 1885 write(nfecra,2011) & 1886 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1887 tetmin, tetmax, tplumn, tplumx, iuiptn,nsubla,nsubla+nlogla 1888 else 1889 write(nfecra,2010) & 1890 uiptmn,uiptmx,uetmin,uetmax,ukmin,ukmax,yplumn,yplumx, & 1891 iuiptn,nsubla,nsubla+nlogla 1892 endif 1893 endif 1894 1895endif 1896 1897!=============================================================================== 1898! 10. Formats 1899!=============================================================================== 1900 1901 2010 format(/, & 1902 3X,'** BOUNDARY CONDITIONS FOR SMOOTH WALLS',/, & 1903 ' ---------------------------------------',/, & 1904 '------------------------------------------------------------',/,& 1905 ' Minimum Maximum',/,& 1906 '------------------------------------------------------------',/,& 1907 ' Rel velocity at the wall uiptn : ',2E12.5 ,/,& 1908 ' Friction velocity uet : ',2E12.5 ,/,& 1909 ' Friction velocity uk : ',2E12.5 ,/,& 1910 ' Dimensionless distance yplus : ',2E12.5 ,/,& 1911 ' ------------------------------------------------------' ,/,& 1912 ' Nb of reversal of the velocity at the wall : ',I10 ,/,& 1913 ' Nb of faces within the viscous sub-layer : ',I10 ,/,& 1914 ' Total number of wall faces : ',I10 ,/,& 1915 '------------------------------------------------------------', & 1916 /,/) 1917 1918 2011 format(/, & 1919 3X,'** BOUNDARY CONDITIONS FOR SMOOTH WALLS',/, & 1920 ' ---------------------------------------',/, & 1921 '------------------------------------------------------------',/,& 1922 ' Minimum Maximum',/,& 1923 '------------------------------------------------------------',/,& 1924 ' Rel velocity at the wall uiptn : ',2E12.5 ,/,& 1925 ' Friction velocity uet : ',2E12.5 ,/,& 1926 ' Friction velocity uk : ',2E12.5 ,/,& 1927 ' Dimensionless distance yplus : ',2E12.5 ,/,& 1928 ' Friction thermal sca. tstar : ',2E12.5 ,/,& 1929 ' Dim-less thermal sca. tplus : ',2E12.5 ,/,& 1930 ' ------------------------------------------------------' ,/,& 1931 ' Nb of reversal of the velocity at the wall : ',I10 ,/,& 1932 ' Nb of faces within the viscous sub-layer : ',I10 ,/,& 1933 ' Total number of wall faces : ',I10 ,/,& 1934 '------------------------------------------------------------', & 1935 /,/) 1936 1937 2020 format( & 1938'@ ',/,& 1939'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1940'@ ',/,& 1941'@ @@ WARNING: MESH NOT ENOUGH REFINED AT THE WALL ',/,& 1942'@ ======== ',/,& 1943'@ The mesh does not seem to be enough refined at the wall ',/,& 1944'@ to be able to run a laminar simulation. ',/,& 1945'@ ',/,& 1946'@ The last time step at which too large values for the ',/,& 1947'@ dimensionless distance to the wall (yplus) have been ',/,& 1948'@ observed is the time step ',I10 ,/,& 1949'@ ',/,& 1950'@ The minimum value for yplus must be lower than the ',/,& 1951'@ limit value YPLULI = ',E14.5 ,/,& 1952'@ ',/,& 1953'@ Have a look at the distribution of yplus at the wall ',/,& 1954'@ (with EnSight or ParaView for example) to conclude on ',/,& 1955'@ the way the results quality might be affected. ') 1956 1957 2030 format( & 1958'@ ',/,& 1959'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1960'@ ',/,& 1961'@ @@ WARNING: MESH NOT ENOUGH REFINED AT THE WALL ',/,& 1962'@ ======== ',/,& 1963'@ The mesh does not seem to be enough refined at the wall ',/,& 1964'@ to be able to run a v2f simulation ',/,& 1965'@ (phi-fbar or BL-v2/k) ',/,& 1966'@ ',/,& 1967'@ The last time step at which too large values for the ',/,& 1968'@ dimensionless distance to the wall (yplus) have been ',/,& 1969'@ observed is the time step ',I10 ,/,& 1970'@ ',/,& 1971'@ The minimum value for yplus must be lower than the ',/,& 1972'@ limit value YPLULI = ',E14.5 ,/,& 1973'@ ',/,& 1974'@ Have a look at the distribution of yplus at the wall ',/,& 1975'@ (with EnSight or ParaView for example) to conclude on ',/,& 1976'@ the way the results quality might be affected. ') 1977 1978 2040 format( & 1979'@ ',/,& 1980'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1981'@ ',/,& 1982'@ @@ WARNING: MESH TOO REFINED AT THE WALL ',/,& 1983'@ ======== ',/,& 1984'@ The mesh seems to be too refined at the wall to use ',/,& 1985'@ a high-Reynolds turbulence model. ',/,& 1986'@ ',/,& 1987'@ The last time step at which too small values for the ',/,& 1988'@ dimensionless distance to the wall (yplus) have been ',/,& 1989'@ observed is the time step ',I10 ,/,& 1990'@ ',/,& 1991'@ The minimum value for yplus must be greater than the ',/,& 1992'@ limit value YPLULI = ',E14.5 ,/,& 1993'@ ',/,& 1994'@ Have a look at the distribution of yplus at the wall ',/,& 1995'@ (with EnSight or ParaView for example) to conclude on ',/,& 1996'@ the way the results quality might be affected. ') 1997 2050 format( & 1998'@ ',/,& 1999'@ This warning is only printed at the first two ',/,& 2000'@ occurences of the problem and at the last time step ',/,& 2001'@ of the calculation. The vanishing of the message does ',/,& 2002'@ not necessarily mean the vanishing of the problem. ',/,& 2003'@ ',/,& 2004'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 2005'@ ',/) 2006 2060 format( & 2007'@ ',/,& 2008'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 2009'@ ',/) 2010 2011!---- 2012! End 2013!---- 2014 2015return 2016end subroutine 2017 2018!=============================================================================== 2019! Local functions 2020!=============================================================================== 2021 2022!------------------------------------------------------------------------------- 2023! Arguments 2024!______________________________________________________________________________. 2025! mode name role ! 2026!______________________________________________________________________________! 2027!> \param[in] iscal scalar id 2028!> \param[in] isvhb indicator to save exchange coeffient 2029!> \param[in,out] icodcl face boundary condition code: 2030!> - 1 Dirichlet 2031!> - 3 Neumann 2032!> - 4 sliding and 2033!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 2034!> - 5 smooth wall and 2035!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 2036!> - 6 rough wall and 2037!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 2038!> - 9 free inlet/outlet 2039!> (input mass flux blocked to 0) 2040!> \param[in,out] rcodcl boundary condition values: 2041!> - rcodcl(1) value of the dirichlet 2042!> - rcodcl(2) value of the exterior exchange 2043!> coefficient (infinite if no exchange) 2044!> - rcodcl(3) value flux density 2045!> (negative if gain) in w/m2 2046!> -# for the velocity \f$ (\mu+\mu_T) 2047!> \gradv \vect{u} \cdot \vect{n} \f$ 2048!> -# for the pressure \f$ \Delta t 2049!> \grad P \cdot \vect{n} \f$ 2050!> -# for a scalar \f$ cp \left( K + 2051!> \dfrac{K_T}{\sigma_T} \right) 2052!> \grad T \cdot \vect{n} \f$ 2053!> \param[in] byplus dimensionless distance to the wall 2054!> \param[in] bdplus dimensionless shift to the wall 2055!> for scalable wall functions 2056!> \param[in] buk dimensionless velocity 2057!> \param[in,out] hbord exchange coefficient at boundary 2058!> \param[in] theipb value of thermal scalar at \f$ \centip \f$ 2059!> of boundary cells 2060!> \param[out] tetmax maximum local ustar value 2061!> \param[out] tetmin minimum local ustar value 2062!> \param[out] tplumx maximum local tplus value 2063!> \param[out] tplumn minimum local tplus value 2064!_______________________________________________________________________________ 2065 2066subroutine clptur_scalar & 2067 ( iscal , isvhb , icodcl , rcodcl , & 2068 byplus , bdplus , buk , bcfnns , & 2069 hbord , theipb , & 2070 tetmax , tetmin , tplumx , tplumn ) 2071 2072!=============================================================================== 2073! Module files 2074!=============================================================================== 2075 2076use paramx 2077use numvar 2078use optcal 2079use cstphy 2080use cstnum 2081use pointe 2082use entsor 2083use albase 2084use parall 2085use ppppar 2086use ppthch 2087use ppincl 2088use radiat 2089use cplsat 2090use mesh 2091use field 2092use lagran 2093use turbomachinery 2094use cs_c_bindings 2095 2096!=============================================================================== 2097 2098implicit none 2099 2100! Arguments 2101 2102integer iscal, isvhb 2103integer, pointer, dimension(:,:) :: icodcl 2104 2105double precision, pointer, dimension(:,:,:) :: rcodcl 2106double precision, dimension(:) :: byplus, bdplus, buk, bcfnns 2107double precision, pointer, dimension(:) :: hbord, theipb 2108double precision tetmax, tetmin, tplumx, tplumn 2109 2110! Local variables 2111 2112integer ivar, f_id, b_f_id, isvhbl 2113integer f_id_ut, f_id_al 2114integer ifac, iel, isou, jsou 2115integer iscacp, ifcvsl, itplus, itstar 2116integer f_id_rough 2117integer kturt, turb_flux_model, turb_flux_model_type 2118 2119double precision cpp, rkl, prdtl, visclc, romc, tplus, cofimp, cpscv 2120double precision distfi, distbf, fikis, hint_al, heq, hflui, hext 2121double precision yplus, dplus, phit, pimp, pimp_al, rcprod, temp, tet, uk 2122double precision viscis, visctc, xmutlm, ypth, xnuii 2123double precision rinfiv(3), pimpv(3) 2124double precision visci(3,3), hintt(6) 2125double precision turb_schmidt, exchange_coef, visls_0 2126double precision mut_lm_dmut 2127double precision rough_t 2128 2129character(len=80) :: fname 2130 2131double precision, dimension(:), pointer :: val_s, bval_s, crom, viscls 2132double precision, dimension(:), pointer :: viscl, visct, cpro_cp, cpro_cv 2133 2134double precision, dimension(:), pointer :: bpro_rough_t 2135double precision, dimension(:), pointer :: bfconv, bhconv 2136double precision, dimension(:), pointer :: tplusp, tstarp, dist_theipb 2137double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp, hextp 2138double precision, dimension(:), pointer :: a_al, b_al, af_al, bf_al 2139double precision, dimension(:,:), pointer :: coefaut, cofafut, cofarut, visten 2140double precision, dimension(:,:,:), pointer :: coefbut, cofbfut, cofbrut 2141 2142double precision, allocatable, dimension(:) :: hbnd, hint, yptp 2143 2144integer, save :: kbfid = -1 2145 2146type(var_cal_opt) :: vcopt 2147 2148logical(c_bool), dimension(:), pointer :: cpl_faces 2149 2150!=============================================================================== 2151 2152ivar = isca(iscal) 2153f_id = ivarfl(ivar) 2154 2155call field_get_val_s(ivarfl(ivar), val_s) 2156 2157call field_get_val_s(iviscl, viscl) 2158call field_get_val_s(ivisct, visct) 2159 2160call field_get_key_int (f_id, kivisl, ifcvsl) 2161 2162if (ifcvsl .ge. 0) then 2163 call field_get_val_s(ifcvsl, viscls) 2164endif 2165 2166call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 2167 2168! If we have no diffusion, no boundary face should have a wall BC type 2169! (this is ensured in typecl) 2170 2171if (vcopt%idiff .eq. 0) then 2172 tetmax = 0.d0 2173 tetmin = 0.d0 2174 tplumx = 0.d0 2175 tplumn = 0.d0 2176 return 2177endif 2178 2179! Get the turbulent flux model for the scalar 2180call field_get_key_id('turbulent_flux_model', kturt) 2181call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 2182turb_flux_model_type = turb_flux_model / 10 2183 2184if ( iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0 & 2185 .or.turb_flux_model_type.eq.3) then 2186 if (iturb.ne.32.or.turb_flux_model_type.eq.3) then 2187 call field_get_val_v(ivsten, visten) 2188 else ! EBRSM and (GGDH or AFM) 2189 call field_get_val_v(ivstes, visten) 2190 endif 2191endif 2192 2193call field_get_coefa_s(f_id, coefap) 2194call field_get_coefb_s(f_id, coefbp) 2195call field_get_coefaf_s(f_id, cofafp) 2196call field_get_coefbf_s(f_id, cofbfp) 2197 2198call field_get_val_s(icrom, crom) 2199if (icp.ge.0) then 2200 call field_get_val_s(icp, cpro_cp) 2201endif 2202 2203if (ippmod(icompf) .ge. 0) then 2204 if (icv.ge.0) then 2205 call field_get_val_s(icv, cpro_cv) 2206 endif 2207endif 2208 2209isvhbl = 0 2210if (iscal.eq.isvhb) then 2211 isvhbl = isvhb 2212endif 2213 2214if (iscal.eq.iscalt) then 2215 ! min. and max. of wall friction of the thermal scalar 2216 tetmax = -grand 2217 tetmin = grand 2218 ! min. and max. of T+ 2219 tplumx = -grand 2220 tplumn = grand 2221endif 2222 2223rinfiv(1) = rinfin 2224rinfiv(2) = rinfin 2225rinfiv(3) = rinfin 2226 2227ypth = 0.d0 2228 2229if (turb_flux_model_type.eq.3) then 2230 2231 ! Turbulent diffusive flux of the scalar T 2232 ! (blending factor so that the component v'T' have only 2233 ! mu_T/(mu+mu_T)* Phi_T) 2234 2235 ! Name of the scalar ivar 2236 call field_get_name(ivarfl(ivar), fname) 2237 2238 ! Index of the corresponding turbulent flux 2239 call field_get_id(trim(fname)//'_turbulent_flux', f_id_ut) 2240 2241 call field_get_coefa_v(f_id_ut,coefaut) 2242 call field_get_coefb_v(f_id_ut,coefbut) 2243 call field_get_coefaf_v(f_id_ut,cofafut) 2244 call field_get_coefbf_v(f_id_ut,cofbfut) 2245 call field_get_coefad_v(f_id_ut,cofarut) 2246 call field_get_coefbd_v(f_id_ut,cofbrut) 2247 2248endif 2249 2250! EB-GGDH/AFM/DFM alpha boundary conditions 2251if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then 2252 2253 ! Name of the scalar ivar 2254 call field_get_name(ivarfl(ivar), fname) 2255 2256 ! Index of the corresponding turbulent flux 2257 call field_get_id(trim(fname)//'_alpha', f_id_al) 2258 2259 call field_get_coefa_s (f_id_al, a_al) 2260 call field_get_coefb_s (f_id_al, b_al) 2261 call field_get_coefaf_s(f_id_al, af_al) 2262 call field_get_coefbf_s(f_id_al, bf_al) 2263endif 2264 2265! pointers to T+ and T* if saved 2266 2267itplus = -1 2268itstar = -1 2269 2270tplusp => null() 2271tstarp => null() 2272 2273if (iscal.eq.iscalt) then 2274 call field_get_id_try('tplus', itplus) 2275 if (itplus.ge.0) then 2276 call field_get_val_s (itplus, tplusp) 2277 endif 2278 call field_get_id_try('tstar', itstar) 2279 if (itstar.ge.0) then 2280 call field_get_val_s (itstar, tstarp) 2281 endif 2282 if (vcopt%icoupl.gt.0) then 2283 allocate(dist_theipb(nfabor)) 2284 call cs_ic_field_dist_data_by_face_id(f_id, 1, theipb, dist_theipb) 2285 endif 2286endif 2287 2288bpro_rough_t => null() 2289 2290call field_get_id_try("boundary_roughness", f_id_rough) 2291if (f_id_rough.ge.0) then 2292 ! same thermal roughness if not specified 2293 call field_get_val_s(f_id_rough, bpro_rough_t) 2294endif 2295 2296call field_get_id_try("boundary_thermal_roughness", f_id_rough) 2297if (f_id_rough.ge.0) then 2298 call field_get_val_s(f_id_rough, bpro_rough_t) 2299endif 2300 2301if (vcopt%icoupl.gt.0) then 2302 call field_get_coupled_faces(f_id, cpl_faces) 2303endif 2304 2305allocate(hbnd(nfabor),hint(nfabor),yptp(nfabor)) 2306 2307! Pointers to specific fields 2308 2309if (iirayo.ge.1 .and. iscal.eq.iscalt) then 2310 call field_get_val_s_by_name("rad_convective_flux", bfconv) 2311 call field_get_val_s_by_name("rad_exchange_coefficient", bhconv) 2312endif 2313 2314! FIXME not really the BC value 2315if (kbfid.lt.0) call field_get_key_id("boundary_value_id", kbfid) 2316 2317call field_get_key_int(f_id, kbfid, b_f_id) 2318 2319! If thermal variable has no boundary but temperature does, use it 2320if (b_f_id .lt. 0 .and. iscal.eq.iscalt .and. itherm.eq.2) then 2321 b_f_id = itempb 2322endif 2323 2324if (b_f_id .ge. 0) then 2325 call field_get_val_s(b_f_id, bval_s) 2326else 2327 bval_s => null() 2328endif 2329 2330! Does the scalar behave as a temperature ? 2331call field_get_key_int(f_id, kscacp, iscacp) 2332 2333! Retrieve turbulent Schmidt value for current scalar 2334call field_get_key_double(f_id, ksigmas, turb_schmidt) 2335 2336! Reference diffusivity 2337call field_get_key_double(f_id, kvisl0, visls_0) 2338 2339! --- Loop on boundary faces 2340do ifac = 1, nfabor 2341 2342 yplus = byplus(ifac) 2343 dplus = bdplus(ifac) 2344 uk = buk(ifac) 2345 2346 ! Test on the presence of a smooth wall condition (start) 2347 if (icodcl(ifac,iu).eq.5) then 2348 2349 iel = ifabor(ifac) 2350 2351 ! Physical quantities 2352 2353 visclc = viscl(iel) 2354 visctc = visct(iel) 2355 romc = crom(iel) 2356 2357 xnuii = visclc / romc 2358 2359 ! Geometric quantities 2360 distbf = distb(ifac) 2361 2362 cpp = 1.d0 2363 if (iscacp.eq.1) then 2364 if (icp.ge.0) then 2365 cpp = cpro_cp(iel) 2366 else 2367 cpp = cp0 2368 endif 2369 endif 2370 2371 if (ifcvsl.lt.0) then 2372 rkl = visls_0 2373 prdtl = cpp*visclc/rkl 2374 else 2375 rkl = viscls(iel) 2376 prdtl = cpp*visclc/rkl 2377 endif 2378 2379 ! --> Compressible module: 2380 ! On suppose que le nombre de Pr doit etre 2381 ! defini de la meme facon que l'on resolve 2382 ! en enthalpie ou en energie, soit Mu*Cp/Lambda. 2383 ! Si l'on resout en energie, on a calcule ci-dessus 2384 ! Mu*Cv/Lambda. 2385 2386 if (iscal.eq.iscalt .and. itherm.eq.3) then 2387 if (icp.ge.0) then 2388 prdtl = prdtl*cpro_cp(iel) 2389 else 2390 prdtl = prdtl*cp0 2391 endif 2392 if (icv.ge.0) then 2393 prdtl = prdtl/cpro_cv(iel) 2394 else 2395 prdtl = prdtl/cv0 2396 endif 2397 endif 2398 2399 ! Scalar diffusivity 2400 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 2401 ! En compressible, pour l'energie LAMBDA/CV+CP/CV*(MUT/TURB_SCHMIDT) 2402 if (iscal.eq.iscalt .and. itherm.eq.3) then 2403 if (icp.ge.0) then 2404 cpscv = cpro_cp(iel) 2405 else 2406 cpscv = cp0 2407 endif 2408 if (icv.ge.0) then 2409 cpscv = cpscv/cpro_cv(iel) 2410 else 2411 cpscv = cpscv/cv0 2412 endif 2413 hint(ifac) = (rkl+vcopt%idifft*cpscv*visctc/turb_schmidt)/distbf 2414 else 2415 hint(ifac) = (rkl+vcopt%idifft*cpp*visctc/turb_schmidt)/distbf 2416 endif 2417 2418 ! Symmetric tensor diffusivity (GGDH or AFM) 2419 else if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 2420 ! En compressible, pour l'energie LAMBDA/CV+CP/CV*(MUT/SIGMAS) 2421 if (iscal.eq.iscalt .and. itherm.eq.3) then 2422 if (icp.ge.0) then 2423 cpscv = cpro_cp(iel) 2424 else 2425 cpscv = cp0 2426 endif 2427 if (icv.ge.0) then 2428 cpscv = cpscv/cpro_cv(iel) 2429 else 2430 cpscv = cpscv/cv0 2431 endif 2432 temp = vcopt%idifft*cpscv*ctheta(iscal)/csrij 2433 else 2434 temp = vcopt%idifft*cpp*ctheta(iscal)/csrij 2435 endif 2436 visci(1,1) = temp*visten(1,iel) + rkl 2437 visci(2,2) = temp*visten(2,iel) + rkl 2438 visci(3,3) = temp*visten(3,iel) + rkl 2439 visci(1,2) = temp*visten(4,iel) 2440 visci(2,1) = temp*visten(4,iel) 2441 visci(2,3) = temp*visten(5,iel) 2442 visci(3,2) = temp*visten(5,iel) 2443 visci(1,3) = temp*visten(6,iel) 2444 visci(3,1) = temp*visten(6,iel) 2445 2446 ! ||Ki.S||^2 2447 viscis = ( visci(1,1)*surfbo(1,ifac) & 2448 + visci(1,2)*surfbo(2,ifac) & 2449 + visci(1,3)*surfbo(3,ifac))**2 & 2450 + ( visci(2,1)*surfbo(1,ifac) & 2451 + visci(2,2)*surfbo(2,ifac) & 2452 + visci(2,3)*surfbo(3,ifac))**2 & 2453 + ( visci(3,1)*surfbo(1,ifac) & 2454 + visci(3,2)*surfbo(2,ifac) & 2455 + visci(3,3)*surfbo(3,ifac))**2 2456 2457 ! IF.Ki.S 2458 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 2459 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 2460 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 2461 )*surfbo(1,ifac) & 2462 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 2463 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 2464 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 2465 )*surfbo(2,ifac) & 2466 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 2467 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 2468 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 2469 )*surfbo(3,ifac) 2470 2471 distfi = distb(ifac) 2472 2473 ! Take I" so that I"F= eps*||FI||*Ki.n when I" is not in cell i 2474 ! NB: eps =1.d-1 must be consistent with vitens.f90 2475 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 2476 2477 hint(ifac) = viscis/surfbn(ifac)/fikis 2478 endif 2479 if (f_id_rough.ge.0) then 2480 rough_t = bpro_rough_t(ifac) 2481 else 2482 rough_t = 0.d0 2483 endif 2484 2485 ! wall function and Dirichlet or Neumann on the scalar 2486 if ( iturb.ne.0 .and. & 2487 (icodcl(ifac,ivar).eq.5 .or.icodcl(ifac,ivar).eq.6.or. & 2488 icodcl(ifac,ivar).eq.15 .or. & 2489 icodcl(ifac,ivar).eq.3)) then 2490 ! Note: to make things clearer yplus is always 2491 ! "y uk /nu" even for rough modelling. And the roughness correction is 2492 ! multiplied afterwards where needed. 2493 call hturbp(iwalfs,xnuii,prdtl,turb_schmidt,rough_t,uk,yplus,dplus,hflui,ypth) 2494 2495 ! Correction for non-neutral condition in atmospheric flows 2496 hflui = hflui * bcfnns(ifac) 2497 2498 ! Compute yk/T+ *PrT, take stability into account 2499 yptp(ifac) = hflui/prdtl 2500 ! Compute 2501 ! lambda/y * Pr_l *yk/T+ = lambda / nu * Pr_l * uk / T+ = rho cp uk / T+ 2502 ! so "Pr_l * yk/T+" is the correction factor compared to a laminar profile 2503 hflui = rkl/distbf *hflui 2504 2505 ! User exchange coefficient 2506 if (icodcl(ifac,ivar).eq.15) then 2507 hflui = rcodcl(ifac,ivar,2) 2508 yptp(ifac) = hflui/prdtl * distbf/rkl 2509 endif 2510 2511 else 2512 2513 ! y+/T+ *PrT 2514 yptp(ifac) = 1.d0/prdtl 2515 hflui = hint(ifac) 2516 2517 endif 2518 2519 hbnd(ifac) = hflui 2520 2521 endif ! smooth wall condition 2522 2523enddo 2524 2525! internal coupling 2526if (vcopt%icoupl.gt.0) then 2527 ! Update exchange coef. in coupling entity of current scalar 2528 call cs_ic_field_set_exchcoeff(f_id, hbnd) 2529 ! Get external exchange coef. 2530 call field_get_hext(f_id, hextp) 2531endif 2532 2533! --- Loop on boundary faces 2534do ifac = 1, nfabor 2535 2536 yplus = byplus(ifac) 2537 dplus = bdplus(ifac) 2538 uk = buk(ifac) 2539 2540 ! Test on the presence of a smooth wall condition (start) 2541 if (icodcl(ifac,iu).eq.5) then 2542 2543 iel = ifabor(ifac) 2544 2545 ! Physical quantities 2546 2547 visclc = viscl(iel) 2548 visctc = visct(iel) 2549 romc = crom(iel) 2550 2551 xnuii = visclc / romc 2552 2553 ! Geometric quantities 2554 distbf = distb(ifac) 2555 2556 cpp = 1.d0 2557 if (iscacp.eq.1) then 2558 if (icp.ge.0) then 2559 cpp = cpro_cp(iel) 2560 else 2561 cpp = cp0 2562 endif 2563 endif 2564 2565 if (ifcvsl.lt.0) then 2566 rkl = visls_0 2567 else 2568 rkl = viscls(iel) 2569 endif 2570 2571 hext = rcodcl(ifac,ivar,2) 2572 pimp = rcodcl(ifac,ivar,1) 2573 2574 if (vcopt%icoupl.gt.0) then 2575 if (cpl_faces(ifac)) then 2576 hext = hextp(ifac)/surfbn(ifac) 2577 endif 2578 endif 2579 2580 hflui = hbnd(ifac) 2581 2582 if (abs(hext).gt.rinfin*0.5d0.or.icodcl(ifac,ivar).eq.15) then 2583 heq = hflui 2584 if (vcopt%icoupl.gt.0) then 2585 ! ensure correct saving of flux in case of rad coupling 2586 if (cpl_faces(ifac)) then 2587 heq = hflui*hext/(hflui+hext) 2588 endif 2589 endif 2590 else 2591 heq = hflui*hext/(hflui+hext) 2592 endif 2593 2594 ! ---> Dirichlet Boundary condition with a wall function correction 2595 ! with or without an additional exchange coefficient hext 2596 2597 if (icodcl(ifac,ivar).eq.5.or.icodcl(ifac,ivar).eq.6 & 2598 .or.icodcl(ifac,ivar).eq.15) then 2599 ! DFM: the gradient BCs are so that the production term 2600 ! of u'T' is correcty computed 2601 if (turb_flux_model_type.ge.1) then 2602 ! In the log layer 2603 if (yplus.ge.ypth.and.iturb.ne.0) then 2604 xmutlm = xkappa*visclc*yplus 2605 if (cell_is_active(iel).eq.1) then 2606 mut_lm_dmut = xmutlm/max(visctc,1.e-12*visclc) 2607 else 2608 mut_lm_dmut = 0.d0 2609 endif 2610 2611 rcprod = min(xkappa , max(1.0d0,sqrt(mut_lm_dmut))/(yplus+dplus)) 2612 2613 cofimp = 1.d0 - yptp(ifac)*turb_schmidt/xkappa* & 2614 (2.0d0*rcprod - 1.0d0/(2.0d0*yplus+dplus)) 2615 2616 ! In the viscous sub-layer 2617 else 2618 cofimp = 0.d0 2619 endif 2620 else 2621 cofimp = 1.d0 - heq/hint(ifac) 2622 endif 2623 2624 ! To be coherent with a wall function, clip it to 0 2625 cofimp = max(cofimp, 0.d0) 2626 2627 ! Gradient BCs 2628 coefap(ifac) = (1.d0 - cofimp)*pimp 2629 coefbp(ifac) = cofimp 2630 2631 ! Flux BCs 2632 cofafp(ifac) = -heq*pimp 2633 cofbfp(ifac) = heq 2634 2635 ! Set coef for coupled face just to ensure relevant saving 2636 ! of bfconv if rad transfer activated 2637 if (vcopt%icoupl.gt.0) then 2638 if (cpl_faces(ifac)) then 2639 ! Flux BCs 2640 cofafp(ifac) = -heq*dist_theipb(ifac) 2641 cofbfp(ifac) = heq 2642 endif 2643 endif 2644 2645 ! Storage of the thermal exchange coefficient 2646 ! (conversion in case of energy or enthalpy) 2647 ! the exchange coefficient is in W/(m2 K) 2648 ! Useful for thermal coupling or radiative transfer 2649 if (iirayo.ge.1 .and. iscal.eq.iscalt.or.isvhbl.gt.0) then 2650 ! Enthalpy 2651 if (itherm.eq.2) then 2652 ! If Cp is variable 2653 if (icp.ge.0) then 2654 exchange_coef = hflui*cpro_cp(iel) 2655 else 2656 exchange_coef = hflui*cp0 2657 endif 2658 2659 ! Total energy (compressible module) 2660 elseif (itherm.eq.3) then 2661 ! If Cv is variable 2662 if (icv.ge.0) then 2663 exchange_coef = hflui*cpro_cv(iel) 2664 else 2665 exchange_coef = hflui*cv0 2666 endif 2667 2668 ! Temperature 2669 elseif (iscacp.eq.1) then 2670 exchange_coef = hflui 2671 endif 2672 endif 2673 2674 ! ---> Thermal coupling, store h = lambda/d 2675 if (isvhbl.gt.0) hbord(ifac) = exchange_coef 2676 2677 ! ---> Radiative transfer 2678 if (iirayo.ge.1 .and. iscal.eq.iscalt) then 2679 bhconv(ifac) = exchange_coef 2680 2681 ! The outgoing flux is stored (Q = h(Ti'-Tp): negative if 2682 ! gain for the fluid) in W/m2 2683 bfconv(ifac) = cofafp(ifac) + cofbfp(ifac)*theipb(ifac) 2684 endif 2685 2686 ! For the coupled faces with h_user (ie icodcl(ifac,ivar)=15) 2687 ! reset to zero af/bf coeff. 2688 ! By default icodcl(ifac,ivar)=3) for coupled faces 2689 if (vcopt%icoupl.gt.0) then 2690 if (cpl_faces(ifac)) then 2691 ! Flux BCs 2692 cofafp(ifac) = 0.d0 2693 cofbfp(ifac) = 0.d0 2694 endif 2695 endif 2696 2697 endif ! End if icodcl.eq.5 or 6 or 15 2698 2699 !--> Turbulent heat flux 2700 if (turb_flux_model_type.eq.3) then 2701 2702 ! Turbulent diffusive flux of the scalar T 2703 ! (blending factor so that the component v'T' have only 2704 ! mu_T/(mu+mu_T)* Phi_T) 2705 if (icodcl(ifac,ivar).eq.5.or.icodcl(ifac,ivar).eq.6 & 2706 .or.icodcl(ifac,ivar).eq.15) then 2707 phit = cofafp(ifac)+cofbfp(ifac)*val_s(iel) 2708 elseif (icodcl(ifac,ivar).eq.3) then 2709 phit = rcodcl(ifac,ivar,3) 2710 elseif (icodcl(ifac,ivar).eq.1) then 2711 phit = heq *(val_s(iel) - pimp) 2712 else 2713 phit = 0.d0 2714 endif 2715 2716 hintt(1) = 0.5d0*(visclc+rkl)/distbf & 2717 + visten(1,iel)*ctheta(iscal)/distbf/csrij 2718 hintt(2) = 0.5d0*(visclc+rkl)/distbf & 2719 + visten(2,iel)*ctheta(iscal)/distbf/csrij 2720 hintt(3) = 0.5d0*(visclc+rkl)/distbf & 2721 + visten(3,iel)*ctheta(iscal)/distbf/csrij 2722 hintt(4) = visten(4,iel)*ctheta(iscal)/distbf/csrij 2723 hintt(5) = visten(5,iel)*ctheta(iscal)/distbf/csrij 2724 hintt(6) = visten(6,iel)*ctheta(iscal)/distbf/csrij 2725 2726 ! Dirichlet Boundary Condition 2727 !----------------------------- 2728 2729 ! Add rho*uk*Tet to T'v' in High Reynolds 2730 if (yplus.ge.ypth) then 2731 do isou = 1, 3 2732 pimpv(isou) = surfbo(isou,ifac)*phit/(surfbn(ifac)*cpp*romc) 2733 enddo 2734 else 2735 do isou = 1, 3 2736 pimpv(isou) = 0.d0 2737 enddo 2738 endif 2739 2740 call set_dirichlet_vector_aniso & 2741 !======================== 2742 ( coefaut(:,ifac) , cofafut(:,ifac) , & 2743 coefbut(:,:,ifac), cofbfut(:,:,ifac), & 2744 pimpv , hintt , rinfiv ) 2745 2746 ! Boundary conditions used in the temperature equation 2747 do isou = 1, 3 2748 cofarut(isou,ifac) = 0.d0 2749 do jsou = 1, 3 2750 cofbrut(isou,jsou,ifac) = 0.d0 2751 enddo 2752 enddo 2753 2754 endif 2755 2756 ! EB-GGDH/AFM/DFM alpha boundary conditions 2757 if (turb_flux_model.eq.11 .or. & 2758 turb_flux_model.eq.21 .or. & 2759 turb_flux_model.eq.31) then 2760 2761 ! Dirichlet Boundary Condition 2762 !----------------------------- 2763 pimp_al = 0.d0 2764 hint_al = 1.d0/distbf 2765 2766 call set_dirichlet_scalar & 2767 !==================== 2768 ( a_al(ifac), af_al(ifac), & 2769 b_al(ifac), bf_al(ifac), & 2770 pimp_al , hint_al , rinfin ) 2771 2772 endif 2773 2774 ! Save the values of T^star and T^+ for post-processing 2775 2776 if (b_f_id.ge.0 .or. iscal.eq.iscalt) then 2777 2778 ! Wall function 2779 if (icodcl(ifac,ivar).eq.5.or.icodcl(ifac,ivar).eq.6 & 2780 .or.icodcl(ifac,ivar).eq.15) then 2781 if (iscal.eq.iscalt) then 2782 phit = cofafp(ifac)+cofbfp(ifac)*theipb(ifac) 2783 else 2784 phit = cofafp(ifac)+cofbfp(ifac)*bval_s(ifac) 2785 endif 2786 ! Imposed flux with wall function for post-processing 2787 elseif (icodcl(ifac,ivar).eq.3) then 2788 phit = rcodcl(ifac,ivar,3) ! = 0 if current face is coupled 2789 elseif (icodcl(ifac,ivar).eq.1) then 2790 if (iscal.eq.iscalt) then 2791 phit = heq *(theipb(ifac) - pimp) 2792 else 2793 phit = heq *(bval_s(ifac) - pimp) 2794 endif 2795 else 2796 phit = 0.d0 2797 endif 2798 2799 ! if face is coupled 2800 if (vcopt%icoupl.gt.0) then 2801 if (cpl_faces(ifac)) then 2802 phit = heq*(theipb(ifac)-dist_theipb(ifac)) 2803 endif 2804 endif 2805 2806 if (f_id_rough.ge.0) then 2807 rough_t = bpro_rough_t(ifac) 2808 else 2809 rough_t = 0.d0 2810 endif 2811 2812 tet = phit/(romc*cpp*max(uk, epzero)) 2813 ! T+ = (T_I - T_w) / Tet 2814 tplus = max(yplus, epzero)/yptp(ifac) 2815 2816 if (b_f_id.ge.0) bval_s(ifac) = bval_s(ifac) - tplus*tet 2817 2818 if (itplus.ge.0) tplusp(ifac) = tplus 2819 if (itstar.ge.0) tstarp(ifac) = tet 2820 2821 if (iscal.eq.iscalt) then 2822 tetmax = max(tet, tetmax) 2823 tetmin = min(tet, tetmin) 2824 tplumx = max(tplus,tplumx) 2825 tplumn = min(tplus,tplumn) 2826 endif 2827 2828 endif 2829 2830 endif ! smooth wall condition 2831 2832enddo 2833 2834deallocate(hbnd, hint, yptp) 2835 2836if (iscal.eq.iscalt.and.vcopt%icoupl.gt.0) then 2837 deallocate(dist_theipb) 2838endif 2839 2840!---- 2841! End 2842!---- 2843 2844return 2845end subroutine 2846 2847!=============================================================================== 2848! Local functions 2849!=============================================================================== 2850 2851!------------------------------------------------------------------------------- 2852! Arguments 2853!______________________________________________________________________________. 2854! mode name role ! 2855!______________________________________________________________________________! 2856!> \param[in] iscal scalar id 2857!> \param[in] isvhb indicator to save exchange coeffient 2858!> \param[in,out] icodcl face boundary condition code: 2859!> - 1 Dirichlet 2860!> - 3 Neumann 2861!> - 4 sliding and 2862!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 2863!> - 5 smooth wall and 2864!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 2865!> - 6 rough wall and 2866!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 2867!> - 9 free inlet/outlet 2868!> (input mass flux blocked to 0) 2869!> \param[in,out] rcodcl boundary condition values: 2870!> - rcodcl(1) value of the dirichlet 2871!> - rcodcl(2) value of the exterior exchange 2872!> coefficient (infinite if no exchange) 2873!> - rcodcl(3) value flux density 2874!> (negative if gain) in w/m2 2875!> -# for the velocity \f$ (\mu+\mu_T) 2876!> \gradv \vect{u} \cdot \vect{n} \f$ 2877!> -# for the pressure \f$ \Delta t 2878!> \grad P \cdot \vect{n} \f$ 2879!> -# for a scalar \f$ cp \left( K + 2880!> \dfrac{K_T}{\sigma_T} \right) 2881!> \grad T \cdot \vect{n} \f$ 2882!> \param[in] byplus dimensionless distance to the wall 2883!> \param[in] bdplus dimensionless shift to the wall 2884!> for scalable wall functions 2885!> \param[in] buk dimensionless velocity 2886!> \param[in,out] hbord exchange coefficient at boundary 2887!_______________________________________________________________________________ 2888 2889subroutine clptur_vector & 2890 ( iscal , isvhb , icodcl , & 2891 rcodcl , & 2892 byplus , bdplus , buk ) 2893 2894!=============================================================================== 2895! Module files 2896!=============================================================================== 2897 2898use paramx 2899use numvar 2900use optcal 2901use cstphy 2902use cstnum 2903use pointe 2904use entsor 2905use albase 2906use parall 2907use ppppar 2908use ppthch 2909use ppincl 2910use radiat 2911use cplsat 2912use mesh 2913use field 2914use field_operator 2915use lagran 2916use turbomachinery 2917use cs_c_bindings 2918use atincl 2919 2920!=============================================================================== 2921 2922implicit none 2923 2924! Arguments 2925 2926integer iscal, isvhb 2927integer, pointer, dimension(:,:) :: icodcl 2928 2929double precision, pointer, dimension(:,:,:) :: rcodcl 2930double precision, dimension(:) :: byplus, bdplus, buk 2931 2932! Local variables 2933 2934integer ivar, f_id, isvhbl 2935integer ifac, iel 2936integer iscacp, ifcvsl 2937integer f_id_rough 2938integer kturt, turb_flux_model, turb_flux_model_type 2939 2940double precision cpp, rkl, prdtl, visclc, romc, cofimp 2941double precision distbf, heq, yptp, hflui, hext 2942double precision yplus, dplus, rcprod, uk 2943double precision visctc, xmutlm, ypth, xnuii, srfbnf 2944double precision rcodcx, rcodcy, rcodcz, rcodcn, rnx, rny, rnz 2945double precision turb_schmidt, visls_0 2946double precision rough_t 2947 2948double precision, dimension(:), pointer :: bpro_rough_t 2949double precision, dimension(:), pointer :: crom, viscls 2950double precision, dimension(:,:), pointer :: val_p_v 2951double precision, dimension(:), pointer :: viscl, visct, cpro_cp, hextp 2952 2953double precision, dimension(:,:), pointer :: coefav, cofafv 2954double precision, dimension(:,:,:), pointer :: coefbv, cofbfv 2955 2956double precision, allocatable, dimension(:) :: hbnd, hint 2957 2958type(var_cal_opt) :: vcopt 2959 2960logical(c_bool), dimension(:), pointer :: cpl_faces 2961 2962!=============================================================================== 2963 2964ivar = isca(iscal) 2965f_id = ivarfl(ivar) 2966 2967call field_get_val_prev_v(f_id, val_p_v) 2968 2969call field_get_val_s(iviscl, viscl) 2970call field_get_val_s(ivisct, visct) 2971 2972call field_get_key_int (f_id, kivisl, ifcvsl) 2973 2974if (ifcvsl .ge. 0) then 2975 call field_get_val_s(ifcvsl, viscls) 2976endif 2977 2978call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 2979 2980call field_get_coefa_v(f_id, coefav) 2981call field_get_coefb_v(f_id, coefbv) 2982call field_get_coefaf_v(f_id, cofafv) 2983call field_get_coefbf_v(f_id, cofbfv) 2984 2985call field_get_val_s(icrom, crom) 2986if (icp.ge.0) then 2987 call field_get_val_s(icp, cpro_cp) 2988endif 2989 2990! Does the vector behave as a temperature ? 2991call field_get_key_int(f_id, kscacp, iscacp) 2992 2993! retrieve turbulent Schmidt value for current vector 2994call field_get_key_double(f_id, ksigmas, turb_schmidt) 2995 2996! Reference diffusivity 2997call field_get_key_double(f_id, kvisl0, visls_0) 2998 2999! Get the turbulent flux model 3000call field_get_key_id('turbulent_flux_model', kturt) 3001call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 3002turb_flux_model_type = turb_flux_model / 10 3003 3004isvhbl = 0 3005if (iscal.eq.isvhb) then 3006 isvhbl = isvhb 3007endif 3008 3009ypth = 0.d0 3010 3011bpro_rough_t => null() 3012 3013call field_get_id_try("boundary_roughness", f_id_rough) 3014if (f_id_rough.ge.0) then 3015 ! same thermal roughness if not specified 3016 call field_get_val_s(f_id_rough, bpro_rough_t) 3017endif 3018 3019call field_get_id_try("boundary_thermal_roughness", f_id_rough) 3020if (f_id_rough.ge.0) then 3021 call field_get_val_s(f_id_rough, bpro_rough_t) 3022endif 3023 3024if (vcopt%icoupl.gt.0) then 3025 call field_get_coupled_faces(f_id, cpl_faces) 3026endif 3027 3028allocate(hbnd(nfabor),hint(nfabor)) 3029 3030! --- Loop on boundary faces 3031do ifac = 1, nfabor 3032 3033 yplus = byplus(ifac) 3034 dplus = bdplus(ifac) 3035 uk = buk(ifac) 3036 3037 ! Geometric quantities 3038 distbf = distb(ifac) 3039 3040 ! Test on the presence of a smooth wall condition (start) 3041 if (icodcl(ifac,iu).eq.5) then 3042 3043 iel = ifabor(ifac) 3044 3045 ! Physical quantities 3046 3047 visclc = viscl(iel) 3048 visctc = visct(iel) 3049 romc = crom(iel) 3050 3051 xnuii = visclc / romc 3052 3053 cpp = 1.d0 3054 if (iscacp.eq.1) then 3055 if (icp.ge.0) then 3056 cpp = cpro_cp(iel) 3057 else 3058 cpp = cp0 3059 endif 3060 endif 3061 3062 if (ifcvsl.lt.0) then 3063 rkl = visls_0 3064 prdtl = cpp*visclc/rkl 3065 else 3066 rkl = viscls(iel) 3067 prdtl = cpp*visclc/rkl 3068 endif 3069 3070 ! Scalar diffusivity 3071 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 3072 hint(ifac) = (rkl+vcopt%idifft*cpp*visctc/turb_schmidt)/distbf 3073 else 3074 ! TODO if (vcopt%idften.eq.6) 3075 call csexit(1) 3076 endif 3077 3078 ! wall function and Dirichlet or Neumann on the scalar 3079 if (iturb.ne.0.and.(icodcl(ifac,ivar).eq.5.or.icodcl(ifac,ivar).eq.3)) then 3080 3081 if (f_id_rough.ge.0) then 3082 rough_t = bpro_rough_t(ifac) 3083 else 3084 rough_t = 0.d0 3085 endif 3086 3087 !FIXME use Re* = rough_t * uk / nu ? * PrT ? 3088 call hturbp(iwalfs,xnuii,prdtl,turb_schmidt,rough_t,uk,yplus,dplus,hflui,ypth) 3089 ! Compute (y+-d+)/T+ *PrT 3090 yptp = hflui/prdtl 3091 ! Compute lambda/y * (y+-d+)/T+ 3092 hflui = rkl/distbf *hflui 3093 3094 ! User exchange coefficient 3095 else if (icodcl(ifac,ivar).eq.15) then 3096 hflui = rcodcl(ifac,ivar,2) 3097 3098 else 3099 3100 ! y+/T+ *PrT 3101 yptp = 1.d0/prdtl 3102 hflui = hint(ifac) 3103 3104 endif 3105 3106 hbnd(ifac) = hflui 3107 3108 endif ! smooth wall condition 3109 3110enddo 3111 3112! internal coupling 3113if (vcopt%icoupl.gt.0) then 3114 ! Update exchange coef. in coupling entity of current scalar 3115 call cs_ic_field_set_exchcoeff(f_id, hbnd) 3116 ! Get external exchange coef. 3117 call field_get_hext(f_id, hextp) 3118endif 3119 3120! --- Loop on boundary faces 3121do ifac = 1, nfabor 3122 3123 yplus = byplus(ifac) 3124 dplus = bdplus(ifac) 3125 3126 ! Test on the presence of a smooth wall condition (start) 3127 if (icodcl(ifac,iu).eq.5) then 3128 3129 iel = ifabor(ifac) 3130 3131 srfbnf = surfbn(ifac) 3132 3133 ! Local framework 3134 3135 ! Unit normal 3136 3137 rnx = surfbo(1,ifac)/srfbnf 3138 rny = surfbo(2,ifac)/srfbnf 3139 rnz = surfbo(3,ifac)/srfbnf 3140 3141 ! Handle Dirichlet vector values 3142 3143 rcodcx = rcodcl(ifac,ivar ,1) 3144 rcodcy = rcodcl(ifac,ivar+1,1) 3145 rcodcz = rcodcl(ifac,ivar+2,1) 3146 3147 ! Keep tangential part 3148 3149 rcodcn = rcodcx*rnx+rcodcy*rny+rcodcz*rnz 3150 rcodcx = rcodcx -rcodcn*rnx 3151 rcodcy = rcodcy -rcodcn*rny 3152 rcodcz = rcodcz -rcodcn*rnz 3153 3154 ! Physical quantities 3155 3156 visclc = viscl(iel) 3157 visctc = visct(iel) 3158 3159 hext = rcodcl(ifac,ivar,2) 3160 3161 if (vcopt%icoupl.gt.0) then 3162 if (cpl_faces(ifac)) then 3163 hext = hextp(ifac)/surfbn(ifac) 3164 endif 3165 endif 3166 3167 hflui = hbnd(ifac) 3168 3169 if (abs(hext).gt.rinfin*0.5d0.or.icodcl(ifac,ivar).eq.15) then 3170 heq = hflui 3171 else 3172 heq = hflui*hext/(hflui+hext) 3173 endif 3174 3175 ! ---> Dirichlet Boundary condition with a wall function correction 3176 ! with or without an additional exchange coefficient hext 3177 3178 if (icodcl(ifac,ivar).eq.5.or.icodcl(ifac,ivar).eq.15) then 3179 ! DFM: the gradient BCs are so that the production term 3180 ! of u'T' is correcty computed 3181 if (turb_flux_model_type.ge.1) then 3182 ! In the log layer 3183 if (yplus.ge.ypth.and.iturb.ne.0) then 3184 xmutlm = xkappa*visclc*(yplus+dplus) 3185 rcprod = min(xkappa , max(1.0d0,sqrt(xmutlm/visctc))/(yplus+dplus)) 3186 3187 cofimp = 1.d0 - yptp*turb_schmidt/xkappa* & 3188 (2.0d0*rcprod - 1.0d0/(2.0d0*yplus+dplus)) 3189 3190 ! In the viscous sub-layer 3191 else 3192 cofimp = 0.d0 3193 endif 3194 else 3195 cofimp = 1.d0 - heq/hint(ifac) 3196 endif 3197 3198 ! To be coherent with a wall function, clip it to 0 3199 cofimp = max(cofimp, 0.d0) 3200 3201 ! Coupled solving of the velocity components 3202 3203 ! Gradient boundary conditions 3204 !----------------------------- 3205 rcodcn = rcodcx*rnx+rcodcy*rny+rcodcz*rnz 3206 3207 coefav(1,ifac) = (1.d0-cofimp)*(rcodcx - rcodcn*rnx) + rcodcn*rnx 3208 coefav(2,ifac) = (1.d0-cofimp)*(rcodcy - rcodcn*rny) + rcodcn*rny 3209 coefav(3,ifac) = (1.d0-cofimp)*(rcodcz - rcodcn*rnz) + rcodcn*rnz 3210 3211 ! Projection in order to have the vector parallel to the wall 3212 ! B = cofimp * ( IDENTITY - n x n ) 3213 3214 coefbv(1,1,ifac) = cofimp*(1.d0-rnx**2) 3215 coefbv(2,2,ifac) = cofimp*(1.d0-rny**2) 3216 coefbv(3,3,ifac) = cofimp*(1.d0-rnz**2) 3217 coefbv(1,2,ifac) = -cofimp*rnx*rny 3218 coefbv(1,3,ifac) = -cofimp*rnx*rnz 3219 coefbv(2,1,ifac) = -cofimp*rny*rnx 3220 coefbv(2,3,ifac) = -cofimp*rny*rnz 3221 coefbv(3,1,ifac) = -cofimp*rnz*rnx 3222 coefbv(3,2,ifac) = -cofimp*rnz*rny 3223 3224 ! Flux boundary conditions 3225 !------------------------- 3226 3227 cofafv(1,ifac) = -heq*(rcodcx - rcodcn*rnx) - hint(ifac)*rcodcn*rnx 3228 cofafv(2,ifac) = -heq*(rcodcy - rcodcn*rny) - hint(ifac)*rcodcn*rny 3229 cofafv(3,ifac) = -heq*(rcodcz - rcodcn*rnz) - hint(ifac)*rcodcn*rnz 3230 3231 ! Projection 3232 ! B = heq*( IDENTITY - n x n ) 3233 3234 cofbfv(1,1,ifac) = heq*(1.d0-rnx**2) + hint(ifac)*rnx**2 3235 cofbfv(2,2,ifac) = heq*(1.d0-rny**2) + hint(ifac)*rny**2 3236 cofbfv(3,3,ifac) = heq*(1.d0-rnz**2) + hint(ifac)*rnz**2 3237 3238 cofbfv(1,2,ifac) = (hint(ifac) - heq)*rnx*rny 3239 cofbfv(1,3,ifac) = (hint(ifac) - heq)*rnx*rnz 3240 cofbfv(2,1,ifac) = (hint(ifac) - heq)*rny*rnx 3241 cofbfv(2,3,ifac) = (hint(ifac) - heq)*rny*rnz 3242 cofbfv(3,1,ifac) = (hint(ifac) - heq)*rnz*rnx 3243 cofbfv(3,2,ifac) = (hint(ifac) - heq)*rnz*rny 3244 3245 endif ! End if icodcl.eq.5 3246 3247 ! TODO : postprocessing at the boundary 3248 3249 endif ! smooth wall condition 3250 3251enddo 3252 3253deallocate(hbnd,hint) 3254 3255!---- 3256! End 3257!---- 3258 3259return 3260end subroutine 3261