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_libby_williams.f90 29!> \brief Example of cs_f_user_boundary_conditions subroutine.f90 for 30!> Libby-Williams gas 31! 32!------------------------------------------------------------------------------- 33 34!------------------------------------------------------------------------------- 35! Arguments 36!______________________________________________________________________________. 37! mode name role ! 38!______________________________________________________________________________! 39!> \param[in] nvar total number of variables 40!> \param[in] nscal total number of scalars 41!> \param[out] icodcl boundary condition code: 42!> - 1 Dirichlet 43!> - 2 Radiative outlet 44!> - 3 Neumann 45!> - 4 sliding and 46!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 47!> - 5 smooth wall and 48!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 49!> - 6 rough wall and 50!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 51!> - 9 free inlet/outlet 52!> (input mass flux blocked to 0) 53!> \param[in] itrifb indirection for boundary faces ordering 54!> \param[in,out] itypfb boundary face types 55!> \param[in,out] izfppp boundary face zone number 56!> \param[in] dt time step (per cell) 57!> \param[in,out] rcodcl boundary condition values: 58!> - rcodcl(1) value of the dirichlet 59!> - rcodcl(2) value of the exterior exchange 60!> coefficient (infinite if no exchange) 61!> - rcodcl(3) value flux density 62!> (negative if gain) in w/m2 63!> -# for the velocity \f$ (\mu+\mu_T) 64!> \gradt \, \vect{u} \cdot \vect{n} \f$ 65!> -# for the pressure \f$ \Delta t 66!> \grad P \cdot \vect{n} \f$ 67!> -# for a scalar \f$ cp \left( K + 68!> \dfrac{K_T}{\sigma_T} \right) 69!> \grad T \cdot \vect{n} \f$ 70!_______________________________________________________________________________ 71 72subroutine cs_f_user_boundary_conditions & 73 ( nvar , nscal , & 74 icodcl , itrifb , itypfb , izfppp , & 75 dt , & 76 rcodcl ) 77 78!=============================================================================== 79 80!=============================================================================== 81! Module files 82!=============================================================================== 83 84use paramx 85use numvar 86use optcal 87use cstphy 88use cstnum 89use entsor 90use parall 91use period 92use ppppar 93use ppthch 94use coincl 95use cpincl 96use ppincl 97use ppcpfu 98use atincl 99use ctincl 100use cs_fuel_incl 101use mesh 102use field 103 104!=============================================================================== 105 106implicit none 107 108! Arguments 109 110integer nvar , nscal 111 112integer icodcl(nfabor,nvar) 113integer itrifb(nfabor), itypfb(nfabor) 114integer izfppp(nfabor) 115 116double precision dt(ncelet) 117double precision rcodcl(nfabor,nvar,3) 118 119! Local variables 120 121!< [loc_var_dec] 122integer ifac, izone, ii 123integer ilelt, nlelt 124 125double precision uref2, xkent, xeent, d2s3 126 127integer, allocatable, dimension(:) :: lstelt 128!< [loc_var_dec] 129 130!=============================================================================== 131 132!=============================================================================== 133! Initialization 134!=============================================================================== 135 136allocate(lstelt(nfabor)) ! temporary array for boundary faces selection 137 138d2s3 = 2.d0/3.d0 139 140!=============================================================================== 141! Assign boundary conditions to boundary faces here 142 143! For each subset: 144! - use selection criteria to filter boundary faces of a given subset 145! - loop on faces from a subset 146! - set the boundary condition for each face 147!=============================================================================== 148 149! Definition of a burned gas inlet (pilot flame) for each face of color 11 150 151!< [example_1] 152call getfbr('11', nlelt, lstelt) 153!========== 154 155do ilelt = 1, nlelt 156 157 ifac = lstelt(ilelt) 158 159 ! Type of pre-defined boundary condidition (see above) 160 itypfb(ifac) = ientre 161 162 ! Zone number (arbitrary number between 1 and n) 163 izone = 1 164 165 ! Allocation of the actual face to the zone 166 izfppp(ifac) = izone 167 168 ! Indicating the inlet as a burned gas inlet 169 ientgb(izone) = 1 170 ! The incoming burned gas flow refers to: 171 ! a) a massflow rate -> iqimp() = 1 172 iqimp(izone) = 0 173 qimp (izone) = zero 174 175 ! b) an inlet velocity -> iqimp() = 0 176 rcodcl(ifac,iu,1) = 0.d0 177 rcodcl(ifac,iv,1) = 0.d0 178 rcodcl(ifac,iw,1) = 21.47d0 179 ! ATTENTION: If iqimp() = 1 the direction vector of the massfow has 180 ! to be given here. 181 182 ! Mean Mixture Fraction at Inlet 183 fment(izone) = 1.d0*fs(1) 184 ! Inlet Temperature in K 185 tkent(izone) = 2000.d0 186 187 ! Boundary Conditions of Turbulence 188 icalke(izone) = 1 189 190 ! - If ICALKE = 0 the boundary conditions of turbulence at 191 ! the inlet are calculated as follows: 192 193 if(icalke(izone).eq.0) then 194 195 uref2 = rcodcl(ifac,iu,1)**2 & 196 +rcodcl(ifac,iv,1)**2 & 197 +rcodcl(ifac,iw,1)**2 198 uref2 = max(uref2,1.d-12) 199 xkent = epzero 200 xeent = epzero 201 202 if (itytur.eq.2) then 203 204 rcodcl(ifac,ik,1) = xkent 205 rcodcl(ifac,iep,1) = xeent 206 207 elseif(itytur.eq.3) then 208 209 rcodcl(ifac,ir11,1) = d2s3*xkent 210 rcodcl(ifac,ir22,1) = d2s3*xkent 211 rcodcl(ifac,ir33,1) = d2s3*xkent 212 rcodcl(ifac,ir12,1) = 0.d0 213 rcodcl(ifac,ir13,1) = 0.d0 214 rcodcl(ifac,ir23,1) = 0.d0 215 rcodcl(ifac,iep,1) = xeent 216 217 elseif (iturb.eq.50) then 218 219 rcodcl(ifac,ik,1) = xkent 220 rcodcl(ifac,iep,1) = xeent 221 rcodcl(ifac,iphi,1) = d2s3 222 rcodcl(ifac,ifb,1) = 0.d0 223 224 elseif (iturb.eq.60) then 225 226 rcodcl(ifac,ik,1) = xkent 227 rcodcl(ifac,iomg,1) = xeent/cmu/xkent 228 229 elseif (iturb.eq.70) then 230 231 rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent 232 233 endif 234 235 endif 236 237 ! - If ICALKE = 1 the boundary conditions of turbulence at 238 ! the inlet refer to both, a hydraulic diameter and a 239 ! reference velocity. 240 ! 241 dh(izone) = 0.032d0 242 243 ! - If ICALKE = 2 the boundary conditions of turbulence at 244 ! the inlet refer to a turbulence intensity. 245 246 xintur(izone) = 0.d0 247 248enddo 249!< [example_1] 250 251! Definition of an unburned gas inlet for each face of color 12 252 253!< [example_2] 254call getfbr('12', nlelt, lstelt) 255!========== 256 257do ilelt = 1, nlelt 258 259 ifac = lstelt(ilelt) 260 261 ! Type of pre-defined boundary condidition (see above) 262 itypfb(ifac) = ientre 263 264 ! Zone number (arbitrary number between 1 and n) 265 izone = 2 266 267 ! Allocation of the actual face to the zone 268 izfppp(ifac) = izone 269 270 ! Indicating the inlet as an unburned gas inlet 271 ientgf(izone) = 1 272 ! The incoming unburned gas flow refers to: 273 ! a) a massflow rate -> iqimp() = 1 274 iqimp(izone) = 0 275 qimp(izone) = zero 276 277 ! b) an inlet velocity -> iqimp() = 0 278 rcodcl(ifac,iu,1) = 60.d0 279 rcodcl(ifac,iv,1) = 0.d0 280 rcodcl(ifac,iw,1) = 0.d0 281 ! ATTENTION: If iqimp() = 1 the direction vector of the massfow has 282 ! to be given here. 283 284 ! Mean Mixture Fraction at Inlet 285 fment(izone) = 8.d-1*fs(1) 286 287 ! Inlet Temperature in K 288 tkent(izone) = 600.d0 289 290 ! Boundary Conditions of Turbulence 291 icalke(izone) = 1 292 293 ! - If ICALKE = 1 the boundary conditions of turbulence at 294 ! the inlet refer to both, a hydraulic diameter and a 295 ! reference velocity. 296 297 dh(izone) = 0.08d0 298 299 ! - If ICALKE = 2 the boundary conditions of turbulence at 300 ! the inlet refer to a turbulence intensity. 301 302 xintur(izone) = 0.1d0 303 304enddo 305!< [example_2] 306 307! Definition of a wall for each face of color 51 and 5 308 309!< [example_3] 310call getfbr('51 or 5', nlelt, lstelt) 311!========== 312 313do ilelt = 1, nlelt 314 315 ifac = lstelt(ilelt) 316 317 ! Type de condition aux limites pour les variables standard 318 itypfb(ifac) = iparoi 319 320 ! Zone number (arbitrary number between 1 and n) 321 izone = 4 322 323 ! Allocation of the actual face to the zone 324 izfppp(ifac) = izone 325 326enddo 327!< [example_3] 328 329! Definition of an exit for each face of color 91 and 9 330 331!< [example_4] 332call getfbr('91 or 9', nlelt, lstelt) 333!========== 334 335do ilelt = 1, nlelt 336 337 ifac = lstelt(ilelt) 338 339 ! Type de condition aux limites pour les variables standard 340 itypfb(ifac) = isolib 341 342 ! Zone number (arbitrary number between 1 and n) 343 izone = 5 344 345 ! Allocation of the actual face to the zone 9 346 izfppp(ifac) = izone 347 348enddo 349!< [example_4] 350 351! Definition of symmetric boundary conditions for each 352! face of color 41 and 4. 353 354!< [example_5] 355call getfbr('41 or 4', nlelt, lstelt) 356!========== 357 358do ilelt = 1, nlelt 359 360 ifac = lstelt(ilelt) 361 362 ! Type de condition aux limites pour les variables standard 363 itypfb(ifac) = isymet 364 365 ! Zone number (arbitrary number between 1 and n) 366 izone = 6 367 368 ! Allocation of the actual face to the zonec 369 izfppp(ifac) = izone 370 371enddo 372!< [example_5] 373 374!-------- 375! Formats 376!-------- 377 378!---- 379! End 380!---- 381 382deallocate(lstelt) ! temporary array for boundary faces selection 383 384return 385end subroutine cs_f_user_boundary_conditions 386