A Simple Unified Coordinate Method for Compressible
Homogeneous Two-Phase Flows
Keh-Ming Shyue
1. Introduction
Our goal is to present a simple moving grid approach for the efficient numerical resolution of compressible multiphase flows with real materials in more than one space dimension. We consider a model homogeneous two-phase flow problem as an example, where the constitutive law for the thermodynamic property of the fluid component of interest is taken to fulfill a stiffened gas equation of state of the form
(1.1) p(ρ, e) = (γ − 1) ρe − γB
(cf. [10] and references therein). Here p, ρ, e, γ, and B are in turn the pressure, the density, the specific internal energy, the adiabatic constant (γ > 1), and the reference pressure. We use the Euler equations of gas dynamics as a model system for the principal motion of each fluid component in that the conservation form of the equations in an Nd-dimensional Cartesian coordinate ~X = (x1, x2, . . . , xNd) can
be written as (1.2) ∂ ∂t ρ ρui E + Nd X j=1 ∂ ∂xj ρuj ρuiuj+ pδij Euj+ puj = 0 for i = 1, 2, . . . , Nd.
Here ui, E, and δij denote the particle velocity in the xi-direction, the total energy, and the Kronecker delta, respectively. We have E = ρe +PNd
j=1ρu 2
j/2 as usual. To solve this homogeneous two-phase flow problem numerically, we want to use an extension of a shock-capturing method based on unified coordinate systems proposed previously by Hui and coworkers for single-phase flows (cf. [5, 6, 7]). This method belongs to a class of moving grid methods in which the temporal evolution of the physical grid coordinate ~X is assumed to satisfy the Lagrange-like condition
(1.3) ∂ ~X
∂t = h~u,
1991 Mathematics Subject Classification. Primary 65M06, 76L05; Secondary 76M20, 76T05. Key words and phrases. Unified coordinate method, Compressible two-phase flows, Wave-propagation method.
The author was supported in part by National Science Council of Republic of China Grant #96-2115-M-002-008. The final version of this paper will be submitted for publication elsewhere.
c
0000 (copyright holder)
where ~u = (u1, u2, . . . , uN d) is the vector of the particle velocity, and h ∈ [0, 1] is a freely chosen parameter that takes on the value of zero when we have a fixed Eulerian grid, and on the value of unity when we have a Lagrangian grid. We employ a fluid-mixture formulation of equations that are written in the generalized curvilinear coordinate as a basis to the modelling of the numerically induced mixing between two different fluid components within a grid cell. In addition to that, a set of geometric conservation laws is included to describe the basic motion of the grid metrics associated with the coordinate transformation between the physical and computational grid. This will be described further in Section 2.
In the current implementation of the unified coordinate method for two-phase flows, we use a finite volume method based on an f -wave decomposition viewpoint developed by Bale et al. [2] with the dimensional-splitting technique incorporated in the method for multidimensional problems. The method is a variant of the Godunov scheme, and has been widely used in many applications including the single-phase fluid of ideal gases on general quadrilateral grids [9], and hyperbolic systems on curved manifolds [12] and on spheres [11]. Numerical results presented in Section 3 gives some indications of the feasibility of the proposed method for two-phase flow problems in one and two space dimensions.
2. Mathematical models in unified coordinates
We begin our discussion by reviewing briefly the mathematical model for single-phase flow problems in a unified-coordinate moving grid system as proposed in [7]. For simplicity, we consider the two-dimensional case Nd = 2 as an example, and having (1.3) in mind we introduce a coordinate mapping from the physical domain (t, x1, x2) to the computational domain (τ, ξ1, ξ2) as
dt = dτ,
dx1= hu1dτ + a11dξ1+ a12dξ2, dx2= hu2dτ + a21dξ1+ a22dξ2. (2.1)
Here a11, a12, a21, and a22 are the metric terms of the mapping. Note that if we define Dh/Dt as being the material derivative following the pseudo-particle, whose velocity is h(u1, u2), i.e.,
Dh Dt := ∂ ∂t+ hu1 ∂ ∂x1 + hu2 ∂ ∂x2 , from (2.1) we would have
Dhξ1
Dt = 0,
Dhξ2
Dt = 0,
yielding coordinate lines ξ1 and ξ2 to be material functions of the pseudo-particle. Hence, computationally grid cells in this coordinate system would move and change shape with pseudo-particles rather than with fluid particles as in Lagrangian coor-dinates.
Now under the mapping (2.1), Euler Eqs. (1.2) can be transformed into the new coordinate system as
(2.2) ∂ ∂τ ρJ ρJu1 ρJu2 EJ + ∂ ∂ξ1 ρJU1 ρu1JU1+ pa22 ρu2JU1−pa21 EJU1+ pJV1 + ∂ ∂ξ2 ρJU2 ρu1JU2−pa12 ρu2JU2+ pa11 EJU2+ pJV2 = 0,
where J = a11a22−a12a21is the Jacobian of the mapping, V1= (u1a22−u2a21)/J, V2 = (u2a11−u1a12)/J are the velocity component normal to the ξ1- and ξ2 -direction, respectively, and Ui = (1 − h)Vi is the contravariant velocity in the ξi-direction. Note that in deriving (2.2) which is written in a strong conservation law form we have assumed and imposed the following metric identities implicitly,
∂a22 ∂ξ1 −∂a12 ∂ξ2 = 0, (2.3) ∂a21 ∂ξ1 −∂a11 ∂ξ2 = 0, (2.4) ∂J ∂τ − ∂ ∂ξ1 (hV1) − ∂ ∂ξ2 (hV2) = 0. (2.5)
It is known in the literature that the first two identities constitute a differential statement of space conservation for a closed grid cell, and the last identity expresses volume conservation. The combination of these three identities can be referred as the geometric conservation laws for a moving grid system (cf. [16, 18]).
Here as in the original unified coordinate method and its variant (cf. [3, 5, 7]), our approach to model the aforementioned geometric conservation laws, in particular to (2.5), is based on the compatibility condition of the mixed second order partial derivatives ∂2x
i/∂τ ∂ξj and ∂2xi/∂ξj∂τ , yielding a set of equations for the time variation of each of the metrics aij as follows
(2.6) ∂aij
∂τ =
∂ ∂ξj
(hui) for i, j = 1, 2.
It is certain that, with a prescribed region of the physical domain (rectangular or not), the initial grid system and also the initial condition of the above equations can be obtained by a chosen technique from numerical grid generators (cf. [17]) when it is necessary. Once the solutions of (2.6) are known during a time step, it is an easy manner to compute J from an explicit formula as defined above.
In summary, if we have chosen a fixed parameter h in the mathematical for-mulation of model equations, but have not taken the more general case where it is constrained by an additional condition (cf. [5, 7]), the governing equations to be solved for single-phase flows in a unified-coordinate method consist of the physical part (2.2) and the geometrical part (2.6).
To generalize the above single-phase flow model to the present homogeneous two-phase flow, as in many of the diffused-interface type algorithms developed in the literature (cf. [1, 8]), we need to devise a proper model and treatment of numerical mixing so that both the pressure and velocity remain in equilibrium without introducing any spurious oscillations for grid cells near interfaces where two or more fluid components are mixed. The approach we take here follows essentially the same idea as discussed in [14], and is an extension of the fluid-mixture type algorithm from a fixed Cartesian grid to a moving quadrilateral grid. That is to say, we use (2.2) as a model system that describes the motion of the fluid mixtures of the conserved variables in a multicomponent grid cell, and derive a set of effective equations for the problem-dependent material functions in those cells so that the pressure can be computed easily from the equation of state.
To find out the aforementioned effective equations for the mixture of material quantities in the present stiffened gas case, it is usual to start with an interface-only problem where both the pressure and each component of the particle velocities
are constant in the domain, while the other variables such as the density and the material quantities are having jumps across some interfaces. In this instance, from (2.6) we observe easily the time invariance of the grid metrics aij as
∂aij
∂τ = 0 for i, j = 1, 2,
and from (2.2) we obtain equations for the density and total internal energy as ∂ρ ∂τ + U1 ∂ρ ∂ξ1 + U2 ∂ρ ∂ξ2 = 0, ∂ρe ∂τ + U1 ∂ρe ∂ξ1 + U2 ∂ρe ∂ξ2 = 0, (2.7)
in a respective manner. By inserting the equation of state (1.1) into the latter one, we find an alternative description of the energy equation
(2.8) ∂ ∂τ p + γB γ − 1 + U1 ∂ ∂ξ1 p + γB γ − 1 + U2 ∂ ∂ξ2 p + γB γ − 1 = 0
that is in relation to the pressure and also the material quantities γ and B. To maintain the pressure in equilibrium as it should be for this interface-only problem, as before we split (2.8) into the following two equations for the fluid mixtures 1/(γ − 1) and γB/(γ − 1) as ∂ ∂τ 1 γ − 1 + U1 ∂ ∂ξ1 1 γ − 1 + U2 ∂ ∂ξ2 1 γ − 1 = 0, (2.9) ∂ ∂τ γB γ − 1 + U1 ∂ ∂ξ1 γB γ − 1 + U2 ∂ ∂ξ2 γB γ − 1 = 0, (2.10)
respectively. Combining the above two equations to (2.2) and (2.6) yields a so-called γ-based model system that is fundamental in our unified coordinate method for describing the behavior of the numerical mixing between different fluid components near the interface. With that, there is no difficulty to compute the pressure based on the equation of state,
p = " E − P2 j=1(ρuj)2 2ρ − γB γ − 1 # 1 γ − 1 .
We note that analogously in our previous work (cf. [14]), the above γ-based model can be reformulated into a slightly simple volume-fraction model, but to save some space the result of that will not be presented here. It should be noted also that while the approach we have taken here may not be the best way to model the geometric conservation law (2.5), see [15] for more details, it is a simple way to be coped with our current implementation of a dimension-by-dimension splitting method for numerical approximation of multi-dimensional problems, and yet can give reasonable results for some sample problems as considered in Section 3.
3. Numerical examples
To find approximate solutions of our two-phase flow model in unified coordinate systems, i.e., a combined equations with the physical part: (2.2), (2.9), (2.10), and the geometric part: (2.6), we use a flux-based wave decomposition method developed by Bale et al. [2] with the dimensional-splitting technique incorporated in the method for multidimensional problems. This method is a variant of the
standard wave-propagation scheme [9] in that we solve one-dimensional Riemann problems in the direction normal to each cell interface as usual. Rather than using the resulting jumps of the state variables moving at constant speeds to update the solutions in neighboring grid cells, we use the resulting jumps of fluxes for that instead, which has shown to be a much more accurate approach than the former one for a class of hyperbolic problems with spatially varying fluxes.
We now present some sample numerical results obtained using our unified co-ordinate method for single- and two-phase flow problems, see [15] for additional results. Without stated otherwise, we have carried out all the tests using the Courant number ν = 0.5, and the minmod limiter in the high-resolution version of the method. The material-dependent parameters we use are (γ, B) = (1.4, 0) and (4.4, 6 × 108Pa) for the gas- and liquid-phase, respectively.
Example 3.1. To begin with, we consider a one-dimensional complex wave interactions problem studied by Woodward and Colella [19]. That is, for the initial condition, we have three constant states with data
ρ u1 p L = 1 0 103 , ρ u1 p M = 1 0 10−2 , ρ u1 p R = 1 0 102 ,
where L is the state used for x1 ∈[0, 0.1), M is the state used for x1 ∈[0.1, 0.9), and R is the state used for x1∈[0.9, 1]. Here the fluid is a perfect gas in the entire domain, and there are two solid walls at the left and right boundaries.
In this setup, it is known that, after breaking the membranes at x1 = 0.1 and 0.9, a shock wave, contact discontinuity, and rarefaction wave develop at each discontinuity individually. As time progresses, the shock waves move toward each other and then collide, yielding a new contact discontinuity from the collision. Further collisions then occur.
Figure 1 shows numerical results for the density, velocity, and pressure at three different times t = 0.016, 0.032, and 0.038 with a 200 mesh points. From the figure, we observe a sensible improvement of the contact discontinuities, especially, the emergent one after the head-on collision of the initial shock waves, when our unified coordinate results with h = 0.99 are in comparison with the Eulerian results. In addition, we find good resolution of the shock wave and rarefaction as we compare our numerical results with the fine grid solutions obtained using the Eulerian grid approach with a 2000 mesh points. The physical grid coordinates generated by our unified coordinate method are plotted in Fig. 2 at eleven different times t = i × 0.0038, for i = 0, 1, . . . , 10, noticing interesting numerical temporal behavior of the grid movement that follows closely to the main feature of the underlying flow. Example3.2. We are next concerned with a two-dimensional Riemann prob-lem in that the initial condition in the four quadrants of a unit square domain is composed of four shock waves with data
(ρ, u1, u2, p)1= (ρ0, 0, 0, p0) , (ρ, u1, u2, p)2= (¯ρ, ¯u1, 0, ¯p) , (ρ, u1, u2, p)3= (ρ0, ¯u1, ¯u2, p0) , (ρ, u1, u2, p)4= (¯ρ, 0, ¯u2, ¯p) , where ρ0 = 1.1, p0 = 1.1, ¯ρ = 0.5065, ¯u1 = 0.8939, ¯u2 = 0.8939, and ¯p = 0.35 (cf. [13]). Here we again consider a single component fluid with a perfect gas in the whole domain, while there are non-reflecting boundary conditions on all sides.
0 0.5 1 0 2 4 6 8 ρ 0 0.5 1 −10 0 10 20 0 0.5 1 0 100 200 300 400 u1 p t = 0.016 t = 0.016 t = 0.016 0 0.5 1 0 5 10 15 20 ρ Fine grid h=0 h=0.99 0 0.5 1 −5 0 5 10 15 0 0.5 1 0 200 400 600 u1 p t = 0.032 t = 0.032 t = 0.032 0 0.5 1 0 2 4 6 8 ρ 0 0.5 1 −5 0 5 10 15 0 0.5 1 0 200 400 600 u1 p x1 x1 x1 t = 0.038 t = 0.038 t = 0.038
Figure 1. Numerical results for the Woodward-Colella problem at three different times t = 0.016, 0.032, and 0.038. Here the solid line is the fine grid solution obtained using h = 0 with 2000 meshes, and the dotted and triangular points are the computed solutions obtained using h = 0 and 0.99, respectively, with 200 meshes.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.01 0.02 0.03 0.04 x1 ti m e
Figure 2. Physical grid coordinates for the unified coordinate run shown in Fig. 1 at eleven different times t = i × 0.0038 for i = 0, 1, . . . , 10. Each little dashed line displayed in the graph gives a cell-center location of the grid system at a specific output time.
As in Example 3.1 we run this problem using two different grid-movement parameters h = 0 and 0.99, and in Fig. 3 show contour plots of the density and pressure, and the physical grid system at time t = 0.2 with a 200 × 200 grid. From the figure, it is interesting to see that the collisions between the initial shock waves creates an oval shape region bounded by the incident and reflected shock waves.
0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Density Pressure Physical grid
a) h = 0.99 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Density Pressure Physical grid
b) h = 0
Figure 3. Numerical results for a 4-shocks two-dimensional Rie-mann problem. Contours of the density and pressure, and the physical grid coordinates are shown at time t = 0.2 obtained using h = 0 and 0.99 with a 200 × 200 grid. Here the computed velocity field is superimposed into the pressure contours. For clarity, solu-tion coarsening factors of 10 and 4 in each x1- and x2-direction, respectively, are used to graph the velocity and the physical grid.
When we make a comparison of our results with the one appeared in Fig. 6 of [13], for instance, (which was done using a state-of-the-art second order Eulerian method with a 400×400 grids), we find notably a better resolution of our unified coordinate result than the Eulerian one for the slip lines that are situated in the shock-waves bounded oval region, and similar solution behaviors for the location and structure of the reflected shock waves.
Example 3.3. As an example to show how our method works on two-phase
flows, we are interested in a model underwater explosion problem (cf. [4]). In this test, we take a rectangular domain (x1, x2) ∈ [−2, 2] × [−1.5, 1]m2, and employ the initial condition that consists of a stationary horizontal air-water interface at the x2= 0 axis and a circular gas bubble in water with the center (x01, x
0
2) = (0, −0.3)m and of the radius r0 = 0.12m. Here above the air-water interface, the fluid is a perfect gas at the standard atmospheric condition,
(ρ, u1, u2, p) =
1.2 kg/m3, 0, 0, 105 Pa.
Below the air-water interface, in region inside the gas bubble the fluid is modeled as a perfect gas also with the state variables
(ρ, u1, u2, p) =
1250 kg/m3, 0, 0, 109 Pa,
and in region outside the gas bubble the fluid is water with the state variables (ρ, u1, u2, p) =
103kg/m3, 0, 0, 105Pa.
The boundary conditions considered here are solid wall on the left, right, and bottom sides, and non-reflecting on the top side of the domain.
Note that due to the pressure difference between the fluids at r = r0, breaking of this underwater bubble would results in an outward-going shock wave in water, an inward-going rarefaction wave in gas, and a material interface lying in between that separates the gas and water. Soon after this shock wave is diffracted through the nearby flat air-water interface, it is known in the literature (cf. [4]) that the topology of the underwater bubble will undergo a change from the original circular-shape to an oval-like circular-shape. As time evolves, this gas bubble would continue rising upward, causing the subsequent deformation of the horizontal air-water interface.
Contour plots of the density and pressure at four different times t = 0.2, 0.4, 0.8, and 1.2ms are presented in Fig. 4, where we have performed the computation using h = 0 and 0.95 with a 400 × 250 grid. From the density plot, we clearly observe the improvement on the use of the unified coordinate method over the Eulerian method to the sharpness of the solution structure near the interfaces. Moreover, from the pressure plot, we see the smooth variation of the solution near the interface, without introducing any spurious oscillations. In Fig. 5, we present the physical grid system for the unified coordinate run shown in Fig. 4 where the dynamical movement of the grids is seen.
References
[1] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach, J. Comput. Phys. 125 (1996), 150–160.
[2] D. Bale, R. J. LeVeque, S. Mitran, and J. A. Rossmanith, A wave propagation method for conservation laws and balance laws with spatially varying flux functions, SIAM J. Sci. Comput. 24 (2002), 955–978.
[3] B. Despres and C. Mazeran, Lagrangian gas dynamics in two dimensions and Lagrangian systems, Arch. Rat. Mech. Anal. 178 (2005), 327–372.
[4] J. Grove and R. Menikoff, The anomalous reflection of a shock wave at a material interface, J. Fluid Mech. 219 (1990), 313–336.
[5] W. H. Hui, The unified coordinate system in computational fluid dynamics, Comm. Comput. Phys. 2 (2007), 577–610.
[6] W. H. Hui and S. Koudriakov, A unified coordinate system for solving the three-dimensional Euler equtions, J. Comput. Phys. 172 (2001), 235–260.
[7] W. H. Hui, P. Y. Li, and Z. W. Li, A unified coordinate system for solving the two-dimensional Euler equtions, J. Comput. Phys. 153 (1999), 596–637.
[8] S. Karni, Multicomponent flow calculations by a consistent primitive algorithm, J. Comput. Phys. 112 (1994), 31–43.
[9] R. J. LeVeque, Finite volume methods for hyperbolic problems, Cambridge University Press, 2002.
[10] R. Menikoff and B. Plohr, The Riemann problem for fluid flow of real materials, Rev. Mod. Phys. 61 (1989), 75–130.
[11] J. A. Rossmanith, A wave propagation method for hyperbolic systems on the sphere, J. Com-put. Phys. 213 (2006), 629–658.
[12] J. A. Rossmanith, D. Bale, and R. J. LeVeque, A wave propagation algorithm for hyperbolic systems on curved manifolds, J. Comput. Phys. 199 (2004), 631–662.
[13] C. W. Schulz-Rinne, J. P. Collins, and H. M. Glaz, Numerical solution of the Riemann problem for two-dimensional gas dynamics, SIAM J. Sci. Comput. 14 (1993), 1394–1414. [14] K.-M. Shyue, An efficient shock-capturing algorithm for compressible multicomponent
a) Density air water h = 0 h = 0.95 t = 0.2ms t = 0.4ms t = 0.8ms t = 1.2ms
Figure 4. Contour plots for the simulation of a model underwater explosion problem. (a) The density and (b) the pressure are shown at four different times t = 0.2, 0.4, 0.8, and 1.2ms obtained using h = 0 and 0.95 with a 400 × 250 grid.
[15] , An efficient unified coordinate method for compressible multicomponent problems, In preparation (2008).
[16] P. D. Thomas and C. K. Lombard, Geometric conservation law and its application to flow computations on moving grids, AIAA J. 17 (1979), 1030–1037.
[17] J. F. Thompson, B. K. Soni, and N. P. Weatherill, Handbook of grid generation, CRC Press, 1999.
[18] M. R. Visbal and D. V. Gaitonde, On the use of higher-order finite-difference schemes on curvilinear and deforming meshes, J. Comput. Phys. 181 (2002), 155–185.
[19] P. Woodward and P. Colella, The piece-wise parabolic method (PPM) for gas dynamical simulations, J. Comput. Phys. 54 (1984), 174–201.
Department of Mathematics, National Taiwan University, Taipei, Taiwan 10617 E-mail address: [email protected]
b) Pressure air water h = 0 h = 0.95 t = 0.2ms t = 0.4ms t = 0.8ms t = 1.2ms Figure 4. (continued) −1 0 1 −1 −0.5 0 0.5 −1 0 1 −1 −0.5 0 0.5 t = 0.2ms t = 0.4ms −1 0 1 −1 −0.5 0 0.5 −1 0 1 −1 −0.5 0 0.5 t = 0.8ms t = 1.2ms
Figure 5. Physical grid system for the unified coordinate run with h = 0.95 shown in Fig. 4. For clarity, a grid coarsening factor of 5 in each x1- and x2-direction is used to make each of the graphs.
