* 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 ' 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 ' 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' ;