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! Example of cs_f_user_boundary_conditions subroutine.f90 for inlet 30! with automatic inlet profile. 31 32! This example assumes the mesh is orthogonal at the inlet. 33! 34!------------------------------------------------------------------------------- 35 36!------------------------------------------------------------------------------- 37! Arguments 38!______________________________________________________________________________. 39! mode name role ! 40!______________________________________________________________________________! 41!> \param[in] nvar total number of variables 42!> \param[in] nscal total number of scalars 43!> \param[out] icodcl boundary condition code: 44!> - 1 Dirichlet 45!> - 2 Radiative outlet 46!> - 3 Neumann 47!> - 4 sliding and 48!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 49!> - 5 smooth wall and 50!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 51!> - 6 rough wall and 52!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 53!> - 9 free inlet/outlet 54!> (input mass flux blocked to 0) 55!> - 13 Dirichlet for the advection operator and 56!> Neumann for the diffusion operator 57!> \param[in] itrifb indirection for boundary faces ordering 58!> \param[in,out] itypfb boundary face types 59!> \param[in,out] izfppp boundary face zone number 60!> \param[in] dt time step (per cell) 61!> \param[in,out] rcodcl boundary condition values: 62!> - rcodcl(1) value of the dirichlet 63!> - rcodcl(2) value of the exterior exchange 64!> coefficient (infinite if no exchange) 65!> - rcodcl(3) value flux density 66!> (negative if gain) in w/m2 67!> -# for the velocity \f$ (\mu+\mu_T) 68!> \gradt \, \vect{u} \cdot \vect{n} \f$ 69!> -# for the pressure \f$ \Delta t 70!> \grad P \cdot \vect{n} \f$ 71!> -# for a scalar \f$ cp \left( K + 72!> \dfrac{K_T}{\sigma_T} \right) 73!> \grad T \cdot \vect{n} \f$ 74!_______________________________________________________________________________ 75 76subroutine cs_f_user_boundary_conditions & 77 ( nvar , nscal , & 78 icodcl , itrifb , itypfb , izfppp , & 79 dt , & 80 rcodcl ) 81 82!=============================================================================== 83 84!=============================================================================== 85! Module files 86!=============================================================================== 87 88use paramx 89use numvar 90use optcal 91use cstphy 92use cstnum 93use entsor 94use parall 95use period 96use ppppar 97use ppthch 98use coincl 99use cpincl 100use ppincl 101use ppcpfu 102use atincl 103use atsoil 104use ctincl 105use cs_fuel_incl 106use mesh 107use field 108use turbomachinery 109use cs_c_bindings 110 111!=============================================================================== 112 113implicit none 114 115! Arguments 116 117integer nvar , nscal 118 119integer icodcl(nfabor,nvar) 120integer itrifb(nfabor), itypfb(nfabor) 121integer izfppp(nfabor) 122 123double precision dt(ncelet) 124double precision rcodcl(nfabor,nvar,3) 125 126! Local variables 127 128!< [loc_var_dec] 129integer ifac, iel, ii, ilelt, nlelt 130 131double precision xustar2, xdh, d2s3, rhomoy 132double precision acc(2), fmprsc, fmul, uref2, vnrm 133 134integer, allocatable, dimension(:) :: lstelt, mrkcel 135 136double precision, dimension(:), pointer :: bfpro_rom 137double precision, dimension(:), pointer :: cvar_k, cvar_ep, cvar_phi 138double precision, dimension(:), pointer :: cvar_omg 139double precision, dimension(:), pointer :: cvar_al, cvar_fb 140double precision, dimension(:), pointer :: cvar_nusa 141double precision, dimension(:), pointer :: cvar_scal 142double precision, dimension(:,:), pointer :: cvar_vel 143double precision, dimension(:,:), pointer :: cvar_rij 144!< [loc_var_dec] 145 146!=============================================================================== 147! Initialization 148!=============================================================================== 149 150!< [init] 151allocate(lstelt(nfabor)) ! temporary array for boundary faces selection 152 153! Map field arrays 154call field_get_val_v(ivarfl(iu), cvar_vel) 155call field_get_val_s(ibrom, bfpro_rom) 156 157if (itytur.eq.2) then 158 159 call field_get_val_s(ivarfl(ik), cvar_k) 160 call field_get_val_s(ivarfl(iep), cvar_ep) 161 162elseif (itytur.eq.3) then 163 call field_get_val_v(ivarfl(irij), cvar_rij) 164 call field_get_val_s(ivarfl(iep), cvar_ep) 165 166 if (iturb.eq.32) then 167 call field_get_val_s(ivarfl(ial), cvar_al) 168 endif 169 170elseif (itytur.eq.5) then 171 172 call field_get_val_s(ivarfl(ik), cvar_k) 173 call field_get_val_s(ivarfl(iep), cvar_ep) 174 call field_get_val_s(ivarfl(iphi), cvar_phi) 175 176 if (iturb.eq.50) then 177 call field_get_val_s(ivarfl(ifb), cvar_fb) 178 elseif (iturb.eq.51) then 179 call field_get_val_s(ivarfl(ial), cvar_al) 180 endif 181 182elseif (iturb.eq.60) then 183 184 call field_get_val_s(ivarfl(ik), cvar_k) 185 call field_get_val_s(ivarfl(iomg), cvar_omg) 186 187elseif (iturb.eq.70) then 188 189 call field_get_val_s(ivarfl(inusa), cvar_nusa) 190 191endif 192 193d2s3 = 2.d0/3.d0 194!< [init] 195 196!=============================================================================== 197! Assign a pseudo-periodic channel type inlet to a set of boundary faces. 198 199! For each subset: 200! - use selection criteria to filter boundary faces of a given subset 201! - loop on faces from a subset 202! - set the boundary condition for each face 203 204! A feedback loop is used so as to progressively reach a state similar 205! to that of a periodic channel at the inlet. 206!=============================================================================== 207 208!< [example_1] 209call getfbr('INLET', nlelt, lstelt) 210!========== 211 212fmprsc = 1.d0 ! mean prescribed velocity 213 214if (ntcabs.eq.1) then 215 216 ! For the Rij-EBRSM model (and possibly V2f), we need a non-flat profile, 217 ! so as to ensure turbulent production, and avoid laminarization; 218 ! here, we simply divide the initial velocity by 10 for inlet 219 ! faces adjacent to the wall. 220 221 ! The loop below assumes wall conditions have been defined first 222 ! (in the GUI, or in this file, before the current test). 223 224 if (iturb.eq.32 .or. itytur.eq.5) then 225 226 allocate(mrkcel(ncelet)) 227 do iel = 1, ncelet 228 mrkcel(iel) = 0 229 enddo 230 231 do ifac = 1, nfabor 232 if (itypfb(ifac) .eq. iparoi) then 233 iel = ifabor(ifac) 234 mrkcel(iel) = 1 235 endif 236 enddo 237 238 endif 239 240 do ilelt = 1, nlelt 241 242 ifac = lstelt(ilelt) 243 iel = ifabor(ifac) 244 245 itypfb(ifac) = ientre 246 247 rcodcl(ifac,iu,1) = - fmprsc * surfbo(1,ifac) / surfbn(ifac) 248 rcodcl(ifac,iv,1) = - fmprsc * surfbo(2,ifac) / surfbn(ifac) 249 rcodcl(ifac,iw,1) = - fmprsc * surfbo(3,ifac) / surfbn(ifac) 250 251 if (iturb.eq.32 .or. itytur.eq.5) then 252 if (mrkcel(iel) .eq. 1) then 253 rcodcl(ifac,iu,1) = fmprsc/10.d0 254 endif 255 endif 256 257 uref2 = rcodcl(ifac,iu,1)**2 & 258 + rcodcl(ifac,iv,1)**2 & 259 + rcodcl(ifac,iw,1)**2 260 uref2 = max(uref2,1.d-12) 261 262 ! Turbulence example computed using equations valid for a pipe. 263 264 ! We will be careful to specify a hydraulic diameter adapted 265 ! to the current inlet. 266 267 ! We will also be careful if necessary to use a more precise 268 ! formula for the dynamic viscosity use in the calculation of 269 ! the Reynolds number (especially if it is variable, it may be 270 ! useful to take the law from 'usphyv'. Here, we use by default 271 ! the 'viscl0" value. 272 ! Regarding the density, we have access to its value at boundary 273 ! faces (romb) so this value is the one used here (specifically, 274 ! it is consistent with the processing in 'usphyv', in case of 275 ! variable density) 276 277 ! Hydraulic diameter 278 xdh = 1.d0 279 280 ! Calculation of turbulent inlet conditions using 281 ! the turbulence intensity and standard laws for a circular pipe 282 ! (their initialization is not needed here but is good practice) 283 284 rhomoy = bfpro_rom(ifac) 285 286 call turbulence_bc_inlet_hyd_diam(ifac, uref2, xdh, rhomoy, viscl0, & 287 rcodcl) 288 289 ! Handle scalars 290 if (nscal.gt.0) then 291 do ii = 1, nscal 292 rcodcl(ifac,isca(ii),1) = 1.d0 293 enddo 294 endif 295 296 enddo 297 298 if (iturb.eq.32 .or. itytur.eq.5) then 299 deallocate(mrkcel) 300 endif 301 302else 303 304! Subsequent time steps 305!---------------------- 306 307 acc(1) = 0.d0 308 acc(2) = 0.d0 309 310 ! Estimate multiplier 311 312 do ilelt = 1, nlelt 313 314 ifac = lstelt(ilelt) 315 iel = ifabor(ifac) 316 317 vnrm = sqrt(cvar_vel(1,iel)**2 + cvar_vel(2,iel)**2 + cvar_vel(3,iel)**2) 318 acc(1) = acc(1) + vnrm*surfbn(ifac) 319 acc(2) = acc(2) + surfbn(ifac) 320 321 enddo 322 323 if (irangp.ge.0) then 324 call parrsm(2, acc) 325 endif 326 327 if (acc(1).le.epzero) then 328 fmul = 0.d0 ! zero velocity in bulk domain 329 else 330 fmul = fmprsc/(acc(1)/acc(2)) ! 1 / estimate flow multiplier 331 endif 332 333 ! Apply BC 334 335 do ilelt = 1, nlelt 336 337 ifac = lstelt(ilelt) 338 iel = ifabor(ifac) 339 340 itypfb(ifac) = ientre 341 342 vnrm = sqrt(cvar_vel(1,iel)**2 + cvar_vel(2,iel)**2 + cvar_vel(3,iel)**2) 343 344 rcodcl(ifac,iu,1) = - fmul * vnrm * surfbo(1,ifac) / surfbn(ifac) 345 rcodcl(ifac,iv,1) = - fmul * vnrm * surfbo(2,ifac) / surfbn(ifac) 346 rcodcl(ifac,iw,1) = - fmul * vnrm * surfbo(3,ifac) / surfbn(ifac) 347 348 if (itytur.eq.2) then 349 350 rcodcl(ifac,ik,1) = cvar_k(iel) 351 rcodcl(ifac,iep,1) = cvar_ep(iel) 352 353 elseif (itytur.eq.3) then 354 355 rcodcl(ifac,ir11,1) = cvar_rij(1,iel) 356 rcodcl(ifac,ir22,1) = cvar_rij(2,iel) 357 rcodcl(ifac,ir33,1) = cvar_rij(3,iel) 358 rcodcl(ifac,ir12,1) = cvar_rij(4,iel) 359 rcodcl(ifac,ir13,1) = cvar_rij(6,iel) 360 rcodcl(ifac,ir23,1) = cvar_rij(5,iel) 361 362 rcodcl(ifac,iep,1) = cvar_ep(iel) 363 364 if (iturb.eq.32) then 365 rcodcl(ifac,ial,1) = cvar_al(iel) 366 endif 367 368 elseif (itytur.eq.5) then 369 370 rcodcl(ifac,ik,1) = cvar_k(iel) 371 rcodcl(ifac,iep,1) = cvar_ep(iel) 372 rcodcl(ifac,iphi,1) = cvar_phi(iel) 373 374 if (iturb.eq.50) then 375 rcodcl(ifac,ifb,1) = cvar_fb(iel) 376 elseif (iturb.eq.51) then 377 rcodcl(ifac,ial,1) = cvar_al(iel) 378 endif 379 380 elseif (iturb.eq.60) then 381 382 rcodcl(ifac,ik,1) = cvar_k(iel) 383 rcodcl(ifac,iomg,1) = cvar_omg(iel) 384 385 elseif (iturb.eq.70) then 386 387 rcodcl(ifac,inusa,1) = cvar_nusa(iel) 388 389 endif 390 391 ! Handle scalars (a correction similar to that of velocity is suggested 392 ! rather than the simpler code below) 393 if (nscal.gt.0) then 394 do ii = 1, nscal 395 call field_get_val_s(ivarfl(isca(ii)), cvar_scal) 396 rcodcl(ifac,isca(ii),1) = cvar_scal(iel) 397 enddo 398 endif 399 400 enddo 401 402endif 403!< [example_1] 404 405!-------- 406! Formats 407!-------- 408 409!---- 410! End 411!---- 412 413deallocate(lstelt) ! temporary array for boundary faces selection 414 415return 416end subroutine cs_f_user_boundary_conditions 417