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 clsyvt.f90 28!> 29!> \brief Symmetry boundary conditions for vectors and tensors. 30!> 31!> Correspond to the code icodcl(ivar) = 4. 32!> 33!> Please refer to the 34!> <a href="../../theory.pdf#clsyvt"><b>clsyvt</b></a> section of the 35!> theory guide for more informations. 36!------------------------------------------------------------------------------- 37 38!------------------------------------------------------------------------------- 39! Arguments 40!______________________________________________________________________________. 41! mode name role ! 42!______________________________________________________________________________! 43!> \param[in] nscal total number of scalars 44!> \param[in,out] icodcl face boundary condition code: 45!> - 1 Dirichlet 46!> - 3 Neumann 47!> - 4 sliding and 48!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 49!> - 5 smooth wall and 50!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 51!> - 6 rough wall and 52!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 53!> - 9 free inlet/outlet 54!> (input mass flux blocked to 0) 55!> \param[in,out] rcodcl boundary condition values: 56!> - rcodcl(1) value of the dirichlet 57!> - rcodcl(2) value of the exterior exchange 58!> coefficient (infinite if no exchange) 59!> - rcodcl(3) value flux density 60!> (negative if gain) in w/m2 or roughness 61!> in m if icodcl=6 62!> -# for the velocity \f$ (\mu+\mu_T) 63!> \gradv \vect{u} \cdot \vect{n} \f$ 64!> -# for the pressure \f$ \Delta t 65!> \grad P \cdot \vect{n} \f$ 66!> -# for a scalar \f$ cp \left( K + 67!> \dfrac{K_T}{\sigma_T} \right) 68!> \grad T \cdot \vect{n} \f$ 69!> \param[in] velipb value of the velocity at \f$ \centip \f$ 70!> of boundary cells 71!> \param[in] rijipb value of \f$ R_{ij} \f$ at \f$ \centip \f$ 72!> of boundary cells 73!_______________________________________________________________________________ 74 75 76subroutine clsyvt & 77 ( nscal , icodcl , & 78 rcodcl , & 79 velipb , rijipb ) 80 81!=============================================================================== 82 83!=============================================================================== 84! Module files 85!=============================================================================== 86 87use paramx 88use numvar 89use optcal 90use cstphy 91use cstnum 92use pointe 93use entsor 94use albase 95use field 96use mesh 97use cs_c_bindings 98 99!=============================================================================== 100 101implicit none 102 103! Arguments 104 105integer nscal 106 107integer, pointer, dimension(:,:) :: icodcl 108 109double precision, pointer, dimension(:,:,:) :: rcodcl 110double precision, dimension(:,:) :: velipb 111double precision, pointer, dimension(:,:) :: rijipb 112 113! Local variables 114 115integer ifac, ii, isou 116integer iel, f_dim, clsyme 117integer iscal , ivar, idftnp 118integer kturt, turb_flux_model, turb_flux_model_type 119 120double precision rnx, rny, rnz 121double precision upx, upy, upz, usn 122double precision tx, ty, tz, txn, t2x, t2y, t2z 123double precision eloglo(3,3), alpha(6,6) 124double precision srfbnf, rcodcn, hint, visclc, visctc, distbf 125double precision cpp 126double precision hintv(6), rnn(6), htnn(6) 127double precision visci(3,3), fikis, viscis, distfi 128double precision fcoefa(6), fcoefb(6), fcofaf(6), fcofbf(6), fcofad(6), fcofbd(6) 129 130double precision, dimension(:,:), pointer :: coefau, cofafu, cfaale, claale 131double precision, dimension(:,:), pointer :: visten 132double precision, dimension(:,:,:), pointer :: coefbu, cofbfu, cfbale, clbale 133double precision, dimension(:,:), pointer :: coefa_rij, coefaf_rij, coefad_rij 134double precision, dimension(:,:,:), pointer :: coefb_rij, coefbf_rij, coefbd_rij 135 136double precision, dimension(:), pointer :: viscl, visct 137double precision, dimension(:), pointer :: cpro_visma_s 138double precision, dimension(:,:), pointer :: cpro_visma_v 139 140type(var_cal_opt) :: vcopt_uma, vcopt_rij 141 142!=============================================================================== 143 144!=============================================================================== 145! 1. Initializations 146!=============================================================================== 147 148cpp = 0.d0 149 150if (itytur.eq.3 .and. idirsm.eq.1) call field_get_val_v(ivsten, visten) 151 152call field_get_val_s(iviscl, viscl) 153call field_get_val_s(ivisct, visct) 154 155! Boundary Conditions 156 157call field_get_coefa_v(ivarfl(iu), coefau) 158call field_get_coefb_v(ivarfl(iu), coefbu) 159call field_get_coefaf_v(ivarfl(iu), cofafu) 160call field_get_coefbf_v(ivarfl(iu), cofbfu) 161 162if (itytur.eq.3) then 163 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt_rij) 164 call field_get_coefa_v(ivarfl(irij), coefa_rij) 165 call field_get_coefb_v(ivarfl(irij), coefb_rij) 166 call field_get_coefaf_v(ivarfl(irij), coefaf_rij) 167 call field_get_coefbf_v(ivarfl(irij), coefbf_rij) 168 call field_get_coefad_v(ivarfl(irij), coefad_rij) 169 call field_get_coefbd_v(ivarfl(irij), coefbd_rij) 170endif 171 172! Get the turbulent flux model key id 173call field_get_key_id('turbulent_flux_model', kturt) 174 175! --- Begin the loop over boundary faces 176do ifac = 1, nfabor 177 178 ! --- Test sur la presence d'une condition de symetrie vitesse : debut 179 if (icodcl(ifac,iu).eq.4) then 180 181 ! --- To cancel the mass flux 182 isympa(ifac) = 0 183 184 ! Geometric quantities 185 srfbnf = surfbn(ifac) 186 187 !=========================================================================== 188 ! 1. Local framework 189 !=========================================================================== 190 191 ! Unit normal 192 193 rnx = surfbo(1,ifac)/srfbnf 194 rny = surfbo(2,ifac)/srfbnf 195 rnz = surfbo(3,ifac)/srfbnf 196 197 ! En ALE, on a eventuellement une vitesse de deplacement de la face 198 ! donc seule la composante normale importe (on continue a determiner 199 ! TX a partir de la vitesse tangentielle absolue car l'orientation 200 ! de TX et T2X est sans importance pour les symetries) 201 rcodcn = 0.d0 202 if (iale.ge.1) then 203 rcodcn = rcodcl(ifac,iu,1)*rnx & 204 + rcodcl(ifac,iv,1)*rny & 205 + rcodcl(ifac,iw,1)*rnz 206 endif 207 208 upx = velipb(ifac,1) 209 upy = velipb(ifac,2) 210 upz = velipb(ifac,3) 211 212 if (itytur.eq.3) then 213 214 ! Relative tangential velocity 215 216 usn = upx*rnx+upy*rny+upz*rnz 217 tx = upx -usn*rnx 218 ty = upy -usn*rny 219 tz = upz -usn*rnz 220 txn = sqrt( tx**2 +ty**2 +tz**2 ) 221 222 ! Unit tangent 223 224 if( txn.ge.epzero) then 225 226 tx = tx/txn 227 ty = ty/txn 228 tz = tz/txn 229 230 else 231 232 ! If the velocity is zero, 233 ! Tx, Ty, Tz is not used (we cancel the velocity), so we assign any 234 ! value (zero for example) 235 236 tx = 0.d0 237 ty = 0.d0 238 tz = 0.d0 239 240 endif 241 242 ! --> T2 = RN X T (where X is the cross product) 243 244 t2x = rny*tz - rnz*ty 245 t2y = rnz*tx - rnx*tz 246 t2z = rnx*ty - rny*tx 247 248 ! --> Orthogonal matrix for change of reference frame ELOGLOij 249 ! (from local to global reference frame) 250 251 ! |TX -RNX T2X| 252 ! ELOGLO = |TY -RNY T2Y| 253 ! |TZ -RNZ T2Z| 254 255 ! Its transpose ELOGLOt is its inverse 256 257 eloglo(1,1) = tx 258 eloglo(1,2) = -rnx 259 eloglo(1,3) = t2x 260 eloglo(2,1) = ty 261 eloglo(2,2) = -rny 262 eloglo(2,3) = t2y 263 eloglo(3,1) = tz 264 eloglo(3,2) = -rnz 265 eloglo(3,3) = t2z 266 267 ! Compute Reynolds stress transformation matrix 268 269 clsyme = 1 270 call turbulence_bc_rij_transform(clsyme, eloglo, alpha) 271 272 endif 273 274 !=========================================================================== 275 ! 2. Boundary conditions on the velocity 276 ! (totaly implicit) 277 ! The condition is a zero (except in ALE) Dirichlet on the normal component 278 ! a homogenous Neumann on the other components 279 !=========================================================================== 280 281 ! Coupled solving of the velocity components 282 283 iel = ifabor(ifac) 284 ! --- Physical properties 285 visclc = viscl(iel) 286 visctc = visct(iel) 287 288 ! --- Geometrical quantity 289 distbf = distb(ifac) 290 291 if (itytur.eq.3) then 292 hint = visclc /distbf 293 else 294 hint = ( visclc+visctc )/distbf 295 endif 296 297 ! Gradient BCs 298 coefau(1,ifac) = rcodcn*rnx 299 coefau(2,ifac) = rcodcn*rny 300 coefau(3,ifac) = rcodcn*rnz 301 302 coefbu(1,1,ifac) = 1.d0-rnx**2 303 coefbu(2,2,ifac) = 1.d0-rny**2 304 coefbu(3,3,ifac) = 1.d0-rnz**2 305 306 coefbu(1,2,ifac) = -rnx*rny 307 coefbu(1,3,ifac) = -rnx*rnz 308 coefbu(2,1,ifac) = -rny*rnx 309 coefbu(2,3,ifac) = -rny*rnz 310 coefbu(3,1,ifac) = -rnz*rnx 311 coefbu(3,2,ifac) = -rnz*rny 312 313 ! Flux BCs 314 cofafu(1,ifac) = -hint*rcodcn*rnx 315 cofafu(2,ifac) = -hint*rcodcn*rny 316 cofafu(3,ifac) = -hint*rcodcn*rnz 317 318 cofbfu(1,1,ifac) = hint*rnx**2 319 cofbfu(2,2,ifac) = hint*rny**2 320 cofbfu(3,3,ifac) = hint*rnz**2 321 322 cofbfu(1,2,ifac) = hint*rnx*rny 323 cofbfu(1,3,ifac) = hint*rnx*rnz 324 cofbfu(2,1,ifac) = hint*rny*rnx 325 cofbfu(2,3,ifac) = hint*rny*rnz 326 cofbfu(3,1,ifac) = hint*rnz*rnx 327 cofbfu(3,2,ifac) = hint*rnz*rny 328 329 !=========================================================================== 330 ! 3. Boundary conditions on Rij (partially implicited) 331 !=========================================================================== 332 333 if (itytur.eq.3) then 334 335 ! Symmetric tensor diffusivity (Daly Harlow -- GGDH) 336 if (iand(vcopt_rij%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then 337 338 visci(1,1) = visclc + visten(1,iel) 339 visci(2,2) = visclc + visten(2,iel) 340 visci(3,3) = visclc + visten(3,iel) 341 visci(1,2) = visten(4,iel) 342 visci(2,1) = visten(4,iel) 343 visci(2,3) = visten(5,iel) 344 visci(3,2) = visten(5,iel) 345 visci(1,3) = visten(6,iel) 346 visci(3,1) = visten(6,iel) 347 348 ! ||Ki.S||^2 349 viscis = ( visci(1,1)*surfbo(1,ifac) & 350 + visci(1,2)*surfbo(2,ifac) & 351 + visci(1,3)*surfbo(3,ifac))**2 & 352 + ( visci(2,1)*surfbo(1,ifac) & 353 + visci(2,2)*surfbo(2,ifac) & 354 + visci(2,3)*surfbo(3,ifac))**2 & 355 + ( visci(3,1)*surfbo(1,ifac) & 356 + visci(3,2)*surfbo(2,ifac) & 357 + visci(3,3)*surfbo(3,ifac))**2 358 359 ! IF.Ki.S 360 fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1) & 361 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1) & 362 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1) & 363 )*surfbo(1,ifac) & 364 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2) & 365 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2) & 366 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2) & 367 )*surfbo(2,ifac) & 368 + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3) & 369 + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3) & 370 + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3) & 371 )*surfbo(3,ifac) 372 373 distfi = distb(ifac) 374 375 ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji 376 ! NB: eps =1.d-1 must be consistent with vitens.f90 377 fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi) 378 379 hint = viscis/surfbn(ifac)/fikis 380 381 ! Scalar diffusivity 382 else 383 hint = (visclc+visctc*csrij/cmu)/distbf 384 endif 385 386 do isou = 1, 6 387 fcoefa(isou) = 0.d0 388 fcofad(isou) = 0.d0 389 fcoefb(isou) = 0.d0 390 fcofbd(isou) = 0.d0 391 enddo 392 393 do isou = 1,6 394 395 if (isou.eq.1) then 396 ivar = ir11 397 elseif (isou.eq.2) then 398 ivar = ir22 399 elseif (isou.eq.3) then 400 ivar = ir33 401 elseif (isou.eq.4) then 402 ivar = ir12 403 elseif (isou.eq.5) then 404 ivar = ir23 405 elseif (isou.eq.6) then 406 ivar = ir13 407 endif 408 409 ! Partial (or total if coupled) implicitation 410 if (irijco.eq.1) then 411 coefa_rij(isou, ifac) = 0.d0 412 coefaf_rij(isou, ifac) = 0.d0 413 coefad_rij(isou, ifac) = 0.d0 414 do ii = 1, 6 415 coefb_rij(isou,ii, ifac) = alpha(ii, isou) 416 if (ii.eq.isou) then 417 coefbf_rij(isou,ii, ifac) = hint * (1.d0 - coefb_rij(isou,ii, ifac)) 418 else 419 coefbf_rij(isou,ii, ifac) = - hint * coefb_rij(isou,ii, ifac) 420 endif 421 coefbd_rij(isou,ii, ifac) = coefb_rij(isou,ii, ifac) 422 enddo 423 424 else if (iclsyr.eq.1) then 425 do ii = 1, 6 426 if (ii.ne.isou) then 427 fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii) 428 endif 429 enddo 430 fcoefb(isou) = alpha(isou,isou) 431 else 432 do ii = 1, 6 433 fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii) 434 enddo 435 fcoefb(isou) = 0.d0 436 endif 437 438 ! Boundary conditions for the momentum equation 439 fcofad(isou) = fcoefa(isou) 440 fcofbd(isou) = fcoefb(isou) 441 442 ! Translate into Diffusive flux BCs 443 fcofaf(isou) = -hint*fcoefa(isou) 444 fcofbf(isou) = hint*(1.d0-fcoefb(isou)) 445 enddo 446 447 if (irijco.ne.1) then 448 do isou = 1, 6 449 coefa_rij(isou,ifac) = fcoefa(isou) 450 coefaf_rij(isou,ifac) = fcofaf(isou) 451 coefad_rij(isou,ifac) = fcofad(isou) 452 do ii = 1,6 453 coefb_rij(isou,ii,ifac) = 0 454 coefbf_rij(isou,ii,ifac) = 0 455 coefbd_rij(isou,ii,ifac) = 0 456 enddo 457 coefb_rij(isou,isou,ifac) = fcoefb(isou) 458 coefbf_rij(isou,isou,ifac) = fcofbf(isou) 459 coefbd_rij(isou,isou,ifac) = fcofbd(isou) 460 enddo 461 endif 462 463 endif 464 465 endif 466! --- Test sur la presence d'une condition de symetrie vitesse : fin 467 468enddo 469! --- End of loop over boundary faces 470 471!=========================================================================== 472! 3.bis Boundary conditions on transported vectors 473!=========================================================================== 474 475do iscal = 1, nscal 476 477 ! Get the associated turbulent flux model 478 call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 479 turb_flux_model_type = turb_flux_model / 10 480 481 ! u'T' 482 if (turb_flux_model_type.eq.3) then 483 call clsyvt_scalar(iscal, icodcl) 484 endif 485 486 ! additional transported vectors 487 call field_get_dim(ivarfl(isca(iscal)), f_dim) 488 if (f_dim.gt.1) then 489 call clsyvt_vector(iscal, icodcl) 490 endif 491enddo 492 493!=============================================================================== 494! 4. Symmetry boundary conditions for mesh velocity (ALE module) 495!=============================================================================== 496 497if (iale.eq.1) then 498 499 call field_get_coefa_v(ivarfl(iuma), claale) 500 call field_get_coefb_v(ivarfl(iuma), clbale) 501 call field_get_coefaf_v(ivarfl(iuma), cfaale) 502 call field_get_coefbf_v(ivarfl(iuma), cfbale) 503 504 call field_get_key_struct_var_cal_opt(ivarfl(iuma), vcopt_uma) 505 idftnp = vcopt_uma%idften 506 507 if (iand(idftnp, ISOTROPIC_DIFFUSION).ne.0) then 508 call field_get_val_s(ivisma, cpro_visma_s) 509 else if (iand(idftnp, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then 510 call field_get_val_v(ivisma, cpro_visma_v) 511 endif 512 513 do ifac = 1, nfabor 514 if (icodcl(ifac,iuma).eq.4) then 515 516 iel = ifabor(ifac) 517 518 ! For a sliding boundary, the normal velocity is enforced to zero 519 ! whereas the other components have an Homogenous Neumann 520 ! NB: no recontruction in I' here 521 522 ! --- Geometrical quantity 523 distbf = distb(ifac) 524 srfbnf = surfbn(ifac) 525 526 ! --- Physical properties 527 if (iand(idftnp, ISOTROPIC_DIFFUSION).ne.0) then 528 do ii = 1, 3 529 hintv(ii) = cpro_visma_s(iel)/distbf 530 hintv(ii+3) = 0.d0 531 enddo 532 else if (iand(idftnp, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then 533 do ii = 1, 6 534 hintv(ii) = cpro_visma_v(ii,iel)/distbf 535 enddo 536 endif 537 538 ! Unit normal 539 rnx = surfbo(1,ifac)/srfbnf 540 rny = surfbo(2,ifac)/srfbnf 541 rnz = surfbo(3,ifac)/srfbnf 542 543 ! Coupled solving of the velocity components 544 545 ! Gradient BCs 546 claale(1,ifac) = 0.d0 547 claale(2,ifac) = 0.d0 548 claale(3,ifac) = 0.d0 549 550 clbale(1,1,ifac) = 1.d0-rnx**2 551 clbale(2,2,ifac) = 1.d0-rny**2 552 clbale(3,3,ifac) = 1.d0-rnz**2 553 554 clbale(1,2,ifac) = -rnx*rny 555 clbale(2,1,ifac) = -rny*rnx 556 clbale(1,3,ifac) = -rnx*rnz 557 clbale(3,1,ifac) = -rnz*rnx 558 clbale(2,3,ifac) = -rny*rnz 559 clbale(3,2,ifac) = -rnz*rny 560 561 ! Flux BCs 562 cfaale(1,ifac) = 0.d0 563 cfaale(2,ifac) = 0.d0 564 cfaale(3,ifac) = 0.d0 565 566 rnn(1) = rnx**2 567 rnn(2) = rny**2 568 rnn(3) = rnz**2 569 rnn(4) = rnx*rny 570 rnn(5) = rny*rnz 571 rnn(6) = rnx*rnz 572 573 call symmetric_matrix_product(hintv, rnn, htnn) 574 cfbale(1,1,ifac) = htnn(1) 575 cfbale(2,2,ifac) = htnn(2) 576 cfbale(3,3,ifac) = htnn(3) 577 cfbale(1,2,ifac) = htnn(4) 578 cfbale(2,1,ifac) = htnn(4) 579 cfbale(1,3,ifac) = htnn(6) 580 cfbale(3,1,ifac) = htnn(6) 581 cfbale(2,3,ifac) = htnn(5) 582 cfbale(3,2,ifac) = htnn(5) 583 584 endif 585 enddo 586 587endif 588 589!---- 590! End 591!---- 592 593return 594end subroutine 595 596!=============================================================================== 597! Local functions 598!=============================================================================== 599 600!------------------------------------------------------------------------------- 601! Arguments 602!______________________________________________________________________________. 603! mode name role ! 604!______________________________________________________________________________! 605!> \param[in] iscal scalar id 606!> \param[in,out] icodcl face boundary condition code: 607!> - 1 Dirichlet 608!> - 3 Neumann 609!> - 4 sliding and 610!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 611!> - 5 smooth wall and 612!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 613!> - 6 rough wall and 614!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 615!> - 9 free inlet/outlet 616!> (input mass flux blocked to 0) 617!_______________________________________________________________________________ 618 619subroutine clsyvt_scalar & 620 ( iscal , icodcl ) 621 622!=============================================================================== 623! Module files 624!=============================================================================== 625 626use paramx 627use numvar 628use optcal 629use cstphy 630use cstnum 631use dimens, only: nvar 632use pointe 633use entsor 634use albase 635use parall 636use ppppar 637use ppthch 638use ppincl 639use cplsat 640use mesh 641use field 642 643!=============================================================================== 644 645implicit none 646 647! Arguments 648 649integer iscal 650 651integer icodcl(nfabor,nvar) 652 653! Local variables 654 655integer ivar, f_id 656integer ifac, iel, isou, jsou 657integer iscacp, ifcvsl 658integer kturt, turb_flux_model, turb_flux_model_type 659 660double precision cpp, rkl, visclc, visls_0 661double precision distbf, srfbnf 662double precision rnx, rny, rnz 663double precision hintt(6) 664double precision hint, qimp 665 666character(len=80) :: fname 667 668double precision, dimension(:), pointer :: val_s, crom 669double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp 670double precision, dimension(:,:), pointer :: coefaut, cofafut, cofarut, visten 671double precision, dimension(:,:,:), pointer :: coefbut, cofbfut, cofbrut 672double precision, dimension(:), pointer :: a_al, b_al, af_al, bf_al 673 674double precision, dimension(:), pointer :: viscl, viscls, cpro_cp 675 676!=============================================================================== 677 678! Get the turbulent flux model for the scalar 679call field_get_key_id('turbulent_flux_model', kturt) 680call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 681turb_flux_model_type = turb_flux_model / 10 682 683if (turb_flux_model_type.ne.3) return 684 685ivar = isca(iscal) 686f_id = ivarfl(ivar) 687 688call field_get_val_s(ivarfl(ivar), val_s) 689call field_get_val_v(ivsten, visten) 690call field_get_val_s(iviscl, viscl) 691 692call field_get_coefa_s(f_id, coefap) 693call field_get_coefb_s(f_id, coefbp) 694call field_get_coefaf_s(f_id, cofafp) 695call field_get_coefbf_s(f_id, cofbfp) 696 697call field_get_val_s(icrom, crom) 698 699if (icp.ge.0) then 700 call field_get_val_s(icp, cpro_cp) 701endif 702 703! Reference diffusivity of the primal scalar 704call field_get_key_double(f_id, kvisl0, visls_0) 705 706call field_get_key_int (f_id, kivisl, ifcvsl) 707if (ifcvsl .ge. 0) then 708 call field_get_val_s(ifcvsl, viscls) 709endif 710 711! Does the scalar behave as a temperature ? 712call field_get_key_int(f_id, kscacp, iscacp) 713 714! Turbulent diffusive flux of the scalar T 715! (blending factor so that the component v'T' have only 716! mu_T/(mu+mu_T)* Phi_T) 717 718! Name of the scalar ivar 719call field_get_name(ivarfl(ivar), fname) 720 721! Index of the corresponding turbulent flux 722call field_get_id(trim(fname)//'_turbulent_flux', f_id) 723 724call field_get_coefa_v(f_id,coefaut) 725call field_get_coefb_v(f_id,coefbut) 726call field_get_coefaf_v(f_id,cofafut) 727call field_get_coefbf_v(f_id,cofbfut) 728call field_get_coefad_v(f_id,cofarut) 729call field_get_coefbd_v(f_id,cofbrut) 730 731! EB-GGDH/AFM/DFM alpha boundary conditions 732if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then 733 734 ! Name of the scalar ivar 735 call field_get_name(ivarfl(ivar), fname) 736 737 ! Index of the corresponding turbulent flux 738 call field_get_id(trim(fname)//'_alpha', f_id) 739 740 call field_get_coefa_s (f_id, a_al) 741 call field_get_coefb_s (f_id, b_al) 742 call field_get_coefaf_s(f_id, af_al) 743 call field_get_coefbf_s(f_id, bf_al) 744endif 745 746! --- Loop on boundary faces 747do ifac = 1, nfabor 748 749 ! Test on symmetry boundary condition: start 750 if (icodcl(ifac,iu).eq.4) then 751 752 iel = ifabor(ifac) 753 ! --- Physical Properties 754 visclc = viscl(iel) 755 cpp = 1.d0 756 if (iscacp.eq.1) then 757 if (icp.ge.0) then 758 cpp = cpro_cp(iel) 759 else 760 cpp = cp0 761 endif 762 endif 763 764 ! --- Geometrical quantities 765 distbf = distb(ifac) 766 srfbnf = surfbn(ifac) 767 768 rnx = surfbo(1,ifac)/srfbnf 769 rny = surfbo(2,ifac)/srfbnf 770 rnz = surfbo(3,ifac)/srfbnf 771 772 if (ifcvsl .lt. 0) then 773 rkl = visls_0/cpp 774 else 775 rkl = viscls(iel)/cpp 776 endif 777 778 !FIXME for EB DFM 779 do isou = 1, 6 780 if (isou.le.3) then 781 hintt(isou) = (0.5d0*(visclc + rkl) & 782 + ctheta(iscal)*visten(isou,iel)/csrij) / distbf 783 else 784 hintt(isou) = ctheta(iscal)*visten(isou,iel) / csrij / distbf 785 endif 786 enddo 787 788 ! Gradient BCs 789 coefaut(1,ifac) = 0.d0 790 coefaut(2,ifac) = 0.d0 791 coefaut(3,ifac) = 0.d0 792 793 coefbut(1,1,ifac) = 1.d0-rnx**2 794 coefbut(2,2,ifac) = 1.d0-rny**2 795 coefbut(3,3,ifac) = 1.d0-rnz**2 796 797 coefbut(1,2,ifac) = -rnx*rny 798 coefbut(1,3,ifac) = -rnx*rnz 799 coefbut(2,1,ifac) = -rny*rnx 800 coefbut(2,3,ifac) = -rny*rnz 801 coefbut(3,1,ifac) = -rnz*rnx 802 coefbut(3,2,ifac) = -rnz*rny 803 804 ! Flux BCs 805 cofafut(1,ifac) = 0.d0 806 cofafut(2,ifac) = 0.d0 807 cofafut(3,ifac) = 0.d0 808 809 cofbfut(1,1,ifac) = hintt(1)*rnx**2 + hintt(4)*rnx*rny + hintt(6)*rnx*rnz 810 cofbfut(2,2,ifac) = hintt(4)*rnx*rny + hintt(2)*rny**2 + hintt(5)*rny*rnz 811 cofbfut(3,3,ifac) = hintt(6)*rnx*rnz + hintt(5)*rny*rnz + hintt(3)*rnz**2 812 813 cofbfut(1,2,ifac) = hintt(1)*rnx*rny + hintt(4)*rny**2 + hintt(6)*rny*rnz 814 cofbfut(2,1,ifac) = hintt(4)*rnx**2 + hintt(2)*rny*rnx + hintt(5)*rnx*rnz 815 cofbfut(1,3,ifac) = hintt(1)*rnx*rnz + hintt(4)*rny*rnz + hintt(6)*rnz**2 816 cofbfut(3,1,ifac) = hintt(6)*rnx**2 + hintt(5)*rny*rnx + hintt(3)*rnz*rnx 817 cofbfut(2,3,ifac) = hintt(4)*rnx*rnz + hintt(2)*rny*rnz + hintt(5)*rnz**2 818 cofbfut(3,2,ifac) = hintt(6)*rnx*rny + hintt(5)*rny**2 + hintt(3)*rnz*rny 819 820 ! Boundary conditions for thermal transport equation 821 do isou = 1, 3 822 cofarut(isou,ifac) = coefaut(isou,ifac) 823 do jsou =1, 3 824 cofbrut(isou,jsou,ifac) = coefbut(isou,jsou,ifac) 825 enddo 826 enddo 827 828 ! EB-GGDH/AFM/DFM alpha boundary conditions 829 if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then 830 831 ! Dirichlet Boundary Condition 832 !----------------------------- 833 834 qimp = 0.d0 835 836 hint = 1.d0/distbf 837 838 call set_neumann_scalar & 839 !==================== 840 ( a_al(ifac), af_al(ifac), & 841 b_al(ifac), bf_al(ifac), & 842 qimp , hint ) 843 844 endif 845 846 847 endif 848 ! Test on velocity symmetry condition: end 849 850enddo 851 852!---- 853! End 854!---- 855 856return 857end subroutine 858 859!=============================================================================== 860! Local functions 861!=============================================================================== 862 863!------------------------------------------------------------------------------- 864! Arguments 865!______________________________________________________________________________. 866! mode name role ! 867!______________________________________________________________________________! 868!> \param[in] iscal scalar id 869!> \param[in,out] icodcl face boundary condition code: 870!> - 1 Dirichlet 871!> - 3 Neumann 872!> - 4 sliding and 873!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 874!> - 5 smooth wall and 875!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 876!> - 6 rough wall and 877!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 878!> - 9 free inlet/outlet 879!> (input mass flux blocked to 0) 880!_______________________________________________________________________________ 881 882subroutine clsyvt_vector & 883 ( iscal , icodcl ) 884 885!=============================================================================== 886! Module files 887!=============================================================================== 888 889use paramx 890use numvar 891use optcal 892use cstphy 893use cstnum 894use dimens, only: nvar 895use pointe 896use entsor 897use albase 898use parall 899use ppppar 900use ppthch 901use ppincl 902use cplsat 903use mesh 904use field 905use cs_c_bindings 906 907!=============================================================================== 908 909implicit none 910 911! Arguments 912 913integer iscal 914 915integer icodcl(nfabor,nvar) 916 917! Local variables 918 919integer ivar, f_id 920integer ifac, iel 921integer ifcvsl 922 923double precision rkl, visclc, visls_0 924double precision distbf, srfbnf 925double precision rnx, rny, rnz, temp 926double precision hintt(6) 927double precision turb_schmidt 928 929double precision, dimension(:), pointer :: crom 930double precision, dimension(:,:), pointer :: coefav, cofafv 931double precision, dimension(:,:,:), pointer :: coefbv, cofbfv 932double precision, dimension(:,:), pointer :: visten 933 934double precision, dimension(:), pointer :: viscl, visct, viscls 935 936type(var_cal_opt) :: vcopt 937 938!=============================================================================== 939 940ivar = isca(iscal) 941f_id = ivarfl(ivar) 942 943call field_get_key_struct_var_cal_opt(f_id, vcopt) 944 945if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 946 if (iturb.ne.32) then 947 call field_get_val_v(ivsten, visten) 948 else ! EBRSM and (GGDH or AFM) 949 call field_get_val_v(ivstes, visten) 950 endif 951endif 952 953call field_get_val_s(iviscl, viscl) 954call field_get_val_s(ivisct, visct) 955 956call field_get_coefa_v(f_id, coefav) 957call field_get_coefb_v(f_id, coefbv) 958call field_get_coefaf_v(f_id, cofafv) 959call field_get_coefbf_v(f_id, cofbfv) 960 961call field_get_val_s(icrom, crom) 962 963call field_get_key_int (f_id, kivisl, ifcvsl) 964if (ifcvsl .ge. 0) then 965 call field_get_val_s(ifcvsl, viscls) 966endif 967 968! retrieve turbulent Schmidt value for current scalar 969call field_get_key_double(ivarfl(ivar), ksigmas, turb_schmidt) 970 971! Reference diffusivity 972call field_get_key_double(f_id, kvisl0, visls_0) 973 974! --- Loop on boundary faces 975do ifac = 1, nfabor 976 977 ! Test on symmetry boundary condition: start 978 if (icodcl(ifac,ivar).eq.4) then 979 980 iel = ifabor(ifac) 981 ! --- Physical Properties 982 visclc = viscl(iel) 983 984 ! --- Geometrical quantities 985 distbf = distb(ifac) 986 srfbnf = surfbn(ifac) 987 988 rnx = surfbo(1,ifac)/srfbnf 989 rny = surfbo(2,ifac)/srfbnf 990 rnz = surfbo(3,ifac)/srfbnf 991 992 if (ifcvsl .lt. 0) then 993 rkl = visls_0 994 else 995 rkl = viscls(iel) 996 endif 997 998 ! Isotropic diffusivity 999 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 1000 hintt(1) = (vcopt%idifft*max(visct(iel),zero)/turb_schmidt + rkl)/distbf 1001 hintt(2) = hintt(1) 1002 hintt(3) = hintt(1) 1003 hintt(4) = 0.d0 1004 hintt(5) = 0.d0 1005 hintt(6) = 0.d0 1006 ! Symmetric tensor diffusivity 1007 elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 1008 temp = vcopt%idifft*ctheta(iscal)/csrij 1009 hintt(1) = (temp*visten(1,iel) + rkl)/distbf 1010 hintt(2) = (temp*visten(2,iel) + rkl)/distbf 1011 hintt(3) = (temp*visten(3,iel) + rkl)/distbf 1012 hintt(4) = temp*visten(4,iel) /distbf 1013 hintt(5) = temp*visten(5,iel) /distbf 1014 hintt(6) = temp*visten(6,iel) /distbf 1015 endif 1016 1017 ! Gradient BCs 1018 coefav(1,ifac) = 0.d0 1019 coefav(2,ifac) = 0.d0 1020 coefav(3,ifac) = 0.d0 1021 1022 coefbv(1,1,ifac) = 1.d0-rnx**2 1023 coefbv(2,2,ifac) = 1.d0-rny**2 1024 coefbv(3,3,ifac) = 1.d0-rnz**2 1025 1026 coefbv(1,2,ifac) = -rnx*rny 1027 coefbv(1,3,ifac) = -rnx*rnz 1028 coefbv(2,1,ifac) = -rny*rnx 1029 coefbv(2,3,ifac) = -rny*rnz 1030 coefbv(3,1,ifac) = -rnz*rnx 1031 coefbv(3,2,ifac) = -rnz*rny 1032 1033 ! Flux BCs 1034 cofafv(1,ifac) = 0.d0 1035 cofafv(2,ifac) = 0.d0 1036 cofafv(3,ifac) = 0.d0 1037 1038 cofbfv(1,1,ifac) = hintt(1)*rnx**2 + hintt(4)*rnx*rny + hintt(6)*rnx*rnz 1039 cofbfv(2,2,ifac) = hintt(4)*rnx*rny + hintt(2)*rny**2 + hintt(5)*rny*rnz 1040 cofbfv(3,3,ifac) = hintt(6)*rnx*rnz + hintt(5)*rny*rnz + hintt(3)*rnz**2 1041 1042 cofbfv(1,2,ifac) = hintt(1)*rnx*rny + hintt(4)*rny**2 + hintt(6)*rny*rnz 1043 cofbfv(2,1,ifac) = hintt(1)*rnx*rny + hintt(4)*rny**2 + hintt(6)*rny*rnz 1044 cofbfv(1,3,ifac) = hintt(1)*rnx*rnz + hintt(4)*rny*rnz + hintt(6)*rnz**2 1045 cofbfv(3,1,ifac) = hintt(1)*rnx*rnz + hintt(4)*rny*rnz + hintt(6)*rnz**2 1046 cofbfv(2,3,ifac) = hintt(4)*rnx*rnz + hintt(2)*rny*rnz + hintt(5)*rnz**2 1047 cofbfv(3,2,ifac) = hintt(4)*rnx*rnz + hintt(2)*rny*rnz + hintt(5)*rnz**2 1048 1049 endif 1050 ! Test on velocity symmetry condition: end 1051 1052enddo 1053 1054!---- 1055! End 1056!---- 1057 1058return 1059end subroutine 1060