Télécharger 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' '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'; '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 ; * * Computation of the conservative variables * GNX1 = RN1 * UN1 ; GNY1 = RN1 * UT1 ; ((und1 * und1) '+' (utd1 * utd1)) ); RETN1 = RN1 '*' (EIN1 '+' ECIN1) ; * * Computation of the resi contribution in (n,t) * 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'; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMTOT = TDOMTOT . 'QUAF' ; * * **** Redefinition de P1FAC dans $DOMTOT 'FACE' * TX = -1 * NY ; TY = NX ; P1FAC ; AA = 'CHANGER' GEOPC 'POI1' ; * 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; 'TITRE' 'Domaine et FACEL'; 'FINSI' ; *********************** **** Les CHPOINTs **** *********************** 'UY' uyg1) '+' 'UY' uyd1) ; CHGN1 = CHRN1 '*' CHVN1 ; 'ESP2' Y2) ; 'ESP2' ALPHA2) ; CHRPHI = CHPHI1 * CHRN1 ; CHRY1 = CHY1 * CHRN1 ; TABG = 'TABLE' ; TABG . 'ESPNEULE' = 'ESP1' ; TABG1 = 'TABLE' ; TABG1 . 'ESPNEULE' = 'ESP3' ; PH1 = CHRPHI '/' CHRN1 ; 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' ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ; 'P1' 0.0 'P2' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 ; 'P1' 0.0 ; 'P1DX' 0.0 'P1DY' 0.0 'P2DX' 0.0 'P2DY' 0.0 ; 'P1' 0.0 'P2' 0.0 ; * **** L'operateur 'PRET' * CHFPH1 CHFRN1 CHFVN1 CHFPN1 = 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 = 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) ; 'RY1' 'RY2' 'ALP1' 'ALP2') ; TABG $DOMTOT LISTINC1 CHPHI1 CHFPH1 CHFRN1 CHFVN1 CHFPN1 MAILIM VRAI 0.0 ; TABG1 $DOMTOT LISTINC2 CHPHI1 CHFPH1 CHFRN1 CHFVN1 CHFPN1 CHFYNB CHFALB CHAL1 MAILIM VRAI 0.0 ; ERRO = 'MAXIMUM' (ERRO '-' CHPRES) 'ABS' ; ERRO = ERRO '/' ('MAXIMUM' CHPRES 'ABS') ; SI (ERRO > 1.0D-12) ; '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 ; 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 ; 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 - Tous droits réservés.
Mentions légales