1!-------------------------------------------------------------------------------
2
3! This file is part of Code_Saturne, a general-purpose CFD tool.
4!
5! Copyright (C) 1998-2021 EDF S.A.
6!
7! This program is free software; you can redistribute it and/or modify it under
8! the terms of the GNU General Public License as published by the Free Software
9! Foundation; either version 2 of the License, or (at your option) any later
10! version.
11!
12! This program is distributed in the hope that it will be useful, but WITHOUT
13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
15! details.
16!
17! You should have received a copy of the GNU General Public License along with
18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
19! Street, Fifth Floor, Boston, MA 02110-1301, USA.
20
21!-------------------------------------------------------------------------------
22
23!===============================================================================
24! Function :
25! --------
26
27!> \file d3ptcl.f90
28!> \brief Automatic boundary conditions for 3 PTHEM gas diffusion flame model
29
30!-------------------------------------------------------------------------------
31! Arguments
32!______________________________________________________________________________.
33!  mode           name          role                                           !
34!______________________________________________________________________________!
35!> \param[in]     nvar          total number of variables
36!> \param[in]     itypfb        boundary face types
37!> \param[out]    izfppp        boundary face zone number
38!> \param[out]    icodcl        boundary condition code:
39!>                               - 1 Dirichlet
40!>                               - 2 Radiative outlet
41!>                               - 3 Neumann
42!>                               - 4 sliding and
43!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
44!>                               - 5 smooth wall and
45!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
46!>                               - 6 rough wall and
47!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
48!>                               - 9 free inlet/outlet
49!>                                 (input mass flux blocked to 0)
50!> \param[out]    rcodcl        boundary condition values:
51!>                               - rcodcl(1) value of the dirichlet
52!>                               - rcodcl(2) value of the exterior exchange
53!>                                 coefficient (infinite if no exchange)
54!>                               - rcodcl(3) value flux density
55!>                                 (negative if gain) in w/m2 or roughness
56!>                                 in m if icodcl=6
57!>                                 -# for the velocity \f$ (\mu+\mu_T)
58!>                                    \gradt \, \vect{u} \cdot \vect{n}  \f$
59!>                                 -# for the pressure \f$ \Delta t
60!>                                    \grad P \cdot \vect{n}  \f$
61!>                                 -# for a scalar \f$ cp \left( K +
62!>                                     \dfrac{K_T}{\sigma_T} \right)
63!>                                     \grad T \cdot \vect{n} \f$
64!_______________________________________________________________________________
65
66subroutine d3ptcl &
67 ( itypfb , izfppp , icodcl,                                             &
68   rcodcl )
69
70!===============================================================================
71
72!===============================================================================
73! Module files
74!===============================================================================
75
76use paramx
77use dimens, only : nscal, nvar
78use numvar
79use optcal
80use cstphy
81use cstnum
82use entsor
83use parall
84use ppppar
85use ppthch
86use coincl
87use cpincl
88use ppincl
89use mesh
90use field
91use cs_c_bindings
92
93!===============================================================================
94
95implicit none
96
97! Arguments
98
99integer          itypfb(nfabor)
100integer          izfppp(nfabor)
101integer          icodcl(nfabor,nvar)
102
103double precision rcodcl(nfabor,nvar,3)
104
105! Local variables
106
107integer          igg, ifac, izone, mode, iscal
108integer          ii, iel, ifue, ioxy
109integer          icke
110double precision qimabs, qisqc, viscla, d2s3, uref2, rhomoy, dhy, xiturb
111double precision ustar2, xkent, xeent
112double precision qcalc(nozppm)
113double precision coefg(ngazgm)
114double precision, dimension(:), pointer ::  brom
115double precision, dimension(:), pointer :: viscl
116
117!===============================================================================
118! 0.  INITIALISATIONS
119!===============================================================================
120
121call field_get_val_s(ibrom, brom)
122call field_get_val_s(iviscl, viscl)
123
124d2s3 = 2.d0/3.d0
125
126do igg = 1, ngazgm
127  coefg(igg) = zero
128enddo
129
130!===============================================================================
131! 1.  ECHANGES EN PARALLELE POUR LES DONNEES UTILISATEUR
132!===============================================================================
133
134!  En realite on pourrait eviter cet echange en modifiant usd3pc et en
135!    demandant a l'utilisateur de donner les grandeurs dependant de la
136!    zone hors de la boucle sur les faces de bord : les grandeurs
137!    seraient ainsi disponibles sur tous les processeurs. Cependant,
138!    ca rend le sous programme utilisateur un peu plus complique et
139!    surtout, si l'utilisateur le modifie de travers, ca ne marche pas.
140!  On suppose que toutes les gandeurs fournies sont positives, ce qui
141!    permet d'utiliser un max pour que tous les procs les connaissent.
142!    Si ce n'est pas le cas, c'est plus complique mais on pourrait
143!    s'en tirer avec un max quand meme, sauf qimp qui peut etre negatif.
144
145if(irangp.ge.0) then
146  call parmax(tinfue)
147  call parmax(tinoxy)
148  ! qimabs  = |qimp|, or 0 if no boundary faces on the rank
149  do ii = 1, nozapm
150    qimabs = abs(qimp(ii))
151    ! the rank owning the max of qimabs, meaning |qimp|, share
152    ! qimp with the others
153    call parmxl(1,qimabs,qimp(ii))
154  enddo
155  call parimx(nozapm,iqimp )
156  call parimx(nozapm,ientox)
157  call parimx(nozapm,ientfu)
158endif
159
160!===============================================================================
161! 2.  REMPLISSAGE DU TABLEAU DES CONDITIONS LIMITES POUR LES SCALAIRES
162!
163!       ON BOUCLE SUR TOUTES LES FACES D'ENTREE
164!                     =========================
165
166!       ON DETERMINE LA FAMILLE ET SES PROPRIETES
167!       ON IMPOSE LES CONDITIONS AUX LIMITES
168!===============================================================================
169
170! Fuel and Oxidant inlet are checked
171ifue = 0
172ioxy = 0
173do ii = 1, nzfppp
174  izone = ilzppp(ii)
175  if (ientfu(izone).eq.1) then
176    ifue = 1
177  elseif(ientox(izone).eq.1) then
178    ioxy = 1
179  endif
180enddo
181if(irangp.ge.0) then
182  call parcmx(ifue)
183  call parcmx(ioxy)
184endif
185
186! Fuel inlet at TINFUE: HINFUE calculation
187if (ifue.eq.1) then
188  coefg(1) = 1.d0
189  coefg(2) = zero
190  coefg(3) = zero
191  mode    = -1
192  call cothht                                                   &
193  !==========
194      ( mode   , ngazg , ngazgm  , coefg  ,                     &
195        npo    , npot   , th     , ehgazg ,                     &
196        hinfue , tinfue )
197endif
198
199! Oxidant inlet at TINOXY: HINOXY calculation
200if (ioxy.eq.1) then
201  coefg(1) = zero
202  coefg(2) = 1.d0
203  coefg(3) = zero
204  mode    = -1
205  call cothht                                                   &
206  !==========
207      ( mode   , ngazg , ngazgm  , coefg  ,                     &
208        npo    , npot   , th     , ehgazg ,                     &
209        hinoxy , tinoxy  )
210endif
211
212
213do ifac = 1, nfabor
214
215  izone = izfppp(ifac)
216
217  if (itypfb(ifac).eq.ientre.or.itypfb(ifac).eq.i_convective_inlet) then
218
219    ! Fuel inlet at TINFUE
220    if (ientfu(izone).eq.1) then
221
222      ! Neumann for outflow
223      if (qimp(izone).lt.0.d0) then
224
225        do iscal = 1, nscal
226          icodcl(ifac,isca(iscal)) = 3
227          rcodcl(ifac,isca(iscal),3) = 0.d0
228        enddo
229
230      ! Dirichlet for inflow
231      else
232        ! Mean mixture fraction
233        rcodcl(ifac,isca(ifm),1)   = 1.d0
234
235        ! Mixture fraction variance
236        rcodcl(ifac,isca(ifp2m),1) = 0.d0
237
238        ! Mixture enthalpy
239        if (ippmod(icod3p).eq.1) then
240          rcodcl(ifac,isca(iscalt),1) = hinfue
241        endif
242
243        ! Soot model
244        if (isoot.ge.1) then
245          rcodcl(ifac,isca(ifsm),1) = 0.d0
246          rcodcl(ifac,isca(inpm),1) = 0.d0
247        endif
248
249      endif
250
251      ! Density
252      brom(ifac) = pther/(cs_physical_constants_r*tinfue/wmolg(1))
253
254    ! Oxydant inlet at TINOXY
255    elseif (ientox(izone).eq.1) then
256
257      ! Neumann for outflow
258      if (qimp(izone).lt.0.d0) then
259
260        do iscal = 1, nscal
261          icodcl(ifac,isca(iscal)) = 3
262          rcodcl(ifac,isca(iscal),3) = 0.d0
263        enddo
264
265      ! Dirichlet for inflow
266      else
267
268        ! Mean mixture fraction
269        rcodcl(ifac,isca(ifm),1)   = 0.d0
270
271        ! Mixture fraction variance
272        rcodcl(ifac,isca(ifp2m),1) = 0.d0
273
274        ! Mixture enthalpy
275        if ( ippmod(icod3p).eq.1 ) then
276          rcodcl(ifac,isca(iscalt),1) = hinoxy
277        endif
278
279        ! Soot model
280        if (isoot.ge.1) then
281          rcodcl(ifac,isca(ifsm),1) = 0.d0
282          rcodcl(ifac,isca(inpm),1) = 0.d0
283        endif
284
285        ! Density
286        brom(ifac) = pther/(cs_physical_constants_r*tinoxy/wmolg(2))
287
288      endif
289
290    endif
291
292  endif
293
294enddo
295
296!===============================================================================
297! 3.  SI IQIMP = 1 : CORRECTION DES VITESSES (EN NORME) POUR CONTROLER
298!                    LES DEBITS IMPOSES
299!     SI IQIMP = 0 : CALCUL DE QIMP
300
301!     ON BOUCLE SUR TOUTES LES FACES D'ENTREE
302!                   =========================
303!===============================================================================
304
305! --- Debit calcule
306
307do izone = 1, nozppm
308  qcalc(izone) = 0.d0
309enddo
310do ifac = 1, nfabor
311  izone = izfppp(ifac)
312  qcalc(izone) = qcalc(izone) - brom(ifac) *             &
313     ( rcodcl(ifac,iu,1)*surfbo(1,ifac) +                  &
314       rcodcl(ifac,iv,1)*surfbo(2,ifac) +                  &
315       rcodcl(ifac,iw,1)*surfbo(3,ifac) )
316enddo
317
318if(irangp.ge.0) then
319  call parrsm(nozapm,qcalc)
320endif
321
322do izone = 1, nozapm
323  if (iqimp(izone).eq.0) then
324    qimp(izone) = qcalc(izone)
325  endif
326enddo
327
328! --- Correction des vitesses en norme
329
330do ifac = 1, nfabor
331  izone = izfppp(ifac)
332  if (iqimp(izone).eq.1) then
333    if (abs(qcalc(izone)).gt.epzero) then
334      qisqc = qimp(izone)/qcalc(izone)
335    else
336      qisqc = 0.d0
337    endif
338    rcodcl(ifac,iu,1) = rcodcl(ifac,iu,1)*qisqc
339    rcodcl(ifac,iv,1) = rcodcl(ifac,iv,1)*qisqc
340    rcodcl(ifac,iw,1) = rcodcl(ifac,iw,1)*qisqc
341  endif
342enddo
343
344!===============================================================================
345! 4.  REMPLISSAGE DU TABLEAU DES CONDITIONS LIMITES POUR LA TURBULENCE
346!       ON BOUCLE SUR TOUTES LES FACES D'ENTREE
347!                     =========================
348!         ON DETERMINE LA FAMILLE ET SES PROPRIETES
349!         ON IMPOSE LES CONDITIONS AUX LIMITES POUR LA TURBULENCE
350!         (pour n'importe quel modele)
351!===============================================================================
352
353do ifac = 1, nfabor
354
355  izone = izfppp(ifac)
356
357  if (itypfb(ifac).eq.ientre.or.itypfb(ifac).eq.i_convective_inlet) then
358
359    ! Neumann for outflow
360    if (qimp(izone).lt.0.d0) then
361
362      if (itytur.eq.2) then
363
364        icodcl(ifac,ik ) = 3
365        icodcl(ifac,iep) = 3
366        rcodcl(ifac,ik ,3) = 0.d0
367        rcodcl(ifac,iep,3) = 0.d0
368
369      elseif (itytur.eq.3) then
370
371        icodcl(ifac,ir11) = 3
372        icodcl(ifac,ir22) = 3
373        icodcl(ifac,ir33) = 3
374        icodcl(ifac,ir12) = 3
375        icodcl(ifac,ir13) = 3
376        icodcl(ifac,ir23) = 3
377        icodcl(ifac,iep ) = 3
378        rcodcl(ifac,ir11,3) = 0.d0
379        rcodcl(ifac,ir22,3) = 0.d0
380        rcodcl(ifac,ir33,3) = 0.d0
381        rcodcl(ifac,ir12,3) = 0.d0
382        rcodcl(ifac,ir13,3) = 0.d0
383        rcodcl(ifac,ir23,3) = 0.d0
384        rcodcl(ifac,iep, 3) = 0.d0
385
386      elseif (iturb.eq.50) then
387
388        icodcl(ifac,ik  ) = 3
389        icodcl(ifac,iep ) = 3
390        icodcl(ifac,iphi) = 3
391        icodcl(ifac,ifb ) = 3
392        rcodcl(ifac,ik,  3) = 0.d0
393        rcodcl(ifac,iep, 3) = 0.d0
394        rcodcl(ifac,iphi,3) = 0.d0
395        rcodcl(ifac,ifb, 3) = 0.d0
396
397      elseif (iturb.eq.60) then
398
399        icodcl(ifac,ik  ) = 3
400        icodcl(ifac,iomg) = 3
401        rcodcl(ifac,ik,  3) = 0.d0
402        rcodcl(ifac,iomg,3) = 0.d0
403
404      elseif(iturb.eq.70) then
405
406        icodcl(ifac,inusa) = 3
407        rcodcl(ifac,inusa,3) = 0.d0
408
409      endif
410
411    ! Dirichlet for inflow
412    else
413
414      ! La turbulence est calculee par defaut si ICALKE different de 0
415      !    - soit a partir du diametre hydraulique, d'une vitesse
416      !      de reference adaptes a l'entree courante si ICALKE = 1
417      !    - soit a partir du diametre hydraulique, d'une vitesse
418      !      de reference et de l'intensite turvulente
419      !      adaptes a l'entree courante si ICALKE = 2
420
421      if ( icalke(izone).ne.0 ) then
422
423        uref2 = rcodcl(ifac,iu,1)**2                         &
424              + rcodcl(ifac,iv,1)**2                         &
425              + rcodcl(ifac,iw,1)**2
426        uref2 = max(uref2,epzero)
427        rhomoy = brom(ifac)
428        iel    = ifabor(ifac)
429        viscla = viscl(iel)
430        icke   = icalke(izone)
431        dhy    = dh(izone)
432        xiturb = xintur(izone)
433        ustar2 = 0.d0
434        xkent = epzero
435        xeent = epzero
436
437        if (icke.eq.1) then
438          !   Calculation of turbulent inlet conditions using
439          !     standard laws for a circular pipe
440          !     (their initialization is not needed here but is good practice).
441          call turbulence_bc_inlet_hyd_diam(ifac, uref2, dhy, rhomoy, viscla,  &
442                                            rcodcl)
443        else if (icke.eq.2) then
444
445          ! Calculation of turbulent inlet conditions using
446          !   the turbulence intensity and standard laws for a circular pipe
447          !   (their initialization is not needed here but is good practice)
448
449          call turbulence_bc_inlet_turb_intensity(ifac, uref2, xiturb, dhy,  &
450                                                  rcodcl)
451
452
453        endif
454
455      endif ! icalke
456
457    endif ! qimp
458
459  endif ! itypfb
460
461enddo
462
463!----
464! FORMATS
465!----
466
467
468!----
469! FIN
470!----
471
472return
473end subroutine
474