* fichier : konv_cons.dgibi
***********************************************************
***********************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la ****
**** solution des ****
**** Equations d Euler pour un gaz parfait ****
**** OPERATEURS PRIM, PRET, KONV ****
**** ****
**** Cas gaz monoespece, "calorically perfect" ****
**** Consistence ****
**** ****
**** Methodes: VANLEER ****
**** VLH ****
**** HUSVL ****
**** HUSVLH ****
**** GODUNOV ****
**** AUSMPLUS ****
**** ROE ****
**** SS ****
**** RUSANOV ****
**** AUSMPLM ****
**** RUSANOLM ****
**** CENTERED ****
**** HLLC ****
**** HLLCLM ****
**** ****
**** A. BECCANTINI DRN/DMT/SEMT/LTMF DECEMBRE 2003 ****
***********************************************************
***********************************************************
'OPTION' 'DIME' 2 'ELEM' QUA4 'ECHO' 1 'TRAC' 'X' ;
*
*** GRAPH
*
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 la consistence
*
P10 = 1.0 0.5;
*
*** Etats gauche et droite eguals
*
rogd = 1.11;
ungd = 1.3;
utgd = 18.2;
pgd = 1234.7;
gamgd = 1.4D0;
*
*** retgd
*
ecingd = 0.5D0 '*' rogd '*' ((ungd '*' ungd) '+' (utgd '*' utgd));
retgd = (pgd '/' (gamgd '-' 1.0)) '+' ecingd;
*
*** flux en (n,t)
*
f1gd = ungd '*' rogd;
f2gd = (rogd '*' (ungd '*' ungd)) '+' pgd;
f3gd = rogd '*' ungd '*' utgd ;
f4gd = ungd * (retgd '+' pgd);
****************************************************
****************************************************
******** 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' ('CHAINE' '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 ****
***********************
uxgd = (ungd '*' ('COS' ANGLE)) '-' (utgd '*' ('SIN' ANGLE));
uygd = (ungd '*' ('SIN' ANGLE)) '+' (utgd '*' ('COS' ANGLE));
GAMMAN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gamgd;
RN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rogd;
RN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rogd;;
RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RN1 'ET' RN2);
GN1 = 'KCHT' $DOM1 'VECT' 'CENTRE' ((rogd '*' uxgd)
(rogd '*' uygd));
GN2 = 'KCHT' $DOM2 'VECT' 'CENTRE' ((rogd '*' uxgd)
(rogd '*' uygd));
GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN1 'ET' GN2);
PN1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' pgd;
PN2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' pgd;
PN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (PN1 'ET' PN2);
*
*** La vitesse
*
VITX = ('EXCO' 'UX ' GN) '/' RN;
VITY = ('EXCO' 'UY ' GN) '/' RN;
*
*** L'energie totale (RETN)
*
GM1 = GAMMAN '-' ('KCHT' $DOMTOT 'SCAL' 'CENTRE' 1.0D0) ;
RETHER = PN '/' GM1 ;
RECIN = (0.5D0 '*' ('KOPS' GN 'PSCA' GN )) '/' RN;
RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RECIN '+' RETHER);
***************************
**** L'operateur PRIM****
***************************
VITESSE PRES = 'PRIM' 'PERFMONO' RN GN RETN GAMMAN ;
***************************
**** L'operateur PRET****
***************************
ORDESP = 1;
ORDTEM = 1;
ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORDESP ORDTEM
$DOMTOT RN VITESSE PRES GAMMAN ;
*********************************************************
*** 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;
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', premier ordre en espace et en temps.
*
ERRO = 1D-13 '*' rogd;
LOGI1 = (('ABS' (rog '-' rod)) < ERRO) 'ET'
(('ABS' (rog '-' rogd)) < ERRO);
ERRO = 1D-13 '*' pgd;
LOGI2 = (('ABS' (pg '-' pd )) < ERRO) 'ET'
(('ABS' (pg '-' pgd)) < ERRO) ;
LOGI1 = LOGI1 'ET' LOGI2;
ERRO = 1D-13 '*' ungd;
LOGI2 = (('ABS' (ung '-' und)) < ERRO) 'ET'
(('ABS' (ung '-' (ORIENT '*'ungd))) < ERRO) ;
LOGI1 = LOGI1 'ET' LOGI2;
ERRO = 1D-13 '*' utgd;
LOGI2 = (( 'ABS' (utg '-' utd)) < ERRO) 'ET'
(('ABS' (utg '-' (ORIENT '*'utgd))) < ERRO) ;
LOGI1 = LOGI1 'ET' LOGI2;
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'Probleme en PRET' ;
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'dRO ' ('ABS' (rog '-' rod))
' ' ('ABS' (rog '-' rogd)));
'MESSAGE' ('CHAINE' 'dp ' ('ABS' (pg '-' pd))
' ' ('ABS' (pg '-' pgd)));
'MESSAGE' ('CHAINE' 'dun ' ('ABS' (ung '-' und))
' ' ('ABS' (ung '-' (ORIENT '*'ungd))));
'MESSAGE' ('CHAINE' 'dRO ' ('ABS' (utg '-' utd))
' ' ('ABS' (utg '-' (ORIENT '*'utgd))));
'ERREUR' 5;
'FINSI' ;
ORDESP = 1 ;
ORDTEM = 1 ;
'REPETER' BLMETO 16 ;
'SI' ('EGA' &BLMETO 1) ;
METO = 'VANLEER' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 2) ;
METO = 'VLH' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 3) ;
METO = 'HUSVL' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 4) ;
METO = 'HUSVLH' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 5) ;
METO = 'GODUNOV' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 6) ;
METO = 'AUSMPLUS' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 7) ;
METO = 'ROE' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 8) ;
METO = 'SS' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 9) ;
METO = 'RUSANOV' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 10) ;
METO = 'AUSMPLM' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 11) ;
METO = 'RUSANOLM' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 12) ;
METO = 'CENTERED' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 13) ;
METO = 'ROELM' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 14) ;
METO = 'HLLC' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 15) ;
METO = 'HLLCLM' ;
'FINSI' ;
'SI' ('EGA' &BLMETO 16) ;
METO = 'AUSMPUP' ;
'FINSI' ;
'MESSAGE' ('CHAINE' 'METO = ' METO) ;
LISTINCO = ('MOTS' 'RN' 'RUNX' 'RUNY' 'RETN') ;
'SI' ((&BLMETO 'EGA' 10) 'OU' (&BLMETO 'EGA' 11)
'OU' (&BLMETO 'EGA' 13) 'OU' (&BLMETO 'EGA' 15)
'OU' (&BLMETO 'EGA' 16)) ;
VCUTOFF = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL'
0.0001 ;
CHPFLU DT = 'KONV' 'VF' 'PERFMONO' 'FLUX' METO
$DOMTOT ROF VITF PF GAMF LISTINCO VCUTOFF VCUTOFF ;
'SINON' ;
CHPFLU DT = 'KONV' 'VF' 'PERFMONO' 'FLUX' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
'FINSI' ;
*
**** Les flux aux interfaces sont dans de
* CHPOINT FACE
*
FLUX1 = 'EXCO' ('EXTRAIRE' LISTINCO 1) CHPFLU ;
FLUX2X = 'EXCO' ('EXTRAIRE' LISTINCO 2) CHPFLU ;
FLUX2Y = 'EXCO' ('EXTRAIRE' LISTINCO 3) CHPFLU ;
FLUX3 = 'EXCO' ('EXTRAIRE' LISTINCO 4) CHPFLU ;
FLUX2N = (FLUX2X '*' ('COS' ANGLE)) '+' (FLUX2Y * ('SIN' ANGLE));
FLUX2T = (FLUX2Y '*' ('COS' ANGLE)) '-' (FLUX2X * ('SIN' ANGLE));
f1 = 'EXTRAIRE' FLUX1 'SCAL' P1;
ERRO = 1D-13 '*' f1gd;
LOGI1 = ('ABS' ((f1 '*' ORIENT) '-' f1gd)) < ERRO ;
f2 = 'EXTRAIRE' FLUX2N 'SCAL' P1;
ERRO = 1D-13 '*' f2gd;
LOGI2 = ('ABS' ((f2 '*' ORIENT) '-' f2gd)) < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
f3 = 'EXTRAIRE' FLUX2T 'SCAL' P1;
ERRO = 1D-13 '*' f3gd;
LOGI2 = ('ABS' ((f3 '*' ORIENT) '-' f3gd)) < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
f4 = 'EXTRAIRE' FLUX3 'SCAL' P1;
ERRO = 1D-13 '*' f4gd;
LOGI2 = ('ABS' ((f4 '*' ORIENT) '-' f4gd)) < ERRO;
LOGI1 = LOGI1 'ET' LOGI2;
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'OPERATEUR KONV';
'MESSAGE' ('CHAINE' METO);
'MESSAGE' ;
'MESSAGE' ('CHAINE' 'df1 ' ('ABS' (f1gd '-' (ORIENT '*'f1)))
' f1 ' f1);
'MESSAGE' ('CHAINE' 'df2 ' ('ABS' (f2gd '-' (ORIENT '*'f2)))
' f2 ' f2);
'MESSAGE' ('CHAINE' 'df3 ' ('ABS' (f3gd '-' (ORIENT '*'f3)))
' f3 ' f3);
'MESSAGE' ('CHAINE' 'df4 ' ('ABS' (f4gd '-' (ORIENT '*'f4)))
' f4 ' f4);
'ERREUR' 5;
'FINSI' ;
*
****** Le residu
*
*
****** CHPDIR = chpoint sur les face de chaque elts
*
*
* NB Dans ce cas particulier, DOMTOT n'a que de QUA4,
* donc ($DOMTOT . 'XXNORMAE') contient une seule
* zona elementaire!
RESI MATCEL = 'KOPS' 'MATRIK' ;
MAICHA = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAE') 'MAILLAGE' ;
NBELEM = 'NBEL' MAICHA ;
'REPETER' BLELE NBELEM ;
CHPDIR MATCEL = 'KOPS' 'MATRIK' ;
ELEM1 = &BLELE 'ELEM' MAICHA ;
ELEMC1 = 'CHANGER' ELEM1 'POI1' ;
NBPOIN = 'NBNO' ELEMC1 ;
PCEN = ('DOMA' $DOMTOT 'CENTRE') 'POIN' &BLELE ;
XVOL = ('DOMA' $DOMTOT 'XXVOLUM') 'EXTRAIRE' PCEN 'SCAL' ;
MAICEN = 'MANU' 'POI1'PCEN ;
'REPETER' BLEPO NBPOIN ;
P1 = &BLEPO 'POIN' ELEMC1 ;
ELEP1 = 'MANUEL' 'POI1' P1 ;
COEFNO = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAE')
'SCAL' 1 &BLELE &BLEPO ;
CHPDIR = CHPDIR 'ET' ('MANUEL' 'CHPO' ELEP1 1 'SCAL'
COEFNO 'NATU' 'DISCRET') ;
'FIN' BLEPO ;
MAIDIR = 'EXTRAIRE' CHPDIR 'MAILLAGE' ;
*
********* Le flux qui sort
*
FLUSOR = CHPDIR '*' ('REDU' CHPFLU MAIDIR)
('MOTS' 'SCAL' 'SCAL' 'SCAL' 'SCAL')
LISTINCO LISTINCO ;
CHPUN = 'MANUEL' 'CHPO' MAIDIR 1 'SCAL' (-1./XVOL) ;
NBCOMP = 'DIME' LISTINCO ;
'REPETER' BLCOM NBCOMP ;
NOMCOM = 'EXTRAIRE' LISTINCO &BLCOM ;
CELSCA = 'XTY' FLUSOR CHPUN ('MOTS' NOMCOM)
('MOTS' 'SCAL') ;
RESI = RESI 'ET' ('MANUEL' 'CHPO' MAICEN 1
NOMCOM CELSCA 'NATU' 'DISCRET') ;
'FIN' BLCOM ;
'FIN' BLELE ;
'SI' ((&BLMETO 'EGA' 10) 'OU' (&BLMETO 'EGA' 11)
'OU' (&BLMETO 'EGA' 13) 'OU' (&BLMETO 'EGA' 15)
'OU' (&BLMETO 'EGA' 16)) ;
VCUTOFF = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL'
0.00001 ;
CHPRES DT1 = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO VCUTOFF VCUTOFF ;
'SINON' ;
CHPRES DT1 = 'KONV' 'VF' 'PERFMONO' 'RESI' METO
$DOMTOT ROF VITF PF GAMF LISTINCO ;
'FINSI' ;
*
ERRO = 'MAXIMUM' (RESI '-' CHPRES) 'ABS' ;
'SI' (ERRO > 1.0D-11) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dans KONV' ;
'ERREUR' 5 ;
'FINSI' ;
ERRO = 'ABS' (DT '-' DT1) ;
'SI' (ERRO > 1.0D-11) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dans KONV' ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' BLMETO ;
****************************************************
****************************************************
******** Fin boucle sur les angles *********
****************************************************
****************************************************
'FIN' BLOC;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer