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_lagrangian.f90
29!> \brief Example of cs_f_user_boundary_conditions subroutine.f90 for
30!>        lagrangian 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
124
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! ---- Face de type entree correspondant a une entree d'air
153!        Par exemple : Air primaire , secondaire ou Air tertiaire
154
155!< [example_1]
156call getfbr('12', nlelt, lstelt)
157!==========
158
159do ilelt = 1, nlelt
160
161  ifac = lstelt(ilelt)
162
163!   Type de condition aux limites pour les variables standard
164  itypfb(ifac) = ientre
165
166!   Numero de zone (on les numerote de 1 a n)
167  izone = 1
168
169!      - Reperage de la zone a laquelle appartient la face
170  izfppp(ifac) = izone
171
172! ------ Pour ces faces d'entree , on est a debit impose
173
174  ientat(izone) = 1
175  iqimp(izone)  = 1
176!      - Debit en kg/s pour l'air
177  qimpat(izone) = 1.46d-03
178!      - Temperature en K pour l'air
179  timpat(izone) = 400.d0 + tkelvi
180
181
182! ------ On impose en couleur 12 une entree a debit impose
183!        L'utilisateur donne donc ici uniquement
184!          la direction du vecteur vitesse
185
186  rcodcl(ifac,iu,1) = 0.d0
187  rcodcl(ifac,iv,1) = 0.d0
188  rcodcl(ifac,iw,1) = 5.d0
189
190! ------ Traitement de la turbulence
191
192!        La turbulence est calculee par defaut si ICALKE different de 0
193!          - soit a partir du diametre hydraulique, d'une vitesse
194!            de reference adaptes a l'entree courante si ICALKE = 1
195!          - soit a partir du diametre hydraulique, d'une vitesse
196!            de reference et de l'intensite turvulente
197!            adaptes a l'entree courante si ICALKE = 2
198
199!      Choix pour le calcul automatique ICALKE = 1 ou 2
200  icalke(izone) = 1
201!      Saisie des donnees
202  dh(izone)     = 0.1d0
203  xintur(izone) = 0.1d0
204
205
206
207! Exemple de cas ou ICALKE(IZONE) = 0 : DEBUT
208!    Eliminer ces lignes pour la clarte si on a fait le choix ICALKE(IZONE) = 1
209
210  if(icalke(izone).eq.0) then
211
212!         Calcul de k et epsilon en entree (XKENT et XEENT) a partir
213!           l'intensite turbulente et de lois standards en conduite
214!           circulaire (leur initialisation est inutile mais plus
215!           propre)
216    uref2 = rcodcl(ifac,iu,1)**2                           &
217           +rcodcl(ifac,iv,1)**2                           &
218           +rcodcl(ifac,iw,1)**2
219    uref2 = max(uref2,1.d-12)
220    xkent  = epzero
221    xeent  = epzero
222
223!     (ITYTUR est un indicateur qui vaut ITURB/10)
224    if    (itytur.eq.2) then
225
226      rcodcl(ifac,ik,1)  = xkent
227      rcodcl(ifac,iep,1) = xeent
228
229    elseif(itytur.eq.3) then
230
231      rcodcl(ifac,ir11,1) = d2s3*xkent
232      rcodcl(ifac,ir22,1) = d2s3*xkent
233      rcodcl(ifac,ir33,1) = d2s3*xkent
234      rcodcl(ifac,ir12,1) = 0.d0
235      rcodcl(ifac,ir13,1) = 0.d0
236      rcodcl(ifac,ir23,1) = 0.d0
237      rcodcl(ifac,iep,1)  = xeent
238
239    elseif (iturb.eq.50) then
240
241      rcodcl(ifac,ik,1)   = xkent
242      rcodcl(ifac,iep,1)  = xeent
243      rcodcl(ifac,iphi,1) = d2s3
244      rcodcl(ifac,ifb,1)  = 0.d0
245
246    elseif (iturb.eq.60) then
247
248      rcodcl(ifac,ik,1)   = xkent
249      rcodcl(ifac,iomg,1) = xeent/cmu/xkent
250
251    elseif (iturb.eq.70) then
252
253      rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent
254
255    endif
256
257  endif
258! Exemple de cas ou ICALKE(IZONE) = 0 : FIN
259
260! ------ Traitement des scalaires physiques particulieres
261!        Ils sont traites automatiquement
262
263
264! ------ Traitement des scalaires utilisateurs
265
266  if ( (nscal-nscapp).gt.0 ) then
267    do ii = 1, (nscal-nscapp)
268      rcodcl(ifac,isca(ii),1) = 1.d0
269    enddo
270  endif
271
272enddo
273!< [example_1]
274
275! --- On impose en couleur 15 une paroi laterale
276
277!< [example_2]
278call getfbr('15', nlelt, lstelt)
279!==========
280
281do ilelt = 1, nlelt
282
283  ifac = lstelt(ilelt)
284
285!          PAROI : DEBIT NUL (FLUX NUL POUR LA PRESSION)
286!                  FROTTEMENT POUR LES VITESSES (+GRANDEURS TURB)
287!                  FLUX NUL SUR LES SCALAIRES
288
289!   Type de condition aux limites pour les variables standard
290  itypfb(ifac)   = iparoi
291
292
293!   Numero de zone (on les numerote de 1 a n)
294  izone = 2
295
296!      - Reperage de la zone a laquelle appartient la face
297  izfppp(ifac) = izone
298
299enddo
300!< [example_2]
301
302! --- On impose en couleur 19 une sortie
303
304!< [example_3]
305call getfbr('19', nlelt, lstelt)
306!==========
307
308do ilelt = 1, nlelt
309
310  ifac = lstelt(ilelt)
311
312!          SORTIE : FLUX NUL VITESSE ET TEMPERATURE, PRESSION IMPOSEE
313!            Noter que la pression sera recalee a P0
314!                sur la premiere face de sortie libre (ISOLIB)
315
316!   Type de condition aux limites pour les variables standard
317  itypfb(ifac)   = isolib
318
319!   Numero de zone (on les numerote de 1 a n)
320  izone = 3
321
322!      - Reperage de la zone a laquelle appartient la face
323  izfppp(ifac) = izone
324
325enddo
326!< [example_3]
327
328
329! --- On impose en couleur 14 et 4 une symetrie
330
331!< [example_4]
332call getfbr('14 or 4', nlelt, lstelt)
333!==========
334
335do ilelt = 1, nlelt
336
337  ifac = lstelt(ilelt)
338
339!          SYMETRIES
340
341!   Type de condition aux limites pour les variables standard
342  itypfb(ifac)   = isymet
343
344!   Numero de zone (on les numerote de 1 a n)
345  izone = 4
346
347!      - Reperage de la zone a laquelle appartient la face
348  izfppp(ifac) = izone
349
350enddo
351!< [example_4]
352
353!--------
354! Formats
355!--------
356
357!----
358! End
359!----
360
361deallocate(lstelt)  ! temporary array for boundary faces selection
362
363return
364end subroutine cs_f_user_boundary_conditions
365