* fichier : primtest3_3D.dgibi *********************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la **** **** solution des **** **** Equations d'Euler pour un gaz parfait **** **** OPERATEUR PRIM **** **** Cas: gaz multiespece "calorically perfect" **** **** Cas 3D **** **** **** **** A. BECCANTINI DRN/DMT/SEMT/TTMF MARS 1998 **** *********************************************************** 'OPTION' 'DIME' 3 ; 'OPTION' 'ELEM' CUB8 ; 'OPTION' 'ECHO' 0; 'OPTION' 'TRAC' 'X' ; * *** GRAPHIQUES * * GRAPH = VRAI ; GRAPH = FAUX ; ************ * 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) ; DOM2 = S11 'VOLU' NZ 'TRANS' (0.0 (-1.0*H) 0.0) ; DOMTOT = DOM1 ET DOM2; 'ELIMINATION' DOMTOT 1D-6; $DOMTOT = 'DOMA' DOMTOT ; $DOM1 = 'DOMA' DOM1 'INCL' $DOMTOT ; $DOM2 = 'DOMA' DOM2 'INCL' $DOMTOT ; 'SI' GRAPH; 'TRACER' ($DOMTOT . 'MAILLAGE' ) 'TITRE' 'Maillage'; 'FINSI' ; ************************ **** MODELE DU GAZ **** ************************ NESP = 4; * *** GAS: H_2, O_2, H_2O, N_2 * * CP, CV en J/Kg/K @ T = 3000 * PGAS = 'TABLE' ; PGAS . 'CP' = 'TABLE' ; PGAS . 'CP' . 'H2 ' = .18729066D+05 ; PGAS . 'CP' . 'O2 ' = .11886820D+04 ; PGAS . 'CP' . 'H2O ' = .31209047D+04 ; PGAS . 'CP' . 'N2 ' = .12993995D+04 ; PGAS . 'CV' = 'TABLE' ; PGAS . 'CV' . 'H2 ' = .14571861D+05 ; PGAS . 'CV' . 'O2 ' = .92885670D+03 ; PGAS . 'CV' . 'H2O ' = .26589930D+04 ; PGAS . 'CV' . 'N2 ' = .10024563D+04; * **** Especes qui sont dans les equations d'Euler * PGAS . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ; * **** Espece qui n'y est pas * PGAS . 'ESPNEULE' = 'N2 '; * *********************** **** Les CHPOINTs **** *********************** RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 1.9D0; RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 2.9D0; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2); GN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (4.9D0 7.9D0 12.3); GN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (0.9D0 11.7D0 1.54); GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN1 'ET' GN2); PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 1.19D2; PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 2.92D3; PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2); * **** Les fractions massiques YN 'ET' les masses RYN * * N.B. YN n'as pas le meme supporte geometrique * que RN, PN; en general la numerotation * des noeuds est diferente, mais ce n'est * pas grave: 'PRIM' change les masse des * especes RYN dans la subroutine QUEPOI!!! * UNCH = 'KCHT' $DOM1 'SCAL' 'CENTRE' 1.0D0 ; Y1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 0.1 ; RY1 = 'NOMC' 'H2' (RN1 * Y1) 'NATU' 'DISCRET'; YH2 = 'NOMC' 'H2' Y1 'NATU' 'DISCRET'; Y2 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 0.2 ; RY2 = 'NOMC' 'O2' (RN1 * Y2) 'NATU' 'DISCRET'; YO2 = 'NOMC' 'O2' Y2 'NATU' 'DISCRET'; Y3 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 0.3 ; RY3 = 'NOMC' 'H2O' (RN1 * Y3) 'NATU' 'DISCRET'; YH2O = 'NOMC' 'H2O' Y3 'NATU' 'DISCRET'; YTOT = Y1 '+' Y2 '+' Y3 ; Y4 = 'KCHT' $DOM1 'SCAL' 'CENTRE' (UNCH '-' YTOT); YN2 = 'NOMC' 'N2' Y4 'NATU' 'DISCRET'; YND1 = YH2 'ET' YO2 'ET' YH2O 'ET' YN2; RYND1 = RY2 'ET' RY3 'ET' RY1; UNCH = 'KCHT' $DOM2 'SCAL' 'CENTRE' 1.0D0 ; Y1 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 0.1 ; RY1 = 'NOMC' 'H2' (RN2 * Y1) 'NATU' 'DISCRET'; YH2 = 'NOMC' 'H2' Y1 'NATU' 'DISCRET'; Y2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 0.2 ; RY2 = 'NOMC' 'O2' (RN2 * Y2) 'NATU' 'DISCRET'; YO2 = 'NOMC' 'O2' Y2 'NATU' 'DISCRET'; Y3 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 0.3 ; RY3 = 'NOMC' 'H2O' (RN2 * Y3) 'NATU' 'DISCRET'; YH2O = 'NOMC' 'H2O' Y3 'NATU' 'DISCRET'; YTOT = Y1 '+' Y2 '+' Y3 ; Y4 = 'KCHT' $DOM2 'SCAL' 'CENTRE' (UNCH '-' YTOT); YN2 = 'NOMC' 'N2' Y4 'NATU' 'DISCRET'; YND2 = YH2 'ET' YO2 'ET' YH2O 'ET' YN2; RYND2 = RY3 'ET' RY2 'ET' RY1; YN = YND1 'ET' YND2; RYN = RYND1 'ET' RYND2 ; * **** Control de la Table PGAS * * **** La vitesse * VITX = ('EXCO' 'UX ' GN) '/' RN; VITY = ('EXCO' 'UY ' GN) '/' RN; VITZ = ('EXCO' 'UZ ' GN) '/' RN; * **** Les GAMMAs * NOMCOM = 'EXTRAIRE' YN 'COMP'; CPTOT = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 0.0D0; CVTOT = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 0.0D0; 'REPETER' BL1 ('DIME' NOMCOM); NOMCEL = 'EXTRAIRE' NOMCOM &BL1; YCEL = 'EXCO' NOMCEL YN ; CPCEL = (PGAS . 'CP' . NOMCEL); CVCEL = (PGAS . 'CV' . NOMCEL); CPTOT = CPTOT '+' (YCEL * CPCEL); CVTOT = CVTOT '+' (YCEL * CVCEL); 'FIN' BL1; GAMMA = CPTOT '/' CVTOT; RGAS = (CPTOT '-' CVTOT); * *** L'energie totale (ROET) et la temperature * GM1 = GAMMA '-' ('KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.0D0); ETHER = PN '/' GM1; lcg=mots 'UX' 'UY' 'UZ'; ECIN = (0.5D0 * (GN lcg 'PSCA' GN lcg)) '/' RN; EN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (ECIN '+' ETHER); TEMPE = (PN '/' RN) '/' RGAS; ************************ **** L'operateur ***** ************************ VITESSE PRES TEMPNEW YNNEW GAMNEW = 'PRIM' 'PERFMULT' PGAS RN GN EN RYN ; * **** L'erreur (pas d'integral pour semplifier le calcul); * VIT1X = 'EXCO' 'UX ' VITESSE ; VIT1Y = 'EXCO' 'UY ' VITESSE ; VIT1Z = 'EXCO' 'UZ ' VITESSE ; ERRVX = 'MAXIMUM' (VIT1X '-' VITX) 'ABS' ; ERRVY = 'MAXIMUM' (VIT1Y '-' VITY) 'ABS' ; ERRVZ = 'MAXIMUM' (VIT1Z '-' VITZ) 'ABS'; VCELL = ('MAXIMUM' VITX 'ABS') '+' ('MAXIMUM' VITY 'ABS' )'+' ('MAXIMUM' VITZ 'ABS' ); 'SI' (VCELL > 0); ERRVX = ERRVX '/' VCELL; ERRVY = ERRVY '/' VCELL; ERRVZ = ERRVZ '/' VCELL; 'FINSI' ; ERRP = ('MAXIMUM' (PRES '-' PN) 'ABS') '/' ('MAXIMUM' PN); ERRT = ('MAXIMUM' (TEMPNEW '-' TEMPE) 'ABS') '/' ('MAXIMUM' TEMPE) ; ERRG = ('MAXIMUM' (GAMMA '-' GAMNEW) 'ABS') '/' ('MAXIMUM' GAMMA); LCOMP = 'EXTRAIRE' YNNEW 'COMP' ; NCOMP = 'DIME' LCOMP ; ERRY = 0.0D0 ; 'REPETER' BL1 NCOMP ; NOMCEL = 'EXTRAIRE' LCOMP &BL1 ; YN1 = 'EXCO' NOMCEL YN 'SCAL' ; YN2 = 'EXCO' NOMCEL YNNEW 'SCAL' ; ERRY = ERRY '+' ('MAXIMUM' (YN1 '-' YN2) 'ABS') ; 'FINSI' ; ERRY = ERRY '/' NCOMP; 'SI' (('MAXIMUM' ('PROG' ERRVX ERRVY ERRVZ ERRP ERRT ERRG ERRY ) '>' 1.0D-12)); 'MESSAGE' ('CHAINE' 'Erreur maximum'); 'LISTE' ('MAXIMUM' ('PROG' ERRVX ERRVY ERRVZ ERRP ERRT ERRG ERRY )); 'ERREUR' 5; 'FINSI' ; 'FIN' ;