* fichier : tubedeto2d1.dgibi ************************************************************************ * Section : Chimie Combustion ************************************************************************ **************************************************************** ***** 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' ;