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-fuel.f90
29!> \brief Example of cs_f_user_boundary_conditions subroutine.f90 for fuel
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, ii
122integer          izone
123integer          iclafu
124integer          ilelt, nlelt
125
126double precision uref2, d2s3
127double precision xkent, xeent
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
140
141d2s3 = 2.d0/3.d0
142
143!===============================================================================
144! Assign boundary conditions to boundary faces here
145
146! For each subset:
147! - use selection criteria to filter boundary faces of a given subset
148! - loop on faces from a subset
149!   - set the boundary condition for each face
150!===============================================================================
151
152!  Each kind of condition for extended physic is allocated with a number
153!  IZONE ( 0<IZONE<= NOZPPM ; NOZPPM allocated in ppppar.h)
154
155! ---- The 12 color is a pure air inlet
156
157!< [example_1]
158call getfbr('12', nlelt, lstelt)
159!==========
160
161do ilelt = 1, nlelt
162
163  ifac = lstelt(ilelt)
164
165!   kind of boundary condition for standard variables
166  itypfb(ifac) = ientre
167
168!   Zone number allocation
169  izone = 1
170
171!   Zone number storage
172  izfppp(ifac) = izone
173
174! ------ This inlet have a fixed mass flux
175
176  ientat(izone) = 1
177  iqimp(izone)  = 1
178  inmoxy(izone) = 1
179!     - Air mass flow rate kg/s
180  qimpat(izone) = 1.46d-03
181!     - Air inlet temperature K
182  timpat(izone) = 400.d0 + tkelvi
183!     - Fuel mass flow rate kg/s
184  qimpfl(izone) = 0.d0
185
186! ----- The 12 color is a fixed mass flow rate inlet
187!        user gives only the speed vector direction
188!        (spedd vector norm is irrelevent)
189!
190  rcodcl(ifac,iu,1) = 0.d0
191  rcodcl(ifac,iv,1) = 0.d0
192  rcodcl(ifac,iw,1) = 5.d0
193
194!   Boundary conditions of turbulence
195  icalke(izone) = 1
196!
197!    - If ICALKE = 0 the boundary conditions of turbulence at
198!      the inlet are calculated as follows:
199
200  if (icalke(izone).eq.0) then
201
202    uref2 = rcodcl(ifac,iu,1)**2                           &
203           +rcodcl(ifac,iv,1)**2                           &
204           +rcodcl(ifac,iw,1)**2
205    uref2 = max(uref2,1.d-12)
206    xkent  = epzero
207    xeent  = epzero
208
209    if    (itytur.eq.2) then
210
211      rcodcl(ifac,ik,1)  = xkent
212      rcodcl(ifac,iep,1) = xeent
213
214    elseif(itytur.eq.3) then
215
216      rcodcl(ifac,ir11,1) = d2s3*xkent
217      rcodcl(ifac,ir22,1) = d2s3*xkent
218      rcodcl(ifac,ir33,1) = d2s3*xkent
219      rcodcl(ifac,ir12,1) = 0.d0
220      rcodcl(ifac,ir13,1) = 0.d0
221      rcodcl(ifac,ir23,1) = 0.d0
222      rcodcl(ifac,iep,1)  = xeent
223
224    elseif (iturb.eq.50) then
225
226      rcodcl(ifac,ik,1)   = xkent
227      rcodcl(ifac,iep,1)  = xeent
228      rcodcl(ifac,iphi,1) = d2s3
229      rcodcl(ifac,ifb,1)  = 0.d0
230
231    elseif (iturb.eq.60) then
232
233      rcodcl(ifac,ik,1)   = xkent
234      rcodcl(ifac,iomg,1) = xeent/cmu/xkent
235
236    endif
237
238  endif
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
249! ------ Treatment of user's scalars
250
251  if ( (nscal-nscapp).gt.0 ) then
252    do ii = 1, (nscal-nscapp)
253      rcodcl(ifac,isca(ii),1) = 1.d0
254    enddo
255  endif
256
257enddo
258!< [example_1]
259
260! ---- Inlet of both primary Air and Fuel
261
262!< [example_2]
263call getfbr('11', nlelt, lstelt)
264!==========
265
266do ilelt = 1, nlelt
267
268  ifac = lstelt(ilelt)
269
270!   kind of boundary condition for standard variables
271  itypfb(ifac) = ientre
272
273!   Zone number allocation
274  izone = 2
275
276!   Zone number storage
277  izfppp(ifac) = izone
278
279! ------ This inlet have a fixed mass flux
280
281  ientfl(izone) = 1
282  iqimp(izone)  = 1
283  inmoxy(izone) = 1
284!     - Air mass flow rate in kg/s
285  qimpat(izone) = 1.46d-03
286!     - Air Temperature at inlet in K
287  timpat(izone) = 800.d0  + tkelvi
288
289!     - Fuel mass flow rate in kg/s
290  qimpfl(izone) = 1.46d-04/360.d0
291
292!     - PERCENTAGE mass fraction of each granulometric class
293!       ICLAFU (1 < ICLAFU < NCLAFU )
294  iclafu = 1
295  distfu(izone,iclafu) = 100.d0
296
297!     - Fuel Temperature at inlet in K
298  timpfl(izone) = 100.d0  + tkelvi
299
300! ----- The 11 color is a fixed mass flow rate inlet
301!        user gives only the speed vector direction
302!        (spedd vector norm is irrelevent)
303!
304
305  rcodcl(ifac,iu,1) = 0.d0
306  rcodcl(ifac,iv,1) = 0.d0
307  rcodcl(ifac,iw,1) = 5.d0
308
309!   Boundary conditions of turbulence
310  icalke(izone) = 1
311!
312!    - If ICALKE = 0 the boundary conditions of turbulence at
313!      the inlet are calculated as follows:
314
315  if(icalke(izone).eq.0) then
316
317    uref2 = rcodcl(ifac,iu,1)**2                           &
318           +rcodcl(ifac,iv,1)**2                           &
319           +rcodcl(ifac,iw,1)**2
320    uref2 = max(uref2,1.d-12)
321    xkent  = epzero
322    xeent  = epzero
323
324    if    (itytur.eq.2) then
325
326      rcodcl(ifac,ik,1)  = xkent
327      rcodcl(ifac,iep,1) = xeent
328
329    elseif(itytur.eq.3) then
330
331      rcodcl(ifac,ir11,1) = d2s3*xkent
332      rcodcl(ifac,ir22,1) = d2s3*xkent
333      rcodcl(ifac,ir33,1) = d2s3*xkent
334      rcodcl(ifac,ir12,1) = 0.d0
335      rcodcl(ifac,ir13,1) = 0.d0
336      rcodcl(ifac,ir23,1) = 0.d0
337      rcodcl(ifac,iep,1)  = xeent
338
339    elseif (iturb.eq.50) then
340
341      rcodcl(ifac,ik,1)   = xkent
342      rcodcl(ifac,iep,1)  = xeent
343      rcodcl(ifac,iphi,1) = d2s3
344      rcodcl(ifac,ifb,1)  = 0.d0
345
346    elseif (iturb.eq.60) then
347
348      rcodcl(ifac,ik,1)   = xkent
349      rcodcl(ifac,iomg,1) = xeent/cmu/xkent
350
351    endif
352
353  endif
354!
355!    - If ICALKE = 1 the boundary conditions of turbulence at
356!      the inlet refer to both, a hydraulic diameter and a
357!      reference velocity.
358!
359  dh(izone)     = 0.032d0
360!
361!    - If ICALKE = 2 the boundary conditions of turbulence at
362!      the inlet refer to a turbulence intensity.
363!
364  xintur(izone) = 0.d0
365enddo
366!< [example_2]
367
368! --- Color 15 is a wall
369
370!< [example_3]
371call getfbr('15', nlelt, lstelt)
372!==========
373
374do ilelt = 1, nlelt
375
376  ifac = lstelt(ilelt)
377
378!          WALL  : nul mass flow rate (nul pressure flux)
379!                  rubbing for speed (and turbulence)
380!                  nul scalar fluxes
381
382!   kind of boundary condition for standard variables
383  itypfb(ifac)   = iparoi
384
385
386!   Zone number allocation
387  izone = 3
388
389!   Zone number storage
390  izfppp(ifac) = izone
391
392enddo
393!< [example_3]
394
395! --- Color 19 is an outlet
396
397!< [example_4]
398call getfbr('19', nlelt, lstelt)
399!==========
400
401do ilelt = 1, nlelt
402
403  ifac = lstelt(ilelt)
404
405!          OUTLET : nul fluxes for speed and scalar
406!                   pressure fixed
407
408!   kind of boundary condition for standard variables
409    itypfb(ifac)   = isolib
410
411!   Zone number allocation
412    izone = 4
413
414!   Zone number storage
415    izfppp(ifac) = izone
416
417  enddo
418!< [example_4]
419
420! --- 14 and 4 are symetry
421
422!< [example_5]
423call getfbr('14 or 4', nlelt, lstelt)
424!==========
425
426do ilelt = 1, nlelt
427
428  ifac = lstelt(ilelt)
429
430!          SYMETRY
431
432!   kind of boundary condition for standard variables
433  itypfb(ifac)   = isymet
434
435!   Zone number allocation
436  izone = 5
437
438!   Zone number storage
439  izfppp(ifac) = izone
440
441enddo
442!< [example_5]
443
444!--------
445! Formats
446!--------
447
448!----
449! End
450!----
451
452deallocate(lstelt)  ! temporary array for boundary faces selection
453
454return
455end subroutine cs_f_user_boundary_conditions
456