* fichier : prim_ther_2es.dgibi
**********************************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la solution des ****
**** Equations d'Euler pour un gaz parfait ****
**** OPERATEUR PRIM ****
**** Gaz multi-especes (deux especes) "thermally perfect" ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/LTMF JANVIER 2000 ****
**********************************************************************
*
'OPTION' 'DIME' 2 'ELEM' QUA4
'ECHO' 0 'TRAC' 'X';
*
*
**** GRAPH
*
*
GRAPH = FAUX ;
* GRAPH = VRAI ;
*
*
*** On considere une melange H2, O2
*
* De la table JANAF on a, pour les cp (J/Kg/K)
*
LTEMP = 'PROG' .2000D+03 .3000D+03 .4000D+03 .5000D+03
.6000D+03
.7000D+03 .8000D+03 .9000D+03 .1000D+04 .1100D+04
.1200D+04 .1300D+04 .1400D+04 .1500D+04 .1600D+04
.1700D+04 .1800D+04 .1900D+04 .2000D+04 .2100D+04
.2200D+04 .2300D+04 .2400D+04 .2500D+04 .2600D+04 ;
LTEMP = LTEMP 'ET' ('PROG'
.2700D+04 .2800D+04 .2900D+04 .3000D+04 .3100D+04
.3200D+04 .3300D+04 .3400D+04 .3500D+04 .3600D+04
.4000D+04 .4500D+04 .5000D+04 .5500D+04 .6000D+04) ;
LCVH2 = 'PROG' .9944D+04 .1018D+05 .1035D+05 .1039D+05
.1042D+05
.1048D+05 .1057D+05 .1070D+05 .1086D+05 .1104D+05
.1125D+05 .1147D+05 .1168D+05 .1190D+05 .1211D+05
.1231D+05 .1251D+05 .1270D+05 .1288D+05 .1305D+05
.1321D+05 .1337D+05 .1351D+05 .1365D+05 .1378D+05 ;
LCVH2 = LCVH2 'ET' ('PROG'
.1391D+05 .1403D+05 .1415D+05 .1426D+05 .1438D+05
.1449D+05 .1460D+05 .1471D+05 .1480D+05 .1490D+05
.1528D+05 .1573D+05 .1613D+05 .1646D+05 .1669D+05 ) ;
LCVO2 = 'PROG' .6500D+03 .6580D+03 .6810D+03 .7120D+03
.7430D+03
.7710D+03 .7940D+03 .8140D+03 .8300D+03 .8430D+03
.8550D+03 .8650D+03 .8740D+03 .8820D+03 .8900D+03
.8970D+03 .9040D+03 .9110D+03 .9180D+03 .9250D+03
.9320D+03 .9390D+03 .9450D+03 .9520D+03 .9580D+03 ;
LCVO2 = LCVO2 'ET' ('PROG'
.9640D+03 .9700D+03 .9760D+03 .9810D+03 .9870D+03
.9920D+03 .9970D+03 .1001D+04 .1006D+04 .1010D+04
.1034D+04 .1054D+04 .1074D+04 .1097D+04 .1127D+04 ) ;
'MESSAGE' ;
'MESSAGE' ('cv_i = polynoms de 4-eme degree');
'MESSAGE' ;
ERRMAX = 1.0D-6 ;
*) Degre de polynome = 4
***********************
**** LA TABLE PGAZ ****
***********************
PGAZ = 'TABLE' ;
*
**** Ordre des polynoms
*
PGAZ . 'NORD' = 4 ;
*
**** Especes qui sont dans les equations d'Euler
*
PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' ;
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'O2 ' ;
*
PGAZ . 'H2 ' = 'TABLE' ;
PGAZ . 'O2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'H2 ' . 'R' = 4130.0 ;
PGAZ . 'O2 ' . 'R' = 259.8 ;
*
**** Regressions polynomials
*
PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
*
**** "Enthalpies" (ou energies) de formations a OK (J/Kg)
*
PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
*
*** Fin PGAZ
*
A0H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 1 ;
A1H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 2 ;
A2H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 3 ;
A3H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 4 ;
A4H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 5 ;
A0O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 1 ;
A1O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 2 ;
A2O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 3 ;
A3O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 4 ;
A4O2 = 'EXTRAIRE' (PGAZ . 'O2 ' . 'A') 5 ;
LCV1H2 = 'PROG' ;
LCV1O2 = 'PROG' ;
'REPETER' BL1 ('DIME' LTEMP ) ;
T = 'EXTRAIRE' LTEMP &BL1 ;
T2 = T * T ;
T3 = T2 * T ;
T4 = T3 * T;
LCV1H2 = LCV1H2 'ET' ('PROG'
(A0H2 '+' (A1H2 * T) '+' (A2H2 * T2) '+' (A3H2 * T3) '+'
(A4H2 * T4))) ;
LCV1O2 = LCV1O2 'ET' ('PROG'
(A0O2 '+' (A1O2 * T) '+' (A2O2 * T2) '+' (A3O2 * T3) '+'
(A4O2 * T4))) ;
'FIN' BL1 ;
EVCVH2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCVH2 ;
EVCV1H2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCV1H2 ;
EVCVO2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCVO2 ;
EVCV1O2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCV1O2 ;
TAB1 = 'TABLE' ;
TAB1 . 'TITRE' = 'TABLE' ;
TAB1 . 1 ='TIRR ';
TAB1 . 'TITRE' . 1 = 'JANAF DATA';
TAB1 . 2 = 'MARQ CROI';
TAB1 . 'TITRE' . 2 = 'POLYN. REGR.';
'SI' GRAPH ;
'DESSIN' (EVCVH2 'ET' EVCV1H2) 'LEGE' TAB1 'TITRE'
'H2, polynome de 4-eme degre' ;
'DESSIN' (EVCVO2 'ET' EVCV1O2) 'LEGE' TAB1 'TITRE'
'O2, polynome de 4-eme degre' ;
'FINSI' ;
***************************
***** DOMAINE SPATIAL ****
***************************
A1 = 0.0D0 0.0D0;
A2 = 2.0D0 0.0D0;
A3 = 3.0D0 0.0D0;
A4 = 4.0D0 1.0D0;
A5 = 1.0D0 1.0D0;
A6 = 0.0D0 1.0D0;
A7 = 7.0D0 8.0D0 ;
MAI1 = ('MANUEL' 'POI1' A1) 'ET'
('MANUEL' 'POI1' A2) 'ET'
('MANUEL' 'POI1' A3) 'ET'
('MANUEL' 'POI1' A4) 'ET'
('MANUEL' 'POI1' A5) 'ET'
('MANUEL' 'POI1' A6) ;
*
MAI2 = ('MANUEL' 'POI1' A2) 'ET'
('MANUEL' 'POI1' A1) 'ET'
('MANUEL' 'POI1' A3) 'ET'
('MANUEL' 'POI1' A4) 'ET'
('MANUEL' 'POI1' A5) 'ET'
('MANUEL' 'POI1' A6) ;
*
MAI3 = ('MANUEL' 'POI1' A3) 'ET'
('MANUEL' 'POI1' A2) 'ET'
('MANUEL' 'POI1' A1) 'ET'
('MANUEL' 'POI1' A4) 'ET'
('MANUEL' 'POI1' A5) 'ET'
('MANUEL' 'POI1' A6) ;
*
MAI4 = ('MANUEL' 'POI1' A4) 'ET'
('MANUEL' 'POI1' A3) 'ET'
('MANUEL' 'POI1' A2) 'ET'
('MANUEL' 'POI1' A1) 'ET'
('MANUEL' 'POI1' A5) 'ET'
('MANUEL' 'POI1' A6) ;
*
**** Pression, densite, temperature et vitesse
*
RN = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAI1 ;
TN = 'BRUI' 'BLAN' 'UNIF' 5000 1000 MAI1 ;
YN = 'NOMC' 'H2 '
('BRUI' 'BLAN' 'UNIF' 0.3 0.1 MAI1)
'NATU' 'DISCRET' ;
UNCH = 'MANUEL' 'CHPO' MAI2 1 'SCAL' 1.0D0 'NATU' 'DISCRET' ;
YO2 = 'NOMC' 'O2 ' (UNCH '-' ('PSCAL' UNCH YN ('MOTS' 'SCAL')
(PGAZ . 'ESPEULE'))) 'NATU' 'DISCRET' ;
YNTOT = YN 'ET' YO2 ;
RYN = YN '*' RN (PGAZ . 'ESPEULE') ('MOTS' 'SCAL')
(PGAZ . 'ESPEULE') ;
LMOT1 = ('MOTS' (PGAZ . 'ESPNEULE')) 'ET' (PGAZ . 'ESPEULE') ;
RYNTOT = YNTOT * RN LMOT1 ('MOTS' 'SCAL' 'SCAL') LMOT1 ;
VNX = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAI1 ;
VNY = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAI3 ;
VN = ('NOMC' VNX 'UX' 'NATU' 'DISCRET' ) 'ET'
('NOMC' VNY 'UY' 'NATU' 'DISCRET' ) ;
GN = VN '*' RN ('MOTS' 'UY' 'UX') ('MOTS' 'SCAL' 'SCAL')
('MOTS' 'UY' 'UX') ;
GN1 = 'COPIER' GN ;
*
* Calcul de l'energie thermique et du gamma
*
* Ether = \sum_{i=1,4} Y_i \int_{0}^{T} cv_i(x) dx
* gamma = (\sum_{i=1,4} Y_i cp_i) / (\sum_{i=1,4} Y_i cv_i)
*
T2 = TN * TN ;
T3 = T2 * TN ;
T4 = T3 * TN ;
T5 = T4 * TN ;
ETERH2 = 'NOMC' ((A0H2 * TN) '+' ((A1H2 '/' 2.0) * T2) '+'
((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
((A4H2 '/' 5.0) * T5)) 'H2 ' 'NATU' 'DISCRET' ;
ETERO2 = 'NOMC' ((A0O2 * TN) '+' ((A1O2 '/' 2.0) * T2) '+'
((A2O2 '/' 3.0) * T3) '+' ((A3O2 '/' 4.0) * T4) '+'
((A4O2 '/' 5.0) * T5)) 'O2 ' 'NATU' 'DISCRET' ;
ETHERSIN = ETERH2 'ET' ETERO2 ;
RETHER = 'PSCAL' ETHERSIN RYNTOT LMOT1 LMOT1 ;
LMOT2 = 'MOTS' 'UX ' 'UY ';
RECIN = 'PSCAL' (0.5 '*' VN) GN LMOT2 LMOT2 ;
RET = RECIN '+' RETHER ;
CVH2 = 'NOMC' (A0H2 '+' (A1H2 * TN) '+'
(A2H2 * T2) '+' (A3H2 * T3) '+'
(A4H2 * T4)) 'H2 ' 'NATU' 'DISCRET' ;
CVO2 = 'NOMC' (A0O2 '+' (A1O2 * TN) '+'
(A2O2 * T2) '+' (A3O2 * T3) '+'
(A4O2 * T4)) 'O2 ' 'NATU' 'DISCRET' ;
CPH2 = CVH2 '+' (PGAZ . 'H2 ' . 'R') ;
CPO2 = CVO2 '+' (PGAZ . 'O2 ' . 'R') ;
CPSIN = CPH2 'ET' CPO2 ;
CPTOT = 'PSCAL' CPSIN YNTOT LMOT1 LMOT1 ;
CVSIN = CVH2 'ET' CVO2 ;
CVTOT = 'PSCAL' CVSIN YNTOT LMOT1 LMOT1 ;
GAMN = CPTOT '/' CVTOT ;
RTOT = CPTOT '-' CVTOT ;
PN = (RTOT '*' TN) '*' RN ;
*
**** PRIM
*
VITESSE PRES TEMPNEW YNNEW GAMNEW
= 'PRIM' 'PERFTEMP' PGAZ RN GN RET RYN ;
*
**** TEST
*
ERRGN = ('MAXIMUM' (GN '-' GN1) 'ABS') '/'
('MAXIMUM' GN 'ABS' ) ;
ERRV = ('MAXIMUM' (VN '-' VITESSE) 'ABS')
'/' ('MAXIMUM' VN 'ABS' ) ;
ERRP = 'MAXIMUM' ((PRES '-' PN) '/' PN) 'ABS';
ERRT = 'MAXIMUM' ((TEMPNEW '-' TN) '/' TN) 'ABS';
ERRG = 'MAXIMUM' ((GAMN '-' GAMNEW) '/' GAMN) 'ABS';
ERRY = ('MAXIMUM' (YN '-' YNNEW) 'ABS')
'/' ('MAXIMUM' YN) ;
'SI' (('MAXIMUM' ('PROG' ERRGN ERRV ERRP ERRT ERRG ERRY )) '>'
ERRMAX);
'MESSAGE' ('CHAINE' 'Erreur maximum');
'LISTE' ('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG ERRY
));
'ERREUR' 5;
'FINSI' ;
**** Avec une temperature de premier tentative
VITESSE PRES TEMPNEW2 YNNEW GAMNEW
= 'PRIM' 'PERFTEMP' PGAZ RN GN RET RYN TEMPNEW ;
ERRT = 'MAXIMUM' ((TEMPNEW '-' TEMPNEW2) '/' TEMPNEW) 'ABS';
'SI' (ERRT > ERRMAX);
'MESSAGE' ('CHAINE' 'Erreur maximum' ERRT) ;
'ERREUR' 5;
'FINSI' ;
**** On teste aussi la difference entre l'energie initiale et celle la.
T2 = TEMPNEW * TEMPNEW ;
T3 = T2 * TEMPNEW ;
T4 = T3 * TEMPNEW ;
T5 = T4 * TEMPNEW ;
ETERH2 = 'NOMC' ((A0H2 * TEMPNEW) '+' ((A1H2 '/' 2.0) * T2) '+'
((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
((A4H2 '/' 5.0) * T5)) 'H2 ' 'NATU' 'DISCRET' ;
ETERO2 = 'NOMC' ((A0O2 * TEMPNEW) '+' ((A1O2 '/' 2.0) * T2) '+'
((A2O2 '/' 3.0) * T3) '+' ((A3O2 '/' 4.0) * T4) '+'
((A4O2 '/' 5.0) * T5)) 'O2 ' 'NATU' 'DISCRET' ;
ETHERSIN = ETERH2 'ET' ETERO2 ;
RETHER_N = 'PSCAL' ETHERSIN RYNTOT LMOT1 LMOT1 ;
ERRET = 'MAXIMUM' ((RETHER '-' RETHER_N) '/' RETHER) 'ABS' ;
'SI' (ERRET > ERRMAX);
'MESSAGE' ('CHAINE' 'Erreur maximum' ERRET) ;
'ERREUR' 5;
'FINSI' ;
*
**** Scalaires passifs
*
PGAZ . 'SCALPASS' = 'MOTS' 'REPS' 'RTAU' ;
SN = 'MANUEL' 'CHPO' MAI1 2 'RTAU' ('PROG' 12.5 11.5 10.5
9.5 12.9 11.0) 'REPS' ('PROG' 11.5 10.5
9.5 12.9 11.0 12.5) ;
RSN = SN '*' RN ('MOTS' 'RTAU' 'REPS') ('MOTS' 'SCAL' 'SCAL')
('MOTS' 'RTAU' 'REPS') ;
RSN1 = 'COPIER' RSN ;
VITESSE PRES TEMPNEW2 YNNEW SNNEW GAMNEW
= 'PRIM' 'PERFTEMP' PGAZ RN GN RET RYN RSN TEMPNEW ;
ERRSN = ('MAXIMUM' (SN '-' SNNEW) 'ABS') '/'
('MAXIMUM' SN 'ABS' ) ;
ERRRSN = ('MAXIMUM' (RSN '-' RSN1) 'ABS') '/'
('MAXIMUM' RSN 'ABS' ) ;
ERRGN = ('MAXIMUM' (GN '-' GN1) 'ABS') '/'
('MAXIMUM' GN 'ABS' ) ;
ERRV = ('MAXIMUM' (VN '-' VITESSE) 'ABS')
'/' ('MAXIMUM' VN 'ABS' ) ;
ERRP = 'MAXIMUM' ((PRES '-' PN) '/' PN) 'ABS';
ERRT = 'MAXIMUM' ((TEMPNEW '-' TN) '/' TN) 'ABS';
ERRG = 'MAXIMUM' ((GAMN '-' GAMNEW) '/' GAMN) 'ABS';
ERRY = ('MAXIMUM' (YN '-' YNNEW) 'ABS')
'/' ('MAXIMUM' YN) ;
'SI' (('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG ERRY ERRSN ERRRSN )) '>'
ERRMAX);
'MESSAGE' ('CHAINE' 'Erreur maximum');
'LISTE' ('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG ERRY
));
'ERREUR' 5;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer