* fichier : pret_ther4.dgibi ************************************************************************ * Section : Fluides Darcy ************************************************************************ ********************************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la solution des **** **** Equations d'Euler pour un gaz parfait. **** **** **** **** OPERATEUR PRET **** **** Operateur qui 'recontruit les variables primitives aux faces **** **** Cas gaz "thermally perfect" multi-especes **** **** Interieur et Mur **** **** Splitting des scalaires passifs **** **** **** **** A. BECCANTINI DRN/DMT/SEMT/TTMF FEVRIER 2000 **** ********************************************************************** 'OPTION' 'DIME' 2 ; 'OPTION' 'ELEM' QUA4 ; 'OPTION' 'ECHO' 0 ; 'OPTION' 'TRAC' 'X' ; * *** GRAPH * * GRAPH = VRAI ; GRAPH = FAUX ; *************************** *** PROCEDURE RUO1 **** *************************** * *** Precedure pour la rotation d'un tenseur de premier * ordre * 'DEBPROC' RUO1 ; 'ARGUMENT' ALPHA*'FLOTTANT' UN*'FLOTTANT' UT*'FLOTTANT'; SINA = 'SIN' ALPHA ; COSA = 'COS' ALPHA ; UX = (UN * COSA ) '-' (UT * SINA) ; UY = (UN * SINA ) '+' (UT * COSA) ; 'FINPROC' UX UY; *************************** *** PROCEDURE RUO2 **** *************************** * *** Precedure pour la rotation d'un tenseur de deuxieme * ordre * 'DEBPROC' RUO2 ; 'ARGUMENT' ALPHA*'FLOTTANT' UNN*'FLOTTANT' UNT*'FLOTTANT' UTN*'FLOTTANT' UTT*'FLOTTANT'; * **** (n,t) -> (x,y) * * n = CA x '+' SA y ; * t = -SA x '+' CA y ; * * UNT = DUN/DT * SA = 'SIN' ALPHA ; CA = 'COS' ALPHA ; CA2 = CA * CA ; CASA = CA * SA ; SA2 = SA * SA; UXX = (CA2 * UNN) '-' (CASA * (UNT '+' UTN)) '+' (SA2 *UTT) ; UYX = (CASA * (UNN '-' UTT)) '-' (SA2 * UNT) '+' (CA2 * UTN ) ; UXY = (CASA * (UNN '-' UTT)) '+' (CA2 * UNT) '-' (SA2 * UTN ) ; UYY = (SA2 * UNN) '+' (CASA * (UNT '+' UTN)) '+' (CA2 *UTT) ; 'FINPROC' UXX UXY UYX UYY; ******************************************************* **** Cas multi-especes : la table proprieté de gaz **** ******************************************************* PGAZ = 'TABLE' ; PGAZ2= 'TABLE' ; * **** Especes qui sont dans les equations d'Euler * PGAZ . 'ESPEULE' = 'MOTS' 'H2 ' 'N2 ' ; PGAZ2. 'ESPEULE' = 'MOTS' 'H2 ' 'N2 ' ; * **** Espece qui n'y est pas * PGAZ . 'ESPNEULE' = 'O2 '; PGAZ2. 'ESPNEULE' = 'O2 '; * PGAZ . 'H2 ' = 'TABLE' ; PGAZ . 'N2 ' = 'TABLE' ; PGAZ . 'O2 ' = 'TABLE' ; PGAZ2. 'H2 ' = 'TABLE' ; PGAZ2. 'N2 ' = 'TABLE' ; PGAZ2. 'O2 ' = 'TABLE' ; * **** R (J/Kg/K) * PGAZ . 'H2 ' . 'R' = 4130.0 ; PGAZ . 'N2 ' . 'R' = 296.8 ; PGAZ . 'O2 ' . 'R' = 259.8 ; PGAZ2. 'H2 ' . 'R' = 4130.0 ; PGAZ2. 'N2 ' . 'R' = 296.8 ; PGAZ2. 'O2 ' . 'R' = 259.8 ; * **** Regressions polynomials * PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.0 0.0 0.0 0.0 ; PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 0.0 0.0 0.0 0. ; PGAZ . 'N2 ' . 'A' = 'PROG' 652.940766 0.0 0.0 0.0 0.0 ; PGAZ2. 'H2 ' . 'A' = 'PROG' 9834.91866 0.0 0.0 0.0 0.0 ; PGAZ2. 'O2 ' . 'A' = 'PROG' 575.012333 0.0 0.0 0.0 0. ; PGAZ2. 'N2 ' . 'A' = 'PROG' 652.940766 0.0 0.0 0.0 0.0 ; * **** "Enthalpies" (ou energies) de formations a OK (J/Kg) * Note: ce sont des entites numeriques * h_i = h_i(T0) '-' \int_0^{T0} cp_i(x) dx = * h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0) * PGAZ . 'H2 ' . 'H0K' = -4.195D6 ; PGAZ . 'O2 ' . 'H0K' = -2.634D5 ; PGAZ . 'N2 ' . 'H0K' = -2.953D5 ; PGAZ2. 'H2 ' . 'H0K' = -4.195D6 ; PGAZ2. 'O2 ' . 'H0K' = -2.634D5 ; PGAZ2. 'N2 ' . 'H0K' = -2.953D5 ; PGAZ2 . 'SCALPASS' = 'MOTS' 'TITI' 'TUTU' 'TATA' ; *************************** ***** DOMAINE SPATIAL **** *************************** * **** Deux carre * A1 = 0.0D0 0.0D0 ; A2 = 1.0D0 0.0D0 ; A3 = 2.0D0 0.0D0 ; A4 = 2.0D0 1.0D0 ; A5 = 1.0D0 1.0D0 ; A6 = 0.0D0 1.0D0 ; L12 = A1 'DROIT' 1 A2 ; L23 = A2 'DROIT' 1 A3 ; L34 = A3 'DROIT' 1 A4 ; L45 = A4 'DROIT' 1 A5 ; L56 = A5 'DROIT' 1 A6 ; L61 = A6 'DROIT' 1 A1 ; L25 = A2 'DROIT' 1 A5 ; DOM10 = 'DALL' L12 L25 L56 L61 'PLANE'; DOM20 = 'DALL' L23 L34 L45 ('INVERSE' L25) 'PLANE'; * **** P10: Point face entre le deux carre, ou on fait les controlles **** P20: Point face sur le "mur" ** P10 = 1.0 0.5 ; P20 = 2.0 0.5 ; DOM1 = DOM10 ; DOM2 = DOM20 ; DOMTOT = DOM1 ET DOM2; 'ELIMINATION' (DOM1 ET DOM2) 1D-6; DOMTOT = DOM1 ET DOM2; $DOMTOT = 'MODELISER' DOMTOT 'EULER'; $DOM1 = 'MODELISER' DOM1 'EULER'; $DOM2 = 'MODELISER' DOM2 'EULER'; TDOMTOT = 'DOMA' $DOMTOT 'VF'; TDOM1 = 'DOMA' $DOM1 'VF'; TDOM2 = 'DOMA' $DOM2 'VF'; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMTOT = TDOMTOT . 'QUAF' ; 'ELIMINATION' (MDOMTOT ET MDOM1) 0.0001 ; 'ELIMINATION' (MDOMTOT ET MDOM2) 0.0001 ; 'SI' GRAPH; 'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL') 'ET' P10 'ET' P20) 'TITRE' 'Domaine et FACEL'; 'FINSI' ; * **** TEST1: si les limiteurs sont nuls ou les gradients son nuls * le deuxieme ordre en espace 'et premier ordre * en temps degenere en premier ordre en espace * RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.0D0; PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 2.0D0; VIT = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (4.0 5.0); GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.5D0; * **** Limiteurs nuls, gradients non-nuls * GRADR = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (1.0D0 1.0D0); GRADP = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (1.0D0 1.0D0); GRADVX = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (1.0D0 1.0D0); GRADVY = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (1.0D0 1.0D0); GRADV = GRADVX 'ET' GRADVY ; ALR = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.0D0 ; ALP = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.0D0 ; ALV = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (0.0 0.0D0) ; * YN YG = 'PROG' 0.2 0.1 ; YD = 'PROG' 0.0 0.13 ; YN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' ('MOTS' 'H2 ' 'N2 ' ) (('EXTRAIRE' 1 YG) ('EXTRAIRE' 2 YG)) ; YN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' ('MOTS' 'H2 ' 'N2 ') (('EXTRAIRE' 1 YD) ('EXTRAIRE' 2 YD)) ; YN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' ('MOTS' 'H2 ' 'N2 ') (YN1 'ET' YN2) ; GRADYH2 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (1.0D0 1.0D0); GRADYN2 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (1.0D0 1.0D0); GRADYN = GRADYH2 'ET' GRADYN2 ; ALYN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (0.0 0.0) ; * scag = 11.342 ; scad = 34.2; SCAC11 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'TITI' scag ; SCAC12 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'TITI' scad ; SCAC1 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'TITI' (SCAC11 'ET' SCAC12) ; SCAC21 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'TUTU' scag ; SCAC22 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'TUTU' scad ; SCAC2 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'TUTU' (SCAC21 'ET' SCAC22) ; SCAC31 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'TATA' scag ; SCAC32 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'TATA' scad ; SCAC3 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'TATA' (SCAC31 'ET' SCAC32) ; SCAC = SCAC1 'ET' SCAC2 'ET' SCAC3 ; GRADS1 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (1.0D0 1.0D0); GRADS2 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (1.0D0 1.0D0); GRADS3 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P3DX' 'P3DY' (1.0D0 1.0D0); GRADS = GRADS1 'ET' GRADS2 'ET' GRADS3 ; ALS1 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.0 ; ALS2 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P2' 0.0 ; ALS3 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P3' 0.0 ; ALS = ALS1 'ET' ALS2 'ET' ALS3 ; * ORDESP = 1; ORDTEM = 1; ROF1 VITF1 PF1 YF1 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN VIT PN YN ; ROF11 VITF11 PF11 YF11 SCAF11 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN VIT PN YN SCAC ; ORDESP = 2; ORDTEM = 1; ROF VITF PF YF = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN; ROF22 VITF22 PF22 YF22 SCAF22 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN SCAC GRADS ALS; ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ROF '-' ROF1) 'ABS') ('MAXIMUM' (ROF '-' ROF11) 'ABS') ('MAXIMUM' (ROF '-' ROF22) 'ABS') ('MAXIMUM' (PF '-' PF1) 'ABS') ('MAXIMUM' (PF '-' PF11) 'ABS') ('MAXIMUM' (PF '-' PF22) 'ABS') ('MAXIMUM' (VITF '-' VITF1) 'ABS') ('MAXIMUM'(VITF '-' VITF11)'ABS') ('MAXIMUM' (VITF '-' VITF22) 'ABS') ) ; ERRO = 'MAXIMUM' ('PROG' ERRO ('MAXIMUM' (YF '-' YF1) 'ABS') ('MAXIMUM' (YF '-' YF11) 'ABS') ('MAXIMUM' (YF '-' YF22) 'ABS') ('MAXIMUM' (SCAF11 '-' SCAF22) 'ABS')); 'SI' (ERRO > 1.0D-12); 'ERREUR' 5; 'FINSI' ; * **** Limiteurs non-nuls, gradients nuls * GRADR = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (0.0D0 0.0D0); GRADP = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (0.0D0 0.0D0); GRADVX = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (0.0D0 0.0D0); GRADVY = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (0.0D0 0.0D0); GRADV = GRADVX 'ET' GRADVY ; ALR = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.4D0 ; ALP = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 0.8D0 ; ALV = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (0.9 0.9D0) ; GRADYN = ('KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (0.0D0 0.0D0)) 'ET' ('KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (0.0D0 0.0D0)); ALYN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2 ' (1.0 1.0) ; GRADS1 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (0.0D0 0.0D0); GRADS2 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (0.0D0 0.0D0); GRADS3 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P3DX' 'P3DY' (0.0D0 0.0D0); GRADS = GRADS1 'ET' GRADS2 'ET' GRADS3 ; ALS1 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' 1.0 ; ALS2 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P2' 1.0 ; ALS3 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P3' 1.0 ; ALS = ALS1 'ET' ALS2 'ET' ALS3 ; * ORDESP = 1; ORDTEM = 1; ROF1 VITF1 PF1 YF1 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN VIT PN YN ; ROF11 VITF11 PF11 YF11 SCAF11 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN VIT PN YN SCAC ; ORDESP = 2; ORDTEM = 1; ROF VITF PF YF = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN; ROF22 VITF22 PF22 YF22 SCAF22 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN SCAC GRADS ALS; ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ROF '-' ROF1) 'ABS') ('MAXIMUM' (ROF '-' ROF11) 'ABS') ('MAXIMUM' (ROF '-' ROF22) 'ABS') ('MAXIMUM' (PF '-' PF1) 'ABS') ('MAXIMUM' (PF '-' PF11) 'ABS') ('MAXIMUM' (PF '-' PF22) 'ABS') ('MAXIMUM' (VITF '-' VITF1) 'ABS') ('MAXIMUM'(VITF '-' VITF11)'ABS') ('MAXIMUM' (VITF '-' VITF22) 'ABS') ) ; ERRO = 'MAXIMUM' ('PROG' ERRO ('MAXIMUM' (YF '-' YF1) 'ABS') ('MAXIMUM' (YF '-' YF11) 'ABS') ('MAXIMUM' (YF '-' YF22) 'ABS') ('MAXIMUM' (SCAF11 '-' SCAF22) 'ABS')); 'SI' (ERRO > 1.0D-12); 'ERREUR' 5; 'FINSI' ; * **** TEST2 : Deuxieme ordre en espace, deuxieme ordre en * temps: invariance par rotation * * Il faut remarquer que: * densite, pression sont des tenseurs d'ordre 0, * donc leurs gradients sont des tenseurs d'ordre 1 * la vitesse est un tenseur d'ordre 1, * donc sont gradient est un tenseur d'ordre 2 * **** On considere un repere solidal avec la surface, (n,t), * et un repere (x,y) * On consider une probleme 1D dans le repaire (n,t) * * Dans le cas de transort de scalaires passifs * (PGAZ2 . 'SCALPASS' = 'MOTS' 'TITI' 'TUTU' 'TATA') * on impose TITI = N2, TOTO = H2, TATA = N2 * *** Etats gauche et droite * gamg = 1.4D0; rog = 1.00 ; grong = -1.0D0 ; grotg = 0.0D0; alrog = 0.1 ; pg = 2.0D0 ; gpng = -3.0D0 ; gptg = 0.0D0; alpg = 0.2 ; ung = 3.0D0 ; gunng = -5.0D0 ; guntg = 0.0D0 ; utg = 0.0D0 ; gutng = 0.0D0 ; guttg = 0.0D0; alug = 0.3 ; yh2g = 0.1 ; gyh2ng = -0.0111 ; gyh2tg = 0.0 ; alyh2g = 0.123 ; yn2g = 0.34 ; gyn2ng = -0.0311 ; gyn2tg = 0.0 ; alyn2g = 0.323 ; c2g = gamg '*' pg '/' rog ; gamd = 1.8 ; rod = 1.50 ; grond = -1.5D0 ; grotd = 0.0D0 ; alrod = 0.25 ; pd = 2.5D0 ; gpnd = -3.5D0 ; gptd = 0.0D0 ; alpd = 0.3 ; und = 3.5D0 ; gunnd = -5.5D0 ; guntd = 0.0D0 ; utd = 0.0D0 ; gutnd = 0.0D0 ; guttd = 0.0D0; alud = 0.15 ; yh2d = 0.1131 ; gyh2nd = -0.03111 ; gyh2td = 0.0 ; alyh2d = 0.4123 ; yn2d = 0.134 ; gyn2nd = -0.01311 ; gyn2td = 0.0 ; alyn2d = 0.1323 ; c2d = gamd '*' pd '/' rod ; deltat = 1.0 ; * Encrement predictive drog = -1.0 '*' ((ung '*' grong '*' alrog) '+' (rog '*' gunng '*' alug)) ; drog = drog '*' deltat ; drod = -1.0 '*' ((und '*' grond '*' alrod) '+' (rod '*' gunnd '*' alud)) ; drod = drod '*' deltat ; dung = -1.0 '*' ((ung '*' gunng '*' alug) '+' (alpg '*' gpng '/' rog)) ; dung = dung '*' deltat ; dund = -1.0 '*' ((und '*' gunnd '*' alud) '+' (alpd '*' gpnd '/' rod)) ; dund = dund '*' deltat ; dpg = -1.0 '*' ((rog '*' c2g '*' gunng '*' alug) '+' (ung '*' gpng '*' alpg)) ; dpg = dpg '*' deltat ; dpd = -1.0 '*' ((rod '*' c2d '*' gunnd '*' alud) '+' (und '*' gpnd '*' alpd)) ; dpd = dpd '*' deltat ; dyh2g = -1.0 '*' ((ung '*' gyh2ng '*' alyh2g)) ; dyh2g = dyh2g '*' deltat ; dyh2d = -1.0 '*' ((und '*' gyh2nd '*' alyh2d)) ; dyh2d = dyh2d '*' deltat ; dyn2g = -1.0 '*' ((ung '*' gyn2ng '*' alyn2g)) ; dyn2g = dyn2g '*' deltat ; dyn2d = -1.0 '*' ((und '*' gyn2nd '*' alyn2d)) ; dyn2d = dyn2d '*' deltat ; * Encrement corrective; dxg = -dxd = 0.5 (point P1) dsrog = 0.5 '*' grong '*' alrog ; dsrod = -0.5 '*' grond '*' alrod ; dspg = 0.5 '*' gpng '*' alpg ; dspd = -0.5 '*' gpnd '*' alpd ; dsung = 0.5 '*' gunng '*' alug ; dsund = -0.5 '*' gunnd '*' alud ; dsyh2g = 0.5 '*' gyh2ng '*' alyh2g ; dsyh2d = -0.5 '*' gyh2nd '*' alyh2d ; dsyn2g = 0.5 '*' gyn2ng '*' alyn2g ; dsyn2d = -0.5 '*' gyn2nd '*' alyn2d ; * Encrement corrective; dxd = 0.5 (point P2) dsrod2 = -1.0 '*' dsrod ; dspd2 = -1.0 '*' dspd ; dsund2 = -1.0 '*' dsund ; dsyh2d2 = -1.0 '*' dsyh2d ; dsyn2d2 = -1.0 '*' dsyn2d ; * Les etats (point P1) etrog = rog '+' drog '+' dsrog ; etrod = rod '+' drod '+' dsrod ; etung = ung '+' dung '+' dsung ; etund = und '+' dund '+' dsund ; etpg = pg '+' dpg '+' dspg ; etpd = pd '+' dpd '+' dspd ; etyh2g = yh2g '+' dyh2g '+' dsyh2g ; etyh2d = yh2d '+' dyh2d '+' dsyh2d ; etyn2g = yn2g '+' dyn2g '+' dsyn2g ; etyn2d = yn2d '+' dyn2d '+' dsyn2d ; * Les etats (point P2) etrod2 = rod '+' drod '+' dsrod2 ; etund2 = und '+' dund '+' dsund2 ; etpd2 = pd '+' dpd '+' dspd2 ; etyh2d2 = yh2d '+' dyh2d '+' dsyh2d2 ; etyn2d2 = yn2d '+' dyn2d '+' dsyn2d2 ; * *** Etats a gauche et a droite du point face P1 * * ro, un, ut, p, yh2, yn2 ETATG = 'PROG' etrog etung 0.0 etpg etyh2g etyn2g ; ETATD = 'PROG' etrod etund 0.0 etpd etyh2d etyn2d ; * *** Etats a gauche et a droite du point face P1 * * ro, un, ut, p, yh2, yn2 ETATG2 = 'PROG' etrod2 etund2 0.0 etpd2 etyh2d2 etyn2d2 ; ETATD2 = 'PROG' etrod2 (-1.0 '*'etund2) 0.0 etpd2 etyh2d2 etyn2d2 ; * *** Rotation * DANGLE = 85; ANGLE = -1.0 '*' DANGLE ; ORIG = 0.0D0 0.0D0; 'REPETER' BL1 5; ANGLE = ANGLE '+' DANGLE ; uxg uyg = RUO1 ANGLE UNG UTG ; uxd uyd = RUO1 ANGLE UND UTD ; groxg groyg = RUO1 ANGLE grong grotg; groxd groyd = RUO1 ANGLE grond grotd; gpxg gpyg = RUO1 ANGLE gpng gptg; gpxd gpyd = RUO1 ANGLE gpnd gptd; guxxg guxyg guyxg guyyg = RUO2 ANGLE gunng guntg gutng guttg; guxxd guxyd guyxd guyyd = RUO2 ANGLE gunnd guntd gutnd guttd; gyh2xg gyh2yg = RUO1 ANGLE gyh2ng gyh2tg ; gyh2xd gyh2yd = RUO1 ANGLE gyh2nd gyh2td ; gyn2xg gyn2yg = RUO1 ANGLE gyn2ng gyn2tg ; gyn2xd gyn2yd = RUO1 ANGLE gyn2nd gyn2td ; 'MESSAGE' ; 'MESSAGE' (CHAIN 'Angle de rotation= ' ANGLE); 'MESSAGE' ; DOM1 = DOM10 'TOURNER' ANGLE ORIG; DOM2 = DOM20 'TOURNER' ANGLE ORIG; P1 = P10 'TOURNER' ANGLE ORIG; P2 = P20 'TOURNER' ANGLE ORIG; DOMTOT = DOM1 ET DOM2; 'ELIMINATION' (DOM1 ET DOM2) 1D-6; DOMTOT = DOM1 ET DOM2; $DOMTOT = 'MODELISER' DOMTOT 'EULER'; $DOM1 = 'MODELISER' DOM1 'EULER'; $DOM2 = 'MODELISER' DOM2 'EULER'; TDOMTOT = 'DOMA' $DOMTOT 'VF'; TDOM1 = 'DOMA' $DOM1 'VF'; TDOM2 = 'DOMA' $DOM2 'VF'; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMTOT = TDOMTOT . 'QUAF' ; 'ELIMINATION' (MDOMTOT ET MDOM1) 0.0001 ; 'ELIMINATION' (MDOMTOT ET MDOM2) 0.0001 ; 'SI' GRAPH; 'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET' ('DOMA' $DOMTOT 'FACEL') 'ET' P1 'ET' P2) 'TITRE' 'Domaine et FACEL'; 'FINSI' ; * **** Redefinition de P1 dans $DOMTOT . 'FACE' * P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1; P2 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P2; *********************** **** Les CHPOINTs **** *********************** GAM1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' gamg ; GAM2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' gamd ; GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (GAM1 'ET' GAM2) ; RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog; RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod; RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2); VIT1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (uxg uyg); VIT2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (uxd uyd); VIT = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (VIT1 'ET' VIT2); PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pg; PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pd; PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2); YN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'H2 ' 'N2 ' (yh2g yn2g) ; YN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'H2 ' 'N2 ' (yh2d yn2d) ; YN = YN1 'ET' YN2 ; * **** On impose les gradients et le limiteurs * ALR1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'P1' alrog ; ALR2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'P1' alrod ; ALR = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' (ALR1 'ET' ALR2) ; ALP1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'P1' alpg ; ALP2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'P1' alpd ; ALP = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'P1' (ALP1 'ET' ALP2) ; ALV1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (alug alug) ; ALV2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (alud alud) ; ALV = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (ALV1 'ET' ALV2) ; ALY1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (alyh2g alyn2g) ; ALY2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (alyh2d alyn2d) ; ALYN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (ALY1 'ET' ALY2) ; GRADR1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (groxg groyg) ; GRADR2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (groxd groyd) ; GRADR = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (GRADR1 'ET' GRADR2); GRADP1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (gpxg gpyg) ; GRADP2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (gpxd gpyd) ; GRADP = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (GRADP1 'ET' GRADP2); GRADVX1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (guxxg guxyg); GRADVX2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (guxxd guxyd); GRADVY1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (guyxg guyyg); GRADVY2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (guyxd guyyd); GRADVX = GRADVX1 'ET' GRADVX2 ; GRADVY = GRADVY1 'ET' GRADVY2 ; GRADV = GRADVX 'ET' GRADVY ; GRADYN = ('KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (gyh2xg gyh2yg )) 'ET' ('KCHT' $DOM1 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (gyn2xg gyn2yg )) 'ET' ('KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P1DX' 'P1DY' (gyh2xd gyh2yd )) 'ET' ('KCHT' $DOM2 'VECT' 'CENTRE' 'COMP' 'P2DX' 'P2DY' (gyn2xd gyn2yd )); *N.B: PGAZ2 . 'SCALPASS' = 'MOTS' 'TITI' 'TUTU' 'TATA' * PGAZ2 . 'ESPEULE' = 'MOTS' 'H2' 'N2' SCAC = 'EXCO' ('MOTS' 'N2' 'H2' 'N2') YN ('MOTS' 'TITI' 'TUTU' 'TATA') 'NATURE' 'DISCRET' ; GRADS = 'EXCO' ('MOTS' 'P2DX' 'P2DY' 'P1DX' 'P1DY' 'P2DX' 'P2DY') GRADYN ('MOTS' 'P1DX' 'P1DY' 'P2DX' 'P2DY' 'P3DX' 'P3DY' ) 'NATURE' 'DISCRET' ; ALS = 'EXCO' ('MOTS' 'P2' 'P1' 'P2') ALYN ('MOTS' 'P1' 'P2' 'P3') 'NATURE' 'DISCRET' ; ORDESP = 2; ORDTEM = 1; ROF VITF PF YF = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN ; ROF1 VITF1 PF1 YF1 SCAF1 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN SCAC GRADS ALS ; YF11 = 'EXCO' 'H2' YF1 'SCAL' ; YF12 = 'EXCO' 'N2' YF1 'SCAL' ; SCAF11 = 'EXCO' 'TITI' SCAF1 'SCAL' ; SCAF12 = 'EXCO' 'TUTU' SCAF1 'SCAL' ; SCAF13 = 'EXCO' 'TATA' SCAF1 'SCAL' ; ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ROF '-' ROF) 'ABS') ('MAXIMUM' (PF '-' PF1) 'ABS') ('MAXIMUM' (VITF '-' VITF1) 'ABS') ('MAXIMUM' (YF '-' YF1) 'ABS') ); ERRO = 'MAXIMUM' ('PROG' ERRO ('MAXIMUM' (YF11 '-' SCAF12) 'ABS') ('MAXIMUM' (YF12 '-' SCAF11) 'ABS') ('MAXIMUM' (YF12 '-' SCAF13) 'ABS')) ; 'SI' (ERRO > 1.0D-12); 'ERREUR' 5; 'FINSI' ; * ORDESP = 2; ORDTEM = 2; ROF VITF PF YF = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN GAMMA deltat ; ROF1 VITF1 PF1 YF1 SCAF1 = 'PRET' 'PERFTEMP' ORDESP ORDTEM $DOMTOT PGAZ2 RN GRADR ALR VIT GRADV ALV PN GRADP ALP YN GRADYN ALYN SCAC GRADS ALS GAMMA deltat ; YF11 = 'EXCO' 'H2' YF1 'SCAL' ; YF12 = 'EXCO' 'N2' YF1 'SCAL' ; SCAF11 = 'EXCO' 'TITI' SCAF1 'SCAL' ; SCAF12 = 'EXCO' 'TUTU' SCAF1 'SCAL' ; SCAF13 = 'EXCO' 'TATA' SCAF1 'SCAL' ; ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ROF '-' ROF) 'ABS') ('MAXIMUM' (PF '-' PF1) 'ABS') ('MAXIMUM' (VITF '-' VITF1) 'ABS') ('MAXIMUM' (YF '-' YF1) 'ABS') ); ERRO = 'MAXIMUM' ('PROG' ERRO ('MAXIMUM' (YF11 '-' SCAF12) 'ABS') ('MAXIMUM' (YF12 '-' SCAF11) 'ABS') ('MAXIMUM' (YF12 '-' SCAF13) 'ABS')) ; 'SI' (ERRO > 1.0D-12); 'ERREUR' 5; 'FINSI' ; ********************************************************* *** Control des etats sur la surface qui contient P1 **** ********************************************************* GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P1 ; GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P1 ; GEOP12 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P2 ; ROGEOP1 = 'REDU' ROF GEOP1; VGEOP1 = 'REDU' VITF GEOP1; PGEOP1 = 'REDU' PF GEOP1; REFGEOP1 = 'REDU' VITF GEOP2; YNGEOP1 = 'REDU' YF GEOP1 ; ROGEOP2 = 'REDU' ROF GEOP12 ; VGEOP2 = 'REDU' VITF GEOP12 ; PGEOP2 = 'REDU' PF GEOP12 ; YNGEOP2 = 'REDU' YF GEOP12 ; rofag = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 1; rofad = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 3; unfag = 'EXTRAIRE' VGEOP1 'UN ' 1 1 1; unfad = 'EXTRAIRE' VGEOP1 'UN ' 1 1 3; utfag = 'EXTRAIRE' VGEOP1 'UT ' 1 1 1; utfad = 'EXTRAIRE' VGEOP1 'UT ' 1 1 3; pfag = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1; pfad = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 3; yh2fag = 'EXTRAIRE' YNGEOP1 'H2 ' 1 1 1; yh2fad= 'EXTRAIRE' YNGEOP1 'H2 ' 1 1 3; yn2fag = 'EXTRAIRE' YNGEOP1 'N2 ' 1 1 1; yn2fad= 'EXTRAIRE' YNGEOP1 'N2 ' 1 1 3; rofag2 = 'EXTRAIRE' ROGEOP2 'SCAL' 1 1 1; rofad2 = 'EXTRAIRE' ROGEOP2 'SCAL' 1 1 3; unfag2 = 'EXTRAIRE' VGEOP2 'UN ' 1 1 1; unfad2 = 'EXTRAIRE' VGEOP2 'UN ' 1 1 3; utfag2 = 'EXTRAIRE' VGEOP2 'UT ' 1 1 1; utfad2 = 'EXTRAIRE' VGEOP2 'UT ' 1 1 3; pfag2 = 'EXTRAIRE' PGEOP2 'SCAL' 1 1 1; pfad2 = 'EXTRAIRE' PGEOP2 'SCAL' 1 1 3; yh2fag2 = 'EXTRAIRE' YNGEOP2 'H2 ' 1 1 1; yh2fad2 = 'EXTRAIRE' YNGEOP2 'H2 ' 1 1 3; yn2fag2 = 'EXTRAIRE' YNGEOP2 'N2 ' 1 1 1; yn2fad2 = 'EXTRAIRE' YNGEOP2 'N2 ' 1 1 3; * **** Orientation de la normal n de castem par raport a la * notre; t est par consequence (sur P1) * NCOS = 'EXTRAIRE' REFGEOP1 'NX' 1 1 1; NSIN = 'EXTRAIRE' REFGEOP1 'NY' 1 1 1; 'SI' (('ABS' NCOS) > ('ABS' NSIN)); ORIENT = ('COS' ANGLE) '/' NCOS; 'SINON'; ORIENT = ('SIN' ANGLE) '/' NSIN; 'FINSI' ; ORIENT = 'SIGN' ORIENT; * **** ORIENT = -1 -> Mon etat gauche est son etat droite * 'SI' (ORIENT > 0); ERRLIG = 'PROG' rofag (unfag '*' ORIENT) (utfag '*' ORIENT) pfag yh2fag yn2fag ; ERRLID = 'PROG' rofad (unfad '*' ORIENT) (utfad '*' ORIENT) pfad yh2fad yn2fad ; 'SINON' ; ERRLID = 'PROG' rofag (unfag '*' ORIENT) (utfag '*' ORIENT) pfag yh2fag yn2fag; ERRLIG = 'PROG' rofad (unfad '*' ORIENT) (utfad '*' ORIENT) pfad yh2fad yn2fad ; 'FINSI' ; ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ETATG '-' ERRLIG) 'ABS') ('MAXIMUM' (ETATD '-' ERRLID) 'ABS') ); 'SI' (ERRO > 1.0D-14) 'MESSAGE' 'Ordre en espace = 2'; 'MESSAGE' 'Ordre en temps = 2'; 'ERREUR' 5 ; 'FINSI' ; * **** Sir le mur (P2) * ERRLIG = 'PROG' rofag2 unfag2 utfag2 pfag2 yh2fag2 yn2fag2 ; ERRLID = 'PROG' rofad2 unfad2 utfad2 pfad2 yh2fad2 yn2fad2 ; ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ETATG2 '-' ERRLIG) 'ABS') ('MAXIMUM' (ETATD2 '-' ERRLID) 'ABS') ); 'SI' (ERRO > 1.0D-14) 'Message' 'Mur' ; 'MESSAGE' 'Ordre en espace = 2'; 'MESSAGE' 'Ordre en temps = 2'; 'ERREUR' 5 ; 'FINSI' ; 'FIN' BL1; 'FIN' ;