← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/comb/cs_fuel_scast.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-24 00:00:08 UTC
  • mfrom: (6.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20111124000008-2vo99e38267942q5
Tags: 2.1.0-3
Install a missing file

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/comb/cs_fuel_scast.f90
 
1
!-------------------------------------------------------------------------------
 
2
 
 
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
4
!
 
5
! Copyright (C) 1998-2011 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
subroutine cs_fuel_scast &
 
24
!=======================
 
25
 
 
26
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
 
27
   iscal  ,                                                       &
 
28
   itypfb ,                                                       &
 
29
   icepdc , icetsm , itypsm ,                                     &
 
30
   izfppp ,                                                       &
 
31
   dt     , rtpa   , rtp    , propce , propfa , propfb ,          &
 
32
   coefa  , coefb  , ckupdc , smacel ,                            &
 
33
   smbrs  , rovsdt )
 
34
 
 
35
!===============================================================================
 
36
! FONCTION :
 
37
! ----------
 
38
 
 
39
! ROUTINE PHYSIQUE PARTICULIERE : FLAMME FUEL
 
40
!   ON PRECISE LES TERMES SOURCES POUR UN SCALAIRE PP
 
41
!   SUR UN PAS DE TEMPS
 
42
 
 
43
 
 
44
! ATTENTION : LE TRAITEMENT DES TERMES SOURCES EST DIFFERENT
 
45
! ---------   DE CELUI DE USTSSC.F
 
46
 
 
47
! ON RESOUT ROVSDT*D(VAR) = SMBRS
 
48
 
 
49
! ROVSDT ET SMBRS CONTIENNENT DEJA D'EVENTUELS TERMES SOURCES
 
50
!  UTILISATEUR. IL FAUT DONC LES INCREMENTER ET PAS LES
 
51
!  ECRASER
 
52
 
 
53
! POUR DES QUESTIONS DE STABILITE, ON NE RAJOUTE DANS ROVSDT
 
54
!  QUE DES TERMES POSITIFS. IL N'Y A PAS DE CONTRAINTE POUR
 
55
!  SMBRS
 
56
 
 
57
! DANS LE CAS D'UN TERME SOURCE EN CEXP + CIMP*VAR ON DOIT
 
58
! ECRIRE :
 
59
!          SMBRS  = SMBRS  + CEXP + CIMP*VAR
 
60
!          ROVSDT = ROVSDT + MAX(-CIMP,ZERO)
 
61
 
 
62
! ON FOURNIT ICI ROVSDT ET SMBRS (ILS CONTIENNENT RHO*VOLUME)
 
63
!    SMBRS en kg variable/s :
 
64
!     ex : pour la vitesse            kg m/s2
 
65
!          pour les temperatures      kg degres/s
 
66
!          pour les enthalpies        Joules/s
 
67
!    ROVSDT en kg /s
 
68
 
 
69
!-------------------------------------------------------------------------------
 
70
!ARGU                             ARGUMENTS
 
71
!__________________.____._____.________________________________________________.
 
72
! name             !type!mode ! role                                           !
 
73
!__________________!____!_____!________________________________________________!
 
74
! nvar             ! i  ! <-- ! total number of variables                      !
 
75
! nscal            ! i  ! <-- ! total number of scalars                        !
 
76
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
 
77
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
 
78
! iscal            ! i  ! <-- ! scalar number                                  !
 
79
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
 
80
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
 
81
! itypsm           ! te ! <-- ! type de source de masse pour les               !
 
82
! (ncesmp,nvar)    !    !     !  variables (cf. ustsma)                        !
 
83
! izfppp           ! te ! --> ! numero de zone de la face de bord              !
 
84
! (nfabor)         !    !     !  pour le module phys. part.                    !
 
85
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
 
86
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
 
87
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
 
88
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
 
89
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
 
90
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
 
91
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
 
92
!  (nfabor, *)     !    !     !                                                !
 
93
! ckupdc           ! tr ! <-- ! tableau de travail pour pdc                    !
 
94
!  (ncepdp,6)      !    !     !                                                !
 
95
! smacel           ! tr ! <-- ! valeur des variables associee a la             !
 
96
! (ncesmp,*   )    !    !     !  source de masse                               !
 
97
!                  !    !     !  pour ivar=ipr, smacel=flux de masse           !
 
98
! smbrs(ncelet)    ! tr ! --> ! second membre explicite                        !
 
99
! rovsdt(ncelet    ! tr ! --> ! partie diagonale implicite                     !
 
100
!__________________!____!_____!________________________________________________!
 
101
 
 
102
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
103
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
104
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
105
!            --- tableau de travail
 
106
!===============================================================================
 
107
 
 
108
!===============================================================================
 
109
! Module files
 
110
!===============================================================================
 
111
 
 
112
use dimens, only: ndimfb
 
113
use paramx
 
114
use numvar
 
115
use entsor
 
116
use optcal
 
117
use cstphy
 
118
use cstnum
 
119
use parall
 
120
use period
 
121
use ppppar
 
122
use ppthch
 
123
use coincl
 
124
use cpincl
 
125
use cs_fuel_incl
 
126
use ppincl
 
127
use ppcpfu
 
128
use mesh
 
129
 
 
130
!===============================================================================
 
131
 
 
132
implicit none
 
133
 
 
134
! Arguments
 
135
 
 
136
integer          nvar   , nscal
 
137
integer          ncepdp , ncesmp
 
138
integer          iscal
 
139
 
 
140
integer          itypfb(nfabor)
 
141
integer          icepdc(ncepdp)
 
142
integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
 
143
integer          izfppp(nfabor)
 
144
 
 
145
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
 
146
double precision propce(ncelet,*)
 
147
double precision propfa(nfac,*), propfb(ndimfb,*)
 
148
double precision coefa(ndimfb,*), coefb(ndimfb,*)
 
149
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
 
150
double precision smbrs(ncelet), rovsdt(ncelet)
 
151
 
 
152
! Local variables
 
153
 
 
154
character*80     chaine
 
155
integer          ivar , ipcrom, iel, icla , numcla
 
156
integer          iexp1 , iexp2 , iexp3 , ifac
 
157
integer          iscala
 
158
integer          ipcro2 , ipcte1 , ipcte2
 
159
integer          ipcdia , ipcvst
 
160
integer          imode  , iesp
 
161
integer          ipcgev , ipcght , ipchgl
 
162
integer          itermx,nbpauv,nbrich,nbepau,nberic
 
163
integer          nbarre,nbimax,nbpass
 
164
 
 
165
double precision aux, rhovst
 
166
double precision t1, h1, t2, h2
 
167
double precision xng,rhofol , rom
 
168
double precision gameva,fdev,fsd,ftrac,hfov
 
169
double precision ho2,hco,xesp(ngazem),xcoke,t2mt1
 
170
double precision gmech,gmvap,gmhet
 
171
 
 
172
double precision xxco,xxo2,xxco2,xxh2o,xco2mx
 
173
double precision xkp,xk0p,xkm,xk0m,wplus,wmoins,t0p,t0m
 
174
double precision auxp,auxm, aux1 , aux2 , aux3 , w1
 
175
 
 
176
double precision xeq,anmr,xcot,xo2t,xco2e,xo2e,xcoe,tauchi,tautur
 
177
double precision sqh2o , x2 , wmhcn , wmno ,wmo2
 
178
 
 
179
integer          iterch
 
180
double precision err1mx,err2mx
 
181
 
 
182
double precision errch,ter1,ddelta,fn,qpr
 
183
double precision auxmax,auxmin
 
184
double precision ymoy,volm,volmp,dmp
 
185
!
 
186
double precision fn0,fn1,fn2,anmr0,anmr1,anmr2
 
187
double precision lnk0p,l10k0e,lnk0m,t0e,xco2eq,xcoeq,xo2eq
 
188
double precision xcom,xo2m,xkcequ,xkpequ,xden
 
189
!===============================================================================
 
190
! 1. INITIALISATION
 
191
!===============================================================================
 
192
 
 
193
! --- Numero du scalaire a traiter : ISCAL
 
194
 
 
195
! --- Numero de la variable associee au scalaire a traiter ISCAL
 
196
ivar = isca(iscal)
 
197
 
 
198
! --- Nom de la variable associee au scalaire a traiter ISCAL
 
199
chaine = nomvar(ipprtp(ivar))
 
200
 
 
201
! --- Numero des grandeurs physiques (voir usclim)
 
202
ipcrom = ipproc(irom)
 
203
ipcvst = ipproc(ivisct)
 
204
 
 
205
! --- Temperature phase gaz
 
206
 
 
207
ipcte1 = ipproc(itemp1)
 
208
 
 
209
!===============================================================================
 
210
! 2. PRISE EN COMPTE DES TERMES SOURCES POUR LES VARIABLES RELATIVES
 
211
!    AUX CLASSES DE PARTICULES
 
212
!===============================================================================
 
213
 
 
214
! --> Terme source pour l'enthalpie du liquide
 
215
 
 
216
if ( ivar .ge. isca(ih2(1))     .and.                            &
 
217
     ivar .le. isca(ih2(nclafu))      ) then
 
218
 
 
219
  if (iwarni(ivar).ge.1) then
 
220
    write(nfecra,1000) chaine(1:8)
 
221
  endif
 
222
 
 
223
  numcla = ivar-isca(ih2(1))+1
 
224
 
 
225
  ipcro2 = ipproc(irom2 (numcla))
 
226
  ipcdia = ipproc(idiam2(numcla))
 
227
  ipcte2 = ipproc(itemp2(numcla))
 
228
  ipcgev = ipproc(igmeva(numcla))
 
229
  ipcght = ipproc(igmhtf(numcla))
 
230
  ipchgl = ipproc(ih1hlf(numcla))
 
231
 
 
232
!       La variable est l'enthalpie du liquide ramen�e �
 
233
!       la masse de m�lange
 
234
!       Les flux interfaciaux contribuent � faire varier l'enthalpie
 
235
!       du liquide
 
236
!       La vapeur emporte son enthalpie
 
237
!       flux = PROPCE(IEL,IPPROC(IGMEVA))
 
238
!       enthalpie massique reconstitu�e � partir de EHGAZE(IFOV )
 
239
!       � la temp�rature de la goutte
 
240
!       L'oxydation h�terog�ne comporte un flux "entrant" d'O2
 
241
!       un flux sortant de CO
 
242
!       le flux net est celui du carbone
 
243
!       fluxIN  = 16/12 * PROPCE(IEL,IPPROC(IGMHTF))
 
244
!       fluxOUT = 28/12 * PROPCE(IEL,IPPROC(IGMHTF))
 
245
!       Enthalpie entrante reconstitu�e � partir de EHGAZE(IO2 )
 
246
!       � la temp�rature du gaz environnant
 
247
!       Enthalpie sortante reconstitu�e � partir de EHGAZE(ICO )
 
248
!       � la temp�rature du grain
 
249
 
 
250
  imode = -1
 
251
  do iel = 1, ncel
 
252
!
 
253
    if ( rtpa(iel,isca(iyfol(numcla))) .gt. epsifl ) then
 
254
!
 
255
      rom = propce(iel,ipcrom)
 
256
 
 
257
      do iesp = 1, ngazem
 
258
        xesp(iesp) = zero
 
259
      enddo
 
260
 
 
261
      xesp(ifov) = 1.d0
 
262
      call cs_fuel_htconvers1(imode,hfov,xesp,propce(iel,ipcte2))
 
263
 
 
264
      xesp(ifov) = zero
 
265
      xesp(io2)  = 1.d0
 
266
      call cs_fuel_htconvers1(imode,ho2 ,xesp,propce(iel,ipcte1))
 
267
 
 
268
      xesp(io2)  = zero
 
269
      xesp(ico)  = 1.d0
 
270
      call cs_fuel_htconvers1(imode,hco,xesp,propce(iel,ipcte2))
 
271
 
 
272
      t2mt1 = propce(iel,ipcte2)-propce(iel,ipcte1)
 
273
 
 
274
      gmech = -propce(iel,ipchgl)*t2mt1
 
275
      gmvap = propce(iel,ipcgev)*hfov*t2mt1
 
276
      gmhet = 16.d0/12.d0*propce(iel,ipcght)*ho2                    &
 
277
             -28.d0/12.d0*propce(iel,ipcght)*hco
 
278
 
 
279
      smbrs(iel) = smbrs(iel) +                                     &
 
280
           ( gmech+gmvap+gmhet )*rom*volume(iel)
 
281
      rhovst = ( propce(iel,ipchgl)                                 &
 
282
                -propce(iel,ipcgev)*hfov )/cp2fol                   &
 
283
              *rom*volume(iel)
 
284
      rovsdt(iel) = rovsdt(iel) +  max(zero,rhovst)
 
285
 
 
286
    endif
 
287
 
 
288
  enddo
 
289
 
 
290
! --> T.S. pour la masse de liquide
 
291
 
 
292
elseif ( ivar .ge. isca(iyfol(1))     .and.                       &
 
293
         ivar .le. isca(iyfol(nclafu))        ) then
 
294
 
 
295
  if (iwarni(ivar).ge.1) then
 
296
    write(nfecra,1000) chaine(1:8)
 
297
  endif
 
298
 
 
299
  numcla = ivar-isca(iyfol(1))+1
 
300
 
 
301
  ipcro2 = ipproc(irom2 (numcla))
 
302
  ipcdia = ipproc(idiam2(numcla))
 
303
  ipcte2 = ipproc(itemp2(numcla))
 
304
  ipcgev = ipproc(igmeva(numcla))
 
305
  ipcght = ipproc(igmhtf(numcla))
 
306
  ipchgl = ipproc(ih1hlf(numcla))
 
307
 
 
308
  do iel = 1, ncel
 
309
 
 
310
    t2mt1 =  propce(iel,ipcte2)-propce(iel,ipcte1)
 
311
    gmvap = -propce(iel,ipcgev)*t2mt1
 
312
    gmhet = -propce(iel,ipcght)
 
313
 
 
314
    smbrs(iel) = smbrs(iel)                                       &
 
315
         - propce(iel,ipcrom)*volume(iel)*(gmvap+gmhet)
 
316
    if ( rtpa(iel,ivar).gt.epsifl ) then
 
317
      rhovst = propce(iel,ipcrom)*volume(iel)*(gmvap + gmhet)       &
 
318
              / rtpa(iel,ivar)
 
319
    else
 
320
      rhovst = 0.d0
 
321
    endif
 
322
    rovsdt(iel) = rovsdt(iel) + max(zero,rhovst)
 
323
 
 
324
  enddo
 
325
 
 
326
! --> T.S. pour le traceur de la vapeur
 
327
 
 
328
elseif ( ivar .eq. isca(ifvap) ) then
 
329
 
 
330
  if (iwarni(ivar).ge.1) then
 
331
    write(nfecra,1000) chaine(1:8)
 
332
  endif
 
333
 
 
334
  do icla = 1, nclafu
 
335
 
 
336
    ipcte2 = ipproc(itemp2(icla))
 
337
    ipcgev = ipproc(igmeva(icla))
 
338
 
 
339
    do iel = 1, ncel
 
340
 
 
341
      t2mt1 = propce(iel,ipcte2)-propce(iel,ipcte1)
 
342
      if ( rtpa(iel,isca(iyfol(icla))) .gt. epsifl ) then
 
343
        gmvap = -propce(iel,ipcgev)*t2mt1*rtp(iel,isca(iyfol(icla)))    &
 
344
                / rtpa(iel,isca(iyfol(icla)))
 
345
      else
 
346
        gmvap = -propce(iel,ipcgev)*t2mt1
 
347
      endif
 
348
 
 
349
      smbrs(iel) = smbrs(iel)                                     &
 
350
                 + gmvap*propce(iel,ipcrom)*volume(iel)
 
351
    enddo
 
352
 
 
353
  enddo
 
354
 
 
355
! --> T.S. pour le traceur du C ex r�action heterogene
 
356
 
 
357
elseif ( ivar .eq. isca(if7m) ) then
 
358
 
 
359
  if (iwarni(ivar).ge.1) then
 
360
    write(nfecra,1000) chaine(1:8)
 
361
  endif
 
362
 
 
363
  do icla = 1, nclafu
 
364
 
 
365
    ipcght = ipproc(igmhtf(icla))
 
366
 
 
367
    do iel = 1, ncel
 
368
      if (rtpa(iel,isca(iyfol(icla))) .gt. epsifl) then
 
369
        smbrs(iel) = smbrs(iel)                                        &
 
370
             -propce(iel,ipcrom)*propce(iel,ipcght)*volume(iel)        &
 
371
                                *rtp(iel,isca(iyfol(icla)))            &
 
372
                                /rtpa(iel,isca(iyfol(icla)))
 
373
      else
 
374
        smbrs(iel) = smbrs(iel)                                        &
 
375
                    -propce(iel,ipcrom)*propce(iel,ipcght)*volume(iel)
 
376
      endif
 
377
 
 
378
    enddo
 
379
 
 
380
  enddo
 
381
 
 
382
endif
 
383
 
 
384
! --> Terme source pour la variance du traceur 4 (Air)
 
385
 
 
386
if ( ivar.eq.isca(ifvp2m) ) then
 
387
 
 
388
  if (iwarni(ivar).ge.1) then
 
389
    write(nfecra,1000) chaine(1:8)
 
390
  endif
 
391
 
 
392
! ---- Calcul des termes sources explicite et implicite
 
393
!      relatif aux echanges interfaciaux entre phases
 
394
 
 
395
  call cs_fuel_fp2st &
 
396
 !==================
 
397
 ( nvar   , nscal  , ncepdp , ncesmp ,                             &
 
398
   iscal  ,                                                        &
 
399
   itypfb ,                                                        &
 
400
   icepdc , icetsm , itypsm ,                                      &
 
401
   dt     , rtpa   , rtp    , propce , propfa , propfb ,           &
 
402
   smbrs  , rovsdt )
 
403
 
 
404
endif
 
405
 
 
406
 
 
407
! --> Terme source pour CO2
 
408
 
 
409
if ( ieqco2 .ge. 1 ) then
 
410
 
 
411
  if ( ivar.eq.isca(iyco2) ) then
 
412
 
 
413
 
 
414
    if (iwarni(ivar).ge.1) then
 
415
      write(nfecra,1000) chaine(1:8)
 
416
    endif
 
417
 
 
418
! ---- Contribution du TS interfacial aux bilans explicite et implicite
 
419
 
 
420
! Oxydation du CO
 
421
! ===============
 
422
 
 
423
!  Dryer Glassman : XK0P en (moles/m3)**(-0.75) s-1
 
424
!          XK0P = 1.26D10
 
425
!           XK0P = 1.26D7 * (1.1)**(NTCABS)
 
426
!           IF ( XK0P .GT. 1.26D10 ) XK0P=1.26D10
 
427
!           T0P  = 4807.D0
 
428
!  Howard : XK0P en (moles/m3)**(-0.75) s-1
 
429
!             XK0P = 4.11D9
 
430
!             T0P  = 15090.D0
 
431
!  Westbrook & Dryer
 
432
 
 
433
    lnk0p = 23.256d0
 
434
    t0p  = 20096.d0
 
435
!
 
436
!  Hawkin et Smith Purdue University Engeneering Bulletin, i
 
437
!  Research series 108 vol 33, n 3n 1949
 
438
!  Kp = 10**(4.6-14833/T)
 
439
!  Constante d'equilibre en pression partielle (atm           !)
 
440
!  XKOE est le log decimal de la constante pre-exponentielle
 
441
!  TOE  n'est PAS une temerature d'activation  ... il reste un lg(e)
 
442
!  pour repasser en Kc et utiliser des concetrations (moles/m3)
 
443
!  Kc = (1/RT)**variation nb moles * Kp
 
444
!  ici Kc = sqrt(0.082*T)*Kp
 
445
 
 
446
    l10k0e = 4.6d0
 
447
    t0e  = 14833.d0
 
448
! Dissociation du CO2 (Trinh Minh Chinh)
 
449
! ===================
 
450
!          XK0M = 5.D8
 
451
!          T0M  = 4807.D0
 
452
!          XK0M = 0.D0
 
453
!  Westbrook & Dryer
 
454
 
 
455
    lnk0m = 20.03d0
 
456
    t0m  = 20096.d0
 
457
 
 
458
    err1mx = 0.d0
 
459
    err2mx = 0.d0
 
460
 
 
461
! Nombre d'iterations
 
462
    itermx = 500
 
463
! Nombre de points converges
 
464
 
 
465
   nbpauv = 0
 
466
   nbepau = 0
 
467
   nbrich = 0
 
468
   nberic = 0
 
469
   nbpass = 0
 
470
   nbarre = 0
 
471
   nbimax = 0
 
472
! Precision pour la convergence
 
473
   errch = 1.d-8
 
474
 
 
475
   do iel = 1, ncel
 
476
!
 
477
     xxco  = propce(iel,ipproc(iym1(ico  )))/wmole(ico)           &
 
478
            *propce(iel,ipproc(irom1))
 
479
     xxo2  = propce(iel,ipproc(iym1(io2  )))/wmole(io2)           &
 
480
            *propce(iel,ipproc(irom1))
 
481
     xxco2 = propce(iel,ipproc(iym1(ico2 )))/wmole(ico2)          &
 
482
            *propce(iel,ipproc(irom1))
 
483
     xxh2o = propce(iel,ipproc(iym1(ih2o )))/wmole(ih2o)          &
 
484
            *propce(iel,ipproc(irom1))
 
485
!
 
486
     xxco  = max(xxco ,zero)
 
487
     xxo2  = max(xxo2 ,zero)
 
488
     xxco2 = max(xxco2,zero)
 
489
     xxh2o = max(xxh2o,zero)
 
490
     sqh2o = sqrt(xxh2o)
 
491
!
 
492
     xkp = exp(lnk0p-t0p/propce(iel,ipproc(itemp1)))
 
493
     xkm = exp(lnk0m-t0m/propce(iel,ipproc(itemp1)))
 
494
!
 
495
     xkpequ = 10.d0**(l10k0e-t0e/propce(iel,ipproc(itemp1)))
 
496
     xkcequ = xkpequ                                              &
 
497
             /sqrt(8.32d0*propce(iel,ipproc(itemp1))/1.015d5)
 
498
 
 
499
!        initialisation par l'�tat transport�
 
500
 
 
501
     anmr  = xxco2
 
502
     xcom  = xxco + xxco2
 
503
     xo2m  = xxo2 + 0.5d0*xxco2
 
504
!
 
505
     if ( propce(iel,ipproc(itemp1)) .gt. 1200.d0 ) then
 
506
 
 
507
!           Recherche de l'�tat d'�quilibre
 
508
!           Recerche it�rative sans controle de convergence
 
509
!            (pour conserver la parallelisation sur les mailles)
 
510
!           sur le nombre de moles de reaction s�parant
 
511
!           l'etat avant r�action (tel que calcul� par Cpcym)
 
512
!           de l'�tat d'�quilibre
 
513
!          ANMR doit etre borne entre 0 et Min(XCOM,2.*XO2M)
 
514
!          on recherche la solution par dichotomie
 
515
 
 
516
       anmr0 = 0.d0
 
517
       anmr1 = min(xcom,2.d0*xo2m)
 
518
       iterch = 0
 
519
       fn2    = 1.d0
 
520
       fn0  = -0.5d0                           * anmr0**3         &
 
521
            + (     xcom    + xo2m - xkcequ**2) * anmr0**2        &
 
522
            - (.5d0*xcom    +2.d0*xo2m)*xcom   * anmr0            &
 
523
            +       xcom**2 * xo2m
 
524
       fn1  = -0.5d0                           * anmr1**3         &
 
525
            + (     xcom    + xo2m - xkcequ**2) * anmr1**2        &
 
526
            - (.5d0*xcom    +2.d0*xo2m)*xcom   * anmr1            &
 
527
            +       xcom**2 * xo2m
 
528
 
 
529
       if ( xo2m.gt.1.d-6) then
 
530
         do while ( iterch.lt.itermx .and. fn2.gt.errch )
 
531
           anmr2 = 0.5d0*(anmr0+anmr1)
 
532
           fn2  = -0.5d0                            * anmr2**3    &
 
533
                + (     xcom    + xo2m - xkcequ**2) * anmr2**2    &
 
534
                - (.5d0*xcom    +2.d0*xo2m)*xcom    * anmr2       &
 
535
                +       xcom**2 * xo2m
 
536
           if(fn0*fn2 .gt. 0.d0) then
 
537
             anmr0 = anmr2
 
538
             fn0 = fn2
 
539
           elseif(fn1*fn2 .gt. 0.d0) then
 
540
             anmr1 = anmr2
 
541
             fn1 = fn2
 
542
           elseif(fn0*fn1 .gt. 0.d0) then
 
543
             iterch = itermx
 
544
             anmr2 = min(xcom,2.d0*xo2m)
 
545
             nbarre = nbarre + 1
 
546
           endif
 
547
           iterch = iterch + 1
 
548
         enddo
 
549
!
 
550
         if ( iterch .ge. itermx) then
 
551
           nberic = nberic + 1
 
552
         else
 
553
           nbimax = max(nbimax,iterch)
 
554
         endif
 
555
         err1mx = max(err1mx,fn2)
 
556
 
 
557
         xco2eq = anmr2
 
558
         xcoeq  = xcom - anmr2
 
559
         xo2eq  = xo2m - 0.5d0 * anmr2
 
560
       else
 
561
         xo2eq  = 0.d0
 
562
         xcoeq  = xxco
 
563
         xco2eq = 0.d0
 
564
       endif
 
565
 
 
566
     else
 
567
 
 
568
       xco2eq = min(xcom,2.d0*xo2m)
 
569
       xo2eq  = xo2m - 0.5d0*xco2eq
 
570
       xcoeq  = xcom - xco2eq
 
571
 
 
572
     endif
 
573
!
 
574
     if ( xco2eq.gt.xxco2 ) then
 
575
!           oxydation
 
576
       xden = xkp*sqh2o*(xxo2)**0.25d0
 
577
     else
 
578
!           dissociation
 
579
       xden = xkm
 
580
     endif
 
581
     if ( xden .ne. 0.d0 ) then
 
582
 
 
583
       tauchi = 1.d0/xden
 
584
       tautur = rtpa(iel,ik)/rtpa(iel,iep)
 
585
 
 
586
       x2 = 0.d0
 
587
       do icla = 1, nclafu
 
588
         x2 = x2 + rtpa(iel,isca(iyfol(icla)))
 
589
       enddo
 
590
 
 
591
!    On transporte CO2
 
592
 
 
593
       smbrs(iel)  = smbrs(iel)                                   &
 
594
                    +wmole(ico2)/propce(iel,ipproc(irom1))        &
 
595
         * (xco2eq-xxco2)/(tauchi+tautur)                         &
 
596
         * (1.d0-x2)                                              &
 
597
         * volume(iel) * propce(iel,ipcrom)
 
598
!
 
599
       w1 = volume(iel)*propce(iel,ipcrom)/(tauchi+tautur)
 
600
       rovsdt(iel) = rovsdt(iel) +   max(w1,zero)
 
601
 
 
602
     else
 
603
       rovsdt(iel) = rovsdt(iel) + 0.d0
 
604
       smbrs(iel)  = smbrs(iel)  + 0.d0
 
605
     endif
 
606
 
 
607
   enddo
 
608
!
 
609
   if(irangp.ge.0) then
 
610
     call parcpt(nberic)
 
611
     call parmax(err1mx)
 
612
     call parcpt(nbpass)
 
613
     call parcpt(nbarre)
 
614
     call parcpt(nbarre)
 
615
     call parcmx(nbimax)
 
616
   endif
 
617
 
 
618
   write(nfecra,*) ' Max Error = ', err1mx
 
619
   write(nfecra,*) ' no Points   ', nberic, nbarre, nbpass
 
620
   write(nfecra,*) ' Iter max number ', nbimax
 
621
 
 
622
 
 
623
  endif
 
624
 
 
625
endif
 
626
 
 
627
 
 
628
! --> Terme source pour HCN et NO : uniquement a partir de la 2eme
 
629
!                                   iter
 
630
 
 
631
if ( ieqnox .eq. 1 .and. ntcabs .gt. 1) then
 
632
 
 
633
  if ( ivar.eq.isca(iyhcn) .or. ivar.eq.isca(iyno) ) then
 
634
 
 
635
    iexp1  = ipproc(ighcn1)
 
636
    iexp2  = ipproc(ighcn2)
 
637
    iexp3  = ipproc(ignoth)
 
638
 
 
639
! QPR= %N lib�r� pendant l'evaporation/taux de matieres volatiles
 
640
!          moyen
 
641
 
 
642
    qpr = 1.3d0
 
643
 
 
644
! YMOY = % vapeur en sorties
 
645
 
 
646
    ymoy = 0.7d0
 
647
 
 
648
! Azote dans le fuel
 
649
 
 
650
    fn = 0.015
 
651
 
 
652
! Masse molaire
 
653
 
 
654
    wmhcn = wmole(ihcn)
 
655
    wmno  = 0.030d0
 
656
    wmo2  = wmole(io2)
 
657
 
 
658
    if ( ivar.eq.isca(iyhcn) ) then
 
659
 
 
660
!        Terme source HCN
 
661
 
 
662
      if (iwarni(ivar).ge.1) then
 
663
        write(nfecra,1000) chaine(1:8)
 
664
      endif
 
665
 
 
666
      auxmin = 1.d+20
 
667
      auxmax =-1.d+20
 
668
 
 
669
      do iel=1,ncel
 
670
!
 
671
        xxo2 = propce(iel,ipproc(iym1(io2)))                        &
 
672
              *propce(iel,ipproc(immel))/wmo2
 
673
 
 
674
        aux = volume(iel)*propce(iel,ipcrom)                       &
 
675
             *( propce(iel,iexp2)                                  &
 
676
               +propce(iel,iexp1)*rtpa(iel,isca(iyno))             &
 
677
                                 *propce(iel,ipproc(immel))/wmno )
 
678
 
 
679
        smbrs(iel)  = smbrs(iel)  - aux*rtpa(iel,ivar)
 
680
        rovsdt(iel) = rovsdt(iel) + aux
 
681
 
 
682
        gmvap = 0.d0
 
683
        gmhet = 0.d0
 
684
        do icla=1,nclafu
 
685
 
 
686
          ipcgev = ipproc(igmeva(icla))
 
687
          ipcght = ipproc(igmhtf(icla))
 
688
          ipcte2 = ipproc(itemp2(icla))
 
689
          ipcte1 = ipproc(itemp1)
 
690
 
 
691
          gmvap = gmvap                                           &
 
692
                 + propce(iel,ipcrom)*propce(iel,ipcgev)          &
 
693
                  *(propce(iel,ipcte2)-propce(iel,ipcte1))
 
694
 
 
695
          gmhet = gmhet                                           &
 
696
                 +propce(iel,ipcrom)*propce(iel,ipcght)
 
697
 
 
698
        enddo
 
699
        if ( xxo2 .gt. 0.03d0 ) then
 
700
          aux = -volume(iel)*fn*wmhcn/(wmole(in2)/2.d0)           &
 
701
                *( qpr*gmvap+(1.d0-qpr*ymoy)/(1.d0-ymoy)*gmhet )
 
702
        else
 
703
          aux = -volume(iel)*fn*wmhcn/(wmole(in2)/2.d0)           &
 
704
                            *(qpr*gmvap)
 
705
        endif
 
706
        smbrs(iel)  = smbrs(iel) + aux
 
707
 
 
708
      enddo
 
709
 
 
710
    endif
 
711
 
 
712
    if ( ivar.eq.isca(iyno) ) then
 
713
 
 
714
!        Terme source NO
 
715
 
 
716
      if (iwarni(ivar).ge.1) then
 
717
        write(nfecra,1000) chaine(1:8)
 
718
      endif
 
719
 
 
720
      do iel=1,ncel
 
721
 
 
722
        aux1 = volume(iel)*propce(iel,ipcrom)                     &
 
723
              *propce(iel,iexp1)*rtpa(iel,isca(iyhcn))            &
 
724
              *propce(iel,ipproc(immel))/wmhcn
 
725
        aux2 = volume(iel)*propce(iel,ipcrom)                     &
 
726
              *propce(iel,iexp2)*rtpa(iel,isca(iyhcn))            &
 
727
              *wmno/wmhcn
 
728
        aux3 = volume(iel)*propce(iel,ipcrom)**1.5d0              &
 
729
              *propce(iel,iexp3)                                  &
 
730
              *propce(iel,ipproc(iym1(in2)))
 
731
 
 
732
        smbrs(iel)  = smbrs(iel) - aux1*rtpa(iel,ivar)            &
 
733
                               + aux2 + aux3
 
734
        rovsdt(iel) = rovsdt(iel) + aux1
 
735
      enddo
 
736
 
 
737
    endif
 
738
 
 
739
  endif
 
740
 
 
741
endif
 
742
 
 
743
!--------
 
744
! Formats
 
745
!--------
 
746
 
 
747
 1000 format(' TERMES SOURCES PHYSIQUE PARTICULIERE POUR LA VARIABLE '  &
 
748
       ,a8,/)
 
749
 2000 format(' ATTENTION : LE TAUX DE VAPEUR EN SORTIE',                &
 
750
       ' EST NULLE',/,                                                  &
 
751
       ' ON PREND QPR = 0.5')
 
752
 
 
753
!----
 
754
! End
 
755
!----
 
756
 
 
757
return
 
758
 
 
759
end subroutine