* fichier : consistence1_HUSVL.dgibi *********************************************************** **** APPROCHE VF "Cell-Centred Formulation" pour la **** **** solution des **** **** Equations d'Euler pour un gaz parfait **** **** **** **** Gaz monoespece "calorically perfect" **** **** Consistence, methode HUSVL **** **** **** **** A. BECCANTINI DRN/DMT/SEMT/TTMF NOVEMBRE 1998 **** *********************************************************** 'OPTION' 'DIME' 2 ; 'OPTION' 'ELEM' 'QUA4' ; 'OPTION' 'ISOV' 'SULI' ; 'OPTION' 'ECHO' 0 ; 'OPTION' 'TRAC' 'X'; GRAPH = FAUX ; * GRAPH = VRAI ; * *** Methodes possibles : * * 'GODUNOV' * 'VANLEER' * 'VLH' * 'HUSVL' * 'HUSVLH' * METO = 'HUSVL' ; ************ * MAILLAGE * ************ NX = 5 ; NY = 1 ; L = 1.0D0 ; DX = L/NX/2.0D0 ; H = DX ; xD = 0.5D0*L ; xG = -1.0D0*xD ; yH = 0.5D0*H ; yB = -1.0D0*yH ; A1 = xG yB ; A2 = 0.0D0 yB ; A3 = xD yB ; A4 = xD yH ; A5 = 0.0D0 yH ; A6 = xG yH ; VECTG = XG 0.0D0 ; VECTD = XD 0.0D0 ; xBG = xG-DX ; XBD = xD+DX ; B1 = xBG yB ; B2 = xBG yH ; B3 = xBD yB ; B4 = xBD yH ; BAS1 = A1 'DROI' NX A2 ; BAS2 = A2 'DROI' NX A3 ; HAU2 = A4 'DROI' NX A5 ; HAU1 = A5 'DROI' NX A6 ; LATG = B1 'DROI' NY B2 ; LAT1 = A1 'DROI' NY A6 ; LAT12 = A2 'DROI' NY A5 ; LAT2 = A3 'DROI' NY A4 ; LATD = B3 'DROI' NY B4 ; DOM1 = LAT12 'TRANSLATION' NX VECTG ; DOM2 = LAT12 'TRANSLATION' NX VECTD ; VECTFG = (-1.0D0*DX) 0.0D0; VECTFD = DX 0.0D0; VECTBH = 0.0D0 (yH '-' yB); VECTHB = 0.0D0 (yB '-' yH); FRONTG1 = LAT1 'TRANSLATION' 1 VECTFG ; DCEL1 = FRONTG1 'ET' DOM1; FRONTG2 = DCEL1 'PLUS' VECTBH; FRONTG3 = DCEL1 'PLUS' VECTHB; FRONTG = (FRONTG1 'ET' FRONTG2 'ET' FRONTG3) 'COULEUR' 'ROUG'; 'ELIMINATION' 0.0001 FRONTG ; FRONTD1 = LAT2 'TRANSLATION' 1 VECTFD; DCEL2 = FRONTD1 'ET' DOM2; FRONTD2 = DCEL2 'PLUS' VECTBH; FRONTD3 = DCEL2 'PLUS' VECTHB; FRONTD = (FRONTD1 'ET' FRONTD2 'ET' FRONTD3) 'COULEUR' 'VERT'; 'ELIMINATION' 0.0001 FRONTD ; * *** Rotation * ANGLE = 30.0D0; ORIG = 0.0D0 0.0D0; 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'ANGLE = ' ANGLE); 'MESSAGE' ; DOM1 = DOM1 'TOURNER' ANGLE ORIG; DOM2 = DOM2 'TOURNER' ANGLE ORIG; FRONTG = FRONTG 'TOURNER' ANGLE ORIG; FRONTD = FRONTD 'TOURNER' ANGLE ORIG; DOMINT = DOM1 'ET' DOM2 ; 'ELIMINATION' DOMINT 1D-6 ; FRONT = FRONTG ET FRONTD ; 'ELIMINATION' FRONT 1D-6 ; DOMTOT = DOMINT ET FRONT ; 'ELIMINATION' DOMTOT 1D-6 ; 'SI' GRAPH; 'TRACER' DOMTOT 'TITRE' 'Maillage' ; 'FINSI' ; $DOMTOT = 'MODE' DOMTOT 'EULER' ; $DOMINT = 'MODE' DOMINT 'EULER' ; $DOM1 = 'MODE' DOM1 'EULER' ; $DOM2 = 'MODE' DOM2 'EULER' ; $FRONTG = 'MODE' FRONTG 'EULER' ; $FRONTD = 'MODE' FRONTD 'EULER' ; $FRONT = 'MODE' FRONT 'EULER' ; TDOMTOT = 'DOMA' $DOMTOT 'VF' ; TDOMINT = 'DOMA' $DOMINT 'VF' ; TDOM1 = 'DOMA' $DOM1 'VF' ; TDOM2 = 'DOMA' $DOM2 'VF' ; TFRONTG = 'DOMA' $FRONTG 'VF' ; TFRONTD = 'DOMA' $FRONTD 'VF' ; TFRONT = 'DOMA' $FRONT 'VF' ; MDOMTOT = TDOMTOT . 'QUAF' ; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MDOMINT = TDOMINT . 'QUAF' ; MFRONTG = TFRONTG . 'QUAF' ; MFRONTD = TFRONTD . 'QUAF' ; MFRONT = TFRONT . 'QUAF' ; 'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2 'ET' MDOMINT 'ET' MFRONTG 'ET' MFRONTD 'ET' MFRONT) 1.E-5 ; ************************************************************ * CONDITIONS INITIALES ET LIMITES. * ************************************************************ gamgd = 1.4D0; * *** Etat gauche * rog = 1.11; ung = 21.3; utg = 11.0; pg = 1234.7; rouxg = ((ung '*' ('COS' ANGLE)) '-' (utg '*' ('SIN' ANGLE))) '*' rog ; rouyg = ((ung '*' ('SIN' ANGLE)) '+' (utg '*' ('COS' ANGLE))) '*' rog; recing = 0.5D0 '*' rog '*' ((ung '*' ung) '+' (utg '*' utg)); retg = (pg '/' (gamgd '-' 1.0)) '+' recing; * *** Etat droite * rod = 1.11; und = 21.3; utd = 11.0; pd = 1234.7; rouxd = ((und '*' ('COS' ANGLE)) '-' (utd '*' ('SIN' ANGLE))) '*' rod; rouyd = ((und '*' ('SIN' ANGLE)) '+' (utd '*' ('COS' ANGLE))) '*' rod; recind = 0.5D0 '*' rod '*' ((und '*' und) '+' (utd '*' utd)); retd = (pd '/' (gamgd '-' 1.0)) '+' recind; * *** gamma * GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' gamgd; * *** ro * RO_f1 = 'KCHT' $FRONTG 'SCAL' 'CENTRE' rog ; RO_f2 = 'KCHT' $FRONTD 'SCAL' 'CENTRE' rod ; RO_f = 'KCHT' $FRONT 'SCAL' 'CENTRE' (RO_f1 'ET' RO_f2) ; RO_i1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' rog; RO_i2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' rod; RO_i = 'KCHT' $DOMINT 'SCAL' 'CENTRE' (RO_i1 'ET' RO_i2); RN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RO_i 'ET' RO_f) ; * *** ro u, ro v * GN_f1 = 'KCHT' $FRONTG 'VECT' 'CENTRE' (rouxg rouyg); GN_f2 = 'KCHT' $FRONTD 'VECT' 'CENTRE' (rouxd rouyd); GN_f = 'KCHT' $FRONT 'VECT' 'CENTRE' (GN_f1 'ET' GN_f2); GN_i1 = 'KCHT' $DOM1 'VECT' 'CENTRE' (rouxg rouyg); GN_i2 = 'KCHT' $DOM2 'VECT' 'CENTRE' (rouxd rouyd); GN_i = 'KCHT' $DOMINT 'VECT' 'CENTRE' (GN_i1 'ET' GN_i2); GN = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (GN_i 'ET' GN_f) ; * *** ro e * RE_f1 = 'KCHT' $FRONTG 'SCAL' 'CENTRE' retg ; RE_f2 = 'KCHT' $FRONTD 'SCAL' 'CENTRE' retd ; RE_f = 'KCHT' $FRONT 'SCAL' 'CENTRE' (RE_f1 'ET' RE_f2) ; RE_i1 = 'KCHT' $DOM1 'SCAL' 'CENTRE' retg ; RE_i2 = 'KCHT' $DOM2 'SCAL' 'CENTRE' retd ; RE_i = 'KCHT' $DOMINT 'SCAL' 'CENTRE' (RE_i1 'ET' RE_i2); RETN = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' (RE_i ET RE_f) ; ******************************************************** *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM *** ******************************************************** MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE'; * *** GRAPHIQUE DES C.I. * 'SI' GRAPH; * *** CREATION DE CHAMELEM * CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN; CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ; CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN; TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RHO at t=' 0.0); TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'RHOET at t=' 0.0); TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'RHOU at t=' 0.0); 'FINSI' ; ZERO =1.0D-12 ; TFINAL = 100. ; NITER = 50 ; SAFFAC = 0.5 ; * RN0 = 'COPIER' RN ; GN0 = 'COPIER' GN ; RETN0 = 'COPIER' RETN ; * TPS = 0. ; 'REPETER' BL1 NITER ; * **** Primitive variables * VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMMA ; ROF VITF PF GAMF = 'PRET' 'PERFMONO' 1 1 $DOMTOT RN VN PN GAMMA; RESIDU DELTAT = 'KONV' 'VF' 'PERFMONO' 'RESI' METO $DOMTOT ('MOTS' 'RN' 'RUX' 'RUY' 'RETN') ROF VITF PF GAMF ; DT_CON = SAFFAC '*' DELTAT ; * **** The time step linked to tps * DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ; * **** Total time step * DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ; * **** Increment of the variables (convection) * RESIDU = DTMIN '*' RESIDU ; * DRN = 'REDU' ('EXCO' 'RN' RESIDU 'SCAL') ('DOMA' $DOMINT 'CENTRE'); DGN = 'REDU' ('EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY')) ('DOMA' $DOMINT 'CENTRE') ; DRETN = 'REDU' ('EXCO' 'RETN' RESIDU 'SCAL') ('DOMA' $DOMINT 'CENTRE') ; * TPS = TPS '+' DTMIN ; RN = RN '+' DRN ; GN = GN '+' DGN ; RETN = RETN '+' DRETN ; * 'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ; 'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ; 'FINSI' ; 'SI' (TPS > TFINAL) ; 'QUITTER' BL1 ; 'FINSI' ; 'FIN' BL1 ; * *** GRAPHIQUE DES SOLUTIONS * 'SI' (GRAPH); * *** CREATION DE CHAMELEM * CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ; CHM_GN = 'KCHA' $DOMTOT 'CHAM' GN ; CHM_RETN = 'KCHA' $DOMTOT 'CHAM' RETN ; TRAC CHM_RN MOD1 'TITR' ('CHAINE' 'RHO at t=' TPS); TRAC CHM_RETN MOD1 'TITR' ('CHAINE' 'RHOET at t=' TPS); TRAC CHM_GN MOD1 'TITR' ('CHAINE' 'RHOU at t=' TPS); 'FINSI' ; * *** Critere de control * ERRO = 'MAXIMUM' (RN '-' RN0) 'ABS' ; 'SI' (ERRO '>' 1.0D-6) ; 'ERREUR' 5 ; 'FINSI' ; 'FIN' ;