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_ebu.f90
29!> \brief Example of cs_f_user_boundary_conditions subroutine.f90 for EBU 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          ifac, izone, ii
122integer          ilelt, nlelt
123
124double precision uref2, xkent, xeent, d2s3
125
126integer, allocatable, dimension(:) :: lstelt
127!< [loc_var_dec]
128
129!===============================================================================
130
131!===============================================================================
132! Initialization
133!===============================================================================
134
135allocate(lstelt(nfabor))  ! temporary array for boundary faces selection
136
137d2s3 = 2.d0/3.d0
138
139!===============================================================================
140! Assign boundary conditions to boundary faces here
141
142! For each subset:
143! - use selection criteria to filter boundary faces of a given subset
144! - loop on faces from a subset
145!   - set the boundary condition for each face
146!===============================================================================
147
148! Definition of a burned gas inlet (pilot flame) for each face of color 11
149
150!< [example_1]
151call getfbr('11', nlelt, lstelt)
152!==========
153
154do ilelt = 1, nlelt
155
156  ifac = lstelt(ilelt)
157
158!   Type of pre-defined boundary condidition (see above)
159  itypfb(ifac) = ientre
160
161!   Zone number (arbitrary number between 1 and n)
162  izone = 1
163
164!   Allocation of the actual face to the zone
165  izfppp(ifac) = izone
166
167!   Indicating the inlet as a burned gas inlet
168  ientgb(izone) = 1
169!   The incoming burned gas flow refers to:
170!   a) a massflow rate   -> iqimp()  = 1
171  iqimp(izone) = 0
172  qimp (izone) = zero
173!
174!   b) an inlet velocity -> iqimp()  = 0
175  rcodcl(ifac,iu,1) = 0.d0
176  rcodcl(ifac,iv,1) = 0.d0
177  rcodcl(ifac,iw,1) = 21.47d0
178! ATTENTION: If iqimp()  = 1 the direction vector of the massfow has
179!            to be given here.
180!
181!   Mean Mixture Fraction at Inlet
182  fment(izone) = 1.d0*fs(1)
183!   Inlet Temperature in K
184  tkent(izone) = 2000.d0
185
186!   Boundary Conditions of Turbulence
187  icalke(izone) = 1
188!
189!    - If ICALKE = 0 the boundary conditions of turbulence at
190!      the inlet are calculated as follows:
191
192  if(icalke(izone).eq.0) then
193
194    uref2 = rcodcl(ifac,iu,1)**2                           &
195           +rcodcl(ifac,iv,1)**2                           &
196           +rcodcl(ifac,iw,1)**2
197    uref2 = max(uref2,1.d-12)
198    xkent  = epzero
199    xeent  = epzero
200
201    if    (itytur.eq.2) then
202
203      rcodcl(ifac,ik,1)  = xkent
204      rcodcl(ifac,iep,1) = xeent
205
206    elseif(itytur.eq.3) then
207
208      rcodcl(ifac,ir11,1) = d2s3*xkent
209      rcodcl(ifac,ir22,1) = d2s3*xkent
210      rcodcl(ifac,ir33,1) = d2s3*xkent
211      rcodcl(ifac,ir12,1) = 0.d0
212      rcodcl(ifac,ir13,1) = 0.d0
213      rcodcl(ifac,ir23,1) = 0.d0
214      rcodcl(ifac,iep,1)  = xeent
215
216    elseif (iturb.eq.50) then
217
218      rcodcl(ifac,ik,1)   = xkent
219      rcodcl(ifac,iep,1)  = xeent
220      rcodcl(ifac,iphi,1) = d2s3
221      rcodcl(ifac,ifb,1)  = 0.d0
222
223    elseif (iturb.eq.60) then
224
225      rcodcl(ifac,ik,1)   = xkent
226      rcodcl(ifac,iomg,1) = xeent/cmu/xkent
227
228    elseif (iturb.eq.70) then
229
230      rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent
231
232    endif
233
234  endif
235!
236!    - If ICALKE = 1 the boundary conditions of turbulence at
237!      the inlet refer to both, a hydraulic diameter and a
238!      reference velocity.
239!
240  dh(izone)     = 0.032d0
241!
242!    - If ICALKE = 2 the boundary conditions of turbulence at
243!      the inlet refer to a turbulence intensity.
244!
245  xintur(izone) = 0.d0
246
247enddo
248!< [example_1]
249
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