* fichier : pret_ther4.dgibi
**********************************************************************
**** 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 = 'SIGNE' 'FLOTTANT' 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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer