* fichier : tubedeto2d1.dgibi
****************************************************************
***** PROPAGATION D UNE DETONATION DANS UN TUBE *****
***** MODELE CREBCOM *****
***** A. BECCANTINI, SFME/LTMF, 15.01.02 *****
****************************************************************
*
* The file is divided into 5 parts
*
* 1) mesh
* 2) initial conditions and gas properties
* 3) parameters used in the computation
* 4) the main part (where the computation is realised)
* 5) the post-treatment
*
'OPTION' 'ECHO' 0 'DIME' 2
'TRAC' 'X' ;
****************************************************************
***** 1D mesh *****
****************************************************************
*
*
* A6 ---- A5 ------------------- A4 ----------------------- A8
* | | | |
* | AR | ZONE1 | ZONE2 |
* | | | |
* A1 ---- A2 ------------------- A3 ----------------------- A7
* | | | |
* | L1 | L2 | L3 |
* |<---->|<--------------------->|<----------------------->|
*
* AR = activation region
* ZONE1
* ZONE2
*
RAF = 1 ;
NL2SL1 = 49 ;
NL3SL1 = 100 ;
L1 = 0.4 ;
L2 = NL2SL1 '*' L1 ;
L3 = NL3SL1 '*' L1 ;
* If NL2SL1 = 49 and NL3SL1 = 50, then L2 = 19.6, L3 = 20.0 ;
N1 = RAF ;
N2 = NL2SL1 '*' N1 ;
N3 = NL3SL1 '*' N1 ;
* For sake of simplicity, we will only consider one mesh in y-direction
DX = L1 '/' N1 ;
NY = 4 ;
NZ = NY ;
H1 = NY '*' DX ;
A1 = 0.0 0.0 ;
A2 = L1 0.0 ;
A3 = (L1 '+' L2) 0.0 ;
A4 = (L1 '+' L2) H1 ;
A5 = L1 H1 ;
A6 = 0.0 H1 ;
A7 = (L1 '+' L2 '+' L3) 0.0 ;
A8 = (L1 '+' L2 '+' L3) H1 ;
'OPTION' 'ELEM' 'SEG2' ;
A1A2 = A1 'DROIT' N1 A2 ;
A2A3 = A2 'DROIT' N2 A3 ;
A3A4 = A3 'DROIT' NY A4 ;
A4A5 = A4 'DROIT' N2 A5 ;
A5A2 = A5 'DROIT' NY A2 ;
A2A5 = 'INVERSE' A5A2 ;
A5A6 = A5 'DROIT' N1 A6 ;
A6A1 = A6 'DROIT' NY A1 ;
A3A7 = A3 'DROIT' N3 A7 ;
A7A8 = A7 'DROIT' NY A8 ;
A8A4 = A8 'DROIT' N3 A4 ;
A4A3 = 'INVERSE' A3A4 ;
*
**** 2D mesh
*
* DOM1 is the activation region
* DOM2 is zone 1
* DOM3 is zone 2
'OPTION' 'ELEM' 'QUA4' ;
DOM1 = 'TRANSLATION' A2A5 N1 ((-1.0 '*' L1) 0.0) ;
DOM2 = 'TRANSLATION' A2A5 N2 (L2 0.0) ;
DOM3 = 'TRANSLATION' A3A4 N3 (L3 0.0) ;
* DOMTOT = total region
DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
'ELIMINATION' (1.0D-3 '/' RAF) DOMTOT ;
*
**** The tables 'DOMAINE'
*
$DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
$DOM1 = 'MODELISER' DOM1 'EULER' ;
$DOM2 = 'MODELISER' DOM2 'EULER' ;
$DOM3 = 'MODELISER' DOM3 'EULER' ;
TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
TDOM1 = 'DOMA' $DOM1 'VF' ;
TDOM2 = 'DOMA' $DOM2 'VF' ;
TDOM3 = 'DOMA' $DOM3 'VF' ;
QDOMTOT = TDOMTOT . 'QUAF' ;
QDOM1 = TDOM1 . 'QUAF' ;
QDOM2 = TDOM2 . 'QUAF' ;
QDOM3 = TDOM3 . 'QUAF' ;
'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM1 ;
'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM2 ;
'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM3 ;
*
**** MOD1 = model (created just to display the CHAMELEMs)
*
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
*
**** Line for graphics
*
'OPTION' 'ELEM' 'SEG2' ;
P1 = (DX '/' 2) (DX '/' 2);
P2 = ((L1 '+' L2 '+' L3) '-' (DX '/' 2)) (DX '/' 2);
LIGEVO = (P1 'DROIT' (N1 '+' N2 '+' N3 '-' 1) P2) 'COULEUR' 'ROUG' ;
'ELIMINATION' LIGEVO ('DOMA' $DOMTOT 'CENTRE') (0.001 '/' RAF) ;
* 'OPTION' 'SAUVER' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;
****************************************************************
****************************************************************
***** Initial conditions and gas properties. *****
****************************************************************
****************************************************************
*
**** 0 from a numerical point of view
*
zero = 1.0d-9 ;
*
**** (Non-Uniform) initial conditions.(DOM1 et DOM2)
* We compute the conservative variables.
*
*
* Contact discontinuity between DOM2 and DOM3 (same pressure,
* different temperature)
*
PINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 99700 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 99700
'NATU' 'DISCRET') ;
TINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 285. 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 400
'NATU' 'DISCRET') ;
WINI = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 2 'UX' 0.0 'UY' 0.0;
* We suppose that the combustion is complete everywhere, i.e. the
* progress variable is 1 after the detonation transition.
CSIMAX = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 1.0 ;
*
**** From initial conditions on molar fractions, we duduce
* the data for the virtual mixture of burned/unburned
* gases
*
**** Molar fractions (before combustion)
*
XH21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
('DOMA' $DOM2 'CENTRE')) 1 'H2' 0.3 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'H2' 0.2
'NATU' 'DISCRET') ;
XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' zero
'NATURE' 'DISCRET' ;
XO21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
('DOMA' $DOM2 'CENTRE')) 1 'O2' 0.147 'NATU' 'DISCRET') 'ET'
('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'O2' 0.2
'NATU' 'DISCRET') ;
XN = XH21 'ET' XH2O1 'ET' XO21 ;
CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT '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 ;
*
***** 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 ';
*
* 'OPTION' 'SAUVER' 'gas_properties.sauv' ;
* 'SAUVER' PGAZ ;
*
**** We compute the mass fractions
*
MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2' Mh2
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'O2' Mo2
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' Mh2o
'NATURE' 'DISCRET') 'ET'
('MANUEL' 'CHPO' ('DOMA' $DOMTOT '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 ;
*
**** Burned gas
* Here x refers to the mole producted per 1 mole of unburned gas
* N.B. 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' $DOMTOT 'CENTRE') 1 'O2' zero ;
XH22a = XH21 '-' ('NOMC' 'H2' (2. '*' (XO21 '-' XO22a))) ;
XO22a = XO22a ;
XH22a = XH22a ;
XH22b = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT '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' 'Probleme dan le calcul du gaz brule' ;
'ERREUR' 21 ;
'FINSI' ;
*
**** xtot2 = mole of burnet 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 ;
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' 'Probleme dan le calcul du gaz brule' ;
'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 ;
*
**** The conservative variables
*
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'))
;
RN = PINI '/' (RINI '*' TINI) ;
* Internal energy
eini = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT '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 ;
RETN = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY')))) ;
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' $DOMTOT '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')) ;
'SI' VRAI ;
* Test
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' 'Probleme calcul energie totale' ;
'ERREUR' 21 ;
'FINSI' ;
'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
'MESSAGE' 'Probleme calcul energie totale' ;
'ERREUR' 21 ;
'FINSI' ;
'FINSI' ;
****************************************************************
****************************************************************
***** Parameters for the computations ******
****************************************************************
****************************************************************
* Upwind scheme
* METO = 'VLH' ;
METO = 'SS' ;
* Iterations
* Final time
* Safety Factor fot the time step
* Second order reconstruction?
NITER = 10000 ;
TFINAL = 10.0D-3 ;
SAFFAC = 0.7 ;
LOGSO = VRAI ;
*
* EPS_CC linked to CREBCOM criterion
*
EPS_CC = 0.8 ;
* 'OPTION' 'SAUVER' 'parameters.sauv' ;
* 'SAUVER' METO NITER TFINAL SAFFAC EPS_CC LOGSO ;
****************************************************************
****************************************************************
***** The computation ******
****************************************************************
****************************************************************
* Mesh dimension (in 2D!)
DELTAXN = ('DOMA' $DOMTOT 'VOLUME') '**' (1. '/' 2.) ;
YN = RYN '/' RN ;
*
**** We impose we have combustion in the ignition region (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.
* Since DELTAXN is constant, DELTATC does not change during the
* computation.
*
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
('REDU' RN ('DOMA' $DOM1 'CENTRE'))
('REDU' YN ('DOMA' $DOM1 'CENTRE'))
('REDU' YININ ('DOMA' $DOM1 'CENTRE'))
('REDU' YFINN ('DOMA' $DOM1 'CENTRE'))
('REDU' K0N ('DOMA' $DOM1 'CENTRE'))
('REDU' DELTAXN ('DOMA' $DOM1 'CENTRE'))
0.0 (1D4 '*' DELTATC) 0.3 ;
RYN = RYN '+' DELTARY ;
YN = RYN '/' RN ;
RETN = RETN '+' DELTARE ;
RN0 = 'COPIER' RN ;
GN0 = 'COPIER' GN ;
RETN0 = 'COPIER' RETN ;
RYN0 = 'COPIER' RYN ;
RSN0 = 'COPIER' RSN ;
*
**** Parameter for the time loop
*
* Names of the residuum components
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' 'H2IN'
'H2FI' 'K0' ;
*
**** Geometrical coefficient to compute gradients
*
GRADR CACCA COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
('MOTS' 'SCAL') RN ;
GRADV CACCA COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
('MOTS' 'UX' 'UY') GN ;
TPS = 0.0 ;
*
**** The chemistry time step
*
*
DT_CHEM = SAFFAC '*' DELTATC ;
*
**** Temporal loop
*
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
TNM1 = 'COPIER' TN ;
'TEMPS' 'ZERO' ;
'REPETER' BL1 NITER ;
*
**** Primitive variables
*
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN
TNM1 ;
TNM1 = 'COPIER' TN ;
'SI' LOGSO ;
GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') RN 'GRADGEO' COEFSCAL ;
GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'SCAL') PN 'GRADGEO' COEFSCAL ;
GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
('MOTS' 'UX' 'UY') VN 'GRADGEO' COEFVECT ;
GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'H2' 'O2' 'H2O') YN 'GRADGEO' COEFSCAL ;
GRADS ALS = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'H2IN' 'H2FI' 'K0') SN 'GRADGEO' COEFSCAL ;
ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 2 1
$DOMTOT PGAZ
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
YN GRADY ALY
SN GRADS ALS ;
'SINON' ;
ROF VITF PF YF SF = 'PRET' 'PERFTEMP' 1 1
$DOMTOT PGAZ
RN
VN
PN
YN
SN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
$DOMTOT PGAZ LISTINCO ROF VITF PF YF SF ;
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
*
**** Total time step
*
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM) ;
*
**** Increment of the variables (convection)
*
RESIDU = DTMIN '*' RESIDU ;
DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
DRYN = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
('MOTS' 'H2' 'O2' 'H2O') ;
DRSN = 'EXCO' ('MOTS' 'H2IN' 'H2FI' 'K0') RESIDU
('MOTS' 'H2IN' 'H2FI' 'K0') ;
TPS = TPS '+' DTMIN ;
RN = RN '+' DRN ;
GN = GN '+' DGN ;
RETN = RETN '+' DRETN ;
RYN = RYN '+' DRYN ;
RSN = RSN '+' DRSN ;
YN = RYN '/' RN ;
SN = RSN '/' RN ;
*
**** Increment of the variables (chemical source)
*
K0N = 'EXCO' 'K0' SN 'SCAL' ;
YININ = 'EXCO' 'H2IN' SN 'H2' ;
YFINN = 'EXCO' 'H2FI' SN 'H2' ;
DRETN DRYN = 'FLAM' 'CREBCOM2' $DOMTOT PGAZ LMOT1 LCOEF RN YN
YININ YFINN K0N DELTAXN EPS_CC DTMIN 0.3 ;
RYN = RYN '+' DRYN ;
RETN = RETN '+' DRETN ;
'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ;
'FINSI' ;
'SI' (TPS > TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
'TEMPS' ;
* 'OPTION' 'SAUVER' ('CHAINE' 'result.sauv_' RAF 'tps_' TPS ) ;
* 'SAUVER' ;
****************************************************************
****************************************************************
***** The post-treatment ******
****************************************************************
****************************************************************
* 'OPTI' 'REST' 'result.sauv_1tps_5' ;
* 'REST' ;
GRAPH = VRAI ;
GRAPH = FAUX ;
*
**** The mesh
*
'SI' GRAPH ;
'TRACER' ('DOMA' $DOMTOT 'MAILLAGE')
'TITR' ('CHAINE' 'Maillage: ' ('NBEL' DOMTOT) ' elts');
'FINSI' ;
*
**** Initial conditions
*
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
'SI' GRAPH ;
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP'
PGAZ RN0 GN0 RETN0 RYN0 RSN0 ;
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN0 ;
CHM_VN = 'KCHA' $DOMTOT 'CHAM' VN ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
CHM_YN = 'KCHA' $DOMTOT 'CHAM' YN ;
CHM_SN = 'KCHA' $DOMTOT 'CHAM' SN ;
'TRAC' CHM_RN MOD1
'TITR' ('CHAINE' 'rho at t=' 0.0) ;
'TRAC' CHM_VN MOD1
'TITR' ('CHAINE' 'v at t= ' 0.0) ;
'TRAC' CHM_PN MOD1
'TITR' ('CHAINE' 'p at t= ' 0.0) ;
'TRAC' CHM_TN MOD1
'TITR' ('CHAINE' 'T at t= ' 0.0) ;
'TRAC' CHM_YN MOD1
'TITR' ('CHAINE' 'y at t= ' 0.0) ;
'TRAC' CHM_SN MOD1
'TITR' ('CHAINE' 's at t= ' 0.0) ;
'FINSI' ;
*
**** Results
*
VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;
'SI' GRAPH ;
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_VN = 'KCHA' $DOMTOT 'CHAM' VN ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
CHM_YN = 'KCHA' $DOMTOT 'CHAM' YN ;
CHM_SN = 'KCHA' $DOMTOT 'CHAM' SN ;
'TRAC' CHM_RN MOD1
'TITR' ('CHAINE' 'rho at t=' TPS) ;
'TRAC' CHM_VN MOD1
'TITR' ('CHAINE' 'v at t= ' TPS) ;
'TRAC' CHM_PN MOD1
'TITR' ('CHAINE' 'p at t= ' TPS) ;
'TRAC' CHM_TN MOD1
'TITR' ('CHAINE' 'T at t= ' TPS) ;
'TRAC' CHM_YN MOD1
'TITR' ('CHAINE' 'y at t= ' TPS) ;
'TRAC' CHM_SN MOD1
'TITR' ('CHAINE' 's at t= ' TPS) ;
'FINSI' ;
*
*** Evolution Objects
*
SS = 'COORDONNEE' 1 LIGEVO;
evxx = 'EVOL' 'CHPO' ss LIGEVO ;
lxx = 'EXTRAIRE' evxx 'ORDO' ;
x0 = 'MINIMUM' lxx;
x1 = 'MAXIMUM' lxx ;
titolo = 'CHAINE' ' at t = ' TPS ' Methode = ' METO
' SF = ' SAFFAC ;
VNN = 'EXCO' 'UX' VN ;
VNT = 'EXCO' 'UY' VN ;
* rho
evro = 'EVOL' 'CHPO' RN LIGEVO ;
lro = 'EXTRAIRE' evro 'ORDO' ;
evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro ;
tro = 'CHAINE' 'ro ' titolo ;
* u
evu = 'EVOL' 'CHPO' VNN LIGEVO ;
lu = 'EXTRAIRE' evu 'ORDO' ;
evu = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
tu = 'CHAINE' 'u ' titolo ;
* yini
YININ = 'EXCO' 'H2IN' SN ;
evyini = 'EVOL' 'CHPO' YININ LIGEVO ;
lyini = 'EXTRAIRE' evyini 'ORDO' ;
evyini = 'EVOL' 'MANU' 'x' lxx 'yini' lyini ;
tyini = 'CHAINE' 'yini ' titolo ;
* yfin
YFINN = 'EXCO' 'H2FI' SN ;
evyfin = 'EVOL' 'CHPO' YFINN LIGEVO ;
lyfin = 'EXTRAIRE' evyfin 'ORDO' ;
evyfin = 'EVOL' 'MANU' 'x' lxx 'yfin' lyfin ;
tyfin = 'CHAINE' 'yfin ' titolo ;
* K0
K0N = 'EXCO' 'K0' SN ;
evk0 = 'EVOL' 'CHPO' K0N LIGEVO ;
lk0 = 'EXTRAIRE' evk0 'ORDO' ;
evk0 = 'EVOL' 'MANU' 'x' lxx 'k0' lk0 ;
tk0 = 'CHAINE' 'k0 ' titolo ;
* csi
CSIN = (('EXCO' 'H2' YN) '-' YININ) '/'
(YFINN '-' YININ) ;
evcsi = 'EVOL' 'CHPO' CSIN LIGEVO ;
lcsi = 'EXTRAIRE' evcsi 'ORDO' ;
evcsi = 'EVOL' 'MANU' 'x' lxx 'csi' lcsi ;
tcsi = 'CHAINE' 'csi ' titolo ;
* v
evv = 'EVOL' 'CHPO' VNT LIGEVO ;
lv = 'EXTRAIRE' evv 'ORDO' ;
evv = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
tv = 'CHAINE' 'v ' titolo ;
* p
evp = 'EVOL' 'CHPO' PN LIGEVO ;
lp = 'EXTRAIRE' evp 'ORDO' ;
evp = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
tp = 'CHAINE' 'p ' titolo ;
* T
evT = 'EVOL' 'CHPO' TN LIGEVO ;
lT = 'EXTRAIRE' evT 'ORDO' ;
evT = 'EVOL' 'MANU' 'x' lxx 'T' lT ;
tT = 'CHAINE' 'T ' titolo ;
'SI' GRAPH ;
'DESSIN' evro 'TITRE' tro 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evp 'TITRE' tp 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evu 'TITRE' tu 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evv 'TITRE' tv 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evT 'TITRE' tT 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evcsi 'TITRE' tcsi 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evk0 'TITRE' tk0 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evyini 'TITRE' tyini 'XBOR' x0 x1 'MIMA' ;
'DESSIN' evyfin 'TITRE' tyfin 'XBOR' x0 x1 'MIMA' ;
'FINSI' ;
* Test
PVN = 3.2E6 ;
PCJ = 1.1E6 ;
RHOVN = 3.7 ;
RHOCJ = 1.6 ;
PMAX = 'MAXIMUM' PN ;
'SI' ((PMAX '>' PVN) 'OU' (PMAX '<' PCJ)) ;
'MESSAGE' 'Pression = ???' ;
'ERREUR' 5 ;
'FINSI' ;
RHOMAX = 'MAXIMUM' RN ;
'SI' ((RHOMAX '>' RHOVN) 'OU' (RHOMAX '<' RHOCJ)) ;
'MESSAGE' 'Densite = ???' ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer