* fichier : flamarrh.dgibi
**********************************************************************
**** ****
**** OPERATEUR FLAM ****
**** Combustion Hydrogene-Oxygene, cinetique de type Arrhenius ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/LTMF FEVRIER 1999 ****
**********************************************************************
'OPTION' 'DIME' 2 ;
'OPTION' 'ELEM' QUA4 ;
'OPTION' 'TRAC' 'X' ;
'OPTION' 'ECHO' 0 ;
*
*** GRAPH
*
* GRAPH = VRAI ;
GRAPH = FAUX ;
***********************************************************************
***************** CALCUL DU TAUX DE REACTION *************************
***********************************************************************
'DEBPROC' TAUX ;
*
**** La meme chose que FLAM ARRHENIUS
*
ARGUMENT TS*CHPOINT A*FLOTTANT b*FLOTTANT con*FLOTTANT Ta*FLOTTANT
H0H2*FLOTTANT H0O2*FLOTTANT H0H2O*FLOTTANT DT*FLOTTANT
RYH2*CHPOINT RYO2*CHPOINT RYH2O*CHPOINT T*CHPOINT ;
MH2 = 2.016D-3 ;
MO2 = 31.998D-3 ;
MH2O = 18.015D-3 ;
MAIL = 'EXTRAIRE' TS 'MAILLAGE' ;
NPOI = 'NBEL' MAIL ;
LIRYH2 = 'PROG' NPOI * 0.0D0 ;
LIRYO2 = 'PROG' NPOI * 0.0D0 ;
LIRYH2O = 'PROG' NPOI * 0.0D0 ;
LIDE = 'PROG' NPOI * 0.0D0 ;
'REPETER' BL1 NPOI ;
P1 = 'POIN' &BL1 MAIL ;
RYH2I = 'EXTRAIRE' RYH2 'H2 ' P1 ;
RYO2I = 'EXTRAIRE' RYO2 'O2 ' P1 ;
RYH2OI = 'EXTRAIRE'RYH2O 'H2O ' P1 ;
TS0 = 'EXTRAIRE' TS 'SCAL' P1 ;
T0 = 'EXTRAIRE' T 'SCAL' P1 ;
XH2 = RYH2I '/' MH2 ;
XO2 = RYO2I '/' MO2 ;
'SI' (T0 > TS0) ;
'SI' (XH2 > (2.0D0 * XO2)) ;
* H2 majoritaire
LAMBDA = EXP(-1.0D0*TA/T0);
LAMBDA = LAMBDA * (T0 ** (-1.0D0*B));
LAMBDA = LAMBDA * (RYH2I ** CON);
LAMBDA = LAMBDA * MO2 * A;
NRYO2I = RYO2I / (1.0D0 + (LAMBDA * DT));
DRYO2 = NRYO2I '-' RYO2I;
DMOL = DRYO2 / MO2;
DRYH2 = 2.0D0 * DMOL * MH2;
DRYH2O = -2.0D0 * DMOL * MH2O;
NRYH2I = RYH2I '+' DRYH2 ;
NRYH2OI = RYH2OI '+' DRYH2O ;
'MESSAGE' ('CHAINE' 'POINT ' &BL1) ;
'MESSAGE' 'O2 minoritaire' ;
'MESSAGE' ('CHAINE' 'DTc = ' (1.0 '/' LAMBDA) ' DT = ' DT) ;
'SINON' ;
* O2 majoritaire
LAMBDA = EXP(-1.0D0*TA/T0);
LAMBDA = LAMBDA * (T0 ** (-1.0D0*B));
LAMBDA = LAMBDA * (RYH2I ** (CON-1.0D0)) * RYO2I;
LAMBDA = 2.0D0 * LAMBDA * MH2 * A;
NRYH2I = RYH2I / (1.0D0 + (LAMBDA * DT)) ;
DRYH2 = NRYH2I '-' RYH2I;
DMOL = DRYH2 / MH2 / 2.0D0;
DRYO2 = DMOL * MO2 ;
DRYH2O = -2.0D0 * DMOL * MH2O;
NRYO2I = RYO2I '+' DRYO2 ;
NRYH2OI = RYH2OI '+' DRYH2O ;
'MESSAGE' ('CHAINE' 'POINT ' &BL1) ;
'MESSAGE' 'H2 minoritaire ' ;
'MESSAGE' ('CHAINE' 'DTc = ' (1.0 '/' LAMBDA) ' DT = ' DT) ;
'FINSI' ;
DEI = (DRYH2 * H0H2) + (DRYO2 * H0O2) + (DRYH2O * H0H2O) ;
DEI = -1.0D0 * DEI ;
'SINON' ;
NRYO2I = RYO2I ;
NRYH2I = RYH2I ;
NRYH2OI = RYH2OI ;
DEI = 0.0D0 ;
'FINSI' ;
'REMPLACER' LIRYH2 &BL1 NRYH2I ;
'REMPLACER' LIRYO2 &BL1 NRYO2I ;
'REMPLACER' LIRYH2O &BL1 NRYH2OI ;
'REMPLACER' LIDE &BL1 DEI ;
'FIN' BL1 ;
NRYH2 = 'MANUEL' 'CHPO' MAIL 1 'H2 ' LIRYH2 ;
NRYO2 = 'MANUEL' 'CHPO' MAIL 1 'O2 ' LIRYO2 ;
NRYH2O = 'MANUEL' 'CHPO' MAIL 1 'H2O ' LIRYH2O ;
DE = 'MANUEL' 'CHPO' MAIL 1 'SCAL' LIDE ;
'FINPROC' NRYH2 NRYO2 NRYH2O DE ;
***********************
**** LA TABLE PGAZ ****
***********************
PGAZ = 'TABLE' ;
*
**** Ordre des polynoms
*
PGAZ . 'NORD' = 4 ;
*
**** Especes qui sont dans les equations d'Euler
*
PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'N2 ';
*
PGAZ . 'H2 ' = 'TABLE' ;
PGAZ . 'H2O ' = 'TABLE' ;
PGAZ . 'N2 ' = 'TABLE' ;
PGAZ . 'O2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'H2 ' . 'R' = 4130.0 ;
PGAZ . 'H2O ' . 'R' = 461.4 ;
PGAZ . 'N2 ' . 'R' = 296.8 ;
PGAZ . 'O2 ' . 'R' = 259.8 ;
*
**** Regressions polynomials pour le cvs
*
PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
-1.82753232E-08 2.44485692E-12 ;
PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
8.78233606E-09 -3.05514485E-13 ;
PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
'SI' FAUX ;
*
**** Calcul des enthalpies de formation a 0K
*
* A 298.15K
*
* HOH2 = 0 ;
* HOO2 = 0 ;
* HON2 = 0 ;
* HOH2O = -1.3428D7 J/kg ;
*
RH2 = PGAZ . 'H2 ' . 'R' ;
RO2 = PGAZ . 'O2 ' . 'R' ;
RN2 = PGAZ . 'N2 ' . 'R' ;
RH2O = PGAZ . 'H2O ' . 'R' ;
A0H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 1 ;
A1H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 2 ;
A2H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 3 ;
A3H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 4 ;
A4H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 5 ;
A0H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 1 ;
A1H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 2 ;
A2H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 3 ;
A3H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 4 ;
A4H2O= 'EXTRAIRE' (PGAZ . 'H2O ' . 'A') 5 ;
A0N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 1 ;
A1N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 2 ;
A2N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 3 ;
A3N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 4 ;
A4N2 = 'EXTRAIRE' (PGAZ . 'N2 ' . 'A') 5 ;
A0O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 1 ;
A1O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 2 ;
A2O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 3 ;
A3O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 4 ;
A4O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 5 ;
* On calcule \int_0^298.15 cv(y) dy pour chaque especes
T = 298.15 ;
T2 = T * T ;
T3 = T2 * T ;
T4 = T3 * T ;
T5 = T4 * T ;
ETHERH2 = ((A0H2 * T) '+' ((A1H2 '/' 2.0) * T2) '+'
((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
((A4H2 '/' 5.0) * T5)) ;
ETHERO2 = ((A0O2 * T) '+' ((A1O2 '/' 2.0) * T2) '+'
((A2O2 '/' 3.0) * T3) '+' ((A3O2 '/' 4.0) * T4) '+'
((A4O2 '/' 5.0) * T5)) ;
ETHERH2O= ((A0H2O * T) '+' ((A1H2O '/' 2.0) * T2) '+'
((A2H2O '/' 3.0) * T3) '+' ((A3H2O '/' 4.0) * T4) '+'
((A4H2O '/' 5.0) * T5)) ;
ETHERN2 = ((A0N2 * T) '+' ((A1N2 '/' 2.0) * T2) '+'
((A2N2 '/' 3.0) * T3) '+' ((A3N2 '/' 4.0) * T4) '+'
((A4N2 '/' 5.0) * T5)) ;
*
ENTHH2 = -1.0D0 * (ETHERH2 '+' (RH2 '*' 298.15)) ;
ENTHO2 = -1.0D0 * (ETHERO2 '+' (RO2 '*' 298.15)) ;
ENTHN2 = -1.0D0 * (ETHERN2 '+' (RN2 '*' 298.15)) ;
ENTHH2O = -1.0D0 * (1.3428D7 '+' ETHERH2O '+' (RH2O '*' 298.15)) ;
'FINSI' ;
*
**** "Enthalpies" (ou energies) de formations a OK (J/Kg)
*
PGAZ . 'H2 ' . 'H0K' = -4.194D6 ;
PGAZ . 'H2O ' . 'H0K' = -1.395D7 ;
PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
***************************
***** DOMAINE SPATIAL ****
***************************
*
**** Trois points
*
A1 = 0.0D0 0.0D0 ;
DOM1 = 'MANUEL' 'POI1' A1 ;
A2 = 1.0D0 0.0D0 ;
DOM2 = 'MANUEL' 'POI1' A2 ;
A3 = 2.0D0 0.0D0 ;
DOM3 = 'MANUEL' 'POI1' A3 ;
DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
*
**** FLAM, 'ARRHENIU'
*
* 2 H2 + O2 -> 2 H2O
*
* w = A H(T - Ts) RH2^c RO2 * T^(-b) exp(-Ta/T)
*
* ou
*
* w = taux de reaction MOLAIRE pour la reaction
* chimie, qui est suppose' etre elementaire
*
* H(T - Ts) est la fontion d'Heavyside
*
* N.B. L'energie thermique de chaque epsece est definie par
*
* e = \int_{0}^T cv(y) dy
*
Ta = 8310.0;
A = 1.1725D14;
b = 0.91;
con = 2.0;
DT = 0.1;
TS = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' 1000 'NATURE' 'DISCRET' ;
********************************************
************* Donnes initiales *************
********************************************
RH2 = PGAZ . 'H2 ' . 'R' ;
RO2 = PGAZ . 'O2 ' . 'R' ;
RN2 = PGAZ . 'N2 ' . 'R' ;
RH2O = PGAZ . 'H2O ' . 'R' ;
H0H2 = PGAZ . 'H2 ' . 'H0K' ;
H0O2 = PGAZ . 'O2 ' . 'H0K' ;
H0H2O = PGAZ . 'H2O ' . 'H0K' ;
* H0N2 = PGAZ . 'N2 ' . 'H0K' ;
MH2 = 2.016D-3 ;
MO2 = 31.998D-3 ;
MH2O = 18.015D-3 ;
MN2 = 28.013D-3 ;
XH2 = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 0.3 0.25 0.3) ;
XO2 = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 0.147 0.147 0.147) ;
XH2O = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 1.0D-6 1.0D-6 1.0D-6) ;
XN2 = ('MANUEL' 'CHPO' DOMTOT 1 'SCAL' 1.0 'NATURE' 'DISCRET' ) '-'
(XH2 '+' XO2 '+' XH2O) ;
* Mass des especes dans une mole
MH2SMOL = XH2 * MH2 ;
MO2SMOL = XO2 * MO2 ;
MH2OSMOL = XH2O * MH2O ;
MN2SMOL = XN2 * MN2 ;
MTOT = MH2SMOL '+' MO2SMOL '+' MH2OSMOL '+' MN2SMOL ;
YH2 = MH2SMOL '/' MTOT ;
YO2 = MO2SMOL '/' MTOT ;
YH2O = MH2OSMOL '/' MTOT ;
YN2 = MN2SMOL '/' MTOT ;
RTOT = (YH2 * RH2) '+' (YO2 * RO2) '+' (YH2O * RH2O) '+' (YN2 * RN2) ;
* Pression et temperature
P = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 2.83D6 2.0D6 2.5D6) ;
T = 'MANUEL' 'CHPO' DOMTOT 1 'SCAL' ('PROG' 1527D0 2500 1650 ) ;
RO = P '/' (RTOT '*' T) ;
RYH2 = 'NOMC' 'H2 ' (RO * YH2) ;
RYO2 = 'NOMC' 'O2 ' (RO * YO2) ;
RYH2O = 'NOMC' 'H2O ' (RO * YH2O) ;
* Dans le cas de detonation dans un tube ou C1 = zone von Neumann et
* C3 = reste du tube, on a
* u ~ 1000 m/s
* gamma ~ 1.25
* donc
* ason ~ 875 m/s -> (u + cson) ~ 1875 m/s
*
* DT ~ DX '/' (u + cson)
*
* DX ~ 0.1 m -> DT ~ 5.00000E-05
*
DT = 0.5D-5 ;
NRYH2 NRYO2 NRYH2O DE = 'FLAM' 'ARRHENIU' TS A b con Ta H0H2 H0O2
H0H2O DT RYH2 RYO2 RYH2O T ;
*
**** Conservativite'
*
ERRRO = (('EXCO' 'H2 ' NRYH2) '+' ('EXCO' 'O2 'NRYO2) '+'
('EXCO' 'H2O ' NRYH2O)) '-'
(('EXCO' 'H2 ' RYH2) '+' ('EXCO' 'O2 ' RYO2) '+'
('EXCO' 'H2O ' RYH2O));
ERRRO = ERRRO '/' ('MAXIMUM' RO) ;
'SI' (('MAXIMUM' ERRRO 'ABS') > 1.0D-10) ;
'ERREUR' 5 ;
'FINSI' ;
NRYH21 NRYO21 NRYH2O1 DE1 = TAUX TS A b con Ta H0H2 H0O2
H0H2O DT RYH2 RYO2 RYH2O T ;
ERRH2 = ('MAXIMUM' (NRYH2 '-' NRYH21) 'ABS')
'/' ('MAXIMUM' NRYH2 ) ;
ERRO2 = 'MAXIMUM' (NRYO2 '-' NRYO21) 'ABS'
'/' ('MAXIMUM' NRYO2 ) ;
ERRH2O = 'MAXIMUM' (NRYH2O '-' NRYH2O1) 'ABS'
'/' ('MAXIMUM' NRYH2O ) ;
ERRDE = 'MAXIMUM' (DE '-' DE1) 'ABS' '/'
('MAXIMUM' DE) ;
'SI' (('MAXIMUM' ('PROG' ERRH2 ERRO2 ERRH2O ERRDE)) > 1.0D-10) ;
'ERREUR' 5;
'FINSI' ;
*
**** Option HEAVYSIDE: les especes minoritaires brule completement
*
NRYH2H NRYO2H NRYH2OH DEH = 'FLAM' 'HEAVYSID' TS H0H2 H0O2
H0H2O RYH2 RYO2 RYH2O T ;
*
**** Conservativite'
*
ERRRO = (('EXCO' 'H2 ' NRYH2H) '+' ('EXCO' 'O2 'NRYO2H) '+'
('EXCO' 'H2O ' NRYH2OH)) '-'
(('EXCO' 'H2 ' RYH2) '+' ('EXCO' 'O2 ' RYO2) '+'
('EXCO' 'H2O ' RYH2O));
ERRRO = ERRRO '/' ('MAXIMUM' RO) ;
'SI' (('MAXIMUM' ERRRO 'ABS') > 1.0D-10) ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer