* fichier : pret_ther3.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" mono-espece ****
**** Differents cas tests ****
**** Interieur ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/TTMF DECEMBRE 1998 ****
**********************************************************************
'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 mono-especes : la table proprieté de gaz ****
*******************************************************
PGAZ = 'TABLE' ;
*
**** Especes qui sont dans les equations d'Euler
*
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'O2 ';
*
PGAZ . 'O2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'O2 ' . 'R' = 259.8 ;
*
**** Regressions polynomials
*
PGAZ . 'O2 ' . 'A' = 'PROG' 575.012333 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 . 'O2 ' . 'H0K' = -2.634D5 ;
***************************
***** 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';
*
*** Point face entre le deux carre, ou on fait les controlles
*
P10 = 1.0 0.5;
DOM1 = DOM10 ;
DOM2 = DOM20 ;
DOMTOT = DOM1 'ET' DOM2 ;
'ELIMINATION' DOMTOT 0.0001 ;
$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) '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) ;
*
ORDESP = 1;
ORDTEM = 1;
ROF1 VITF1 PF1 = 'PRET' 'PERFTEMP' ORDESP ORDTEM
$DOMTOT PGAZ RN VIT PN ;
ORDESP = 2;
ORDTEM = 1;
ROF VITF PF = 'PRET' 'PERFTEMP' ORDESP ORDTEM
$DOMTOT PGAZ RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP ;
ERRO = 'MAXIMUM' ('PROG'
('MAXIMUM' (ROF '-' ROF1) 'ABS')
('MAXIMUM' (PF '-' PF1) 'ABS')
('MAXIMUM' (VITF '-' VITF1) 'ABS') );
'SI' (ERRO > 1.0D-15);
'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) ;
ORDESP = 1;
ORDTEM = 1;
ROF1 VITF1 PF1 = 'PRET' 'PERFTEMP' ORDESP ORDTEM
$DOMTOT PGAZ RN VIT PN ;
ORDESP = 2;
ORDTEM = 1;
ROF VITF PF = 'PRET' 'PERFTEMP' ORDESP ORDTEM
$DOMTOT PGAZ RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP ;
ERRO = 'MAXIMUM' ('PROG'
('MAXIMUM' (ROF '-' ROF1) 'ABS')
('MAXIMUM' (PF '-' PF1) 'ABS')
('MAXIMUM' (VITF '-' VITF1) 'ABS')
) ;
'SI' (ERRO > 1.0D-15);
'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)
*
*** 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 ;
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 ;
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 ;
* Encrement corrective; dxg = -dxd = 0.5
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 ;
* Les etats
etrog = rog '+' drog '+' dsrog ;
etrod = rod '+' drod '+' dsrod ;
etung = ung '+' dung '+' dsung ;
etund = und '+' dund '+' dsund ;
etpg = pg '+' dpg '+' dspg ;
etpd = pd '+' dpd '+' dspd ;
*
*** Etats a gauche et a droite du point face P10
*
* ro, un, ut, p
ETATG = 'PROG' etrog etung 0.0 etpg ;
ETATD = 'PROG' etrod etund 0.0 etpd ;
*
*** 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;
'MESSAGE' ;
'MESSAGE' (CHAIN 'Angle de rotation= ' ANGLE);
'MESSAGE' ;
DOM1 = DOM10 'TOURNER' ANGLE ORIG;
DOM2 = DOM20 'TOURNER' ANGLE ORIG;
P1 = P10 '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) 'TITRE' 'Domaine et FACEL';
'FINSI' ;
*
**** Redefinition de P1 dans $DOMTOT . 'FACE'
*
P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
***********************
**** 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);
*
**** 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) ;
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 ;
ORDESP = 2;
ORDTEM = 2;
ROF VITF PF = 'PRET' 'PERFTEMP' ORDESP ORDTEM
$DOMTOT PGAZ RN GRADR ALR
VIT GRADV ALV
PN GRADP ALP
GAMMA deltat ;
*********************************************************
*** 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;
ROGEOP1 = 'REDU' ROF GEOP1;
VGEOP1 = 'REDU' VITF GEOP1;
PGEOP1 = 'REDU' PF GEOP1;
REFGEOP1 = 'REDU' VITF GEOP2;
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;
*
**** Orientation de la normal n de castem par raport a la
* notre; t est par consequence
*
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 ;
ERRLID = 'PROG' rofad (unfad '*' ORIENT)
(utfad '*' ORIENT) pfad ;
'SINON' ;
ERRLID = 'PROG' rofag (unfag '*' ORIENT)
(utfag '*' ORIENT) pfag ;
ERRLIG = 'PROG' rofad (unfad '*' ORIENT)
(utfad '*' ORIENT) pfad ;
'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 = 1';
'ERREUR' 5 ;
'FINSI' ;
'FIN' BL1;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer