for Compressible Two-Phase Flows
Keh-MingShyue
Abstract.
We present a simple approach to the computation of a simpli ed
two-phase ow model involving gases and liquids separated by interfaces in
multiple space dimensions. In contrast to the many popular techniques which
are mainly concerned with the incompressible ow, we consider a
compress-ible version of the model equations without the e ect of surface tension and
viscosity. We use the
-law and the Tait equation of state for approximating
the material property of the gas and the liquid in a respective manner. The
algorithm uses a volume-of- uid formulation of the equations together with a
sti ened gas equation of state that is derived to give an approximate model
for the mixture of more than one phase of the uids within grid cells. A
stan-dard high-resolution shock capturing method based on a wave-propagation
view-point is employed to solve the proposed model. We show results of some
preliminary calculations that illustrate the viability of the method to practical
application without the occurrence of any spurious oscillation in the pressure
near the interfaces. This includes results of a planar shock wave in water over
a bubble of air.
1. Introduction
We are interested in a two-phase ow problem with interfaces that separate regions
of two di erent uid components consisting of gases and liquids. We consider a
two-dimensional ow as an example, and use the compressible Euler equations to
modeling the motion of the gases,
@ @t 0 B B @ u v E 1 C C A
+
@ @x 0 B B @ u u 2+
p uv(
E+
p)
u 1 C C A+
@ @y 0 B B @ v uv v 2+
p(
E+
p)
v 1 C C A= 0
:(1)
Here
is the density,
uand
vare the velocities in the
x- and
y-direction
respec-tively,
pis the pressure, and
Eis the total energy per unit mass. We assume that
ThisworkwassupportedinpartbyNationalScienceCouncilofRepublicofChinaGrants NSC-86-2115-M-002-005andNSC-87-2115-M-002-016.
the equation of the state for the gas satis es the
-law, where
is the ratio of
speci c heats (
>1). Then the internal energy is
e= 1
1
p ;(2)
and
E=
e+ (
u 2+
v 2)
=
2. We note that the four components of Eqs. (1) express
the conservation of mass, momenta in the
x- and
y-direction, and energy,
respec-tively [5].
As it is often the case, we take the liquids to be baratropic that the pressure
p
is a function of the density
only, and in particular it ful lls the Tait equation
of state of the form,
p
(
) =
AB;
(3)
where
Aand
Bare material-dependent constants; these two parameters along with
give a fundamental characterization of the uid properties of interest and can
be obtained from a tting procedure of laboratory data [26]. Typical set of values
are, for water,
= 7,
A= 3001 atm, and
B= 3000 atm [5], and for human
blood,
= 5
:527,
A= 614
:6 MPa, and
B= 614
:6 MPa [17], approximately. It is
important to note that because of the mere dependency of the pressure
pon the
density
, the energy equation of the Euler Eqs. (1) plays no active role to the
determination of this ow behavior, and hence can be neglected. For completeness,
we write down the equations of motion for the liquid as follows,
@ @t 0 @ u v 1 A
+
@ @x 0 @ u u 2+
p(
)
uv 1 A+
@ @y 0 @ v uv v 2+
p(
)
1 A= 0
;(4)
where the pressure
pis governed by (3).
We want to use a state-of-the-art shock-capturing method on a uniform
rect-angular grid for the computation. For convenience, let us suppose that for each
grid cell we have a volume-fraction function
Zrepresenting the type of uid within
the cell. For example, for the liquid only cells we may take
Z= 1, and therefore
for the gas only cells we may set
Z= 0. In case there are some cells cut by the
interfaces where
Z2(0
;1), we then have both of the gas and liquid components
in the cells in which the liquid and the gas are occupied by the volume fractions
Z
and 1
Z, respectively.
It is clear that, when
Z= 0 or
Z= 1, there is no problem in describing
the motion of each of the gas and liquid ows individually, see Eqs. (1){(4). The
principle problem in the current application and also in the other multicomponent
problems (cf. [1, 4, 7, 14, 15]) is however directed to the development of an ecient
and yet accurate way that is capable of dealing with the uid-mixture case when
Z 2
(0
;1). Motivated by the previous work by the author [20], the algorithm we
employed uses a mixture model based on a fraction (or called a
volume-of- uid) formulation of the equations that is devised to model the mixture of the
two di erent equations of state (2) and (3) by the so-called sti ened gas equation
of state (6) (see Section 2). Then from the energy equation of (1), we choose
conditions that should be satis ed to ensure the pressure in equilibrium for these
cells (see Section 3). Numerical results to be presented in Section 4 show that this
is a viable approach for practical problems without any spurious oscillations in
the pressure near the interfaces when we approximate the model equations in a
consistent manner by a uctuation-and-signal type of shock capturing method [11,
12, 19].
It is true that for real applications the e ect of surface tension as well as
viscosity are two important elements to the solutions of a two-phase ow problem
under studied [18]. Among many of the approaches developed over the years to deal
with these situations (cf. [10, 25]), the method based on the level set formulation [2,
22, 23] has shown to be quite e ective for an incompressible version of the problem
where the uids are mostly in a low Mach number regime. Here in contrast to
the work just mentioned, we consider a class of problems where the in uence of
compressibility of the uids to the solutions is vital, but not the surface tension
and viscosity. Examples of this kind cover a family of shock wave problems with
interfaces [8, 9]. Our goal here is to establish a basic solution strategy for the
problem (cf. [3, 6, 7, 24] for other approaches). Realization of this approach that
couples with adaptive grid procedures such as front tracking and adaptive mesh
re nement is in progress (cf. [21]). In a future work, we will take more physical
factors into account, and also look at methods that are ecient for problems with
low and high Mach number ow where the time step restriction is an important
issue to be addressed.
This paper is organized as follows. In Section 2, we remark and discuss the
equations of state that are used in our model two-phase ow problems with gases
and liquids. In Section 3, we review the derivation of the model system based on
a volume-of- uid formulation of the equations. Numerical results for some sample
problems are shown in Section 4.
2. Equations of State
It is known that in gases the speci c entropy, denoted by
S, has great in uence
to the behavior of the pressure
p, while in liquids the in uence of its changes is
negligible to the pressure [5, 26]. In fact, because of this little dependency of the
pressure on the entropy and a consequence of the rst law of thermodynamics, we
may write the internal energy of a liquid as a separable energy of the form
e
=
e (1)(
) +
e (2)(
S
)
;where with the Tait equation of state (3) we nd
e (1)
(
) = 1
1
p(
) +
B :(5)
From this, it is easy to observe strong resemblance of the equation of state of a
liquid to that of an ideal gas (2). With this in mind, it should be sensible to use
a generalized version of Eq. (5), namely, the sti ened gas equation of state,
e= 1
1
p+
B ;(6)
without the explicitly dependence on the density
to the pressure
pto modeling
mixtures of the components of gases and liquids within a grid cell. Note that
a sti ened gas reduces to a
-law gas when
B= 0, and to a baratropic gas of
the form (3) when the pressure
pis independent of the entropy
S. Oftentimes,
the equation of state (6) has been used in other applications to model materials
including compressible liquids and solids [16]. Here we nd it is suitable for the
current application, when combining (6) with a two-phase ow algorithm described
in the next section, see results shown in Section 4 for numerical veri cation of this
statement.
3. Volume-Of-Fluid Algorithm
We now discuss the basic idea of our approach to modeling cells that contain both
of the gas and liquid phases. We use the Euler Eqs. (1) as a model system that
describes the motion of the mixtures of conserved variables such as
,
u,
v, and
Ein a gas-liquid coexistent grid cell. Then with the use of the equation of state (6)
the pressure
pin the equations (1) can be set by
p= (
1)
E(
u)
2+ (
v)
22
B;(7)
when the mixture of material-dependent parameters
and
Bare de ned and
known a priori. It is important to note that in our approach the conditions for
and
Bare chosen so that the pressure
pin (7) retains in equilibrium for a
gas-liquid mixture cell. The derivation of the evolution equations for
and
Bin the
current case follows directly from procedures developed in [20]. However, to make
the paper self-contained, a brief overview of this approach is undertaken below.
We begin by considering a model two-phase ow problem that the pressure
p
and the velocity eld (
u;v) are constant in the domain, while the other
vari-ables such as the density
and the equation of state parameters:
and
B, are
having jumps across a gas-liquid interface. We write Eqs. (1) in the following
non-conservative form,
@ @t+
u @ @x+
v @ @y+
@u @x+
@v @y= 0
; @u @t+
u @u @x+
v @u @y+ 1
@p @x= 0
; @v @t+
u @v @x+
v @v @y+ 1
@p @y= 0
; @ @t(
e) +
@ @x(
eu) +
@ @y(
ev) +
p @u @x+
@v @y= 0
;and obtain easily equations describing the motion of the gas-liquid interface as
@ @t+
u @ @x+
v @ @y= 0
; @ @t(
e) +
u @ @x(
e) +
v @ @y(
e) = 0
:(8)
To see how the pressure
pwould keep in equilibrium as it should be for this
model problem, we insert the equation of state (6) into (8), and have
@ @t p
+
B1
+
u @ @x p+
B1
+
v @ @y p+
B1
= 0
:(9)
Then with the volume-fraction function
Zde ned previously in Section 1, the
total internal energy
ecan be expressed as a function of the volume fraction by
e
=
p+
B1 =
Z p (l)+
(l) B (l) (l)1
+ (1
Z)
p (g) (g)1
;(10)
where the superscript \
l" and \
g" of the states
p,
, and
Brepresent the values
for the liquid and gas phases, respectively. When substituting (10) to (9), and
after some algebraic manipulation, it is not dicult to show the satisfaction of
the required pressure equilibrium,
p=
p(l)
=
p(g)
, when the following evolution
equation for the volume-fraction function is hold,
@Z @t
+
u @Z @x+
v @Z @y= 0
:(11)
Note that in this instance, the mixture of
and
Bare computed by
= 1 + 1
Z (l)1 +
1
Z (g)1
;(12a)
and
B=
Z (l) B (l) (l)1
1 +
Z (l)1 +
1
Z (g)1
:(12b)
In summary, the model equations we proposed to solve two-phase ow
prob-lems with gases and liquids consist of the following equations,
8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > : @ @t
+
@ @x(
u) +
@ @y(
v) = 0
@ @t(
u) +
@ @x(
u 2+
p) +
@ @y(
uv) = 0
@ @t(
v) +
@ @x(
uv) +
@ @y(
v 2+
p) = 0
@ @t(
E) +
@ @x[(
E+
p)
u] +
@ @y[(
E+
p)
v] = 0
if 0
Z<1
@Z @t+
u @Z @x+
v @Z @y= 0
;(13)
where in case
Z= 1 the Tait equation of state (3) is used to compute the pressure
pfor a liquid, while Eq. (7) is employed when 0
Z <1 with
and
Bde ned
by (12a) and (12b), respectively. We note that the model Eqs. (13) is not written in
the full conservation form, but is rather a quasi-conservative system of equations.
Numerical approximation based on the discretization of the model equations via
the wave propagation approach [11, 12] has shown to be quite robust for a wide
variety of problems of practical interests, see results presented in Section 4. The
algorithmic details as well as many other numerical aspects of this approach may
be found in the papers [20, 21].
4. Numerical Results
We now present two sample calculations to demonstrate the feasibility of the
pro-posed approach described in Section 3 for practical applications (cf. [21] for more
numerical results).
Example 4.1.
To begin, we consider a radially symmetric problem that
the computed solutions can be compared with the one-dimensional results for
numerical validation. We use the same setup as done by Davis [6] that inside
a circular membrane of radius
r0
= 0
:
2 the uid is air with density
= 1
:25,
pressure
p= 2
:75, and the adiabatic gas constant
= 1
:4, while outside the
circular membrane the uid is liquid with density
= 1 and the Tait equation of
state (3) of parameter values:
A= 1,
B= 1, and
= 7. Initially both of the air and
the liquid are in a stationary position, but due to the pressure di erence between
the uids, breaking of the membrane occurs instantaneously. For this problem, the
resulting solution consists of an outward-going shock wave in liquid, an
inward-going rarefaction wave in air, and a contact discontinuity lying in between that
separates the air and the liquid.
In Fig. 1, we show numerical results for the density
as well as the pressure
pat time
t= 0
:07. We use a 100
100 grid with the high resolution version of the
wave propagation method (cf. [11, 12]) for the computation (the computational
domain is a square with dimensions [0
;0
:5]
[0
;0
:5]). From the contours of the
plots, it is easy to see the wave pattern as just mentioned after the breaking of the
membrane, where the dashed line in the pressure contours is the approximate
loca-tion of the interface. From the scatter plots, we nd good agreement of the results
as compared with the \true" solution obtained from solving the one-dimensional
multicomponent model with appropriate source terms for the radial symmetry
using the high-resolution method [20] (with mesh side
h= 1
=2000). Noticeably,
here the pressure near the interface behaves in a satisfactory manner without any
spurious oscillations, and the shock wave and contact discontinuity appear to be
very well located.
Example 4.2.
Our next example concerns a planar shock wave in water
over a bubble of air. We take an initial condition that is similar to the one used
by Grove and Meniko [9] where anomalous re ection of a shock wave at a uid
Density
a)
Pressure
0.0 0.2 0.4 0.6 0.8 0.9 1.0 1.1 1.2r
.... . . . . . . . . . . . . . ....... ... . . . . ... . . ... ... . . . . . . . . . . . . . . ....... ... . . . . ... . . ... ... . . . . . . . . . . . . . . ....... ... . . . . ... . . ... ... . . . . . . . . . . . . . . ....... . . . . ... . . ... ..... . . . . . . . . . . . . . . ....... . . . . ... . . ... ... . . . . . . . . . . . . . . .. ....... . . . . ... . . ... .... . . . . . . . . . . . . . . ...... ... . . . . . ... . . ... ... . . . . . . . . . . . . . . . ........ . . . . ... . . ... ..... . . . . . . . . . . . . . . . ....... . . . . .... . . ... ..... . . . . . . . . . . . . . . . ....... ... . . . . ... . . ... .... . . . . . . . . . . . . . . .. ....... . . . . ... . . ... .... . . . . . . . . . . . . . . . ...... ... . . . . . ... . . . ... ..... . . . . . . . . . . . . . . . .. ....... . . . . ... . . .... ..... . . . . . . . . . . . . . . . . ... ... . . . . ... . . ... .... ... . . . . . . . . . . . . . . . ...... ... . . . . ... . . ... .... .. . . . . . . . . . . . . . . . .. ....... . . . ... ... . . . ... ... ... . . . . . . . . . . . . . . .. ... ... . . . . ... . . ... ..... . . . . . . . . . . . . . . .. ....... . . . . ... . . . ... .... . . . . . . . . . . . . . . .. ....... . . . .... . ... ..... . . . . . . . . . . . . .. ... ... . . . . ... . . ... .... . . . . . . . . . . . .. ...... ... . . . . . .... . . ... .... ... . . . . . . .. ....... ... . . . . .... . . ... .... ... .. .. .. ... ....... . . . ... . ... ...... .. .... ...... ... . . . . ... . . ... ... ... ... ...... ... . . . . ... . . ... ... ... ....... . . . . ... . . .... ... ....... . . . . ... . . ... ....... . . . . .... . . ... ....... . . . . ... . . .... ... . . . . . ... . . ... ... . . . . ... . . ... ... . . . . . ... . . ... ... . . . . . ... . . ... ... . . . . ... ... . . . ... ... . . . . . ... . . ... ... . . . . . ... . . ... ... . . . . . ... .... . . ... ... . . . . . ... ... . . . ... ... . . . . . . . ... . . ... ... . . . . . . . . ... . . ... ... .. . . . . . . ... . . ... ... .. . . . . . . . ... . . ... ... .. . . . . . . .. ... . . ... ... .. . . . . . . . . . ... . . ... ... .. . . . . . . . . . . ... . . ... ... .. . . . . . . . . . . .. .... . .... ... .. .. . . . . . . . .. . ... . . ... .... .. .. . . .. . .. .. ... . . ... ... .. .. . . . .... ... . . ... ... .. ... ... . . . ... ... .... . . ... ... . . . ... ... . . . ... .... . . ... .... . . ... ... . . .... ... . . . ... .... . . ... ... . . ... .... . . ... .... . . ... ... . . ... ... . . ... .... . . ... .... . . ... .... . . ... ... . . . ... ... . . . ... ... . . . ... ... . . . ... .... . . . ... ... . . . ... ... . . . ... ... . . . . ... ... . . . ... ... . . . . ... .... . . . ... .... . . . ... ... . . . . ... ... . . . . ... ... . . . . ... ... . . . . ... .... . . . ..... ..... . . . . ... ..... . . . . ..... ..... . . . . ..... ..... .. . . . . ... ..... .. . .. .. ... ... ... .... ....... ... ....................................Density
b)
0.0 0.2 0.4 0.6 0.0 0.5 1.0 1.5 2.0 2.5r
..... . . . . . . . . . . .. ...... ...... . . . . ... ..... . . . . . . . . . . .. ...... ... . . . . ... .... . . . . . . . . . . ...... ...... .. . . . .... ... . . . . . . . . . . .. ...... ...... . . . .... ... . . . . . . . . . .. ...... ...... . . . ... ..... . . . . . . . . . . .. ...... ...... . . . ... .... . . . . . . . . . . ... ...... ... . . . ... ..... . . . . . . . . . . .. ...... ...... . . . . ... ..... . . . . . . . . . . .. ...... ...... . . . . ... ..... . . . . . . . . . . .. ...... ...... . . . . ... ..... . . . . . . . . . . . .. ...... ...... . . . . ... ..... . . . . . . . . . . . .. ...... ...... . . . . ... ..... .. . . . . . . . . . . .. ...... ...... . . . . ... ..... ... . . . . . . . . . .. ...... ...... . . . .... ..... . . . . . . . . . . . .. ...... ...... . . . ... ..... ... . . . . . . .. .. ...... ...... . . . . ... ..... . . . . . . . . .. .. ... ...... ... . . . ... ..... ... . .. . .. .. ...... ...... . . . . ... ..... .. .. .. .... ...... ...... . . . .... ..... .... .. ...... ...... ..... . . . ... ... ... .... ... ... ...... ... . . . ... ...... ... ...... ...... ..... . . . ... ...... .... ...... ...... . . . .... ...... ...... ...... . . . ... ...... ...... .... . . . ... ...... ...... . . . . ... ...... ...... . . . ... ...... ...... . . . ... ...... ... . . . . ... ...... ... . . . ... ...... ... . . . ... ...... .... . . . ... ...... ... . . . ... ...... .... . . . ... ...... ... . . . ... ...... ... . . . ... ...... ... . . . ... ... ...... . . . ... ...... .... . . .... ...... ... . . . ... ...... ... . . . .... ... ....... . . . ... ... ... . . . ... ...... ... . . . ... ...... .. . . . ... ...... ... . . . ... ...... .. . . . ... ...... . . . ... ...... .. . . . ... ...... .. . . . ... ...... . . . ... ...... . . . . ... ...... . . . . ... ...... . . . ... ...... . . . .... ...... . . . . ... ...... . . . . ... ...... . . . ... ...... . . . . ... ...... . . . ... ...... . . . ... ...... . . . .... ...... . . . . ... ...... . . . . ... ...... . . . . ... ...... . . . ... ...... . . . ... ... . . . . ... ... . . . . ... ..... . . . . ... ..... . . . . ... ... . . . . ... ... . . . . ... ... . . . . . ... ... . . . . ... ... . . . . . ... ... . . . . . ... ... . . . . . ... ... . . . . . ... ... . . . . . ... ..... . . . . . ... ... . . . . . . ... ..... . . . . . . ... .... . . . . . . ..... ..... . . . . . . . ... ..... . . . . . . . . ... ... ... . . . . . . . . . ... ..... .. . . . . . . . .. ... ..... . . . . . . .. ....... ... ... ......... ...Pressure
Figure 1.
High resolution results for a radially symmetric
prob-lem at time
t= 0
:07. (a) Contours of the density
and the
pres-sure
p. (b) Scatter plots of the density
and the pressure
pwith
locations measured as a distance from a point of the solution to
the center (0
;0). The solid line in the scatter plot is the \true"
solution obtained from solving the one-dimensional
multicompo-nent model with appropriate source terms for the radial symmetry
using the high-resolution method. The dotted points are the
two-dimensional result. The dashed line in the pressure contour is the
approximate location of the interface.
interface is examined closely there. Here to test the proposed method, we consider
a downward-moving Mach 1
:587 shock wave in water with data in the
and data in the post-shock state by
(
; u; v; p)post-shock = (1
:233
;0
;43
:467
;10
4)
:
The parameters we used for the equation of state of the water are:
A= 3001 atm,
B= 3000 atm, and
= 7. In addition to this, we assume that there is a stationary
bubble of air of radius 0
:2 located just below the shock with density
= 0
:0012
(di-mensional unit g/cc) and the adiabatic gas constant
= 1
:4 for the
-law equation
of state (2). We note that the ratio of acoustic impedances (
c)
water=
