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