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!=============================================================================== 25! Function: 26! -------- 27!> \file visdyn.f90 28!> \brief Calculation of turbulent viscosity for 29!> a dynamic Smagorinsky LES model 30!> 31!> \f[ smago = \dfrac{L_{ij}M_{ij}}{M_{ij}M_{ij}} \f] 32!> 33!> \f[ \mu_T = \rho smago L^2 \sqrt{2 S_{ij}S_{ij}} \f] 34!> \f[ S_{ij} = \dfrac{\der{u_i}{x_j} + \der{u_j}{x_i}}{2}\f] 35!> 36!> We have at edge faces types at previous time step 37!> (except at first time step, when tables itypfb and itrifb 38!> have not been filled). 39!> 40!> Please refer to the 41!> <a href="../../theory.pdf#dynsmago"><b>dynamic Smagorinsky model</b></a> 42!> section of the theory guide for more informations. 43!------------------------------------------------------------------------------- 44 45!------------------------------------------------------------------------------- 46! Arguments 47!______________________________________________________________________________. 48! mode name role 49!______________________________________________________________________________! 50!> \param[in] nvar total number of variables 51!> \param[in] nscal total number of scalars 52!> \param[in] ncepdp number of cells with head loss 53!> \param[in] ncesmp number of cells with mass source term 54!> \param[in] icepdc number of ncepdp cells with losses 55!> \param[in] icetsm number of cells with mass source 56!> \param[in] itypsm type of mass source for the variable 57!> (cf. cs_user_mass_source_terms) 58!> \param[in] dt time step (per cell) 59!> \param[in] ckupdc work array for head losses 60!> \param[in] smacel value of variables associated to the 61!> mass source 62!> for ivar = ipr, smacel = mass flux 63!> \param[out] gradv the computed velocity gradients 64!______________________________________________________________________________! 65 66subroutine visdyn & 67 ( nvar , nscal , ncepdp , ncesmp , & 68 icepdc , icetsm , itypsm , & 69 dt , & 70 ckupdc , smacel, gradv ) 71 72!=============================================================================== 73! Module files 74!=============================================================================== 75 76use paramx 77use numvar 78use cstnum 79use optcal 80use cstphy 81use ppincl 82use entsor 83use parall 84use period 85use mesh 86use field 87use field_operator 88use cs_c_bindings 89 90!=============================================================================== 91 92implicit none 93 94! Arguments 95 96integer nvar , nscal 97integer ncepdp , ncesmp 98 99integer icepdc(ncepdp) 100integer icetsm(ncesmp), itypsm(ncesmp,nvar) 101 102double precision dt(ncelet) 103double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) 104double precision, intent(inout) :: gradv(3,3,ncelet) 105 106! Local variables 107 108integer ii, jj, iel, inc, ivar 109integer iprev 110integer iclipc 111integer iccocg, iclip 112integer key_turb_diff, key_sgs_sca_coef 113integer t_dif_id, sca_dync_id 114integer nswrgp, imligp, iwarnp 115 116double precision coef, radeux, deux, delta, deltaf 117double precision s11, s22, s33, s11f, s22f, s33f 118double precision dudy, dudz, dvdx, dvdz, dwdx, dwdy 119double precision dudyf, dudzf, dvdxf, dvdzf, dwdxf, dwdyf 120double precision xfil, xa, xb, xfil2, xsmgmx, xsmgmn 121double precision xl11, xl22, xl33, xl12, xl13, xl23 122double precision xm11, xm22, xm33, xm12, xm13, xm23 123double precision smagma, smagmi, smagmy 124double precision scal1, scal2, scal3 125double precision epsrgp, climgp 126 127double precision, allocatable, dimension(:) :: w1, w2, w3 128double precision, allocatable, dimension(:) :: w4, w5, w6 129double precision, allocatable, dimension(:) :: w7, w8, w9 130double precision, allocatable, dimension(:) :: s_n, sf_n 131double precision, allocatable, dimension(:) :: w0, xrof, xro 132double precision, allocatable, dimension(:,:) :: grads, scami, scamif 133double precision, allocatable, dimension(:,:) :: xmij, w61, w62, gradsf 134double precision, dimension(:,:), pointer :: coefau 135double precision, dimension(:,:,:), pointer :: coefbu 136double precision, dimension(:), pointer :: crom, coefas, coefbs, cvar_sca 137double precision, dimension(:), pointer :: cpro_turb_diff 138double precision, dimension(:,:), pointer :: vel 139double precision, dimension(:), pointer :: visct, cpro_smago, cpro_sca_dync 140 141type(var_cal_opt) :: vcopt 142 143!=============================================================================== 144 145!=============================================================================== 146! 1. Iniatilization 147!=============================================================================== 148 149! Map field arrays 150call field_get_val_v(ivarfl(iu), vel) 151 152call field_get_coefa_v(ivarfl(iu), coefau) 153call field_get_coefb_v(ivarfl(iu), coefbu) 154 155call field_get_val_s(ivisct, visct) 156call field_get_val_s(icrom, crom) 157 158call field_get_val_s(ismago, cpro_smago) 159 160! For the calculation of the viscosity of the sub-mesh 161xfil = xlesfl 162xfil2 = xlesfd 163xa = ales 164xb = bles 165deux = 2.d0 166radeux = sqrt(deux) 167xsmgmx = smagmx 168xsmgmn = smagmn 169 170! Allocate some work arrays 171 172allocate(w0(ncelet), w1(ncelet)) 173allocate(xmij(6,ncelet)) 174allocate(xro(ncelet), xrof(ncelet)) 175 176! Take into account variable density case: Favre filtering 177! Constant density case: Reynolds filtering 178if(irovar.eq.1) then 179 do iel = 1, ncel 180 xro(iel) = crom(iel) 181 enddo 182else 183 xro(:) = 1.d0 184endif 185 186! In case of constant density, xrof always 1.d0 187call les_filter(1, xro, xrof) 188!=============================================================================== 189! 2. Calculation of velocity gradient and of 190! S11**2+S22**2+S33**2+2*(S12**2+S13**2+S23**2) 191!=============================================================================== 192 193! Allocate temporary arrays for gradients calculation 194allocate(s_n(ncelet), sf_n(ncelet)) 195allocate(w61(6,ncelet), w62(6,ncelet)) 196 197inc = 1 198iprev = 0 199 200call field_gradient_vector(ivarfl(iu), iprev, 0, inc, gradv) 201 202do iel = 1, ncel 203 204 ! gradv(iel, xyz, uvw) 205 s11 = gradv(1, 1, iel) 206 s22 = gradv(2, 2, iel) 207 s33 = gradv(3, 3, iel) 208 dudy = gradv(2, 1, iel) 209 dudz = gradv(3, 1, iel) 210 dvdx = gradv(1, 2, iel) 211 dvdz = gradv(3, 2, iel) 212 dwdx = gradv(1, 3, iel) 213 dwdy = gradv(2, 3, iel) 214 215 216 217! In the case of constant density, s11+s22+s33 is zero 218 xmij(1,iel) = s11 - irovar*1.0d0/3.0d0*(s11+s22+s33) 219 xmij(2,iel) = s22 - irovar*1.0d0/3.0d0*(s11+s22+s33) 220 xmij(3,iel) = s33 - irovar*1.0d0/3.0d0*(s11+s22+s33) 221 xmij(4,iel) = 0.5d0*(dudy+dvdx) 222 xmij(5,iel) = 0.5d0*(dudz+dwdx) 223 xmij(6,iel) = 0.5d0*(dvdz+dwdy) 224 225 s_n(iel) = radeux*sqrt( & 226 xmij(1,iel)**2 + xmij(2,iel)**2 + xmij(3,iel)**2 & 227 + 2.d0*(xmij(4,iel)**2 + xmij(5,iel)**2 & 228 + xmij(6,iel)**2) ) 229 w62(:,iel) = xro(iel)*xmij(:,iel) 230enddo 231 232! w62 temperarily contains rho*S 233call les_filter(6, w62, w61) 234 235! w61 <rho*S>/<rho>, sf_n is ||<rho*S>/<rho>|| 236do iel = 1, ncel 237 w61(:,iel) = w61(:,iel)/xrof(iel) 238 sf_n(iel) = radeux*sqrt( & 239 w61(1,iel)**2 + w61(2,iel)**2 + w61(3,iel)**2 & 240 + 2.d0*(w61(4,iel)**2 + w61(5,iel)**2 & 241 + w61(6,iel)**2) ) 242enddo 243 244! Here XMIJ contains Sij 245! S_n contains ||S|| 246! sqrt(2)*sqrt(S11^2+S22^2+S33^2+2(S12^2+S13^2+S23^2)) 247! Sf_n contains ||SF|| 248! sqrt(2)*sqrt(S11F^2+S22F^2+S33F^2+2(S12F^2+S13F^2+S23F^2)) 249 250!=============================================================================== 251! 3. Calculation of Mij 252!=============================================================================== 253 254do iel = 1, ncel 255 w0(iel) = xfil *(xa*volume(iel))**xb 256enddo 257 258! Reuse xmij as temporary array 259 260do iel = 1, ncel 261 delta = w0(iel) 262 do ii = 1, 6 263 xmij(ii,iel) = -deux*xro(iel)*delta**2*s_n(iel)*xmij(ii,iel) 264 enddo 265enddo 266 267! w62 now contains <-2*rho*delta**2*||S||*S> 268call les_filter(6, xmij, w62) 269 270! Now compute final xmij value: M_ij = alpha_ij - beta_ij 271 272do iel = 1, ncel 273 delta = w0(iel) 274 deltaf = xfil2*delta 275 do ii = 1, 6 276 xmij(ii,iel) = -deux*xrof(iel)*deltaf**2*sf_n(iel)*w61(ii,iel) - w62(ii,iel) 277 enddo 278enddo 279 280deallocate(w61, w62) 281 282!=============================================================================== 283! 4. Calculation of the dynamic Smagorinsky constant 284!=============================================================================== 285 286! Allocate work arrays 287allocate(w2(ncelet), w3(ncelet), w4(ncelet)) 288allocate(w5(ncelet), w6(ncelet), w7(ncelet)) 289allocate(w8(ncelet), w9(ncelet)) 290 291! Filtering the velocity and its square 292 293! U**2 294do iel = 1,ncel 295 w0(iel) = xro(iel)*vel(1,iel)*vel(1,iel) 296enddo 297call les_filter(1, w0, w1) 298 299! V**2 300do iel = 1,ncel 301 w0(iel) = xro(iel)*vel(2,iel)*vel(2,iel) 302enddo 303call les_filter(1, w0, w2) 304 305! W**2 306do iel = 1,ncel 307 w0(iel) = xro(iel)*vel(3,iel)*vel(3,iel) 308enddo 309call les_filter(1, w0, w3) 310 311! UV 312do iel = 1,ncel 313 w0(iel) = xro(iel)*vel(1,iel)*vel(2,iel) 314enddo 315call les_filter(1, w0, w4) 316 317! UW 318do iel = 1,ncel 319 w0(iel) = xro(iel)*vel(1,iel)*vel(3,iel) 320enddo 321call les_filter(1, w0, w5) 322 323! VW 324do iel = 1,ncel 325 w0(iel) = xro(iel)*vel(2,iel)*vel(3,iel) 326enddo 327call les_filter(1, w0, w6) 328 329! U 330do iel = 1,ncel 331 w0(iel) = xro(iel)*vel(1,iel) 332enddo 333call les_filter(1, w0, w7) 334do iel = 1, ncel 335 w7(iel) = w7(iel)/xrof(iel) 336enddo 337 338! V 339do iel = 1,ncel 340 w0(iel) = xro(iel)*vel(2,iel) 341enddo 342call les_filter(1, w0, w8) 343do iel = 1, ncel 344 w8(iel) = w8(iel)/xrof(iel) 345enddo 346 347! W 348do iel = 1,ncel 349 w0(iel) = xro(iel)*vel(3,iel) 350enddo 351call les_filter(1, w0, w9) 352do iel = 1, ncel 353 w9(iel) = w9(iel)/xrof(iel) 354enddo 355 356do iel = 1, ncel 357 358 ! Calculation of Lij 359 xl11 = w1(iel) - xrof(iel) * w7(iel) * w7(iel) 360 xl22 = w2(iel) - xrof(iel) * w8(iel) * w8(iel) 361 xl33 = w3(iel) - xrof(iel) * w9(iel) * w9(iel) 362 xl12 = w4(iel) - xrof(iel) * w7(iel) * w8(iel) 363 xl13 = w5(iel) - xrof(iel) * w7(iel) * w9(iel) 364 xl23 = w6(iel) - xrof(iel) * w8(iel) * w9(iel) 365 366 xm11 = xmij(1,iel) 367 xm22 = xmij(2,iel) 368 xm33 = xmij(3,iel) 369 xm12 = xmij(4,iel) 370 xm13 = xmij(5,iel) 371 xm23 = xmij(6,iel) 372 ! Calculation of Mij :: Lij 373 w1(iel) = xm11 * xl11 + 2.d0* xm12 * xl12 + 2.d0* xm13 * xl13 & 374 + xm22 * xl22 + 2.d0* xm23 * xl23 & 375 + xm33 * xl33 376 ! Calculation of Mij :: Mij 377 w2(iel) = xm11 * xm11 + 2.d0* xm12 * xm12 + 2.d0* xm13 * xm13 & 378 + xm22 * xm22 + 2.d0* xm23 * xm23 & 379 + xm33 * xm33 380 381enddo 382 383deallocate(xmij) 384 385if (irangp.ge.0.or.iperio.eq.1) then 386 call synsca(w1) 387 call synsca(w2) 388endif 389 390! By default we make a local average of numerator and of 391! denominator, then only we make the quotient. 392! The user can do otherwise in ussmag. 393 394call les_filter(1, w1, w3) 395 396call les_filter(1, w2, w4) 397 398do iel = 1, ncel 399 if(abs(w4(iel)).le.epzero) then 400 cpro_smago(iel) = xsmgmx 401 else 402 cpro_smago(iel) = w3(iel)/w4(iel) 403 endif 404enddo 405 406call ussmag & 407 ( nvar , nscal , ncepdp , ncesmp , & 408 icepdc , icetsm , itypsm , & 409 dt , & 410 ckupdc , smacel , & 411 w1 , w2 ) 412 413iclipc = 0 414do iel = 1, ncel 415 if(cpro_smago(iel).ge.xsmgmx) then 416 cpro_smago(iel) = xsmgmx 417 iclipc = iclipc + 1 418 elseif(cpro_smago(iel).le.xsmgmn) then 419 cpro_smago(iel) = xsmgmn 420 iclipc = iclipc + 1 421 endif 422enddo 423 424!=============================================================================== 425! 5. Calculation of (dynamic) viscosity 426!=============================================================================== 427 428do iel = 1, ncel 429 coef = cpro_smago(iel) 430 delta = xfil * (xa*volume(iel))**xb 431 visct(iel) = crom(iel) * coef * delta**2 * s_n(iel) 432enddo 433 434call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 435 436! Some printings 437if (vcopt%iwarni.ge.1) then 438 439 smagma = -1.0d12 440 smagmi = 1.0d12 441 smagmy = 0.d0 442 do iel = 1, ncel 443 smagma = max(smagma,cpro_smago(iel)) 444 smagmi = min(smagmi,cpro_smago(iel)) 445 smagmy = smagmy + cpro_smago(iel)*volume(iel) 446 enddo 447 if (irangp.ge.0) then 448 call parmax(smagma) 449 call parmin(smagmi) 450 call parsom(smagmy) 451 call parcpt(iclipc) 452 endif 453 smagmy = smagmy / voltot 454 write(nfecra,1000) iclipc 455 write(nfecra,2001) 456 write(nfecra,2002) smagmy, smagmi, smagma 457 write(nfecra,2003) 458 459endif 460 461!=============================================================================== 462! 6. Clipping of the turbulent viscosity 463!=============================================================================== 464 465! Pour la LES en modele dynamique on clippe la viscosite turbulente de maniere 466! a ce que mu_t soit positif, 467 468iclipc = 0 469do iel = 1, ncel 470 if (visct(iel).lt.0.d0) then 471 visct(iel) = 0.d0 472 iclipc = iclipc + 1 473 endif 474enddo 475if (vcopt%iwarni.ge.1) then 476 if (irangp.ge.0) then 477 call parcpt(iclipc) 478 endif 479 write(nfecra,2004) iclipc 480endif 481 482!=============================================================================== 483! 7. Scalar turbulent model 484!=============================================================================== 485! In case of gaz combustion, the SGS scalar flux constant and the turbulent 486! diffusivity are only evaluated with the mixture fraction, then applied 487! automatically to the other scalar equations 488 489call field_get_key_id("turbulent_diffusivity_id", key_turb_diff) 490call field_get_key_id("sgs_scalar_flux_coef_id", key_sgs_sca_coef) 491 492do jj = 1, nscal 493 ivar = isca(jj) 494 495 ! For variance of a scalar, the turbulent diffsivity must not be computed 496 if (iscavr(jj).gt.0) cycle 497 498 ! For any scalar other than the mixture fraction in diffusion flames, 499 ! Dt is not computed either. TODO Soot may be an exception 500 if (ippmod(icod3p).ge.0) then 501 if (jj.ne.ifm) cycle 502 endif 503 504 call field_get_key_int(ivarfl(ivar), key_turb_diff, t_dif_id) 505 call field_get_key_int(ivarfl(ivar), key_sgs_sca_coef, sca_dync_id) 506 507 if (t_dif_id.ge.0.and.sca_dync_id.ge.0) then 508 call field_get_val_s(t_dif_id, cpro_turb_diff) 509 call field_get_val_s(sca_dync_id, cpro_sca_dync) 510 511 !================================================================ 512 ! 7.1. Compute the Mi for scalar 513 !================================================================ 514 515 allocate(grads(3, ncelet)) 516 allocate(gradsf(3, ncelet)) 517 518 inc = 1 519 iprev = 0 520 call field_get_coefa_s(ivarfl(ivar), coefas) 521 call field_get_coefb_s(ivarfl(ivar), coefbs) 522 523 call field_get_val_s(ivarfl(ivar),cvar_sca) 524 call field_gradient_scalar(ivarfl(ivar), iprev, imrgra, inc, 0, grads) 525 526 ! compute grad (<rho.Y>/<rho>) 527 allocate(scami(3,ncelet)) 528 allocate(scamif(3,ncelet)) 529 530 do iel = 1, ncel 531 w0(iel) = cvar_sca(iel)*xro(iel) 532 enddo 533 call les_filter(1, w0, w4) 534 do iel = 1, ncel 535 w4(iel) = w4(iel)/xrof(iel) 536 enddo 537 538 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 539 inc = 1 540 nswrgp = vcopt%nswrgr 541 epsrgp = vcopt%epsrgr 542 imligp = vcopt%imligr 543 iwarnp = vcopt%iwarni 544 climgp = vcopt%climgr 545 iccocg = 1 546 547 call gradient_s & 548 (-1 , imrgra , inc , iccocg , nswrgp , imligp , & 549 iwarnp , epsrgp , climgp , & 550 w4 , coefas , coefbs , & 551 gradsf ) 552 553 do iel = 1, ncel 554 delta = xfil * (xa*volume(iel))**xb 555 do ii = 1, 3 556 scami(ii,iel) = -xro(iel)*delta**2*s_n(iel)*grads(ii,iel) 557 enddo 558 enddo 559 560 call les_filter(3, scami, scamif) 561 do iel = 1, ncel 562 deltaf = xfil2*xfil * (xa*volume(iel))**xb 563 do ii = 1, 3 564 scami(ii,iel) = -deltaf**2*xrof(iel)*sf_n(iel)*gradsf(ii,iel) & 565 - scaMif(ii,iel) 566 enddo 567 enddo 568 569 !================================================================ 570 ! 7.2. Compute the Li for scalar 571 !================================================================ 572 ! rho*U*Y 573 do iel = 1, ncel 574 w0(iel) = xro(iel)*vel(1,iel)*cvar_sca(iel) 575 enddo 576 call les_filter(1, w0, w1) 577 578 ! rho*V*Y 579 do iel = 1, ncel 580 w0(iel) = xro(iel)*vel(2,iel)*cvar_sca(iel) 581 enddo 582 call les_filter(1, w0, w2) 583 584 ! rho*W*Y 585 do iel = 1, ncel 586 w0(iel) = xro(iel)*vel(3,iel)*cvar_sca(iel) 587 enddo 588 call les_filter(1, w0, w3) 589 590 do iel = 1, ncel 591 scal1 = w1(iel) - xrof(iel)*w7(iel)*w4(iel) 592 scal2 = w2(iel) - xrof(iel)*w8(iel)*w4(iel) 593 scal3 = w3(iel) - xrof(iel)*w9(iel)*w4(iel) 594 595 w1(iel) = scal1*scami(1,iel) + scal2*scami(2,iel) + scal3*scami(3,iel) 596 w2(iel) = scami(1,iel)**2 + scami(2,iel)**2 + scami(3,iel)**2 597 enddo 598 599 if (irangp.ge.0.or.iperio.eq.1) then 600 call synsca(w1) 601 call synsca(w2) 602 endif 603 604 call les_filter(1, w1, w3) 605 call les_filter(1, w2, w4) 606 607 !================================================================ 608 ! 7.3. Compute the SGS flux coefficient and SGS diffusivity 609 ! Cs >= 0, Dt >=0 610 !================================================================ 611 612 do iel = 1, ncel 613 if(abs(w4(iel)).le.epzero) then 614 cpro_sca_dync(iel) = 0.d0 615 else 616 cpro_sca_dync(iel) = max(w3(iel)/w4(iel), 0.d0) 617 endif 618 619 delta = xfil * (xa*volume(iel))**xb 620 cpro_turb_diff(iel) = crom(iel) * cpro_sca_dync(iel) * delta**2 * s_n(iel) 621 enddo 622 623 deallocate(scami, scamif) 624 deallocate(grads, gradsf) 625 endif 626 627enddo 628 629 630! Free memory 631deallocate(s_n, sf_n) 632deallocate(w9, w8) 633deallocate(w7, w6, w5) 634deallocate(w4, w3, w2) 635deallocate(w1, w0) 636deallocate(xro, xrof) 637 638!---- 639! Formats 640!---- 641 642 1000 format( & 643' Nb of clipping of the Smagorinsky constant', I10,/) 644 2001 format( & 645' --- Informations on the squared Smagorinsky constant' ,/,& 646' --------------------------------' ,/,& 647' Mean value Min value Max value' ,/,& 648' --------------------------------' ) 649 2002 format( & 650 e12.4 , e12.4, e12.4 ) 651 2003 format( & 652' --------------------------------' ,/) 653 2004 format( & 654' Nb of clippings for the turbulent viscosity (mu_t>0):', i10,/) 655 656!---- 657! End 658!---- 659 660return 661end subroutine 662