Mach-Uniform Methods for Compressible Two-Phase Flow
Problems
Keh-Ming Shyue
Department of Mathematics, National Taiwan University, Taipei 106, Taiwan E-mail address: [email protected]
1
Introduction
Our aim in this project is to develop a Mach-uniform method for nearly compressible two-phase flow problems where the flow speed is much less than the speed of sound of the underlying fluid components. The accuracy and efficiency of the method under development should depend weakly on the Mach number in the full range from zero to supersonic values. Representative applications of this kind are such as the rising of gas bubbles in liquids, the falling of liquid drops in the air under gravitational force field, and bubbly flow in liquids, see [12, 15] for more examples. To get a quick impression of some of the difficulties that may occur when one attempts to extend a state-of-the-art compressible solver to a low speed flow scenario, we begin by considering the one-dimensional Euler equations of gas dynamics of the form
∂q
∂t +
∂f (q)
∂x = 0, (1)
where the vector for the conservative variables q and the flux vector f are defined by
q = ρ ρu ρE and f (q) = ρu ρu2+ p ρEu + pu .
Here ρ, u, p, and E, denote the density, the particle velocity, the pressure, and the specific total energy, respectively. For simplicity, we assume that the thermodynamical behavior of the flow satisfies a γ-law gas as
p = (γ − 1)ρe, (2)
where e is the specific internal energy, γ is the ratio of specific heats (1 < γ ≤ 5/3), and as usual E = e + ρu2
/2.
Now by introducing the following set of the dimensionless variables into the above equations: ˜
ρ = ρ/ρ0, u = u/u˜ 0, p = p/(ρ˜ 0c20), E = E/c˜ 2
0, ˜x = x/x0, ˜t = (u0t)/x0,
we find easily a dimensionless form of (1) and (2) as ∂ ∂˜t ˜ ρ ˜ ρ˜u ˜ ρ ˜E + ∂ ∂ ˜x ˜ ρ˜u ˜ ρ˜u2 ˜ ρ ˜E ˜u + ˜p˜u + 0 1 M2 0 ∂ ˜p ∂ ˜x 0 = 0, (3) and ˜ p = (γ − 1) ˜ ρ ˜E −M 2 0 2 ρ˜˜u 2 , (4)
respectively, where M0= u0/c0is a dimensionless Mach number, and the quantity κ0represents a reference state for the variable κ, for κ = ρ, u, and so on. For M0→ 0 (incompressible flow limit), it is easy to see that (3) develops a singularity, because ∂x˜p/M˜ 0→ ∞. Moreover, it can be shown that the eigenvalues of the Jacobian matrix of the flux function in (3) are ˜u and ˜u ± ˜c/M0, where c2= γp/ρ. Clearly, the acoustic wave speeds ˜u ± ˜c/M
0 degenerate as M0→ 0.
Surely, one possible way to simulate low speed flow problems is to treat such a flow as a fully compressible flow and use an explicit method. It is known that with the use of an explicit method we would have a time step restriction, the CFL condition,
νa = ∆t maxj∆x(|uj| + cj) ≤ 1, (5)
which states that, with a given mesh spacing ∆x, for stability the time step size must be inversely proportional to the maximum of the sum of flow speed and the sound speed. Because of this, explicit methods are well-suited only for problems where the flow speed is on the same order of or larger than the sound speed. In the current interests of low speed problems, however, the sound speed could be orders of magnitude larger than the flow speed, thus grossly over resolving in time features of the fluid flow, see [10, 11] for more details.
The central idea of many numerical solver for low-speed flow is then somehow to separate out the “ill-conditioned” part of the flow and treat it implicitly, while hoping that the rest of the flow may be advanced at an acceptable time step with an economical and accurate method. Such is the approach first taken by Harlow and Amsden [5] and later expanded by many others [2, 8, 18] for examples. A common theme in all these method is splitting. That is, we separate the flow variable into various parts, and identify the part of the flow that needs implicit treatment. When this is done properly, we may bypass the stringent CFL condition dominated by the sound speed as in (5), and use a more robust flow speed-based CFL condition as
νf =∆t max∆xj|uj| ≤ 1. (6)
It should be remarked that a somewhat different approach is proposed by Colella and Pao [3] where by using the Hodge decomposition, they write the velocity field as the sum of a divergence-free vector field and a gradient of a scalar function. On the divergence-divergence-free part, the solution of the equations would vary on a time scale determined by the flow speed, while on the other part, the solution of the equations would depend on fast sound waves. Sensible results were reported by their method for some sample problems. Besides this projection method, Klein and co-workers have also proposed a numerical method based on asymptotic [7, 9], which were an extension of the low Mach number asymptotic of Klainermann and Majda [6]. As to the other possible approaches to low speed flow in the literature, see [1, 4, 16] for an example.
It is important to note that almost all the low speed flow methods mentioned above is to single component flow problems. The objective of the current work is devoted to the generalization of many of these methods to problems with more than one fluid component, and with real fluids characterized by non-convex equations of state.
2
Pressure-Correction Methods
2.1
Conservative formulation
We begin by generalizing the barely implicit correction method of Patnaik et al. [8] originally devised for single phase flows to two phases. To do this, analogously to [8], there are two basic steps in the algorithm. In the first step, we solve the continuity and momentum equations of (1) explicitly as ρ∗− ρn ∆t = −Dx(ρu) n, (7) (ρu)∗− (ρu)n ∆t = −Dx ρu 2 + pn , (8)
over a time step ∆t (which is chosen based on the local flow speed νf, but is not based on the local acoustic wave speed νa) for an intermediate states of ρ∗ and (ρu)∗ using the initial data at time tn. Here Dxrepresents a finite-difference operator to the differential operator ∂/∂x, and for simplicity we have omitted the subscript for discretizing the spatial variables. Then, in the second step, we consider an implicit finite-difference scheme of the form for the momentum and energy equations as (ρu)n+1− (ρu)n ∆t = −Dx n ρu2n +ωpn+1+ (1 − ω)pno , (9) (ρE)n+1− (ρE)n ∆t = −Dx(ρE + p) n ωun+1+ (1 − ω)un , (10)
where ω is an implicitness parameter, typically 0.5 ≤ ω ≤ 1 for a stable calculations. Note that we have only treated u and p as implicit variables, but leave ρ as an explicit one.
If we now define the change in pressure as
φ = ω pn+1− pn , (11)
then the correction equation for momentum can be obtained in terms of φ by subtracting (8) from (9), yielding (ρu)n+1− (ρu)∗ ∆t = −Dxω p n+1− pn = −D xφ, (12) and so un+1= u∗−∆t ρ∗Dxφ, (13)
after rearranging (12) and letting ρn+1= ρ∗. Now, with (13), we may derive a single equation for φ by eliminating un+1from (9) and (10).
To accomplish this, we first use the equation of state
(ρe)n+1= pn+1 γn+1− 1 = γn− 1 γn+1− 1 (ρe)n+ φ ω(γn+1− 1) (14)
to obtain a correction equation for the internal energy. Note that for single phase flow problems the ratio of specific heats γ is a fixed constant at all time, but this is no longer true for a multiphase application concerned here. For this reason, we have to come up with a way to compute γn+1at all space and time locations; this can be done in the same way as in [13]. For the moment, suppose that we have already gotten the value, the procedure for finding φ proceeded by substituting (13) and (14) into (10), leading to
1 ∆t " γn− γn+1 γn+1− 1 (ρe)n+ φ ω(γn+1− 1)+ (∆tDxφ)2 2ρ∗ + u ∗∆tD xφ + ρu2∗ 2 − ρu2n 2 # = ω∆tDx (ρE + p)n ρ∗ Dxφ
− Dx{(ρE + p)n[ωu∗+ (1 − ω)un]} .
For convenience, we define the quantity (ρE)∗ as (ρE)∗− (ρE)n
∆t ≡ −Dx{(ρE + p)
n[ωu∗+ (1 − ω)un]} , (16)
and this allows us to rewrite (15) in the form φ ω(γn+1− 1)∆t+ ∆t(Dxφ)2 2ρ∗ + u ∗D xφ − ω∆tDx (ρE + p)n ρ∗ Dxφ = (ρE)∗− (ρE)n ∆t − ρu2∗ − ρu2n 2∆t − 1 ∆t γn− γn+1 γn+1− 1 (ρe)n, (17)
which provides us with a nonlinear elliptic equation for φ. Note that in the original work of [2, 8], the second and third terms from the left-hand side of (17) are not considered in the method. Surely, the effect of these terms to the solution of the method will be one of the major things to look at. No matter what form of the elliptic equation is being used, its solution can be employed to correct un+1 and (ρe)n+1by (12) and (14), respectively.
Up to this point, it should be clear that as compared to the traditional compressible flow solver, the complexity of this algorithm is increased only by an additional solve for the elliptic equation (17), where the right-hand size of the equation as well as the coefficients of the equation can be obtained basically by a predictor solve of (1) and the initial condition of the problem. The proper discretization of (17) as well as the performance of the method to achieve the Mach uniformity and the divergence-free constraint as the Mach number is approaching to zero will be studied and tested carefully in this project.
2.2
Non-conservative formulation
Note that if we start with the basic thermodynamic relation
∆p = ∂p ∂ρ e ∆ρ + ∂p ∂e ρ ∆e and with the perfect gas equations of state, we would have
∆p = p
ρ∆ρ + (γ − 1)ρ∆e, (18)
where ∆p = pn+1− p∗. Suppose that ∆ρ and ∆e can be well-approximated by
∆ρ = −ρ∗∆tD xun+1 and ∆e = −p ρ∗∆tDxu n+1,
respectively. We may then rewrite (18) as
∆p = γp∆tDxun+1= ρc2∆tDxun+1, where ρc2
is assumed to be evaluated at some prescribed state, say at the intermediate state ∗, but surely it can treated in an implicit manner also. Having that, the velocity un+1in above can be eliminated by using the equation of motion:
∆u = un+1− u∗= −Dxp n+1
and thus we obtain the equation for pn+1 Dx Dxp∗ ρn = p∗− p n (ρc2)n(∆t)2 + Dxun ∆t . (19)
The approach just described is basically the method proposed by Yabe and Wang [17] for a universal treatment of compressible and incompressible flows. The applicability of this method to general multicomponet problem has yet to be investigated, though some promising results has been shown for a wide variety of problems [18].
It should be mentioned that, in the method of Yabe and Wang, the Euler equations are rewritten in the following non-conservative form,
∂q ∂t + u ∂q ∂x = ψ, (20) where q = ρ u p and ψ = −ρux −px −γpux .
In each time step, the method consists of the following steps: (1) Non-advection step
With the initial condition at the current time and a time step ∆t chosen by (6), solve the following set of equations,
∂q
∂t = ψ, (21)
to obtain the intermediate states ρ∗, u∗, and p∗. In practice, this is done by first solving the Poisson equation (19) for the intermediate state of the pressure p∗, and then it is used in the discretization of the equations for the density and velocity, yielding ρ∗ and u∗.
(2) Advection step
After ρ∗, u∗, and p∗have been gotten in the non-advection step, we apply a numerical method to update the solutions of the transport equations in (21),
∂q ∂t + u
∂q
∂x = 0,
at the next time step ρn+1, un+1, and pn+1.
3
Numerical Examples
We now present sample one-dimensional results obtained using a preliminary version of the Mach-uniform method based on the non-conservative formulation of Yabe and Wang, see [14] for more numerical examples.
Example 3.1. We begin by considering an interface only problem that the solution of a
Riemann problem consists of a single contact discontinuity in gas dynamics. The initial condition we use is ρ u p γ L = 1 1 1 1.4 and ρ u p γ R = 0.125 1 1 1.2 ,
Here L is the state used for x ∈ [0, 0.5) and R is the state used for x ∈ [0.5, 1].
Numerical results obtained using the uniform method as well as the classical non Mach-uniform method are shown in Fig. 1 at time t = 0.2 with a 100 grid. From the graphs, it is clear that in each of these methods there is no spurious oscillation in the pressure and also the particle velocity near the diffused material interface. We have gotten a slightly sharp resolution of the interface when the Mach-uniform method is in use as it is compared with the non Mach-uniform method case. Note that in the former case, we have chosen a time step based on the flow speed CFL condition (6), while in the latter case, we have used a time step based on the acoustic wave speed CFL condition (5).
Example 3.2. We are next concerned with a two-phase version of the Sod’s Riemann problem.
On the left when x ∈ [0, 0.5), we have the data
(ρ, u, p, γ)L= (1, 0, 1, 1.4), and on the right when x ∈ [0.5, 1], we have the data
(ρ, u, p, γ)R= (0.125, 0, 0.1, 1.2).
We run the problem in a shock tube with 100 mesh points, and show the results in Fig. 2 at time t = 0.14 using both the Mach-uniform and non Mach-uniform methods as in the previous example. From the figure, we observe some spurious oscillations near the shock wave for the Mach-uniform computation, while there is not any oscillation near the shock wave for the non Mach-uniform computation. In many non-conservative numerical methods alike, this unexpected oscillation may be eliminated by adding numerical artificial viscosity in the method formulation.
Example 3.3. Finally, we consider a shock-contact interaction problem that investigates con-vergence of the computed solutions to the correct weak ones in a multicomponent case. The initial condition we use consists of a stationary interface at x = 0.3 separating two fluids of different equation of states, and a rightward going Mach 1.22 shock wave at x = 0.2 traveling from the left to right. The state on the right of the interface is a helium with
(ρ, u, p, γ)R= (0.138, 0, 1, 1.67),
and the state on the right of the interface, (i.e., on the middle and the preshock state), is an air with
(ρ, u, p, γ)M = (1, 0, 1, 1.4). The state behind the shock is
(ρ, u, p, γ)L= (1.3765, 0.3948, 1.57, 1.4).
Figure 3 shows numerical results of this problem at time t = 0.2. In view of the global structure of the solution is concerned, it is clear that the Mach-uniform method gives the correct solution behavior, but the result is too diffusive as compared to the result obtained using the non Mach-uniform method. It is somewhat surprising that (at least in my perspective) the larger the CFL number is, the more numerical diffusion is in the computed solution (not shown here).
4
Conclusion
The work present here is only the preliminary one. Extension of the basic methodology of the non-conservation formulation to more than one space dimension is straightforward. However, the convergence behavior of the non-conservative method to the correct weak solution should be looked into carefully. In this regards, a conservative formulation should be more favorable for shock-wave problems.
0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.5 1 1.5 2 0 0.5 1 0 0.5 1 1.5 2 0 0.5 1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.5 1 1.5 2 0 0.5 1 0 0.5 1 1.5 2 0 0.5 1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x x x x
Density Velocity Pressure Mach number
a)
b)
Figure 1: Numerical results for an interface only problem at time t = 0.2. a) Mach-uniform computation with a time step chosen by the flow-based CFL condition (6), and b) Classical non Mach-uniform computation with a time step chosen by the acoustic-based CFL condition (5). The solid line is the exact solution and the points shows the computed solution with 100 mesh points.
0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 x x x x
Density Velocity Pressure Mach number
a)
b)
Figure 2: Numerical results for a two-phase Riemann problem at time t = 0.14. The graphs of the solutions are displayed in the same manner as in Fig. 1.
0 0.5 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.5 1 1 1.2 1.4 1.6 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.5 1 1 1.2 1.4 1.6 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 x x x x
Density Velocity Pressure Mach number
a)
b)
Figure 3: Numerical results for a shock-contact interaction problem at time t = 0.2. The graphs of the solutions are displayed in the same manner as in Fig. 1.
References
[1] H. Bijl and P. Wesseling. A unified method for computing incompressible and compressible flows in boundary-fitted coordinates. J. Comput. Phys., 141:153–173, 1998.
[2] V. Casulli and D. Greenspan. Pressure method for the numerical solution of transient, com-pressible fluid flows. Int. J. Numer. Meth. Fluids, 4:1001–1012, 1984.
[3] P. Colella and K. Pao. A projection method for low speed flows. J. Comput. Phys., 149:245– 269, 1999.
[4] H. Guillard and C. Viozat. On the behavior of upwind schemes in the low Mach number limit. Computers and Fluids, 28:63–86, 1999.
[5] F. H. Harlow and A. A. Amsden. Numerical calculation of almost incompressible flow. J. Comput. Phys., 3:80–93, 1968.
[6] S. Klainermann and A. Majda. Compressible and incompressible fluids. Comm. Pure Appl. Math., 35:629–651, 1982.
[7] R. Klein. Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics. J. Comput. Phys., 121:213–237, 1995.
[8] G. Patnaik, R. H. Guirguis, J. P. Boris, and E. S. Oran. A barely implicit correction for flux-correct transport. J. Comput. Phys., 71:1–20, 1987.
[9] T. Schneider, N. Botta, K. J. Geratz, and R. Klein. Extension of finite volume compressible flow solvers to multi-dimensional, variable density zero Mach number flows. J. Comput. Phys., 155:248–286, 1999.
[10] J. Sesterhenn, B. Mueller, and H. Thomann. Flux-vector splitting for compressible low Mach number flow. Computers and Fluids, 22:441–451, 1993.
[11] J. Sesterhenn, B. Mueller, and H. Thomann. On the cancellation problem in calculating compressible low Mach number flows. J. Comput. Phys., 151:597–615, 1999.
[12] B. R. Shin, S. Yamamoto, and X. Yuan. Application of preconditioning method to gas-liquid two-phase flow computations. J. Fluids Eng., 126:605–612, 2004.
[13] K.-M. Shyue. An efficient shock-capturing algorithm for compressible multicomponent prob-lems. J. Comput. Phys., 142:208–242, 1998.
[14] K.-M. Shyue. A pressure-correction method for low speed compressible two-phase flow. In preparation, 2006.
[15] D. R. van der Heul, C. Vuik, and P. Wesseling. A staggered scheme for hyperbolic conservation laws applied to unsteady sheet cavitation. Comput. Visual Sci., 2:63–68, 1999.
[16] J. M. Weiss and W. A. Smith. Preconditioning applied to variable and constant density flows. AIAA J., 33:2050–2057, 1995.
[17] T. Yabe and P.-Y. Wang. Unified numerical procedure for compressible and incompressible fluid. J. Physical Soc. Japan, 60:2105–2108, 1991.
[18] T. Yabe, F. Xiao, and T. Utsumi. The constrained interpolation profile method for multiphase analysis. J. Comput. Phys., 169:556–593, 2001.
