JOURNAL OF COMPUTATIONAL PHYSICS**142, 208–242 (1998)**
ARTICLE NO. CP985930

## An Efficient Shock-Capturing Algorithm for Compressible Multicomponent Problems

Keh-Ming Shyue

*Department of Mathematics, National Taiwan University, Taipei, Taiwan, Republic of China*
E-mail: shyue@math.ntu.edu.tw

Received January 6, 1997; revised January 22, 1998

A simple shock-capturing approach to multicomponent flow problems is devel- oped for the compressible Euler equations with a stiffened gas equation of state in multiple space dimensions. The algorithm uses a quasi-conservative formulation of the equations that is derived to ensure the correct fluid mixing when approximating the equations numerically with interfaces. Aγ -based model and a volume-fraction model have been described, and both of them are solved using the standard high-resolution wave propagation method for general hyperbolic systems of partial differential equa- tions. Several calculations are presented with a Roe approximate Riemann solver that show accurate results obtained using the method without any spurious oscillations in the pressure near the interfaces. Convergence of the computed solutions to the correct weak ones has been verified for a two-dimensional Richtmyer–Meshkov un- stable interface problem where we have performed a mesh-refinement study and also shown front-tracking results for comparison. ° 1998 Academic Pressc

**1. INTRODUCTION**

Our goal is to present a simple approach to multicomponent flow of general compressible materials in more than one dimension. We use the Euler equations of gas dynamics as a model system, and consider problems with the so-called “stiffened” gas equation of state for approximating materials including compressible liquids and solids [29, 50]. The algorithm uses a quasi-conservative formulation of the equations proposed by Abgrall [1] to ensure a consistent approximation of the energy equation near the interfaces where regions of different fluid components are separated. Aγ -based model is therefore derived that extends the work of Abgrall [1] from polytropic gases in one dimension to a stiffened gas, and also to multiple dimensions. We give a new formulation of the resulting model in expression of the volume fraction that is more robust for two-component flow problems. This will be discussed further in Sections 2 and 5.

208 0021-9991/98 $25.00

Copyright c° 1998 by Academic Press All rights of reproduction in any form reserved.

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 209

We use the high-resolution wave propagation method developed by LeVeque [40, 41, 43]

*to solve the proposed multicomponent models. This method is a variant of the fluctuation-*
*and-signal method of Roe [59, 60], and has been widely used in many applications including*
the single-component fluid of ideal gases. The main idea behind the method has recently
been implemented in the software package CLAWPACK (Conservation LAWs PACKage)
as the underlying integration routine [46, 47]. The current use of the method is just an easy
extension of the previous one from single-component to multicomponent problems. It is an
efficient and yet accurate scheme without any spurious oscillations in the pressure near an
interface as illustrated by numerical results presented in this paper.

We will only briefly review and describe the method in a shock-capturing framework
in one dimension, see Section 3. Extensions of the method to front-tracking and to two
dimensions are straightforward also by following the procedures outlined in [44, 45],
for example. We will not discuss a front-tracking version of the method here, but fo-
*cus on the more fundamental shock-capturing algorithm and validate its use via numeri-*
cal experimentations. Some preliminary results obtained using the front tracking method
may be found in [65, 66] for two-dimensional unstable fluid interface problems such as
Rayleigh–Taylor and Richtmyer–Meshkov instabilities, see Section 6 for an example also.

Generalization of the approach to three dimensions can be made in a similar manner, but requires more programming work, especially in regard to a front-tracking method [25].

Clearly for real applications the use of a stiffened gas equation of state that appears in an analytical formula (i.e., the constitutive relation (2) in Section 2) represents only a limited number of materials of practical importance [29, 50]. However, there are some problems of sufficient interest and difficulty that the development of a multicomponent algorithm for this equation of state is worthwhile, particularly since in some cases it is relatively easy to compute the exact solutions and check accuracy of the method.

Numerous numerical methods have been developed over the years to handle multicom- ponent flow problems. Consider a non-reacting ideal gas flow, for example. One popular approach among them is to solve an extended system of equations in which additional conservation equations are introduced to the original Euler equations to describe the con- servation of each fluid component separately. Methods of this type, in particular a shock- capturing version of the method, often fail to maintain pressure equilibrium for grid cells near interfaces where two or more fluid components are mixed. Some representative exam- ples that exhibit this erroneous phenomenon are given in [34] for the use of aγ -based (see Fig. 1 also) and a level-set model, and in [9, 12, 71] for the use of a mass-fraction model.

*The exception is the method explained by Jenny et al. [33] that the fluxes obtained using*
conservative Euler solvers are modified in a suitable way to avoid the occurrence of the
pressure errors generated near the interfaces. It is unclear how to extend the method to a
situation other than the ideal gas flows, however.

Another approach introduced by Karni [35, 36] is to solve the Euler equations separately
on each side of the interface using a method designed for a single-component flow, while
the interface is dealt with in a different manner using a pressure evolution equation derived
from the energy equation. Despite the fact that the method is not exactly conservative at
the interface, reasonable results are obtained using this approach in conjunction with either
standard level-set or mass-fraction formulation of ideal gases. Extension of the method to
*a thermally perfect gas was done recently by Fedkiw et al. [21] in one dimension. It should*
*be noted that Cocchi et al. [12] devised a rather similar method of this kind that employs a*

210 KEH-MING SHYUE

linear interpolation technique near the interface instead; some one-dimensional results are shown for a stiffened gas flow.

Our method to modeling multicomponent flow of general compressible materials is mo- tivated by the work of Abgrall [1] in that, based on the physical principles and from the energy equation, we derive the effective equations (see Section 2) for the mixture of material- dependent quantities near the interface. In this method, with the stiffened gas equation of state, we take these equations to be of the form that do not vary their solutions across the shock and rarefaction waves as well. Combining the resulting effective equations (i.e., Eqs. (9)) to the Euler equations yields a model system that is not written in the full con- servation form, but is rather a quasi-conservative system of equations. Abgrall [1] solves a system of this kind using a predictor-corrector method, while we use the high-resolution wave propagation method that gives an efficient implementation of the algorithm. In prin- ciple, when properly modified, it is possible to employ the state-of-the-art shock-capturing methods for hyperbolic systems of conservation laws for solving the model equations as well; see [67] for an example that generalizes the MUSCL scheme [13, 72] to this specific application and also to a van der Waals gas.

There are many other multicomponent approaches available in the literature. Some typical ones are the level-set methods [19, 52], front-tracking methods [25, 28], volume-of-fluid methods [14, 51, 55, 71], and the BGK-based method [74].

This paper is organized as follows. In Section 2, we begin by discussing the basic com- putational models in one dimension that govern the motion of materials characterized by a stiffened gas equation of state in a multicomponent environment. We will give two different formulations of the model equations that are both applicable for practical computations of multicomponent problems. In Section 3, we briefly review an approximate Riemann solver of Roe, and analyze a first order numerical methods based on a wave-propagation approach with an application to our multicomponent models proposed in Section 2. The results of some one-dimensional tests are given in Section 4 that validate this approach. The basic idea of the algorithm is then extended to multiple dimensions in Section 5, and some two-dimensional numerical results are presented. Section 6 shows results for a Richtmyer–

Meshkov unstable interface problem where a mesh-refinement study is performed to check convergence of the computed solutions.

**2. DERIVATION OF MODEL EQUATIONS**

In one dimension, the single-component Euler equations of gas dynamics take the form

∂

*∂t*

ρ
*ρu*
*ρE*

+ ∂*∂x*

*ρu*

*ρu*^{2}*+ p*
*(ρE + p)u*

= 0, (1)

where*ρ is the density, u is the velocity, p is the pressure, and E is the total energy per unit*
mass. We assume a compressible material that the internal energy per unit mass, denoted
*by e, satisfies the “stiffened” gas equation of state,*

*ρe =* *p+ γ p*_{∞}

γ − 1 , (2)

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 211

*and E= e + u*^{2}*/2. Here γ is the usual ratio of specific heats (γ > 1), and p*∞ is a pre-
scribed pressure-like constant; these values can be used to describe the material property of
*interests and can be determined from laboratory experiments via an empirical fit [31, 49].*

For example, for water we haveγ ∼*= 5.5, p*∞= 4921.15 bars [12], and for tungsten we
getγ ∼*= 3.14, p*∞= 1.0 Mbar [54]. Note a stiffened gas reduces to a polytropic gas when
*p*_{∞}= 0. The three components of Eqs. (1) express the conservation of mass, momentum,
and energy, respectively [18].

We are interested in the simulation of multicomponent flow problems. For the equations, we take a popular approach by considering the Euler Eqs. (1) as a model system, see [33, 74] for the use of other governing equations. Our goal here is to derive computational models that may prevent pressure oscillations near the interfaces, when solving the problem numerically with standard shock-capturing methods.

*2.1. Preliminary. To begin, suppose that there are m different fluid components in a*
grid cell, and each of them occupies a distinct region with a volume-fraction function
*Y*^{(i)}*in relation to it, for i= 1, 2, . . . , m. Here by the standard assumptions we have Y*^{(i)}

∈ [0, 1] andP*m*

*i*=1*Y*^{(i)}*= 1. Suppose that for each component i the state variables such as*
ρ^{(i)}*, u*^{(i)}*, p** ^{(i)}*, γ

^{(i)}*, and p*

^{(i)}_{∞}are known a priori. The objective is to define the mixture of the

*pressure p as well as the conserved variablesρ, ρu, and ρE in a consistent manner within*the cell for the Euler Eqs. (1). Note this step is necessary when we initialize the data for computations.

To accomplish this, we follow a common practice by setting*ρ, ρu, and ρE as a volume-*
weighted sum over the set of componentsρ* ^{(i)}*, ρ

^{(i)}*u*

*, andρ*

^{(i)}

^{(i)}*E*

*, for each separately,*

^{(i)}ρ =
X*m*

*i*=1

*Y** ^{(i)}*ρ

*,*

^{(i)}*ρu =*X

*m*

*i*=1

*Y** ^{(i)}*ρ

^{(i)}*u*

*,*

^{(i)}*ρE =*X

*m*

*i*=1

*Y** ^{(i)}*ρ

^{(i)}*e*

*+1*

^{(i)}2*Y** ^{(i)}*ρ

^{(i)}*u*

^{2}, (3)

*where u is the velocity mixture that can be computed easily by*

*u* =*ρu*
ρ =

P*m*

*i*=1*Y** ^{(i)}*ρ

^{(i)}*u*

*P*

^{(i)}*m*

*i*=1*Y** ^{(i)}*ρ

*.*

^{(i)}*To find the pressure mixture p, we need to use the equation of state (2). From (3), after*
some simple algebra, we find

*p+ γ p*∞

γ − 1 *= ρe =*
X*m*

*i*=1

*Y** ^{(i)}*ρ

^{(i)}*e*

*= X*

^{(i)}*m*

*i*=1

*Y** ^{(i)}*
Ã

*p** ^{(i)}*+ γ

^{(i)}*p*

_{∞}

*γ*

^{(i)}*− 1*

^{(i)}!

. (4)

*Note this gives us one equation that is not only for the mixture of p but also forγ and p*_{∞}.
*Obviously, we need to choose two supplementary conditions so as to have p,γ , and p*∞

defined well, leading to the full agreement of (4).

The basic idea of our approach is quite simple. We first split (4) into two parts by setting the terms

*p*
γ − 1 =

X*m*
*i*=1

*Y*^{(i)}*p*^{(i)}

γ* ^{(i)}*− 1 and

*γ p*

_{∞}γ − 1 =

X*m*
*i*=1

*Y** ^{(i)}*γ

^{(i)}*p*

^{(i)}_{∞}

γ* ^{(i)}*− 1 ; (5)

212 KEH-MING SHYUE

this provides us with one condition right away. We then impose the following condition to the computation ofγ ,

1 γ − 1=

X*m*
*i*=1

*Y*^{(i)}

γ* ^{(i)}*− 1. (6a)

Clearly when*γ is known, from (5), it is an easy matter to determine the unknowns p and*
*p*_{∞}. The results are

*p*=
Ã _{m}

X

*i*=1

*Y*^{(i)}*p** ^{(i)}*
γ

*− 1*

^{(i)}!,Ã * _{m}*
X

*i*=1

*Y** ^{(i)}*
γ

*− 1*

^{(i)}!

(6b)

and

*p*_{∞}=
Ã _{m}

X

*i*=1

*Y** ^{(i)}*γ

^{(i)}*p*

^{(i)}_{∞}γ

*− 1*

^{(i)}!,Ã 1+

X*m*
*i*=1

*Y** ^{(i)}*
γ

*− 1*

^{(i)}!

. (6c)

*Notice that in case each of the partial pressures p** ^{(i)}* is in equilibrium within a grid cell

*the pressure p acquired from (6b) would remain in equilibrium also, i.e., p= p*

*, for*

^{(i)}*i= 1, 2, . . . , m. This nice property on the mixture of p is in fact the main reason why we*use (6a), but not other ways to computeγ , see [14, 55] for a more rigorous derivation. In

*addition, we get the total pressure p*=P

*m*

*i*=1*Y*^{(i)}*p** ^{(i)}*when the grid cell is composed of the
same set ofγ

^{(i)}*but with different p*

^{(i)}*. We use (6c) for the mixture of p*

_{∞}in order to make sure that (4) is handled in a consistent manner.

*It is worthwhile to mention that in case Y** ^{(i)}*represents a mass-fraction function of the

*i th component with*ρ

^{(i)}*= ρY*

*, we would have Eq. (4) replaced by*

^{(i)}*p+ γ p*_{∞}

γ − 1 *= ρe =*
X*m*

*i*=1

ρ^{(i)}*e** ^{(i)}*=
X

*m*

*i*=1

Ã

*p** ^{(i)}*+ γ

^{(i)}*p*

^{(i)}_{∞}γ

*− 1*

^{(i)}! .

In this instance, it is a standard approach that setsγ according to (6a) (cf. [26, 39]).

Analogously, by following the same procedures as for the volume-fraction case, we find
*the results of p and p*_{∞},

*p* =
Ã _{m}

X

*i*=1

*p** ^{(i)}*
γ

*− 1*

^{(i)}!,Ã * _{m}*
X

*i*=1

*Y** ^{(i)}*
γ

*− 1*

^{(i)}! ,

*p*_{∞}=
Ã _{m}

X

*i*=1

γ^{(i)}*p*^{(i)}_{∞}
γ* ^{(i)}*− 1

!,Ã 1+

X*m*
*i*=1

*Y** ^{(i)}*
γ

*− 1*

^{(i)}! .

2.2. *γ -based model. It is clear that the motion of the fluid mixtures, obtained from*
above, is governed by the Euler Eqs. (1). In the development of our multicomponent model,
*it is of great importance to first consider the case of an interface only problem where both*
*of the pressure p and velocity u are constants in the domain, while the other variables such*
as*ρ, γ , and P*_{∞}are having jumps across some interfaces.

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 213

To do so, we start out to write (1) in the following non-conservative form,

∂ρ

*∂t* *+ u*∂ρ

*∂x* + ρ*∂u*

*∂x* = 0,

*∂u*

*∂t* *+ u∂u*

*∂x* +1
ρ

*∂p*

*∂x* = 0,

∂

*∂t(ρe) +* ∂

*∂x(ρeu) + p∂u*

*∂x* = 0,
and obtain easily equations describing the motion of the interfaces as

∂ρ

*∂t* *+ u*∂ρ

*∂x* = 0, (7a)

∂

*∂t(ρe) + u* ∂

*∂x(ρe) = 0.* (7b)

With this, it is clear that for an interface the densityρ as well as the total internal energy
*ρe is evolved by the linear transport Eqs. (7a) and (7b).*

*To see how the pressure p would retain in equilibrium as it should be for this problem,*
we insert the equation of state (2) into (7b), and have

∂

*∂t*

µ*p+ γ p*_{∞}
γ − 1

¶
*+ u* ∂

*∂x*

µ*p+ γ p*_{∞}
γ − 1

¶

= 0. (8)

By expanding (8), we may therefore write that equation as µ 1

γ − 1

¶·*∂p*

*∂t* *+ u∂p*

*∂x*

¸
*+ p*

·∂

*∂t*
µ 1

γ − 1

¶
*+ u* ∂

*∂x*
µ 1

γ − 1

¶¸

+

·∂

*∂t*
µ*γ p*_{∞}

γ − 1

¶

*+ u* ∂

*∂x*
µ*γ p*_{∞}

γ − 1

¶¸

= 0.

*Now the requirement that p be in equilibrium leads to the equation of the form*

*p*

·∂

*∂t*
µ 1

γ − 1

¶
*+ u* ∂

*∂x*
µ 1

γ − 1

¶¸

+

·∂

*∂t*
µ*γ p*_{∞}

γ − 1

¶
*+ u* ∂

*∂x*
µ*γ p*_{∞}

γ − 1

¶¸

= 0.

*Since this equation should hold for any p in the physical space, it implies that the terms in*
bracket of the above equation should be vanished simultaneously, yielding a system of two
equations

∂

*∂t*
µ 1

γ − 1

¶
*+ u* ∂

*∂x*
µ 1

γ − 1

¶

= 0

∂

*∂t*
µ*γ p*_{∞}

γ − 1

¶
*+ u* ∂

*∂x*
µ*γ p*_{∞}

γ − 1

¶

= 0.

(9)

Note that these are the evolution equations that should be satisfied for the material-dependent
variables*γ and p*_{∞}in order to have the correct pressure behavior in (8) for the interface. For
convenience, we call (9) the effective equations of the problem, where the initial condition
of the equations are provided, for example, by (6a) and (6c) in a respective manner.

214 KEH-MING SHYUE

Of course, intuitively, there are many other ways that the effective equations can be
*rewritten, while still getting the correct pressure from (8) for this interface only problem.*

One simple example among them is to write (9) as

∂

*∂t*
µ 1

γ − 1

¶ + ∂

*∂x*
µ 1

γ − 1*u*

¶

= 0

∂

*∂t*
µ*γ p*∞

γ − 1

¶ + ∂

*∂x*
µ*γ p*∞

γ − 1*u*

¶

= 0 ;

(10)

*this is a legitimate one to use, for the velocity u is a constant for the problem (see*
*Eqs. (25a)–(25c) also for other examples when p*_{∞}= 0). But since in general we are inter-
ested in shock wave problems as well, with the stiffened gas equation of state (2), we should
only take the effective equations in a form that do not vary their solutions across both of
shocks and rarefaction waves. For this reason, it rules out immediately the use of (10) as the
effective equations. When further taking the numerical aspect of the model equations into
consideration, it turns out that with the full Euler Eqs. (1) the effective equations in a form
of (9) are the proper ones to use for practical multicomponent problems; see Section 3 for
the analysis of a numerical method that approximates the model equations, and Section 4
for numerical examples.

In summary, the model equations we propose to solve multicomponent problems with the stiffened gas equation of state are

∂ρ

*∂t* + ∂

*∂x(ρu) = 0*

∂

*∂t(ρu) +* ∂

*∂x(ρu*^{2}*+ p) = 0*

∂

*∂t(ρE) +* ∂

*∂x*[*(ρE + p)u] = 0*

∂

*∂t*
µ 1

γ − 1

¶
*+ u* ∂

*∂x*
µ 1

γ − 1

¶

= 0

∂

*∂t*
µ*γ p*∞

γ − 1

¶
*+ u* ∂

*∂x*
µ*γ p*∞

γ − 1

¶

= 0.

(11)

Here in this system the first three equations are simply the Euler Eqs. (1) that are used
to make certain the conservation of*ρ, ρu, and ρE, while the last two equations are the*
effective Eqs. (9) of the problem, that are introduced to ensure the correct mixing behavior
of the variables*γ and p*∞on the interfaces.

Note the model Eqs. (11) is not written in the full conservation form, but is rather a
quasi-conservative system of equations. For shock wave problems, this poses no problem
at all, because for stiffened gases*γ and p*_{∞}remain unchanged across genuinely nonlinear
waves such as shocks or rarefactions, and hence we have the usual Riemann invariants and
Rankine–Hugoniot jump conditions across the rarefaction and shock waves, respectively
(cf. [26, 68]). At linearly degenerate waves such as interfaces where there may be jumps in
*γ and p*∞, this again gives no problem since as we have seen in the previous discussion we
would have the desired pressure equilibrium when the model is in use.

Of particular importance is the case of shock wave and interface interaction. It is known that, due to the nonlinearity of the problem, the pressure across the interface will be quite

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 215

different before and right after the wave interaction [7, 18]. Somewhat surprisingly, we find
no major problem in using the model as well, see Section 4 for a representative numerical
test. We emphasize that in this case*γ and p*_{∞}are transported in a passive manner according
to (9) along with the interface, and the variables such as*ρ, u, and p are dealt with in a*
conservative way from the conservation part of Eqs. (11) as usual.

With all these in mind, it should be sensible to use the proposed model as for practical computations. The numerical method to be described in Section 3 is a consistent approxi- mation of the model that gives excellent results for a wide variety of problems as illustrated in Section 4. For convenience, we call this model aγ -based model to be distinct from the other one presented below. Note in [1] Abgrall used a slightly different starting point, but nevertheless he obtained the same set of model equations for polytropic gas problems when

*p*_{∞}= 0.

*2.3. Volume-fraction model. It should be noted that we may reformulate the above*γ -
based model as a volume-fraction or a mass-fraction model that is also applicable to many
multicomponent problems. To demonstrate the basic idea, we consider the case in using
volume-fraction functions as an example.

Similar to the approach used in the construction of our γ -based model, we look for
*effective equations that may preserve the pressure equilibrium for an interface only problem.*

Here with the volume-fraction notion of the states 1*/(γ −1) and γ p*∞/(γ −1) being defined
by (6a) and (6c), the key step is to replace (9) by

∂

*∂t*
Ã _{m}

X

*i*=1

*Y** ^{(i)}*
γ

*− 1*

^{(i)}!
*+ u* ∂

*∂x*
Ã _{m}

X

*i*=1

*Y** ^{(i)}*
γ

*− 1*

^{(i)}!

= 0,

∂

*∂t*
Ã _{m}

X

*i*=1

*Y** ^{(i)}*γ

^{(i)}*p*

_{∞}

*γ*

^{(i)}*− 1*

^{(i)}!
*+ u* ∂

*∂x*
Ã _{m}

X

*i*=1

*Y** ^{(i)}*γ

^{(i)}*p*

^{(i)}_{∞}γ

*− 1*

^{(i)}!

= 0.

*After regrouping terms, we find the transport equation for each volume fraction Y** ^{(i)}*,

*∂Y*^{(i)}

*∂t* *+ u∂Y*^{(i)}

*∂x* = 0 *for i* *= 1, 2, . . . , m.* (12)
*As before when the set of Y** ^{(i)}*is known, we may therefore compute

*γ and p*

_{∞}from (6a)

*and (6c). In effect, in a volume-fraction model, instead of using (9) we use (12) (m of them)*

*as the effective equations of the problem. Note that in case Y*

*stands for the mass-fraction*

^{(i)}*of the component i , we find the effective equations of a mass-fraction model that is of the*same form as in (12). Because of the close connection between the two models, we devote our discussion to one of the models only, namely to the volume-fraction model.

For completeness, we write down the full set of equations for this volume-fraction model,

∂ρ

*∂t* + ∂

*∂x(ρu) = 0*

∂

*∂t(ρu) +* ∂

*∂x(ρu*^{2}*+ p) = 0*

∂

*∂t(ρE) +* ∂

*∂x*[*(ρE + p)u] = 0*

*∂Y*^{(i)}

*∂t* *+ u∂Y*^{(i)}

*∂x* = 0 *for i* *= 1, 2, . . . , m;*

(13)

216 KEH-MING SHYUE

*this gives us totally m*+ 3 equations to be solved. Since the derivation of the volume-
fraction model comes closely out of theγ -based model, it can be shown that this model is
as effective as theγ -based model. But for general multicomponent problems, the γ -based
model is the preferred one to use, because the basic equations for the model stay as five,
see (11), irrespective of the number of components involved in the problem.

**3. NUMERICAL APPROXIMATION OF MODEL EQUATIONS**

We use the high-resolution wave propagation method to compute solutions for our mul-
*ticomponent models introduced in Section 2. The method is a variant of the fluctuation-*
*and-signal formulation of Roe [59, 60] in that we solve the Riemann problem at each cell*
interface, and use the resulting waves (i.e., discontinuities propagating at constant speeds)
to update the solutions in neighboring grid cells (cf. [40, 41]).

For simplicity, we describe the method on a uniform grid with fixed mesh spacing*1x,*
but the method can be extended quite easily to a nonuniform and time-varying grid as well
*[44]. We use a standard finite-volume formulation that the value U*^{n}* _{j}* ∈ R

*approximates the*

^{m}*cell average of the solution over the grid cell [x*

*j*

*, x*

*j*+1

*] at time t*

*n*. The time step is denoted by

*1t.*

In this setup, a first order accurate version of the method in wave-propagation form is a Godunov-type scheme that can be written as

*U*^{n}_{j}^{+1}*= U*^{n}*j* − *1t*
*1x*

X*m*
*k*=1

£(λ^{−}*k*W*k*)^{n}*j*+1+ (λ^{+}*k*W*k*)^{n}*j*

¤, (14)

whereλ*k*∈ R and W*k*∈ R* ^{m}*are solutions obtained from solving Riemann problems at cell

*interfaces x*

*j*

*and x*

*j*+1, λ

^{−}= min(λ, 0), and λ

^{+}= max(λ, 0) (cf. [26]). It is easy to see that the method belongs to a class of upwind schemes, and as it will be shown below that the method is quasi-conservative in the sense that when applying the method to our multicomponent models not only the conservation equations but also the transport equations are approximated in a consistent manner by the method with the chosen Riemann solver.

Concerning stability, it is observed numerically that the method is stable and convergent under mesh refinement provided that the waves in the method affect only the cells adjacent to the interface during the time step; see Section 4 for numerical examples and also the results shown in [43].

*3.1. Method withγ -based model. Consider the γ -based model as an example. We solve*
*the Riemann problem at each cell interface x** _{j}*that consists of (11) with piecewise constant

*data U*

^{n}

_{j}_{−1}

*and U*

^{n}*on the left and on the right of the interface. Rather than computing the exact solution to this Riemann problem, which can be done by iterative procedures (cf. [12, 15, 54]) but is rather expensive, we use a generalized version of the approximate Riemann solver of Roe (cf. [58] and below) in most instances. This is much more efficient to compute than the exact Riemann solution, and provides a very accurate approximation of solution for smooth flows and also for moderate-strength shock waves. As long as the equation of state is not too stiff across the interfaces (which is the application considered here), the solution of the Roe Riemann solver gives a proper resolution of the contact discontinuity to be used in the method (14) for numerical approximation also. It is important to mention that in case we have a stringent set of data that involves strong shock waves and (or) stiff equation of states, we find it is advisable to use the exact Riemann solver so as to properly deal with*

_{j}MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 217

the nonlinear effect of the solution structure; see [67] for an example and also [20, 45, 56]

for other examples that the exact Riemann solver should be in use.

To implement Roe’s approximate Riemann solver, we first write (11) as a quasi-linear system of equations

*∂q*

*∂t* *+ A(q)∂q*

*∂x* = 0 (15)

with

*q*=

ρ
*ρu*
*ρE*

1
γ −1
*γ p*∞

γ −1

, *A(q) =*

0 1 0 0 0

³γ −3 2

´

*u*^{2} *(3 − γ )u* γ − 1 χ 1− γ

³γ −1 2

´

*u*^{3}*− uH H − (γ − 1)u*^{2} *γ u* *χu (1 − γ )u*

0 0 0 *u* 0

0 0 0 0 *u*

,

*where H* *= E + (p/ρ) is the enthalpy, and χ = −ι(γ − 1)*^{2}with*ι = p/(γ − 1). We then*
solve a linear problem

*∂q*

*∂t* *+ ˆA(q**L**, q**R*)*∂q*

*∂x* = 0 (16)

with initial data

*q(x, 0) =*

½*q**L* *for x left to the interface*
*q**R* *for x right to the interface*,

where ˆ*A(q**L**, q**R*) is a constant matrix that depends on the initial data and is a local lineariza-
*tion of the matrix A about an average state. Here as it is often done in many other Roe*
solvers (cf. [10, 23, 26]), we want to seek an average state that the difference of the fluxes
in the conservation part of Eqs. (11) are equal to the respective first order approximation of
the flux differences. That is,

1F* ^{(i)}*= (F

*R*− F

*L*)

^{(i)}*= [ ˆA(q*

*L*

*, q*

*R*

*)(q*

*R*

*− q*

*L*)]

^{(i)}*= [ ˆA(q*

*L*

*, q*

*R*

*)1q]*

*, (17)*

^{(i)}*for i*= 1, 2, 3, where F ∈ R

^{3}is the usual definition of the fluxes for the Euler Eqs. (1), and1F

^{(i)}*is the i th component of*1F.

To accomplish the relation in (17), by taking a similar approach employed in [26] for
real gases (cf. [23] also), we find it is sufficient to get the average states for variables such
as ˆ*u, ˆH, ˆγ , ˆι, and set the matrix ˆA(q**L**, q**R**) = A(ˆu, ˆH, ˆγ , ˆι). The results are*

ˆ

*u* =√ρ^{L}*u**L*+ √ρ*R**u**R*

√ρ*L*+ √ρ*R*

, *H*ˆ =√ρ^{L}*H**L*+ √ρ*R**H**R*

√ρ*L*+ √ρ*R*

,

(18) 1

γ − 1ˆ =√ρ*L*

³ 1
γ*L*−1

´+ √ρ*R*

³ 1
γ*R*−1

´

√ρ*L*+ √ρ*R*

, ˆι = √ρ* ^{L}*ι

*L*+ √ρ

*R*ι

*R*

√ρ*L*+ √ρ*R*

.

Note that in the current derivation of the average states ˆγ and ˆι are chosen so that the expression

*1p = ( ˆγ − 1)*^{2}

· 1

γ − 1ˆ 1ι − ˆι1 µ 1

γ − 1

¶¸

218 KEH-MING SHYUE

is satisfied approximately, and ˆ*p is computed by*( ˆγ −1)ˆι (see [23, 24] for a related discussion
to the case of real gases). As shown in Section 4, we find good results with the use of this
set of average states defined in (18).

The solution of the linear problem (16) consists of five discontinuities propagating at
constant speeds. The jump across each discontinuity is a multiple of the eigenvector of the
matrix ˆ*A, and the propagating speed is the corresponding eigenvalue. We thus have*

*1q = q**R**− q**L* =
X5
*k*=1

αˆ*k**ˆr**k*, (19)

*where ˆr**k**is the kth eigenvector of ˆA with*

*ˆr*1=

1
ˆ
*u− ˆc*
*H*ˆ *− ˆuˆc*

0 0

*, ˆr*2=

1
ˆ
*u*

1
2*u*ˆ^{2}

0 0

*, ˆr*3=

1
ˆ
*u+ ˆc*
*H*ˆ *+ ˆuˆc*

0 0

*, ˆr*4=

0
0
*ˆp*
1
0

*, ˆr*5=

0 0 1 0 1

,

(20)
and ˆ*A ˆr**k*= ˆλ*k**ˆr**k*with ˆλ*k*the corresponding eigenvalue,

λˆ1*= ˆu − ˆc,* λˆ2*= ˆu,* λˆ3*= ˆu + ˆc,* λˆ4= ˆλ5*= ˆu,* (21)
*where ˆc*=p

*( ˆγ − 1)( ˆH − ˆu*^{2}/2) is the speed of sound. The scalar ˆα*k*gives the strength
across the discontinuity that can be determined easily from (19). We find

αˆ2= γ − 1ˆ
*ˆc*^{2}

£*( ˆH − ˆu*^{2}*)1q*^{(1)}*+ ˆu1q*^{(2)}*− 1q*^{(3)}*+ ˆp1q*^{(4)}*+ 1q*^{(5)}¤

= 1ρ − *1p*
*ˆc*^{2}
αˆ3= 1

*2ˆc*

£*(ˆc − ˆu)1q*^{(1)}*+ 1q*^{(2)}*− ˆc ˆα*2

¤=*1p + ˆρˆc1u*
*2ˆc*^{2}
αˆ1*= 1q*^{(1)}− ˆα2− ˆα3= *1p − ˆρˆc1u*

*2ˆc*^{2}
αˆ4*= 1q*^{(4)}= 1

µ 1 γ − 1

¶

αˆ5*= 1q*^{(5)}= 1
µ*γ p*_{∞}

γ − 1

¶ ,

(22)

where ˆρ = √ρ*L*ρ*R*.

Notice that in this Riemann solution there exist three discontinuities propagating at the
same speed, ˆλ2= ˆλ4= ˆλ5*= ˆu. For practical purposes, we may view these discontinuities*
as a single one with the operator ˆW2defined by combining all the jumps across the ˆλ2wave
family, i.e., set ˆW2 = ˆα2*ˆr*2+ ˆα4*ˆr*4+ ˆα5*ˆr*5. Clearly, doing so removes the effect of the
λˆ4 and ˆλ5wave families to the solution. With this notation, we also write ˆW*k* = ˆα*k**ˆr**k*to
*represent the jump across the k*= 1 and 3 waves. Wave propagation methods are based on
using these propagating discontinuities to update the cells averages in the cells neighboring
each interface.

*We now show that in the case of an interface only problem the numerical solution obtained*
using the method would be free of oscillations, and in particular the pressure would remain

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 219

in equilibrium. Without loss of generality, we consider a single Riemann problem that at
*cell interface x**j*the initial data are picked up so that there are no jumps in both the pressure
and velocity, i.e.,*1p*^{n}*j* *= 0 and 1u*^{n}*j* = 0, but is otherwise for the other variables such as
*ρ, γ , and p*∞. With this initial data, the Riemann solution consists of only a single contact
discontinuity with the wave strength ˆα*k j*computed by (22) as

αˆ*1 j*= ˆα*3 j* = 0, αˆ*2 j* = 1ρ^{n}*j*, αˆ*4 j* = 1
µ 1

γ − 1

¶*n*

*j*

, αˆ*5 j*= 1
µ*γ p*_{∞}

γ − 1

¶*n*

*j*

,

*and the associated eigen-structure ˆr**k j*, ˆλ*k j*computed by (20) and (21), respectively.

Now suppose that the wave moves to the right of the interface ˆ*u**j* > 0. To take account
*of the effect of this wave, according to (14), the cell average U*^{n}_{j}^{+1}should be evaluated by

ρ^{n}*j*^{+1} = ρ^{n}*j* − *1t*

*1xu*ˆ*j*1ρ^{n}*j*, (23a)

*(ρu)*^{n}*j*^{+1} *= (ρu)*^{n}*j*− *1t*

*1xu*ˆ*j*[ ˆ*u*1ρ* ^{n}*]

*j*, (23b)

*(ρE)*^{n}*j*^{+1} *= (ρE)*^{n}*j*− *1t*
*1xu*ˆ_{j}

·1

2*u*ˆ^{2}1ρ^{n}*+ ˆp1*
µ 1

γ − 1

¶*n*

+ 1
µ*γ p*∞

γ − 1

¶*n*¸

*j*

, (23c) µ 1

γ − 1

¶*n*+1

*j*

= µ 1

γ − 1

¶*n*

*j*

− *1t*
*1xu*ˆ*j*1

µ 1 γ − 1

¶*n*

*j*

, (23d)

µ *γ p*∞

γ − 1

¶*n*+1

*j*

=
µ*γ p*∞

γ − 1

¶*n*

*j*

− *1t*
*1xu*ˆ*j*1

µ*γ p*∞

γ − 1

¶*n*

*j*

. (23e)

*It follows from (23a) and (23b) that we have the expected state of u*^{n}_{j}^{+1}*= u*^{n}*j*where by (18),
ˆ

*u*_{j}*= u*^{n}*j*. With this result, Eq. (23c) can be simplified to

*(ρe)*^{n}*j*^{+1}*= (ρe)*^{n}*j*− *1t*
*1xu*ˆ_{j}

·
*ˆp1*

µ 1 γ − 1

¶*n*

+ 1
µ*γ p*∞

γ − 1

¶*n*¸

*j*

,

or alternatively to
µ*p+ γ p*∞

γ − 1

¶*n*+1

*j*

=

µ*p+ γ p*∞

γ − 1

¶*n*

*j*

− *1t*
*1xu*ˆ_{j}

·
*ˆp1*

µ 1 γ − 1

¶*n*

+ 1
µ*γ p*∞

γ − 1

¶*n*¸

*j*

,

when using the equation of state (2). Notice that in case (23e) is applied to the above equation, we have further simplification

µ *p*

γ − 1

¶*n*+1
*j*

=

µ *p*

γ − 1

¶*n*
*j*

− *1t*
*1xu*ˆ*j*

·
*ˆp1*

µ 1 γ − 1

¶*n*¸

*j*

,

*and find the pressure equilibrium of the computed solution p*^{n}_{j}^{+1} *= p*^{n}*j*, when Eq. (23d)
*is satisfied along with the fulfillment of the condition ˆp*_{j}*= p*^{n}*j* (see (18) and the comment
thereafter).

As to the behavior of the solutionsρ^{n}*j*^{+1}, γ*j*^{n}^{+1}, and*(p*∞)^{n}*j*^{+1}, from (23a), (23d), (23e), it is
easy to derive the modified equations for each of them individually, and show monotonicity

220 KEH-MING SHYUE

of the results obtained using the method (cf. [42]). See results shown in Section 4 for numerical verification of this statement also.

*In practical applications, there may be some other waves which come into the j th cell*
*and affect the cell average U*^{n}_{j}^{+1}as well. Suppose that we are taking Riemann data at cell
*interface x*_{j}_{+1}so that there is a 1-wave (shock or rarefaction) propagating to the left of the
*interface. In this method (14), we update the cell average U*^{n}_{j}^{+1}by

ρ^{n}*j*^{+1} := ρ^{n}*j*^{+1}− *1t*

*1x(ˆu − ˆc)**j*+1

·*1p*^{n}*− ˆρˆc1u*^{n}*2ˆc*^{2}

¸

*j*+1,
*(ρu)*^{n}*j*^{+1} :*= (ρu)*^{n}*j*^{+1}− *1t*

*1x(ˆu − ˆc)**j*+1

·
*(ˆu − ˆc)*

µ*1p*^{n}*− ˆρˆc1u*^{n}*2ˆc*^{2}

¶¸

*j*+1,
*(ρE)*^{n}*j*^{+1} :*= (ρE)*^{n}*j*^{+1}− *1t*

*1x(ˆu − ˆc)**j*+1

·

*( ˆH − ˆuˆc)*

µ*1p*^{n}*− ˆρˆc1u*^{n}*2ˆc*^{2}

¶¸

*j*+1.
We note that due to a fundamental property of the Roe solver, i.e., the relation in (17), this
is a conservative update of the numerical solutionsρ^{n}*j*^{+1}*, (ρu)*^{n}*j*^{+1}, and*(ρE)*^{n}*j*^{+1}. We get a
*new pressure p*^{n}_{j}^{+1}by using the equation of state (2),

*p*^{n}_{j}^{+1}=¡

γ*j*^{n}^{+1}− 1¢µ
*ρE −*1

2*ρu*^{2}

¶*n*+1

*j*

*− (γ p*∞)^{n}*j*^{+1}.

One advantage of using the wave propagation form is that we are able to handle each wave in turn, and there is no need to compute fluxes and make distinction between the waves (cf. [67]). Extension of the method to higher order accuracy, and in particular to a high-resolution version of the wave propagation scheme, follows easily by incorporating limited slopes and by constructing piecewise linear profiles to the method; see [43, 44]

*for the detail. It is not difficult to show that for the interface only problem we again have*
the required pressure equilibrium that is independent of the limiter being employed to the
high-resolution method. Moreover, we obtain a better resolution of the result as compared
to the first order result; see Section 4 for numerical examples.

*3.2. Method with volume-fraction model. We now consider the use of the Riemann*
solutions associated with the volume-fraction model (13) to the wave propagation method
(14). For simplicity in description, we are concerned with a two-component flow problem
that the quasi-linear system of equations in (15) are defined with

*q*=

ρ
*ρu*
*ρE*
*Y*^{(1)}
*Y*^{(2)}

, *A(q) =*

0 1 0 0 0

³γ −3 2

´

*u*^{2} *(3 − γ )u* γ − 1 χ^{(1)} χ^{(2)}

³_{γ −1}

2

´

*u*^{3}*− uH H − (γ − 1)u*^{2} *γ u* χ^{(1)}*u* χ^{(2)}*u*

0 0 0 *u* 0

0 0 0 0 *u*

,

whereχ^{(i)}*= (1 − γ )(p + γ*^{(i)}*p*^{(i)}_{∞})/(γ* ^{(i)}*− 1). We use the Roe solver as usual in that we
solve the linear problem (16) with the Roe-type matrix ˆ

*A(q*

*L*

*, q*

*R*

*) = A(ˆu, ˆH, ˆι, ˆY*

^{(1)}

*, ˆY*

^{(2)})

*depending on the initial data q*

*L*

*and q*

*R*. Here as in the case of the γ -based model we

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 221

compute the average states ˆ*u, ˆH, ˆι by (18), and set the average state ˆY** ^{(i)}*by the standard

“Roe-averaging” approach also with

*Y*ˆ* ^{(i)}*= √ρ

^{L}*Y*

_{L}*+ √ρ*

^{(i)}*R*

*Y*

_{R}

^{(i)}√ρ*L*+ √ρ*R*

*for i* = 1, 2.

When defining ˆ*A in this way, it is easy to check the satisfaction of the fundamental relation*
in (17).

The result of this linearized Riemann problem takes a rather similar wave structure as to the form shown in (19)–(22) of theγ -based model. Without writing the full solution here, we mention that the only difference between the representation of these two solutions is to the jump across the ˆλ4and ˆλ5families. Here we use the following expressions instead,

αˆ4*= 1Y*^{(1)}, αˆ5 *= 1Y*^{(2)}, *ˆr*4 =

0 0

ˆ
*p*+ γ^{(1)}*p*^{(1)}_{∞}

γ^{(1)}− 1

1 0

, *ˆr*5=

0 0

*ˆp*+ γ^{(2)}*p*^{(2)}_{∞}
γ^{(2)}− 1

0 1

.

As in the case performed for theγ -based model, we now analyze the pressure obtained
*using the method (14) with the volume-fraction model for the interface only problem con-*
sidered there. In this case, it is enough to begin looking at the update of the internal energy
in the form,

µ*p+ γ p*_{∞}
γ − 1

¶*n*+1
*j*

=

µ*p+ γ p*_{∞}
γ − 1

¶*n*
*j*

− *1t*
*1xu*ˆ*j*

" _{2}
X

*i*=1

*ˆp*+ γ^{(i)}*p*^{(i)}_{∞}
γ* ^{(i)}*− 1 1¡

*Y** ^{(i)}*¢

*n*

#

*j*

.

Note the above equation is just a simplified version of the update of the total energy by taking
*into account of the fact that u*^{n}_{j}^{+1}*= u*^{n}*j*. When substituting the volume-fraction relations in
(6a) and (6c) for*γ and p*_{∞}, we may write the equation as

Ã _{2}
X

*i*=1

*Y*^{(i)}*p*

γ* ^{(i)}*− 1+

*Y*

*γ*

^{(i)}

^{(i)}*p*

_{∞}

*γ*

^{(i)}*− 1*

^{(i)}!*n*+1

*j*

=
Ã _{2}

X

*i*=1

*Y*^{(i)}*p*

γ* ^{(i)}*− 1 +

*Y*

*γ*

^{(i)}

^{(i)}*p*

^{(i)}_{∞}γ

*− 1*

^{(i)}!*n*

*j*

−*1t*
*1xu*ˆ*j*

" _{2}
X

*i*=1

*ˆp*+ γ^{(i)}*p*^{(i)}_{∞}
γ* ^{(i)}*− 1 1¡

*Y** ^{(i)}*¢

*n*

#

*j*

.

*After a simple manipulation, we find p*^{n}_{j}^{+1}*= p*^{n}*j*, only if ˆ*p*_{j}*= p*^{n}*j*and the establishment of
*the difference equations for the volume fraction Y** ^{(i)}*,

¡*Y** ^{(i)}*¢

*n*+1

*j* =¡

*Y** ^{(i)}*¢

*n*

*j*−

*1t*

*1xu*ˆ* _{j}*1¡

*Y*

*¢*

^{(i)}*n*

*j* *for i* = 1, 2.

It is easy to see that the latter two requirements are guaranteed by the method.

Note for this two-component flow problem we have described the method uses two
*different volume-fraction functions Y*^{(1)}*and Y*^{(2)}*. In practice, due to the fact that Y*^{(1)}+
*Y*^{(2)}*= 1, we may simply use a volume-fraction function, say Y*^{(1)}, to the model, and set the

222 KEH-MING SHYUE

*value of the other one, Y*^{(2)}*= 1 − Y*^{(1)}, for example. When modeling in this way, it gives
a more robust method to the computation of two-component flow problems than the use of
theγ -based model to the method.

**4. NUMERICAL RESULTS IN ONE DIMENSION**

We now present results to validate our multicomponent algorithm with the Roe solver described in Section 3 in one dimension.

EXAMPLE4.1. *We consider an interface only problem that the solution of a Riemann*
problem consists of a single contact discontinuity in gas dynamics. We take two sets of data
for numerical experiments. In the first set, we have a polytropic gas and use two constant
states as

ρ
*u*
*p*
γ
*p*_{∞}

*L*

=

1 1 1 1.4

0

and

ρ
*u*
*p*
γ
*p*_{∞}

*R*

=

0.125

1 1 1.2

0

, (24)

while in the second set we use the same data as in the first set with the exception that on
*the state R a stiffened gas with*γ*R* *= 4 and (p*∞)*R* *= 1 is employed instead. Here L is the*
*state used for x∈ [0, 0.2) and R is the state used for x ∈ [0.2, 1].*

Results for the first set of data are shown in Fig. 1 where we have employed theγ -based model together with both the first order and high-resolution wave propagation methods to the computations. For comparison, we also include results obtained using three different effective equations to the simulations. They are the conservation equations with either

∂

*∂t*(ργ ) + ∂

*∂x(ργ u) = 0,* (25a)

or

∂

*∂t*
µ ρ

γ − 1

¶ + ∂

*∂x*
µ *ρu*

γ − 1

¶

= 0, (25b)

and the primitive equation

∂γ

*∂t* *+ u*∂γ

*∂x* = 0. (25c)

Note in the cases of (25a) and (25b), we have two model systems that are in the full conservation form; see [34] for a similar consideration.

From the figure, we clearly observe pressure fluctuations in the solutions when employing
any of the equations in (25a)–(25c), but not theγ -based model where (9) is considered, to
*the method. By following a similar analysis conducted in Section 3 for an interface only*
problem, we may explain the observed error behavior in pressure as being the failure to
approximate the energy Eq. (7b) in a consistent manner when those equations are in use
(cf. [34]). Note other variables in the solution are affected by this error also as time evolves.

Numerical evidences suggest however that these errors decrease as the mesh is refined, and

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 223

**FIG. 1.** Comparison plots of four different multicomponent models for the Euler equations with data (24) at
*time t* = 0.12. (a) Results using the first order wave propagation method. (b) Results using the high-resolution
wave propagation method with the “minmod” slope limiter. In each figure the solid line is the exact solution and
the points show the computed solution with 100 mesh points.

224 KEH-MING SHYUE

**FIG. 2.** *High-resolution results for an interface only problem using the same initial data as in Fig. 1 with the*
exception thatγ*R**= 4 and (p*∞)*R*= 1 are used in the current computation. The solid line is the exact solution and
the points show the computed solution with 100 mesh points.

the rate of convergence of the error in the 1-norm is about the order of accuracy of the method that is employed to the computation.

Figure 2 shows results of a run using the second set of data where only the high-resolution solutions with the γ -based model are presented. Notice that the pressure and also the velocity remain at constant states for this stiffened gas simulation. Because the errors become too erroneous in most cases, we do not show results using the type of models given in (25a)–(25c) for this test.

Note in the above computation, we use 100 mesh points and plot the results at time
*t*= 0.12. For the purposes in illustration of the basic solution structure obtained from
using the high-resolution version of the method, we only present results using the sim-
pler “minmod” slope limiter (see [42, 70]) for the runs. Of course, other more sophisticated
limiters, such as “superbee” for example (cf. [2, 70]), can be employed to the methods
for computations also. As far as the global structure of the solution is concerned, we
observe quite similar behavior of the solutions when different limiters are in use to the
method. We further remark that the Courant number (see [42]) in choosing the time step
that maintains stability of the method is 0.9 in the tests. The nonreflecting boundaries
are used on the left and right of the computational domain. Without further notice, we
use the same limiter and the Courant number for all other experiments performed in this
paper.

EXAMPLE4.2. We next are concerned with a two-phase gas-liquid Riemann problem.

*On the left when x*∈ [0, 0.5), we have the gas phase with data
*(ρ, u, p, γ, p*_{∞})*L* = (1.241, 0, 2.753, 1.4, 0),
*and on the right when x*∈ [0.5, 1], we have the liquid phase with data

*(ρ, u, p, γ, p*∞)*R*= (0.991, 0, 3.059 × 10^{−4}, 5.5, 1.505).

We note that the above variables have been nondimensionalized as in the work done by
*Cocchi et al. [12] and Cooke and Chen [17] to simulate underwater explosions in a spheri-*
cally symmetric geometry. We run the problem in a shock tube with 100 mesh points, and
*show the high-resolution results in Fig. 3 at time t* = 0.1 using the γ -based model. From
the figure, we again see the correct behavior of the computed contact discontinuity, and

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 225

**FIG. 3.** *High-resolution results for a two-phase gas-liquid Riemann problem at time t*= 0.1. The solid line
is the exact solution and the points show the computed solution with 100 mesh points.

also the rarefaction and shock waves as in comparison with the exact solution (the solid line shown in the plot). The results of this problem with a source term that accounts for the simplification of a cylindrically symmetric flow will be presented in Section 5 together with two-dimensional results.

EXAMPLE 4.3. Finally, we consider a shock-contact interaction problem studied by
Abgrall [1] and Karni [35] that verifies convergence 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.5 separating two fluids of different equation of states, and a leftward
*going Mach 1.95 shock wave at x* = 0.6 traveling from the right to left. The gas on the left
of the interface is a polytropic gas with

*(ρ, u, p, γ, p*∞)*L* = (1, 0, 1, 1.4, 0),

and the gas on the right of the interface (i.e., on the middle and the preshock state), is a stiffened gas with

*(ρ, u, p, γ, p*_{∞})*M* = (5, 0, 1, 4, 1).

The state behind the shock is

*(ρ, u, p, γ, p*∞)*R* = (7.093, −0.7288, 10, 4, 1).

*The exact solution for this problem in the x-t plane up to time t*= 0.2 is illustrated in Fig. 4
where the density contours are presented.

A snap shot of the computed total internal energy, velocity, and pressure are shown in
*Fig. 5 at time t*= 0.2 where we solve the problem using the high-resolution wave propa-
gation method with 200 mesh points. We can easily see that the shock wave and contact
discontinuity are very well located, and the rarefaction wave moves at the correct speed
with the correct shape. A two-dimensional version of the problem will be considered in
Section 6.

We note that in this section we have only present numerical solutions obtained using theγ -based model to the method. It is interesting to mention that we find little difference between the results as compared to the ones obtained using the volume-fraction model to the method for simulations. Because of this, we omit the presentation of the volume-fraction based numerical results here.