* fichier : shearlayer.dgibi
*********************************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la solution des ****
**** Equations d'Euler pour un gaz parfait ****
**** OPERATEURS PRET, KONV ****
**** Gaz monoespece "calorically perfect" ****
**** HUSVL, HUSVLH, Godunov on "shear layer" ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/TTMF MARS 1998 ****
*********************************************************************
'OPTION' 'DIME' 2 ;
'OPTION' 'ELEM' QUA4 ;
'OPTION' 'ECHO' 0 ;
'OPTION' 'TRAC' 'PS';
GRAPH = FAUX ;
* GRAPH = VRAI ;
***************************
***** DOMAINE SPATIAL ****
***************************
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 ou on controlle les flux
*
P10 = 1.0 0.5;
*
*** Etats gauche "shear layer"
*
rogsl = 1.0;
ungsl = 1.0D0;
utgsl = 10.0;
pgsl = 1.;
gamgsl = 1.4D0;
*
*** retgsl
*
ecingsl = 0.5D0 '*' rogsl '*' ((ungsl '*' ungsl) '+' (utgsl '*' utgsl));
retgsl = (pgsl '/' (gamgsl '-' 1.0)) '+' ecingsl;
*
*** Etats gauche et droite "shear layer"
*
rodsl = 100.0;
undsl = 1.0;
utdsl = 1.;
pdsl = 1.;
gamdsl = 1.4D0;
*
*** retdsl
*
ecindsl = 0.5D0 '*' rodsl '*' ((undsl '*' undsl) '+' (utdsl '*' utdsl));
retdsl = (pdsl '/' (gamdsl '-' 1.0)) '+' ecindsl;
*
*** flux en (n,t),g
*
f1gsl = ungsl '*' rogsl;
f2gsl = (rogsl '*' (ungsl '*' ungsl)) '+' pgsl;
f3gsl = rogsl '*' ungsl '*' utgsl ;
f4gsl = ungsl * (retgsl '+' pgsl);
*
*** flux en (n,t),d
*
f1dsl = undsl '*' rodsl;
f2dsl = (rodsl '*' (undsl '*' undsl)) '+' pdsl;
f3dsl = rodsl '*' undsl '*' utdsl ;
f4dsl = undsl * (retdsl '+' pdsl);
*
*** N.B: f2gsl = f2dsl = pgsl (=pdsl)
* figsl = fidsl = 0.0 \forall i != 2
*
****************************************************
****************************************************
******** Boucle sur les angles *********
****************************************************
****************************************************
DANGLE = 360 '/' 7.15;
ANGLE = 11.3 ;
'REPETER' BLOC 8;
*
*** Rotation
*
ANGLE = ANGLE '+' DANGLE;
ORIG = 0.0D0 0.0D0;
MESSAGE;
MESSAGE (CHAIN 'Angle de rotation= ' ANGLE);
MESSAGE;
DOM1 = DOM10 'TOURNER' ANGLE ORIG;
DOM2 = DOM20 'TOURNER' ANGLE ORIG;
P1 = P10 'TOURNER' ANGLE ORIG;
'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 ****
***********************
uxgsl = (ungsl '*' ('COS' ANGLE)) '-' (utgsl '*' ('SIN' ANGLE));
uygsl = (ungsl '*' ('SIN' ANGLE)) '+' (utgsl '*' ('COS' ANGLE));
uxdsl = (undsl '*' ('COS' ANGLE)) '-' (utdsl '*' ('SIN' ANGLE));
uydsl = (undsl '*' ('SIN' ANGLE)) '+' (utdsl '*' ('COS' ANGLE));
GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gamgsl;
RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rogsl;
RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rodsl;;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
GN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' ((rogsl*uxgsl) (rogsl*uygsl));
GN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' ((rodsl*uxdsl) (rodsl*uydsl));
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN1 'ET' GN2);
PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pgsl;
PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pdsl;
PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
*
*** La vitesse
*
VITX = ('EXCO' 'UX ' GN) '/' RN;
VITY = ('EXCO' 'UY ' GN) '/' RN;
*
*** L'energie totale (ROET)
*
GM1 = GAMMA '-' ('KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.0D0);
ETHER = PN '/' GM1;
ECIN = (0.5D0 * ('KOPS' GN 'PSCA' GN )) '/' RN;
EN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (ECIN '+' ETHER);
***************************
**** L'operateur PRIM****
***************************
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN EN GAMMA;
***************************
**** L'operateur PRET****
***************************
ORDESP = 1;
ORDTEM = 1;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMA;
*********************************************************
*** Control des etats sur la surface qui contient P1 ****
*********************************************************
GEOP1 = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'CONTENANT' P1;
GEOP2 = ('DOMA' $DOMTOT 'FACE') 'ELEM' 'CONTENANT' P1;
ROGEOP1 = 'REDU' ROF GEOP1;
VGEOP1 = 'REDU' VITF GEOP1;
PGEOP1 = 'REDU' PF GEOP1;
GAMGEOP1 = 'REDU' GAMF GEOP1;
REFGEOP1 = 'REDU' VITF GEOP2;
rog = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 1;
rod = 'EXTRAIRE' ROGEOP1 'SCAL' 1 1 3;
ung = 'EXTRAIRE' VGEOP1 'UN ' 1 1 1;
und = 'EXTRAIRE' VGEOP1 'UN ' 1 1 3;
utg = 'EXTRAIRE' VGEOP1 'UT ' 1 1 1;
utd = 'EXTRAIRE' VGEOP1 'UT ' 1 1 3;
pg = 'EXTRAIRE' PGEOP1 'SCAL' 1 1 1;
pd = '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;
*
**** TEST 'PRET' 'ET' 'PRIM'
*
'SI' (ORIENT '>' 0);
*
*** MyGauche = CASTEM Gauche
*
ERRO = 1D-10 '*' ('MAXIMUM' ('PROG' rogsl rodsl));
ERRORO = 'MAXIMUM' ('PROG' ('ABS' (rog '-' rogsl))
('ABS' (rod '-' rodsl)));
LOGI1 = ERRORO < ERRO;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR RO ' ERRORO);
* 'MESSAGE' ;
ERRO = 1D-10 '*' pgsl ;
ERROP = 'MAXIMUM' ('PROG' ('ABS' (pg '-' pgsl))
('ABS' (pd '-' pdsl)));
LOGI2 = ERROP < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR P ' ERROP);
* 'MESSAGE' ;
ERROUN = 'MAXIMUM' ('PROG' ('ABS' (ung '-' ungsl))
('ABS' (und '-' undsl)));
ERRO = 1D-10 ;
LOGI2 = ERROUN < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR UN ' ERROUN);
* 'MESSAGE' ;
ERRO = 1D-10 '*' ('MAXIMUM' ('PROG' utgsl utdsl));
ERROUT = 'MAXIMUM' ('PROG' ('ABS' (utg '-' utgsl))
('ABS' (utd '-' utdsl)));
LOGI2 = ERROUT < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR UT ' ERROUT);
* 'MESSAGE' ;
'SINON' ;
*
*** MyGauche = CASTEM Droite
*
* ung = -undsl
* utg = -utdsl
*
ERRO = 1D-10 '*' ('MAXIMUM' ('PROG' rogsl rodsl));
ERRORO = 'MAXIMUM' ('PROG' ('ABS' (rog '-' rodsl))
('ABS' (rod '-' rogsl)));
LOGI1 = ERRORO < ERRO;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR RO ' ERRORO);
* 'MESSAGE' ;
ERRO = 1D-10 '*' pgsl ;
ERROP = 'MAXIMUM' ('PROG' ('ABS' (pg '-' pdsl))
('ABS' (pd '-' pgsl)));
LOGI2 = ERROP < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR P ' ERROP);
* 'MESSAGE' ;
ERROUN = 'MAXIMUM' ('PROG' ('ABS' (ung '+' undsl))
('ABS' (und '+' ungsl)));
ERRO = 1D-10 ;
LOGI2 = ERROUN < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR UN ' ERROUN);
* 'MESSAGE' ;
ERRO = 1D-10 '*' ('MAXIMUM' ('PROG' utgsl utdsl));
ERROUT = 'MAXIMUM' ('PROG' ('ABS' (utd '+' utgsl))
('ABS' (utg '+' utdsl)));
LOGI2 = ERROUT < ERRO;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR UT ' ERROUT);
* 'MESSAGE' ;
'FINSI';
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'Probleme en PRET ou en PRIM ' ;
'MESSAGE' ;
'ERREUR' 5;
'FINSI' ;
********************************
**** 'KONV', TEST HUSVL ****
********************************
METO = 'HUSVL' ;
TFLUX DT = 'KONV' 'VF' 'PERFMONO' 'FLUX' METO
$DOMTOT ROF VITF PF GAMF
('MOTS' 'RN' 'UX' 'UY' 'RETN') ;
*
**** Les flux aux interfaces sont dans de
* CHPOINT FACE
*
FLUX1 = 'EXCO' 'RN' TFLUX ;
FLUX2X = 'EXCO' 'UX' TFLUX ;
FLUX2Y = 'EXCO' 'UY' TFLUX ;
FLUX3 = 'EXCO' 'RETN' TFLUX ;
FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));
*
**** Test 'KONV'
*
f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
ERRO = 1D-10;
ERROF1 = 'ABS' ((f1 '*' ORIENT) '-' f1gsl);
LOGI1 = ERROF1 < ERRO;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F1 ' ERROF1);
* 'MESSAGE' ;
f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
ERRO = 1D-10 '*' f2gsl;
ERROF2 = 'ABS' ((f2 '*' ORIENT) '-' f2gsl);
LOGI2 = ERROF2 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F2 ' ERROF2);
* 'MESSAGE' ;
f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
ERRO = 1D-10;
ERROF3 = 'ABS' ((f3 '*' ORIENT) '-' f3gsl);
LOGI2 = ERROF3 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F3 ' ERROF3);
* 'MESSAGE' ;
f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
ERRO = 1D-10;
ERROF4 = 'ABS' ((f4 '*' ORIENT) '-' f4gsl);
LOGI2 = ERROF4 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F4 ' ERROF4);
* 'MESSAGE' ;
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'OPERATEUR KONV';
'MESSAGE' ('CHAINE' METO);
'ERREUR' 5;
'FINSI' ;
********************************
**** 'KONV', TEST HUSVLH ****
********************************
METO = 'HUSVLH' ;
TFLUX DT = 'KONV' 'VF' 'PERFMONO' 'FLUX' METO
$DOMTOT ROF VITF PF GAMF
('MOTS' 'RN' 'UX' 'UY' 'RETN') ;
*
**** Les flux aux interfaces sont dans de
* CHPOINT FACE
*
FLUX1 = 'EXCO' 'RN' TFLUX ;
FLUX2X = 'EXCO' 'UX' TFLUX ;
FLUX2Y = 'EXCO' 'UY' TFLUX ;
FLUX3 = 'EXCO' 'RETN' TFLUX ;
FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));
*
**** Test 'KONV'
*
f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
ERRO = 1D-10;
ERROF1 = 'ABS' ((f1 '*' ORIENT) '-' f1gsl);
LOGI1 = ERROF1 < ERRO;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F1 ' ERROF1);
* 'MESSAGE' ;
f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
ERRO = 1D-10 '*' f2gsl;
ERROF2 = 'ABS' ((f2 '*' ORIENT) '-' f2gsl);
LOGI2 = ERROF2 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F2 ' ERROF2);
* 'MESSAGE' ;
f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
ERRO = 1D-10;
ERROF3 = 'ABS' ((f3 '*' ORIENT) '-' f3gsl);
LOGI2 = ERROF3 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F3 ' ERROF3);
* 'MESSAGE' ;
f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
ERRO = 1D-10;
ERROF4 = 'ABS' ((f4 '*' ORIENT) '-' f4gsl);
LOGI2 = ERROF4 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F4 ' ERROF4);
* 'MESSAGE' ;
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'OPERATEUR KONV';
'MESSAGE' ('CHAINE' METO);
'ERREUR' 5;
'FINSI' ;
********************************
**** 'KONV', TEST GODUNOV ****
********************************
METO = 'GODUNOV' ;
TFLUX DT = 'KONV' 'VF' 'PERFMONO' 'FLUX' METO
$DOMTOT ROF VITF PF GAMF
('MOTS' 'RN' 'UX' 'UY' 'RETN') ;
*
**** Les flux aux interfaces sont dans de
* CHPOINT FACE
*
FLUX1 = 'EXCO' 'RN' TFLUX ;
FLUX2X = 'EXCO' 'UX' TFLUX ;
FLUX2Y = 'EXCO' 'UY' TFLUX ;
FLUX3 = 'EXCO' 'RETN' TFLUX ;
FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));
*
**** Test 'KONV'
*
f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
ERRO = 1D-10;
ERROF1 = 'ABS' ((f1 '*' ORIENT) '-' f1gsl);
LOGI1 = ERROF1 < ERRO;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F1 ' ERROF1);
* 'MESSAGE' ;
f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
ERRO = 1D-10 '*' f2gsl;
ERROF2 = 'ABS' ((f2 '*' ORIENT) '-' f2gsl);
LOGI2 = ERROF2 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F2 ' ERROF2);
* 'MESSAGE' ;
f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
ERRO = 1D-10;
ERROF3 = 'ABS' ((f3 '*' ORIENT) '-' f3gsl);
LOGI2 = ERROF3 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F3 ' ERROF3);
* 'MESSAGE' ;
f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
ERRO = 1D-10;
ERROF4 = 'ABS' ((f4 '*' ORIENT) '-' f4gsl);
LOGI2 = ERROF4 < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
* 'MESSAGE' ;
* 'MESSAGE' ('CHAINE' 'ERR F4 ' ERROF4);
* 'MESSAGE' ;
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'OPERATEUR KONV';
'MESSAGE' ('CHAINE' METO);
'ERREUR' 5;
'FINSI' ;
****************************************************
****************************************************
******** Fin boucle sur les angles *********
****************************************************
****************************************************
'FIN' BLOC;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer