Télécharger konv_resi_dem_contact_discontinuty_22.dgibi
* fichier : konv_resi_dem_contact_discontinuity_22.dgibi ************************************************************************ ************************************************************************ *********************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la **** **** solution des **** **** Equations d'Euler pour un gaz parfait **** **** Approche DEM pour la combustion **** **** OPERATEURS PRET, KONV **** **** **** **** Consistency in the case of contact discontinuity. **** **** Left state = phase 2 **** **** Right state = phase 2 **** **** **** **** Methodes: SS, VLH **** **** **** **** A. BECCANTINI DM2S/SEMT/LTMF DECEMBRE 2009 **** *********************************************************** 'OPTION' 'ELEM' QUA4 ; 'OPTION' 'ECHO' 1 ; 'OPTION' 'TRAC' 'X' ; * *** GRAPH * GRAPH = FAUX ; * GRAPH = VRAI ; **** Cas homogene : la table proprieté de gaz ************************************************* **** 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 ; 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 ; *************************** ***** DOMAINE SPATIAL **** *************************** A1 = 0.0D0 0.01D0; A2 = 1.011D0 0.21D0; A3 = 2.12D0 0.021D0; A4 = 2.11D0 1.09D0; A5 = 1.023D0 1.132D0; A6 = 0.098D0 1.199D0; L12 = A1 'DROIT' 1 A2; L23 = A2 'DROIT' 1 A3; L34 = A3 'DROIT' 1 A4; L45 = A4 'DROIT' 1 A5; L56 = A5 'DROIT' 1 A6; L61 = A6 'DROIT' 1 A1; L25 = A2 'DROIT' 1 A5; DOM10 = 'DALL' L12 L25 L56 L61 'PLANE'; 'PLANE'; * *** Point ou on controlle la consistence * P10 = 1.0 0.5; * EPS = 1.0D-16 ; K0 = 10.0 ; * ******************************************* *** Etats gauche et droite **************** ******************************************* * * There are 3 fluxes to control : * 1 -> 2 can be checked by considering an exact reactive shock * 2 -> 1 can be checked by considering an exact reactive shock * 1 -> 1 can be checked via consistency (or via an exact non-reactive * or by comparing the results with a non-reactive one) * 2 -> 2 can be checked via consistency (or via an exact non-reactive * or by comparing the results with a non-reactive one) * * Case 1 -> 1 * Constant state * tg1 = 593.16 ; td1 = 293.16 ; pg1 = 1.023e5 ; pd1 = 1.023e5 ; ung1 = -122.0 ; und1 = -122.0 ; utg1 = 221.0 ; utd1 = 101.0 ; * utg1 = 0.0 ; * utd1 = 0.0 ; tg2 = 2800.15 ; td2 = 1800.15 ; pg2 = 1.023e5 ; pd2 = 1.023e5 ; ung2 = -150.0 ; und2 = -150.0 ; utg2 = 100.0 ; utd2 = 160.0 ; * * 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 ; * * 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) * '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 ; * * Computation of the conservative variables * RN1 = PN1 '/' (RGAS1 '*' TN1) ; RN2 = PN2 '/' (RGAS2 '*' TN2) ; GNX1 = RN1 * UN1 ; GNY1 = RN1 * UT1 ; GNX2 = RN2 * UN2 ; GNY2 = RN2 * UT2 ; ((und1 * und1) '+' (utd1 * utd1))); ((und2 * und2) '+' (utd2 * utd2))); RETN1 = RN1 '*' (ETHER1 '+' ECIN1 '+' EFORM1) ; RETN2 = RN2 '*' (ETHER2 '+' ECIN2 '+' EFORM2) ; * * Computation of the resi contribution in (n,t) * 'SI' (ung2 > 0) ; retg2 = 'EXTRAIRE' RETN2 1 ; rhog2 = 'EXTRAIRE' RN2 1 ; f1gd = ung2 '*' rhog2 ; f2gd = (f1gd '*' ung2) '+' pg2 ; f3gd = f1gd '*' utg2 ; f4gd = ung2 '*' (retg2 '+' pg2); 'SINON' ; retd2 = 'EXTRAIRE' RETN2 2 ; rhod2 = 'EXTRAIRE' RN2 2 ; f1gd = und2 '*' rhod2 ; f2gd = (f1gd '*' und2) '+' pd2 ; f3gd = f1gd '*' utd2 ; f4gd = und2 '*' (retd2 '+' pd2); 'FINSI' ; **************************************************** **************************************************** ******** Boucle sur les angles ********* **************************************************** **************************************************** DANGLE = 360 '/' 7.15; ANGLE = 11.3 ; 'REPETER' BLOC 8; * *** Rotation * ANGLE = ANGLE '+' DANGLE; ORIG = 0.0D0 0.0D0; 'MESSAGE' ; 'MESSAGE' (CHAIN 'Angle de rotation= ' ANGLE); 'MESSAGE' ; DOM1 = DOM10 'TOURNER' ANGLE ORIG; DOM2 = DOM20 'TOURNER' ANGLE ORIG; P1FAC = P10 'TOURNER' ANGLE ORIG; DOMTOT = DOM1 ET DOM2; 'ELIMINATION' DOMTOT 1D-6; $DOMTOT = 'MODELISER' DOMTOT 'EULER'; $DOM1 = 'MODELISER' DOM1 'EULER'; $DOM2 = 'MODELISER' DOM2 'EULER'; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMTOT = TDOMTOT . 'QUAF' ; * * **** Redefinition de P1FAC dans $DOMTOT 'FACE' * TX = -1 * NY ; TY = NX ; P1FAC ; AA = 'CHANGER' GEOPC 'POI1' ; * Si P3 n'existe pas, probleme en FACEL. 'SI' ('NEG' PC1 P2) ; ORIENT = -1 ; PCD = PC1 ; PCG = PC2 ; XVOLD = 'MAXIMUM' (TDOM1 . 'XXVOLUM') ; XVOLG = 'MAXIMUM' (TDOM2 . 'XXVOLUM') ; 'SINON' ; ORIENT = 1 ; PCD = PC2 ; PCG = PC1 ; XVOLG = 'MAXIMUM' (TDOM1 . 'XXVOLUM') ; XVOLD = 'MAXIMUM' (TDOM2 . 'XXVOLUM') ; 'FINSI' ; uxg1 = ((ung1 '*' NX) '+' (utg1 '*' TX)) '*' ORIENT ; uyg1 = ((ung1 '*' NY) '+' (utg1 '*' TY)) '*' ORIENT ; uxg2 = ((ung2 '*' NX) '+' (utg2 '*' TX)) '*' ORIENT ; uyg2 = ((ung2 '*' NY) '+' (utg2 '*' TY)) '*' ORIENT ; uxd1 = ((und1 '*' NX) '+' (utd1 '*' TX)) '*' ORIENT ; uyd1 = ((und1 '*' NY) '+' (utd1 '*' TY)) '*' ORIENT ; uxd2 = ((und2 '*' NX) '+' (utd2 '*' TX)) '*' ORIENT ; uyd2 = ((und2 '*' NY) '+' (utd2 '*' TY)) '*' ORIENT ; 'SI' GRAPH; 'TITRE' 'Domaine et FACEL'; 'FINSI' ; *********************** **** Les CHPOINTs **** *********************** 'UY' uyg1) '+' 'UY' uyd1) ; 'UY' uyg2) '+' 'UY' uyd2) ; CHGN1 = CHRN1 '*' CHVN1 ; CHGN2 = CHRN2 '*' CHVN2 ; R1 R2 V1 V2 P1 P2 T1 T2 = 'PRIM' 'DEM' PGAS CHAL1 CHAL2 (CHAL1 * CHRN1) (CHAL2 * CHRN2) (CHAL1 * CHGN1) (CHAL2 * CHGN2) (CHAL1 * CHRET1) (CHAL2 * CHRET2) CHTN1 CHTN2 EPS ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1' 0.0 ; 'P1' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1' 0.0 ; 'P1' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ; 'P1' 0.0 'P2' 0.0 ; 'P1' 0.0 'P2' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1' 0.0 ; 'P1' 0.0 ; * **** L'operateur 'PRET' * CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 = 'PRET' 'DEM' $DOMTOT CHAL1 (0.0 * GRADAL1) LIMAL1 CHAL2 (0.0 * GRADAL2) LIMAL2 R1 (0.0 * GRADR1) LIMR1 R2 (0.0 * GRADR2) LIMR2 V1 (0.0 * GRADV1) LIMV1 V2 (0.0 * GRADV2) LIMV2 P1 (0.0 * GRADP1) LIMP1 P2 (0.0 * GRADP2) LIMP2 ; * **** L'operateur 'KONV' * 'REPETER' BLMETO 2 ; 'SI' ('EGA' &BLMETO 1) ; METO = 'SS' ; 'FINSI' ; 'SI' ('EGA' &BLMETO 2) ; METO = 'AUSMPUP' ; 'FINSI' ; 'SI' ('EGA' &BLMETO 3) ; METO = 'VLH' ; * NB VLH n'est pas exact sur la discontinuité de contact 'FINSI' ; 'MESSAGE' ('CHAINE' 'METO = ' METO) ; 'ALF2' 'RN2' 'RUX2' 'RUY2' 'RET2') ; SI ('EGA' METO 'AUSMPUP') ; $DOMTOT PGAS LISTINC1 CHAL1 CHAL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS MAILIM VINF VINF ; 'SINON' ; $DOMTOT PGAS LISTINC1 CHAL1 CHAL2 CHFAL1 CHFAL2 CHFRN1 CHFRN2 CHFVN1 CHFVN2 CHFPN1 CHFPN2 K0 GRALP1 EPS MAILIM ; 'FINSI' ; * 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' ETHER1 1)) ; * 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' ECIN1 1)) ; * 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' EFORM1 1)) ; * 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' (pg1)) ; * 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' RN1 1)) ; * 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' (ung1)) ; * 'OPTION' DONN 5 ; RESX2N = (RESX2X '*' NX) '+' (RESX2Y * NY); RESX2T = (RESX2X '*' TX) '+' (RESX2Y * TY); f1 = 'EXTRAIRE' RESX1 'SCAL' PCD ; f1bis = 'EXTRAIRE' RESX1 'SCAL' PCG ; ERRO = 1D-8 '*' f1gd 'ABS' ; LOGI1 = ('ABS' ((f1 * XVOLD) + (f1bis * XVOLG))) < ERRO ; HORROR = ('ABS' ((f1gd * XSURF) + (f1bis * XVOLG))) ; LOGI2 = HORROR < ERRO ; 'MESSAGE' ('CHAINE' 'error = ' 'FORMAT' '(E16.10)' HORROR) ; 'MESSAGE' ('CHAINE' 'val = ' 'FORMAT' '(E16.10)' (f1gd 'ABS')) ; LOGI1 = LOGI1 'ET' LOGI2; f2 = 'EXTRAIRE' RESX2N 'SCAL' PCD ; f2bis = 'EXTRAIRE' RESX2N 'SCAL' PCG ; ERRO = 1D-8 '*' f2gd 'ABS' ; LOGI2 = ('ABS' ((f2 * XVOLD) + (f2bis * XVOLG))) < ERRO ; LOGI1 = LOGI1 'ET' LOGI2; HORROR = ('ABS' ((f2gd * XSURF) + (f2bis * XVOLG))) ; LOGI2 = HORROR < ERRO ; 'MESSAGE' ('CHAINE' 'error = ' 'FORMAT' '(E16.10)' HORROR) ; 'MESSAGE' ('CHAINE' 'val = ' 'FORMAT' '(E16.10)' (f2gd 'ABS')) ; LOGI1 = LOGI1 'ET' LOGI2; f3 = 'EXTRAIRE' RESX2T 'SCAL' PCD ; f3bis = 'EXTRAIRE' RESX2T 'SCAL' PCG ; ERRO = 1D-8 '*' f3gd 'ABS' ; LOGI2 = ('ABS' ((f3 * XVOLD) + (f3bis * XVOLG))) < ERRO ; LOGI1 = LOGI1 'ET' LOGI2; HORROR = ('ABS' ((f3gd * XSURF) + (f3bis * XVOLG))) ; LOGI2 = HORROR < ERRO ; 'MESSAGE' ('CHAINE' 'error = ' 'FORMAT' '(E16.10)' HORROR) ; 'MESSAGE' ('CHAINE' 'val = ' 'FORMAT' '(E16.10)' (f3gd 'ABS')) ; LOGI1 = LOGI1 'ET' LOGI2; f4 = 'EXTRAIRE' RESX3 'SCAL' PCD ; f4bis = 'EXTRAIRE' RESX3 'SCAL' PCG ; ERRO = 5D-7 '*' f4gd 'ABS' ; LOGI2 = ('ABS' ((f4 * XVOLD) + (f4bis * XVOLG))) < ERRO ; LOGI1 = LOGI1 'ET' LOGI2; HORROR = ('ABS' ((f4gd * XSURF) + (f4bis * XVOLG))) ; LOGI2 = HORROR < ERRO ; 'MESSAGE' ('CHAINE' 'error = ' 'FORMAT' '(E16.10)' HORROR) ; 'MESSAGE' ('CHAINE' 'val = ' 'FORMAT' '(E16.10)' (f4gd 'ABS')) ; LOGI1 = LOGI1 'ET' LOGI2; 'SI' ('NON' LOGI1); 'MESSAGE' ; 'MESSAGE' 'OPERATEUR KONV'; 'MESSAGE' ('CHAINE' METO); 'MESSAGE' ; 'ERREUR' 5 ; 'FINSI' ; 'FIN' BLMETO ; **************************************************** **************************************************** ******** Fin boucle sur les angles ********* **************************************************** **************************************************** 'FIN' BLOC; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales