Download carre3D_therper.dgibi
* fichier : carre3D_therper.dgibi
******************************************************************
* CALCUL DU CUBE A CHOC *
* Modele: gaz multi-especes "thermally perfect" *
* Gaz multi-especes "thermally perfect" *
* Controle de la symetrie *
* *
* A. BECCANTINI, SEMT/LTMF, MARS 1999 *
* Modif, 25/09/02, syntaxe de KONV changée *
******************************************************************
'OPTION' 'DIME' 3 'ELEM' 'CUB8' 'ISOV' 'SULI' 'ECHO' 1
'TRAC' 'PS' ;
GRAPH = FAUX ;
GRAPH = VRAI ;
*
*** Methodes possibles :
*
* 'VLH', 'SS'
*
METO = 'VLH' ;
* SAFFAC = safety facton on the CFL
SAFFAC = .5 ;
* LOGSO = second order reconstruction
LOGSO = VRAI ;
TFINAL = 10. ;
NITER = 40 ;
ZERO = 1.0D-12 ;
************
* MAILLAGE *
************
A1 = 0.0 0.0 0.0 ;
A2 = 0.0 0.0 5.0 ;
LAT12 = A1 'DROIT' 5 A2 ;
SUR1 = LAT12 'TRANSLATION' 5 (0.0 5.0 0.0) ;
DOMTOT = 'VOLUME' SUR1 'TRANSLATION' 5 (5.0 0.0 0.0) ;
MDNS = 'EULER' ;
$DOMTOT = 'MODE' DOMTOT MDNS ;
TABM = 'DOMA' $DOMTOT 'VF' ;
MDOMTOT = TABM . 'QUAF' ;
* Comparaison EULER '-' 'NS'
* Element au centre
PCEN = (TABM . 'CENTRE') 'POIN' 'PROC' (2.5 2.5 2.5);
ELCEN = 'MANUEL' 'POI1' PCEN ;
'SI' GRAPH ;
'TRACER' (TABM . '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 ;
retg = .4291145555695540D+04 ;
yh2g = 0.01 ;
yo2g = 0.2 ;
yh2og = 0.3 ;
rouxg = 0.0 ;
rouyg = 0.0 ;
rouzg = 0.0 ;
ROC1 = 'MANUEL' 'CHPO' ELCEN 1 'SCAL' rog ;
ROVC1 = 'MANUEL' 'CHPO' ELCEN 3 'UX'rouxg 'UY' rouyg
'UZ' rouzg ;
RETC1 = 'MANUEL' 'CHPO' ELCEN 1 'SCAL' retg ;
RYH2C1 = 'MANUEL' 'CHPO' ELCEN 1 'H2 ' (rog*yh2g) ;
RYO2C1 = 'MANUEL' 'CHPO' ELCEN 1 'O2 ' (rog*yo2g) ;
RYH2OC1 = 'MANUEL' 'CHPO' ELCEN 1 'H2O ' (rog*yh2og) ;
*
*** Le reste
*
rod = 1.250D-1 ;
retd = .3598345082089522D+01 ;
yh2d = 0.2 ;
yo2d = 0.4 ;
yh2od = 0.1 ;
rouxd = 0.0 ;
rouyd = 0.0 ;
rouzd = 0.0 ;
ROC2 = KCHT $DOMTOT 'SCAL' 'CENTRE' rod ;
ROVC2 = KCHT $DOMTOT 'VECT' 'CENTRE' (rouxd rouyd rouzd);
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' (TABM . '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 ;
* Tests: invariance par rotation de la CI
RN1 = COPIER RN0 ;
GN1 = COPI GN0 ;
RETN1 = COPIER RETN0 ;
RYN1 = COPIER RYN0 ;
MAIL0 = EXTR RN1 'MAILLAGE' ;
* Axe de rotation: PCEN = 2.5 2.5 2.5 '+' VECX ;
VECX = 1.0 0.0 0.0 ;
'REPETER' BL1 3 ;
RN1 = RN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
GN1 = GN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
RETN1 = RETN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
RYN1 = RYN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
MAIL1 = 'EXTRAIRE' RN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RN1 '-' RN0) / RN0) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' RETN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RETN1 '-' RETN0) / RETN0) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' RYN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RYN1 '-' RYN0) / RN0) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
'FIN' BL1 ;
*
**** Parameter for the time loop
*
* Names of the residuum components
LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RUZ' 'RETN' 'H2' 'O2' 'H2O' ;
NOMRN = 'MOTS' 'SCAL' ;
NOMVN = 'MOTS' 'UX' 'UY' 'UZ' ;
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' 'RUZ') RESIDU
('MOTS' 'UX' 'UY' 'UZ') ;
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' ;
* Axe de rotation: PCEN = 2.5 2.5 2.5 '+' VECX ;
VECX = 1.0 0.0 0.0 ;
'REPETER' BL1 3 ;
RN1 = RN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
GN1 = GN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
RETN1 = RETN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
RYN1 = RYN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
MAIL1 = 'EXTRAIRE' RN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RN1 '-' RN) / RN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' GN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((GN1 '-' GN) / (RN * CN)) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' RETN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RETN1 '-' RETN) / RETN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' RYN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RYN1 '-' RYN) / RN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
'FIN' BL1 ;
VECX = 0.0 0.0 1.0 ;
'REPETER' BL1 3 ;
RN1 = RN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
GN1 = GN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
RETN1 = RETN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
RYN1 = RYN1 'TOURNER' 90 PCEN (PCEN 'PLUS' VECX) ;
MAIL1 = 'EXTRAIRE' RN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RN1 '-' RN) / RN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' GN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((GN1 '-' GN) / (RN * CN)) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' RETN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RETN1 '-' RETN) / RETN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
MAIL1 = 'EXTRAIRE' RYN1 'MAILLAGE' ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RYN1 '-' RYN) / RN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
'FIN' BL1 ;
* Symetrie
MAIL1 RN1 = MAIL0 RN 'SYMETRIE'
'PLAN' (2.5 2.5 2.5) ((2.5 2.5 2.5) 'PLUS' (0.0 1.0 0.0))
((2.5 2.5 2.5) 'PLUS' (0.0 0.0 1.0)) ;
'ELIMINATION' 0.001 MAIL0 MAIL1 ;
ERR1 = 'MAXIMUM' ((RN1 - RN) / RN) 'ABS' ;
'SI' (ERR1 > 1.0D-10);
'ERREUR' 5;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer