Download konv_resi_gfmp_consist.dgibi
* fichier : konv_resi_gfmp_constant_state_11.dgibi
***********************************************************
**** APPROCHE VF "Cell-Centred Formulation" pour la ****
**** solution des ****
**** Equations d'Euler pour un gaz parfait ****
**** Approche GFMP ****
**** OPERATEURS 'PRIM', PRET, KONV ****
**** ****
**** Consistency ****
**** ****
**** Methodes: GODUNOV ****
**** ****
**** A. BECCANTINI DM2S/SFME/LTMF DECEMBRE 2010 ****
***********************************************************
'OPTION' 'DIME' 2 ;
'OPTION' 'ELEM' QUA4 ;
'OPTION' 'ECHO' 1 ;
'OPTION' 'TRAC' 'X' ;
*
*** GRAPH
*
GRAPH = FAUX ;
* GRAPH = VRAI ;
***************************
***** DOMAINE SPATIAL ****
***************************
A1 = 0.0D0 0.01D0;
A2 = 1.011D0 0.21D0;
A3 = 2.12D0 0.021D0;
A4 = 2.11D0 1.09D0;
A5 = 1.023D0 1.132D0;
A6 = 0.098D0 1.199D0;
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 ****************
*******************************************
*
* Case 1 -> 1
* Constant state
*
phg1 = 0.77 ;
phd1 = 1.56 ;
yg1 = 0.11;
yd1 = 0.13;
yg2 = 0.11;
yd2 = 0.13;
ag1 = 0.45;
ad1 = 0.27;
ag2 = 0.45;
ad2 = 0.27;
gam = 1.2D0 ;
pinf = 1.011d5 ;
rhog1 = 1.17e0 ;
rhod1 = 1.17e0 ;
pg1 = 1.023e5 ;
pd1 = 1.023e5 ;
ung1 = 122.0 ;
und1 = 122.0 ;
utg1 = 101.0 ;
utd1 = 101.0 ;
* utg1 = 0.0 ;
* utd1 = 0.0 ;
eing1 = (1. '/' (gam '-' 1.0D0)) * (pg1 '+' (gam * pinf)) ;
eing1 = eing1 '/' rhog1 ;
eind1 = (1. '/' (gam '-' 1.0D0)) * (pd1 '+' (gam * pinf)) ;
eind1 = eind1 '/' rhod1 ;
etg1 = 0.0 ;
utd1 = 0.0 ;
PN1 = 'PROG' pg1 pd1 ;
RN1 = 'PROG' rhog1 rhod1 ;
UN1 = 'PROG' ung1 und1 ;
UT1 = 'PROG' utg1 utd1 ;
EIN1 = 'PROG' eing1 eind1 ;
PHI1 = 'PROG' phg1 phd1 ;
Y1 = 'PROG' yg1 yd1 ;
Y2 = 'PROG' yg2 yd2 ;
ALPHA1 = 'PROG' ag1 ad1 ;
ALPHA2 = 'PROG' ag2 ad2 ;
*
* Computation of the conservative variables
*
GNX1 = RN1 * UN1 ;
GNY1 = RN1 * UT1 ;
ECIN1 = 0.5D0 '*' ('PROG' ((ung1 * ung1) '+' (utg1 * utg1))
((und1 * und1) '+' (utd1 * utd1)) );
RETN1 = RN1 '*' (EIN1 '+' ECIN1) ;
*
* Computation of the resi contribution in (n,t)
*
retg1 = MAXI RETN1 ;
rhog1 = 'MAXIMUM' RN1 ;
f0gd = ung1 * rhog1 * phg1 ;
f1gd = ung1 '*' rhog1 ;
f2gd = (f1gd '*' ung1) '+' pg1 ;
f3gd = f1gd '*' utg1 ;
f4gd = ung1 '*' (retg1 '+' pg1) ;
f5gd = ung1 * rhog1 * yg1 ;
f6g = 0.0 ;
f6d = ung1 * (ag1 '-' ad1) ;
****************************************************
****************************************************
******** 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;
P1FAC = P10 'TOURNER' ANGLE ORIG;
DOMTOT = DOM1 ET DOM2;
'ELIMINATION' DOMTOT 1D-6;
$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' ;
GRALP1 = 'NOMC' ('MOTS' 'UX' 'UY') ('DOMA' $DOMTOT 'XXNORMAF')
('MOTS' 'P1DX' 'P1DY') ;
'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) 1.E-6 ;
*
VINF = 'MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' 100. ;
*
**** Redefinition de P1FAC dans $DOMTOT 'FACE'
*
P1FAC = ('DOMA' $DOMTOT 'FACE') 'POIN' 'PROC' P1FAC ;
XSURF = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXSURFAC') P1FAC 'SCAL' ;
NX = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAF') P1FAC 'UX' ;
NY = 'EXTRAIRE' ('DOMA' $DOMTOT 'XXNORMAF') P1FAC 'UY' ;
TX = -1 * NY ;
TY = NX ;
GEOPC = ('DOMA' $DOMTOT 'FACEL') 'ELEM' 'APPUYE' 'LARGEMENT'
P1FAC ;
AA = 'CHANGER' GEOPC 'POI1' ;
PC1 = 'POIN' 1 AA ;
P2 = 'POIN' 1 ('DOMA' $DOMTOT 'CENTRE') ;
PC2 = 'POIN' 3 AA ;
* Si P3 n'existe pas, probleme en FACEL.
'SI' ('NEG' PC1 P2) ;
ORIENT = -1 ;
PCD = PC1 ;
PCG = PC2 ;
XVOLD = 'MAXIMUM' (TDOM1 . 'XXVOLUM') ;
XVOLG = 'MAXIMUM' (TDOM2 . 'XXVOLUM') ;
'SINON' ;
ORIENT = 1 ;
PCD = PC2 ;
PCG = PC1 ;
XVOLG = 'MAXIMUM' (TDOM1 . 'XXVOLUM') ;
XVOLD = 'MAXIMUM' (TDOM2 . 'XXVOLUM') ;
'FINSI' ;
uxg1 = ((ung1 '*' NX) '+' (utg1 '*' TX)) '*' ORIENT ;
uyg1 = ((ung1 '*' NY) '+' (utg1 '*' TY)) '*' ORIENT ;
uxd1 = ((und1 '*' NX) '+' (utd1 '*' TX)) '*' ORIENT ;
uyd1 = ((und1 '*' NY) '+' (utd1 '*' TY)) '*' ORIENT ;
'SI' GRAPH;
'TRACER' (('DOMA' $DOMTOT 'MAILLAGE') 'ET'
('COULEUR' ('DOMA' $DOMTOT 'FACEL') 'ROUGE') 'ET' P1FAC)
'TITRE' 'Domaine et FACEL';
'FINSI' ;
MAILIM = 'DIFF' ('DOMA' $DOMTOT 'FACE') ('MANUEL' 'POI1' P1FAC) ;
***********************
**** Les CHPOINTs ****
***********************
CHRN1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' RN1) ;
CHVN1 = ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 2 'UX' uxg1
'UY' uyg1) '+'
('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 2 'UX' uxd1
'UY' uyd1) ;
CHGN1 = CHRN1 '*' CHVN1 ;
CHRET1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' RETN1) ;
CHPHI1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 1 'SCAL' PHI1) ;
CHY1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'ESP1' Y1
'ESP2' Y2) ;
CHAL1 = ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'ESP1' ALPHA1
'ESP2' ALPHA2) ;
CHRPHI = CHPHI1 * CHRN1 ;
CHRY1 = CHY1 * CHRN1 ;
TABG = 'TABLE' ;
TABG . 'ESPEULE' = 'MOTS' ;
TABG . 'ESPNEULE' = 'ESP1' ;
TABG . 'MGAM' = 'PROG' GAM ;
TABG . 'PGAM' = 'PROG' GAM ;
TABG . 'MPIN' = 'PROG' PINF ;
TABG . 'PPIN' = 'PROG' PINF ;
TABG1 = 'TABLE' ;
TABG1 . 'ESPEULE' = 'MOTS' 'ESP1' 'ESP2' ;
TABG1 . 'ESPNEULE' = 'ESP3' ;
TABG1 . 'MGAM' = 'PROG' GAM GAM GAM ;
TABG1 . 'PGAM' = 'PROG' GAM GAM GAM ;
TABG1 . 'MPIN' = 'PROG' PINF PINF PINF ;
TABG1 . 'PPIN' = 'PROG' PINF PINF PINF ;
PH1 = CHRPHI '/' CHRN1 ;
V1 P1 = 'PRIM' 'GFMP' TABG CHRPHI CHRN1 CHGN1 CHRET1 ;
V1B P1B Y1B = 'PRIM' 'GFMP' TABG1 CHRPHI CHRN1 CHGN1 CHRET1
CHRY1 CHAL1 ;
ERRO = 'MAXIMUM' (Y1B '-' CHY1) 'ABS' ;
ERRO = ERRO '+' (
'MAXIMUM' ((P1 '-' P1B) '/' P1) 'ABS'
);
ERRO = ERRO '+' (
'MAXIMUM' (V1 '-' V1B) 'ABS'
);
SI (ERRO > 1.0D-12) ;
'MESSAGE' 'Problem in PRIM' ;
'MESSAGE' ('CHAINE' 'Error = ', ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
GRADPH1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 2
'P1DX' 0.0 'P1DY' 0.0 ;
LIMPH1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 1
'P1' 0.0 ;
GRADR1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 2
'P1DX' 0.0 'P1DY' 0.0 ;
LIMR1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 1
'P1' 0.0 ;
GRADV1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 4
'P1DX' 0.0 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ;
LIMV1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 2
'P1' 0.0 'P2' 0.0 ;
GRADP1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 2
'P1DX' 0.0 'P1DY' 0.0 ;
LIMP1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 1
'P1' 0.0 ;
GRADY1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 4
'P1DX' 0.0 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ;
LIMY1 = 'MANU' 'CHPO' (TDOMTOT . 'CENTRE') 2
'P1' 0.0 'P2' 0.0 ;
*
**** L'operateur 'PRET'
*
CHFPH1 CHFRN1 CHFVN1 CHFPN1 =
'PRET' 'GFMP' (TABG . 'ESPEULE') $DOMTOT
PH1 (0.0 * GRADPH1) LIMPH1
CHRN1 (0.0 * GRADR1) LIMR1
V1 (0.0 * GRADV1) LIMV1
P1 (0.0 * GRADP1) LIMP1 ;
CHFPH1B CHFRN1B CHFVN1B CHFPN1B
CHFYNB CHFALB =
'PRET' 'GFMP' (TABG1 . 'ESPEULE') $DOMTOT
PH1 (0.0 * GRADPH1) LIMPH1
CHRN1 (0.0 * GRADR1) LIMR1
V1 (0.0 * GRADV1) LIMV1
P1 (0.0 * GRADP1) LIMP1
CHY1 (0.0 * GRADY1) LIMY1
CHAL1 (0.0 * GRADY1) LIMY1;
*
**** L'operateur 'KONV'
*
'REPETER' BLMETO 1 ;
'SI' ('EGA' &BLMETO 1) ;
METO = 'GODUNOV' ;
'FINSI' ;
'MESSAGE' ('CHAINE' 'METO = ' METO) ;
LISTINC1 = ('MOTS' 'PHI1' 'RN1' 'RNX1' 'RUY1' 'RET1') ;
LISTINC2 = ('MOTS' 'PHI1' 'RN1' 'RNX1' 'RUY1' 'RET1'
'RY1' 'RY2' 'ALP1' 'ALP2') ;
CHPRES DT = 'KONV' 'VF' 'GFMP' 'RESI' METO
TABG $DOMTOT
LISTINC1 CHPHI1 CHFPH1 CHFRN1 CHFVN1 CHFPN1
MAILIM VRAI 0.0 ;
CHPRES1 DT1 = 'KONV' 'VF' 'GFMP' 'RESI' METO
TABG1 $DOMTOT
LISTINC2 CHPHI1 CHFPH1 CHFRN1 CHFVN1 CHFPN1
CHFYNB CHFALB CHAL1
MAILIM VRAI 0.0 ;
ERRO = 'EXCO' CHPRES1 LISTINC1 ;
ERRO = 'MAXIMUM' (ERRO '-' CHPRES) 'ABS' ;
ERRO = ERRO '/' ('MAXIMUM' CHPRES 'ABS') ;
SI (ERRO > 1.0D-12) ;
'MESSAGE' 'Problem in KONV 1' ;
'MESSAGE' ('CHAINE' 'Error = ', ERRO) ;
'ERREUR' 5 ;
'FINSI' ;
* 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' ETHER1 1)) ;
* 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' ECIN1 1)) ;
* 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' EFORM1 1)) ;
* 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' (pg1)) ;
* 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' ('EXTRAIRE' RN1 1)) ;
* 'LISTE' ('CHAINE' 'FORMAT' '(E16.10)' (ung1)) ;
* 'OPTION' DONN 5 ;
RESX0 = 'EXCO' ('EXTRAIRE' LISTINC1 1) CHPRES ;
RESX1 = 'EXCO' ('EXTRAIRE' LISTINC1 2) CHPRES ;
RESX2X = 'EXCO' ('EXTRAIRE' LISTINC1 3) CHPRES ;
RESX2Y = 'EXCO' ('EXTRAIRE' LISTINC1 4) CHPRES ;
RESX3 = 'EXCO' ('EXTRAIRE' LISTINC1 5) CHPRES ;
RESX4 = 'EXCO' ('EXTRAIRE' LISTINC2 6) CHPRES1 ;
RESX5 = 'EXCO' ('EXTRAIRE' LISTINC2 8) CHPRES1 ;
RESX2N = (RESX2X '*' NX) '+' (RESX2Y * NY);
RESX2T = (RESX2X '*' TX) '+' (RESX2Y * TY);
f0 = 'EXTRAIRE' RESX0 'SCAL' PCD ;
f0bis = 'EXTRAIRE' RESX0 'SCAL' PCG ;
ERRO = 1D-8 '*' f0gd 'ABS' ;
LOGI1 = ('ABS' ((f0 * XVOLD) + (f0bis * XVOLG))) < ERRO ;
LOGI2 = ('ABS' ((f0gd * XSURF) + (f0bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
f1 = 'EXTRAIRE' RESX1 'SCAL' PCD ;
f1bis = 'EXTRAIRE' RESX1 'SCAL' PCG ;
ERRO = 1D-8 '*' f1gd 'ABS' ;
LOGI2 = ('ABS' ((f1 * XVOLD) + (f1bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI2 = ('ABS' ((f1gd * XSURF) + (f1bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI1 = LOGI1 'ET' LOGI2;
f2 = 'EXTRAIRE' RESX2N 'SCAL' PCD ;
f2bis = 'EXTRAIRE' RESX2N 'SCAL' PCG ;
ERRO = 1D-8 '*' f2gd 'ABS' ;
LOGI2 = ('ABS' ((f2 * XVOLD) + (f2bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI2 = ('ABS' ((f2gd * XSURF) + (f2bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
f3 = 'EXTRAIRE' RESX2T 'SCAL' PCD ;
f3bis = 'EXTRAIRE' RESX2T 'SCAL' PCG ;
ERRO = 1D-8 '*' f3gd 'ABS' ;
LOGI2 = ('ABS' ((f3 * XVOLD) + (f3bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI2 = ('ABS' ((f3gd * XSURF) + (f3bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
f4 = 'EXTRAIRE' RESX3 'SCAL' PCD ;
f4bis = 'EXTRAIRE' RESX3 'SCAL' PCG ;
ERRO = 5D-7 '*' f4gd 'ABS' ;
LOGI2 = ('ABS' ((f4 * XVOLD) + (f4bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI2 = ('ABS' ((f4gd * XSURF) + (f4bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
* RHOY
f5 = 'EXTRAIRE' RESX4 'SCAL' PCD ;
f5bis = 'EXTRAIRE' RESX4 'SCAL' PCG ;
ERRO = 1D-8 '*' f5gd 'ABS' ;
LOGI2 = ('ABS' ((f5 * XVOLD) + (f5bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI2 = ('ABS' ((f5gd * XSURF) + (f5bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI1 = LOGI1 'ET' LOGI2;
* ALPHA ('NON' conservative form...)
f6 = 'EXTRAIRE' RESX5 'SCAL' PCD ;
f6bis = 'EXTRAIRE' RESX5 'SCAL' PCG ;
ERRO = 1D-8 '*' ('MAXIMUM' ('PROG' f6 f6bis) 'ABS') ;
LOGI2 = ('ABS' ((f6d * XSURF) - (f6 * XVOLD))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI2 = ('ABS' ((f6g * XSURF) + (f6bis * XVOLG))) < ERRO ;
LOGI1 = LOGI1 'ET' LOGI2;
LOGI1 = LOGI1 'ET' LOGI2;
'SI' ('NON' LOGI1);
'MESSAGE' ;
'MESSAGE' 'OPERATEUR KONV';
'MESSAGE' ('CHAINE' METO);
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' BLMETO ;
****************************************************
****************************************************
******** Fin boucle sur les angles *********
****************************************************
****************************************************
'FIN' BLOC;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer