* fichier : tubedeto3d2.dgibi
****************************************************************
***** PROPAGATION D UNE DETONATION DANS UN TUBE *****
***** MODELE COMBUSTION H2-air de PLEXUS *****
***** 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' 3
'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 0.0 ;
A2 = L1 0.0 0.0 ;
A3 = (L1 '+' L2) 0.0 0.0 ;
A4 = (L1 '+' L2) H1 0.0 ;
A5 = L1 H1 0.0 ;
A6 = 0.0 H1 0.0 ;
A7 = (L1 '+' L2 '+' L3) 0.0 0.0 ;
A8 = (L1 '+' L2 '+' L3) H1 0.0 ;
'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' ;
DOM1P = 'TRANSLATION' A2A5 N1 ((-1.0 '*' L1) 0.0 0.0) ;
DOM2P = 'TRANSLATION' A2A5 N2 (L2 0.0 0.0) ;
DOM3P = 'TRANSLATION' A3A4 N3 (L3 0.0 0.0) ;
* DOMTOT = total region
DOMTOTP = DOM1P 'ET' DOM2P 'ET' DOM3P ;
'ELIMINATION' (1.0D-3 '/' RAF) DOMTOTP ;
'ELIMINATION' DOMTOTP (1.0D-3 '/' RAF) (A1A2 'ET' A2A3 'ET'
A3A4 'ET' A4A5 'ET' A5A2 'ET' A5A6 'ET' A6A1 'ET' A3A7
'ET' A7A8 'ET' A8A4 'ET' A4A3) ;
'OPTION' 'ELEM' 'CUB8' ;
DOM1 = 'VOLUME' DOM1P 'TRANSLATION' NZ (0.0 0.0 H1) ;
DOM2 = 'VOLUME' DOM2P 'TRANSLATION' NZ (0.0 0.0 H1) ;
DOM3 = 'VOLUME' DOM3P 'TRANSLATION' NZ (0.0 0.0 H1) ;
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) (DX '/' 2) ;
P2 = ((L1 '+' L2 '+' L3) '-' (DX '/' 2)) (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. *****
****************************************************************
****************************************************************
* 'OPTION' 'ECHO' 1 'TRAC' 'X' ;
* RAF = 1 ;
* 'OPTION' 'REST' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'RESTITUER' ;
*
**** 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') 3 'UX' 0.0
'UY' 0.0 'UZ' 0.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 ;
* '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 ;
*
**** 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' 'UZ') ('MOTS' 'UX' 'UY' 'UZ')))) ;
GN = RN '*' WINI ;
YN = YH21 '+' YO21 '+' YH2O1 ;
RYN = RN '*' YN ;
'SI' VRAI ;
* Test
VN P T Y GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
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' ;
* 'OPTION' 'SAUV' 'ic.sauv';
* 'SAUVER' RN RYN GN RETN zero ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;
****************************************************************
****************************************************************
***** 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 = 5.0D-3 ;
SAFFAC = 0.7 ;
LOGSO = VRAI ;
* 'OPTION' 'SAUVER' 'parameters.sauv' ;
* 'SAUVER' METO NITER TFINAL SAFFAC LOGSO ;
****************************************************************
****************************************************************
***** The computation ******
****************************************************************
****************************************************************
************************
**** PROCEDURE DETO ****
************************
'DEBPROC' DETO ;
'ARGUMENT' PGAZ*'TABLE' TN*'CHPOINT' RYN*'CHPOINT' DTI*'FLOTTANT' ;
*
* dRYO2/dT = MO2 * A * (T ** -B) EXP( - TA / T)
* (RYH2 ** C) * (RYO2)
*
**** Les parametres physiques
*
* Ts = temperature de seuil
*
TS = 'MANUEL' 'CHPO' ('EXTRAIRE' TN 'MAILLAGE') 1 'SCAL' 600 ;
TA = 8310.;
A = 1.1725D14;
B = 0.91;
CON = 2;
H0H2 = PGAZ . 'H2 ' . 'H0K' ;
H0O2 = PGAZ . 'O2 ' . 'H0K' ;
H0H2O = PGAZ . 'H2O ' . 'H0K' ;
RYH2 = 'EXCO' 'H2' RYN 'H2' 'NATU' 'DISCRET' ;
RYO2 = 'EXCO' 'O2' RYN 'O2' 'NATU' 'DISCRET' ;
RYH2O = 'EXCO' 'H2O' RYN 'H2O' 'NATU' 'DISCRET' ;
RYH2 RYO2 RYH2O DELTAE = 'FLAM' 'ARRHENIU' TS A B CON TA
H0H2 H0O2 H0H2O DTI RYH2 RYO2 RYH2O
TN ;
RY = RYH2 'ET' RYO2 'ET' RYH2O ;
'FINPROC' RY DELTAE ;
****************************
**** FIN PROCEDURE DETO ****
****************************
*
**** The table of gas properties
* (PGAZ)
* 'OPTI' 'REST' 'gas_properties.sauv' ;
* 'REST' ;
*
**** The mesh
* RAF = refinement parameter
* DX = mesh dimension
* DOM1 = activation region
* DOM1 and DOM2 = zone 1
* DOM3 = zone 2
* DOMTOT = total region
* LIGEVO (line for evolutions)
* The MMODEL objects
* $DOMTOT
* $DOM1
* $DOM2
* $DOM3
* MOD1 = model object to plot 2D graphics
*
**** Initial conditions
*
* RYN : densities for H2, O2, H2O (kg/m^3)
* RN : total density (kg/m^3)
* GN : total momentum (kg/m^3 m/s)
* RETN : total energy (J/m^3)
*
* zero : "the numerical value of zero"
*
*
* 'OPTI' 'REST' 'ic.sauv' ;
* 'REST' ;
*
**** Parameters for the computation
*
* METO : upwind scheme
* NITER : max iterations number
* TFINAL : final time
* SAFFAC : safety factor for the time step
* LOGSO : second order reconstruction ?
* 'OPTION' 'REST' 'parameters.sauv' ;
* 'RESTITUER' ;
*
**** Chemical time step
*
* We compute d(RHY2)/dT at the von Neumann state
*
* Molar reaction rate
*
TVN = 2000 ;
TA = 8310.;
A = 1.1725D14;
B = 0.91;
CON = 2;
H0H2 = PGAZ . 'H2 ' . 'H0K' ;
H0O2 = PGAZ . 'O2 ' . 'H0K' ;
H0H2O = PGAZ . 'H2O ' . 'H0K' ;
RYNRED = 'REDU' RYN ('DOMA' $DOM2 'CENTRE') ;
RYH2 = 'MAXIMUM' ('EXCO' 'H2' RYNRED 'SCAL') ;
RYO2 = 'MAXIMUM' ('EXCO' 'O2' RYNRED 'SCAL') ;
OMEGA = A '*' (RYH2 '**' CON) '*' RYO2 '*' (TVN '**' (-1.0 '*' B))
'*' ('EXP' (-1.0 '*' TA '*' (TVN '**' -1))) ;
* OMEGA in (mol/m3/s)
* OMEGAH2 in (kg/m3/s)
OMEGAH2 = OMEGA '*' 2 '*' 2E-3 ;
DELTATC = RYH2 '/' OMEGAH2 ;
* NB DELTATC << DELTAT_CON
*
**** Zone d'activation = etat AICC
*
VN PN TN YN GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
RYH2 = ('EXCO' 'H2' RYN 'H2' 'NATU' 'DISCRET') ;
RYO2 = ('EXCO' 'O2' RYN 'O2' 'NATU' 'DISCRET') ;
RYH2O = ('EXCO' 'H2O' RYN 'H2O' 'NATU' 'DISCRET') ;
TN = TN '-' ('REDU' TN ('DOMA' $DOM1 'CENTRE')) '+'
('MANUEL' 'CHPO' ('DOMA' $DOM1 'CENTRE') 1 'SCAL' TVN) ;
TS = 'MANUEL' 'CHPO' ('EXTRAIRE' TN 'MAILLAGE') 1 'SCAL' (0.9 '*' TVN)
;
RYH2 RYO2 RYH2O DELTARE = 'FLAM' 'ARRHENIU' TS A B CON TA
H0H2 H0O2 H0H2O (1.0D5 '*' DELTATC)
RYH2 RYO2 RYH2O TN ;
RYN = RYH2 'ET' RYO2 'ET' RYH2O ;
RETN = RETN '+' DELTARE ;
RN0 = 'COPIER' RN ;
GN0 = 'COPIER' GN ;
RETN0 = 'COPIER' RETN ;
RYN0 = 'COPIER' RYN ;
*
**** Parameter for the time loop
*
* Names of the residuum components
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RUZ' 'RETN' 'H2' 'O2' 'H2O' ;
*
**** 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' 'UZ') GN ;
TPS = 0.0 ;
*
**** The chemistry time step is too little
* We take DT_CHEM > DELTATC
*
*
DT_CHEM = 50 '*' 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 GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
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' 'UZ') VN 'GRADGEO' COEFVECT ;
GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
('MOTS' 'H2' 'O2' 'H2O') YN 'GRADGEO' COEFSCAL ;
ROF VITF PF YF = 'PRET' 'PERFTEMP' 2 1
$DOMTOT PGAZ
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
YN GRADY ALY ;
'SINON' ;
ROF VITF PF YF = 'PRET' 'PERFTEMP' 1 1
$DOMTOT PGAZ
RN
VN
PN
YN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
$DOMTOT PGAZ LISTINCO ROF VITF PF YF ;
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' 'RUZ') RESIDU
('MOTS' 'UX' 'UY' 'UZ') ;
DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
DRYN = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
('MOTS' 'H2' 'O2' 'H2O') ;
TPS = TPS '+' DTMIN ;
RN = RN '+' DRN ;
GN = GN '+' DGN ;
RETN = RETN '+' DRETN ;
RYN = RYN '+' DRYN ;
*
**** Increment of the variables (chemical source)
*
VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
TNM1 ;
RYN DRETN = DETO PGAZ TN RYN DTMIN ;
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' DOMTOT 'TITR' ('CHAINE' 'Maillage: '
('NBEL' DOMTOT) ' elts');
'FINSI' ;
*
**** Initial conditions
*
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
'SI' GRAPH ;
VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP'
PGAZ RN0 GN0 RETN0 RYN0 ;
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 ;
'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) ;
'FINSI' ;
*
**** Results
*
VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
'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 ;
'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) ;
'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 ;
* 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' ;
'FINSI' ;
* Test
PVN = 3.2E6 ;
PCJ = 1.166 ;
RHOVN = 3.7 ;
RHOCJ = 1.4 ;
PMAX = 'MAXIMUM' PN ;
'SI' ((PMAX '>' PVN) 'OU' (PMAX '<' PCJ)) ;
'MESSAGE' 'Pression = ???' ;
'ERREUR' 5 ;
'FINSI' ;
RHOMAX = 'MAXIMUM' RN ;
'SI' ((RHOMAX '>' RHOVN) 'OU' (RHOMAX '<' (0.85 '*' RHOCJ))) ;
'MESSAGE' 'Densite = ???' ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer