Télécharger prim_ther_dem.dgibi
* fichier : prim_ther_dem.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. **** **** **** **** A. BECCANTINI DEN/DM2S/SFME/LTMF JUILLET 2007 **** *********************************************************** * * Conservative variables Primitive variables * * u(1,*) = alpha1 w(1,*) = alpha1 * u(2,*) = alpha1 rho1 w(2,*) = rho1 * u(3,*) = alpha1 rhou1 w(3,*) = u1 * u(4,*) = alpha1 rhov1 w(4,*) = v1 * u(5,*) = alpha1 rhoet1 w(5,*) = p1 * u(6,*) = alpha2 w(6,*) = alpha2 * u(7,*) = alpha2 rho2 w(7,*) = rho2 * u(8,*) = alpha2 rhou2 w(8,*) = u2 * u(9,*) = alpha2 rhov2 w(9,*) = v2 * u(10,*) = alpha2 rhoet2 w(10,*) = p2 * * 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 ; 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' ; YFINPH1 = YFINPH1 '-' YPH1 ; YFINPH2 = YFINPH2 '-' YPH2 ; 'FIN' BLESP ; '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) ; 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 ; 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 ; 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' 'ELEM' 'QUA4' ; '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) * * * **** 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 * * **** 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. * 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 * -2.37281455E-07 1.84701105E-11 ; -1.82753232E-08 2.44485692E-12 ; 8.78233606E-09 -3.05514485E-13 ; 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; A2 = 0.0D0 1.0D0; A3 = 1.0D0 0.0D0; A4 = 1.0D0 1.0D0; A5 = 2.0D0 0.0D0; A6 = 2.0D0 1.0D0; L12 = A1 'DROIT' 4 A2 ; L34 = A3 'DROIT' 4 A4 ; L56 = A5 'DROIT' 4 A6 ; DOM1 = L12 'REGLER' 5 L34 'COULEUR' 'ROUG' ; DOM2 = L34 'REGLER' 5 L56 'COULEUR' 'VERT' ; DOMTOT = DOM1 'ET' DOM2 ; 'SI' GRAPH ; 'TRACER' DOMTOT 'TITRE' 'Domaine' ; 'FINSI' ; MDOMTOT = TDOMTOT . 'QUAF' ; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MAIL1 = TDOM1 . 'CENTRE' ; MAIL2 = TDOM2 . 'CENTRE' ; ************************************************* **** Initial condition on primitive variables *** ************************************************* 'UY' 101. ) '+' 'UY' 102. ) ; 'UY' 160. ) '+' 'UY' 160.) ; 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' ; 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) * * * **** 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 * * **** 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' ; * * **** Polynomial coefficients * -2.37281455E-07 1.84701105E-11 ; -1.82753232E-08 2.44485692E-12 ; 2.33636971E-08 -1.53304905E-12; 8.78233606E-09 -3.05514485E-13 ; * **** 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 = 0.0 ; un2 = 0.0 ; ut1 = 0.0 ; ut2 = 0.0 ; pre1 = 6.909D5 ; pre2 = 6.909D5 ; * 'UY' ut1) ; 'UY' ut2) ; * * * DOMAL = TDOM1 . 'CENTRE' ; alph11) '+' alph12) ; alph21) '+' 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 = ERRO '/' (pre1) ; 'SI' (ERRO > 1.0D-6) ; 'ERREUR' 21 ; 'FINSI' ; ERRO = ERRO '/' (pre2) ; 'SI' (ERRO > 1.0D-6) ; 'ERREUR' 21 ; 'FINSI' ; ERRO = ERRO '/' (T1) ; 'SI' (ERRO > 1.0D-6) ; 'ERREUR' 21 ; 'FINSI' ; ERRO = ERRO '/' (T2) ; 'SI' (ERRO > 1.0D-6) ; 'ERREUR' 21 ; 'FINSI' ; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales