Download stationary_discontinuity.dgibi
* fichier : stationary_discontinuity.dgibi
***********************************************************************
***********************************************************************
***********************************************************************
***** CALCUL DU TUBE A CHOC; CAS MULTIESPECE *****
***** GAZ multi-especes "THERMALLY PERFECT" *****
***** *****
***** FORMULATION VF COMPRESSIBLE EXPLICITE *****
***** Colella-Glaz, discontinuité de contact stationnaire *****
***** *****
***** A. BECCANTINI DRN/DMT/SEMT/LTMF FEVRIER 2000 *****
***********************************************************************
***********************************************************************
***********************************************************************
'MESSAGE' 'A mettre a jours' ;
'FIN' ;
'OPTION' 'DIME' 2 ;
'OPTION' 'ELEM' 'QUA4' ;
'OPTION' 'ISOV' 'SULI' ;
'OPTION' 'ECHO' 0 ;
'OPTION' 'TRAC' 'X';
GRAPH = FAUX ;
* GRAPH = VRAI ;
*
*** Methodes possibles :
*
* 'VLH'
* 'CG'
METO = 'CG' ;
*******************************************************
***** PROCEDURE POUR TESTER CONVERGENCE *****
*******************************************************
'DEBPROC' CALCUL ;
'ARGUMENT' RVX*'TABLE' ;
RV = RVX . 'EQEX' ;
*KDOMA = RV . 'DOMAINE' ;
KDOMA = RV . 'MODTOT' ;
RNi = RV . 'INCO' . 'RNI' ;
RN0 = RV . 'INCO' . 'RN0' ;
*
*** RN0 = solution à t = tn_M1;
*
DD = RV . 'PASDETPS' . 'NUPASDT' ;
NINTE = 20 ;
NN = DD/NINTE ;
LO = (DD-(NINTE*NN)) EGA 0 ;
'SI' ( LO ) ;
ELI = 'MAXIMUM' (RNi '-' RN0) 'ABS' ;
ELI = ('LOG' (ELI + 1.0E-20))/('LOG' 10.) ;
'MESSAGE' ('CHAINE'
'ITER ' DD ' ERREUR LINF ' ELI );
IT = 'PROG' DD ;
ER = 'PROG' ELI ;
RV . 'INCO' . 'IT' = (RV . 'INCO' . 'IT') 'ET' IT ;
RV . 'INCO' . 'ER' = (RV . 'INCO' . 'ER') 'ET' ER ;
'FINSI' ;
RV . 'INCO' . 'RN0' = 'COPIER' RNi;
'FINPROC' ;
*****************************************************
*****************************************************
** PROCEDURE EXEX POUR FORMULATION VF COMPRESSIBLE **
** CAS MULTIESPECES "THERMALLY PERFECT" **
*****************************************************
*****************************************************
'DEBPROC' EXEX ;
'ARGUMENT' RV*TABLE ;
*******************************************
**** Recherche de RV . *KONV . IDEUL ****
*******************************************
*
**** Nom de la table RV . *'KONV' -> NOMT
*
NBOP = 'DIME' (RV . 'LISTOPER' ) ;
'REPETER' BL1 NBOP;
MCEL = 'EXTRAIRE' &BL1 RV . 'LISTOPER';
'SI' ( 'EGA' MCEL 'KONV ');
NOMT = 'MOT' ('TEXTE' ('CHAINE' &BL1 MCEL));
'QUITTER' BL1;
'FINSI' ;
'FIN' BL1;
IEUL = RV . NOMT . 'KOPT' . 'IDEUL';
'SI' ('NON' (IEUL 'EGA' 3));
'MESSAGE' 'EULERMST ???';
'ERREUR' 21;
'FINSI' ;
* Mono-espece ou multi-especes
'SI' ('EXISTE' (RV . NOMT) 'ARG5') ;
LOGMUL = VRAI ;
'SINON' ;
LOGMUL = FAUX ;
'FINSI' ;
*
**** CL
*
LOGLIM = RV . 'INCO' . 'CLIM' ;
******************************************
**** Ordre en espace, ordre en temps ****
******************************************
ORD_ESP = RV . 'ORDREESP' ;
ORD_TPS = RV . 'ORDRETPS' ;
'MESSAGE' '--------------------------';
'MESSAGE' 'Ordre en Espace :' ord_esp;
'MESSAGE' 'Ordre en Temps :' ord_tps;
'MESSAGE' '--------------------------';
'SI' ((ORD_ESP 'EGA' 1) 'ET' (ORD_TPS 'EGA' 2));
'MESSAGE' ;
'MESSAGE' (CHAINE 'Ordre en Espace doit etre 2');
'MESSAGE' (CHAINE 'On impose ça.');
'MESSAGE' ;
RV . 'ORDREESP' = 2 ;
'MESSAGE' ;
'MESSAGE' '--------------------------';
'MESSAGE' 'Ordre en Espace :' ord_esp;
'MESSAGE' 'Ordre en Temps :' ord_tps;
'MESSAGE' '--------------------------';
'FINSI' ;
******************************
**** La table 'PASDETPS' ****
******************************
TPSI = RV . 'TPSI' ;
TFIN = RV . 'TFINAL';
RV . 'PASDETPS' . 'TPS' = TPSI;
*
**** DELTAT-1 est un argument de PRET (prediction)
* Donc on doit l'initialiser.
*
RV . 'PASDETPS' . 'DELTAT-1' = 0.0D0;
CFL = rv.'ALFA' ;
*********************
**** Les TABLES *****
*********************
*
**** RV . 'INCO'
* RV . 'DOMAINE'
* RV . 'KIZD'
* RV . 'KIZG'
*
**** RV . 'INCO' -> KINCO
*
KINCO = (RV . 'INCO') ;
*
**** RV . 'DOMAINE' -> KDOMA
*
*KDOMA = (RV . 'DOMAINE') ;
KDOMA = (RV . 'MODTOT') ;
KDOMA2 = (RV . 'DOMAINE') ;
*
**** RV . 'KIZD' contient les volumes des ELTs
*
'SI' ('NON' ('EXISTE' RV 'KIZD')) ;
'KDIA' RV ;
'FINSI' ;
*
***** RV . 'KIZG' contient les flux aux interfaces.
*
'SI' ('NON' ('EXISTE' RV 'KIZG')) ;
RV . 'KIZG' = 'TABLE' 'KIZG' ;
'FINSI' ;
'SI' LOGMUL ;
*********************************************************
**** Multi-especes, boucle Sur les Pas de Temps ****
*********************************************************
*
**** Evaluation de coeff pour le calcule des pentes
*
KINCO . 'V' KINCO . 'P' KINCO . 'T' KINCO . 'Y'
KINCO . 'GAMMA' = 'PRIM' 'PERFTEMP' (KINCO . 'IPGAZ')
(KINCO . 'RNI') (KINCO . 'GNI') (KINCO . 'ENI')
(KINCO. 'RYNI');
GRADR ALR COEFR = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'RNI');
GRADP ALP COEFP = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'P');
GRADV ALV COEFV = 'PENT' KDOMA 'CENTRE' 'EULEVECT' 'LIMITEUR'
(KINCO . 'V');
GRADY ALY COEFY = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'Y');
I = 0 ;
'REPETER' BLOC1 (RV . 'ITMA') ;
I = I + 1 ;
*
***** Les variables primitives
*
KINCO . 'V' KINCO . 'P' KINCO . 'T' KINCO . 'Y'
KINCO . 'GAMMA' = 'PRIM' 'PERFTEMP' (KINCO . 'IPGAZ')
(KINCO . 'RNI') (KINCO . 'GNI') (KINCO . 'ENI')
(KINCO. 'RYNI');
'SI' (ORD_ESP 'EGA' 1) ;
ROF VITF PF YF = 'PRET' 'PERFTEMP'
ORD_ESP ORD_TPS KDOMA (KINCO . 'IPGAZ') (KINCO . 'RNI')
(KINCO . 'V') (KINCO . 'P') (KINCO . 'Y') ;
'SINON';
*
***** Ordre 2 en espace => calcul des pentes
*
GRADR ALR = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'RNI') 'GRADGEO' COEFR ;
GRADP ALP = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'P') 'GRADGEO' COEFP ;
GRADV ALV = 'PENT' KDOMA 'CENTRE' 'EULEVECT' 'LIMITEUR'
(KINCO . 'V') 'GRADGEO' COEFV ;
GRADY ALY = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'Y') 'GRADGEO' COEFY ;
'SI' (ORD_TPS 'EGA' 1);
ROF VITF PF YF = 'PRET' 'PERFTEMP' ORD_ESP ORD_TPS
KDOMA (KINCO . 'IPGAZ')
(KINCO . 'RNI') GRADR ALR
(KINCO . 'V') GRADV ALV
(KINCO . 'P') GRADP ALP
(KINCO . 'Y') GRADY ALY ;
'SINON' ;
*
********* Ordre 2 en temps
*
ROF VITF PF YF = 'PRET' 'PERFTEMP' ORD_ESP ORD_TPS
KDOMA (KINCO . 'IPGAZ')
(KINCO . 'RNI') GRADR ALR
(KINCO . 'V') GRADV ALV
(KINCO . 'P') GRADP ALP
(KINCO . 'Y') GRADY ALY
(KINCO . 'GAMMA')
((RV . 'PASDETPS' . 'DELTAT-1')/2.0);
'FINSI' ;
'FINSI' ;
*
*********** Creation de MCHAML de type 'FACEL' pour les
* calcul de flux aux interfaces
KINCO . 'RNF' = ROF ;
KINCO . 'VITNF' = VITF ;
KINCO . 'PNF' = PF ;
KINCO . 'YF' = YF ;
*
********* Boucle sur les operateurs
*
'REPETER' BLOC2 NBOP ;
NOMPER = 'EXTRAIRE' &BLOC2 (RV . 'LISTOPER');
NOTABLE = 'MOT' ('TEXTE' ('CHAINE' &BLOC2 NOMPER) ) ;
('TEXTE' NOMPER) (RV . NOTABLE) ;
'FIN' BLOC2 ;
*
********* Mise a jour de la table RV . 'PASDETPS'
*
'SI' ('EXISTE' RV 'DTI');
DTI = 'MINIMUM'
('PROG' ((RV . 'DTI') '/' CFL)
(RV . 'PASDETPS' . 'DTCONV') );
'SINON';
DTI = (RV . 'PASDETPS' . 'DTCONV');
'FINSI';
RV . 'PASDETPS' . 'DELTAT' = DTI ;
TMPS = RV . 'PASDETPS' . 'TPS';
DTI0 = TFIN '-' TMPS;
DTI0 = DTI0 '/' CFL;
'SI' (DTI0 '<EG' DTI);
DTI = DTI0;
RV . 'PASDETPS' . 'DELTAT' = DTI ;
LOGQUIT = VRAI;
'SINON' ;
LOGQUIT = FAUX ;
'FINSI' ;
*
******** On avance au pas de temps suivant
*
* N.B. Tn+1 = Tn + (CFL * RV . 'PASDETPS' . 'DELTAT');
*
'AVCT' RV CFL 'IMPR' (RV.'FIDT') ;
*
******** On detrui les choses qui ne servent plus
*
*
* Les variables primitives
*
'DETR' ( KINCO . 'V');
'DETR' ( KINCO . 'P');
'DETR' ( KINCO . 'T');
'DETR' ( KINCO . 'Y');
'DETR' ( KINCO . 'GAMMA');
'OUBL' KINCO 'V';
'OUBL' KINCO 'P';
'OUBL' KINCO 'T';
'OUBL' KINCO 'Y';
'OUBL' KINCO 'GAMMA';
*
* Les MCHAML faces
*
'DETR' ROF ;
'DETR' VITF ;
'DETR' PF ;
'DETR' YF;
'OUBL' KINCO 'RNF';
'OUBL' KINCO 'VITNF';
'OUBL' KINCO 'PNF';
'OUBL' KINCO 'YF';
*
* Les pentes 'ET' les limiteurs
*
'SI' (ORD_ESP 'EGA' 2);
'DETR' GRADR;
'DETR' GRADP;
'DETR' GRADV;
'DETR' GRADY;
'DETR' ALR;
'DETR' ALP;
'DETR' ALV;
'DETR' ALY;
'FINSI' ;
*
******** On impose le CL
*
*
'SI' LOGLIM;
PROLIM RV ;
'FINSI' ;
'SI' LOGQUIT;
'QUITTER' BLOC1;
'FINSI';
'MENAGE' ;
'FIN' BLOC1 ;
********************************************
**** FIN de Boucle Sur les Pas de Temps ***
********************************************
'SINON' ;
*********************************************************
**** Mono-espece, boucle sur les pas de temps ****
*********************************************************
KINCO . 'V' KINCO . 'P' KINCO . 'T'
KINCO . 'GAMMA' = 'PRIM' 'PERFTEMP' (KINCO . 'IPGAZ')
(KINCO . 'RNI') (KINCO . 'GNI') (KINCO . 'ENI') ;
GRADR ALR COEFR = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'RNI');
GRADP ALP COEFP = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'P');
GRADV ALV COEFV = 'PENT' KDOMA 'CENTRE' 'EULEVECT' 'LIMITEUR'
(KINCO . 'V');
I = 0 ;
'REPETER' BLOC1 (RV . 'ITMA') ;
I = I + 1 ;
*
***** Les variables primitives
*
KINCO . 'V' KINCO . 'P' KINCO . 'T'
KINCO . 'GAMMA' = 'PRIM' 'PERFTEMP' (KINCO . 'IPGAZ')
(KINCO . 'RNI') (KINCO . 'GNI') (KINCO . 'ENI') ;
'SI' (ORD_ESP 'EGA' 1) ;
ROF VITF PF = 'PRET' 'PERFTEMP'
ORD_ESP ORD_TPS KDOMA (KINCO . 'IPGAZ') (KINCO . 'RNI')
(KINCO . 'V') (KINCO . 'P') ;
'SINON';
*
***** Ordre 2 en espace => calcul des pentes
*
GRADR ALR = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'RNI') 'GRADGEO' COEFR ;
GRADP ALP = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'LIMITEUR'
(KINCO . 'P') 'GRADGEO' COEFP ;
GRADV ALV = 'PENT' KDOMA 'CENTRE' 'EULEVECT' 'LIMITEUR'
(KINCO . 'V') 'GRADGEO' COEFV ;
'SI' (ORD_TPS 'EGA' 1);
ROF VITF PF = 'PRET' 'PERFTEMP' ORD_ESP ORD_TPS
KDOMA (KINCO . 'IPGAZ')
(KINCO . 'RNI') GRADR ALR
(KINCO . 'V') GRADV ALV
(KINCO . 'P') GRADP ALP ;
'SINON' ;
*
********* Ordre 2 en temps
*
ROF VITF PF = 'PRET' 'PERFTEMP' ORD_ESP ORD_TPS
KDOMA (KINCO . 'IPGAZ')
(KINCO . 'RNI') GRADR ALR
(KINCO . 'V') GRADV ALV
(KINCO . 'P') GRADP ALP
(KINCO . 'GAMMA')
((RV . 'PASDETPS' . 'DELTAT-1')/2.0);
'FINSI' ;
'FINSI' ;
*
*********** Creation de MCHAML de type 'FACEL' pour les
* calcul de flux aux interfaces
KINCO . 'RNF' = ROF ;
KINCO . 'VITNF' = VITF ;
KINCO . 'PNF' = PF ;
*
********* Boucle sur les operateurs
*
'REPETER' BLOC2 NBOP ;
NOMPER = 'EXTRAIRE' &BLOC2 (RV . 'LISTOPER');
NOTABLE = 'MOT' ('TEXTE' ('CHAINE' &BLOC2 NOMPER) ) ;
('TEXTE' NOMPER) (RV . NOTABLE) ;
'FIN' BLOC2 ;
*
********* Mise a jour de la table RV . 'PASDETPS'
*
'SI' ('EXISTE' RV 'DTI');
DTI = 'MINIMUM'
('PROG' ((RV . 'DTI') '/' CFL)
(RV . 'PASDETPS' . 'DTCONV') );
'SINON';
DTI = (RV . 'PASDETPS' . 'DTCONV');
'FINSI';
RV . 'PASDETPS' . 'DELTAT' = DTI ;
TMPS = RV . 'PASDETPS' . 'TPS';
DTI0 = TFIN '-' TMPS;
DTI0 = DTI0 '/' CFL;
'SI' (DTI0 '<EG' DTI);
DTI = DTI0;
RV . 'PASDETPS' . 'DELTAT' = DTI ;
LOGQUIT = VRAI;
'SINON' ;
LOGQUIT = FAUX ;
'FINSI' ;
*
******** On avance au pas de temps suivant
*
* N.B. Tn+1 = Tn + (CFL * RV . 'PASDETPS' . 'DELTAT');
*
'AVCT' RV CFL 'IMPR' (RV.'FIDT') ;
*
******** On detrui les choses qui ne servent plus
*
*
* Les variables primitives
*
'DETR' ( KINCO . 'V');
'DETR' ( KINCO . 'P');
'DETR' ( KINCO . 'T');
'DETR' ( KINCO . 'GAMMA');
'OUBL' KINCO 'V';
'OUBL' KINCO 'P';
'OUBL' KINCO 'T';
'OUBL' KINCO 'GAMMA';
*
* Les MCHAML faces
*
'DETR' ROF ;
'DETR' VITF ;
'DETR' PF ;
'OUBL' KINCO 'RNF';
'OUBL' KINCO 'VITNF';
'OUBL' KINCO 'PNF';
*
* Les pentes 'ET' les limiteurs
*
'SI' (ORD_ESP 'EGA' 2);
'DETR' GRADR;
'DETR' GRADP;
'DETR' GRADV;
'DETR' ALR;
'DETR' ALP;
'DETR' ALV;
'FINSI' ;
*
******** On impose le CL
*
*
'SI' LOGLIM;
PROLIM RV ;
'FINSI' ;
'SI' LOGQUIT;
'QUITTER' BLOC1;
'FINSI';
'MENAGE' ;
'FIN' BLOC1 ;
********************************************
**** FIN de Boucle Sur les Pas de Temps ***
********************************************
'FINSI' ;
'FINPROC' ;
*****************************************************
*****************************************************
** FIN PROCEDURE EXEX **
*****************************************************
*****************************************************
*****************************************************
*****************************************************
***** PROCEDURE PROLIM *****
*****************************************************
*****************************************************
*
*
**** Cas Euler
*
'DEBPROC' PROLIM ;
'ARGUMENT' RV*TABLE ;
*
**** RV . 'DOMAINE' -> KDOMA
*
*KDOMA = RV . 'DOMAINE' ;
KDOMA = RV . 'MODTOT' ;
KDOMA2 = RV . 'DOMAINE' ;
KINCO = RV . 'INCO' ;
KINCO . 'RNI' = 'KCHT' KDOMA 'SCAL' 'CENTRE'
(KINCO . 'RNI') (KINCO . 'RNLIM');
KINCO . 'GNI' = 'KCHT' KDOMA 'VECT' 'CENTRE'
(KINCO . 'GNI') (KINCO . 'GNLIM');
KINCO . 'ENI' = 'KCHT' KDOMA 'SCAL' 'CENTRE'
(KINCO . 'ENI') (KINCO . 'ENLIM');
KINCO . 'RYNI' = 'KCHT' KDOMA 'SCAL' 'CENTRE'
'COMP' 'N2 ' (KINCO . 'RYNI')
(KINCO . 'RYNLIM') ;
'FINPROC' ;
*****************************************************
*****************************************************
***** FIN PROCEDURE PROLIM *****
*****************************************************
*****************************************************
************
* MAILLAGE *
************
RAF = 2 ;
NY = 4 ;
NX = 3 '*' RAF ;
L = 1.0D0 ;
DX = L '/' NX '/' 2.0D0;
H = NY '*' DX ;
xD = 0.5D0 '*' L ;
xG = -1.0D0 '*' xD ;
yH = 0.5D0 '*' H ;
yB = -1.0D0 '*' yH ;
A1 = xG yB ;
A2 = 0.0D0 yB ;
A3 = xD yB ;
A4 = xD yH ;
A5 = 0.0D0 yH ;
A6 = xG yH ;
VECTG = XG 0.0D0 ;
VECTD = XD 0.0D0 ;
xBG = xG '-' DX;
XBD = xD '+' DX;
B1 = xBG yB;
B2 = xBG yH;
B3 = xBD yB;
B4 = xBD yH;
BAS1 = A1 'DROI' NX A2 ;
BAS2 = A2 'DROI' NX A3 ;
HAU2 = A4 'DROI' NX A5 ;
HAU1 = A5 'DROI' NX A6 ;
LATG = B1 'DROI' NY B2 ;
LAT1 = A1 'DROI' NY A6 ;
LAT12 = A2 'DROI' NY A5 ;
LAT2 = A3 'DROI' NY A4 ;
LATD = B3 'DROI' NY B4 ;
DOM1 = LAT12 'TRANSLATION' NX VECTG ;
DOM2 = LAT12 'TRANSLATION' NX VECTD ;
VECTFG = (-1.0D0 '*' DX) 0.0D0;
VECTFD = DX 0.0D0;
FRONTG = LAT1 'TRANSLATION' 1 VECTFG;
FRONTG = FRONTG 'COUL' ROUG ;
FRONTD = LAT2 'TRANSLATION' 1 VECTFD;
FRONTD = FRONTD 'COUL' VERT ;
*
*** Rotation
*
ANGLE = 30.0D0;
ORIG = 0.0D0 0.0D0;
'MESSAGE';
'MESSAGE' (CHAIN 'ANGLE = ' ANGLE);
'MESSAGE';
DOM1 = DOM1 'TOURNER' ANGLE ORIG;
DOM2 = DOM2 'TOURNER' ANGLE ORIG;
FRONTG = FRONTG 'TOURNER' ANGLE ORIG;
FRONTD = FRONTD 'TOURNER' ANGLE ORIG;
DOMINT = DOM1 'ET' DOM2 ;
'ELIMINATION' DOMINT 1D-6;
FRONT = FRONTG 'ET' FRONTD ;
'ELIMINATION' FRONT 1D-6;
DOMTOT = DOMINT 'ET' FRONT;
'ELIMINATION' DOMTOT 1D-6;
**********************
*** OBJETS MODELES ***
**********************
MDOMTOT = 'CHANGER' DOMTOT 'QUAF' ;
MDOMINT = 'CHANGER' DOMINT 'QUAF' ;
MDOM1 = 'CHANGER' DOM1 'QUAF' ;
MDOM2 = 'CHANGER' DOM2 'QUAF' ;
MFRONTG = 'CHANGER' FRONTG 'QUAF' ;
MFRONTD = 'CHANGER' FRONTD 'QUAF' ;
MFRONT = 'CHANGER' FRONT 'QUAF' ;
'ELIM' (MDOMTOT 'ET' MDOMINT 'ET' MDOM1 'ET' MDOM2 'ET'
MFRONTG 'ET' MFRONTD 'ET' MFRONT) 1.E-5;
MDNS = 'EULER' ;
$DOMTOT = 'MODE' MDOMTOT MDNS ;
$DOMINT = 'MODE' MDOMINT MDNS ;
$DOM1 = 'MODE' MDOM1 MDNS ;
$DOM2 = 'MODE' MDOM2 MDNS ;
$FRONTG = 'MODE' MFRONTG MDNS ;
$FRONTD = 'MODE' MFRONTD MDNS ;
$FRONT = 'MODE' MFRONT MDNS ;
*
******* Creation de la ligne Utilisée pour le Post-Traitement
* reliant les points centres
*
XINIT = XG '+' (0.5D0 '*' DX) ;
YINIT = YB '+' (0.5D0 '*' DX) ;
XFIN = XD '-' (0.5D0 '*' DX) ;
YFIN = YINIT ;
PINI = XINIT YINIT;
PFIN = XFIN YFIN ;
IAUX = (2 '*' NX) '-' 1 ;
COURB = PINI 'DROIT' IAUX PFIN;
COURB = COURB 'TOURNER' ANGLE ORIG;
COURB = COURB 'COULEUR' 'VERT';
'ELIMINATION' 0.001 Courb ('DOMA' $DOMTOT 'CENTRE') ;
'SI' GRAPH ;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' COURB)
'TITRE' ('CHAINE' 'Maillage ');
'FINSI' ;
***********************
**** LA TABLE PGAZ ****
***********************
PGAZ = 'TABLE' ;
*
**** Ordre des polynoms cv_i
*
PGAZ . 'NORD' = 4 ;
*
**** Especes qui sont dans les equations d'Euler
*
PGAZ . 'ESPEULE' = 'MOTS' 'N2 ' ;
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'O2 ';
*
PGAZ . 'N2 ' = 'TABLE' ;
PGAZ . 'O2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'N2 ' . 'R' = 296.8 ;
PGAZ . 'O2 ' . 'R' = 259.8 ;
*
**** Regressions polynomials
*
PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.350522002 -0.000128294865
2.33636971E-08 -1.53304905E-12;
PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
8.78233606E-09 -3.05514485E-13 ;
*
**** "Enthalpies" (ou energies) de formations a OK (J/Kg)
* Note: ce sont des entites numeriques
*
PGAZ . 'O2 ' . 'H0K' = -2.634D5 ;
PGAZ . 'N2 ' . 'H0K' = -2.953D5 ;
**********************************************************
*** Fin PGAZ *********************************************
**********************************************************
*
*** Etatg
*
*
s = 0 ;
YN2g = 0.77785 ;
UMYN2g = 1.0 '-' YN2g ;
RTOTg = ((PGAZ . 'O2 ' . 'R') * UMYN2g) '+'
((PGAZ . 'N2 ' . 'R') * YN2G ) ;
pg = 1.0D6;
rog = 1.0D2 ;
Tg = pg '/' (rog '*' RTOTg) ;
ung = s ;
utg = 0.0 ;
*
A1 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 1) '*' UMYN2g) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 1) '*' YN2g);
A2 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 2) '*' UMYN2g) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 2) '*' YN2g);
A3 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 3) '*' UMYN2g) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 3) '*' YN2g);
A4 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 4) '*' UMYN2g) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 4) '*' YN2g);
A5 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 5) '*' UMYN2g) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 5) '*' YN2g);
T2 = Tg '*' Tg ;
T3 = T2 '*' Tg ;
T4 = T3 '*' Tg ;
T5 = T4 '*' Tg ;
ETHERg = ((A1 * Tg) '+' ((A2 '/' 2.0) * T2) '+'
((A3 '/' 3.0) * T3) '+' ((A4 '/' 4.0) * T4) '+'
((A5 '/' 5.0) * T5)) ;
retg = rog '*' (ETHERg '+' (0.5 '*' ung '*' ung)) ;
rouxg = ((ung '*' ('COS' ANGLE)) '-'
(utg '*' ('SIN' ANGLE))) '*' rog ;
rouyg = ((ung '*' ('SIN' ANGLE)) '+'
(utg '*' ('COS' ANGLE))) '*' rog;
*
*** Etatd
*
YN2d = 0.77785 ;
UMYN2d = 1.0 '-' YN2d ;
RTOTd = ((PGAZ . 'O2 ' . 'R') * UMYN2d) '+'
((PGAZ . 'N2 ' . 'R') * YN2D ) ;
pd = pg ;
rod = rog '/' 5.0 ;
Td = pd '/' (rod '*' RTOTd) ;
und = ung ;
utd = utg ;
*
A1 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 1) '*' UMYN2d) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 1) '*' YN2d);
A2 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 2) '*' UMYN2d) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 2) '*' YN2d);
A3 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 3) '*' UMYN2d) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 3) '*' YN2d);
A4 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 4) '*' UMYN2d) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 4) '*' YN2d);
A5 = (('EXTRAIRE' (PGAZ . 'O2 ' . 'A') 5) '*' UMYN2d) '+'
(('EXTRAIRE' (PGAZ . 'N2 ' . 'A') 5) '*' YN2d);
T2 = Td '*' Td ;
T3 = T2 '*' Td ;
T4 = T3 '*' Td ;
T5 = T4 '*' Td ;
ETHERd = ((A1 * Td) '+' ((A2 '/' 2.0) * T2) '+'
((A3 '/' 3.0) * T3) '+' ((A4 '/' 4.0) * T4) '+'
((A5 '/' 5.0) * T5)) ;
rouxd = ((und '*' ('COS' ANGLE)) '-'
(utd '*' ('SIN' ANGLE))) '*' rod;
rouyd = ((und '*' ('SIN' ANGLE)) '+'
(utd '*' ('COS' ANGLE))) '*' rod;
retd = rod '*' (ETHERd '+' (0.5 '*' und '*' und)) ;
*
*** ro
*
RO_f1 = 'KCHT' $FRONTG 'SCAL' 'CENTRE' rog ;
RO_f2 = 'KCHT' $FRONTD 'SCAL' 'CENTRE' rod ;
RO_f = 'KCHT' $FRONT 'SCAL' 'CENTRE' (RO_f1 'ET' RO_f2) ;
RO_i1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog;
RO_i2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod;
RO_i = 'KCHT' $DOMINT 'SCAL' 'CENTRE' (RO_i1 'ET' RO_i2);
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RO_i 'ET' RO_f) ;
*
*** ro u, ro v
*
GN_f1 = 'KCHT' $FRONTG 'VECT' 'CENTRE' (rouxg rouyg);
GN_f2 = 'KCHT' $FRONTD 'VECT' 'CENTRE' (rouxd rouyd);
GN_f = 'KCHT' $FRONT 'VECT' 'CENTRE' (GN_f1 'ET' GN_f2);
GN_i1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (rouxg rouyg);
GN_i2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (rouxd rouyd);
GN_i = 'KCHT' $DOMINT 'VECT' 'CENTRE' (GN_i1 'ET' GN_i2);
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN_i 'ET' GN_f) ;
*
*** ro e
*
RE_f1 = 'KCHT' $FRONTG 'SCAL' 'CENTRE' retg ;
RE_f2 = 'KCHT' $FRONTD 'SCAL' 'CENTRE' retd ;
RE_f = 'KCHT' $FRONT 'SCAL' 'CENTRE' (RE_f1 'ET' RE_f2) ;
RE_i1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' retg ;
RE_i2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' retd ;
RE_i = 'KCHT' $DOMINT 'SCAL' 'CENTRE' (RE_i1 'ET' RE_i2);
REN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RE_i ET RE_f) ;
*
*** ro y
*
RYN =('KCHT' $FRONTG 'SCAL' 'CENTRE' 'COMP' 'N2 ' (rog '*' YN2g)) 'ET'
('KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'N2 ' (rog '*' YN2g)) 'ET'
('KCHT' $FRONTD 'SCAL' 'CENTRE' 'COMP' 'N2 ' (rod '*' YN2d)) 'ET'
('KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'N2 ' (rod '*' YN2d));
********************************************************
*** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM ***
********************************************************
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
*
**** Les debits dans le repaire tube
*
GNX = 'EXCO' 'UX' GN;
GNY = 'EXCO' 'UY' GN;
GNN = (GNX * ('COS' ANGLE)) '+'
(GNY * ('SIN' ANGLE));
GNT = (GNY * ('COS' ANGLE)) '-'
(GNX * ('SIN' ANGLE));
*
*** Les variables primitives
*
VN PN TN YN GAMN = 'PRIM' 'PERFTEMP' PGAZ
RN GN REN RYN ;
RN0 = 'COPIER' RN ;
*
*** GRAPHIQUE DES C.I.
*
'SI' GRAPH ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
CHM_GNN = 'KCHA' $DOMTOT 'CHAM' GNN ;
CHM_GNT = 'KCHA' $DOMTOT 'CHAM' GNT ;
CHM_REN = 'KCHA' $DOMTOT 'CHAM' REN ;
CHM_RYN = 'KCHA' $DOMTOT 'CHAM' RYN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' 0.0);
TRAC CHM_TN MOD1 'TITR' ('CHAINE' 'T at t=' 0.0);
TRAC CHM_PN MOD1 'TITR' ('CHAINE' 'P at t=' 0.0);
TRAC CHM_REN MOD1 'TITR' ('CHAINE' 'ROET at t=' 0.0);
TRAC CHM_GNN MOD1 'TITR' ('CHAINE' 'ROUN at t=' 0.0);
TRAC CHM_GNT MOD1 'TITR' ('CHAINE' 'ROUT at t=' 0.0);
TRAC CHM_RYN MOD1 'TITR' ('CHAINE' 'RYN at t=' 0.0);
'FINSI' ;
***********************
**** La table EQEX ****
***********************
CFL = 0.5D0 ;
NITER = 100 ;
* TFIN = 0.0005 ;
TFIN = 100 ;
RV = 'EQEX' ('DOMA' $DOMTOT 'TABLE') 'ITMA' NITER 'ALFA' CFL
'TPSI' 0. 'TFINAL' TFIN
*
***** Option VF = volumes finis
* CONS = conservative
* EULER = euler monoespece
* VANLEER = methode
*
'OPTI' 'VF' 'CONS' 'EULERMST' METO
*
***** Procedure CALCUL
*
'ZONE' $DOMTOT 'OPER' 'CALCUL'
*
***** Operateur 'KONV'
*
'ZONE' $DOMTOT 'OPER' 'KONV'
*
***** Arguments de 'KONV' (maximum 8 chatacters)
*
'IPGAZ' 'RNF' 'VITNF' 'PNF' 'YF'
*
***** LIST des inconnues
*
'INCO' 'RNI' 'GNI' 'ENI' 'RYNI' ;
*
*** La table RV . INCO (de soustype INCO);
*
RV . 'INCO' = TABLE 'INCO' ;
*
*** Stokage des inconnues des equations d'Euler.
*
RV . 'INCO' . 'RNI' = 'COPIER' RN ;
RV . 'INCO' . 'GNI' = 'COPIER' GN ;
RV . 'INCO' . 'ENI' = 'COPIER' REN ;
RV . 'INCO' . 'RYNI' = 'COPIER' RYN ;
*
*** Le gaz
*
RV . 'INCO' . 'IPGAZ' = PGAZ ;
*
*** CONDITIONS LIMITS
*
RV . 'INCO' . 'CLIM' = VRAI ;
RV . 'INCO' . 'RNLIM' = 'KCHT' $FRONT 'SCAL' 'CENTRE' RN ;
RV . 'INCO' . 'GNLIM' = 'KCHT' $FRONT 'VECT' 'CENTRE' GN ;
RV . 'INCO' . 'ENLIM' = 'KCHT' $FRONT 'SCAL' 'CENTRE' REN ;
RV . 'INCO' . 'RYNLIM' = 'KCHT' $FRONT 'SCAL' 'CENTRE'
'COMP' 'N2 ' RYN ;
*
**** Ordre en espace
* Ordre en temps
*
IE = 2 ;
IT = 2 ;
RV . 'ORDREESP' = IE ;
RV . 'ORDRETPS' = IT ;
*
**** Test de la convergence
*
RV . 'INCO' . 'RN0' = 'COPIER' RN ;
RV . 'INCO' . 'IT' = 'PROG' ;
RV . 'INCO' . 'ER' = 'PROG' ;
*???
RV . 'MODTOT' = $DOMTOT ;
*
***********************************
*** Execution EXEX ***
***********************************
*
'MESSAGE' ;
'MESSAGE' '**************************';
'MESSAGE' ('CHAINE' 'METHODE : ' METO) ;
'MESSAGE' '**************************';
'MESSAGE' ;
'TEMPS' ;
EXEX RV ;
'TEMPS' ;
TFIN = RV. 'PASDETPS'. 'TPS';
*
**** SOLUTIONS
*
*
**** Les variables conservatives
*
RN = 'COPIER' (RV . 'INCO' . 'RNI');
GN = 'COPIER' (RV . 'INCO'. 'GNI');
REN = 'COPIER' (RV . 'INCO'. 'ENI');
RYN = 'COPIER' (RV . 'INCO'. 'RYNI');
*
**** Les variables primitives
*
VN PN TN YN GAMN = 'PRIM' 'PERFTEMP' PGAZ
RN GN REN RYN ;
*
**** La vitesse dans le repaire tube
*
VNX = 'EXCO' 'UX' VN;
VNY = 'EXCO' 'UY' VN;
VNN = (VNX * ('COS' ANGLE)) '+'
(VNY * ('SIN' ANGLE));
VNT = (VNY * ('COS' ANGLE)) '-'
(VNX * ('SIN' ANGLE));
GNN = VNN * RN;
GNT = VNT * RN;
*
*** GRAPHIQUE DES SOLUTIONS
*
'SI' GRAPH ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN;
CHM_GNN = 'KCHA' $DOMTOT 'CHAM' GNN ;
CHM_GNT = 'KCHA' $DOMTOT 'CHAM' GNT ;
CHM_REN = 'KCHA' $DOMTOT 'CHAM' REN;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' TFIN);
TRAC CHM_REN MOD1 'TITR' ('CHAINE' 'ROET at t=' TFIN);
TRAC CHM_GNN MOD1 'TITR' ('CHAINE' 'ROUN at t=' TFIN);
TRAC CHM_GNT MOD1 'TITR' ('CHAINE' 'ROUT at t=' TFIN);
'FINSI' ;
*
*** Objects evolutions
*
xx yy = 'COORDONNEE' Courb;
ss = 'KOPS' ('KOPS' ('COS' ANGLE) '*' xx) '+'
('KOPS' ('SIN' ANGLE) '*' yy);
evxx = 'EVOL' 'CHPO' ss Courb ;
lxx = 'EXTRAIRE' evxx 'ORDO' ;
x0 = 'MINIMUM' lxx;
x1 = 'MAXIMUM' lxx ;
* ro
evro = 'EVOL' 'CHPO' RN Courb ;
lro = 'EXTRAIRE' evro 'ORDO' ;
evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro;
tro = CHAINE '1D ' METO ' : RO IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
* u
evu = 'EVOL' 'CHPO' VNN Courb ;
lu = 'EXTRAIRE' evu 'ORDO' ;
evu = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
tu = CHAINE '1D ' METO ' : U IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
* v
evv = 'EVOL' 'CHPO' VNT Courb ;
lv = 'EXTRAIRE' evv 'ORDO' ;
evv = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
tv = CHAINE '1D ' METO ' : V IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
'SI' GRAPH ;
'DESSIN' evv 'TITRE' tv 'XBOR' x0 x1;
'FINSI' ;
* p
evp = 'EVOL' 'CHPO' PN Courb ;
lp = 'EXTRAIRE' evp 'ORDO' ;
evp = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
tp = CHAINE '1D ' METO ' : P IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
* gamma
evgam = 'EVOL' 'CHPO' GAMN Courb ;
lgam = 'EXTRAIRE' evgam 'ORDO' ;
evgam = 'EVOL' 'MANU' 'x' lxx 'gamma' lgam ;
tgam = CHAINE '1D ' METO ' : GAMMA IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
* T
evT = 'EVOL' 'CHPO' TN Courb ;
lT = 'EXTRAIRE' evT 'ORDO' ;
evT = 'EVOL' 'MANU' 'x' lxx 'T' lT ;
tT = 'CHAINE' '1D ' METO ' : T IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
* Fractions massiques
EVY = 'TABLE' ;
TY = CHAINE '1D ' METO ' : Y IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
YOLD = 'COPIER' YN ;
NOMECO = 'EXTR' YOLD 'COMP';
NUMECO = 'DIME' NOMECO ;
'REPETER' BLOCCO NUMECO;
NOMCEL = 'EXTRAIRE' NOMECO &BLOCCO ;
Y0 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE'
'COMP' NOMCEL
('EXCO' NOMCEL YOLD NOMCEL) ;
N0 = NOMCEL ;
EVY0 = 'EVOL' 'CHPO' Y0 Courb ;
LY0 = 'EXTRAIRE' EVY0 'ORDO' ;
EVY0 = 'EVOL' 'MANU' 'x' lxx N0 LY0 ;
EVY . &BLOCCO = EVY0;
'FIN' BLOCCO;
EVERR = 'EVOL' 'MANU' 'IT' (RV . 'INCO' . 'IT')
'LOG|E|inf' (RV . 'INCO' . 'ER') ;
'SI' GRAPH ;
'DESSIN' evro 'TITRE' tro 'XBOR' x0 x1 ;
'DESSIN' evT 'TITRE' tT 'XBOR' x0 x1 ;
'DESSIN' evp 'TITRE' tp 'XBOR' x0 x1 ;
'DESSIN' evu 'TITRE' tu 'XBOR' x0 x1 ;
'DESSIN' evgam 'TITRE' tgam 'XBOR' x0 x1 ;
'DESSIN' EVERR 'TITRE' ('CHAINE' METO ' IE =' IE ' IT =' IT) ;
'FINSI' ;
* TESTS
ERRO = ('MAXIMUM' (RN '-' RN0) 'ABS') '/'
('MAXIMUM' RN0);
'SI' (ERRO > 1.0D-4) ;
'MESSAGE' 'CG perde la discontinuite de contact stationnaire' ;
'ERREUR' 5;
'FINSI';
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer