Bethancourt L., Angel M.

Egami, Masahiro

Kyoto University Research Information Repository

TitlePARALLEL COMPUTATION OF COMPRESSIBLEMULTI-FLUID FLOW (Analysis and Numerical Simulation ofComplex Fluid Flow)Author(s)Bethancourt L., Angel M.; Egami, MasahiroCitation数理解析研究所講究録 (2011), 1724: 167-175Issue Date2011-01URL http://hdl.handle.net/2433/170460RightType Departmental Bulletin PaperTextversionpublisherKyoto UniversityPARALLEL COMPUTATION OF COMPRESSIBLE MULTI-FLUID FLOWAngel M. Bethancourt L.Institute ofComputational Fluid DynamicsMasahiro EgamiHonda R&D, Co. Ltd.Abstract. A collapse of a helium bubble due to its interactions with an incoming shockwave is numerically examined. The flow is modeled using the two-dimensionalcompressible Euler equations with an additional equation of motion for the propagatinginterface. The solution is obtained in a two-step manner, first, the Euler’s system ofequationsis solved using Roe’s flux difference splitting scheme; second, the moving interface (affinction of the specific heat ratio is used to mark the fluids interface) is captured using athird-order WENO scheme. For both, a third-order Runge-Kutta method is used for timeintegration. The solver is parallelized, allowing us to conduct fme mesh computations, andenabling the numerical reproduction of the mechanisms observed in previous publishedexperiments.1. IntroductionEarly attempts to compute multi-fluid flow, consisted in incorporating a conservative formofthe advection equation for the material interface into the Euler system of equation [1,2,3],and solve them using standard methods developed for the classical Euler equations.However, it was observed that spurious oscillations develop at the material interfaces due todifferences in fluid properties at the interface, as well as, errors in the location ofthe interface.In order to solve some ofthe problems, a quasi-conservative is proposed in [1], in which theadvection equation is decoupled $\theta om$ the Euler equations and written in non-conservativeform.Additional modifications to the goveming equations are made in [4,5], which render thealgorithm “single-fluid” when calculating the interface fluxes at the expense of conservation,to improve the coupling of the advection-Euler equations and to account for pressureequilibrium across the interface.In this paper, the quasi-conservative method described in [1,4] is implemented. Becauseour interest is on high-resolution computations, this implementation is extended using MPI(Message Passing Interface) libraries for its parallelization. As application, the interactionsbetween a helium bubble and a Mach 1.22 shock wave in air is studied.1724 2011 167-175 1672. Equations ofmotionFor the problem at hand, we limit our analysis to two immiscible ideal gases,$p=(\gamma-1)p$ , characterized by their constant ratio of specific heats, $\gamma$ . Thetwo-dimensional Euler equations for such multi-component flow in conservative form are$\frac{\partial Q}{\partial t}+\frac{\partial E}{\ }+ \frac{\partial F}{\partial y}=0$ , (1)$Q=\{\begin{array}{l}\rho\rho u\rho ve\end{array}\};E=\{\begin{array}{l}\rho upu^{2}+p\rho uv(e+p)u\end{array}\};F=\{\begin{array}{l}\rho v\rho uv\rho v^{2}+p(e+p)v\end{array}\}$where $p$ is the density, $u$ and $v$ are the velocity components, $p$ is the pressure, and $e$ is the totalenergy ofthe fluid.The moving interface is captured using the discontinuity in the ratio of the specific heatsfor the different components. Since material interfaces are advected by the flow, anyhnction of $\gamma$ obeys the advection equation [1],$\frac{\partial f(\gamma)}{\partial t}+u\frac{\partial f(\gamma)}{\ }+v \frac{\partial f(\gamma)}{\Phi}=0$ , . $\cdot$ . $f( \gamma)=\frac{1}{\gamma-1}$ (2)Spatial discretization of (1) is done using Roe’s flux difference splitting scheme. Thenumerical flux is given by:$E_{i+1/2}=E(Q_{L},Q_{R})= \frac{1}{2}[E(Q_{L})+E(Q_{R})]-\frac{1}{2}|A|[Q_{R}-Q_{L}]$ , where $A= \frac{\partial E}{\partial Q}$ (3)$A$ is the Jacobian matrix based on Roe’s average of $Q_{L}$ and $Q_{R}$ . The modification introducedin [4] consists in calculating two fluxes at the cell interface, by “freezing” the values of $\gamma$ ineither side ofthe interface:$E_{i+1/2}^{L}=E(Q_{L},Q_{R};\gamma_{L})$$E_{i+1/2}^{R}=E(Q_{L},Q_{R};\gamma_{R})$Then, use $E_{i+1/2}^{L}$ to update the flux in cell $i$ and $E_{i+1/2}^{R}$ to update the flux in cell $i+l$ .The spatial discretization of (2) is carried out using a third-order WENO scheme [6]. ATVD third order Runge-Kutta method is used to solve both the Euler and advection equations[6]. The equations are written as, $Q_{t}=L(u)$ , where $L(u)$ is the discretization of the spatialderivatives, and advanced in time according to the algorithm:168$Q^{(1)}=Q^{n}+\Delta tL(Q^{n})$$Q^{(2)}= \frac{3}{4}Q^{n}+\frac{1}{4}Q^{(1)}+\frac{1}{4}\Delta tL(Q^{(1)})$$Q^{n+1}= \frac{1}{3}Q^{n}+\frac{2}{3}Q^{(2)}+\frac{2}{3}\Delta tL(Q^{(2)})$(4)The parallelization of the code is achieved by using Domain Decomposition inconjunction with a SPMD (Single Program Multiple Data) technique. The MPI libraries areused to communicate between the different processors. Our system consists of 8 nodes, eachnode contains of $2cpu$-IntelXeon (8 cores). The maximum number ofcores available is 64.3. Problem DefinitionOur focus is on the dynamics ofthe interaction ofa Mach 1.22 shock and a helium bubble,as presented in [7,8]. The initial conditions for this problem are given in Fig. 1, the upperand lower boundaries are considered solid walls, the inflow is $fi\cdot om$ the right and propertiesare specified, at outflow (left) a zero gradient is imposed. A uniform grid with spatialresolution of 0.05 mm is used for the present computation, resulting in a grid size of712lxl781. It takes approximately $12 hr30\min$ to obtain a 528 ps sample (time step:$3.3333x10^{-8}s)$ using 32 processors.4. Computational ResultsFigure 2 shows a sequence of idealized Schlieren images [8], showing the overall wavestmcture. At 24 ps, Figure $2a$ shows a reffacted shock inside of the bubble, moving fasterthan incident shock, since the helium has a higher sound speed than the air. Upstream ofthebubble, a reflected wave is shown moving away $\theta om$ the bubble. At 36 ps, Figure $2b$ showsa four-shock configuration which it is usually termed twin regular reflection-reffaction (TRR),also shown in Figure 3. At 48 ps, Figure $2c$ shows the re$\theta acted$ wave emerging $\theta om$ thebubble and becoming the transmitted wave, also two cusps appear in the left-side of thebubble indicating intemally reflected waves that converge toward the axis of the bubble(Figure $2d,$ $56$ ps). At 64 ps, Figure $2e$ shows the intemal reflected waves diverging aftercrossing, and the two branches ofthe transmitted shock cross downstream ofthe bubble. At88 ps, Figure 2fthe initial reflected shock has reflected from the upper and bottom walls. Atlarger times (228 ps), Figures $2g$ shows that the bubble has ffirther deformed into a kidneyshaped object. Finally, Figures $2h$ and $2i$ show an induced jet of air along the axis of thebubble, resulting in the development oftwo vortical structures.Figure 3 shows the enlargement of the twin regular reflection-refraction (TRR) shockconfiguration as well as the locations used to compare the shock waves velocities. Table 1169shows the comparison of the present computation wave velocities with the experiment ofHass and Sturtevant [7] and Quick and Kami [8] computation. Overall agreement isexcellent, except for the values of the reflected and upstream velocities; this disagreement isdue to the existence of bubble contamination which affects the values of its physicalproperties.Figure 4 uses the density, vorticity and density gradient ffinction to visualize the flow.Clearly, the bulk of the vorticity is produced along the bubble interface [9], following theprogress of the reffacted shock inside the bubble. This is consistent with the source term$(S \equiv\frac{\nabla p\cross\nabla p}{\rho^{2}})$ appearing in vorticity equation ofan inviscid fluid. The motion produced bythe generation ofvorticity is the driven force inducing the jet along the axis ofthe bubble.5. ConclusionsThe goal of simulating high resolution compressible multi-fluid flow was achieved. Aquasi-conservative method is successffilly implemented in a parallel environment to modelthe interactions between a helium bubble and a shock wave. The main features ofthe floware captured, and compare quite well with previous numerical and experimental results.Further investigation will be aimed to consider flows with more stiff conditions.References[1] Abgrall R., 1996 How to Prevent Pressure Oscillations in Multicomponent FlowCalculations: A Quasi Conservative Approach, J Comp. Phys. 125, 150-160[2] Larrouturou, B., 1991, How to Preserve the Mass Fractions Positivity when ComputingCompressible Multi-Component Flows, J Comp. Phys. 95, 59-84[3] Mulder, W., Osher, S., and Sethian, J.A., 1992, Computing Interffice Motion inCompressible Gas Dynamics, J Comp. Phys. 100, 209-228[4] Abgrall, R., and Karni, S., 2001, Computations of Compressible Multifluids, J Comp.Phys. 169, 594-623[5] Fedkiw, R.P., Aslam, T., Berriman, B., and Osher, S., 1999, A non-oscillatory Eulerianapproach to interfaces in multimaterial flows (the Ghost Fluid Method), J Comp. Phys. 152,457-492[6] Shu, C.W., 1997, Essentially Non-Oscillatory and Weighted Essentially Non-OscillatorySchemes for Hyperbolic Conservation Laws, ICASE Report No. 97-65[7] Haas, J.-F., and Sturtevant, B. 1987 Interaction ofweak shock waves with cylindrical andspherical gas inhomogeneities. JFluidMech 181, 41-76.[8] Quirk, J. J., and Kami, S. 1996, On the dynamics of a shock-bubble interaction J FluidMech. 4129-163.170[9] Picone, J.M., and Boris, J.P., 1988, Vorticity generation by shock propagation throughbubbles in a gas, J FluidMech. 189, 23-51(a) Initial Conditions (356mm,89mm)$t\dot{m}\text{\’{e}} 000000\triangleleft p0\mathfrak{g}ld’ 121x179t$Figure 1. Computational domain and initial parameters.171(d) $56ps$ (g) $228\mu s$(e) $64ps$ (h) $412ps$$—–\vee\cdot\cdot.\cdot--$(i) $464ps$(f) $88ps$—-$\cdots$Figure 2. Idealized Schlieren images generated following the methodology outlined in [8](using the magnitude ofthe gradient ofthe density). a) 24 ps, b) 36 ps, c) 48 ps, d) 56 ps, e)64 ps, f) 88 ps, g) 228 ps, h) 412 ps, i) 464 ps. Times are estimates after the shock impactinto the bubble.172Figure 3. Twin regular reflection-reffaction (TRR) shock configuration as well as thelocations used to compare the shock waves velocities. $V_{S}$ : velocity ofincident shock. $V_{R}$ :velocity of refracted shock. $V_{T}$ : velocity of transmitted shock. $V_{ui}$ : velocity of upstreamedge of bubble, initial times. $V_{di}$ : velocity of downstream edge of bubble, initial times.$V_{j}$ : velocity ofjet head.Table 1. Comparison of shock waves velocities173$)$— $-$$\S$–$\overline{\overline{\underline{-}}\underline{g}\cdot.--}-\overline{-}-$$g^{\frac{1}{a^{m}\cdot-\frac{\sim}{}-}}\underline{-}i\underline{=.}-\cdots.--$$-\overline{\ovalbox{\tt\small REJECT}}$ $\overline{\ovalbox{\tt\small REJECT}}$ $—-\overline{\ovalbox{\tt\small REJECT}}_{j}^{*}$$\overline{\underline{-}}\underline{-\cdot\cdot}$. $\overline{-}---$ $\underline{-}-\cdot$ $-\overline{-}-$ $–\overline{g\mu\underline{\cdot}\cdot.\cdot}$ $—’-$$\overline{\underline{\epsilon}_{\ddot{r}^{-}-}}$ $rightarrow-$ $\underline{\epsilon}:\cdot-\cdot\cdot$ $-.-$ $—-\underline{\underline{-}}ia^{\infty}-$ $’-$ $\cdot\cdot$–Figure 4 (a-1). For caption see next page.174$-7...:\overline{\overline{-}}a^{m}..\cdot\cdot\wedge.-\cdot..\cdot\cdot\cdot\underline{}.\underline{-}j1^{*\cdot\cdot\cdots\cdot-\infty-},’.ae^{:}.\cdot..\cdot\cdot\cdot _{\dot{\mathfrak{t}}}\overline{-}\cdot-\overline{r-}arrow$ $|\underline{r_{\overline{\wedgearrow-\cdot-}\cdot:.-}\overline{g}}$.. –arrow\sim-$Figure 4. Flow visualization. Density, Vorticity and Density Function.175

PARALLEL COMPUTATION OF COMPRESSIBLE MULTI-FLUID FLOW (Analysis and Numerical Simulation of Complex Fluid Flow)

http://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/170460/1/1724-20.pdf

PARALLEL COMPUTATION OF COMPRESSIBLE MULTI-FLUID FLOW (Analysis and Numerical Simulation of Complex Fluid Flow)

Abstract

Similar works

Full text

Available Versions

Kyoto University Research Information Repository