Download stationary_shock.dgibi
* fichier : stationary_shock.dgibi
***********************************************************************
***********************************************************************
***********************************************************************
*** CALCUL DU TUBE A CHOC; CAS MULTIESPECE ****
*** GAZ mono-especes "THERMALLY PERFECT" ****
*** ****
*** FORMULATION VF COMPRESSIBLE EXPLICITE ****
*** DIFFERENTS SOLVEURS ****
*** ****
*** Choc 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' ;
***********************
**** LA TABLE PGAZ ****
***********************
PGAZ = 'TABLE' ;
*
**** Ordre des polynoms cv_i
*
PGAZ . 'NORD' = 4 ;
*
**** La seul espece
*
PGAZ . 'ESPNEULE' = 'H2 ';
*
PGAZ . 'H2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'H2 ' . 'R' = 4130.0 ;
*
**** Regressions polynomials
*
PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
*
**** "Enthalpies" (ou energies) de formations a OK (J/Kg)
* Note: ce sont des entites numeriques
*
PGAZ . 'H2 ' . 'H0K' = -4.195D6 ;
*
*** Fin PGAZ
*
*
**** Choc stationnaire
*
pg = 1.0D5 ;
rog = 1.0 ;
Tg = (pg '*' rog) '/' (PGAZ . 'H2 ' . 'R') ;
* sstren > 1.0
sstren = 4.01 ;
Td = sstren '*' Tg ;
* Det. ug, rod, ud
A0H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 1 ;
A1H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 2 ;
A2H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 3 ;
A3H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 4 ;
A4H2 = 'EXTRAIRE' (PGAZ . 'H2 ' . 'A') 5 ;
R = PGAZ . 'H2 ' . 'R' ;
T2 = Tg '*' Tg ;
T3 = T2 '*' Tg ;
T4 = T3 '*' Tg ;
T5 = T4 '*' Tg ;
ETHERg = ((A0H2 * Tg) '+' ((A1H2 '/' 2.0) * T2) '+'
((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
((A4H2 '/' 5.0) * T5)) ;
T2 = Td '*' Td ;
T3 = T2 '*' Td ;
T4 = T3 '*' Td ;
T5 = T4 '*' Td ;
ETHERd = ((A0H2 * Td) '+' ((A1H2 '/' 2.0) * T2) '+'
((A2H2 '/' 3.0) * T3) '+' ((A3H2 '/' 4.0) * T4) '+'
((A4H2 '/' 5.0) * T5)) ;
RTg = R '*' Tg ;
RTd = R '*' Td ;
* x = rog '/' rod `in (0,1]
NDEF = 100 ;
DX = 1.0 /NDEF ;
X = 0.0 ;
AEQ2 = RTg ;
AEQ1 = (2.0 '*' (ETHERd '-' ETHERg)) '+' RTd '-' RTg ;
AEQ0 = -1.0 '*' RTd ;
LIX = 'PROG' X ;
LIF = 'PROG' AEQ0 ;
'REPETER' BL1 NDEF ;
X = X '+' DX ;
X2 = X '*' X ;
FUN = (AEQ2 '*' X2) '+' (AEQ1 '*' X) '+' AEQ0 ;
LIX = LIX 'ET' ('PROG' X) ;
LIF = LIF 'ET' ('PROG' FUN) ;
'FIN' BL1 ;
EVFUN = 'EVOL' 'MANU' 'rog/rod' LIX 'fun' LIF ;
'SI' FAUX ;
'DESSIN' EVFUN 'MIMA' ;
'FINSI' ;
* Selection de la racine:
* XSOL \in (0,1]
XSOL = ((((AEQ1 '**' 2) '-' (4.0 '*' AEQ0 * AEQ2)) '**' 0.5) '-' AEQ1)
'/' (2.0 * AEQ2) ;
*
rod = rog '/' XSOL ;
pd = RTd '*' rod;
RC1 = 2.0 '*' (ETHERd '+' (R '*' Td) '-' (ETHERg '+' (R '*' Tg))) ;
ud = RC1 '*' (XSOL '**' 2) '/' (1.0 '-' (XSOL '**' 2)) ;
ud = ud '**' 0.5 ;
ug = (ud '*' rod) '/' rog ;
* Control de Rankine-Hugoniot
err1 = 'ABS' (((rog '*' ug) '-' (rod '*' ud)) '/' (rog '*' ug)) ;
err2 = 'ABS' (((rog '*' ug '*' ug) '+' pg '-'
(rod '*' ud '*' ud) '-' pd)
'/' ((rog '*' ug '*' ug) '+' pg)) ;
err3 = 'ABS' (((ETHERg '+' RTg '+' (0.5 '*' ug '*' ug)) '-'
(ETHERd '+' RTd '+' (0.5 '*' ud '*' ud))) '/'
(ETHERd '+' RTd '+' (0.5 '*' ud '*' ud))) ;
ERMAX = 'MAXIMUM' ( 'PROG' err1 err2 err3) ;
'SI' (ERMAX > 1.0d-10);
'ERREUR' 'No Rankine-Hugoniot';
'FINSI' ;
*******************************************************
***** 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 = 5;
NN = DD/NINTE ;
LO = (DD-(NINTE*NN)) EGA 0 ;
'SI' ( LO ) ;
RNi = 'KCHT' KDOMA 'SCAL' 'CENTRE' RNi ;
RN0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' RN0 ;
ERR = 'KOPS' RNi '-' RN0 ;
ERR = 'KOPS' ERR '/' RN0;
ELI = 'MAXIMUM' ERR 'ABS' ;
'MESSAGE' eli ;
ELI = ('LOG' (ELI + 1.0E-10))/('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 mono-espece
*
'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');
'FINPROC' ;
*****************************************************
*****************************************************
***** FIN PROCEDURE PROLIM *****
*****************************************************
*****************************************************
************
* MAILLAGE *
************
RAF = 4 ;
NY = RAF ;
NX = 6 '*' 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 = (-2.0D0 '*' DX) 0.0D0;
VECTFD = (2.0 '*' DX) 0.0D0;
FRONTG = LAT1 'TRANSLATION' 2 VECTFG;
FRONTG = FRONTG 'COUL' ROUG ;
FRONTD = LAT2 'TRANSLATION' 2 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' ;
*
*** Etats exacte
*
ro1 = rog ;
u1 = ug ;
p1 = pg ;
ETHER1 = ETHERg ;
T1 = Tg ;
ro2 = rod ;
u2 = ud ;
p2 = pd ;
ETHER2 = ETHERd ;
T2 = Td ;
*
**** On inverse les etats, pour controller qu'on a le meme resultats
*
'SI' VRAI ;
rog = ro2 ;
ug = -1.0 '*' u2 ;
pg = p2 ;
ETHERg = ETHER2 ;
Tg = T2 ;
rod = ro1 ;
ud = -1.0 '*' u1 ;
pd = p1 ;
ETHERd = ETHER1 ;
Td = T1 ;
'FINSI' ;
*
*** Etat gauche (rog, ug, ETHERg, Tg , Pg deja definis)
*
ung = ug ;
utg = 0.0 ;
retg = rog '*' (ETHERg '+' (0.5 '*' ung '*' ung)) ;
*
rouxg = ((ung '*' ('COS' ANGLE)) '-'
(utg '*' ('SIN' ANGLE))) '*' rog ;
rouyg = ((ung '*' ('SIN' ANGLE)) '+'
(utg '*' ('COS' ANGLE))) '*' rog;
*
*** Etat droite
*
und = ud ;
utd = 0.0D0 ;
retd = rod '*' (ETHERd '+' (0.5 '*' und '*' und)) ;
rouxd = ((und '*' ('COS' ANGLE)) '-'
(utd '*' ('SIN' ANGLE))) '*' rod;
rouyd = ((und '*' ('SIN' ANGLE)) '+'
(utd '*' ('COS' ANGLE))) '*' rod;
*
*** 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) ;
********************************************************
*** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM ***
********************************************************
MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
*
**** Les variables primitives
*
VN PN TN GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN REN ;
*
**** Enthalpie totale
*
HTN = ((REN '+' PN) '/' RN) ;
*
**** Les vitesses dans le repaire tube
*
VX = 'EXCO' 'UX ' VN ;
VY = 'EXCO' 'UY ' VN ;
VNN = (VX * ('COS' ANGLE)) '+'
(VY * ('SIN' ANGLE)) ;
VNT = (VY * ('COS' ANGLE)) '-'
(VX * ('SIN' ANGLE)) ;
*
*** GRAPHIQUE DES C.I.
*
'SI' GRAPH ;
* 'SI' FAUX ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_VNN = 'KCHA' $DOMTOT 'CHAM' VNN ;
CHM_VNT = 'KCHA' $DOMTOT 'CHAM' VNT ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
CHM_HTN = 'KCHA' $DOMTOT 'CHAM' HTN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' 0.0);
TRAC CHM_PN MOD1 'TITR' ('CHAINE' 'P at t=' 0.0);
TRAC CHM_VNN MOD1 'TITR' ('CHAINE' 'VNN at t=' 0.0);
TRAC CHM_VNT MOD1 'TITR' ('CHAINE' 'VNT at t=' 0.0);
TRAC CHM_TN MOD1 'TITR' ('CHAINE' 'T at t=' 0.0);
TRAC CHM_HTN MOD1 'TITR' ('CHAINE' 'ht at t=' 0.0);
'FINSI' ;
***********************
**** La table EQEX ****
***********************
CFL = 0.8D0 ;
NITER = 50 ;
TFIN = 10000 ;
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'
*
***** LIST des inconnues
*
'INCO' 'RNI' 'GNI' 'ENI' ;
*
*** 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 ;
*
*** Le gaz
*
RV . 'INCO' . 'IPGAZ' = PGAZ ;
*
*** CONDITIONS LIMITS
*
RV . 'INCO' . 'CLIM' = VRAI ;
RV . 'INCO' . 'RNLIM' = 'COPIER' RO_f ;
RV . 'INCO' . 'GNLIM' = 'COPIER' GN_f ;
RV . 'INCO' . 'ENLIM' = 'COPIER' RE_f ;
*
**** Ordre en espace
* Ordre en temps
*
IE = 2 ;
IT = 2 ;
RV . 'ORDREESP' = IE ;
RV . 'ORDRETPS' = IT ;
*
**** Test de la convergence
*
RN0 = 'COPIER' RN ;
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');
*
**** Les variables primitives
*
VN PN TN GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN REN ;
*
**** Enthalpie totale
*
HTN = ((REN '+' PN) '/' RN) ;
*
**** Les vitesses dans le repaire tube
*
VX = 'EXCO' 'UX ' VN ;
VY = 'EXCO' 'UY ' VN ;
VNN = (VX * ('COS' ANGLE)) '+'
(VY * ('SIN' ANGLE)) ;
VNT = (VY * ('COS' ANGLE)) '-'
(VX * ('SIN' ANGLE)) ;
*
*** GRAPHIQUE DES SOLUTIONS
*
'SI' GRAPH ;
* 'SI' FAUX ;
*
*** CREATION DE CHAMELEM
*
CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
CHM_VNN = 'KCHA' $DOMTOT 'CHAM' VNN ;
CHM_VNT = 'KCHA' $DOMTOT 'CHAM' VNT ;
CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
CHM_TN = 'KCHA' $DOMTOT 'CHAM' TN ;
CHM_HTN = 'KCHA' $DOMTOT 'CHAM' HTN ;
TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RO at t=' TFIN);
TRAC CHM_PN MOD1 'TITR' ('CHAINE' 'P at t=' TFIN);
TRAC CHM_VNN MOD1 'TITR' ('CHAINE' 'VNN at t=' TFIN);
TRAC CHM_VNT MOD1 'TITR' ('CHAINE' 'VNT at t=' TFIN);
TRAC CHM_TN MOD1 'TITR' ('CHAINE' 'T at t=' TFIN);
TRAC CHM_HTN MOD1 'TITR' ('CHAINE' 'ht 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' ;
* ht
evht = 'EVOL' 'CHPO' HTN Courb ;
lht = 'EXTRAIRE' evht 'ORDO' ;
evht = 'EVOL' 'MANU' 'x' lxx 'ht' lht ;
tht = CHAINE '1D ' METO ' : ht IT ' IT ' IE ' IE
' tmps ' TFIN ' elem ' 'QUA4' ;
*
*** Solution analitique
lpan = 'PROG' ;
lroan = 'PROG' ;
luan = 'PROG' ;
lhtan = 'PROG' ;
hg = ETHERg '+' (pg '/' rog) '+' (0.5 '*' ug '*' ug) ;
'REPETER' BLOC2 ('DIME' lxx);
x = 'EXTRAIRE' lxx &BLOC2;
'SI' (x < 0) ;
lroan = lroan 'ET' ('PROG' rog);
lpan = lpan 'ET' ('PROG' pg);
luan = luan 'ET' ('PROG' ug);
lhtan = lhtan 'ET' ('PROG' hg);
'SINON' ;
lroan = lroan 'ET' ('PROG' rod);
lpan = lpan 'ET' ('PROG' pd);
luan = luan 'ET' ('PROG' ud);
lhtan = lhtan 'ET' ('PROG' hg);
'FINSI' ;
'FIN' BLOC2;
evroa = 'EVOL' 'MANU' 'x' lxx 'ro' lroan;
evpa = 'EVOL' 'MANU' 'x' lxx 'p' lpan;
evua = 'EVOL' 'MANU' 'x' lxx 'ua' luan;
evha = 'EVOL' 'MANU' 'x' lxx 'ht' lhtan;
'SI' GRAPH ;
TAB1=TABLE;
TAB1.'TITRE'= TABLE ;
TAB1.1='MARQ TRIB REGU';
TAB1.'TITRE' . 1 = MOT 'Numerique' ;
TAB1.2='MARQ CROI REGU';
TAB1.'TITRE' . 2 = MOT 'Exacte' ;
'DESSIN' (evro 'ET' evroa) 'TITRE' tro 'XBOR' x0 x1
'LEGE' TAB1;
'DESSIN' (evp 'ET' evpa) 'TITRE' tp 'XBOR' x0 x1
'LEGE' TAB1;
'DESSIN' (evu 'ET' evua) 'TITRE' tu 'XBOR' x0 x1
'LEGE' TAB1;
'DESSIN' (evht 'ET' evha) 'TITRE' tht 'XBOR' x0 x1
'LEGE' TAB1;
'FINSI' ;
*
**** Tests
*
* Convergence
ERCON = 'EXTRAIRE' ('DIME' (RV . 'INCO' . 'ER')) (RV . 'INCO' . 'ER') ;
'SI' (ERCON > -3) ;
'MESSAGE' 'Probleme de convergence' ;
'ERREUR' 5;
'FINSI' ;
'SI' GRAPH ;
EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
(RV.INCO.'ER') ;
'DESSIN' EVOL4 'XBOR' 0. (2. '*' niter) 'YBOR' -12.0 0.0
'TITRE' ('CHAINE' METO ' IE =' IE ' IT =' IT) ;
'FINSI' ;
ERRO = 'ABS' (LRO '-' LROAN) ;
ERP = 'ABS' (LP '-' LPAN) ;
ERU = 'ABS' (LU '-' LUAN) ;
ERH = 'ABS' (LHT '-' LHTAN) ;
L1RO = 0.0 ;
L1P = 0.0 ;
L1U = 0.0 ;
L1H = 0.0 ;
NDIM = 'DIME' erro ;
'REPETER' BL1 NDIM ;
L1RO = L1RO '+' ('EXTRAIRE' &BL1 ERRO) ;
L1P = L1P '+' ('EXTRAIRE' &BL1 ERP) ;
L1U = L1U '+' ('EXTRAIRE' &BL1 ERU) ;
L1H = L1H '+' ('EXTRAIRE' &BL1 ERH) ;
'FIN' BL1 ;
L1RO = L1RO '/' (NDIM * ('MAXIMUM' lroan)) ;
L1P = L1P '/' (NDIM * ('MAXIMUM' lpan)) ;
L1U = L1U '/' (NDIM * ('MAXIMUM' luan 'ABS' )) ;
L1H = L1H '/' (NDIM * ('MAXIMUM' lhtan)) ;
ERRMAX = 'MAXIMUM' ('PROG' L1RO L1P L1U L1H) ;
'SI' (ERRMAX > 1.0D-4);
'ERREUR' 5;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer