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 29!> \file cs_user_boundary_conditions-atmospheric.f90 30!> 31!> Atmospheric example of cs_user_boundary_conditions.f90 subroutine 32!> 33!------------------------------------------------------------------------------- 34 35!------------------------------------------------------------------------------- 36! Arguments 37!______________________________________________________________________________. 38! mode name role ! 39!______________________________________________________________________________! 40!> \param[in] nvar total number of variables 41!> \param[in] nscal total number of scalars 42!> \param[out] icodcl boundary condition code: 43!> - 1 Dirichlet 44!> - 2 Radiative outlet 45!> - 3 Neumann 46!> - 4 sliding and 47!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 48!> - 5 smooth wall and 49!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 50!> - 6 rough wall and 51!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 52!> - 9 free inlet/outlet 53!> (input mass flux blocked to 0) 54!> - 13 Dirichlet for the advection operator and 55!> Neumann for the diffusion operator 56!> \param[in] itrifb indirection for boundary faces ordering 57!> \param[in,out] itypfb boundary face types 58!> \param[in,out] izfppp boundary face zone number 59!> \param[in] dt time step (per cell) 60!> \param[in,out] rcodcl boundary condition values: 61!> - rcodcl(1) value of the dirichlet 62!> - rcodcl(2) value of the exterior exchange 63!> coefficient (infinite if no exchange) 64!> - rcodcl(3) value flux density 65!> (negative if gain) in w/m2 66!> -# for the velocity \f$ (\mu+\mu_T) 67!> \gradt \, \vect{u} \cdot \vect{n} \f$ 68!> -# for the pressure \f$ \Delta t 69!> \grad P \cdot \vect{n} \f$ 70!> -# for a scalar \f$ cp \left( K + 71!> \dfrac{K_T}{\sigma_T} \right) 72!> \grad T \cdot \vect{n} \f$ 73!_______________________________________________________________________________ 74 75subroutine cs_f_user_boundary_conditions & 76 ( nvar , nscal , & 77 icodcl , itrifb , itypfb , izfppp , & 78 dt , & 79 rcodcl ) 80 81!=============================================================================== 82 83!=============================================================================== 84! Module files 85!=============================================================================== 86 87use paramx 88use numvar 89use optcal 90use cstphy 91use cstnum 92use entsor 93use parall 94use period 95use ppppar 96use ppthch 97use coincl 98use cpincl 99use ppincl 100use ppcpfu 101use atchem 102use atincl 103use atsoil 104use ctincl 105use cs_fuel_incl 106use mesh 107use field 108use turbomachinery 109use iso_c_binding 110use cs_c_bindings 111 112!=============================================================================== 113 114implicit none 115 116! Arguments 117 118integer nvar , nscal 119 120integer icodcl(nfabor,nvar) 121integer itrifb(nfabor), itypfb(nfabor) 122integer izfppp(nfabor) 123 124double precision dt(ncelet) 125double precision rcodcl(nfabor,nvar,3) 126 127! Local variables 128 129!< [loc_var_dec] 130integer ifac, ii 131integer izone 132integer ilelt, nlelt 133integer f_id_rough, f_id_t_rough 134double precision d2s3 135double precision zref, xuref 136double precision ustar, rugd, rugt 137double precision zent, xuent, xvent 138double precision xkent, xeent 139 140integer, allocatable, dimension(:) :: lstelt 141double precision, dimension(:), pointer :: bpro_roughness 142double precision, dimension(:), pointer :: bpro_roughness_t 143!< [loc_var_dec] 144 145!=============================================================================== 146 147!=============================================================================== 148! Initialization 149!=============================================================================== 150 151allocate(lstelt(nfabor)) ! temporary array for boundary faces selection 152 153d2s3 = 2.d0/3.d0 154 155! Paremeters for the analytical rough wall law (neutral) 156zref = 10.d0 157xuref = 10.d0 158rugd = 0.1d0 159rugt = rugd 160 161!=============================================================================== 162! Assign boundary conditions to boundary faces here 163 164! For each subset: 165! - use selection criteria to filter boundary faces of a given subset 166! - loop on faces from a subset 167! - set the boundary condition for each face 168!=============================================================================== 169 170! --- For boundary faces of color 11, 171! assign an inlet boundary condition for all phases prescribed from the 172! meteo profile with automatic choice between inlet/ outlet according to 173! the meteo profile 174 175!< [example_1] 176call getfbr('11',nlelt,lstelt) 177 178! Get a new zone number (1 <= izone <= nozppm) 179izone = maxval(izfppp) + 1 180 181do ilelt = 1, nlelt 182 183 ifac = lstelt(ilelt) 184 185 ! - Zone to which the face belongs 186 izfppp(ifac) = izone 187 188 ! - Boundary conditions are prescribed from the meteo profile 189 iprofm(izone) = 1 190 191 ! - Chemical boundary conditions are prescribed from the chemistry profile 192 iprofc(izone) = 1 193 194 ! - boundary condition type can be set to ientre or i_convective_inlet 195 196 itypfb(ifac) = ientre 197 198 ! - automatic determination of type (inlet/outlet) according to sign of 199 ! mass flux 200 201 iautom(ifac) = 1 202 203enddo 204!< [example_1] 205 206 207! ---For boundary faces of color 21, 208! assign an inlet boundary condition for all phases prescribed from the 209! meteo profile 210 211!< [example_2] 212call getfbr('21',nlelt,lstelt) 213 214! Get a new zone number (1 <= izone <= nozppm) 215izone = maxval(izfppp) + 1 216 217do ilelt = 1, nlelt 218 219 ifac = lstelt(ilelt) 220 221 ! - Zone to which the face belongs 222 izfppp(ifac) = izone 223 224 ! - Boundary conditions are prescribed from the meteo profile 225 iprofm(izone) = 1 226 227 ! - Chemical boundary conditions are prescribed from the chemistry profile 228 iprofc(izone) = 1 229 230 ! - Assign inlet boundary conditions 231 itypfb(ifac) = ientre 232 233enddo 234!< [example_2] 235 236! --- For boundary faces of color 31, 237! assign an inlet boundary condition for all phases prescribed from the 238! meteo profile except for dynamical variables which are prescribed with 239! a rough log law. 240 241!< [example_3] 242call getfbr('31',nlelt,lstelt) 243 244! Get a new zone number (1 <= izone <= nozppm) 245izone = maxval(izfppp) + 1 246 247do ilelt = 1, nlelt 248 249 ifac = lstelt(ilelt) 250 251 ! - Zone to which the face belongs 252 izfppp(ifac) = izone 253 254 ! - Boundary conditions are prescribed from the meteo profile 255 iprofm(izone) = 1 256 257 ! - Chemical boundary conditions are prescribed from the chemistry profile 258 iprofc(izone) = 1 259 260 ! - Dynamical variables are prescribed with a rough log law 261 zent=cdgfbo(3,ifac) 262 263 ustar = xkappa*xuref/log((zref+rugd)/rugd) 264 xuent = ustar/xkappa*log((zent+rugd)/rugd) 265 xvent = 0.d0 266 xkent = ustar**2/sqrt(cmu) 267 xeent = ustar**3/xkappa/(zent+rugd) 268 269 itypfb(ifac) = ientre 270 271 rcodcl(ifac,iu,1) = xuent 272 rcodcl(ifac,iv,1) = xvent 273 rcodcl(ifac,iw,1) = 0.d0 274 275 ! itytur is a flag equal to iturb/10 276 if (itytur.eq.2) then 277 278 rcodcl(ifac,ik,1) = xkent 279 rcodcl(ifac,iep,1) = xeent 280 281 elseif(itytur.eq.3) then 282 283 rcodcl(ifac,ir11,1) = d2s3*xkent 284 rcodcl(ifac,ir22,1) = d2s3*xkent 285 rcodcl(ifac,ir33,1) = d2s3*xkent 286 rcodcl(ifac,ir12,1) = 0.d0 287 rcodcl(ifac,ir13,1) = 0.d0 288 rcodcl(ifac,ir23,1) = 0.d0 289 rcodcl(ifac,iep,1) = xeent 290 291 elseif(iturb.eq.50) then 292 293 rcodcl(ifac,ik,1) = xkent 294 rcodcl(ifac,iep,1) = xeent 295 rcodcl(ifac,iphi,1) = d2s3 296 rcodcl(ifac,ifb,1) = 0.d0 297 298 elseif(iturb.eq.60) then 299 300 rcodcl(ifac,ik,1) = xkent 301 rcodcl(ifac,iomg,1) = xeent/cmu/xkent 302 303 elseif(iturb.eq.70) then 304 305 rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent 306 307 endif 308 309enddo 310 311!< [example_3] 312! --- Prescribe at boundary faces of color '12' an outlet for all phases 313call getfbr('12', nlelt, lstelt) 314 315! Get a new zone number (1 <= izone <= nozppm) 316izone = maxval(izfppp) + 1 317 318do ilelt = 1, nlelt 319 320 ifac = lstelt(ilelt) 321 322 ! - Zone to which the zone belongs 323 izfppp(ifac) = izone 324 325 ! Outlet: zero flux for velocity and temperature, prescribed pressure 326 ! Note that the pressure will be set to P0 at the first 327 ! free outlet face (isolib) 328 329 itypfb(ifac) = isolib 330 331enddo 332 333! --- Prescribe at boundary faces of color 15 a rough wall 334!< [example_4] 335call field_get_id_try("boundary_roughness", f_id_rough) 336if (f_id_rough.ge.0) call field_get_val_s(f_id_rough, bpro_roughness) 337call field_get_id_try("boundary_thermal_roughness", f_id_t_rough) 338if (f_id_t_rough.ge.0) call field_get_val_s(f_id_t_rough, bpro_roughness_t) 339 340call getfbr('15', nlelt, lstelt) 341 342! Get a new zone number (1 <= izone <= nozppm) 343izone = maxval(izfppp) + 1 344 345do ilelt = 1, nlelt 346 347 ifac = lstelt(ilelt) 348 349 ! Wall: zero flow (zero flux for pressure) 350 ! rough friction for velocities (+ turbulent variables) 351 352 ! - Zone to which the zone belongs 353 izfppp(ifac) = izone 354 355 itypfb(ifac) = iparug 356 357 ! Roughness for velocity: rugd 358 bpro_roughness(ifac) = rugd 359 360 ! Roughness for scalars (if required) 361 if (f_id_rough.ge.0) & 362 bpro_roughness_t(ifac) = 0.01d0 363 364 ! By default zero flux for scalars 365 366enddo 367!< [example_4] 368 369! --- Prescribe at boundary faces of color 4 a symmetry for all phases 370!< [example_5] 371call getfbr('4', nlelt, lstelt) 372 373! Get a new zone number (1 <= izone <= nozppm) 374izone = maxval(izfppp) + 1 375 376do ilelt = 1, nlelt 377 378 ifac = lstelt(ilelt) 379 380 ! - Zone to which the zone belongs 381 izfppp(ifac) = izone 382 383 itypfb(ifac) = isymet 384 385enddo 386!< [example_5] 387 388!-------- 389! Formats 390!-------- 391 392!---- 393! End 394!---- 395 396deallocate(lstelt) ! temporary array for boundary faces selection 397 398return 399end subroutine cs_f_user_boundary_conditions 400