Download prim_ther_dem3D.dgibi
* fichier : prim_ther_dem3D.dgibi
***********************************************************
**** FV "Cell-Centred Formulation" for the solution of ****
**** the Euler equations for a thermally perfect gas. ****
**** Discrete Equation Method for the propagation of ****
**** infinitely thin flames in initially homogeneous ****
**** media. ****
**** PRIM operator ****
**** It computes primitive variables from conservative ****
**** ones. ****
**** 3D case ****
**** ****
**** A. BECCANTINI DEN/DM2S/SFME/LTMF DECEMBER 2010 ****
***********************************************************
*
* alpha1 = 1 -> unburnt phase
* alpha2 = 1 -> burnt phase
* In exact algebra, alpha2 = 1 - alpha1
* alpha2 can be seen as a progress variable
*
****************************
**** PRODECURE PRIMCONS ****
****************************
*
'DEBPROC' PRIMCONS ;
'ARGUMENT' PGAS*TABLE TN1*'CHPOINT ' TN2*'CHPOINT '
PN1*'CHPOINT ' PN2*'CHPOINT ' VN1*'CHPOINT ' VN2*'CHPOINT ' ;
*
* ETHER = int_0^T cv(T') dT' T < TMAX
* = int_0^TMAX cv(T') dT' '+'
* cv(TMAX) T >= TMAX
ESP1 = 'EXTRAIRE' (PGAS . 'SPECIES') 1 ;
* DY1 = y_i - y_f for the species 1
DY1 = (('EXTRAIRE' (PGAS . 'MASSFRA') 1) '-'
('EXTRAIRE' (PGAS . 'MASSFRA') 2)) ;
COEF1 = ('EXTRAIRE' (PGAS . 'CHEMCOEF') 1) '*'
(PGAS . ESP1 . 'W') ;
YFINPH1 = 1.0 ;
YFINPH2 = 1.0 ;
'SI' (COEF1 > 0) ;
YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
YPH1 = YPH2 '+' DY1 ;
'SINON' ;
YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') 2 ;
YPH2 = YPH1 '-' DY1 ;
'FINSI' ;
YFINPH1 = YFINPH1 '-' YPH1 ;
YFINPH2 = YFINPH2 '-' YPH2 ;
PRYPH1 = 'PROG' YPH1 ;
PRYPH2 = 'PROG' YPH2 ;
'REPETER' BLESP (('DIME' (PGAS . 'SPECIES')) '-' 2) ;
ESP = 'EXTRAIRE' (PGAS . 'SPECIES') (&BLESP '+' 1) ;
COEF = ('EXTRAIRE' (PGAS . 'CHEMCOEF') (&BLESP '+' 1))
'*' (PGAS . ESP . 'W') ;
DY = (DY1 * (COEF '/' COEF1)) ;
'SI' (COEF > 0) ;
YPH2 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
YPH1 = YPH2 '+' DY ;
'SINON' ;
YPH1 = 'EXTRAIRE' (PGAS . 'MASSFRA') (&BLESP '+' 2) ;
YPH2 = YPH1 '-' DY ;
'FINSI' ;
PRYPH1 = PRYPH1 'ET' ('PROG' YPH1) ;
PRYPH2 = PRYPH2 'ET' ('PROG' YPH2) ;
YFINPH1 = YFINPH1 '-' YPH1 ;
YFINPH2 = YFINPH2 '-' YPH2 ;
'FIN' BLESP ;
PRYPH1 = PRYPH1 'ET' ('PROG' YFINPH1) ;
PRYPH2 = PRYPH2 'ET' ('PROG' YFINPH2) ;
'LISTE' PRYPH1 ;
'LISTE' PRYPH2 ;
*
TMAX = (PGAS . 'TMAX') ;
* TCAL1 = MIN TN1, TMAX
TCAL1 = 0.5D0 '*' ((TMAX '+' TN1) '-' ('ABS' (TN1 '-' TMAX))) ;
DTN1 = TN1 '-' TCAL1 ;
* TCAL1 = MIN TN1, TMAX
TCAL2 = 0.5D0 '*' ((TMAX '+' TN2) '-' ('ABS' (TN2 '-' TMAX))) ;
DTN2 = TN2 '-' TCAL2 ;
*
* Internal energy (J/kg in SI)
*
ETHER1 = 0.0 ;
CV1 = 0.0 ;
ETHER2 = 0.0 ;
CV2 = 0.0 ;
FUNTN1 = 1.0 ;
FUNTN2 = 1.0 ;
'REPETER' BLPO ((PGAS . 'NORD') '+' 1) ;
'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
AA = 'EXTRAIRE' (PGAS . ESP . 'A') &BLPO ;
DCV1 = (AA * YCEL1 * FUNTN1) ;
DCV2 = (AA * YCEL2 * FUNTN2) ;
CV1 = CV1 '+' DCV1 ;
CV2 = CV2 '+' DCV2 ;
ETHER1 = ETHER1 '+' (DCV1 * TCAL1 '/' (&BLPO)) ;
ETHER2 = ETHER2 '+' (DCV2 * TCAL2 '/' (&BLPO)) ;
'FIN' BLESP ;
FUNTN1 = FUNTN1 '*' TCAL1 ;
FUNTN2 = FUNTN2 '*' TCAL2 ;
'FIN' BLPO ;
ETHER1 = ETHER1 '+' (CV1 '*' DTN1) ;
ETHER2 = ETHER2 '+' (CV2 '*' DTN2) ;
*
* Formation energy/enthalpy (J/kg in SI) and gas constant (J/kg/K)
*
EFORM1 = 0.0 ;
EFORM2 = 0.0 ;
RGAS1 = 0.0 ;
RGAS2 = 0.0 ;
'REPETER' BLESP ('DIME' (PGAS . 'SPECIES')) ;
ESP = 'EXTRAIRE' (PGAS . 'SPECIES') &BLESP ;
YCEL1 = 'EXTRAIRE' PRYPH1 &BLESP ;
YCEL2 = 'EXTRAIRE' PRYPH2 &BLESP ;
EFORM1 = EFORM1 '+' (YCEL1 * (PGAS . ESP . 'H0K')) ;
EFORM2 = EFORM2 '+' (YCEL2 * (PGAS . ESP . 'H0K')) ;
RGAS1 = RGAS1 '+' (YCEL1 * (PGAS . 'RUNIV') '/'
(PGAS . ESP . 'W')) ;
RGAS2 = RGAS2 '+' (YCEL2 * (PGAS . 'RUNIV') '/'
(PGAS . ESP . 'W')) ;
'FIN' BLESP ;
'LISTE' RGAS1 ;
'LISTE' RGAS2 ;
*
* Computation of the conservative variables
*
RN1 = PN1 '/' (RGAS1 '*' TN1) ;
RN2 = PN2 '/' (RGAS2 '*' TN2) ;
GN1 = RN1 * VN1 ;
GN2 = RN2 * VN2 ;
LVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
ECIN1 = 0.5D0 '*' ('PSCAL' VN1 VN1 LVEL LVEL) ;
ECIN2 = 0.5D0 '*' ('PSCAL' VN2 VN2 LVEL LVEL) ;
RETN1 = RN1 '*' (ETHER1 '+' ECIN1 '+' EFORM1) ;
RETN2 = RN2 '*' (ETHER2 '+' ECIN2 '+' EFORM2) ;
*
'RESPRO' RN1 RN2 GN1 GN2 RETN1 RETN2 ;
'FINPROC' ;
*
********************************
**** FIN PRODECURE PRIMCONS ****
********************************
*
*
'OPTION' 'DIME' 3 ;
'OPTION' 'ELEM' 'CUB8' ;
'OPTION' 'ECHO' 1 ;
'OPTION' 'TRAC' 'X';
*
**** GRAPH
*
GRAPH = FAUX ;
* GRAPH = VRAI ;
*
*** We consider a mixture of H2, O2, H2O, N2
*
*************************************************
**** The table for the properties of the gas ****
*************************************************
*
PGAS = 'TABLE' ;
*
**** Order of the polynomial order for cv = cv(T)
* For T > TMAX, cv(T) = cv(Tmax)
*
PGAS . 'TMAX' = 6000.0 ;
PGAS . 'NORD' = 4 ;
*
**** Species involved in the mixture (before or after
* the chemical reaction)
*
PGAS . 'SPECIES' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' 'N2 ' ;
*
*
**** Coefficient of the chemical reaction.
* Note that for the first species this coefficient should be positive
* Normal, we take it equal to 1.
*
* H2 '+' 0.5 O2 ---> H2O
*
PGAS . 'CHEMCOEF' = 'PROG' 1.0 0.5 -1.0 0.0 ;
*
**** Mass fraction of the first species before and after the combustion
* Final mass fractions of the species with positive coefficients.
* Final mass fractions of the species with non-positive coefficient.
* The mass fraction of the last species is not given.
*
PGAS . 'MASSFRA' = 'PROG' 0.285219E-01 0.964039E-11 0.765104E-10
0.127442E-10 ;
*
**** Coef with the gas properties
*
PGAS . 'H2 ' = 'TABLE' ;
PGAS . 'H2O ' = 'TABLE' ;
PGAS . 'N2 ' = 'TABLE' ;
PGAS . 'O2 ' = 'TABLE' ;
*
**** Runiv (J/mole/K)
*
PGAS . 'RUNIV' = 8.31441 ;
*
**** W (kg/mole). Gas constant (J/kg/K = Runiv/W)
*
PGAS . 'H2 ' . 'W' = 2.016E-3 ;
PGAS . 'O2 ' . 'W' = 31.999E-3 ;
PGAS . 'H2O ' . 'W' = 18.0155E-3 ;
PGAS . 'N2 ' . 'W' = 28.013E-3 ;
*
**** Polynomial coefficients
*
PGAS . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAS . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
-1.82753232E-08 2.44485692E-12 ;
PGAS . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
8.78233606E-09 -3.05514485E-13 ;
PGAS . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
*
**** Formation enthalpies (energies) at 0K (J/Kg)
*
PGAS . 'H2 ' . 'H0K' = -4.195D6 ;
PGAS . 'H2O ' . 'H0K' = -1.395D7 ;
PGAS . 'N2 ' . 'H0K' = -2.953D5 ;
PGAS . 'O2 ' . 'H0K' = -2.634D5 ;
*************************************************
**** The mesh ****
*************************************************
A1 = 0.0D0 0.0D0 0.0D0 ;
A2 = 0.0D0 1.0D0 0.0D0 ;
A3 = 1.0D0 0.0D0 0.0D0 ;
A4 = 1.0D0 1.0D0 0.0D0 ;
A5 = 2.0D0 0.0D0 0.0D0 ;
A6 = 2.0D0 1.0D0 0.0D0 ;
MAIL1 = 'MANUEL' 'POI1' A1 ;
MAIL1 = (MAIL1 'ET' A2) 'ET' A3 ;
MAIL2 = 'MANUEL' 'POI1' A4 ;
MAIL2 = (MAIL2 'ET' A5) 'ET' A6 ;
MAIL = MAIL1 'ET' MAIL2 ;
'SI' GRAPH ;
'TRACER' MAIL 'TITRE' 'Domaine' ;
'FINSI' ;
*************************************************
**** Initial condition on primitive variables ***
*************************************************
TN1 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 293.16) '+'
('MANUEL' 'CHPO' MAIL2 1 'SCAL' 283.16) ;
TN2 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 2800.15) '+'
('MANUEL' 'CHPO' MAIL2 1 'SCAL' 2910.15) ;
PN1 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 1.013E5) '+'
('MANUEL' 'CHPO' MAIL2 1 'SCAL' 1.011E5) ;
PN2 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 1.023E5) '+'
('MANUEL' 'CHPO' MAIL2 1 'SCAL' 1.033E5) ;
VN1 = ('MANUEL' 'CHPO' MAIL1 3 'UX' 122.
'UY' 101. 'UZ' 0.0 ) '+'
('MANUEL' 'CHPO' MAIL2 3 'UX' 120.
'UY' 102. 'UZ' 0.0 ) ;
VN2 = ('MANUEL' 'CHPO' MAIL1 3 'UX' 150.
'UY' 160. 'UZ' 0.0 ) '+'
('MANUEL' 'CHPO' MAIL2 3 'UX' 140.
'UY' 160. 'UZ' 0.0 ) ;
ALPHA1 = ('MANUEL' 'CHPO' MAIL1 1 'SCAL' 0.9) '+'
('MANUEL' 'CHPO' MAIL2 1 'SCAL' 0.1) ;
ALPHA2 = -1 '*' (ALPHA1 '-' 1.0) ;
*
* Computation of the conservative variables
*
RN1 RN2 GN1 GN2 RETN1 RETN2 = PRIMCONS
PGAS TN1 TN2 PN1 PN2 VN1 VN2 ;
*
ARN1 = ALPHA1 * RN1 ;
ARN2 = ALPHA2 * RN2 ;
AGN1 = ALPHA1 * GN1 ;
AGN2 = ALPHA2 * GN2 ;
ARETN1 = ALPHA1 * RETN1 ;
ARETN2 = ALPHA2 * RETN2 ;
*
* K0 = 22.5 ;
EPS = 1.0D-16 ;
*
TG1 = TN1 '*' 10. ;
TG2 = TN2 '/' 10. ;
R1 R2 V1 V2 P1 P2 T1 T2
= 'PRIM' 'DEM' PGAS ALPHA1 ALPHA2
ARN1 ARN2 AGN1 AGN2
ARETN1 ARETN2 TG1 TG2 EPS ;
*
ERRO1 = 'MAXIMUM' ((TN1 '-' T1) '/' TN1) 'ABS' ;
ERRO2 = 'MAXIMUM' ((TN2 '-' T2) '/' TN2) 'ABS' ;
ERRO3 = 'MAXIMUM' ((PN2 '-' P2) '/' PN2) 'ABS' ;
ERRO4 = 'MAXIMUM' ((PN1 '-' P1) '/' PN1) 'ABS' ;
ERRO5 = 'MAXIMUM' ((RN1 '-' R1) '/' RN1) 'ABS' ;
ERRO6 = 'MAXIMUM' ((RN2 '-' R2) '/' RN2) 'ABS' ;
ERRO = 'MAXIMUM' ('PROG' ERRO1 ERRO2 ERRO3 ERRO4
ERRO5 ERRO6) ;
'SI' (ERRO > 1.0D-12) ;
'ERREUR' 5 ;
'FINSI' ;
*
*** We consider a mixture of H2, O2, H2O
*
*************************************************
**** The table for the properties of the gas ****
*************************************************
*
PGAS2 = 'TABLE' ;
*
**** Order of the polynomial order for cv = cv(T)
* For T > TMAX, cv(T) = cv(Tmax)
*
PGAS2 . 'TMAX' = 6000.0 ;
PGAS2 . 'NORD' = 4 ;
*
**** Species involved in the mixture (before or after
* the chemical reaction)
*
PGAS2 . 'SPECIES' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
*
*
**** Coefficient of the chemical reaction.
* Note that for the first species this coefficient should be positive
* Normal, we take it equal to 1.
*
* H2 '+' 0.5 O2 ---> H2O
*
PGAS2 . 'CHEMCOEF' = 'PROG' 1.0 0.5 -1.0 ;
*
**** Coef with the gas properties
*
PGAS2 . 'H2 ' = 'TABLE' ;
PGAS2 . 'H2O ' = 'TABLE' ;
PGAS2 . 'O2 ' = 'TABLE' ;
*
**** Runiv (J/mole/K)
*
PGAS2 . 'RUNIV' = 8.31441 ;
*
**** W (kg/mole). Gas constant (J/kg/K = Runiv/W)
*
PGAS2 . 'H2 ' . 'W' = 2. * 1.00797E-3 ;
PGAS2 . 'O2 ' . 'W' = 2. * 15.9994E-3 ;
PGAS2 . 'H2O ' . 'W' = (PGAS2 . 'H2 ' . 'W' ) '+'
(0.5 * (PGAS2 . 'O2 ' . 'W' )) ;
*
**** Mass fraction of the first species before and after the combustion
* Final mass fractions of the species with positive coefficients.
* Final mass fractions of the species with non-positive coefficient.
* The mass fraction of the last species is not given.
*
NH2 = 1.;
NO2 = 0.4999999999 * NH2 ;
NH2O = 1.0D-12 ;
MH2 = NH2 * (PGAS2 . 'H2 ' . 'W') ;
MO2 = NO2 * (PGAS2 . 'O2 ' . 'W') ;
MH2O = NH2O * (PGAS2 . 'H2O ' . 'W') ;
MTOT = MH2 '+' MO2 '+' MH2O ;
YH2 = MH2 '/' MTOT ;
YO2 = MO2 '/' MTOT ;
YH2O = MH2O '/' MTOT ;
SI FAUX ;
*
**** Complete combustion
*
SI (NH2 > (2.0 * NO2)) ;
YO2F = 1.0D-12 ;
YH2F = YH2 '-' ((YO2 '-' YO2F) * ((PGAS2 . 'H2 ' . 'W') '/'
(0.5 * (PGAS2 . 'O2 ' . 'W')))) ;
'SINON' ;
YH2F = 1.0D-12 ;
YO2F = YO2 '-' (0.5 * (YH2 '-' YH2F) * ((PGAS2 . 'O2 ' . 'W')
'/'(PGAS2 . 'H2 ' . 'W'))) ;
'FINSI' ;
'SINON' ;
YH2F = 1. '/' 3. * YH2 ;
YO2F = YO2 '-' (0.5 * (YH2 '-' YH2F) * ((PGAS2 . 'O2 ' . 'W') '/'
(PGAS2 . 'H2 ' . 'W'))) ;
'FINSI' ;
YH2OF = YH2O '+' (YH2 '-' YH2F) '+' (YO2 '-' YO2F) ;
ERRO = ((YH2F '+' YO2F '+' YH2OF) '-' 1.0) ;
SI (ERRO > 1.0D-12) ;
'ERREUR' 21 ;
'FINSI' ;
*
PGAS2 . 'MASSFRA' = 'PROG' YH2 YH2F YO2F ;
*
**** Polynomial coefficients
*
PGAS2 . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAS2 . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
-1.82753232E-08 2.44485692E-12 ;
PGAS2 . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
*
**** Formation enthalpies (energies) at 0K (J/Kg)
*
PGAS2 . 'H2 ' . 'H0K' = -4.195D6 ;
PGAS2 . 'H2O ' . 'H0K' = -1.395D7 ;
PGAS2 . 'O2 ' . 'H0K' = -2.634D5 ;
*
*
*************************************************
**** The initial conditions *********************
*************************************************
*
eps = 1.0D-64 ;
K0 = 3000. ;
*
T1 = 473. ;
* alph11 = 1.0 '-' 1.0D-5 ;
alph11 = 1.0 '-' 5.0D-3 ;
* alph12 = 1.0D-5 ;
alph12 = 5.0D-3 ;
alph21 = 1.0 '-' alph11 ;
alph22 = 1.0 - alph12 ;
T2 = 4000. ;
un1 = 1900. ;
un2 = 1200. ;
ut1 = 1221. ;
ut2 = 1432. ;
utt1 = 1231. ;
utt2 = 1439. ;
pre1 = 6.909D5 ;
pre2 = 6.909D5 ;
*
CHVN1 = ('MANUEL' 'CHPO' MAIL 3 'UX' un1
'UY' ut1 'UZ' utt1) ;
CHVN2 = ('MANUEL' 'CHPO' MAIL 3 'UX' un2
'UY' ut2 'UZ' utt2) ;
*
CHTN1 = ('MANUEL' 'CHPO' MAIL 1 'SCAL' T1) ;
CHTN2 = ('MANUEL' 'CHPO' MAIL 1 'SCAL' T2) ;
*
CHPN1 = ('MANUEL' 'CHPO' MAIL 1 'SCAL' pre1) ;
CHPN2 = ('MANUEL' 'CHPO' MAIL 1 'SCAL' pre2) ;
*
DOMAL = MAIL1 ;
DOMRES = 'DIFF' MAIL DOMAL ;
CHAL1 = ('MANUEL' 'CHPO' DOMRES 1 'SCAL'
alph11) '+'
('MANUEL' 'CHPO' DOMAL 1 'SCAL'
alph12) ;
CHAL2 = ('MANUEL' 'CHPO' DOMRES 1 'SCAL'
alph21) '+'
('MANUEL' 'CHPO' DOMAL 1 'SCAL'
alph22) ;
*
CHRN1 CHRN2 CHGN1 CHGN2 CHRET1 CHRET2 = PRIMCONS
PGAS2 CHTN1 CHTN2 CHPN1 CHPN2 CHVN1 CHVN2 ;
*
RN1 RN2 VN1 VN2 PN1 PN2 TN1 TN2 = 'PRIM' 'DEM' PGAS2 CHAL1 CHAL2
(CHAL1 * CHRN1) (CHAL2 * CHRN2) (CHAL1 * CHGN1)
(CHAL2 * CHGN2) (CHAL1 * CHRET1) (CHAL2 * CHRET2)
CHTN1 CHTN2 EPS ;
TN1M = COPIER TN1 ;
TN2M = COPIER TN2 ;
*
ERRO = 'MAXIMUM' (PRE1 '-' ('REDU' PN1 DOMAL)) 'ABS' ;
ERRO = ERRO '/' (pre1) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
ERRO = 'MAXIMUM' (PRE2 '-' ('REDU' PN2 DOMRES)) 'ABS' ;
ERRO = ERRO '/' (pre2) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
ERRO = 'MAXIMUM' (T1 '-' ('REDU' TN1 DOMAL)) 'ABS' ;
ERRO = ERRO '/' (T1) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
ERRO = 'MAXIMUM' (T2 '-' ('REDU' TN2 DOMRES)) 'ABS' ;
ERRO = ERRO '/' (T2) ;
'SI' (ERRO > 1.0D-6) ;
'MESSAGE' ('CHAINE' 'Erreur PRIM ' ERRO);
'ERREUR' 21 ;
'FINSI' ;
*
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer