1
!-------------------------------------------------------------------------------
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
5
! Copyright (C) 1998-2011 EDF S.A.
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
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
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.
21
!-------------------------------------------------------------------------------
23
subroutine cs_fuel_noxst &
24
!=======================
28
pdfm1 , pdfm2 , doxyd , dfuel , hrec , &
29
f3m , f4m , f5m , f6m , f7m , f8m , f9m , &
30
fs3no , fs4no , yfs4no , enthox , &
33
!===============================================================================
38
! K1 exp(-E1/RT) (conversion HCN en N2)
39
! K2 exp(-E2/RT) (conversion HCN en NO)
40
! K3 exp(-E3/RT) (NO thermique)
44
!__________________.____._____.________________________________________________.
45
! name !type!mode ! role !
46
!__________________!____!_____!________________________________________________!
47
! ncelet ! i ! <-- ! number of extended (real + ghost) cells !
48
! ncel ! i ! <-- ! number of cells !
49
! indpdf ! te ! <-- ! passage par les pdf !
50
!__________________!____!_____!________________________________________________!
52
! TYPE : E (ENTIER), R (REEL), A (ALPHAMNUMERIQUE), T (TABLEAU)
53
! L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
54
! MODE : <-- donnee, --> resultat, <-> Donnee modifiee
55
! --- tableau de travail
56
!===============================================================================
58
!==============================================================================
60
!==============================================================================
79
!===============================================================================
88
double precision pdfm1(ncel) , pdfm2(ncel) , dfuel(ncel)
89
double precision f3m(ncel) , f4m(ncel) , f5m(ncel)
90
double precision f6m(ncel) , f7m(ncel) , f8m(ncel) ,f9m(ncel)
91
double precision doxyd(ncel) , hrec(ncel)
92
double precision fs3no(ncel) , fs4no(ncel) , yfs4no(ncel,ngazg)
93
double precision enthox(ncel)
95
double precision rtp(ncelet,*), propce(ncelet,*)
99
integer iel , icla , i
100
integer ipdf1 , ipdf2 , ipdf3
101
integer iexp1 , iexp2 , iexp3
104
parameter (npart = 200 )
106
double precision ee1,ee2,ee3,kk1,kk2,kk3
107
double precision wmo2,tg,xo2,bb,xmx2
108
double precision bb1 , bb2 , bb3 , bb4
109
double precision dirac , tfuel , tmpgaz , lrf ,lro
110
double precision qqq , rrr , sss , ttt , uuu
111
double precision gt1,gt2,gt3,gt10,gt20,gt30
112
double precision ts4,ts4num,ts4den
114
double precision val(npart+1),tt(npart+1) , gs(npart+1) , yyo2(npart+1)
116
integer icha,ige,numcha,mode
117
double precision xxf , hhf , hfuel,tfs4ad,hfs4ad
118
double precision xhf1,xhf2 , aa , den
119
double precision xesp(ngazem),coefe(ngazem),f1mc(ncharm),f2mc(ncharm)
120
double precision yo2moy,yo2ox,yo2cb,yo24den,yo24num,yo2s4,yo2
121
double precision toxyd,hoxyd , som , react , deltamol
123
integer i300 ,i000 ,imini,i2500,i2600,i2700,i2800,i3000,i3500
124
integer i4000,i5000,imaxi,inok
125
integer nbpt , nbclip1 , nbclip2
126
integer nbclip30,nbclip31
127
double precision ts4min,ts4max
128
double precision somm , sommax , sommin , ts4admin,ts4admax
129
double precision yo2min,yo2max,yo2min1,yo2max1
130
double precision yo2oxmin,yo2oxmax
131
double precision toxmin,toxmax
134
!===============================================================================
135
! 1. CALCULS PRELIMINAIRES
136
!===============================================================================
138
! --- Masses molaires
144
iexp1 = ipproc(ighcn1)
145
iexp2 = ipproc(ighcn2)
146
iexp3 = ipproc(ignoth)
148
! Parametres des lois d'Arrhenius
157
! Pour les termes, indicateur de calcul par PDF ou non
158
! = 1 --> passage par pdf
159
! = 0 --> on ne passe pas par les pdf
168
propce(iel,iexp1) = zero
169
propce(iel,iexp2) = zero
170
propce(iel,iexp3) = zero
173
!===============================================================================
174
! 2. CALCUL SANS LES PDF
175
!===============================================================================
179
tg = propce(iel,ipproc(itemp1))
180
yo2 = propce(iel,ipproc(iym1(io2)))
181
xo2 = yo2*propce(iel,ipproc(immel))/wmo2
183
! Reaction HCN + NO + 1/4 O2 ---> N2 + 1/2 H2O + CO
185
propce(iel,iexp1) = kk1*exp(-ee1/tg)
187
! Reaction HCN + 5/4 O2 --> NO + 1/2 H2O + CO
189
if ( xo2 .gt. 0.018d0 ) then
191
else if ( xo2 .lt. 0.0025d0 ) then
194
bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
196
propce(iel,iexp2) = kk2*exp(-ee2/tg)*(xo2**bb)
198
! No thermique : Zeldovich
200
propce(iel,iexp3) = kk3*exp(-ee3/tg)*(xo2**0.5d0)
204
!===============================================================================
205
! 3. CALCUL AVEC LES PDF
206
!===============================================================================
208
if ( ipdf1 .eq. 1 .or. ipdf2 .eq. 1 .or. ipdf3 .eq. 1 ) then
246
if ( indpdf(iel) .eq. 1 ) then
248
! Calcul de Yo2 dans l'oxydant
251
bb1 = max(0.d0 ,pdfm1(iel))
252
bb2 = min(fs3no(iel),pdfm2(iel))
253
if ( bb2 .gt. bb1 ) then
258
qqq = bb2**2 - bb1**2
260
gt1 = lro*qqq/(2.d0*fs3no(iel))
263
if ( propce(iel, ipproc(iym1(io2))) .gt. 0.d0 ) then
264
yo2ox = propce(iel, ipproc(iym1(io2)))/(-gt1+gt2+gt3)
266
yo2oxmin = min(yo2oxmin,yo2ox)
267
yo2oxmax = max(yo2oxmax,yo2ox)
269
yo2moy = propce(iel, ipproc(iym1(io2)))
271
dirac = dfuel(iel)*yo2cb + doxyd(iel)*yo2ox
273
bb1 = max(0.D0 ,pdfm1(iel))
274
bb2 = min(fs4no(iel),pdfm2(iel))
275
bb3 = max(fs4no(iel),pdfm1(iel))
276
bb4 = min(fs3no(iel),pdfm2(iel))
277
if ( bb2 .gt. bb1 ) then
282
if ( bb4 .gt. bb3 ) then
288
qqq = bb2**2 - bb1**2
290
sss = bb4**2 - bb3**2
292
uuu = fs4no(iel)-fs3no(iel)
294
gt1 = lro*qqq/(2.d0*fs4no(iel))
297
gt10= lrf*sss/(2.d0*uuu)
298
gt20= lrf*ttt*fs3no(iel)/uuu
300
yo24num = yo2moy - dirac + yo2ox*(gt1 -gt2)
301
yo24den = gt1+gt10-gt20
303
yo2s4 = yo24num/yo24den
310
! Affichage et clipping
312
yo2min = min(yo2s4,yo2min)
313
yo2max = max(yo2s4,yo2max)
314
if ( yo2s4 .lt. 0.d0 ) then
315
nbclip30 = nbclip30+1
318
if ( yo2s4 .gt. yo2ox ) then
319
nbclip31 = nbclip31+1
322
yo2min1 = min(yo2s4,yo2min1)
323
yo2max1 = max(yo2s4,yo2max1)
326
! Calcul de Tfioul moyen
327
! ----------------------
333
xmx2 = xmx2 + rtp(iel,isca(iyfol(icla)))
335
if ( xmx2 .gt. 0.d0 ) then
337
ipctem=ipproc(itemp2(icla))
338
tfuel = tfuel +rtp(iel,isca(iyfol(icla)))*propce(iel,ipctem)
344
tfuel = propce(iel,ipproc(itemp1))
347
! On recupere la valeur de Toxyd a partir de hoxyd
354
coefe(io2) =( af3(io2 )*f3m(iel)+af4(io2 )*f4m(iel)+af5(io2 )*f5m(iel) &
355
+af6(io2 )*f6m(iel)+af7(io2 )*f7m(iel)+af8(io2 )*f8m(iel) &
356
+af9(io2 )*f9m(iel) ) * wmole(io2)
358
coefe(in2) =( af3(in2 )*f3m(iel)+af4(in2 )*f4m(iel)+af5(in2 )*f5m(iel) &
359
+af6(in2 )*f6m(iel)+af7(in2 )*f7m(iel)+af8(in2 )*f8m(iel) &
360
+af9(in2 )*f9m(iel) ) * wmole(in2)
362
coefe(ico2)=( af3(ico2)*f3m(iel)+af4(ico2)*f4m(iel)+af5(ico2)*f5m(iel) &
363
+af6(ico2)*f6m(iel)+af7(ico2)*f7m(iel)+af8(ico2)*f8m(iel) &
364
+af9(ico2)*f9m(iel) ) * wmole(ico2)
366
coefe(ih2o)=( af3(ih2o)*f3m(iel)+af4(ih2o)*f4m(iel)+af5(ih2o)*f5m(iel) &
367
+af6(ih2o)*f6m(iel)+af7(ih2o)*f7m(iel)+af8(ih2o)*f8m(iel) &
368
+af9(ih2o)*f9m(iel) ) * wmole(ih2o)
370
coefe(ico)=( af3(ico)*f3m(iel)+af4(ico)*f4m(iel)+af5(ico)*f5m(iel) &
371
+af6(ico)*f6m(iel)+af7(ico)*f7m(iel)+af8(ico)*f8m(iel) &
372
+af9(ico)*f9m(iel) ) * wmole(ico)
374
coefe(ihy)=( af3(ihy)*f3m(iel)+af4(ihy)*f4m(iel)+af5(ihy)*f5m(iel) &
375
+af6(ihy)*f6m(iel)+af7(ihy)*f7m(iel)+af8(ihy)*f8m(iel) &
376
+af9(ihy)*f9m(iel) ) * wmole(ihy)
378
! Nbre de mole qui a reagis avec CO pour faire du CO2
380
deltamol = (coefe(io2)-yo2ox)/wmole(io2)
381
react = min(2.d0*deltamol,coefe(ico)/wmole(ico))
382
coefe(io2 )= coefe(io2 )-react/2.d0*wmole(io2)
383
coefe(ico )= coefe(ico )-react *wmole(ico)
384
coefe(ico2)= coefe(ico2)+react *wmole(ico2)
387
som = coefe(io2)+coefe(in2)+coefe(ico2) &
388
+coefe(ih2o)+coefe(ico)+coefe(ihy)
389
coefe(io2 ) = coefe(io2 ) /som
390
coefe(in2 ) = coefe(in2 ) /som
391
coefe(ico2) = coefe(ico2) /som
392
coefe(ih2o) = coefe(ih2o) /som
393
coefe(ico ) = coefe(ico ) /som
394
coefe(ihy ) = coefe(ihy ) /som
402
call cs_fuel_htconvers1 ( mode , hoxyd , coefe , toxyd )
403
!======================
405
toxmin = min(toxmin,toxyd)
406
toxmax = max(toxmax,toxyd)
408
if ( toxyd .gt. propce(iel,ipproc(itemp1)) ) then
409
toxyd = propce(iel,ipproc(itemp1))
412
! On initialise par les temperatures Toxy et Tfuel aux extremites
414
dirac = dfuel(iel)*tfuel + doxyd(iel)*toxyd
416
! On recupere la valeur de la temperature moyenne
418
tmpgaz = propce(iel,ipproc(itemp1))
420
bb1 = max(0.D0 ,pdfm1(iel))
421
bb2 = min(fs4no(iel),pdfm2(iel))
422
bb3 = max(fs4no(iel),pdfm1(iel))
423
bb4 = min(1.d0 ,pdfm2(iel))
425
if ( bb2 .gt. bb1 ) then
430
if ( bb4 .gt. bb3 ) then
436
qqq = bb2**2 - bb1**2
438
sss = bb4**2 - bb3**2
440
uuu = 1.d0 - fs4no(iel)
442
gt1 = lro*qqq/(2.d0*fs4no(iel))
445
gt10= lrf*sss/(2.d0*uuu)
449
ts4num = tmpgaz - dirac + toxyd*(gt1 -gt2 ) &
450
- tfuel*(gt10+gt20-gt30)
451
ts4den = gt1-gt10+gt30
455
! Calcul de l'enthalpie des vapeur a Tfuel
459
hhf = hhf + rtp(iel,isca(ih2(icla)))
462
if ( xmx2 .gt. epsifl ) then
466
hfs4ad = fs4no(iel)*hfuel+(1.d0-fs4no(iel))*hoxyd
472
coefe(ige) = yfs4no(iel,ige)
476
call cs_fuel_htconvers1 ( mode , hfs4ad , coefe , tfs4ad )
477
!======================
480
! Calcul pour affichage
482
ts4min = min(ts4min,ts4)
483
ts4max = max(ts4max,ts4)
485
ts4admin= min(ts4admin,tfs4ad)
486
ts4admax= max(ts4admax,tfs4ad)
490
somm = somm + yfs4no(iel,ige)
492
sommin=min(sommin,somm)
493
sommax=max(sommax,somm)
495
if ( (ts4 .gt. min(toxyd,tmpgaz,tfuel)) .and. &
496
(ts4 .lt. 2400.d0) ) inok = inok + 1
497
if ( ts4 .lt. min(toxyd,tmpgaz,tfuel) ) then
498
if ( ts4 .ge. 300.d0 ) then
500
else if ( ts4 .gt. 0 ) then
507
if ( ts4 .gt. 2400.d0 ) then
508
if ( ts4 .lt. 2500 ) then
510
else if ( ts4 .lt. 2600 ) then
512
else if ( ts4 .lt. 2700 ) then
514
else if ( ts4 .lt. 2800 ) then
516
else if ( ts4 .lt. 3000 ) then
518
else if ( ts4 .lt. 3500 ) then
520
else if ( ts4 .lt. 4000 ) then
522
else if ( ts4 .lt. 5000 ) then
529
! Fin calcul pour affichage
532
! Clipping de Ts4 : a min(toxyd,tmpgaz,tfuel) en min
536
if ( ts4 .lt. min(toxyd,tmpgaz,tfuel) ) then
537
nbclip1 = nbclip1 + 1
538
ts4 = min(toxyd,tmpgaz,tfuel)
541
if ( ts4 .gt. tfs4ad ) then
542
nbclip2 = nbclip2 + 1
546
! Concentration oxygene
548
xo2 = propce(iel,ipproc(iym1(io2 ))) &
549
*propce(iel,ipproc(immel))/wmole(io2)
554
gs(i) = pdfm1(iel)+dble(i-1)/dble(npart)*(pdfm2(iel)-pdfm1(iel))
556
if( gs(i) .lt. fs4no(iel) ) then
557
tt(i) = (ts4-toxyd)/fs4no(iel)* gs(i) + toxyd
559
tt(i) = (tfuel-ts4)/(1.d0-fs4no(iel))*gs(i) &
560
+ tfuel - (tfuel-ts4)/(1.d0-fs4no(iel))
563
if ( gs(i) .lt. fs4no(iel) ) then
564
yyo2(i) = (yo2s4-yo2ox)/fs4no(iel) * gs(i) + yo2ox
565
else if ( gs(i) .lt. fs3no(iel) ) then
566
aa = yo2s4/(fs4no(iel)-fs3no(iel))
567
yyo2(i) = aa * ( gs(i) -fs3no(iel) )
575
dgs = ( pdfm2(iel) - pdfm1(iel) ) / dble(npart)
577
! Calcul de K1*EXP(-E1/T)
579
if ( ipdf1 .eq. 1 ) then
581
propce(iel,iexp1)= kk1*exp(-ee1/toxyd)*doxyd(iel) &
582
+kk1*exp(-ee1/tfuel)*dfuel(iel)
585
val(i) = kk1*exp(-ee1/tt(i))*hrec(iel)
589
propce(iel,iexp1) = propce(iel,iexp1) &
590
+0.5d0*dgs*(val(i)+val(i+1))
595
! Calcul de K2*EXP(-E2/T)
597
if ( ipdf2 .eq. 1 ) then
599
if ( xo2 .gt. 0.d0 ) then
601
if(xo2.gt.0.018d0) then
603
else if(xo2 .lt. 0.0025d0) then
606
bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
609
propce(iel,iexp2) = kk2*exp(-ee2/toxyd)*doxyd(iel)*(xo2**bb) &
610
+kk2*exp(-ee2/tfuel)*dfuel(iel)*(xo2**bb)
613
val(i) = kk2*exp(-ee2/tt(i))*hrec(iel)
617
propce(iel,iexp2) = propce(iel,iexp2) &
618
+0.5d0*dgs*(val(i)+val(i+1))*(xo2**bb)
621
propce(iel,iexp2) = 0.d0
626
! Calcul de K3*EXP(-E3/T)
628
if ( ipdf3 .eq. 1 ) then
630
if ( xo2 .gt. 0.d0 ) then
632
propce(iel,iexp3) = kk3*exp(-ee3/toxyd) &
633
*doxyd(iel)*(yo2ox**0.5d0) &
634
+kk3*exp(-ee3/tfuel) &
635
*dfuel(iel)*(yo2cb**0.5d0)
638
if ( gs(i) .le. fs3no(iel) ) then
639
val(i) = kk3*exp(-ee3/tt(i))*hrec(iel)*(yyo2(i)**0.5d0)
646
propce(iel,iexp3) = propce(iel,iexp3) &
647
+0.5d0*dgs*(val(i)+val(i+1))
651
propce(iel,iexp3)= 0.d0
660
if ( irangp .ge. 0 ) then
679
call parmin(ts4admin)
680
call parmax(ts4admax)
684
call parcpt(nbclip30)
685
call parcpt(nbclip31)
693
call parmin(yo2oxmin)
694
call parmax(yo2oxmax)
696
!===============================================================================
698
write(nfecra,*) ' Min max de TSox ',toxmin,toxmax
699
write(nfecra,*) ' Min max de TS4 ',ts4min,ts4max
700
write(nfecra,*) ' nombre de pts a Tmini < ts4 < 2400 ',inok
701
write(nfecra,*) ' nombre de pts a ts4 < 0 ',imini
702
write(nfecra,*) ' nombre de pts a 0 < ts4 < 300 ',i000
703
write(nfecra,*) ' nombre de pts a 300 < ts4 < Tmini ',i300
704
write(nfecra,*) ' nombre de pts a 2400 < ts4 < 2500 ',i2500
705
write(nfecra,*) ' nombre de pts a 2500 < ts4 < 2600 ',i2600
706
write(nfecra,*) ' nombre de pts a 2600 < ts4 < 2700 ',i2700
707
write(nfecra,*) ' nombre de pts a 2700 < ts4 < 2800 ',i2800
708
write(nfecra,*) ' nombre de pts a 2800 < ts4 < 3000 ',i3000
709
write(nfecra,*) ' nombre de pts a 3000 < ts4 < 3500 ',i3500
710
write(nfecra,*) ' nombre de pts a 3500 < ts4 < 4000 ',i4000
711
write(nfecra,*) ' nombre de pts a 4000 < ts4 < 5000 ',i5000
712
write(nfecra,*) ' nombre de pts a 5000 < ts4 ',imaxi
713
write(nfecra,*) ' Min max de TS4ad ',ts4admin,ts4admax
714
write(nfecra,*) ' Min max concentration en fs4 ',sommin,sommax
715
write(nfecra,*) ' Clipping en min en ',nbclip1,' points sur ',nbpt,' points'
716
write(nfecra,*) ' Clipping en max en ',nbclip2,' points sur ',nbpt,' points'
719
write(nfecra,*) ' Min max de Yo2ox en 0 ',yo2oxmin,yo2oxmax
720
write(nfecra,*) ' Min max de Yo2 en fs4 avant clipping ',yo2min,yo2max
721
write(nfecra,*) ' Clipping en min sur Yo2 en fs4 ',nbclip30
722
write(nfecra,*) ' Clipping en max sur Yo2 en fs4 ',nbclip31
723
write(nfecra,*) ' Min max de Yo2 en fs4 apres clipping ',yo2min1,yo2max1
724
!===============================================================================
added
removed
src/comb/cs_fuel_noxst.f90
