* fichier : pret_gfmp.dgibi ************************************************************************ ************************************************************************ ********************************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la solution des **** **** Equations d'Euler. **** **** "Ghost Fluid Method for the Poor" **** **** **** **** OPERATEUR PRET **** **** Operateur qui 'recontruit les variables primitives aux faces **** **** "Stiffened gas" **** **** Interieur et murs **** **** **** **** A. BECCANTINI DEN/DM2S/SFME/LTMF MAI 2011 **** ********************************************************************** 'TRAC' 'X' ; * *** GRAPH * GRAPH = VRAI ; GRAPH = FAUX ; *************************** ***** DOMAINE SPATIAL **** *************************** * **** Deux carre * A1 = 0.0D0 0.0D0; A2 = 2.0D0 0.0D0; A3 = 2.2D0 0.0D0; A4 = 2.0D0 0.8D0; A5 = 1.0D0 0.90D0; A6 = 0.0D0 1.1D0; 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'; * *** Etats a gauche et a droite * * al1g = 0.9 ; limgal1 = 0.11 ; rn1g = 1.10 ; limgrn1 = 0.21 ; ux1g = 130. ; limgux1 = 0.31 ; uy1g = 111. ; limguy1 = 0.41 ; pn1g = 1.13E5 ; limgpn1 = 0.51 ; xgragal1 = 11.0D-2 ; ygragal1 = 12.0D-2 ; xgragrn1 = 13.0D-1 ; ygragrn1 = 11.0D-1 ; xgragux1 = 15.0 ; ygragux1 = 12.0 ; xgraguy1 = 17.0 ; ygraguy1 = 19.0 ; xgragpn1 = 21.0D2 ; ygragpn1 = 21.2D2 ; al1d = 0.8 ; limdal1 = 0.17 ; rn1d = 1.31 ; limdrn1 = 0.19 ; ux1d = 230. ; limdux1 = 0.22 ; uy1d = 431. ; limduy1 = 0.32 ; pn1d = 1.31E5 ; limdpn1 = 0.42 ; xgradal1 = 11.0D-2 ; ygradal1 = 11.2D-2 ; xgradrn1 = 21.0D-1 ; ygradrn1 = 21.2D-1 ; xgradux1 = 31.0 ; ygradux1 = 31.2 ; xgraduy1 = 41.0 ; ygraduy1 = 41.3 ; xgradpn1 = 51.0D1 ; ygradpn1 = 51.2D1 ; DOM1 = DOM10 ; DOM2 = DOM20 ; DOMTOT = DOM1 ET DOM2; 'ELIMINATION' (DOM1 ET DOM2) 1D-6; $DOMTOT = 'MODELISER' DOMTOT 'EULER'; $DOM1 = 'MODELISER' DOM1 'EULER'; $DOM2 = 'MODELISER' DOM2 'EULER'; $L25 = 'MODELISER' L25 'EULER'; $L12 = 'MODELISER' L12 'EULER'; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMTOT = TDOMTOT . 'QUAF' ; ML25 = TL25 . 'QUAF' ; ML12 = TL12 . 'QUAF' ; 'ELIMINATION' (MDOMTOT ET MDOM1 ET ML25) 0.0001 ; 'ELIMINATION' (MDOMTOT ET MDOM2 ET ML12) 0.0001 ; 'SI' GRAPH; ('MANUEL' 'POI1' P1 'COULEUR' 'ROUG') 'ET' ('MANUEL' 'POI1' P2 'COULEUR' 'ROUG')) 'TITRE' 'Domaine et points'; 'FINSI' ; **** La densité, la pression, la vitesse 'P1DX' xgragal1 'P1DY' ygragal1) '+' 'P1DX' xgradal1 'P1DY' ygradal1) ; 'P1DX' xgragrn1 'P1DY' ygragrn1) '+' 'P1DX' xgradrn1 'P1DY' ygradrn1) ; 'P1DX' xgragpn1 'P1DY' ygragpn1) '+' 'P1DX' xgradpn1 'P1DY' ygradpn1) ; 'P1DX' xgragux1 'P1DY' ygragux1) '+' 'P1DX' xgradux1 'P1DY' ygradux1) ; 'P2DX' xgraguy1 'P2DY' ygraguy1) '+' 'P2DX' xgraduy1 'P2DY' ygraduy1) ; VN1 = UXN1 '+' UYN1 ; GRADVN1 = GRADUX1 '+' GRADUY1 ; LIMVN1 = LIMUX1 '+' LIMUY1 ; ************************************************************************* *** We check if the reconstruction degenerates to first order as the **** *** gradients or the limiters are zero **** ************************************************************************* AL1G0 RN1G0 VN1G0 PN1G0 = 'PRET' 'GFMP' NESP $DOMTOT AL1 (0.0 * GRADAL1) LIMAL1 RN1 (0.0 * GRADRN1) LIMRN1 VN1 (0.0 * GRADVN1) LIMVN1 PN1 (0.0 * GRADPN1) LIMPN1 ; AL1L0 RN1L0 VN1L0 PN1L0 = 'PRET' 'GFMP' NESP $DOMTOT AL1 GRADAL1 (0.0 * LIMAL1) RN1 GRADRN1 (0.0 * LIMRN1) VN1 GRADVN1 (0.0 * LIMVN1) PN1 GRADPN1 (0.0 * LIMPN1) ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-16) ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; * *** Control des etats sur la surface qui contient P1 * al1gp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 1 ; al1dp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 3 ; rn1gp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 1 ; rn1dp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 3 ; pn1gp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 1 ; pn1dp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 3 ; un1gp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 1 ; un1dp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 3 ; ut1gp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 1 ; ut1dp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 3 ; nxp1 = 'EXTRAIRE' TN1GEP1F 'NX ' 1 1 1 ; nyp1 = 'EXTRAIRE' TN1GEP1F 'NY ' 1 1 1 ; txp1 = 'EXTRAIRE' TN1GEP1F 'TX ' 1 1 1 ; typ1 = 'EXTRAIRE' TN1GEP1F 'TY ' 1 1 1 ; ux1gp1 = (un1gp1 * nxp1) '+' (ut1gp1 * txp1) ; uy1gp1 = (un1gp1 * nyp1) '+' (ut1gp1 * typ1) ; ux1dp1 = (un1dp1 * nxp1) '+' (ut1dp1 * txp1) ; uy1dp1 = (un1dp1 * nyp1) '+' (ut1dp1 * typ1) ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-12) ; 'MESSAGE' 'Les etats' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; * **** Wall * al1gp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 1 ; al1dp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 3 ; rn1gp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 1 ; rn1dp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 3 ; pn1gp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 1 ; pn1dp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 3 ; un1gp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 1 ; un1dp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 3 ; ut1gp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 1 ; ut1dp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 3 ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-16) ; 'MESSAGE' 'Les etats (mur)' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; ***************************************************** *** We check the second order reconstruction. ******* ***************************************************** AL1F RN1F VN1F PN1F = 'PRET' 'GFMP' NESP $DOMTOT AL1 GRADAL1 LIMAL1 RN1 GRADRN1 LIMRN1 VN1 GRADVN1 LIMVN1 PN1 GRADPN1 LIMPN1 ; *** Internal point al1gp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 1 ; al1dp1 = 'EXTRAIRE' AL1GEOP1 'SCAL' 1 1 3 ; rn1gp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 1 ; rn1dp1 = 'EXTRAIRE' RN1GEOP1 'SCAL' 1 1 3 ; pn1gp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 1 ; pn1dp1 = 'EXTRAIRE' PN1GEOP1 'SCAL' 1 1 3 ; un1gp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 1 ; un1dp1 = 'EXTRAIRE' VN1GEOP1 'UN ' 1 1 3 ; ut1gp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 1 ; ut1dp1 = 'EXTRAIRE' VN1GEOP1 'UT ' 1 1 3 ; nxp1 = 'EXTRAIRE' TN1GEP1F 'NX ' 1 1 1 ; nyp1 = 'EXTRAIRE' TN1GEP1F 'NY ' 1 1 1 ; txp1 = 'EXTRAIRE' TN1GEP1F 'TX ' 1 1 1 ; typ1 = 'EXTRAIRE' TN1GEP1F 'TY ' 1 1 1 ; ux1gp1 = (un1gp1 * nxp1) '+' (ut1gp1 * txp1) ; uy1gp1 = (un1gp1 * nyp1) '+' (ut1gp1 * typ1) ; ux1dp1 = (un1dp1 * nxp1) '+' (ut1dp1 * txp1) ; uy1dp1 = (un1dp1 * nyp1) '+' (ut1dp1 * typ1) ; * Selection of the left and the right state 'MESSAGE' 'Probleme en maicen' ; 'ERREUR' 5 ; 'FINSI' ; DXPC1 = ('COORDONNEE' 1 P1) '-' ('COORDONNEE' 1 PC1) ; DYPC1 = ('COORDONNEE' 2 P1) '-' ('COORDONNEE' 2 PC1) ; DXPC2 = ('COORDONNEE' 1 P1) '-' ('COORDONNEE' 1 PC2) ; DYPC2 = ('COORDONNEE' 2 P1) '-' ('COORDONNEE' 2 PC2) ; PSPC1 = (DXPC1 * nxp1) '+' (DYPC1 * nyp1) ; PSPC2 = (DXPC2 * nxp1) '+' (DYPC2 * nyp1) ; LOG1 = (PSPC1 > 0.0) 'ET' (PSPC2 < 0.0) ; 'SI' LOG1 ; PG = PC1 ; DXG = DXPC1 ; DYG = DYPC1 ; PD = PC2 ; DXD = DXPC2 ; DYD = DYPC2 ; 'FINSI' ; LOG2 = (PSPC1 < 0.0) 'ET' (PSPC2 > 0.0) ; 'SI' LOG2 ; PG = PC2 ; DXG = DXPC2 ; DYG = DYPC2 ; PD = PC1 ; DXD = DXPC1 ; DYD = DYPC1 ; 'FINSI' ; 'SI' (('NON' LOG1) 'ET' ('NON' LOG2)) ; 'MESSAGE' 'Probleme en facel' ; 'ERREUR' 5 ; 'FINSI' ; al1gn = ('EXTRAIRE' AL1 'SCAL' PG) '+' ( ( (('EXTRAIRE' GRADAL1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADAL1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMAL1 'P1' PG)) ; rn1gn = ('EXTRAIRE' RN1 'SCAL' PG) '+' ( ( (('EXTRAIRE' GRADRN1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADRN1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMRN1 'P1' PG)) ; pn1gn = ('EXTRAIRE' PN1 'SCAL' PG) '+' ( ( (('EXTRAIRE' GRADPN1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADPN1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMPN1 'P1' PG)) ; ux1gn = ('EXTRAIRE' VN1 'UX ' PG) '+' ( ( (('EXTRAIRE' GRADUX1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADUX1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMUX1 'P1' PG)) ; uy1gn = ('EXTRAIRE' VN1 'UY ' PG) '+' ( ( (('EXTRAIRE' GRADUY1 'P2DX' PG) * DXG) '+' (('EXTRAIRE' GRADUY1 'P2DY' PG) * DYG) ) * ('EXTRAIRE' LIMUY1 'P2' PG)) ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-13) ; 'MESSAGE' 'Deuxieme ordre' ; 'MESSAGE' 'Les etats gauches' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; al1dn = ('EXTRAIRE' AL1 'SCAL' PD) '+' ( ( (('EXTRAIRE' GRADAL1 'P1DX' PD) * DXD) '+' (('EXTRAIRE' GRADAL1 'P1DY' PD) * DYD) ) * ('EXTRAIRE' LIMAL1 'P1' PD)) ; rn1dn = ('EXTRAIRE' RN1 'SCAL' PD) '+' ( ( (('EXTRAIRE' GRADRN1 'P1DX' PD) * DXD) '+' (('EXTRAIRE' GRADRN1 'P1DY' PD) * DYD) ) * ('EXTRAIRE' LIMRN1 'P1' PD)) ; pn1dn = ('EXTRAIRE' PN1 'SCAL' PD) '+' ( ( (('EXTRAIRE' GRADPN1 'P1DX' PD) * DXD) '+' (('EXTRAIRE' GRADPN1 'P1DY' PD) * DYD) ) * ('EXTRAIRE' LIMPN1 'P1' PD)) ; ux1dn = ('EXTRAIRE' VN1 'UX ' PD) '+' ( ( (('EXTRAIRE' GRADUX1 'P1DX' PD) * DXD) '+' (('EXTRAIRE' GRADUX1 'P1DY' PD) * DYD) ) * ('EXTRAIRE' LIMUX1 'P1' PD)) ; uy1dn = ('EXTRAIRE' VN1 'UY ' PD) '+' ( ( (('EXTRAIRE' GRADUY1 'P2DX' PD) * DXD) '+' (('EXTRAIRE' GRADUY1 'P2DY' PD) * DYD) ) * ('EXTRAIRE' LIMUY1 'P2' PD)) ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-13) ; 'MESSAGE' 'Deuxieme ordre' ; 'MESSAGE' 'Les etats droites' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; *** Wall point 'SI' GRAPH ; 'TRACER' (DOMTOT 'ET' ('MANUEL' 'POI1' P2 'COULEUR' 'ROUG')) 'TITRE' 'Wall' ; 'FINSI' ; al1gp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 1 ; al1dp2 = 'EXTRAIRE' AL1GEOP2 'SCAL' 1 1 3 ; rn1gp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 1 ; rn1dp2 = 'EXTRAIRE' RN1GEOP2 'SCAL' 1 1 3 ; pn1gp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 1 ; pn1dp2 = 'EXTRAIRE' PN1GEOP2 'SCAL' 1 1 3 ; un1gp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 1 ; un1dp2 = 'EXTRAIRE' VN1GEOP2 'UN ' 1 1 3 ; ut1gp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 1 ; ut1dp2 = 'EXTRAIRE' VN1GEOP2 'UT ' 1 1 3 ; nxp2 = 'EXTRAIRE' TN1GEP2F 'NX ' 1 1 1 ; nyp2 = 'EXTRAIRE' TN1GEP2F 'NY ' 1 1 1 ; txp2 = 'EXTRAIRE' TN1GEP2F 'TX ' 1 1 1 ; typ2 = 'EXTRAIRE' TN1GEP2F 'TY ' 1 1 1 ; ux1gp2 = (un1gp2 * nxp2) '+' (ut1gp2 * txp2) ; uy1gp2 = (un1gp2 * nyp2) '+' (ut1gp2 * typ2) ; * Selection of the left and the right state 'MESSAGE' 'Probleme en maicen' ; 'ERREUR' 5 ; 'FINSI' ; DXPC1 = ('COORDONNEE' 1 P2) '-' ('COORDONNEE' 1 PC1) ; DYPC1 = ('COORDONNEE' 2 P2) '-' ('COORDONNEE' 2 PC1) ; PSPC1 = (DXPC1 * nxp2) '+' (DYPC1 * nyp2) ; LOG1 = (PSPC1 > 0.0) ; 'SI' LOG1 ; PG = PC1 ; DXG = DXPC1 ; DYG = DYPC1 ; 'FINSI' ; 'SI' (('NON' LOG1)) ; 'MESSAGE' 'Probleme en facel' ; 'ERREUR' 5 ; 'FINSI' ; al1gn = ('EXTRAIRE' AL1 'SCAL' PG) '+' ( ( (('EXTRAIRE' GRADAL1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADAL1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMAL1 'P1' PG)) ; rn1gn = ('EXTRAIRE' RN1 'SCAL' PG) '+' ( ( (('EXTRAIRE' GRADRN1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADRN1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMRN1 'P1' PG)) ; pn1gn = ('EXTRAIRE' PN1 'SCAL' PG) '+' ( ( (('EXTRAIRE' GRADPN1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADPN1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMPN1 'P1' PG)) ; ux1gn = ('EXTRAIRE' VN1 'UX ' PG) '+' ( ( (('EXTRAIRE' GRADUX1 'P1DX' PG) * DXG) '+' (('EXTRAIRE' GRADUX1 'P1DY' PG) * DYG) ) * ('EXTRAIRE' LIMUX1 'P1' PG)) ; uy1gn = ('EXTRAIRE' VN1 'UY ' PG) '+' ( ( (('EXTRAIRE' GRADUY1 'P2DX' PG) * DXG) '+' (('EXTRAIRE' GRADUY1 'P2DY' PG) * DYG) ) * ('EXTRAIRE' LIMUY1 'P2' PG)) ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-13) ; 'MESSAGE' 'Deuxieme ordre' ; 'MESSAGE' 'Les etats gauches au mur' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; ERRO = 'MAXIMUM' ERRO 'ABS' ; 'SI' (ERRO > 1.0D-13) ; 'MESSAGE' 'Deuxieme ordre' ; 'MESSAGE' 'Les etats gauches au mur' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; * One species involved * For the sake of simplicity, we put Y = alpha = phi * and we test with respect to the previous case. * PH1F RN1F VN1F PN1F = 'PRET' 'GFMP' NESP $DOMTOT AL1 GRADAL1 LIMAL1 RN1 GRADRN1 LIMRN1 VN1 GRADVN1 LIMVN1 PN1 GRADPN1 LIMPN1 ; PH0F RN1F VN1F PN1F = 'PRET' 'GFMP' NESP $DOMTOT AL1 GRADAL1 (LIMAL1 * 0.0D0) RN1 GRADRN1 LIMRN1 VN1 GRADVN1 LIMVN1 PN1 GRADPN1 LIMPN1 ; PHFM RN1F VN1F PN1F YNF ALF = 'PRET' 'GFMP' NESP1 $DOMTOT AL1 GRADAL1 (LIMAL1 * 0.0D0) RN1 GRADRN1 LIMRN1 VN1 GRADVN1 LIMVN1 PN1 GRADPN1 LIMPN1 ; ERRO = 'MAXIMUM' (PHFM '-' YNF) 'ABS' ; ERRO = ERRO '+' ( 'MAXIMUM' (PH0F '-' YNF) 'ABS' ) ; ERRO = ERRO '+' ( 'MAXIMUM' (ALF '-' PH1F) 'ABS' ) ; PHFM RN1F VN1F PN1F YNF1 ALF1 = 'PRET' 'GFMP' NESP1 $DOMTOT AL1 GRADAL1 (LIMAL1 * 0.0D0) RN1 GRADRN1 LIMRN1 VN1 GRADVN1 LIMVN1 PN1 GRADPN1 LIMPN1 ; ERRO = ERRO '+' ( 'MAXIMUM' (ALF1 '-' YNF) 'ABS' ) ; ERRO = ERRO '+' ( 'MAXIMUM' (ALF '-' YNF1) 'ABS' ) ; 'SI' (ERRO > 1.0D-13) ; 'MESSAGE' 'Multi-espece' ; 'MESSAGE' 'Erreur trop importante' ; 'MESSAGE' ('CHAINE' 'erro = ' ERRO) ; 'ERREUR' 5 ; 'FINSI' ; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales