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-pulverized_coal.f90
29!> \brief Example of cs_f_user_boundary_conditions subroutine.f90
30!>        for pulverized coal
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, ii
123integer          izone
124integer          icha, iclapc
125integer          ilelt, nlelt
126
127double precision uref2, d2s3
128double precision xkent, xeent
129
130integer, allocatable, dimension(:) :: lstelt
131!< [loc_var_dec]
132
133!===============================================================================
134
135!===============================================================================
136! Initialization
137!===============================================================================
138
139allocate(lstelt(nfabor))  ! temporary array for boundary faces selection
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! ---- BOUNDARY FACE corresponding to AIR INLET
153!      e.g. : secondary or tertiary air
154
155!< [example_1]
156call getfbr('12', nlelt, lstelt)
157!==========
158
159do ilelt = 1, nlelt
160
161  ifac = lstelt(ilelt)
162
163!   Kind of boundary conditions for standard variables
164  itypfb(ifac) = ientre
165
166!   zone's number (from 1 to n)
167  izone = 1
168
169!      - Allocation of the zone's number to the face
170  izfppp(ifac) = izone
171
172!      - For theses inlet faces, mass flux is fixed
173
174  ientat(izone) = 1
175  iqimp(izone)  = 1
176!      - Oxidizer's number (1 to 3)
177  inmoxy(izone) = 1
178!      - Oxidizer mass flow rate  in kg/s
179  qimpat(izone) = 1.46d-03
180!      - Oxidizer's Temperature in K
181  timpat(izone) = 400.d0 + tkelvi
182
183!      - The color 12 becomes an fixed flow rate inlet
184!        The user gives speed vector direction
185!        (speed vector norm is irrelevent)
186
187  rcodcl(ifac,iu,1) = 0.d0
188  rcodcl(ifac,iv,1) = 0.d0
189  rcodcl(ifac,iw,1) = 5.d0
190
191! ------ Turbulence treatment
192!   Boundary conditions of turbulence
193  icalke(izone) = 1
194!
195!    - If ICALKE = 0 the boundary conditions of turbulence at
196!      the inlet are calculated as follows:
197
198  if(icalke(izone).eq.0) then
199
200    uref2 = rcodcl(ifac,iu,1)**2                           &
201           +rcodcl(ifac,iv,1)**2                           &
202           +rcodcl(ifac,iw,1)**2
203    uref2 = max(uref2,1.d-12)
204    xkent  = epzero
205    xeent  = epzero
206
207    if    (itytur.eq.2) then
208
209      rcodcl(ifac,ik,1)  = xkent
210      rcodcl(ifac,iep,1) = xeent
211
212    elseif(itytur.eq.3) then
213
214      rcodcl(ifac,ir11,1) = d2s3*xkent
215      rcodcl(ifac,ir22,1) = d2s3*xkent
216      rcodcl(ifac,ir33,1) = d2s3*xkent
217      rcodcl(ifac,ir12,1) = 0.d0
218      rcodcl(ifac,ir13,1) = 0.d0
219      rcodcl(ifac,ir23,1) = 0.d0
220      rcodcl(ifac,iep,1)  = xeent
221
222    elseif (iturb.eq.50) then
223
224      rcodcl(ifac,ik,1)   = xkent
225      rcodcl(ifac,iep,1)  = xeent
226      rcodcl(ifac,iphi,1) = d2s3
227      rcodcl(ifac,ifb,1)  = 0.d0
228
229    elseif (iturb.eq.60) then
230
231      rcodcl(ifac,ik,1)   = xkent
232      rcodcl(ifac,iomg,1) = xeent/cmu/xkent
233
234    elseif (iturb.eq.70) then
235
236      rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent
237
238    endif
239
240  endif
241!
242!    - If ICALKE = 1 the boundary conditions of turbulence at
243!      the inlet refer to both, a hydraulic diameter and a
244!      reference velocity.
245!
246  dh(izone)     = 0.032d0
247!
248!    - If ICALKE = 2 the boundary conditions of turbulence at
249!      the inlet refer to a turbulence intensity.
250!
251  xintur(izone) = 0.d0
252
253! ------ Automatic treatment of scalars for extended physic
254
255
256! ------ treatment of user's scalars
257
258  if ( (nscal-nscapp).gt.0 ) then
259    do ii = 1, (nscal-nscapp)
260      rcodcl(ifac,isca(ii),1) = 1.d0
261    enddo
262  endif
263
264enddo
265!< [example_1]
266
267! ---- BOUNDARY FACE for pulverised COAL & primary air INLET
268
269!< [example_2]
270call getfbr('11', nlelt, lstelt)
271!==========
272
273do ilelt = 1, nlelt
274
275  ifac = lstelt(ilelt)
276
277!   Kind of boundary conditions for standard variables
278  itypfb(ifac) = ientre
279
280!   zone's number (from 1 to n)
281  izone = 2
282
283!      - Allocation of the zone's number to the face
284  izfppp(ifac) = izone
285
286!      - For theses inlet faces, mass flux is fixed
287
288  ientcp(izone) = 1
289  iqimp(izone)  = 1
290!      - Oxidizer's number (1 to 3)
291  inmoxy(izone) = 1
292!      - Oxidizer's mass flow rate in kg/s
293  qimpat(izone) = 1.46d-03
294!      - Oxidizer's Temperature in K
295  timpat(izone) = 800.d0  + tkelvi
296
297!        Coal inlet, initialization
298  do icha = 1, ncharm
299    qimpcp(izone,icha) = zero
300    timpcp(izone,icha) = zero
301    do iclapc = 1, ncpcmx
302      distch(izone,icha,iclapc) = zero
303    enddo
304  enddo
305
306! Code_Saturne deals with NCHA different coals (component of blend)
307!       every coal is described by NCLPCH(icha) class of particles
308!       (each of them described by an inlet diameter)
309!
310!      - Treatment for the first coal
311  icha = 1
312!      - Coal mass flow rate in kg/s
313   qimpcp(izone,icha) = 1.46d-4
314!      - PERCENTAGE mass fraction of each granulometric class
315  do iclapc = 1, nclpch(icha)
316    distch(izone,icha,iclapc) = 100.d0/dble(nclpch(icha))
317  enddo
318!      - Inlet temperature for coal & primary air
319  timpcp(izone,icha) = 800.d0 + tkelvi
320
321!      - The color 11 becomes an fixed flow rate inlet
322!        The user gives speed vector direction
323!        (speed vector norm is irrelevent)
324
325  rcodcl(ifac,iu,1) = 0.d0
326  rcodcl(ifac,iv,1) = 0.d0
327  rcodcl(ifac,iw,1) = 5.d0
328
329! PPl
330! ------ Traitement de la turbulence
331
332!        La turbulence est calculee par defaut si ICALKE different de 0
333!          - soit a partir du diametre hydraulique, d'une vitesse
334!            de reference adaptes a l'entree courante si ICALKE = 1
335!          - soit a partir du diametre hydraulique, d'une vitesse
336!            de reference et de l'intensite turvulente
337!            adaptes a l'entree courante si ICALKE = 2
338
339!      Choix pour le calcul automatique ICALKE = 1 ou 2
340  icalke(izone) = 1
341!      Saisie des donnees
342  dh(izone)     = 0.1d0
343  xintur(izone) = 0.1d0
344
345! PPl
346!
347enddo
348!< [example_2]
349
350!     The color 15 become a WALL
351
352!< [example_3]
353call getfbr('15', nlelt, lstelt)
354!==========
355
356do ilelt = 1, nlelt
357
358  ifac = lstelt(ilelt)
359
360!          WALL : NUL MASS FLUX (PRESSURE FLUX is zero valued)
361!                 FRICTION FOR SPEED (& TURBULENCE)
362!                 NUL SCALAR FLUX
363
364!   Kind of boundary conditions for standard variables
365  itypfb(ifac)   = iparoi
366
367
368!   zone's number (from 1 to n)
369  izone = 3
370
371!      - Allocation of the zone's number to the face
372  izfppp(ifac) = izone
373
374enddo
375!< [example_3]
376
377
378!     The color 19 becomes an OUTLET
379
380!< [example_4]
381call getfbr('19', nlelt, lstelt)
382!==========
383
384do ilelt = 1, nlelt
385
386  ifac = lstelt(ilelt)
387
388!          OUTLET : NUL FLUX for SPEED & SCALARS, FIXED PRESSURE
389
390!   Kind of boundary conditions for standard variables
391  itypfb(ifac)   = isolib
392
393!   zone's number (from 1 to n)
394  izone = 4
395
396!      - Allocation of the zone's number to the face
397  izfppp(ifac) = izone
398
399enddo
400!< [example_4]
401
402!     The color 14 becomes a symetry plane
403
404!< [example_5]
405call getfbr('14 or 4', nlelt, lstelt)
406!==========
407
408do ilelt = 1, nlelt
409
410  ifac = lstelt(ilelt)
411
412!          SYMETRIES
413
414!   Kind of boundary conditions for standard variables
415  itypfb(ifac)   = isymet
416
417!   zone's number (from 1 to n)
418  izone = 5
419
420!      - Allocation of the zone's number to the face
421  izfppp(ifac) = izone
422
423enddo
424!< [example_5]
425
426
427!--------
428! Formats
429!--------
430
431!----
432! End
433!----
434
435deallocate(lstelt)  ! temporary array for boundary faces selection
436
437return
438end subroutine cs_f_user_boundary_conditions
439