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 d3ptcl.f90 28!> \brief Automatic boundary conditions for 3 PTHEM gas diffusion flame model 29 30!------------------------------------------------------------------------------- 31! Arguments 32!______________________________________________________________________________. 33! mode name role ! 34!______________________________________________________________________________! 35!> \param[in] nvar total number of variables 36!> \param[in] itypfb boundary face types 37!> \param[out] izfppp boundary face zone number 38!> \param[out] icodcl boundary condition code: 39!> - 1 Dirichlet 40!> - 2 Radiative outlet 41!> - 3 Neumann 42!> - 4 sliding and 43!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 44!> - 5 smooth wall and 45!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 46!> - 6 rough wall and 47!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 48!> - 9 free inlet/outlet 49!> (input mass flux blocked to 0) 50!> \param[out] rcodcl boundary condition values: 51!> - rcodcl(1) value of the dirichlet 52!> - rcodcl(2) value of the exterior exchange 53!> coefficient (infinite if no exchange) 54!> - rcodcl(3) value flux density 55!> (negative if gain) in w/m2 or roughness 56!> in m if icodcl=6 57!> -# for the velocity \f$ (\mu+\mu_T) 58!> \gradt \, \vect{u} \cdot \vect{n} \f$ 59!> -# for the pressure \f$ \Delta t 60!> \grad P \cdot \vect{n} \f$ 61!> -# for a scalar \f$ cp \left( K + 62!> \dfrac{K_T}{\sigma_T} \right) 63!> \grad T \cdot \vect{n} \f$ 64!_______________________________________________________________________________ 65 66subroutine d3ptcl & 67 ( itypfb , izfppp , icodcl, & 68 rcodcl ) 69 70!=============================================================================== 71 72!=============================================================================== 73! Module files 74!=============================================================================== 75 76use paramx 77use dimens, only : nscal, nvar 78use numvar 79use optcal 80use cstphy 81use cstnum 82use entsor 83use parall 84use ppppar 85use ppthch 86use coincl 87use cpincl 88use ppincl 89use mesh 90use field 91use cs_c_bindings 92 93!=============================================================================== 94 95implicit none 96 97! Arguments 98 99integer itypfb(nfabor) 100integer izfppp(nfabor) 101integer icodcl(nfabor,nvar) 102 103double precision rcodcl(nfabor,nvar,3) 104 105! Local variables 106 107integer igg, ifac, izone, mode, iscal 108integer ii, iel, ifue, ioxy 109integer icke 110double precision qimabs, qisqc, viscla, d2s3, uref2, rhomoy, dhy, xiturb 111double precision ustar2, xkent, xeent 112double precision qcalc(nozppm) 113double precision coefg(ngazgm) 114double precision, dimension(:), pointer :: brom 115double precision, dimension(:), pointer :: viscl 116 117!=============================================================================== 118! 0. INITIALISATIONS 119!=============================================================================== 120 121call field_get_val_s(ibrom, brom) 122call field_get_val_s(iviscl, viscl) 123 124d2s3 = 2.d0/3.d0 125 126do igg = 1, ngazgm 127 coefg(igg) = zero 128enddo 129 130!=============================================================================== 131! 1. ECHANGES EN PARALLELE POUR LES DONNEES UTILISATEUR 132!=============================================================================== 133 134! En realite on pourrait eviter cet echange en modifiant usd3pc et en 135! demandant a l'utilisateur de donner les grandeurs dependant de la 136! zone hors de la boucle sur les faces de bord : les grandeurs 137! seraient ainsi disponibles sur tous les processeurs. Cependant, 138! ca rend le sous programme utilisateur un peu plus complique et 139! surtout, si l'utilisateur le modifie de travers, ca ne marche pas. 140! On suppose que toutes les gandeurs fournies sont positives, ce qui 141! permet d'utiliser un max pour que tous les procs les connaissent. 142! Si ce n'est pas le cas, c'est plus complique mais on pourrait 143! s'en tirer avec un max quand meme, sauf qimp qui peut etre negatif. 144 145if(irangp.ge.0) then 146 call parmax(tinfue) 147 call parmax(tinoxy) 148 ! qimabs = |qimp|, or 0 if no boundary faces on the rank 149 do ii = 1, nozapm 150 qimabs = abs(qimp(ii)) 151 ! the rank owning the max of qimabs, meaning |qimp|, share 152 ! qimp with the others 153 call parmxl(1,qimabs,qimp(ii)) 154 enddo 155 call parimx(nozapm,iqimp ) 156 call parimx(nozapm,ientox) 157 call parimx(nozapm,ientfu) 158endif 159 160!=============================================================================== 161! 2. REMPLISSAGE DU TABLEAU DES CONDITIONS LIMITES POUR LES SCALAIRES 162! 163! ON BOUCLE SUR TOUTES LES FACES D'ENTREE 164! ========================= 165 166! ON DETERMINE LA FAMILLE ET SES PROPRIETES 167! ON IMPOSE LES CONDITIONS AUX LIMITES 168!=============================================================================== 169 170! Fuel and Oxidant inlet are checked 171ifue = 0 172ioxy = 0 173do ii = 1, nzfppp 174 izone = ilzppp(ii) 175 if (ientfu(izone).eq.1) then 176 ifue = 1 177 elseif(ientox(izone).eq.1) then 178 ioxy = 1 179 endif 180enddo 181if(irangp.ge.0) then 182 call parcmx(ifue) 183 call parcmx(ioxy) 184endif 185 186! Fuel inlet at TINFUE: HINFUE calculation 187if (ifue.eq.1) then 188 coefg(1) = 1.d0 189 coefg(2) = zero 190 coefg(3) = zero 191 mode = -1 192 call cothht & 193 !========== 194 ( mode , ngazg , ngazgm , coefg , & 195 npo , npot , th , ehgazg , & 196 hinfue , tinfue ) 197endif 198 199! Oxidant inlet at TINOXY: HINOXY calculation 200if (ioxy.eq.1) then 201 coefg(1) = zero 202 coefg(2) = 1.d0 203 coefg(3) = zero 204 mode = -1 205 call cothht & 206 !========== 207 ( mode , ngazg , ngazgm , coefg , & 208 npo , npot , th , ehgazg , & 209 hinoxy , tinoxy ) 210endif 211 212 213do ifac = 1, nfabor 214 215 izone = izfppp(ifac) 216 217 if (itypfb(ifac).eq.ientre.or.itypfb(ifac).eq.i_convective_inlet) then 218 219 ! Fuel inlet at TINFUE 220 if (ientfu(izone).eq.1) then 221 222 ! Neumann for outflow 223 if (qimp(izone).lt.0.d0) then 224 225 do iscal = 1, nscal 226 icodcl(ifac,isca(iscal)) = 3 227 rcodcl(ifac,isca(iscal),3) = 0.d0 228 enddo 229 230 ! Dirichlet for inflow 231 else 232 ! Mean mixture fraction 233 rcodcl(ifac,isca(ifm),1) = 1.d0 234 235 ! Mixture fraction variance 236 rcodcl(ifac,isca(ifp2m),1) = 0.d0 237 238 ! Mixture enthalpy 239 if (ippmod(icod3p).eq.1) then 240 rcodcl(ifac,isca(iscalt),1) = hinfue 241 endif 242 243 ! Soot model 244 if (isoot.ge.1) then 245 rcodcl(ifac,isca(ifsm),1) = 0.d0 246 rcodcl(ifac,isca(inpm),1) = 0.d0 247 endif 248 249 endif 250 251 ! Density 252 brom(ifac) = pther/(cs_physical_constants_r*tinfue/wmolg(1)) 253 254 ! Oxydant inlet at TINOXY 255 elseif (ientox(izone).eq.1) then 256 257 ! Neumann for outflow 258 if (qimp(izone).lt.0.d0) then 259 260 do iscal = 1, nscal 261 icodcl(ifac,isca(iscal)) = 3 262 rcodcl(ifac,isca(iscal),3) = 0.d0 263 enddo 264 265 ! Dirichlet for inflow 266 else 267 268 ! Mean mixture fraction 269 rcodcl(ifac,isca(ifm),1) = 0.d0 270 271 ! Mixture fraction variance 272 rcodcl(ifac,isca(ifp2m),1) = 0.d0 273 274 ! Mixture enthalpy 275 if ( ippmod(icod3p).eq.1 ) then 276 rcodcl(ifac,isca(iscalt),1) = hinoxy 277 endif 278 279 ! Soot model 280 if (isoot.ge.1) then 281 rcodcl(ifac,isca(ifsm),1) = 0.d0 282 rcodcl(ifac,isca(inpm),1) = 0.d0 283 endif 284 285 ! Density 286 brom(ifac) = pther/(cs_physical_constants_r*tinoxy/wmolg(2)) 287 288 endif 289 290 endif 291 292 endif 293 294enddo 295 296!=============================================================================== 297! 3. SI IQIMP = 1 : CORRECTION DES VITESSES (EN NORME) POUR CONTROLER 298! LES DEBITS IMPOSES 299! SI IQIMP = 0 : CALCUL DE QIMP 300 301! ON BOUCLE SUR TOUTES LES FACES D'ENTREE 302! ========================= 303!=============================================================================== 304 305! --- Debit calcule 306 307do izone = 1, nozppm 308 qcalc(izone) = 0.d0 309enddo 310do ifac = 1, nfabor 311 izone = izfppp(ifac) 312 qcalc(izone) = qcalc(izone) - brom(ifac) * & 313 ( rcodcl(ifac,iu,1)*surfbo(1,ifac) + & 314 rcodcl(ifac,iv,1)*surfbo(2,ifac) + & 315 rcodcl(ifac,iw,1)*surfbo(3,ifac) ) 316enddo 317 318if(irangp.ge.0) then 319 call parrsm(nozapm,qcalc) 320endif 321 322do izone = 1, nozapm 323 if (iqimp(izone).eq.0) then 324 qimp(izone) = qcalc(izone) 325 endif 326enddo 327 328! --- Correction des vitesses en norme 329 330do ifac = 1, nfabor 331 izone = izfppp(ifac) 332 if (iqimp(izone).eq.1) then 333 if (abs(qcalc(izone)).gt.epzero) then 334 qisqc = qimp(izone)/qcalc(izone) 335 else 336 qisqc = 0.d0 337 endif 338 rcodcl(ifac,iu,1) = rcodcl(ifac,iu,1)*qisqc 339 rcodcl(ifac,iv,1) = rcodcl(ifac,iv,1)*qisqc 340 rcodcl(ifac,iw,1) = rcodcl(ifac,iw,1)*qisqc 341 endif 342enddo 343 344!=============================================================================== 345! 4. REMPLISSAGE DU TABLEAU DES CONDITIONS LIMITES POUR LA TURBULENCE 346! ON BOUCLE SUR TOUTES LES FACES D'ENTREE 347! ========================= 348! ON DETERMINE LA FAMILLE ET SES PROPRIETES 349! ON IMPOSE LES CONDITIONS AUX LIMITES POUR LA TURBULENCE 350! (pour n'importe quel modele) 351!=============================================================================== 352 353do ifac = 1, nfabor 354 355 izone = izfppp(ifac) 356 357 if (itypfb(ifac).eq.ientre.or.itypfb(ifac).eq.i_convective_inlet) then 358 359 ! Neumann for outflow 360 if (qimp(izone).lt.0.d0) then 361 362 if (itytur.eq.2) then 363 364 icodcl(ifac,ik ) = 3 365 icodcl(ifac,iep) = 3 366 rcodcl(ifac,ik ,3) = 0.d0 367 rcodcl(ifac,iep,3) = 0.d0 368 369 elseif (itytur.eq.3) then 370 371 icodcl(ifac,ir11) = 3 372 icodcl(ifac,ir22) = 3 373 icodcl(ifac,ir33) = 3 374 icodcl(ifac,ir12) = 3 375 icodcl(ifac,ir13) = 3 376 icodcl(ifac,ir23) = 3 377 icodcl(ifac,iep ) = 3 378 rcodcl(ifac,ir11,3) = 0.d0 379 rcodcl(ifac,ir22,3) = 0.d0 380 rcodcl(ifac,ir33,3) = 0.d0 381 rcodcl(ifac,ir12,3) = 0.d0 382 rcodcl(ifac,ir13,3) = 0.d0 383 rcodcl(ifac,ir23,3) = 0.d0 384 rcodcl(ifac,iep, 3) = 0.d0 385 386 elseif (iturb.eq.50) then 387 388 icodcl(ifac,ik ) = 3 389 icodcl(ifac,iep ) = 3 390 icodcl(ifac,iphi) = 3 391 icodcl(ifac,ifb ) = 3 392 rcodcl(ifac,ik, 3) = 0.d0 393 rcodcl(ifac,iep, 3) = 0.d0 394 rcodcl(ifac,iphi,3) = 0.d0 395 rcodcl(ifac,ifb, 3) = 0.d0 396 397 elseif (iturb.eq.60) then 398 399 icodcl(ifac,ik ) = 3 400 icodcl(ifac,iomg) = 3 401 rcodcl(ifac,ik, 3) = 0.d0 402 rcodcl(ifac,iomg,3) = 0.d0 403 404 elseif(iturb.eq.70) then 405 406 icodcl(ifac,inusa) = 3 407 rcodcl(ifac,inusa,3) = 0.d0 408 409 endif 410 411 ! Dirichlet for inflow 412 else 413 414 ! La turbulence est calculee par defaut si ICALKE different de 0 415 ! - soit a partir du diametre hydraulique, d'une vitesse 416 ! de reference adaptes a l'entree courante si ICALKE = 1 417 ! - soit a partir du diametre hydraulique, d'une vitesse 418 ! de reference et de l'intensite turvulente 419 ! adaptes a l'entree courante si ICALKE = 2 420 421 if ( icalke(izone).ne.0 ) then 422 423 uref2 = rcodcl(ifac,iu,1)**2 & 424 + rcodcl(ifac,iv,1)**2 & 425 + rcodcl(ifac,iw,1)**2 426 uref2 = max(uref2,epzero) 427 rhomoy = brom(ifac) 428 iel = ifabor(ifac) 429 viscla = viscl(iel) 430 icke = icalke(izone) 431 dhy = dh(izone) 432 xiturb = xintur(izone) 433 ustar2 = 0.d0 434 xkent = epzero 435 xeent = epzero 436 437 if (icke.eq.1) then 438 ! Calculation of turbulent inlet conditions using 439 ! standard laws for a circular pipe 440 ! (their initialization is not needed here but is good practice). 441 call turbulence_bc_inlet_hyd_diam(ifac, uref2, dhy, rhomoy, viscla, & 442 rcodcl) 443 else if (icke.eq.2) then 444 445 ! Calculation of turbulent inlet conditions using 446 ! the turbulence intensity and standard laws for a circular pipe 447 ! (their initialization is not needed here but is good practice) 448 449 call turbulence_bc_inlet_turb_intensity(ifac, uref2, xiturb, dhy, & 450 rcodcl) 451 452 453 endif 454 455 endif ! icalke 456 457 endif ! qimp 458 459 endif ! itypfb 460 461enddo 462 463!---- 464! FORMATS 465!---- 466 467 468!---- 469! FIN 470!---- 471 472return 473end subroutine 474