1!------------------------------------------------------------------------------- 2 3!VERS 4 5! This file is part of Code_Saturne, a general-purpose CFD tool. 6! 7! Copyright (C) 1998-2021 EDF S.A. 8! 9! This program is free software; you can redistribute it and/or modify it under 10! the terms of the GNU General Public License as published by the Free Software 11! Foundation; either version 2 of the License, or (at your option) any later 12! version. 13! 14! This program is distributed in the hope that it will be useful, but WITHOUT 15! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 16! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 17! details. 18! 19! You should have received a copy of the GNU General Public License along with 20! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 21! Street, Fifth Floor, Boston, MA 02110-1301, USA. 22 23!------------------------------------------------------------------------------- 24 25!=============================================================================== 26! Function: 27! --------- 28!> \file cs_user_boundary_conditions-gas_3pthem.f90 29!> \brief Example of cs_user_boundary_conditions subroutine.f90 for 3 PTHEM gas 30! 31!------------------------------------------------------------------------------- 32 33!------------------------------------------------------------------------------- 34! Arguments 35!______________________________________________________________________________. 36! mode name role ! 37!______________________________________________________________________________! 38!> \param[in] nvar total number of variables 39!> \param[in] nscal total number of scalars 40!> \param[out] icodcl boundary condition code: 41!> - 1 Dirichlet 42!> - 2 Radiative outlet 43!> - 3 Neumann 44!> - 4 sliding and 45!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 46!> - 5 smooth wall and 47!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 48!> - 6 rough wall and 49!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 50!> - 9 free inlet/outlet 51!> (input mass flux blocked to 0) 52!> \param[in] itrifb indirection for boundary faces ordering 53!> \param[in,out] itypfb boundary face types 54!> \param[in,out] izfppp boundary face zone number 55!> \param[in] dt time step (per cell) 56!> \param[in,out] rcodcl boundary condition values: 57!> - rcodcl(1) value of the dirichlet 58!> - rcodcl(2) value of the exterior exchange 59!> coefficient (infinite if no exchange) 60!> - rcodcl(3) value flux density 61!> (negative if gain) in w/m2 62!> -# for the velocity \f$ (\mu+\mu_T) 63!> \gradt \, \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!_______________________________________________________________________________ 70 71subroutine cs_f_user_boundary_conditions & 72 ( nvar , nscal , & 73 icodcl , itrifb , itypfb , izfppp , & 74 dt , & 75 rcodcl ) 76 77!=============================================================================== 78 79!=============================================================================== 80! Module files 81!=============================================================================== 82 83use paramx 84use numvar 85use optcal 86use cstphy 87use cstnum 88use entsor 89use parall 90use period 91use ppppar 92use ppthch 93use coincl 94use cpincl 95use ppincl 96use ppcpfu 97use atincl 98use ctincl 99use cs_fuel_incl 100use mesh 101use field 102 103!=============================================================================== 104 105implicit none 106 107! Arguments 108 109integer nvar , nscal 110 111integer icodcl(nfabor,nvar) 112integer itrifb(nfabor), itypfb(nfabor) 113integer izfppp(nfabor) 114 115double precision dt(ncelet) 116double precision rcodcl(nfabor,nvar,3) 117 118! Local variables 119 120!< [loc_var_dec] 121integer iel, ifac, izone, ii 122integer ilelt, nlelt 123 124double precision uref2, d2s3 125double precision xkent, xeent 126double precision dp, rho, S0, Kadm, radm, madm, Qadm, Padm, Ploc 127double precision, dimension(:), pointer :: crom 128 129integer, allocatable, dimension(:) :: lstelt 130!< [loc_var_dec] 131 132!=============================================================================== 133 134!=============================================================================== 135! Initialization 136!=============================================================================== 137 138allocate(lstelt(nfabor)) ! temporary array for boundary faces selection 139 140d2s3 = 2.d0/3.d0 141 142!=============================================================================== 143! Assign boundary conditions to boundary faces here 144 145! For each subset: 146! - use selection criteria to filter boundary faces of a given subset 147! - loop on faces from a subset 148! - set the boundary condition for each face 149!=============================================================================== 150 151! Definition of a fuel flow inlet for each face of color 11 152 153!< [example_1] 154call getfbr('11', nlelt, lstelt) 155!========== 156 157do ilelt = 1, nlelt 158 159 ifac = lstelt(ilelt) 160 161! Type of pre-defined boundary condidition (see above) 162 itypfb(ifac) = ientre 163 164! Zone number (arbitrary number between 1 and n) 165 izone = 1 166 167! Allocation of the actual face to the zone 168 izfppp(ifac) = izone 169 170! Indicating the inlet as a fuel flow inlet 171 ientfu(izone) = 1 172 173! Inlet Temperature in K 174 tinfue = 436.d0 175 176! The incoming fuel flow refers to: 177! a) a massflow rate -> iqimp() = 1 178 iqimp(izone) = 1 179 qimp(izone) = 2.62609d-4 / 72.d0 180! 181! b) an inlet velocity -> iqimp() = 0 182 rcodcl(ifac,iu,1) = 0.d0 183 rcodcl(ifac,iv,1) = 0.d0 184 rcodcl(ifac,iw,1) = 21.47d0 185! ATTENTION: If iqimp() = 1 the direction vector of the massfow has 186! to begiven here. 187! 188! Boundary conditions of turbulence 189 icalke(izone) = 1 190! 191! - If ICALKE = 0 the boundary conditions of turbulence at 192! the inlet are calculated as follows: 193 194 if(icalke(izone).eq.0) then 195 196 uref2 = rcodcl(ifac,iu,1)**2 & 197 +rcodcl(ifac,iv,1)**2 & 198 +rcodcl(ifac,iw,1)**2 199 uref2 = max(uref2,1.d-12) 200 xkent = epzero 201 xeent = epzero 202 203 if (itytur.eq.2) then 204 205 rcodcl(ifac,ik,1) = xkent 206 rcodcl(ifac,iep,1) = xeent 207 208 elseif(itytur.eq.3) then 209 210 rcodcl(ifac,ir11,1) = d2s3*xkent 211 rcodcl(ifac,ir22,1) = d2s3*xkent 212 rcodcl(ifac,ir33,1) = d2s3*xkent 213 rcodcl(ifac,ir12,1) = 0.d0 214 rcodcl(ifac,ir13,1) = 0.d0 215 rcodcl(ifac,ir23,1) = 0.d0 216 rcodcl(ifac,iep,1) = xeent 217 218 elseif (iturb.eq.50) then 219 220 rcodcl(ifac,ik,1) = xkent 221 rcodcl(ifac,iep,1) = xeent 222 rcodcl(ifac,iphi,1) = d2s3 223 rcodcl(ifac,ifb,1) = 0.d0 224 225 elseif (iturb.eq.60) then 226 227 rcodcl(ifac,ik,1) = xkent 228 rcodcl(ifac,iomg,1) = xeent/cmu/xkent 229 230 elseif (iturb.eq.70) then 231 232 rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent 233 234 endif 235 236 endif 237! 238! - If ICALKE = 1 the boundary conditions of turbulence at 239! the inlet refer to both, a hydraulic diameter and a 240! reference velocity. 241! 242 dh(izone) = 0.032d0 243! 244! - If ICALKE = 2 the boundary conditions of turbulence at 245! the inlet refer to a turbulence intensity. 246! 247 xintur(izone) = 0.d0 248 249enddo 250!< [example_1] 251 252! Definition of an air flow inlet for each face of color 21 253 254!< [example_2] 255call getfbr('21', nlelt, lstelt) 256!========== 257 258do ilelt = 1, nlelt 259 260 ifac = lstelt(ilelt) 261 262! Type of pre-defined boundary condidition (see above) 263 itypfb(ifac) = ientre 264 265! Zone number (arbitrary number between 1 and n) 266 izone = 2 267 268! Allocation of the actual face to the zone 269 izfppp(ifac) = izone 270 271! Indicating the inlet as a air flow inlet 272 ientox(izone) = 1 273 274! Inlet Temperature in K 275 tinoxy = 353.d0 276 277! The inflowing fuel flow refers to: 278! a) a massflow rate -> iqimp() = 1 279 iqimp(izone) = 1 280 qimp(izone) = 4.282d-3 / 72.d0 281! 282! b) an inlet velocity -> iqimp() = 0 283 rcodcl(ifac,iu,1) = 0.d0 284 rcodcl(ifac,iv,1) = 0.d0 285 rcodcl(ifac,iw,1) = 0.097d0 286! ATTENTION: If iqimp() = 1 the direction vector of the massfow has 287! to be given here. 288 289! Boundary conditions of turbulence 290 icalke(izone) = 1 291! 292! - If ICALKE = 0 the boundary conditions of turbulence at 293! the inlet are calculated as follows: 294 295 if(icalke(izone).eq.0) then 296 297 uref2 = rcodcl(ifac,iu,1)**2 & 298 +rcodcl(ifac,iv,1)**2 & 299 +rcodcl(ifac,iw,1)**2 300 uref2 = max(uref2,1.d-12) 301 xkent = epzero 302 xeent = epzero 303 304 if (itytur.eq.2) then 305 306 rcodcl(ifac,ik,1) = xkent 307 rcodcl(ifac,iep,1) = xeent 308 309 elseif(itytur.eq.3) then 310 311 rcodcl(ifac,ir11,1) = d2s3*xkent 312 rcodcl(ifac,ir22,1) = d2s3*xkent 313 rcodcl(ifac,ir33,1) = d2s3*xkent 314 rcodcl(ifac,ir12,1) = 0.d0 315 rcodcl(ifac,ir13,1) = 0.d0 316 rcodcl(ifac,ir23,1) = 0.d0 317 rcodcl(ifac,iep,1) = xeent 318 319 elseif (iturb.eq.50) then 320 321 rcodcl(ifac,ik,1) = xkent 322 rcodcl(ifac,iep,1) = xeent 323 rcodcl(ifac,iphi,1) = d2s3 324 rcodcl(ifac,ifb,1) = 0.d0 325 326 elseif (iturb.eq.60) then 327 328 rcodcl(ifac,ik,1) = xkent 329 rcodcl(ifac,iomg,1) = xeent/cmu/xkent 330 331 elseif (iturb.eq.70) then 332 333 rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent 334 335 endif 336 337 endif 338! 339! - If ICALKE = 1 the boundary conditions of turbulence at 340! the inlet refer to both, a hydraulic diameter and a 341! reference velocity. 342! 343 dh(izone) = 0.218d0 344! 345! - If ICALKE = 2 the boundary conditions of turbulence at 346! the inlet refer to a turbulence intensity. 347! 348 xintur(izone) = 0.d0 349! 350enddo 351!< [example_2] 352 353! Definition of a wall for each face of color 51 up to 59 354 355!< [example_3] 356call getfbr('51 to 59', nlelt, lstelt) 357!========== 358 359do ilelt = 1, nlelt 360 361 ifac = lstelt(ilelt) 362 363! Type of pre-defined boundary condidition (see above) 364 itypfb(ifac) = iparoi 365 366! Zone number (arbitrary number between 1 and n) 367 izone = 3 368 369! Allocation of the actual face to the zone 370 izfppp(ifac) = izone 371 372enddo 373!< [example_3] 374 375 376! Definition of an exit for each face of color 91 377 378!< [example_4] 379call getfbr('91', nlelt, lstelt) 380!========== 381 382do ilelt = 1, nlelt 383 384 ifac = lstelt(ilelt) 385 386! Type of pre-defined boundary condidition (see above) 387 itypfb(ifac) = isolib 388 389! Zone number (arbitrary number between 1 and n) 390 izone = 4 391 392! Allocation of the actual face to the zone 393 izfppp(ifac) = izone 394 395enddo 396!< [example_4] 397 398! Definition of symmetric boundary conditions for each 399! face of color 41 and 4. 400 401!< [example_5] 402call getfbr('41 or 4', nlelt, lstelt) 403!========== 404 405do ilelt = 1, nlelt 406 407 ifac = lstelt(ilelt) 408 409! Type of pre-defined boundary condidition (see above) 410 itypfb(ifac) = isymet 411 412! Zone number (arbitrary number between 1 and n) 413 izone = 5 414 415! Allocation of the actual face to the zone 416 izfppp(ifac) = izone 417 418enddo 419!< [example_5] 420 421! Definition of an air supply for each face of color 61 422 423!< [example_6] 424call getfbr('61', nlelt, lstelt) 425!========== 426 427call field_get_val_s(icrom, crom) 428 429! Head losses in the ventilation pipe (often from exp data) 430 431! nominal room pressure 432Ploc= 51.d0 ! Pa 433! vent surface 434S0 = acos(-1.d0) * (75.d-3)**2 ! m2 435! density 436radm = (P0+ploc)/(cs_physical_constants_r*tinoxy/wmolg(2)) 437! nominal vent pressure 438Padm = 119.d0 ! Pa 439! nominal flow rate 440Qadm = 504.d0 ! m3/h 441madm = Qadm*radm/3600.d0 442! head loss 443Kadm = 2.d0*radm*(S0/madm)**2 * abs(Ploc-Padm) 444 445! pressure gradient (Pther computed if ipthrm = 1) 446dp = (Pther - P0) - Padm 447 448do ilelt = 1, nlelt 449 450 ifac = lstelt(ilelt) 451 iel = ifabor(ifac) 452 453 ! Density 454 if (dp.gt.0.d0) then 455 rho = crom(iel) 456 else 457 rho = radm 458 endif 459 460 ! Mass flow rate (opposed to the boundary face normal vector) 461 madm = - sign(1.d0,dp) * s0 * sqrt(2.d0*rho/kadm*abs(dp)) 462 463 ! Type of pre-defined boundary condidition (see above) 464 itypfb(ifac) = i_convective_inlet 465 466 ! Zone number (arbitrary number between 1 and n) 467 izone = 6 468 469 ! Allocation of the actual face to the zone 470 izfppp(ifac) = izone 471 472 ! Indicating the inlet as a air flow inlet 473 ientox(izone) = 1 474 475 ! Inlet Temperature in K 476 tinoxy = t0 477 478 ! The inflowing fuel flow refers to: 479 ! a) a massflow rate -> iqimp() = 1 480 iqimp(izone) = 1 481 qimp(izone) = madm 482 483 ! b) an inlet velocity -> iqimp() = 0 484 rcodcl(ifac,iu,1) = 1.d0 485 rcodcl(ifac,iv,1) = 0.d0 486 rcodcl(ifac,iw,1) = 0.d0 487 488 ! Boundary conditions of turbulence 489 icalke(izone) = 1 490 491 ! - If ICALKE = 0 the boundary conditions of turbulence at 492 ! the inlet are calculated as follows: 493 if(icalke(izone).eq.0) then 494 495 uref2 = rcodcl(ifac,iu,1)**2 & 496 +rcodcl(ifac,iv,1)**2 & 497 +rcodcl(ifac,iw,1)**2 498 uref2 = max(uref2,1.d-12) 499 xkent = epzero 500 xeent = epzero 501 502 if (itytur.eq.2) then 503 504 rcodcl(ifac,ik,1) = xkent 505 rcodcl(ifac,iep,1) = xeent 506 507 elseif(itytur.eq.3) then 508 509 rcodcl(ifac,ir11,1) = d2s3*xkent 510 rcodcl(ifac,ir22,1) = d2s3*xkent 511 rcodcl(ifac,ir33,1) = d2s3*xkent 512 rcodcl(ifac,ir12,1) = 0.d0 513 rcodcl(ifac,ir13,1) = 0.d0 514 rcodcl(ifac,ir23,1) = 0.d0 515 rcodcl(ifac,iep,1) = xeent 516 517 elseif (iturb.eq.50) then 518 519 rcodcl(ifac,ik,1) = xkent 520 rcodcl(ifac,iep,1) = xeent 521 rcodcl(ifac,iphi,1) = d2s3 522 rcodcl(ifac,ifb,1) = 0.d0 523 524 elseif (iturb.eq.60) then 525 526 rcodcl(ifac,ik,1) = xkent 527 rcodcl(ifac,iomg,1) = xeent/cmu/xkent 528 529 elseif (iturb.eq.70) then 530 531 rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent 532 533 endif 534 535 endif 536 537 ! - If ICALKE = 1 the boundary conditions of turbulence at 538 ! the inlet refer to both, a hydraulic diameter and a 539 ! reference velocity. 540 dh(izone) = 0.218d0 541 542 ! - If ICALKE = 2 the boundary conditions of turbulence at 543 ! the inlet refer to a turbulence intensity. 544 xintur(izone) = 0.d0 545 546enddo 547!< [example_6] 548 549!-------- 550! Formats 551!-------- 552 553!---- 554! End 555!---- 556 557deallocate(lstelt) ! temporary array for boundary faces selection 558 559return 560end subroutine cs_f_user_boundary_conditions 561