* fichier : pret_ther.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 mono-espece "thermally perfect" ****
**** Premier ordre en espace, premier ordre en temps ****
**** Enterieur et murs ****
**** ****
**** 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;
***************************
***** 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 ;
*
*** Point face sur le mur (a droite)
*
P20 = 2.0 0.5 ;
*
*** Etats a gauche et a droite du point face P10
*
* ro, un, ut, p
rog = 1.100 ;
pg = 2.13D0 ;
ung = 1.17D0 ;
utg = 1.14D0 ;
sca1g = 1.76 ;
sca2g = 311.76 ;
rod = 1.50 ;
pd = 2.5D0 ;
und = 3.5D0 ;
utd = 4.5D0 ;
sca1d = 217.6 ;
sca2d = 110.0 ;
ETATG = 'PROG' rog ung utg pg sca1g sca2g ;
ETATD = 'PROG' rod und utd pd sca1d sca2d ;
**** Cas monoespece : la table proprieté de gaz
*
* PGAZ2 avec 'SCALPASS'
PGAZ = 'TABLE' ;
PGAZ2 = 'TABLE' ;
*
**** Especes qui sont dans les equations d'Euler (rien)
*
*
**** Espece qui n'y est pas
*
PGAZ . 'ESPNEULE' = 'H2 ';
PGAZ2 . 'ESPNEULE' = 'H2 ';
*
PGAZ . 'H2 ' = 'TABLE' ;
PGAZ2 . 'H2 ' = 'TABLE' ;
*
**** R (J/Kg/K)
*
PGAZ . 'H2 ' . 'R' = 4130.0 ;
PGAZ2 . 'H2 ' . 'R' = 4130.0 ;
*
**** Regressions polynomials
*
PGAZ . 'H2 ' . 'A' = 'PROG' 9834.91866 0.54273926 0.000862203836
-2.37281455E-07 1.84701105E-11 ;
PGAZ2 . '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 ;
PGAZ2 . 'H2 ' . 'H0K' = -4.195D6 ;
PGAZ2 . 'SCALPASS' = 'MOTS' 'TITU' 'TOTO' ;
**** Test1: rotation
ANGLE = 7.0D0;
DANGLE = 85;
ORIG = 0.0D0 0.0D0;
'REPETER' BL1 5;
ANGLE = ANGLE '+' DANGLE ;
uxg uyg = RUO1 ANGLE UNG UTG ;
uxd uyd = RUO1 ANGLE UND UTD ;
'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' P10) 'TITRE' 'Domaine et FACEL';
'FINSI' ;
**** La densité, la pression, la vitesse
RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog ;
RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod ;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2) ;
PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pg ;
PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pd ;
PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2) ;
VN1X = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'UX ' uxg ;
VN1Y = 'KCHT' $DOM1 'SCAL' 'CENTRE' 'COMP' 'UY ' uyg ;
VN2X = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'UX ' uxd ;
VN2Y = 'KCHT' $DOM2 'SCAL' 'CENTRE' 'COMP' 'UY ' uyd ;
VN = 'KCHT' $DOMTOT 'VECT' 'CENTRE'
(VN1X 'ET' VN1Y 'ET' VN2X 'ET' VN2Y);
SCAN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' 'COMP'
('MOTS' 'TITU' 'TOTO') (sca1g sca2g) ;
SCAN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' 'COMP'
('MOTS' 'TITU' 'TOTO') (sca1d sca2d) ;
SCAN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' 'COMP'
('MOTS' 'TITU' 'TOTO') (SCAN1 'ET' SCAN2) ;
*
**** Cas monoespece: dans la table PGAZ il n'y a pas de
* PGAZ . 'ESPEULE'
ROF VITF PF = 'PRET' 'PERFTEMP' 1 1 $DOMTOT PGAZ RN VN PN ;
ROF2 VITF2 PF2 SCAF = 'PRET' 'PERFTEMP'
1 1 $DOMTOT PGAZ2 RN VN PN SCAN ;
ERRO = 'MAXIMUM' (ROF '-' ROF2) 'ABS' ;
ERVI = 'MAXIMUM' (VITF '-' VITF2) 'ABS' ;
ERP = 'MAXIMUM' (PF '-' PF2) 'ABS' ;
'SI' (('MAXIMUM' ('PROG' ERRO ERVI ERP )) > 1.0D-10) ;
'MESSAGE' 'Erreur trop importante' ;
'ERREUR' 5 ;
'FINSI' ;
*
**** Redefinition de P1, P2 dans $DOMTOT . 'FACE'
*
P1 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1;
P2 = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P2;
*********************************************************
*** 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;
SGEOP1 = 'REDU' SCAF GEOP1 ;
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;
s1fag = 'EXTRAIRE' SGEOP1 'TITU' 1 1 1;
s1fad = 'EXTRAIRE' SGEOP1 'TITU' 1 1 3;
s2fag = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 1;
s2fad = 'EXTRAIRE' SGEOP1 'TOTO' 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 s1fag s2fag ;
ERRLID = 'PROG' rofad (unfad '*' ORIENT)
(utfad '*' ORIENT) pfad s1fad s2fad ;
'SINON' ;
ERRLID = 'PROG' rofag (unfag '*' ORIENT)
(utfag '*' ORIENT) pfag s1fag s2fag ;
ERRLIG = 'PROG' rofad (unfad '*' ORIENT)
(utfad '*' ORIENT) pfad s1fad s2fad ;
'FINSI' ;
ERRO = 'MAXIMUM' ('PROG' ('MAXIMUM' (ETATG '-' ERRLIG) 'ABS')
('MAXIMUM' (ETATD '-' ERRLID) 'ABS'));
'SI' (ERRO > 1.0D-12) ;
'MESSAGE' 'Interieur' ;
'MESSAGE' 'Ordre en espace = 1';
'MESSAGE' 'Ordre en temps = 1';
'ERREUR' 5 ;
'FINSI' ;
*********************************************************
*** Control des etats sur la surface qui contient P2 ****
*********************************************************
GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'APPUYE' 'LARGEMENT' P2;
ROGEOP1 = 'REDU' ROF GEOP1;
VGEOP1 = 'REDU' VITF GEOP1;
PGEOP1 = 'REDU' PF GEOP1;
REFGEOP1 = 'REDU' VITF GEOP2;
SGEOP1 = 'REDU' SCAF GEOP1 ;
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;
s1fag = 'EXTRAIRE' SGEOP1 'TITU' 1 1 1;
s1fad = 'EXTRAIRE' SGEOP1 'TITU' 1 1 3;
s2fag = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 1;
s2fad = 'EXTRAIRE' SGEOP1 'TOTO' 1 1 3;
*
* Par convention l'etat droite est l'etat mirroir
*
* Sur le mur: Gauche = etat reel
* Droit = etat mirroir
*
LOG1 = ('ABS' (rod '-' rofag)) > 1.0D-12 ;
LOG1 = 'OU' LOG1 (('ABS' (rod '-' rofad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (und '-' unfag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (und '+' unfad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (utd '-' utfag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (utd '-' utfad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (pd '-' pfag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (pd '-' pfad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca1d '-' s1fad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fad)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fag)) > 1.0D-12) ;
LOG1 = 'OU' LOG1 (('ABS' (sca2d '-' s2fad)) > 1.0D-12) ;
'SI' (LOG1) ;
'MESSAGE' 'Mur' ;
'MESSAGE' 'Ordre en espace = 1';
'MESSAGE' 'Ordre en temps = 1';
'ERREUR' 5 ;
'FINSI' ;
'FIN' BL1 ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer