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