* fichier : prim_ther_mono_3D.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ *********************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la **** **** solution des **** **** Equations d'Euler pour un gaz parfait **** **** OPERATEUR PRIM **** **** Gaz monoespece "thermally perfect" **** **** Cas 3D **** **** **** **** A. BECCANTINI DRN/DMT/SEMT/TTMF JUILLET 1998 **** *********************************************************** * 'OPTION' 'DIME' 3 ; 'OPTI' 'ELEM' 'CUB8' ; 'OPTION' 'ECHO' 0 ; 'OPTION' 'TRAC' 'X'; * * **** GRAPH * * GRAPH = FAUX ; * GRAPH = VRAI ; * ERRMAX = 1.0D-6 ; 'MESSAGE' ; 'MESSAGE' ('cv = polynom de 4-eme degree'); 'MESSAGE' ; * *** On considere une melange H2, O2, H2O, N2 * * 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 ) ; *********************** **** LA TABLE PGAZ **** *********************** PGAZ = 'TABLE' ; PGAZ . 'NORD' = 4 ; * **** La seul espece * PGAZ . 'ESPNEULE' = 'H2 '; * PGAZ . 'H2 ' = 'TABLE' ; * **** R (J/Kg/K) * PGAZ . 'H2 ' . 'R' = 4130.0 ; * **** Regressions polynomials * PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836 -2.37281455E-07 1.84701105E-11 ; * **** "Enthalpies" (ou energies) de formations a OK (J/Kg) * Note: ce sont des entites numeriques PGAZ . 'H2 ' . 'H0K' = -4.195D6 ; * *** 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 ; LCV1H2 = 'PROG' ; LCV1H2O = 'PROG' ; LCV1N2 = '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))) ; 'FIN' BL1 ; EVCVH2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCVH2 ; EVCV1H2 = 'EVOL' 'MANU' 'T(K)' LTEMP 'CV (J/Kg/K)' LCV1H2 ; 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' ; 'FINSI' ; *************************** ***** DOMAINE SPATIAL **** *************************** ************ * MAILLAGE * ************ NX = 10 ; NY = 2 ; NZ = 2 ; L = 1.0D0 ; DX = L '/' NX '/' 2.0D0 ; H = NY '*' DX ; P = NZ '*' DX ; xD = 0.5D0 '*' L ; xG = -1.0D0 '*' xD ; yH = 0.5D0 '*' H ; yB = -1.0D0 '*' yH ; zV = 0.5D0 '*' P ; zR = -1.0D0 '*' zV ; A1 = xG yB zR ; A2 = 0.0D0 yB zR ; A7 = xG yB zV ; A8 = 0.0D0 yB zV ; BAS1 = A1 'DROI' NX A2 ; BAS3 = A7 'DROI' NX A8 ; BAS5 = A1 'DROI' NZ A7 ; BAS6 = A2 'DROI' NZ A8 ; S11 = 'DALL' BAS3 BAS6 BAS1 BAS5 ; DOM1 = S11 'VOLU' NZ 'TRANS' (0.0 H 0.0) ; $DOM1 = 'DOMA' DOM1 ; * **** Pression, densite, temperature et vitesse * MAIL = $DOM1 . 'CENTRE' ; RN = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAIL ; TN = 'BRUI' 'BLAN' 'UNIF' 3900 2000 MAIL ; VNX = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAIL ; VNY = 'BRUI' 'BLAN' 'UNIF' 2.0 0.1 MAIL ; VNZ = 'BRUI' 'BLAN' 'UNIF' 3.0 0.1 MAIL ; VN = ('NOMC' VNX 'UX' 'NATU' 'DISCRET' ) 'ET' ('NOMC' VNY 'UY' 'NATU' 'DISCRET' ) 'ET' ('NOMC' VNZ 'UZ' 'NATU' 'DISCRET' ) ; GN = ('NOMC' (VNX * RN) 'UX' 'NATU' 'DISCRET' ) 'ET' ('NOMC' (VNY * RN) 'UY' 'NATU' 'DISCRET' ) 'ET' ('NOMC' (VNZ * RN) 'UZ' 'NATU' 'DISCRET' ) ; * * 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 ; ETHER = ((A0H2 * TN) '+' ((A1H2 '/' 2.0) * T2) '+' ((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+' ((A4H2 '/' 5.0) * T5)) ; LISCOM2 = 'MOTS' 'UX ' 'UY ' 'UZ '; RECIN = 'PSCAL' (0.5 '*' VN) GN LISCOM2 LISCOM2 ; RET = RECIN '+' (ETHER * RN) ; CVTOT = (A0H2 '+' (A1H2 * TN) '+' (A2H2 * T2) '+' (A3H2 * T3) '+' (A4H2 * T4)) ; CPTOT = CVTOT '+' (PGAZ . 'H2 ' . 'R') ; GAMN = CPTOT '/' CVTOT ; RTOT = CPTOT '-' CVTOT ; PN = (RTOT '*' TN) '*' RN ; * **** PRIM * VITESSE PRES TEMPNEW GAMNEW = 'PRIM' 'PERFTEMP' PGAZ RN GN RET ; * **** TEST * ERRV = ('MAXIMUM' (VN '-' VITESSE) 'ABS') '/' ('MAXIMUM' VN) ; ERRP = 'MAXIMUM' ((PRES '-' PN) '/' PN) 'ABS'; ERRT = 'MAXIMUM' ((TEMPNEW '-' TN) '/' TN) 'ABS'; ERRG = 'MAXIMUM' ((GAMN '-' GAMNEW) '/' GAMN) 'ABS'; 'SI' (('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG ) '>' ERRMAX)); 'MESSAGE' ('CHAINE' 'Erreur maximum'); 'LISTE' ('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG )); 'ERREUR' 5; 'FINSI' ; VITESSE PRES TEMPNEW2 GAMNEW = 'PRIM' 'PERFTEMP' PGAZ RN GN RET 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 ; ETHER_N = ((A0H2 * TEMPNEW) '+' ((A1H2 '/' 2.0) * T2) '+' ((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+' ((A4H2 '/' 5.0) * T5)) ; ERRET = 'MAXIMUM' ((ETHER '-' ETHER_N) '/' ETHER) 'ABS' ; 'SI' (ERRET > ERRMAX); 'MESSAGE' ('CHAINE' 'Erreur maximum' ERRET) ; 'ERREUR' 5; 'FINSI' ; *) Degre de polynome = 0 ERRMAX = 1.0D-10 ; 'MESSAGE' ; 'MESSAGE' ('cv = polynom constant'); 'MESSAGE' ; *********************** **** LA TABLE PGAZ **** *********************** PGAZ = 'TABLE' ; * **** Degre des polynoms * PGAZ . 'NORD' = 0 ; * **** Especes qui sont dans les equations d'Euler * * **** Espece qui n'y est pas * PGAZ . 'ESPNEULE' = 'H2 '; * PGAZ . 'H2 ' = 'TABLE' ; * **** R (J/Kg/K) * PGAZ . 'H2 ' . 'R' = 4130.0 ; * **** Regressions polynomials * * * Cas particulier: gaz "calorically perfect" * PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 ; * **** "Enthalpies" (ou energies) de formations a OK (J/Kg) * PGAZ . 'H2 ' . 'H0K' = -4.195D6 ; * *** Fin PGAZ * A0H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 1 ; * **** Pression, densite, temperature et vitesse * RN = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAIL ; TN = 'BRUI' 'BLAN' 'UNIF' 3900 200 MAIL ; YH2 = 'BRUI' 'BLAN' 'UNIF' 0.3 0.1 MAIL ; YO2 = 'BRUI' 'BLAN' 'UNIF' 0.2 0.1 MAIL ; YH2O = 'BRUI' 'BLAN' 'UNIF' 0.1 0.05 MAIL ; YN2 = 1.0 '-' (YH2 '+' YO2 '+' YH2O) ; RYH2 = RN * YH2 ; RYH2O = RN * YH2O ; RYO2 = RN * YO2 ; RYN2 = RN * YN2 ; YH2 = 'NOMC' YH2 'H2 ' 'NATU' 'DISCRET' ; YH2O = 'NOMC' YH2O 'H2O ' 'NATU' 'DISCRET' ; YO2 = 'NOMC' YO2 'O2 ' 'NATU' 'DISCRET' ; YN2 = 'NOMC' YN2 'N2 ' 'NATU' 'DISCRET' ; YN = YH2 'ET' YO2 'ET' YH2O ; RYH2 = 'NOMC' RYH2 'H2 ' 'NATU' 'DISCRET' ; RYH2O = 'NOMC' RYH2O 'H2O ' 'NATU' 'DISCRET' ; RYO2 = 'NOMC' RYO2 'O2 ' 'NATU' 'DISCRET' ; RYN2 = 'NOMC' RYN2 'N2 ' 'NATU' 'DISCRET' ; VNX = 'BRUI' 'BLAN' 'UNIF' 1.0 0.1 MAIL ; VNY = 'BRUI' 'BLAN' 'UNIF' 2.0 0.1 MAIL ; VNZ = 'BRUI' 'BLAN' 'UNIF' 3.0 0.1 MAIL ; VN = ('NOMC' VNX 'UX' 'NATU' 'DISCRET' ) 'ET' ('NOMC' VNY 'UY' 'NATU' 'DISCRET' ) 'ET' ('NOMC' VNZ 'UZ' 'NATU' 'DISCRET' ) ; GN = ('NOMC' (VNX * RN) 'UX' 'NATU' 'DISCRET' ) 'ET' ('NOMC' (VNY * RN) 'UY' 'NATU' 'DISCRET' ) 'ET' ('NOMC' (VNZ * RN) 'UZ' 'NATU' 'DISCRET' ) ; RYN = RYH2 'ET' RYH2O 'ET' RYO2 ; * * Calcul de l'energie thermique et du gamma * RETHER = (A0H2 '*' RN '*' TN) ; LISCOM2 = 'MOTS' 'UX ' 'UY ' 'UZ '; RECIN = 'PSCAL' (0.5 '*' VN) GN LISCOM2 LISCOM2 ; RET = RECIN '+' RETHER ; CHUN = 'MANUEL' 'CHPO' MAIL 1 'SCAL' 1.0 ; CV = 'NOMC' (A0H2 '*' CHUN) 'SCAL' 'NATU' 'DISCRET' ; CP = CV '+' (PGAZ . 'H2 ' . 'R') ; GAMN = CP '/' CV ; RTOT = CP '-' CV ; PN = (RTOT '*' TN) '*' RN ; * **** PRIM * VITESSE PRES TEMPNEW GAMNEW = 'PRIM' 'PERFTEMP' PGAZ RN GN RET ; * **** TEST * ERRV = ('MAXIMUM' (VN '-' VITESSE) 'ABS') '/' ('MAXIMUM' VN) ; ERRP = 'MAXIMUM' ((PRES '-' PN) '/' PN) 'ABS'; ERRT = 'MAXIMUM' ((TEMPNEW '-' TN) '/' TN) 'ABS'; ERRG = 'MAXIMUM' ((GAMN '-' GAMNEW) '/' GAMN) 'ABS'; 'SI' (('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG ) '>' ERRMAX)); 'MESSAGE' ('CHAINE' 'Erreur maximum'); 'LISTE' ('MAXIMUM' ('PROG' ERRV ERRP ERRT ERRG )); 'ERREUR' 5; 'FINSI' ; * Avec une temperature d'essai VITESSE PRES TEMPNEW2 GAMNEW = 'PRIM' 'PERFTEMP' PGAZ RN GN RET 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. RETHER_N = (A0H2 '*' RN '*' TEMPNEW) ; ERRET = 'MAXIMUM' ((RETHER '-' RETHER_N) '/' RETHER) 'ABS' ; 'SI' (ERRET > ERRMAX); 'MESSAGE' ('CHAINE' 'Erreur maximum' ERRET) ; 'ERREUR' 5; 'FINSI' ; 'FIN' ;