* fichier : shock3d.dgibi ************************************************************************ * NAME : main3d.dgibi * DESCRIPTION : file to solve 3D tests: * 2. shock tube * * LANGUAGE : GIBIANE-CAST3M * AUTHOR : José R. Garcia Cascales * Universidad Politecnica de Cartagena, * jr.garcia@upct.es ************************************************************************ * VERSION : v1, 04/04/2002, version initiale * HISTORIQUE : vfinal, 22/11/2005 ************************************************************************ * * MESH DEFINITION * ************************************************************************ GRAPH = VRAI ; GRAPH = FAUX ; * ** number of elements in the z axis (after rotation) = NZ * NZ = 100 ; NX = 2; NY = 2; H = 10.D0 ; DZ = H/NZ ; L = NX*DZ ; P = NY*DZ ; xR = 0.5D0*L ; xL = -1.0D0*xR ; yR = 0.5D0*P ; yL = -1.0D0*yR ; P1 = XR YR 0. ; P2 = XL YR 0. ; P3 = XL YL 0. ; P4 = XR YL 0. ; P1P2 = 'DROIT' NX P1 P2 ; P2P3 = 'DROIT' NY P2 P3 ; P3P4 = 'DROIT' NX P3 P4 ; P4P1 = 'DROIT' NY P4 P1 ; BASEJR = 'DALLER' P1P2 P2P3 P3P4 P4P1 ; POSH = H/2. ; NEGH = (-1.)*POSH ; MIDDIV = NZ/2 ; V1 = 0. 0. POSH ; V2 = 0. 0. NEGH ; DOMTOT = DOM1 'ET' DOM2 ; MDOMTOT = 'MODELISER' DOMTOT 'EULER'; MDOM1 = 'MODELISER' DOM1 'EULER'; MDOM2 = 'MODELISER' DOM2 'EULER'; TDOMTOT = $DOMTOT. 'QUAF'; TDOM1 = $DOM1. 'QUAF'; TDOM2 = $DOM2. 'QUAF'; 'ELIMINATION' (TDOMTOT 'ET' TDOM1 'ET' TDOM2) 1D-6; * ******* Creation de la ligne Utilisée pour le Post-Traitement * reliant les points centres * XINIT = 0.5D0*DZ ; YINIT = 0.5D0*DZ ; ZINIT = NEGH '+' (0.5D0*DZ) ; XFIN = XINIT ; YFIN = YINIT ; ZFIN = POSH '-' (0.5D0*DZ) ; PINI = XINIT YINIT ZINIT; PFIN = XFIN YFIN ZFIN; DIV = NZ '-' 1 ; COURB = PINI 'DROIT' DIV PFIN; COURB = COURB 'COULEUR' 'VERT'; 'ELIMINATION' (1.0D-6) COURB $DOMTOT.CENTRE; 'SI' GRAPH ; 'TRACER' (DOMTOT 'ET' COURB) 'TITRE' 'Domaine total' ; 'FINSI' ; ************************************************************************ * * INITIAL CONDITIONS * ************************************************************************ CP = 2.d0 ; CVM = 0.d0 ; * * 'DISCRET') 'ET' 'DISCRET') ; 'UVZ' 0.d0 'NATURE' 'DISCRET') ; 'ULZ' 0.d0 'NATURE' 'DISCRET') ; 'DISCRET') 'ET' 'DISCRET') ; 'DISCRET') ; 'DISCRET'); 'DISCRET'); 'DISCRET'); 'DISCRET'); 'DISCRET'); 'DISCRET'); 'DISCRET'); * vch7 (rhov) = gammav*p/((gammav - 1.D0)*Cpv*T) * vch8 (rhol) = gammal*(p + pil)/((gammal - 1.D0)*Cpl*T) * IMPROVE, VCH7 & VCH8 ARE NOT ASSOCIATED TO ANY MESH!!!!!!!!! VCH7 = NUM1 '/' NUM2 ; NUM3 = 2.8d0 * (VCH4 '+' 8.5d8) ; VCH8 = NUM3 '/' NUM4 ; *************************************************************** ***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES ***** *************************************************************** 'DEBPROC' CONS ; 'ARGUMENT'VCH1*'CHPOINT' VCH2*'CHPOINT' VCH3*'CHPOINT' VCH4*'CHPOINT' VCH5*'CHPOINT' VCH6*'CHPOINT' VCH7*'CHPOINT' VCH8*'CHPOINT' ; AL = VCH1 '-' 1.0 ; AL = AL * (-1.0) ; 'UVY' 'UVZ') 'NATURE' 'DISCRET' ; 'ULZ')'NATURE' 'DISCRET' ; ECV = 0.5d0 * (UVX2 '+' UVY2 '+' UVY2) ; * ev = Cpv*v(5)/gammav * IEV = ('PSCAL' CPV VCH5 ('MOTS' 'SCAL') ('MOTS' 'SCAL'))'/' 1.4d0 ; IEV = (1008.d0 * VCH5)'/' 1.4d0 ; TIEV = IEV '+' ECV ; WCH5 = WCH1 * TIEV ; ECL = 0.5d0 * (ULX2 '+' ULY2 '+' ULZ2) ; * el = Cpl*v(6)/gammal + pil/rl * NUM9 = ('PSCAL' CPL VCH5 ('MOTS' 'SCAL') ('MOTS' 'SCAL'))'/' 2.8d0 ; NUM9 = (4186.d0 * VCH6)'/' 2.8d0 ; NUM0 = 8.5d8 '*'(VCH8 '**' -1.) ; IEL = NUM9 '+' NUM0 ; TIEL = IEL '+' ECL ; WCH6 = WCH2 * TIEL ; 'RESPRO' WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 ; 'FINPROC' ; *************************************************************** WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 = CONS VCH1 VCH2 VCH3 VCH4 VCH5 VCH6 VCH7 VCH8 ; * * We create old alpha, to satisfy the prim operator * necessities * OALP = VCH1 ; OTV = VCH5 ; OTL = VCH6 ; PINT RL RV TL TV P UL UV ALP = 'PRIM' 'TWOFLUID' WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 OALP OTV OTL CP CVM; ERRO = 'MAXIMUM' (VCH4 '-' P) 'ABS' ; 'SI' (ERRO > 1.0D-6) ; 'MESSAGE' 'Problem in the stic file!!!' 'ERREUR' 5 ; 'FINSI' ; * **** Plot of IC * 'SI' GRAPH ; MOD1 = 'MODELISER' ($DOMTOT . 'MAILLAGE' ) 'THERMIQUE'; 'TRACER' CHM_V1 MOD1 'TRACER' CHM_V2 MOD1 'TRACER' CHM_V3 MOD1 'TRACER' CHM_V4 MOD1 'TRACER' CHM_V5 MOD1 'TRACER' CHM_V6 MOD1 'TRACER' CHM_V7 MOD1 'TRACER' CHM_V8 MOD1 'FINSI' ; ************************************************************************* * * MAIN PROGRAM * ************************************************************************* * Different upwind scheme available * AUSM + * METO = 'AUSMP1' ; * Preconditioned AUSM + METO = 'AUSMP2'; * AUSMDV * METO = 'AUSMDV1'; * Preconditioned AUSMDV * METO = 'AUSMDV2'; * * Iterations * NITER = 10000000 ; * * Gravity, Final time, CFL-like parameter * g = 0.d0 ; TFINAL = 0.005d0 ; MCFL = 0.1D0 ; * * Order in space, order in time and a useful counter * ORDESP = 2; ORDTEM = 2; COUNT = 1; * * TIME IS RESETED TO 0 seconds * TPS = 0. ; * * Names of the residuum components * 'MESSAGE' ; 'MESSAGE' ('CHAINE' 'Methode = ' METO) ; 'MESSAGE' ; * * Some preliminary calculation for the evaluation of gradients * 1st. for the gradient of void fraction * 2nd. for the gradients in 2nd order * DALP LIMCH GRG = 'PENT' MDOMTOT 'CENTRE' 'SI' (ORDESP 'EGA' 2) ; GRADAL CACA COEFSCAL = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' GRADUV CACA COEFVECT = 'PENT' MDOMTOT 'CENTRE' 'EULEVECT' 'FINSI' ; * **** Temporal loop * 'REPETER' BL1 NITER ; * **** Evaluation of variables at each side of the interface (PRET operator) * 'SI' ((ORDESP 'EGA' 1) 'OU' (&BL1 'EGA' 1)) ; RLF RVF TLF TVF PF ULF UVF ALPF = 'PRET' 'TWOFLUID' 1 1 $DOMTOT ALP UV UL P TV TL RV RL ; 'SINON' ; GRADAL LIMALP = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADUV LIMUV = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADUL LIMUL = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADP LIMP = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADTV LIMTV = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADTL LIMTL = 'PENT' MDOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' RLF RVF TLF TVF PF ULF UVF ALPF = 'PRET' 'TWOFLUID' ORDESP ORDTEM $DOMTOT ALP GRADAL LIMALP UV GRADUV LIMUV UL GRADUL LIMUL P GRADP LIMP TV GRADTV LIMTV TL GRADTL LIMTL RV RL DELTAT ; 'FINSI' ; * **** Evaluation of Residual and Dt * $DOMTOT LISTINCO ALPF UVF ULF PF TVF TLF RVF RLF ; DT = MCFL '*' DELTAT ; 'SI' (COUNT 'EGA' 1 ) ; COUNT = 0 ; OALP = ALP ; OLDDT = DT ; 'FINSI' ; * **** Increment of the variables (convection) * RESIDU = DT '*' RESIDU ; 'UVY' 'UVZ') ; 'ULY' 'ULZ') ; * * Analizamos la conservacion de la masa * * **** Previous update of conserved variables * TPS = TPS '+' DT ; WCH1 = WCH1 '+' DRVN ; WCH2 = WCH2 '+' DRLN ; WCH3 = WCH3 '+' DRUVN ; WCH4 = WCH4 '+' DRULN ; WCH5 = WCH5 '+' DRETVN ; WCH6 = WCH6 '+' DRETLN ; ************************************************************************ **** Adittion of the source terms, final update of conserved variables ************************************************************************ DALP1 LIMCH = 'PENT' MDOMTOT 'CENTRE' DADP = DT * PINT * DALP1 ; * * Source term corresponding to the vapour momentum equation * 'NATURE' 'DISCRET'); S32 = DT * ALP * RV * g ; S3 = S31 'ET' S321 ; 'UVZ') 'NATURE' 'DISCRET' ; SO3 = (S3 '+' DADP) ; 'UVZ') 'NATURE' 'DISCRET' ; WCH3 = WCH3 '+' SO3 ; * Source term corresponding to the liquid momentum equation 'NATURE' 'DISCRET'); S42 = DT *(ALP '-' 1.d0) * RL * g ; S42 = S42 * (-1.) ; S4 = S41 'ET' S421 ; 'ULZ') 'NATURE' 'DISCRET' ; SO4 = (S4 '-' DADP) ; 'ULZ') 'NATURE' 'DISCRET' ; WCH4 = WCH4 '+' SO4 ; * Source term corresponding to the vapour total energy equation S5 = S32 * VY ; DADT = DT * PINT * (ALP '-' OALP) '/' OLDDT ; SO5 = (S5 '-' DADT) ; WCH5 = WCH5 '+' SO5 ; * Source term corresponding to the liquid total equation S6 = S42 * LY ; SO6 = (S6 '+' DADT) ; WCH6 = WCH6 '+' SO6 ; OALP = ALP ; OLDDT = DT ; OTV = TV ; OTL = TL ; ************************************************************************ **** Updating primitive variables ************************************************************************ PINT RL RV TL TV P UL UV ALP = 'PRIM' 'TWOFLUID' WCH1 WCH2 WCH3 WCH4 WCH5 WCH6 OALP OTV OTL CP CVM; * ANADIR CVM TRAS COMPILAR ************************************************************************ **** Graphical outputs ************************************************************************ 'SI' (((&BL1 '/' 1000) '*' 1000) 'EGA' &BL1) ; 'MESSAGE' &BL1 TPS ; 'FINSI' ; 'SI' (TPS > TFINAL) ; 'MESSAGE' &BL1 TPS ; 'SI' GRAPH ; 'alpha' ; 'uvx (m/s)' ; 'uvy (m/s)' ; 'uvz (m/s)' ; 'ulx (m/s)' ; 'uly (m/s)' ; 'ulz (m/s)' ; 'p (Pa)' ; 'Tv (K)' ; 'Tl (K)' ; 'rv (kg/m3)' ; 'rl (kg/m3)' ; 'FINSI' ; 'QUITTER' BL1 ; 'FINSI' ; 'FIN' BL1 ; ************************************************************************ **** End of loop ************************************************************************ * 'SI' GRAPH ; MOD1 = 'MODELISER' ($DOMTOT . 'MAILLAGE' ) 'THERMIQUE'; 'TRACER' CHM_V1 MOD1 'TRACER' CHM_V2 MOD1 'TRACER' CHM_V3 MOD1 'TRACER' CHM_V4 MOD1 'TRACER' CHM_V5 MOD1 'TRACER' CHM_V6 MOD1 'TRACER' CHM_V7 MOD1 'TRACER' CHM_V8 MOD1 'FINSI' ; ************************************************************************ *** saving the results in a file ************************************************************************ * 'OPTION' 'SAUVER' ('CHAINE' 'shock3d.sauv') ; * 'SAUVER' ; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales