← Back to branch summary

~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to users/fuel/usfucl.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-01 17:43:32 UTC
  • mto: (6.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 11.
  • Revision ID: package-import@ubuntu.com-20111101174332-tl4vk45no0x3emc3
Tags: upstream-2.1.0
Import upstream version 2.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/fuel/usfucl.f90
1
 
!-------------------------------------------------------------------------------
2
 
 
3
 
!VERS
4
 
 
5
 
 
6
 
!     This file is part of the Code_Saturne Kernel, element of the
7
 
!     Code_Saturne CFD tool.
8
 
 
9
 
!     Copyright (C) 1998-2009 EDF S.A., France
10
 
 
11
 
!     contact: saturne-support@edf.fr
12
 
 
13
 
!     The Code_Saturne Kernel is free software; you can redistribute it
14
 
!     and/or modify it under the terms of the GNU General Public License
15
 
!     as published by the Free Software Foundation; either version 2 of
16
 
!     the License, or (at your option) any later version.
17
 
 
18
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
19
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
20
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21
 
!     GNU General Public License for more details.
22
 
 
23
 
!     You should have received a copy of the GNU General Public License
24
 
!     along with the Code_Saturne Kernel; if not, write to the
25
 
!     Free Software Foundation, Inc.,
26
 
!     51 Franklin St, Fifth Floor,
27
 
!     Boston, MA  02110-1301  USA
28
 
 
29
 
!-------------------------------------------------------------------------------
30
 
 
31
 
subroutine usfucl &
32
 
!================
33
 
 
34
 
 ( idbia0 , idbra0 ,                                              &
35
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
36
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
37
 
   nvar   , nscal  , nphas  ,                                     &
38
 
   nideve , nrdeve , nituse , nrtuse ,                            &
39
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , &
40
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
41
 
   icodcl , itrifb , itypfb , izfppp ,                            &
42
 
   idevel , ituser , ia     ,                                     &
43
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
44
 
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
45
 
   coefa  , coefb  , rcodcl ,                                     &
46
 
   w1     , w2     , w3     , w4     , w5     , w6     , coefu  , &
47
 
   rdevel , rtuser , ra     )
48
 
 
49
 
!===============================================================================
50
 
! PURPOSE  :
51
 
! --------
52
 
 
53
 
!    USER SUBROUTINE for extended physic
54
 
!                Combsution of heavy Fuel oil
55
 
!    Allocation of boundary conditions (ICODCL, RCODCL)
56
 
!    for variables unknowns during USCLIM
57
 
 
58
 
 
59
 
! Introduction
60
 
! ============
61
 
 
62
 
! Here we define boundary conditions on a per-face basis.
63
 
 
64
 
! Boundary faces may be identified using the 'getfbr' subroutine.
65
 
 
66
 
!  getfbr(string, nelts, eltlst) :
67
 
!  - string is a user-supplied character string containing
68
 
!    selection criteria;
69
 
!  - nelts is set by the subroutine. It is an integer value
70
 
!    corresponding to the number of boundary faces verifying the
71
 
!    selection criteria;
72
 
!  - lstelt is set by the subroutine. It is an integer array of
73
 
!    size nelts containing the list of boundary faces verifying
74
 
!    the selection criteria.
75
 
 
76
 
!  string may contain:
77
 
!  - references to colors (ex.: 1, 8, 26, ...
78
 
!  - references to groups (ex.: inlet, group1, ...)
79
 
!  - geometric criteria (ex. x < 0.1, y >= 0.25, ...)
80
 
!  These criteria may be combined using logical operators
81
 
!  ('and', 'or') and parentheses.
82
 
!  Example: '1 and (group2 or group3) and y < 1' will select boundary
83
 
!  faces of color 1, belonging to groups 'group2' or 'group3' and
84
 
!  with face center coordinate y less than 1.
85
 
 
86
 
 
87
 
 
88
 
! Boundary condition types
89
 
! ========================
90
 
 
91
 
! Boundary conditions may be assigned in two ways.
92
 
 
93
 
 
94
 
!    For "standard" boundary conditions:
95
 
!    -----------------------------------
96
 
 
97
 
!     (inlet, free outlet, wall, symmetry), we define a code
98
 
!     in the 'itypfb' array (of dimensions number of boundary faces,
99
 
!     number of phases). This code will then be used by a non-user
100
 
!     subroutine to assign the following conditions (scalars in
101
 
!     particular will receive the conditions of the phase to which
102
 
!     they are assigned). Thus:
103
 
 
104
 
!     Code      |  Boundary type
105
 
!     --------------------------
106
 
!      ientre   |   Inlet
107
 
!      isolib   |   Free outlet
108
 
!      isymet   |   Symmetry
109
 
!      iparoi   |   Wall (smooth)
110
 
!      iparug   |   Rough wall
111
 
 
112
 
!     Integers ientre, isolib, isymet, iparoi, iparug
113
 
!     are defined elsewhere (param.h). Their value is greater than
114
 
!     or equal to 1 and less than or equal to ntypmx
115
 
!     (value fixed in paramx.h)
116
 
 
117
 
 
118
 
!     In addition, some values must be defined:
119
 
 
120
 
 
121
 
!     - Inlet (more precisely, inlet/outlet with prescribed flow, as
122
 
!              the flow may be prescribed as an outflow):
123
 
 
124
 
!       -> Dirichlet conditions on variables
125
 
!         other than pressure are mandatory if the flow is incoming,
126
 
!         optional if the flow is outgoing (the code assigns 0 flux
127
 
!         if no Dirichlet is specified); thus,
128
 
!         at face 'ifac', for the variable 'ivar': rcodcl(ifac, ivar, 1)
129
 
 
130
 
 
131
 
!     - Smooth wall: (= impermeable solid, with smooth friction)
132
 
 
133
 
!       -> Velocity value for sliding wall if applicable
134
 
!         at face ifac, rcodcl(ifac, iu, 1)
135
 
!                       rcodcl(ifac, iv, 1)
136
 
!                       rcodcl(ifac, iw, 1)
137
 
!       -> Specific code and prescribed temperature value
138
 
!         at wall, if applicable:
139
 
!         at face ifac, icodcl(ifac, ivar)    = 5
140
 
!                       rcodcl(ifac, ivar, 1) = prescribed temperature
141
 
!       -> Specific code and prescribed flux value
142
 
!         at wall, if applicable:
143
 
!         at face ifac, icodcl(ifac, ivar)    = 3
144
 
!                       rcodcl(ifac, ivar, 3) = prescribed flux
145
 
!                                        =
146
 
!        Note that the default condition for scalars
147
 
!         (other than k and epsilon) is homogeneous Neumann.
148
 
 
149
 
 
150
 
!     - Rough wall: (= impermeable solid, with rough friction)
151
 
 
152
 
!       -> Velocity value for sliding wall if applicable
153
 
!         at face ifac, rcodcl(ifac, iu, 1)
154
 
!                       rcodcl(ifac, iv, 1)
155
 
!                       rcodcl(ifac, iw, 1)
156
 
!       -> Value of the dynamic roughness height to specify in
157
 
!                       rcodcl(ifac, iu, 3) (value for iv et iw not used)
158
 
!       -> Specific code and prescribed temperature value
159
 
!         at rough wall, if applicable:
160
 
!         at face ifac, icodcl(ifac, ivar)    = 6
161
 
!                       rcodcl(ifac, ivar, 1) = prescribed temperature
162
 
!                       rcodcl(ifac, ivar, 3) = dynamic roughness height
163
 
!       -> Specific code and prescribed flux value
164
 
!         at rough wall, if applicable:
165
 
!         at face ifac, icodcl(ifac, ivar)    = 3
166
 
!                       rcodcl(ifac, ivar, 3) = prescribed flux
167
 
!                                        =
168
 
!        Note that the default condition for scalars
169
 
!         (other than k and epsilon) is homogeneous Neumann.
170
 
 
171
 
!     - Symmetry (= impermeable frictionless wall):
172
 
 
173
 
!       -> Nothing to specify
174
 
 
175
 
 
176
 
!     - Free outlet (more precisely free inlet/outlet with prescribed pressure)
177
 
 
178
 
!       -> Nothing to prescribe for pressure and velocity
179
 
!          For scalars and turbulent values, a Dirichlet value may optionally
180
 
!            be specified. The behavior is as follows:
181
 
!              * pressure is always handled as a Dirichlet condition
182
 
!              * if the mass flow is inflowing:
183
 
!                  we retain the velocity at infinity
184
 
!                  Dirichlet condition for scalars and turbulent values
185
 
!                    (or zero flux if the user has not specified a
186
 
!                    Dirichlet value)
187
 
!                if the mass flow is outflowing:
188
 
!                  we prescribe zero flux on the velocity, the scalars,
189
 
!                  and turbulent values
190
 
 
191
 
!       Note that the pressure will be reset to P0
192
 
!           on the first free outlet face found
193
 
 
194
 
 
195
 
!    For "non-standard" conditions:
196
 
!    ------------------------------
197
 
 
198
 
!     Other than (inlet, free outlet, wall, symmetry), we define
199
 
!      - on one hand, for each face:
200
 
!        -> an admissible 'itypfb' value
201
 
!           (i.e. greater than or equal to 1 and less than or equal to
202
 
!            ntypmx; see its value in paramx.h).
203
 
!           The values predefined in paramx.h:
204
 
!           'ientre', 'isolib', 'isymet', 'iparoi', 'iparug' are in
205
 
!           this range, and it is preferable not to assign one of these
206
 
!           integers to 'itypfb' randomly or in an inconsiderate manner.
207
 
!           To avoid this, we may use 'iindef' if we wish to avoid
208
 
!           checking values in paramx.h. 'iindef' is an admissible
209
 
!           value to which no predefined boundary condition is attached.
210
 
!           Note that the 'itypfb' array is reinitialized at each time
211
 
!           step to the non-admissible value of 0. If we forget to
212
 
!           modify 'typfb' for a given face, the code will stop.
213
 
 
214
 
!      - and on the other hand, for each face and each variable:
215
 
!        -> a code             icodcl(ifac, ivar)
216
 
!        -> three real values  rcodcl(ifac, ivar, 1)
217
 
!                              rcodcl(ifac, ivar, 2)
218
 
!                              rcodcl(ifac, ivar, 3)
219
 
!     The value of 'icodcl' is taken from the following:
220
 
!       1: Dirichlet      (usable for any variable)
221
 
!       3: Neumann        (usable for any variable)
222
 
!       4: Symmetry       (usable only for the velocity and
223
 
!                          components of the Rij tensor)
224
 
!       5: Smooth wall    (usable for any variable except for pressure)
225
 
!       6: Rough wall     (usable for any variable except for pressure)
226
 
!       9: Free outlet    (usable only for velocity)
227
 
!     The values of the 3 'rcodcl' components are
228
 
!      rcodcl(ifac, ivar, 1):
229
 
!         Dirichlet for the variable          if icodcl(ifac, ivar) =  1
230
 
!         wall value (sliding velocity, temp) if icodcl(ifac, ivar) =  5
231
 
!         The dimension of rcodcl(ifac, ivar, 1) is that of the
232
 
!           resolved variable: ex U (velocity in m/s),
233
 
!                                 T (temperature in degrees)
234
 
!                                 H (enthalpy in J/kg)
235
 
!                                 F (passive scalar in -)
236
 
!      rcodcl(ifac, ivar, 2):
237
 
!         "exterior" exchange coefficient (between the prescribed value
238
 
!                          and the value at the domain boundary)
239
 
!                          rinfin = infinite by default
240
 
!         For velocities U,                in kg/(m2 s):
241
 
!           rcodcl(ifac, ivar, 2) =          (viscl+visct) / d
242
 
!         For the pressure P,              in  s/m:
243
 
!           rcodcl(ifac, ivar, 2) =                     dt / d
244
 
!         For temperatures T,              in Watt/(m2 degres):
245
 
!           rcodcl(ifac, ivar, 2) = Cp*(viscls+visct/sigmas) / d
246
 
!         For enthalpies H,                in kg /(m2 s):
247
 
!           rcodcl(ifac, ivar, 2) =    (viscls+visct/sigmas) / d
248
 
!         For other scalars F              in:
249
 
!           rcodcl(ifac, ivar, 2) =    (viscls+visct/sigmas) / d
250
 
!              (d has the dimension of a distance in m)
251
 
!
252
 
!      rcodcl(ifac, ivar, 3) if icodcl(ifac, ivar) <> 6:
253
 
!        Flux density (< 0 if gain, n outwards-facing normal)
254
 
!                         if icodcl(ifac, ivar)= 3
255
 
!         For velocities U,                in kg/(m s2) = J:
256
 
!           rcodcl(ifac, ivar, 3) =         -(viscl+visct) * (grad U).n
257
 
!         For pressure P,                  en kg/(m2 s):
258
 
!           rcodcl(ifac, ivar, 3) =                    -dt * (grad P).n
259
 
!         For temperatures T,              in Watt/m2:
260
 
!           rcodcl(ifac, ivar, 3) = -Cp*(viscls+visct/sigmas) * (grad T).n
261
 
!         For enthalpies H,                in Watt/m2:
262
 
!           rcodcl(ifac, ivar, 3) = -(viscls+visct/sigmas) * (grad H).n
263
 
!         For other scalars F in :
264
 
!           rcodcl(ifac, ivar, 3) = -(viscls+visct/sigmas) * (grad F).n
265
 
 
266
 
!      rcodcl(ifac, ivar, 3) if icodcl(ifac, ivar) = 6:
267
 
!        Roughness for the rough wall law
268
 
!         For velocities U, dynamic roughness
269
 
!           rcodcl(ifac, ivar, 3) = rugd
270
 
!         For other scalars, thermal roughness
271
 
!           rcodcl(ifac, ivar, 3) = rugt
272
 
 
273
 
 
274
 
!      Note that if the user assigns a value to itypfb equal to
275
 
!       ientre, isolib, isymet, iparoi, or iparug
276
 
!       and does not modify icodcl (zero value by default),
277
 
!       itypfb will define the boundary condition type.
278
 
 
279
 
!      To the contrary, if the user prescribes
280
 
!        icodcl(ifac, ivar) (nonzero),
281
 
!        the values assigned to rcodcl will be used for the considered
282
 
!        face and variable (if rcodcl values are not set, the default
283
 
!        values will be used for the face and variable, so:
284
 
!                                 rcodcl(ifac, ivar, 1) = 0.d0
285
 
!                                 rcodcl(ifac, ivar, 2) = rinfin
286
 
!                                 rcodcl(ifac, ivar, 3) = 0.d0)
287
 
!        Especially, we may have for example:
288
 
!        -> set itypfb(ifac, iphas) = iparoi
289
 
!        which prescribes default wall conditions for all variables at
290
 
!        face ifac,
291
 
!        -> and define IN ADDITION for variable ivar on this face
292
 
!        specific conditions by specifying
293
 
!        icodcl(ifac, ivar) and the 3 rcodcl values.
294
 
 
295
 
 
296
 
!      The user may also assign to itypfb a value not equal to
297
 
!       ientre, isolib, isymet, iparoi, iparug, iindef
298
 
!       but greater than or equal to 1 and less than or equal to
299
 
!       ntypmx (see values in param.h) to distinguish
300
 
!       groups or colors in other subroutines which are specific
301
 
!       to the case and in which itypfb is accessible.
302
 
!       In this case though it will be necessary to
303
 
!       prescribe boundary conditions by assigning values to
304
 
!       icodcl and to the 3 rcodcl fields (as the value of itypfb
305
 
!       will not be predefined in the code).
306
 
 
307
 
 
308
 
! Consistency rules
309
 
! =================
310
 
 
311
 
!       A few consistency rules between 'icodcl' codes for
312
 
!         variables with non-standard boundary conditions:
313
 
 
314
 
!           Codes for velocity components must be identical
315
 
!           Codes for Rij components must be identical
316
 
!           If code (velocity or Rij) = 4
317
 
!             we must have code (velocity and Rij) = 4
318
 
!           If code (velocity or turbulence) = 5
319
 
!             we must have code (velocity and turbulence) = 5
320
 
!           If code (velocity or turbulence) = 6
321
 
!             we must have code (velocity and turbulence) = 6
322
 
!           If scalar code (except pressure or fluctuations) = 5
323
 
!             we must have velocity code = 5
324
 
!           If scalar code (except pressure or fluctuations) = 6
325
 
!             we must have velocity code = 6
326
 
 
327
 
 
328
 
! Remarks
329
 
! =======
330
 
 
331
 
!       Caution: to prescribe a flux (nonzero) to Rij,
332
 
!                the viscosity to take into account is viscl
333
 
!                even if visct exists (visct=rho cmu k2/epsilon)
334
 
 
335
 
!       We have the ordering array for boundary faces from the
336
 
!           previous time step (except for the fist time step,
337
 
!           where 'itrifb' has not been set yet).
338
 
!       The array of boundary face types 'itypfb' has been
339
 
!           reset before entering the subroutine.
340
 
 
341
 
 
342
 
!       Note how to access some variables:
343
 
 
344
 
! Cell values
345
 
!               Let         iel = ifabor(ifac)
346
 
 
347
 
! * Density                         phase iphas, cell iel:
348
 
!                  propce(iel, ipproc(irom(iphas)))
349
 
! * Dynamic molecular viscosity     phase iphas, cell iel:
350
 
!                  propce(iel, ipproc(iviscl(iphas)))
351
 
! * Turbulent viscosity   dynamique phase iphas, cell iel:
352
 
!                  propce(iel, ipproc(ivisct(iphas)))
353
 
! * Specific heat                   phase iphas, cell iel:
354
 
!                  propce(iel, ipproc(icp(iphasl))
355
 
! * Diffusivity: lambda          scalaire iscal, cell iel:
356
 
!                  propce(iel, ipproc(ivisls(iscal)))
357
 
 
358
 
! Boundary face values
359
 
 
360
 
! * Density                        phase iphas, boundary face ifac :
361
 
!                  propfb(ifac, ipprob(irom(iphas)))
362
 
! * Mass flow relative to variable ivar, boundary face ifac:
363
 
!      (i.e. the mass flow used for convecting ivar)
364
 
!                  propfb(ifac, pprob(ifluma(ivar )))
365
 
! * For other values                  at boundary face ifac:
366
 
!      take as an approximation the value in the adjacent cell iel
367
 
!      i.e. as above with iel = ifabor(ifac).
368
 
 
369
 
!-------------------------------------------------------------------------------
370
 
!   All cells can be identified by using the subroutine 'getcel'.
371
 
!    Syntax of getcel:
372
 
!     getcel(string, nelts, eltlst) :
373
 
!     - string is a user-supplied character string containing
374
 
!       selection criteria;
375
 
!     - nelts is set by the subroutine. It is an integer value
376
 
!       corresponding to the number of boundary faces verifying the
377
 
!       selection criteria;
378
 
!     - lstelt is set by the subroutine. It is an integer array of
379
 
!       size nelts containing the list of boundary faces verifying
380
 
!       the selection criteria.
381
 
 
382
 
!       string may contain:
383
 
!       - references to colors (ex.: 1, 8, 26, ...
384
 
!       - references to groups (ex.: inlet, group1, ...)
385
 
!       - geometric criteria (ex. x < 0.1, y >= 0.25, ...)
386
 
!
387
 
!       These criteria may be combined using logical operators
388
 
!       ('and', 'or') and parentheses.
389
 
!       Example: '1 and (group2 or group3) and y < 1' will select boundary
390
 
!       faces of color 1, belonging to groups 'group2' or 'group3' and
391
 
!       with face center coordinate y less than 1.
392
 
!
393
 
!   All boundary faces may be identified using the 'getfbr' subroutine.
394
 
!    Syntax of getfbr:
395
 
!     getfbr(string, nelts, eltlst) :
396
 
!     - string is a user-supplied character string containing
397
 
!       selection criteria;
398
 
!     - nelts is set by the subroutine. It is an integer value
399
 
!       corresponding to the number of boundary faces verifying the
400
 
!       selection criteria;
401
 
!     - lstelt is set by the subroutine. It is an integer array of
402
 
!       size nelts containing the list of boundary faces verifying
403
 
!       the selection criteria.
404
 
!
405
 
!     string may contain:
406
 
!     - references to colors (ex.: 1, 8, 26, ...
407
 
!     - references to groups (ex.: inlet, group1, ...)
408
 
!     - geometric criteria (ex. x < 0.1, y >= 0.25, ...)
409
 
!
410
 
!     These criteria may be combined using logical operators
411
 
!     ('and', 'or') and parentheses.
412
 
!
413
 
!   All internam faces may be identified using the 'getfac' subroutine.
414
 
!    Syntax of getfac:
415
 
!     getfac(string, nelts, eltlst) :
416
 
!     - string is a user-supplied character string containing
417
 
!       selection criteria;
418
 
!     - nelts is set by the subroutine. It is an integer value
419
 
!       corresponding to the number of boundary faces verifying the
420
 
!       selection criteria;
421
 
!     - lstelt is set by the subroutine. It is an integer array of
422
 
!       size nelts containing the list of boundary faces verifying
423
 
!       the selection criteria.
424
 
!
425
 
!     string may contain:
426
 
!     - references to colors (ex.: 1, 8, 26, ...
427
 
!     - references to groups (ex.: inlet, group1, ...)
428
 
!     - geometric criteria (ex. x < 0.1, y >= 0.25, ...)
429
 
!
430
 
!     These criteria may be combined using logical operators
431
 
!     ('and', 'or') and parentheses.
432
 
!-------------------------------------------------------------------------------
433
 
! Arguments
434
 
!__________________.____._____.________________________________________________.
435
 
! name             !type!mode ! role                                           !
436
 
!__________________!____!_____!________________________________________________!
437
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
438
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
439
 
! ndim             ! i  ! <-- ! spatial dimension                              !
440
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
441
 
! ncel             ! i  ! <-- ! number of cells                                !
442
 
! nfac             ! i  ! <-- ! number of interior faces                       !
443
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
444
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
445
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
446
 
! nnod             ! i  ! <-- ! number of vertices                             !
447
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
448
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
449
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
450
 
! nvar             ! i  ! <-- ! total number of variables                      !
451
 
! nscal            ! i  ! <-- ! total number of scalars                        !
452
 
! nphas            ! i  ! <-- ! number of phases                               !
453
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
454
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
455
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
456
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
457
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
458
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
459
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
460
 
!  (nfml, nprfml)  !    !     !                                                !
461
 
! maxelt           !  e ! <-- ! max number of cells and faces (int/boundary)   !
462
 
! lstelt(maxelt)   ! ia ! --- ! work array                                     !
463
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
464
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
465
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
466
 
! nodfac(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
467
 
! icodcl           ! ia ! --> ! boundary condition code                        !
468
 
!  (nfabor, nvar)  !    !     ! = 1  -> Dirichlet                              !
469
 
!                  !    !     ! = 2  -> flux density                           !
470
 
!                  !    !     ! = 4  -> sliding wall and u.n=0 (velocity)      !
471
 
!                  !    !     ! = 5  -> friction and u.n=0 (velocity)          !
472
 
!                  !    !     ! = 6  -> roughness and u.n=0 (velocity)         !
473
 
!                  !    !     ! = 9  -> free inlet/outlet (velocity)           !
474
 
!                  !    !     !         inflowing possibly blocked             !
475
 
! itrifb(nfabor    ! ia ! <-- ! indirection for boundary faces ordering)       !
476
 
!  (nfabor, nphas) !    !     !                                                !
477
 
! itypfb           ! ia ! --> ! boundary face types                            !
478
 
!  (nfabor, nphas) !    !     !                                                !
479
 
! idevel(nideve)   ! ia ! <-- ! integer work array for temporary developpement !
480
 
! ituser(nituse    ! ia ! <-- ! user-reserved integer work array               !
481
 
! ia(*)            ! ia ! --- ! main integer work array                        !
482
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
483
 
!  (ndim, ncelet)  !    !     !                                                !
484
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
485
 
!  (ndim, nfac)    !    !     !                                                !
486
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
487
 
!  (ndim, nfavor)  !    !     !                                                !
488
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
489
 
!  (ndim, nfac)    !    !     !                                                !
490
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
491
 
!  (ndim, nfabor)  !    !     !                                                !
492
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
493
 
!  (ndim, nnod)    !    !     !                                                !
494
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
495
 
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
496
 
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
497
 
!  (ncelet, *)     !    !     !  (at current and preceding time steps)         !
498
 
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
499
 
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
500
 
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
501
 
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
502
 
!  (nfabor, *)     !    !     !                                                !
503
 
! rcodcl           ! ra ! --> ! boundary condition values                      !
504
 
!                  !    !     ! rcodcl(1) = Dirichlet value                    !
505
 
!                  !    !     ! rcodcl(2) = exterior exchange coefficient      !
506
 
!                  !    !     !  (infinite if no exchange)                     !
507
 
!                  !    !     ! rcodcl(3) = flux density value                 !
508
 
!                  !    !     !  (negative for gain) in w/m2 or                !
509
 
!                  !    !     !  roughness height (m) if icodcl=6              !
510
 
!                  !    !     ! for velocities           ( vistl+visct)*gradu  !
511
 
!                  !    !     ! for pressure                         dt*gradp  !
512
 
!                  !    !     ! for scalars    cp*(viscls+visct/sigmas)*gradt  !
513
 
! w1,2,3,4,5,6     ! ra ! --- ! work arrays                                    !
514
 
!  (ncelet)        !    !     !  (computation of pressure gradient)            !
515
 
! coefu            ! ra ! --- ! tab de trav                                    !
516
 
!  (nfabor, 3)     !    !     !  (computation of pressure gradient)            !
517
 
! rdevel(nrdeve)   ! ra ! <-> ! tab reel complementaire developemt             !
518
 
! rdevel(nideve)   ! ra ! <-- ! real work array for temporary developpement    !
519
 
! rtuser(nituse    ! ra ! <-- ! user-reserved real work array                  !
520
 
! ra(*)            ! ra ! --- ! main real work array                           !
521
 
!__________________!____!_____!________________________________________________!
522
 
 
523
 
!     Type: i (integer), r (real), s (string), a (array), l (logical),
524
 
!           and composite types (ex: ra real array)
525
 
!     mode: <-- input, --> output, <-> modifies data, --- work array
526
 
!===============================================================================
527
 
implicit none
528
 
 
529
 
!===============================================================================
530
 
!     Common blocks
531
 
!===============================================================================
532
 
 
533
 
include "ihmpre.h"
534
 
include "paramx.h"
535
 
include "pointe.h"
536
 
include "numvar.h"
537
 
include "optcal.h"
538
 
include "cstphy.h"
539
 
include "cstnum.h"
540
 
include "entsor.h"
541
 
include "parall.h"
542
 
include "period.h"
543
 
include "ppppar.h"
544
 
include "ppthch.h"
545
 
include "coincl.h"
546
 
include "cpincl.h"
547
 
include "fuincl.h"
548
 
include "ppincl.h"
549
 
 
550
 
!===============================================================================
551
 
 
552
 
! Arguments
553
 
 
554
 
integer          idbia0 , idbra0
555
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
556
 
integer          nfml   , nprfml
557
 
integer          nnod   , lndfac , lndfbr , ncelbr
558
 
integer          nvar   , nscal  , nphas
559
 
integer          nideve , nrdeve , nituse , nrtuse
560
 
 
561
 
integer          ifacel(2,nfac) , ifabor(nfabor)
562
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
563
 
integer          iprfml(nfml,nprfml)
564
 
integer          maxelt, lstelt(maxelt)
565
 
integer          ipnfac(nfac+1), nodfac(lndfac)
566
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
567
 
integer          icodcl(nfabor,nvar)
568
 
integer          itrifb(nfabor,nphas), itypfb(nfabor,nphas)
569
 
integer          izfppp(nfabor)
570
 
integer          idevel(nideve), ituser(nituse), ia(*)
571
 
 
572
 
double precision xyzcen(ndim,ncelet)
573
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
574
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
575
 
double precision xyznod(ndim,nnod), volume(ncelet)
576
 
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
577
 
double precision propce(ncelet,*)
578
 
double precision propfa(nfac,*), propfb(nfabor,*)
579
 
double precision coefa(nfabor,*), coefb(nfabor,*)
580
 
double precision rcodcl(nfabor,nvar,3)
581
 
double precision w1(ncelet),w2(ncelet),w3(ncelet)
582
 
double precision w4(ncelet),w5(ncelet),w6(ncelet)
583
 
double precision coefu(nfabor,ndim)
584
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
585
 
 
586
 
! LOCAL VARIABLES
587
 
 
588
 
 
589
 
integer          idebia, idebra
590
 
integer          ifac, iphas, ii
591
 
integer          izone
592
 
integer          iclafu
593
 
integer          ilelt, nlelt
594
 
 
595
 
double precision uref2, d2s3
596
 
double precision xkent, xeent
597
 
 
598
 
!===============================================================================
599
 
 
600
 
! TEST_TO_REMOVE_FOR_USE_OF_SUBROUTINE_START
601
 
!===============================================================================
602
 
! 0.  THIS TEST CERTIFY THIS VERY ROUINE IS USED
603
 
!     IN PLACE OF LIBRARY'S ONE
604
 
!===============================================================================
605
 
 
606
 
  if(1.eq.1) then
607
 
    write(nfecra,9001)
608
 
    call csexit (1)
609
 
    !==========
610
 
  endif
611
 
 
612
 
 9001 format(                                                           &
613
 
'@                                                            ',/,&
614
 
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
615
 
'@                                                            ',/,&
616
 
'@ @@ Beware    : Stop during Boundary Conditions Allocation  ',/,&
617
 
'@    =========                                               ',/,&
618
 
'@    FUEL                                                    ',/,&
619
 
'@     User subroutine USFUCL must be completed               ',/, &
620
 
'@                                                            ',/,&
621
 
'@                                                            ',/,&
622
 
'@  Computation will be stopped                               ',/,&
623
 
'@                                                            ',/,&
624
 
'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
625
 
'@                                                            ',/)
626
 
 
627
 
! TEST_TO_REMOVE_FOR_USE_OF_SUBROUTINE_END
628
 
 
629
 
 
630
 
!===============================================================================
631
 
! 1.  INITIALISATIONS
632
 
 
633
 
!===============================================================================
634
 
 
635
 
idebia = idbia0
636
 
idebra = idbra0
637
 
 
638
 
d2s3 = 2.d0/3.d0
639
 
!===============================================================================
640
 
! 2.  ALLOCATION OF BOUNDARY CONDITIONS VECTOR
641
 
!       LOOP ON BOUNDARY FACE
642
 
!         FAMILY AND PROPERTIES ARE DETERMINED
643
 
!         BOUNDARY CONDITION ALLOCATED
644
 
!
645
 
!     BOUNDARY CONDITIONS ON BONDARY FACE HAVE TO BE ALLOCATED HERE
646
 
!
647
 
!     USER'S WORK TO DO
648
 
!
649
 
!===============================================================================
650
 
 
651
 
!  Each kind of condition for extended physic is allocated with a number
652
 
!  IZONE ( 0<IZONE<= NOZPPM ; NOZPPM allocated in ppppar.h)
653
 
 
654
 
iphas = 1
655
 
 
656
 
! ---- The 12 color is a pure air inlet
657
 
 
658
 
CALL GETFBR('12',NLELT,LSTELT)
659
 
!==========
660
 
 
661
 
do ilelt = 1, nlelt
662
 
 
663
 
  ifac = lstelt(ilelt)
664
 
 
665
 
!   kind of boundary condition for standard variables
666
 
  itypfb(ifac,iphas) = ientre
667
 
 
668
 
!   Zone number allocation
669
 
  izone = 1
670
 
 
671
 
!   Zone number storage
672
 
  izfppp(ifac) = izone
673
 
 
674
 
! ------ This inlet have a fixed mass flux
675
 
 
676
 
  ientat(izone) = 1
677
 
  iqimp(izone)  = 1
678
 
!     - Air mass flow rate kg/s
679
 
  qimpat(izone) = 1.46d-03
680
 
!     - Air inlet temperature K
681
 
  timpat(izone) = 400.d0 + tkelvi
682
 
!     - Fuel mass flow rate kg/s
683
 
  qimpfl(izone) = 0.d0
684
 
 
685
 
! ----- The 12 color is a fixed mass flow rate inlet
686
 
!        user gives only the speed vector direction
687
 
!        (spedd vector norm is irrelevent)
688
 
!
689
 
  rcodcl(ifac,iu(iphas),1) = 0.d0
690
 
  rcodcl(ifac,iv(iphas),1) = 0.d0
691
 
  rcodcl(ifac,iw(iphas),1) = 5.d0
692
 
 
693
 
!   Boundary conditions of turbulence
694
 
  icalke(izone) = 1
695
 
!
696
 
!    - If ICALKE = 0 the boundary conditions of turbulence at
697
 
!      the inlet are calculated as follows:
698
 
 
699
 
  if(icalke(izone).eq.0) then
700
 
 
701
 
    uref2 = rcodcl(ifac,iu(iphas),1)**2                           &
702
 
           +rcodcl(ifac,iv(iphas),1)**2                           &
703
 
           +rcodcl(ifac,iw(iphas),1)**2
704
 
    uref2 = max(uref2,1.d-12)
705
 
    xkent  = epzero
706
 
    xeent  = epzero
707
 
 
708
 
    call keenin                                                   &
709
 
    !==========
710
 
      ( uref2, xintur(izone), dh(izone), cmu, xkappa,             &
711
 
        xkent, xeent )
712
 
 
713
 
    if    (itytur(iphas).eq.2) then
714
 
 
715
 
      rcodcl(ifac,ik(iphas),1)  = xkent
716
 
      rcodcl(ifac,iep(iphas),1) = xeent
717
 
 
718
 
    elseif(itytur(iphas).eq.3) then
719
 
 
720
 
      rcodcl(ifac,ir11(iphas),1) = d2s3*xkent
721
 
      rcodcl(ifac,ir22(iphas),1) = d2s3*xkent
722
 
      rcodcl(ifac,ir33(iphas),1) = d2s3*xkent
723
 
      rcodcl(ifac,ir12(iphas),1) = 0.d0
724
 
      rcodcl(ifac,ir13(iphas),1) = 0.d0
725
 
      rcodcl(ifac,ir23(iphas),1) = 0.d0
726
 
      rcodcl(ifac,iep(iphas),1)  = xeent
727
 
 
728
 
    elseif (iturb(iphas).eq.50) then
729
 
 
730
 
      rcodcl(ifac,ik(iphas),1)   = xkent
731
 
      rcodcl(ifac,iep(iphas),1)  = xeent
732
 
      rcodcl(ifac,iphi(iphas),1) = d2s3
733
 
      rcodcl(ifac,ifb(iphas),1)  = 0.d0
734
 
 
735
 
    elseif (iturb(iphas).eq.60) then
736
 
 
737
 
      rcodcl(ifac,ik(iphas),1)   = xkent
738
 
      rcodcl(ifac,iomg(iphas),1) = xeent/cmu/xkent
739
 
 
740
 
    endif
741
 
 
742
 
  endif
743
 
!
744
 
!    - If ICALKE = 1 the boundary conditions of turbulence at
745
 
!      the inlet refer to both, a hydraulic diameter and a
746
 
!      reference velocity given in usini1.f90.
747
 
!
748
 
  dh(izone)     = 0.032d0
749
 
!
750
 
!    - If ICALKE = 2 the boundary conditions of turbulence at
751
 
!      the inlet refer to a turbulence intensity.
752
 
!
753
 
  xintur(izone) = 0.d0
754
 
 
755
 
! ------ Treatment of user's scalars
756
 
 
757
 
  if ( (nscal-nscapp).gt.0 ) then
758
 
    do ii = 1, (nscal-nscapp)
759
 
      if(iphsca(ii).eq.iphas) then
760
 
        rcodcl(ifac,isca(ii),1) = 1.d0
761
 
      endif
762
 
    enddo
763
 
  endif
764
 
 
765
 
enddo
766
 
 
767
 
! ---- Inlet of both primary Air and Fuel
768
 
 
769
 
CALL GETFBR('11',NLELT,LSTELT)
770
 
!==========
771
 
 
772
 
do ilelt = 1, nlelt
773
 
 
774
 
  ifac = lstelt(ilelt)
775
 
 
776
 
!   kind of boundary condition for standard variables
777
 
  itypfb(ifac,iphas) = ientre
778
 
 
779
 
!   Zone number allocation
780
 
  izone = 2
781
 
 
782
 
!   Zone number storage
783
 
  izfppp(ifac) = izone
784
 
 
785
 
! ------ This inlet have a fixed mass flux
786
 
 
787
 
  ientfl(izone) = 1
788
 
  iqimp(izone)  = 1
789
 
!     - Air mass flow rate in kg/s
790
 
  qimpat(izone) = 1.46d-03
791
 
!     - Air Temperature at inlet in K
792
 
  timpat(izone) = 800.d0  + tkelvi
793
 
 
794
 
!     - Fuel mass flow rate in kg/s
795
 
  qimpfl(izone) = 1.46d-04/360.d0
796
 
 
797
 
!     - PERCENTAGE mass fraction of each granulometric class
798
 
!       ICLAFU (1 < ICLAFU < NCLAFU )
799
 
  iclafu = 1
800
 
  distfu(izone,iclafu) = 100.d0
801
 
 
802
 
!     - Fuel Temperature at inlet in K
803
 
  timpfl(izone) = 100.d0  + tkelvi
804
 
 
805
 
! ----- The 11 color is a fixed mass flow rate inlet
806
 
!        user gives only the speed vector direction
807
 
!        (spedd vector norm is irrelevent)
808
 
!
809
 
 
810
 
  rcodcl(ifac,iu(iphas),1) = 0.d0
811
 
  rcodcl(ifac,iv(iphas),1) = 0.d0
812
 
  rcodcl(ifac,iw(iphas),1) = 5.d0
813
 
 
814
 
!   Boundary conditions of turbulence
815
 
  icalke(izone) = 1
816
 
!
817
 
!    - If ICALKE = 0 the boundary conditions of turbulence at
818
 
!      the inlet are calculated as follows:
819
 
 
820
 
  if(icalke(izone).eq.0) then
821
 
 
822
 
    uref2 = rcodcl(ifac,iu(iphas),1)**2                           &
823
 
           +rcodcl(ifac,iv(iphas),1)**2                           &
824
 
           +rcodcl(ifac,iw(iphas),1)**2
825
 
    uref2 = max(uref2,1.d-12)
826
 
    xkent  = epzero
827
 
    xeent  = epzero
828
 
 
829
 
    call keenin                                                   &
830
 
    !==========
831
 
      ( uref2, xintur(izone), dh(izone), cmu, xkappa,             &
832
 
        xkent, xeent )
833
 
 
834
 
    if    (itytur(iphas).eq.2) then
835
 
 
836
 
      rcodcl(ifac,ik(iphas),1)  = xkent
837
 
      rcodcl(ifac,iep(iphas),1) = xeent
838
 
 
839
 
    elseif(itytur(iphas).eq.3) then
840
 
 
841
 
      rcodcl(ifac,ir11(iphas),1) = d2s3*xkent
842
 
      rcodcl(ifac,ir22(iphas),1) = d2s3*xkent
843
 
      rcodcl(ifac,ir33(iphas),1) = d2s3*xkent
844
 
      rcodcl(ifac,ir12(iphas),1) = 0.d0
845
 
      rcodcl(ifac,ir13(iphas),1) = 0.d0
846
 
      rcodcl(ifac,ir23(iphas),1) = 0.d0
847
 
      rcodcl(ifac,iep(iphas),1)  = xeent
848
 
 
849
 
    elseif (iturb(iphas).eq.50) then
850
 
 
851
 
      rcodcl(ifac,ik(iphas),1)   = xkent
852
 
      rcodcl(ifac,iep(iphas),1)  = xeent
853
 
      rcodcl(ifac,iphi(iphas),1) = d2s3
854
 
      rcodcl(ifac,ifb(iphas),1)  = 0.d0
855
 
 
856
 
    elseif (iturb(iphas).eq.60) then
857
 
 
858
 
      rcodcl(ifac,ik(iphas),1)   = xkent
859
 
      rcodcl(ifac,iomg(iphas),1) = xeent/cmu/xkent
860
 
 
861
 
    endif
862
 
 
863
 
  endif
864
 
!
865
 
!    - If ICALKE = 1 the boundary conditions of turbulence at
866
 
!      the inlet refer to both, a hydraulic diameter and a
867
 
!      reference velocity given in usini1.f90.
868
 
!
869
 
  dh(izone)     = 0.032d0
870
 
!
871
 
!    - If ICALKE = 2 the boundary conditions of turbulence at
872
 
!      the inlet refer to a turbulence intensity.
873
 
!
874
 
  xintur(izone) = 0.d0
875
 
enddo
876
 
 
877
 
! --- Color 15 is a wall
878
 
 
879
 
CALL GETFBR('15',NLELT,LSTELT)
880
 
!==========
881
 
 
882
 
do ilelt = 1, nlelt
883
 
 
884
 
  ifac = lstelt(ilelt)
885
 
 
886
 
!          WALL  : nul mass flow rate (nul pressure flux)
887
 
!                  rubbing for speed (and turbulence)
888
 
!                  nul scalar fluxes
889
 
 
890
 
!   kind of boundary condition for standard variables
891
 
  itypfb(ifac,iphas)   = iparoi
892
 
 
893
 
 
894
 
!   Zone number allocation
895
 
  izone = 3
896
 
 
897
 
!   Zone number storage
898
 
  izfppp(ifac) = izone
899
 
 
900
 
enddo
901
 
 
902
 
! --- Color 19 is an outlet
903
 
 
904
 
CALL GETFBR('19',NLELT,LSTELT)
905
 
!==========
906
 
 
907
 
do ilelt = 1, nlelt
908
 
 
909
 
  ifac = lstelt(ilelt)
910
 
 
911
 
!          OUTLET : nul fluxes for speed and scalar
912
 
!                   pressure fixed
913
 
 
914
 
!   kind of boundary condition for standard variables
915
 
    itypfb(ifac,iphas)   = isolib
916
 
 
917
 
!   Zone number allocation
918
 
    izone = 4
919
 
 
920
 
!   Zone number storage
921
 
    izfppp(ifac) = izone
922
 
 
923
 
  enddo
924
 
 
925
 
! --- 14 and 4 are symetry
926
 
 
927
 
CALL GETFBR('14 or 4',NLELT,LSTELT)
928
 
!==========
929
 
 
930
 
do ilelt = 1, nlelt
931
 
 
932
 
  ifac = lstelt(ilelt)
933
 
 
934
 
!          SYMETRY
935
 
 
936
 
!   kind of boundary condition for standard variables
937
 
  itypfb(ifac,iphas)   = isymet
938
 
 
939
 
!   Zone number allocation
940
 
  izone = 5
941
 
 
942
 
!   Zone number storage
943
 
  izfppp(ifac) = izone
944
 
 
945
 
enddo
946
 
 
947
 
 
948
 
!----
949
 
! END
950
 
!----
951
 
 
952
 
return
953
 
end subroutine