* fichier : comb.dgibi
*************************************************************************
*
* NOM : comb.dgibi
*
* DESCRIPTION : We compute the AICC (Adiabatic Isochoric Complete
* Combustion), the AIBCC, (Adiabatic IsoBaric Complete
* Combustion), the CJ (Chapman - Jouguet) and VN (von
* Neumann) states of a ZND detonation for several mixtures
* at NIST-normal conditions involving H2 and air.
* We check that some physical properties are satisfied.
* Comparison of the results obtained using the DETO
* operator and the ones given by the CREBCOM combustion
* model (operators PRIM and FLAM).
*
* LANGAGE : GIBIANE-CAST3M
* AUTEUR : Alberto BECCANTINI (CEA/DEN/DM2S/SFME/LTMF)
* beccantini@cea.fr
* DATE : 28/09/2006
*************************************************************************
*
'OPTION' 'ECHO' 1 'DIME' 2
'TRAC' 'X' 'ELEM' 'QUA4' ;
****************************************************************
***** 1D mesh *****
****************************************************************
*
*
* A4 ---- A3
* | |
* | COM |
* | |
* A1 ---- A2
*
A1 = 0.0 0.0 ;
A2 = 1.0 0.0 ;
A1A2 = A1 'DROIT' 1 A2 ;
DOM1 = A1A2 'TRANSLATION' 1 (0.0 1.0) ;
*
**** The table 'DOMAINE'
*
$DOM1 = 'MODELISER' DOM1 'EULER' ;
TDOM1 = 'DOMA' $DOM1 'VF' ;
QDOM1 = TDOM1 . 'QUAF' ;
*
***** Initial conditions and gas properties.
*
*
**** 0 from a numerical point of view
*
zero = 1.0d-9 ;
*
PINI = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.013D5
'NATU' 'DISCRET' ;
TINI = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 298.15
'NATU' 'DISCRET' ;
WINI = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 2 'UX' 0.0 'UY' 0.0
'NATU' 'DISCRET' ;
* We suppose that the combustion is complete everywhere
CSIMAX = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.0 ;
*
***** Gas properties
*
PGAZ = 'TABLE' ;
* Polynomial degree of specific heats
PGAZ . 'NORD' = 4 ;
* Species explicitly treated in the Euler Equations
PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
* Species non explicitly treated
PGAZ . 'ESPNEULE' = 'N2 ';
* Single gases properties
PGAZ . 'H2 ' = 'TABLE' ;
PGAZ . 'H2O ' = 'TABLE' ;
PGAZ . 'N2 ' = 'TABLE' ;
PGAZ . 'O2 ' = 'TABLE' ;
* R (J/Kg/K)
mH2 = 2. '*' 1.00797E-3 ;
mo2 = 2. '*' 15.9994E-3 ;
mH2O = mh2 '+' (0.5 '*' mo2) ;
mN2 = 2 '*' 14.0067E-3 ;
RGAS = 8.31441 ;
PGAZ . 'H2 ' . 'R' = RGAS '/' mh2 ;
PGAZ . 'H2O ' . 'R' = RGAS '/' mh2o ;
PGAZ . 'N2 ' . 'R' = RGAS '/' mn2 ;
PGAZ . 'O2 ' . 'R' = RGAS '/' mo2 ;
* Polynomials regressions coefficients
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;
* Formation enthalpies at 0K (J/Kg)
*
* h_i(0K) = h_i(T0) '-' \int_0^{T0} cp_i(x) dx
* = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
*
PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
PGAZ . 'H2O ' . 'H0K' = -1.395D7 ;
PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
* Passive scalars names
PGAZ . 'SCALPASS' = 'MOTS' 'H2IN' 'H2FI' 'K0 ';
*
**** Initial conditions on molar fractions of hydrogen.
*
LISTXH2 = PROG 0.08 0.134 0.188 0.242 0.29581 0.350 ;
LISXH2F = PROG ;
LISXO2F = PROG ;
LISXWAF = PROG ;
LISXN2F = PROG ;
LISYH2F = PROG ;
LISYO2F = PROG ;
LISYWAF = PROG ;
LISYN2F = PROG ;
* Inside
LRINSI = PROG ;
LTINSI = PROG ;
LPINSI = PROG ;
LCINSI = PROG ;
LRGINSI = PROG ;
LGINSI = PROG ;
LGAINSI = PROG ;
* VN conditions
LRVN = PROG ;
LTVN = PROG ;
LPVN = PROG ;
LCVN = PROG ;
LRGVN = PROG ;
LGVN = PROG ;
LGAVN = PROG ;
* AICC
LRAICC = PROG ;
LTAICC = PROG ;
LPAICC = PROG ;
LCAICC = PROG ;
LRGAICC = PROG ;
LGAICC = PROG ;
LGAAICC = PROG ;
* AIBCC
LRAIBC = PROG ;
LTAIBC = PROG ;
LPAIBC = PROG ;
LCAIBC = PROG ;
LRGAIBC = PROG ;
LGAIBC = PROG ;
LGAAIBC = PROG ;
* CJ conditions
LRCJ = PROG ;
LTCJ = PROG ;
LPCJ = PROG ;
LCCJ = PROG ;
LRGCJ = PROG ;
LGCJ = PROG ;
LGACJ = PROG ;
* Released energy per mass unit
LQ = PROG ;
*
REPE BLCALC (DIME LISTXH2) ;
************************************************************************
************************************************************************
***** LOOP to determine the initial conditions, AICC, AIBCC, CJ *****
***** anv VN states *****
************************************************************************
************************************************************************
*
**** Molar fractions (before combustion)
*
cell = (EXTR LISTXH2 &BLCALC) ;
cell1 = (1. - cell) / (1. + (0.79/0.21)) ;
XH21 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' cell
'NATU' 'DISCRET' ;
XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2O' zero
'NATURE' 'DISCRET' ;
XO21 = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1
'O2' cell1 'NATU' 'DISCRET' ;
XN = XH21 'ET' XH2O1 'ET' XO21 ;
CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 1.0D0) ;
XN21 = CHUN '-' ('PSCAL' XN CHUN ('MOTS' 'H2' 'O2' 'H2O')
('MOTS' 'SCAL' 'SCAL' 'SCAL') ) ;
XN21 = 'NOMC' 'N2' XN21 'NATURE' 'DISCRET' ;
XN = XN 'ET' XN21 ;
*
**** We compute the mass fractions
*
MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' Mh2
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'O2' Mo2
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2O' Mh2o
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'N2' Mn2
'NATURE' 'DISCRET') ;
MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
YH21 = (XH21 '*' MH2) '/' MTOT ;
YO21 = (XO21 '*' MO2) '/' MTOT ;
YN21 = (XN21 '*' MN2) '/' MTOT ;
YH2O1 = (XH2O1 '*' MH2O) '/' MTOT ;
*
**** Burnt gas
* Here X refers to the mole producted per 1 mole of unburned gas
* N.B. The number of moles changes after the combustion!
*
* H2MAG = 1 in the region where H2 is the majority species.
* Conversely H2MAG = 0
H2MAG = (0.5 '*' ('EXCO' 'H2' XN)) 'MASQUE' 'SUPERIEUR'
('EXCO' 'O2' XN) ;
O2MAG = CHUN '-' H2MAG ;
XO22a = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'O2' zero ;
XH22a = XH21 '-' ('NOMC' 'H2' (2. '*' (XO21 '-' XO22a))) ;
XO22a = XO22a ;
XH22a = XH22a ;
XH22b = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'H2' zero ;
XO22b = XO21 '-' ('NOMC' 'O2' (0.5 '*' (xh21 '-' xh22b))) ;
XH22 = (XH22a '*' H2MAG) '+' (XH22b '*' O2MAG) ;
XO22 = (XO22a '*' H2MAG) '+' (XO22b '*' O2MAG) ;
XH2O2 = XH2O1 '+' ('NOMC' 'H2O' (2 '*' (XO21 '-' XO22))) ;
XN22 = XN21 ;
*
**** mtot = weight of a mole of unburned gas
* mtot2 = weight of a mole of unburned gas after the combustion
* Of course mtot2 = mtot
*
MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
mtotcel = 'MAXIMUM' MTOT ;
'SI' ( ('MAXIMUM' (mtot2 '-' mtot) 'ABS') '>' (zero '*' mtotcel)) ;
'MESSAGE' 'Problem in the computation of the burnt gas (1).' ;
'ERREUR' 21 ;
'FINSI' ;
*
**** xtot2 = mole of burnt gas per mole of unburned gas
*
xtot2 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
*
**** X = real molar fractions
*
xh22 = xh22 '/' xtot2 ;
xo22 = xo22 '/' xtot2 ;
xh2o2 = xh2o2 '/' xtot2 ;
xn22 = xn22 '/' xtot2 ;
LISXH2F = LISXH2F 'ET' (PROG (MAXI XH22)) ;
LISXO2F = LISXO2F 'ET' (PROG (MAXI XO22)) ;
LISXWAF = LISXWAF 'ET' (PROG (MAXI XH2O2)) ;
LISXN2F = LISXN2F 'ET' (PROG (MAXI XN22)) ;
xtot3 = 'PSCAL' (xh22 '+' xo22 '+' xh2O2 '+' xn22)
CHUN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL');
'SI' ( ('MAXIMUM' (xtot3 '-' 1.0D0) 'ABS' ) '>'
zero) ;
'MESSAGE' 'Problem in the computation of the burnt gas (2).' ;
'ERREUR' 21 ;
'FINSI' ;
*
**** mtot2 = weight of 1 mole of burned gas
*
MTOT2 = 'PSCAL' MMOL (XH22 '+' XO22 '+' XH2O2 '+' XN22)
('MOTS' 'H2' 'O2' 'H2O' 'N2')
('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
*
**** Mass fractions of burned gas after the complete combustion
*
YH22 = (XH22 '*' MH2) '/' MTOT2 ;
YO22 = (XO22 '*' MO2) '/' MTOT2 ;
YN22 = (XN22 '*' MN2) '/' MTOT2 ;
YH2O2 = (XH2O2 '*' MH2O) '/' MTOT2 ;
LISYH2F = LISYH2F 'ET' (PROG (MAXI YH22)) ;
LISYO2F = LISYO2F 'ET' (PROG (MAXI YO22)) ;
LISYWAF = LISYWAF 'ET' (PROG (MAXI YH2O2)) ;
LISYN2F = LISYN2F 'ET' (PROG (MAXI YN22)) ;
*
**** The conservative variables
*
* Gas constant
RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ .'H2' . 'R')) '+'
(('NOMC' 'SCAL' yo21) '*' (PGAZ .'O2' . 'R')) '+'
(('NOMC' 'SCAL' yh2o1) '*' (PGAZ .'H2O' . 'R')) '+'
(('NOMC' 'SCAL' yn21) '*' (PGAZ .'N2' . 'R'))
;
* Density
RNINI = PINI '/' (RINI '*' TINI) ;
RN = 'COPIER' RNINI ;
* Internal energy
eini = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
Tfois = Tini ;
'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
&BL1)) ;
eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
Tfois = Tfois * tini ;
'FIN' BL1 ;
RETNINI = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
RETN = RETNINI ;
GN = RN '*' WINI ;
YN = YH21 '+' YO21 '+' YH2O1 ;
RYN = RN '*' YN ;
*
**** K0/DX = chemical source term
*
* i.e. dcsi/dt = K0 '/' DX . {criterium function}
* where
* {criterium function} = 1.0 or 0.0
*
* d(rho etot)/dt = -rho (dcsi/dt) (h_f '-' h_i)
*
* h_f reference enthalpy (0K) of the burned gas
* h_i reference enthalpy (0K) of the unburned gas
*
**** Mass fraction of H2 before and after the combustion
*
* CSIMAX = 1.0 -> complete combustion
*
YININ = yh21 ;
YFINN = YININ '+' ((yh22 '-' yh21) '*' csimax) ;
K0N = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 500 ;
*
**** Passive scalars
*
RSN = RN '*' (('NOMC' 'H2IN' YININ 'NATU' 'DISCRET') 'ET'
('NOMC' 'H2FI' YFINN 'NATU' 'DISCRET') 'ET'
('NOMC' 'K0' K0N 'NATU' 'DISCRET')) ;
* Test of prim
VN P T Y S GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
PCEL = 'MAXIMUM' PINI ;
TCEL = 'MAXIMUM' TINI ;
'SI' (('MAXIMUM' ((P '-' pini) '/' PCEL) 'ABS')> 1.0D-3) ;
'MESSAGE' 'Problem in the computation of the energy (1).'
'ERREUR' 21 ;
'FINSI' ;
'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
'MESSAGE' 'Problem in the computation of the energy (2).'
'ERREUR' 21 ;
'FINSI' ;
PINSI1 = MAXI P ;
TINSI1 = MAXI T ;
RINSI1 = MAXI RN ;
GAMINSI = MAXI GAMN ;
GAMAINSI = 1 + (MAXI (P / RETN)) ;
RGAS = MAXI (P / (RN * T)) ;
LRINSI = LRINSI 'ET' (PROG RINSI1) ;
LTINSI = LTINSI 'ET' (PROG TINSI1) ;
LPINSI = LPINSI 'ET' (PROG PINSI1) ;
LRGINSI = LRGINSI 'ET' (PROG RGAS) ;
LGINSI = LGINSI 'ET' (PROG GAMINSI) ;
LGAINSI = LGAINSI 'ET' (PROG GAMAINSI) ;
LCINSI = LCINSI 'ET' (PROG ((GAMINSI * PINSI1 / RINSI1) ** 0.5)) ;
*
* VN state
*
* VN state via DETO operator (not really complete combustion...)
*
CHPO1 = 'MANUEL''CHPO' ('DOMA' $DOM1 'CENTRE') 6
'H2' (MAXI XH21)
'O2' (MAXI XO21)
'N2' (1. - ((MAXI XH21) + (MAXI XO21)))
'H2O' 0.0
'P' (MAXI PINI) 'T' (MAXI TINI) ;
CHPCJ CHPVN CHPAICC = 'DETO' CHPO1 ;
PVN = (EXCO CHPVN 'PZND') ;
RVN = (EXCO CHPVN 'RZND') ;
TVN = (EXCO CHPVN 'TZND') ;
DVN = (EXCO CHPCJ 'VCJ') ;
* VVN via RH conditions
VVN = DVN '*' (RNINI / RVN) ;
VVN = DVN - VVN ;
VVN = (NOMC VVN 'UX') +
('MANUEL''CHPO' ('DOMA' $DOM1 'CENTRE') 1 'UY' 0.0) ;
* Internal energy
EVN = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
Tfois = TVN ;
'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
&BL1)) ;
eVN = eVN '+' (acel '*' (Tfois '/' &BL1)) ;
Tfois = Tfois * tVN ;
'FIN' BL1 ;
RETNVN = RVN '*' (eVN '+' (0.5 '*' ('PSCAL' VVN VVN
('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
V1VN P1VN T1VN YNVN SNVN GAMVN = 'PRIM' 'PERFTEMP' PGAZ RVN (RVN * VVN)
RETNVN (RVN * YN) RSN ;
ERRO = MAXI ((TVN - T1VN) / T1VN) ABS ;
'SI' (ERRO > 1.0d-3) ;
MESS 'Problem in VN' ;
ERREUR 21 ;
FINS ;
*
* We check that the RH on the energy is satisfied
*
VVN = DVN '*' (RNINI / RVN) ;
HVN = (EVN + (PVN / RVN)) + (0.5 * VVN * VVN) ;
HINI = (EINI + (PINI / RNINI)) + (0.5 * DVN * DVN) ;
ERRO = MAXI ((HVN - HINI) / HVN) ABS ;
'SI' (ERRO > 1.0d-2) ;
MESS 'Problem in VN (2)' ;
ERREUR 21 ;
FINS ;
PVN1 = MAXI PVN ;
TVN1 = MAXI TVN ;
RVN1 = MAXI RVN ;
GAMVN = MAXI GAMVN ;
GAMAVN = 1 + (MAXI (PVN / (RVN * EVN))) ;
RGAS = MAXI (PVN / (RVN * TVN)) ;
LRVN = LRVN 'ET' (PROG RVN1) ;
LTVN = LTVN 'ET' (PROG TVN1) ;
LPVN = LPVN 'ET' (PROG PVN1) ;
LRGVN = LRGVN 'ET' (PROG RGAS) ;
LGVN = LGVN 'ET' (PROG GAMVN) ;
LGAVN = LGAVN 'ET' (PROG GAMAVN) ;
LCVN = LCVN 'ET' (PROG ((GAMVN * PVN1 / RVN1) ** 0.5)) ;
*
***** Combustion via the operator FLAM
*
DELTAXN = ('DOMA' $DOM1 'VOLUME') '**' (1. '/' 2.) ;
YN = RYN '/' RN ;
*
**** We want the complete combustion in DOM1.
* We can use the operator 'FLAM', with epsilon=0 and
* DELTAT >> DELTATC
* where DELTATC is the characteristic time step linked to
* the source term.
*
SN = RSN '/' RN ;
K0N = 'EXCO' 'K0' SN ;
YININ = 'EXCO' 'H2IN' SN 'H2' ;
YFINN = 'EXCO' 'H2FI' SN 'H2' ;
DELTATC = 0.25 '*' ('MINIMUM' (DELTAXN '/' K0N)) ;
LMOT1 = 'MOTS' 'H2' 'O2' 'H2O' ;
LCOEF = 'PROG' 1.0 0.5 -1.0 ;
DELTARE DELTARY = 'FLAM' 'CREBCOM2' $DOM1 PGAZ LMOT1 LCOEF
RN YN YININ YFINN K0N DELTAXN
0.0 (1D4 '*' DELTATC) 0.3 ;
Q = ((MAXI (DELTARE / RN))) ;
Q1 = (('NOMC' 'SCAL' (yh21 '-' yh22)) * (PGAZ .'H2 ' . 'H0K')) +
(('NOMC' 'SCAL' (yo21 '-' yo22)) * (PGAZ .'O2 ' . 'H0K')) +
(('NOMC' 'SCAL' (yn21 '-' yn22)) * (PGAZ .'N2 ' . 'H0K')) +
(('NOMC' 'SCAL' (yh2o1 '-' yh2o2)) * (PGAZ .'H2O ' . 'H0K')) ;
ERRO = MAXI ((Q - Q1) / Q1) ABS ;
'SI' (ERRO > 1.0d-2) ;
MESS 'Problem in Q' ;
ERREUR 21 ;
FINS ;
LQ = LQ 'ET' (PROG Q) ;
RYN = RYN '+' DELTARY ;
YN = RYN '/' RN ;
RETN = RETN '+' DELTARE ;
*
*
*
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
PAICC1 = MAXI PN ;
TAICC1 = MAXI TN ;
RAICC1 = MAXI RN ;
GAMAICC = MAXI GAMMAN ;
GAMAAICC = 1 + (MAXI (PN / RETN)) ;
RGAS = MAXI (PN / (RN * TN)) ;
*
* AICC, CJ via DETO operator (not really complete combustion...)
*
CHPO1 = 'MANUEL''CHPO' ('DOMA' $DOM1 'CENTRE') 6
'H2' (MAXI XH21)
'O2' (MAXI XO21)
'N2' (1. - ((MAXI XH21) + (MAXI XO21)))
'H2O' 0.0
'P' (MAXI PINI) 'T' (MAXI TINI) ;
CHPCJ CHPZND CHPAICC = 'DETO' CHPO1 ;
PAICC2 = MAXI (EXCO CHPAICC 'PAIC') ;
RAICC2 = MAXI (EXCO CHPAICC 'RAIC') ;
TAICC2 = MAXI (EXCO CHPAICC 'TAIC') ;
*
LERRO = ABS (PROG ((PAICC1 - PAICC2) / PAICC1)
((RAICC1 - RAICC2) / RAICC1)
((TAICC1 - TAICC2) / TAICC1)) ;
ERRO = MAXI LERRO ;
'SI' (ERRO > 1.0d-2) ;
MESS 'Problem in AICC' ;
ERREUR 21 ;
FINS ;
LRAICC = LRAICC 'ET' (PROG RAICC1) ;
LTAICC = LTAICC 'ET' (PROG TAICC1) ;
LPAICC = LPAICC 'ET' (PROG PAICC1) ;
LRGAICC = LRGAICC 'ET' (PROG RGAS) ;
LGAICC = LGAICC 'ET' (PROG GAMAICC) ;
LGAAICC = LGAAICC 'ET' (PROG GAMAAICC) ;
LCAICC = LCAICC 'ET' (PROG ((GAMAICC * PAICC1 / RAICC1) ** 0.5)) ;
*
* CJ
*
PCJ = (EXCO CHPCJ 'PCJ') ;
RCJ = (EXCO CHPCJ 'RCJ') ;
TCJ = (EXCO CHPCJ 'TCJ') ;
DCJ = (EXCO CHPCJ 'VCJ') ;
Tfois = TCJ ;
eCJ = 'MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' 0.0 ;
'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
acel = (('NOMC' 'SCAL' yh22) '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yo22) '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yh2o2) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
&BL1)) '+'
(('NOMC' 'SCAL' yn22) '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
&BL1)) ;
eCJ = eCJ '+' (acel '*' (Tfois '/' &BL1)) ;
Tfois = Tfois * TCJ ;
'FIN' BL1 ;
RETNCJ = RCJ '*' eCJ ;
RWCJ = RCJ * VN * 0 ;
V1CJ P1CJ T1CJ YNCJ SNCJ GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RCJ RWCJ
RETNCJ (RCJ * YN) (RCJ * SN) ;
CCJ = (GAMMAN * (P1CJ / RCJ)) '**' 0.5 ;
VCJ = DCJ * (RNINI / RCJ) ;
erro = MAXI ((VCJ / CCJ) - 1.0) ABS ;
* We check weather the speed in the shock reference is equal to the speed of
* sound.
'SI' (erro > 2.0d-2) ;
MESS 'Problem in CJ' ;
ERREUR 21 ;
FINS ;
GAMACJ = 1 + (MAXI (PCJ / RETNCJ)) ;
LRCJ = LRCJ 'ET' (PROG (maxi RCJ)) ;
LTCJ = LTCJ 'ET' (PROG (maxi TCJ)) ;
LPCJ = LPCJ 'ET' (PROG (maxi PCJ)) ;
LRGCJ = LRGCJ 'ET' (PROG RGAS) ;
LGCJ = LGCJ 'ET' (PROG (MAXI GAMMAN)) ;
LGACJ = LGACJ 'ET' (PROG GAMACJ) ;
LCCJ = LCCJ 'ET' (PROG (maxi CCJ)) ;
*
* AIBCC
*
* Gas constant
*
YO2 = 'MAXIMUM' ('EXCO' 'O2' YN) ;
YH2 = 'MAXIMUM' ('EXCO' 'H2' YN) ;
YH2O = 'MAXIMUM' ('EXCO' 'H2O' YN) ;
YN2 = (-1 '*' ((YO2 + YH2 + YH2O) '-' 1)) ;
R = ((PGAZ . 'H2O ' . 'R') '*' YH2O) '+'
((PGAZ . 'O2 ' . 'R') '*' YO2) '+'
((PGAZ . 'H2 ' . 'R') '*' YH2) '+'
((PGAZ . 'N2 ' . 'R') '*' YN2) ;
T0 = 500 ;
T1 = 3500 ;
hfin = 'MAXIMUM' ((RETN '/' RN) '+' (PINI '/' RN)) ;
* hfin = EINI '+' q '+' (PINI '/' RN) ;
*
* Regula falsi
*
'REPETER' BLITER 100 ;
T = (T0 '+' T1) '/' 2. ;
eini = 0.0 ;
Tfois = T ;
'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
acel = (YH2O '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
&BL1)) '+'
(YN2 '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
&BL1)) '+'
(YO2 '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
&BL1)) '+'
(YH2 '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
&BL1)) ;
eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
Tfois = Tfois * T ;
'FIN' BL1 ;
hini = eini '+' (R * T) ;
'MESS' ('CHAINE' 'h, T ' hini ' ' T) ;
'SI' (< ('ABS' (hini '-' hfin)) 1.0D-4) ;
'MESSAGE' OK ;
'QUITTER' BLITER ;
'SINON' ;
'SI' (hini > hfin) ;
T1 = T ;
'SINON' ;
T0 = T ;
'FINSI' ;
'FINSI' ;
'FIN' BLITER ;
*
TAIBC1=T;
PAIBC1=(MAXI PINI) ;
RNAIBCC = PINI '/' (R '*' T) ;
RAIBC1=(MAXI RNAIBCC) ;
eini = 0.0 ;
Tfois = T ;
'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
acel = (YH2O '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
&BL1)) '+'
(YN2 '*' ('EXTRAIRE' (PGAZ .'N2 ' . 'A')
&BL1)) '+'
(YO2 '*' ('EXTRAIRE' (PGAZ .'O2 ' . 'A')
&BL1)) '+'
(YH2 '*' ('EXTRAIRE' (PGAZ .'H2 ' . 'A')
&BL1)) ;
eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
Tfois = Tfois * T ;
'FIN' BL1 ;
RETAIBCC = ('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' eini) '*'
RNAIBCC ;
RYNAIBCC = RNAIBCC '*' YN ;
VN2 P2 T2 Y2 S2 GAMN = 'PRIM' 'PERFTEMP' PGAZ RNAIBCC GN RETAIBCC
RYNAIBCC RSN ;
TAIBC2= MAXI T2 ;
PAIBC2= MAXI P2 ;
LERRO = ABS (PROG ((PAIBC1 - PAIBC2) / PAIBC1)
((TAIBC1 - TAIBC2) / TAIBC1)) ;
ERRO = MAXI LERRO ;
'SI' (ERRO > 1.0d-2) ;
MESS 'Problem in AIBCC' ;
ERREUR 21 ;
FINS ;
LRAIBC = LRAIBC 'ET' (PROG RAIBC1) ;
LTAIBC = LTAIBC 'ET' (PROG TAIBC1) ;
LPAIBC = LPAIBC 'ET' (PROG PAIBC1) ;
GAMAIBC = MAXI GAMN ;
GAMAAIBC = 1 + (MAXI (P2 / RETAIBCC)) ;
RGAS = MAXI (P2 / (RNAIBCC * T2)) ;
LRGAIBC = LRGAIBC 'ET' (PROG RGAS) ;
LGAIBC = LGAIBC 'ET' (PROG GAMAIBC) ;
LGAAIBC = LGAAIBC 'ET' (PROG GAMAAIBC) ;
LCAIBC = LCAIBC 'ET' (PROG ((GAMAIBC * PAIBC1 / RAIBC1) ** 0.5)) ;
*
*************************************************************************
*************************************************************************
*********************** END OF THE LOOP *********************************
*************************************************************************
*************************************************************************
*
FIN BLCALC ;
*
**** Some LaTex tables.
*
SI FAUX ;
OPTI ECHO 0 ;
REPE BLCALC (DIME LISTXH2) ;
*
R = EXTR LRINSI &BLCALC ;
P = EXTR LPINSI &BLCALC ;
T = EXTR LTINSI &BLCALC ;
C = EXTR LCINSI &BLCALC ;
GAM = EXTR LGINSI &BLCALC ;
GAMA = EXTR LGAINSI &BLCALC ;
RGAS = EXTR LRGINSI &BLCALC ;
AA = CHAI &blcalc ' & Inside & ' 'FORMAT' '(E9.3)' R ' & '
'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
MESS AA ;
*
R = EXTR LRAICC &BLCALC ;
P = EXTR LPAICC &BLCALC ;
T = EXTR LTAICC &BLCALC ;
C = EXTR LCAICC &BLCALC ;
GAM = EXTR LGAICC &BLCALC ;
GAMA = EXTR LGAAICC &BLCALC ;
RGAS = EXTR LRGAICC &BLCALC ;
AA = CHAI &blcalc ' & AICC & ' 'FORMAT' '(E9.3)' R ' & '
'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
MESS AA ;
*
R = EXTR LRAIBC &BLCALC ;
P = EXTR LPAIBC &BLCALC ;
T = EXTR LTAIBC &BLCALC ;
C = EXTR LCAIBC &BLCALC ;
GAM = EXTR LGAIBC &BLCALC ;
GAMA = EXTR LGAAIBC &BLCALC ;
RGAS = EXTR LRGAIBC &BLCALC ;
AA = CHAI &blcalc ' & AIBCC & ' 'FORMAT' '(E9.3)' R ' & '
'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
MESS AA ;
*
R = EXTR LRCJ &BLCALC ;
P = EXTR LPCJ &BLCALC ;
T = EXTR LTCJ &BLCALC ;
C = EXTR LCCJ &BLCALC ;
GAM = EXTR LGCJ &BLCALC ;
GAMA = EXTR LGACJ &BLCALC ;
RGAS = EXTR LRGCJ &BLCALC ;
AA = CHAI &blcalc ' & CJ & ' 'FORMAT' '(E9.3)' R ' & '
'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
MESS AA ;
*
R = EXTR LRVN &BLCALC ;
P = EXTR LPVN &BLCALC ;
T = EXTR LTVN &BLCALC ;
C = EXTR LCVN &BLCALC ;
GAM = EXTR LGVN &BLCALC ;
GAMA = EXTR LGAVN &BLCALC ;
RGAS = EXTR LRGVN &BLCALC ;
AA = CHAI &blcalc ' & VN & ' 'FORMAT' '(E9.3)' R ' & '
'FORMAT' '(E9.3)' P ' & ' 'FORMAT' '(E9.3)' T ' & '
'FORMAT' '(E9.3)' C ' & ' 'FORMAT' '(E9.3)' RGAS ' & '
'FORMAT' '(E9.3)' GAM ' & ' 'FORMAT' '(E9.3)' GAMA ' \\';
MESS AA ;
FIN BLCALC ;
OPTI ECHO 1 ;
OPTI ECHO 0 ;
REPE BLCALC (DIME LISTXH2) ;
*
XH2 = EXTR LISXH2F &BLCALC ;
XO2 = EXTR LISXO2F &BLCALC ;
XWA = EXTR LISXWAF &BLCALC ;
XN2 = EXTR LISXN2F &BLCALC ;
YH2 = EXTR LISYH2F &BLCALC ;
YO2 = EXTR LISYO2F &BLCALC ;
YWA = EXTR LISYWAF &BLCALC ;
YN2 = EXTR LISYN2F &BLCALC ;
Q = EXTR LQ &BLCALC ;
* AA = CHAI &blcalc ' & ' 'FORMAT' '(E9.3)' XH2 ' & '
* 'FORMAT' '(E9.3)' XO2 ' & ' 'FORMAT' '(E9.3)' XWA ' & '
* 'FORMAT' '(E9.3)' XN2 ' & ' 'FORMAT' '(E9.3)' YH2 ' & '
* 'FORMAT' '(E9.3)' YO2 ' & ' 'FORMAT' '(E9.3)' YWA ' & '
* 'FORMAT' '(E9.3)' YN2 ' \\';
AA = CHAI &blcalc ' & ' 'FORMAT' '(E9.3)' XH2 ' & '
'FORMAT' '(E9.3)' XO2 ' & ' 'FORMAT' '(E9.3)' XWA ' & '
'FORMAT' '(E9.3)' YH2 ' & '
'FORMAT' '(E9.3)' YO2 ' & ' 'FORMAT' '(E9.3)' YWA ' & '
'FORMAT' '(E9.3)' Q ' \\';
MESS AA ;
MESS '\hline' ;
FIN BLCALC ;
OPTI ECHO 1 ;
FINSI ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer