* fichier : carre_therper.dgibi
******************************************************************
* CALCUL DU QUARRE A CHOC *
* Modele: gaz multi-especes "thermally perfect" *
* Gaz multi-especes "thermally perfect" *
* Controle de la symetrie *
* *
* A. BECCANTINI, SEMT/LTMF, MARS 1999 *
******************************************************************
'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'ISOV' 'SULI' 'ECHO' 1
'TRAC' 'X';
GRAPH = FAUX ;
* GRAPH = VRAI ;
*
*** Methodes possibles :
*
* 'VLH', 'SS'
*
METO = 'VLH' ;
* SAFFAC = safety facton on the CFL
SAFFAC = 1.0 ;
* LOGSO = second order reconstruction
LOGSO = VRAI ;
TFINAL = 10. ;
NITER = 40 ;
ZERO = 1.0D-12 ;
************
* MAILLAGE *
************
PCEN = 0.5 0.5 ;
A1 = 0.75 0.5 ;
A2 = 1.0 0.5 ;
A3 = 1.0 0.75 ;
A4 = 1.0 1.0 ;
A5 = 0.75 1.0 ;
A6 = 0.5 1.0 ;
A7 = 0.5 0.75 ;
MAIL1 = MANU 'QUA4' A1 A3 A5 A7 ;
MAIL1 = MAIL1 ET (MANU 'TRI3' PCEN A1 A7 );
MAIL1 = MAIL1 ET (MANU 'TRI3' A1 A2 A3 );
MAIL1 = MAIL1 ET (MANU 'TRI3' A3 A4 A5 );
MAIL1 = MAIL1 ET (MANU 'TRI3' A5 A6 A7 );
MAIL2 = MAIL1;
REPE BL1 3 ;
MAIL2 = MAIL2 TOUR 90 PCEN ;
MAIL1 = MAIL1 ET MAIL2 ;
FIN BL1 ;
ELIM MAIL1 0.001 ;
VECX = 1.0 0.0 ;
VECY = 0.0 1.0 ;
NELEM = 1 ;
NELEMTOT = (2 * NELEM) + 1 ;
MAILTOT = MAIL1 'COUL' 'BLEU' ;
IX = (-1 * NELEM) - 1 ;
'REPE' BL1 NELEMTOT ;
IX = IX + 1;
IY = (-1 * NELEM) - 1 ;
'REPE' BL2 NELEMTOT ;
IY = IY + 1;
DVECX = IX '*' VECX ;
DVECY = IY '*' VECY ;
DVEC = 'PLUS' DVECX DVECY ;
SI ((IX NEG 0) OU (IY NEG 0)) ;
MAIL2 = MAIL1 PLUS DVEC ;
MAILTOT = MAILTOT 'ET' MAIL2 ;
FINSI ;
'FIN' BL2 ;
'FIN' BL1 ;
'ELIM' MAILTOT 0.01 ;
'ELIM' 0.01 MAIL1 MAILTOT ;
MDNS = 'EULER' ;
$DOMTOT = 'MODE' MAILTOT MDNS ;
$DOM1 = 'MODE' MAIL1 MDNS ;
TAB1 = 'DOMA' $DOMTOT 'VF' ;
TAB2 = 'DOMA' $DOM1 'VF' ;
MDOMTOT = TAB1 . 'QUAF' ;
MDOM1 = TAB2 . 'QUAF' ;
'ELIM' (MDOMTOT 'ET' MDOM1) 1.E-5 ;
'SI' GRAPH ;
'TRACER' (TAB1 . 'MAILLAGE') 'TITRE' ('CHAINE' 'Maillage ');
'FINSI' ;
***********************
**** LA TABLE PGAZ ****
***********************
PGAZ = 'TABLE' ;
*
**** Ordre des polynoms
*
PGAZ . 'NORD' = 4 ;
*
**** Especes qui sont dans les equations d'Euler
*
PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'O2 ' 'H2O ' ;
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'N2 ';
*
PGAZ . 'H2 ' = 'TABLE' ;
PGAZ . 'H2O ' = 'TABLE' ;
PGAZ . 'N2 ' = 'TABLE' ;
PGAZ . 'O2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'H2 ' . 'R' = 4130.0 ;
PGAZ . 'H2O ' . 'R' = 461.4 ;
PGAZ . 'N2 ' . 'R' = 296.8 ;
PGAZ . 'O2 ' . 'R' = 259.8 ;
*
**** Regressions polynomials
*
PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAZ . 'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
-1.82753232E-08 2.44485692E-12 ;
PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
8.78233606E-09 -3.05514485E-13 ;
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 . 'H2O ' . 'H0K' = -1.395D7 ;
PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
*
*** Fin PGAZ
*
*
**** ETAT CENTRE
*
rog = 1.0 ;
ung = 0.0 ;
utg = 0.0 ;
retg = .4291145555695540D+04 ;
yh2g = 0.01 ;
yo2g = 0.2 ;
yh2og = 0.3 ;
rouxg = ung '*' rog ;
rouyg = utg '*' rog ;
ROC1 = KCHT $DOM1 'SCAL' 'CENTRE' rog ;
ROVC1 = KCHT $DOM1 'VECT' 'CENTRE' (rouxg rouyg );
RETC1 = KCHT $DOM1 'SCAL' 'CENTRE' retg ;
RYH2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2 ' (rog*yh2g) ;
RYO2C1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'O2 ' (rog*yo2g) ;
RYH2OC1 = KCHT $DOM1 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rog*yh2og) ;
*
*** Le reste
*
rod = 1.250D-1 ;
und = 0.0D0 ;
utd = 0.0D0 ;
retd = .3598345082089522D+01 ;
yh2d = 0.2 ;
yo2d = 0.4 ;
yh2od = 0.1 ;
rouxd = und '*' rod ;
rouyd = utd '*' rod ;
ROC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' rod ;
ROVC2 = KCHT $DOMTOT 'VECT' 'CENTRE' (rouxd rouyd );
RETC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' retd ;
RYH2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2 ' (rod * yh2d) ;
RYO2C2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2 ' (rod * yo2d) ;
RYH2OC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' (rod * yh2od) ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' ROC2 ROC1 ;
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' ROVC2 ROVC1 ;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' RETC2 RETC1 ;
RYH2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2 ' RYH2C2 RYH2C1 ;
RYO2 = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'O2 ' RYO2C2 RYO2C1 ;
RYH2O = KCHT $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'H2O ' RYH2OC2 RYH2OC1 ;
RYN = RYH2 'ET' RYO2 'ET' RYH2O ;
********************************************************
*** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM ***
********************************************************
MOD1 = 'MODELISER' (TAB1 . 'MAILLAGE') 'THERMIQUE' ;
*
*** GRAPHIQUE DES C.I.
*
'SI' GRAPH ;
* 'SI' FAUX ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ;
CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ;
CHM_RYN = 'KCHA' $DOMTOT 'CHAM' RYN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RN at t=' 0.0);
TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'GN at t=' 0.0);
TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' 0.0);
TRAC CHM_RYN MOD1 'TITR' ('CHAINE' 'RYN at t=' 0.0);
'FINSI' ;
**********************************************************
RN0 = 'COPIER' RN ;
GN0 = 'COPIER' GN ;
RETN0 = 'COPIER' RETN ;
RYN0 = 'COPIER' RYN ;
*
**** Parameter for the time loop
*
* Names of the residuum components
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RETN' 'H2' 'O2' 'H2O' ;
NOMRN = 'MOTS' 'SCAL' ;
NOMVN = 'MOTS' 'UX' 'UY';
NOMY = 'MOTS' 'H2' 'O2' 'H2O' ;
*
**** Geometrical coefficient to compute gradients
*
GRADR CACCA COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
NOMRN RN ;
GRADV CACCA COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
NOMVN GN ;
TPS = 0.0 ;
*
**** Temporal loop
*
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
'MESSAGE' ('CHAINE' 'SAFFAC = ' SAFFAC) ;
'MESSAGE' ;
VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
TNM1 = 'COPIER' TN ;
'TEMPS' 'ZERO' ;
'REPETER' BL1 NITER ;
*
**** Primitive variables
*
VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN
TNM1 ;
TNM1 = 'COPIER' TN ;
'SI' LOGSO ;
GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
NOMRN RN 'GRADGEO' COEFSCAL ;
GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
NOMRN PN 'GRADGEO' COEFSCAL ;
GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
NOMVN VN 'GRADGEO' COEFVECT ;
GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
NOMY YN 'GRADGEO' COEFSCAL ;
ROF VITF PF YF = 'PRET' 'PERFTEMP' 2 1
$DOMTOT PGAZ
RN GRADR ALR
VN GRADV ALV
PN GRADP ALP
YN GRADY ALY ;
'SINON' ;
ROF VITF PF YF = 'PRET' 'PERFTEMP' 1 1
$DOMTOT PGAZ
RN
VN
PN
YN ;
'FINSI' ;
RESIDU DELTAT = 'KONV' 'VF' 'PERFTEMP' 'RESI' METO
$DOMTOT PGAZ LISTINCO ROF VITF PF YF ;
DT_CON = SAFFAC '*' DELTAT ;
*
**** The time step linked to tps
*
DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
*
**** Total time step
*
DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
*
**** Increment of the variables (convection)
*
RESIDU = DTMIN '*' RESIDU ;
DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY') ;
DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
DRYN = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
('MOTS' 'H2' 'O2' 'H2O') ;
TPS = TPS '+' DTMIN ;
RN = RN '+' DRN ;
GN = GN '+' DGN ;
RETN = RETN '+' DRETN ;
RYN = RYN '+' DRYN ;
'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ;
'FINSI' ;
'SI' (TPS > TFINAL) ;
'QUITTER' BL1 ;
'FINSI' ;
'FIN' BL1 ;
'TEMPS' ;
**********************************************************
*
**** Les variables primitives
*
VN PN TN YN GAMN = 'PRIM' 'PERFTEMP' PGAZ
RN GN RETN RYN ;
CN = (GAMN * (PN / RN)) '**' 0.5 ;
*
*** GRAPHIQUE DES SOLUTIONS
*
'SI' GRAPH ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ;
CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ;
CHM_RYN = 'KCHA' $DOMTOT 'CHAM' RYN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' TFIN);
TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'GN at t=' TFIN);
TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'ROET at t=' TFIN);
TRAC CHM_RYN MOD1 'TITR' ('CHAINE' 'RYN at t=' TFIN);
'FINSI' ;
* Tests: invariance par rotation
RN1 = COPIER RN ;
GN1 = COPI GN ;
RETN1 = COPIER RETN ;
RYN1 = COPIER RYN ;
MAIL0 = EXTR RN 'MAILLAGE' ;
REPETER BL1 3 ;
RN1 = RN1 TOUR 90 PCEN ;
GN1 = GN1 TOUR 90 PCEN ;
RETN1 = RETN1 TOUR 90 PCEN ;
RYN1 = RYN1 TOUR 90 PCEN ;
MAIL1 = EXTR RN1 'MAILLAGE' ;
ELIM 0.001 MAIL0 MAIL1 ;
ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ;
SI (ERR1 > 1.0D-10);
ERRE 5;
FINSI ;
MAIL1 = EXTR GN1 'MAILLAGE' ;
ELIM 0.001 MAIL0 MAIL1 ;
ERR1 = MAXI ((GN1 - GN) / (RN * CN)) 'ABS' ;
SI (ERR1 > 1.0D-10);
ERRE 5;
FINSI ;
MAIL1 = EXTR RETN1 'MAILLAGE' ;
ELIM 0.001 MAIL0 MAIL1 ;
ERR1 = MAXI ((RETN1 - RETN) / RETN) 'ABS' ;
SI (ERR1 > 1.0D-10);
ERRE 5;
FINSI ;
MAIL1 = EXTR RYN1 'MAILLAGE' ;
ELIM 0.001 MAIL0 MAIL1 ;
ERR1 = MAXI ((RYN1 - RYN) / RN) 'ABS' ;
SI (ERR1 > 1.0D-10);
ERRE 5;
FINSI ;
FIN BL1 ;
* Symetrie
MAIL1 RN1 = MAIL0 RN SYME 'DROI' (0 0.5) (1.0 0.5) ;
ELIM 0.001 MAIL0 MAIL1 ;
ERR1 = MAXI ((RN1 - RN) / RN) 'ABS' ;
SI (ERR1 > 1.0D-10);
ERRE 5;
FINSI ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer