* 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;
ECIN = (0.5D0 * ('KOPS' GN 'PSCA' GN )) '/' 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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer