Article ID jcph.1999.6349, available online at http://www.idealibrary.com on

### A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with van der Waals

### Equation of State

Keh-Ming Shyue

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

Received September 17, 1998; revised May 19, 1999

In previous work by the author, a simple interface-capturing approach has been developed and validated for compressible multicomponent flows with a stiffened gas equation of state in multiple space dimensions. The algorithm uses a mixture type of the model equations written in a quasi-conservative form to ensure a consistent approximation of the energy equation near the interfaces where two or more fluid components are present in a grid cell. A standard high-resolution wave propagation method is employed to solve the proposed system, giving an efficient implementation of the algorithm. In this paper, the method is extended to a more general two-phase (liquid–gas) flow where the fluid of interests is characterized by a van der Waals- type equation of state. Several numerical results are presented in both one and two space dimensions that show the feasibility of the method with the Roe solver as applied to practical problems without introducing any spurious oscillations in the pressure near the interfaces. This includes a convergence study of a shock wave in liquid over a gas bubble. To deal with a difficult slip line problem where there is a strong shear flow moving along the interface, we implement the method based on the shock-only Riemann solver with an additional update by the scheme to the total kinetic energy. Rather than using solutions from the basic conservation laws for the density and momenta which incurs large errors, the resulting total kinetic energy is used to the computation of the pressure from the equation of state, yielding typically more accurate results than the unmodified method near the slip lines. This is demonstrated by numerical results of some sample two-dimensional Riemann problems. ° 1999 Academic Pressc

*Key Words: fluid-mixture algorithm; interface capturing method; multicomponent*
flow; slip line problems; stiffened gas equation of state; van der Waals gases.

43

0021-9991/99 $30.00 Copyright c° 1999 by Academic Press All rights of reproduction in any form reserved.

**1. INTRODUCTION**

This paper is concerned with the development of a simple interface-capturing approach for efficient numerical resolution of compressible multicomponent flows with a van der Waals-type fluid (i.e., a fluid with the finite size of a molecule and the nonzero cohesive forces between molecules [21]). We consider a simplified two-phase flow where the fluids of interest consist of two different phases, liquid and gas, separated by immiscible interfaces;

see Fig. 1 for a sample setup in two space dimensions. The algorithm uses a Eulerian formulation of the equations in which, on the gas-phase part of the domain, the fluid is governed by the full set of the compressible Euler equations. In two space dimensions, for instance, it takes the form

*∂**t*

*ρ*
*ρu*
*ρv*
*ρE*

* + ∂*^{x}

*ρu*
*ρu*^{2}*+ p*

*ρuv*
*ρEu + pu*

*+ ∂**y*

*ρv*
*ρuv*
*ρv*^{2}*+ p*
*ρEv + pv*

*= 0,* (1)

where*ρ is the density, u and v are the particle velocities in the x- and y-direction, respec-*
*tively, p is the pressure, and E is the specific total energy. We assume that the gas satisfies*
a van der Waals equation of state,

*p(ρ, e) =*

µ*γ − 1*
1*− bρ*

¶

*(ρe + aρ*^{2}*) − aρ*^{2}*,* (2)

so as to deal with the possible real-gas effect (without phase transitions) when both the
*temperature and pressure are high (cf. [21, 74, 78]). Here e denotes the specific internal*
energy,*γ is the ratio of specific heats (γ > 1), and the quantities a, b are the van der Waals*
gas constants for molecular cohesive forces and the finite size of molecules, respectively
*(a≥ 0, 0 ≤ b < 1/ρ, see [44] for numerical values to various gaseous substances). As usual*
*we set E= e + (u*^{2}*+v*^{2}*)/2. Note that a van der Waals gas of the form (2) reduces to a Noble–*

*Abel gas (also called a constant covolume gas) when a*= 0 [54, 73] and to a polytropic gas
*when b*= 0 as well. The four components of (1) express the conservation of mass, momenta
*in the x- and y-direction, and energy, respectively [15].*

**FIG. 1.** A typical example of a two-phase flow setup with interfaces that separate regions of two different
fluids, liquid and gas, into two parts. Note that the gas component in each elliptical-like shape of the domain may
be different from one another. It is those grid cells that are cut by the interfaces requiring special attention for
proper numerical treatments.

On the liquid-phase part of the domain, however, while the motion of a liquid is assumed to be governed by (1), the algorithm uses the stiffened gas equation of state,

*p(ρ, e) = (γ − 1)ρe − γ B,* (3)

for a fundamental characterization of material properties of the liquids. Here*B is a pressure-*
like constant that, together with*γ , can be determined by a fitting procedure from laboratory*
data (cf. [47]). A typical set of parameter values is for water:*γ = 7, B = 3000 atm [15,*
25], and for human blood: *γ = 5.527, B = 614.6 MPa [52], approximately. It should be*
mentioned that in addition to the modeling of a liquid, Eq. (3) is often used to describe other
type of materials, including many compressible solids of practical importance (cf. [47, 55,
60]).

We want to use a state-of-the-art shock-capturing method on a uniform rectangular grid for the multicomponent flow computations. Clearly when grid cells contain only a single phase of the fluid, there is no problem to solve each phase of the equations separately. But in practice due to the presence and the subsequent dynamic-evolution of the interfaces by the solutions of the governing system, it is inevitable to have two or more fluid components staying within a cell (see Fig. 1). Because of this, the need to give a proper modeling and approximation of these mixture grid cells becomes a principle issue in many of the multicomponent algorithms developed in the literature (cf. [1, 12, 20, 29, 30, 46, 60, 63]).

See [60] in particular for a concise survey of the up-to-date multicomponent methods.

The approach we take here is an extension of the work described by Saurel and Abgrall [60] and by the author [63, 64], in that as opposed to a simplier case with the stiffened gas equation of state (3), a modified van der Waals equation of state (5), see Section 2, is introduced as a basic element to the modeling of the mixing between the stiffened and van der Waals gases. With that, assuming uniform pressure equilibrium and constant particle velocity across the interfaces, from the energy equation, we are able to derive the effective equations for the mixture of the material-dependent quantities near the interfaces. As in the previous work [63], we take these equations to be of the form that do not vary their solutions across the shocks and rarefaction waves as well. Combining the resulting set of effective equations to the Euler equations yields a model system that is written in a quasi- conservative form; see (13) and (36) for the one- and two-dimensional models, respectively.

We use the high-resolution wave propagation method developed by LeVeque [34, 39, 43] to solve the proposed system. Numerical results present in Sections 3.4 and 4.3 show that this is a viable approach in both one and two dimensions as the method is applied with the Roe solver to practical problems without introducing any spurious oscillations in the pressure near the interfaces.

To deal with a difficult slip line problem where there is a strong shear flow moving along the interface, we implement the two-dimensional method based on a shock-only Riemann solver with an additional update by the scheme to the total kinetic energy. Rather than using the solutions from the basic conservation laws for the density and momenta, which incurs large errors, see Section 5.1, the resulting total kinetic energy is employed to compute the pressure from the equation of state, yielding typically more accurate results than the unmodified method near slip lines. This will be discussed further in Section 5.2.

It is true that in many multicomponent flow applications physical effects, such as surface tension and viscosity, play an important role to not only controlling the dynamics of the interfaces, but also influencing the structure of the nearby solutions, and hence need to be

taken in the model for realization (cf. [66]). Typical examples of this kind are the popular water-drop problem in air and the rising air bubble problem in water. In these instances, because the fluids are mostly in a low Mach number regime, compressibility of the fluids is often ignored, and so the problem may be formulated into an incompressible form and solved by a state-of-the-art numerical method for incompressible flows; see [7, 68, 69, 75]

for an example.

Here in contrast to the work just mentioned, we consider a class of problems where the influence of compressibility of the fluids to the solutions is vital, but not the surface tension and viscosity. Examples cover a family of shock wave problems with complicated interface patterns [24, 25, 53, 71, 76, 77] and a hydrodynamic model of sonoluminescence (an acoustic-induced light emission phenomenon) [8, 33, 78]. It is interesting to note that the latter problem in fact motivates the current study of two-phase flows with a van der Waals-type fluid where (2) is used for the modeling of real gases and (3) is employed for a simple approximation of liquids. Our goal here is to establish a basic solution strategy and validate its use via numerical experimentation of some sample problems. Direct simulation of sonoluminescence and other important two-phase flow problems such as shock waves in bubbly liquids or liquid–solid suspensions will be considered in the future; see [65] for a preliminary result of the latter problem.

This paper is organized as follows. In Section 2, we begin by discussing the basic equation of state and the associated thermodynamic stability for the mixing of the stiffened and van der Waals gases within a grid cell. In Section 3, we describe the one-dimensional version of the multicomponent algorithm in more details. This includes the construction of the Riemann problem solution by the shock-only or Roe solver (see Section 3.2) and some numerical results that validate the proposed approach (see Section 3.4). Extension of the basic approach to two space dimensions is explained briefly in Section 4, and some sample two-dimensional results are shown in Section 4.3. A study of how the algorithm works for slip line problems is addressed in Section 5.1, and a correction of the algorithm is made in Section 5.2.

**2. EQUATIONS OF STATE**

To begin, we introduce a hybrid version of the equation of state that is necessary in the
algorithm for modeling the mixing between the stiffened and van der Waals gases within a
grid cell. We do this by taking an approach that expresses (2) and (3) in terms of their natural
*variables: the specific entropy S and the specific volume V= 1/ρ, yielding the formulas*

*p(V, S) = A(S)(V − b)*^{−γ}*− aV*^{−2}*,* (4a)
and

*p(V, S) = A(S)V*^{−γ}*− B,* (4b)

in a respective manner, where*A(S) = R exp[(S − S*ref*)/C**V*] takes the same form in both
*cases (cf. [21, 26]). Here C**V* represents the specific heat at constant volume (see below),*R*
*is the universal gas constant, and S*refis the reference state of the specific entropy. Assuming
the fluids under consideration are all in an adiabatic equilibrium with the same entropy, it
is feasible then to define an extension of the equation of state that combines (4a) with (4b)

in the following way:

*p(V, S) = A(S)(V − b)*^{−γ}*− (B + aV*^{−2}*).* (4c)
With this, it is easy to see that not only the limiting van der Waals gas case (4a) is recovered as
*B → 0, but so also is the stiffened gas case (4b) as both a → 0 and b → 0. In the nonlimiting*
*case, however, where none of the material-dependent quantities a, b, and B is close to zero,*
(4c) does give a way to the representation of the cases in between, that is to the mixing of
the stiffened and van der Waals gases.

It should be noted that, as far as the computational efficiency of the method is concerned,
(4c) is of little use in practice, for it requires some extra work to evaluate the specific entropy
from the solutions of the Euler equations. Instead, by employing the first and second laws of
thermodynamics, we may rewrite (4c) in terms of the often-used variables in gas dynamics,
*ρ and e, as*

*p(ρ, e) =*

µ*γ − 1*
1*− bρ*

¶

*(ρe − B + aρ*^{2}*) − (B + aρ*^{2}*).* (5)

It is clear that (5) is a generalization of the equations of state (2) and (3). As we will see in the latter sections (cf. [64, 65] also for a simpler case), with this so-called modified van der Waals equation of state (5), it is very robust to devise an interface-capturing solver for a family of two-phase flow problems considered here.

*As to the computation of the fluid temperature T (which is important to sonolumines-*
cence, for example), with (5), we may simply use one of the formulas,

*p(V, T ) =* *RT*

*V− b* *− B −* *a*
*V*^{2}*,*
*e(V, T ) =* *RT*

*γ − 1+ BV −* *a*
*V,*

for realization. Note that these two equations are easily derived from (5) using the basic thermodynamic principles; see [21] for more details.

It is important to mention that, in this work, the thermodynamic description of the multi-
component flows is limited by the stability requirement that the internal energy defined in
*(5) be a convex function of its dependent variables V and S. Because of this, the immediate*
consequence is the exclusion of the important but difficult problems involving the transition
of phases here. Analogously in [49, 50] for a Mie–Gr¨uneisen-type equation of state (see
below), it can be shown that explicit conditions for the above mentioned thermodynamic
stability of a van der Waals gas are

*(i) The heat capacity coefficients C**V* *and C**p*must be positive,

*C**V* *= (∂**T**e)|**V* = *R*
*γ − 1> 0,*

*C*_{p}*= (∂**T**h)|**p* = *R[γ RT − 2aρ(1 − bρ)*^{2}]
*(γ − 1)[RT − 2aρ(1 − bρ)*^{2}] *> 0,*
*where h= e + (p/ρ) is the specific enthalpy.*

*(ii) The isentropic bulk modulus K**S*must be positive,

*K**S**= ρ(∂*_{ρ}*p)|**S*=
µ *γ*

1*− bρ*

¶

*(p + B + aρ*^{2}*) − 2aρ*^{2}*> 0.*

(iii) The product of the thermal expansion coefficient*ϑ and the Gr¨uneisen coefficient*
*0 must be nonnegative,*

*ϑ0 =*

"

−1
*ρ(∂**T**ρ)*¯¯

¯¯*p*

#"

1
*ρ(∂**e**p)*¯¯

¯¯*ρ*

#

=

· *R(1 − bρ)*
*RT − 2aρ(1 − bρ)*^{2}

¸·*γ − 1*
1*− bρ*

¸

= *R(γ − 1)*

*RT − 2aρ(1 − bρ)*^{2} *≥ 0.*

In summary, for stability, from the above conditions, this amounts to the satisfaction of the inequalities

*γ > 1, RT > 2aρ(1 − bρ)*^{2}*.*

Combining this with the positiveness of the basic thermodynamic states,*ρ, T , and 0, defines*
the domain of the phase space of this van der Waals gas model. Note in particular that from
*item (ii) given above we have K**S**= ρc*^{2}*, and so the positive of K**S*implies that the speed of
*sound c belongs to a set of real numbers.*

For convenience, we write (5) into a more general Mie–Gr¨uneisen-type equation of state of the form

*p(ρ, e) = 0(ρ)[ρe − ρe**H**(ρ)] − p**H**(ρ),* (6)
where*0 is the Gr¨uneisen coefficient of the material of interests, and the density-dependent*
*functions e**H**, p**H* are the reference Hugoniot states. Typically, the specific form of (6) is
rather complicated in general (cf. [27, 28, 48]). Here for fluids described by (5), we have a
relatively simple but nontrivial model to work with, i.e.,

*0(ρ) =* *γ − 1*

1*− bρ, e**H**(ρ) =* *B − aρ*^{2}

*ρ* *, p**H**(ρ) = B + aρ*^{2}*.* (7)
It is without question that experience gained by studying the current van der Waals gas
model will help to the further development of an efficient multicomponent solver for more
general materials as described by (6).

**3. ONE SPACE DIMENSION**

To motivate the basic idea of our method in multiple space dimensions, we begin by
considering one-dimensional problems with both the gas and the liquid governed by the
*Euler Eqs. (1) of the following form in the x-direction,*

*∂**t*

*ρ*
*ρu*
*ρE*

* + ∂**x*

*ρu*
*ρu*^{2}*+ p*
*ρEu + pu*

* = 0,* (8)

and by the equations of state (2) and (3), respectively. The algorithm uses a popular approach
that employs (8) for the motion of the liquid–gas mixtures of the conserved variables*ρ,*
*ρu, and ρE in a multicomponent grid cell (see [1, 61, 79] for the use of other governing*
equations). We compute the pressure based on the equation of state (6), so long as the
mixture of the problem-dependent material quantities appearing in (7), i.e.,*γ, a, b, and B,*
are defined and known a priori as well. In the algorithm, by following a general procedure
proposed in [63], conditions for those material quantities are found so that not only the
pressure is retained in equilibrium for an interface only fluid-mixture cell, but also the
mixture of the total mass remains conservative on the domain (of course, the mass of each
fluid component may not typically be conserved here). For completeness, we next give a
detailed description of that procedure.

*3.1. Derivation of Model Equations*

As in the previous works [63, 64], our starting point is to consider an interface only
*problem where both the pressure p and particle velocity u are constants in the domain*
*(i.e., p and u satisfy the standard surface-tension free dynamic and kinematic boundary*
conditions on the interface [32], respectively), while the other variables such as the density
*ρ and the material-dependent parameters in the Mie–Gr¨uneisen equation of state (6) are*
having jumps across some interfaces. In this setup, we first write (8) in the following
nonconservative form,

*∂**t**ρ + u∂**x**ρ + ρ∂**x**u* *= 0,*

*∂**t**u+ u∂**x**u*+ 1

*ρ∂**x**p* *= 0,*

*∂**t**(ρe) + ∂**x**(ρeu) + p∂**x**u* *= 0,*

and obtain easily two basic transport equations for the motion of*ρ and ρe as*

*∂**t**ρ + u∂**x**ρ = 0,*

*∂**t**(ρe) + u∂**x**(ρe) = 0.*

By inserting the equation of state (6) into the latter one, we find an alternative description of the energy equation

*∂**t*

µ*p+ p**H*

*0* *+ ρe**H*

¶
*+ u∂**x*

µ*p+ p**H*

*0* *+ ρe**H*

¶

*= 0.* (9)

To see how the pressure would retain in equilibrium as it should be for this model problem, we expand (9) into the form

1

*0(∂**t**p+u∂**x**p)+ p*

·

*∂**t*

µ1
*0*

¶
*+u∂**x*

µ1
*0*

¶¸

+

·

*∂**t*

µ*p**H*

*0* *+ρe**H*

¶
*+u∂**x*

µ*p**H*

*0* *+ρe**H*

¶¸

*= 0.*

Applying the assumed state of the pressure equilibrium to the above equation comes to a simplier expression as

*p*

·

*∂**t*

µ1
*0*

¶
*+ u∂**x*

µ1
*0*

¶¸

+

·

*∂**t*

µ*p**H*

*0* *+ ρe**H*

¶
*+ u∂**x*

µ*p**H*

*0* *+ ρe**H*

¶¸

*= 0.*

*Now since this equation should hold for any p in the physical space, it implies right away*
that the terms in the square bracket of the equation should be vanished simultaneously,
yielding a system of the following two equations:

*∂**t*

µ1
*0*

¶
*+ u∂**x*

µ1
*0*

¶

*= 0,* (10a)

*∂**t*

µ*p**H*

*0* *+ ρe**H*

¶
*+ u∂**x*

µ*p**H*

*0* *+ ρe**H*

¶

*= 0.* (10b)

It is important to note that in order to have the correct pressure equilibrium in (9) near the
interfaces, (10a) and (10b) are the two key equations that should be satisfied for any given
expression of*0, p**H**, and e**H*in the Mie–Gr¨uneisen equation of state (6). From them, for a
class of real materials modeled analytically by (6) at least, we may continue to work out
suitable conditions for the further details of the related material parameters.

Consider the current two-phase flow application with*0, p**H**, and e**H*defined by (7), for
example. Equations (10a) and (10b) now have the form

*∂**t*

µ1*− bρ*
*γ − 1*

¶
*+ u∂**x*

µ1*− bρ*
*γ − 1*

¶

*= 0, (11a)*

*∂**t*

µ*γ − bρ*

*γ − 1B +* 2*− γ − bρ*
*γ − 1* *aρ*^{2}

¶
*+ u∂**x*

µ*γ − bρ*

*γ − 1* *B +*2*− γ − bρ*
*γ − 1* *aρ*^{2}

¶

*= 0, (11b)*

in a respective manner. It is obvious that, in addition to (11a) and (11b), we need to impose two supplementary conditions so as to have enough equations for the four material quantities:

*γ , a, b, and B. In our approach, this is done quite easily by simply splitting (11a) into the*
following two parts,

*∂**t*

µ 1
*γ − 1*

¶
*+ u∂**x*

µ 1
*γ − 1*

¶

*= 0,* (11c)

*∂**t*

µ *bρ*
*γ − 1*

¶
*+ u∂**x*

µ *bρ*
*γ − 1*

¶

*= 0,* (11d)

and also (11b) into the terms

*∂**t*

µ*γ − bρ*
*γ − 1* *B*

¶
*+ u∂**x*

µ*γ − bρ*
*γ − 1* *B*

¶

*= 0,* (11e)

*∂**t*

µ2*− γ − bρ*
*γ − 1* *aρ*^{2}

¶
*+ u∂**x*

µ2*− γ − bρ*
*γ − 1* *aρ*^{2}

¶

*= 0.* *(11f )*

Having done so, we arrive at a primitive form of the transport equations (11c)–(11f) for the
variables 1*/(γ − 1), bρ/(γ − 1), B(γ − bρ)/(γ − 1), and aρ*^{2}*(2 − γ − bρ)/(γ − 1). We*
note that with them for this interface only problem it is sufficient to have all the quantities
*γ , a, b, and B determined at all times, provided that the initial condition has been properly*
set for the computation; see Section 3.1.1.

Up to this point, our discussion stresses only on an approach that is capable of maintaining the pressure in equilibrium for a model interface only problem. Since in practice we are

interested in shock wave problems as well, we should thus take the equations, i.e., (11c)–

(11f ), in a form such that*γ , a, b, and B remain unchanged across both shocks and rarefaction*
waves. Concerning this, it is easy to see that with*γ governed by (11c) there is no problem*
*to do so (cf. [1, 63]). For b andB, however, due to the appearance of the linear factor of ρ*
in (11d) and (11e), it turns out that, in a time when such a scenario occurs, for consistency
these two equations should be modified in such a way that each of them reduces to the basic
mass conservation law of the fluid mixture. The derivation of the modification is simple.

Without going into the details, we write down the corrected version of the corresponding equation as follows:

*∂**t*

µ *bρ*
*γ − 1*

¶
*+ ∂**x*

µ *bρ*
*γ − 1u*

¶

*= 0,* (12a)

*∂**t*

µ*γ − bρ*
*γ − 1* *B*

¶
*+ ∂**x*

µ*γ − bρ*
*γ − 1* *Bu*

¶

=
µ *γ*

*γ − 1B*

¶

*∂**x**u.* (12b)

*We next come to the discussion for the quantity a described by (11f). In this situation, if*
we assume the proper smoothness of the nonlinear factors of*ρ in the equation (such as in*
*the case of rarefaction waves), as in (12a) and (12b) for the quantities b andB, respectively,*
we may derive an equation of the form

*∂**t*

µ2*− γ − bρ*
*γ − 1* *aρ*^{2}

¶
*+ ∂**x*

µ2*− γ − bρ*
*γ − 1* *aρ*^{2}*u*

¶

= −

µ2*− γ − 2bρ*
*γ − 1* *aρ*^{2}

¶

*∂**x**u* (12c)

*that admits the desired constancy of a and the mass conservation as well. It is important to*
note that this is not the case, however, if the smoothness assumption on*ρ is violated; this*
happens in the event of shock waves, for there is no way to differentiate the discontinuous
terms*ρ*^{2}and*ρ*^{3}yielding the partial derivatives on*ρ itself. Despite its apparent difficulty*
*to fulfill the requirement of a across both shocks and rarefaction waves, (12c) is the correct*
*form to be used in the model so that the fluid mixture aρ*^{2}*(2 − γ − bρ)/(γ − 1) can be*
solved and applied readily to the computation of the pressure; see (14) below.

Here, rather than using (12c), we introduce a simplier linear advection equation of the form

*∂**t**a+ u∂**x**a*= 0 (12d)

*for the motion of the mixture of a. Clearly, (12d) gives an accurate description of a in all*
the solution regions discussed above. In practice, it is a good one to use in the algorithm for
numerical approximation, see Section 3.4 for some sample results.

To sum up, as in [63], we use the term effective equations to describe the set of equations that govern the motion of the material-dependent mixtures of the problems. In the current case, there consists of a couple sytem of five equations, i.e., (11c) and (12a)–(12d). It should be noted that when further taking the numerical aspect of the model approximation into consideration, with the full Euler equations, these effective equations are the proper ones to use for practical problems, but not the other form of the equations; see Section 3.3 and [60, 63] for the details.

Putting all the things together, with the modified van der Waals equation of state (5), the model equations we propose to solve one-dimensional multicomponent pro- blems are

*∂**t**ρ + ∂**x**(ρu) = 0*

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

*∂**t**(ρE) + ∂**x**(ρEu + pu) = 0*

*∂**t*

¡ _{b}_{ρ}

*γ − 1*

¢*+ ∂**x*

¡ _{b}_{ρ}

*γ − 1**u*¢

= 0

*∂**t*

¡_{γ − bρ}

*γ − 1**B*¢
*+ ∂**x*

¡_{γ − bρ}

*γ − 1**Bu*¢

=¡ _{γ B}

*γ − 1*

¢*∂**x**u*

*∂**t*

¡_{2}_{− γ − bρ}

*γ − 1* *aρ*^{2}¢
*+ ∂**x*

¡_{2}_{− γ − bρ}

*γ − 1* *aρ*^{2}*u*¢

= −¡_{2}_{− γ − 2bρ}

*γ − 1* *aρ*^{2}¢

*∂**x**u*

*∂**t*

¡ _{1}

*γ − 1*

¢*+ u∂**x*

¡ _{1}

*γ − 1*

¢= 0

*∂**t**a+ u∂**x**a*= 0;

(13)

this gives us eight equations in total to be solved in one space dimension that is nicely
independent of the number of fluid components involved in the problem. It is clear that
in this system the first three are the Euler equations which are used to make certain the
conservation of the basic fluid mixtures,*ρ, ρu, and ρE, while the remaining are the effective*
equations that are introduced to ensure the correct mixing of the problem-dependent material
variables near the interfaces. With a system expressed in this way, there is no problem to
compute all the state variables of interest, including the pressure from the equation of
state

*p*=

·

*ρE −(ρu)*^{2}
2*ρ* −

µ*γ − bρ*
*γ − 1B*

¶

−

µ2*− γ − bρ*
*γ − 1* *aρ*^{2}

¶¸Áµ 1

*γ − 1*− *bρ*
*γ − 1*

¶
*. (14)*

Note that (13) reduces to a well-tested*γ -based model for a stiffened gas when a = b = 0*
(cf. [60, 63, 64]) and for a polytropic gas when*B = 0 as well (cf. [1]). As before, the proposed*
system (13) is not written in the full conservation form, but is rather a quasi-conservative
system of equations. In addition, the nonzero terms on the right-hand side of the fifth and
sixth equations of (13) should be viewed as an integrated part of the whole system, but not
be considered as a source term. This fact can be realized easily by formulating the model as
a quasi-linear system of equations (of course, with the assumption of a proper smoothness
of the solutions), namely,

*∂**t**q+ A(q)∂**x**q*= 0 (15a)

with the state vector

*q*=

·

*ρ, ρu, ρE,* *bρ*

*γ − 1,γ − bρ*

*γ − 1* *B,*2*− γ − bρ*

*γ − 1* *aρ*^{2}*,* 1
*γ − 1, a*

¸*T*

(15b)

and the matrix

*A(q) =*

0 1 0 0 0 0 0 0

*K− u*^{2} *u(2 − 0)* *0* *p0* *−0* *−0* *−p0 0*

*u(K − H) H − u*^{2}*0 u(0 + 1) up0 −u0 −u0 −up0 0*

*−ϕu* *ϕ* 0 *u* 0 0 0 0

*ϕuB* *−ϕB* 0 0 *u* 0 0 0

*−χu* *χ* 0 0 0 *u* 0 0

0 0 0 0 0 0 *u* 0

0 0 0 0 0 0 0 *u*

*. (15c)*

*Here K= 0u*^{2}*/2, H = E + (p/ρ), ϕ = b/(γ − 1), and χ = aρ(4 − 2γ − 3bρ)/(γ − 1).*

Recall that*0 is the Gr¨uneisen coefficient, see (7), which gives different values for different*
fluids.

*It is easy to show that for each physically relevant value of the state variables q defined*
*in the region of thermodynamic stability (see Section 2), the eigenstructure of the matrix A*
is possessed of real eigenvalues

*3 = diag(λ*1*, λ*2*, . . . , λ*8*) = diag(u − c, u, u + c, u, . . . , u)* (16a)
and a complete set of eigenvectors of the form

*R= (r*1*, r*2*, . . . , r*8*) =*

1 1 1 0 0 0 0 0

*u− c* *u* *u+ c* 0 0 0 0 0

*H− uc* ^{1}_{2}*u*^{2} *H+ uc −p 1 1 p 0*

*ϕ* 0 *ϕ* 1 0 0 0 0

*ϕB* 0 *ϕB* 0 1 0 0 0

*χ* 0 *χ* 0 0 1 0 0

0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 1

(16b)

*with Ar**k**= λ**k**r**k*. Thus, (15a) is a hyperbolic system of partial differential equations and so
is our multicomponent model (13). Regarding discontinuous solutions of the system, such
as shock waves or contact discontinuities, (13) has the usual form of the Rankine–Hugoniot
jump conditions across the waves; see Section 3.2.1 for more details.

With these comments, it should be sensible to use the proposed model for practical computations. The numerical method to be discussed in Section 3.3 is a consistent approxi- mation of the model that gives excellent results for a wide variety of problems as illustrated in Section 3.4.

*3.1.1. Initialization of fluid-mixture cells. Consider a typical multicomponent setting*
*in which there are m different fluids 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 assumption, we have Y*^{(i)}*∈ [0, 1] and*P*m*

*i*=1*Y** ^{(i)}*= 1. Suppose that for each

*component i the state variables such asρ*

^{(i)}*, u*

^{(i)}*, p*

^{(i)}*, γ*

^{(i)}*, a*

^{(i)}*, b*

*, and*

^{(i)}*B*

*are known a priori. The objective is to give a proper definition of the fluid mixtures so that they can be used as an initial condition for our model equations (13) to the computations.*

^{(i)}In the algorithm, we follow a common practice by evaluating the mixing states*ρ, ρu, and*
*ρe as a volume-weighted sum over the set of components ρ*^{(i)}*, ρ*^{(i)}*u** ^{(i)}*, and

*ρ*

^{(i)}*e*

*given above,*

^{(i)}

*ρ*
*ρu*
*ρe*

=X^{m}

*i*=1

*Y*^{(i)}

*ρ*^{(i)}*ρ*^{(i)}*u*^{(i)}*ρ*^{(i)}*e*^{(i)}

* .* (17)

With this result, we then compute the mixture of the total energy by*ρE = ρe + (ρu)*^{2}*/(2ρ);*

this completes the definition of the conservative variables for the first three equations of the model.

To find the initial fluid mixtures, 1*/(γ − 1), bρ/(γ − 1), B(γ − bρ)/(γ − 1), and*
*aρ*^{2}*(2 − γ − bρ)/(γ − 1), for the next four equations, we use the equation of state (5),*
where written as a function of the volume fraction it reads

µ1*− bρ*
*γ − 1*

¶

*(p + B + aρ*^{2}*) + B − aρ*^{2}*= ρe =*
X*m*

*i*=1

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

=
X*m*

*i*=1

*Y*^{(i)}

·µ1*− b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

¶

*(p*^{(i)}*+ B*^{(i)}*+ a*^{(i)}*(ρ*^{(i)}*)*^{2}*) + B*^{(i)}*− a*^{(i)}*(ρ*^{(i)}*)*^{2}

¸
*.*

By taking a similar approach as employed in Section 3.1 for the derivation of the effective equations, it comes out easily a splitting of the above equation into the form

*γ − 1*1
*b**ρ*
*γ − 1*
*γ − bρ*

*γ − 1**B*

2*− γ − bρ*
*γ − 1* *aρ*^{2}

=
X*m*

*i*=1

*Y*^{(i)}

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

*b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

*γ*

^{(i)}*− b*

^{(i)}*ρ*

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

*B*

^{(i)}2*− γ*^{(i)}*− b*^{(i)}*ρ*^{(i)}

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

*a*

^{(i)}*(ρ*

^{(i)}*)*

^{2}

*,* (18)

*where in the process of splitting the terms the pressure p is chosen to satisfy the relation as*
µ1*− bρ*

*γ − 1*

¶
*p*=

X*m*
*i*=1

*Y*^{(i)}

·µ1*− b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

¶
*p*^{(i)}

¸

(19)

(cf. [63] also for a simplier case of how this is done). With (18), it is easy to see that when
each of the partial pressures is in equilibrium within a grid cell, the pressure p acquired
*from (19) would remain in equilibrium also; i.e., p= p*^{(i)}*, for i= 1, 2, . . . , m. Furthermore,*
from (18), we are able to obtain an explicit expression of the material-dependent parameters
*γ, a, b, and B in terms of the volume-fraction function and the original set of data. The*
results are

*γ = 1 + 1*
,Ã _{m}

X

*i*=1

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

!

*,* (20a)

*b*=

" * _{m}*
X

*i*=1

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

*b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

!#,"Ã * _{m}*
X

*i*=1

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

!Ã * _{m}*
X

*i*=1

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

!#

*,* (20b)

*B =*

" * _{m}*
X

*i*=1

*Y*^{(i)}

Ã*γ*^{(i)}*− b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

*B*

^{(i)}!#,"

1+
X*m*

*i*=1

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

1*− b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

!#

*,* (20c)

*a*=

" * _{m}*
X

*i*=1

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

2*− γ*^{(i)}*− b*^{(i)}*ρ*^{(i)}

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

*a*

^{(i)}*(b*

^{(i)}*)*

^{2}

!# ,( Ã * _{m}*
X

*i*=1

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

!2

×

"

−1 +
X*m*

*i*=1

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

1*− b*^{(i)}*ρ*^{(i)}*γ** ^{(i)}*− 1

!#)

*.* (20d)

We note that as in the continuous counterpart (12c), (20d) is not used in practice for setting
*the initial mixture of a, but is done by the formula a*=P*m*

*i*=1*Y*^{(i)}*a** ^{(i)}*instead.

To end this subsection, it should be mentioned that by following the same approach introduced in [63] we may reformulate the model (13) into its variant form, the so-called volume-fraction model, that is robust when both the set of governing equations and the type of equations of state are different from one fluid component to the others, separated by the interfaces. To keep the presentation simple and clear, we omit the discussion of that model here, but refer the reader to [64, 65] for the details.

*3.2. Approximate Riemann Solvers*

Before describing numerical methods to solve (13), we pause to discuss the construction
of the Riemann problem solutions which is one of the major steps in our multicomponent
algorithm. For comparison purposes, we present two popular approaches for the resolution
*of the Riemann problem with piecewise constant data q*_{L}*and q** _{R}* to the left and right of the
interface.

*3.2.1. Shock-only solver. First, we are concerned with a shock-only approximation of*
the Riemann solver that ignores the possibility of rarefaction waves and simply construct a
solution in which each pair of the states is connected along the Hugoniot locus for a shock
*(cf. [3, 10, 13]). In this approach, the key step is to find the midstate (u*_{∗}*, p*_{∗}*) in the u– p*
*phase plane so that it can connect to (u*_{L}*, p**L*) by a 1-shock and to*(u**R**, p**R**) by a 3-shock.*

It is well known that this is equivalent to solving the following nonlinear equation in an
*iterative manner for the pressure p*_{∗}:

*h(p*∗*) = u**∗R**(p*∗*) − u**∗L**(p*∗*) = 0.* (21)
*Here u*_{∗L}*and u** _{∗R}*are the velocities defined by connecting the states along the 1-shock and
3-shock curves, respectively,

*u*_{∗L}*(p) = u**L*− *p− p**L*

*M*_{L}*(p), u**∗R**(p) = u**R*+ *p− p**R*

*M*_{R}*(p),*

*with M** _{ι}* denoting the Lagrangian shock speed, for

*ι = L or R. In the current application*

*with the modified van der Waals equation of state (5), when a*= 0 (i.e., the vanishing of the

*aρ*

^{2}

*term), we may compute M*

*directly by evaluating the formula*

_{ι}*M*_{ι}^{2}*(p) = C*_{ι}^{2}

· 1+

µ*γ**ι*+ 1
2*γ*_{ι}

¶µ*p+ B**ι*

*p*_{ι}*+ B** _{ι}* − 1

¶ ¸
*,*

*where C*_{ι}*= ρ**ι**c** _{ι}*is the Lagrangian sound speed (cf. [73]). Note that this is as a result derived
from the Rankine–Hugoniot jump condition across the shock waves,

*M*_{ι}^{2}*(e*_{∗ι}*− e*_{ι}*) =*¡

*p*^{2}_{∗}*− p*^{2}* _{ι}*¢±

2*,* (22)

*with e*_{∗ι}*= e(p*_{∗}*, ρ*_{∗ι}*) a function of the midstate density ρ*_{∗ι}*, a quantity that is related to p, p** _{ι}*,

*and M*

*in the following way*

_{ι}*ρ*_{∗ι}*(p) =*

·

*ρ*_{ι}^{−1}− *p− p**ι*

*M*_{ι}^{2}*(p)*

¸_{−1}
*,*

*while in the more general case when a*6= 0, there is not such a close form solution available
*for W*_{ι}*. Instead, we need to solve (22) iteratively for M** _{ι}*; this is a typical thing to do when
using the method for real gases (cf. [13]).

When applying a standard root-finding approach such as the secant method to (21), we have a 2-step iteration scheme as follows,

*p*^{(n+1)}_{∗} *= p*_{∗}* ^{(n)}*− ¯¯

*p*

_{∗}

^{(n)}*− p*

^{(n−1)}_{∗}¯¯

¯¯*u*^{(n)}_{∗L}*− u*^{(n−1)}* _{∗L}* ¯¯ + ¯¯

*u*

^{(n)}

_{∗R}*− u*

^{(n−1)}*¯¯£*

_{∗R}*u*^{(n)}_{∗R}*− u*^{(n)}* _{∗L}*¤

*,* (23)

*where u*^{(n)}_{∗ι}*= u**∗ι**[ p*^{(n)}_{∗} ]*, for ι = L or R, and n = 1, 2, . . . (until convergence). With a suitable*
*choice of the starting values p*_{∗}^{(0)}*and p*_{∗}* ^{(1)}*, method (23) typically converges to the exact

*solution p*

_{∗}

*at a superlinear rate [31]. For gas dynamics, it is a common practice to set p*

_{∗}

^{(0)}*and p*

_{∗}

*by*

^{(1)}*p*_{∗}* ^{(0)}*=

*p*

*R*

*C*

*L*

*+ p*

*L*

*C*

*R*

*− (u*

*R*

*− u*

*L*

*)C*

*L*

*C*

*R*

*C**L**+ C**R*

*,*
*u*^{(0)}_{∗L}*= u**L*− *p*^{(0)}_{∗} *− p**L*

*C**L*

*, u*^{(0)}_{∗R}*= u**R*+ *p*^{(0)}_{∗} *− p**R*

*C**R*

*,* (24)

*p*_{∗}* ^{(1)}*=

*p*

_{R}*M*

^{(0)}

_{L}*+ p*

*L*

*M*

^{(0)}

_{R}*− (u*

*R*

*− u*

*L*

*)M*

^{(0)}*L*

*M*

_{R}

^{(0)}*M*

^{(0)}

_{L}*+ M*

^{(0)}*R*

*,*

*with M*_{ι}^{(0)}*= M*_{ι}*[ p*_{∗}* ^{(0)}*]. After a satisfactory convergence of the scheme, we may then calculate

*u*

_{∗}by

*u*_{∗}= *p**L**− p**R**+ u**L**M**L**(p*∗*) + u**R**M**R**(p*∗*)*
*M**L**(p*∗*) + M**R**(p*∗*)* *.*

We note that alternatively we may use a 1-step Newton method for the solution of (21).

Since the derivation of the scheme is more involved due to the need to compute the derivative
*term d p*_{∗}*/du*∗, we do not discuss the method here, but refer the reader to [13] for more
details.

Figure 2 shows a typical solution structure of the Riemann problem considered here.

Clearly, in a shock-only approximate solver, we replace the leftward-going rarefaction wave by a 1-shock and so the solution consists of three discontinuities moving at constant speeds. Here the propagation speed of each discontinuity is determined by

*λ*1*= u*_{∗}− *M**L**(p*∗*)*

*ρ**∗L**(p*∗*), λ*2*= u*_{∗}*, λ*3*= u*_{∗}+ *M**R**(p*∗*)*

*ρ**∗R**(p*∗*),* (25a)

**FIG. 2.** Typical solution structure of the Riemann problem for our multicomponent model (13). Note that in
a shock-only approximation of the approximate Riemann solver, the rarefaction wave is replaced by an entropy-
violating shock.

with the jumps across each of them computed by the difference between the states to the left and right of the discontinuity,

*W*1 *= q*_{∗L}*− q**L**, W*2*= q*_{∗R}*− q*_{∗L}*, W*3*= q**R**− q*_{∗R}*.* (25b)
Wave propagation methods are based on using these propagating discontinuities to update
the cell averages in the cells neighboring each interface.

It is true that no matter what iterative method is employed to the solution of (21), the approach is quite expansive as compared to the approximate solver of Roe described next.

Nevertheless, there are various situations where this approach is worthwhile and can provide more accurate results than the Roe solver does. This includes some examples shown in [18, 41, 56] and many difficult problems with strong shock waves and stiff equations of state. Moreover, it is a straightforward matter to generalize the approach that covers the case with surface tension effect across the interfaces; results on this aspect will be reported elsewhere.

*3.2.2. Roe solver. In a Roe’s approximate Riemann solver, we replace the nonlinear*
*system (13) with data q**L**and q**R*by a linear system of the form

*∂**t**q+ ˆA(q**L**, q**R**)∂**x**q* *= 0.* (26)
Here ˆ*A(q**L**, q**R**) is a constant matrix that depends on the initial data and is a local linearization*
*of the matrix A in (15c) about an average state. To find that matrix, as it is often done in*
many other Roe solvers (cf. [9, 22, 23]), we want to seek an average state such that the
difference of the fluxes in the conservation part of (13) (i.e., the first four equations of the
system) are equal to the respective first order approximations 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]*^{(i)}*,* (27)
*for i* *= 1, 2, 3, 4, where F ∈ R*^{4}is the usual definition of the fluxes for conservation laws,
and*1F*^{(i)}*is the i th component of1F. With that, it is a straightforward matter to obtain*
*the results for ˆu, ˆH, and ˆϕ by the standard “Roe-averaging” approach; i.e., for a given pair*
*(ρ**L**, ρ**R**), the average state for a quantity z is defined by*

*ˆz*=*√ρ**L**z**L**+ √ρ**R**z**R*

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

*.* (28)

Note that in the process of the derivation, as in [63] we have chosen the averages*(d*1*/0)*
and*(dp/0) based on (28) so that the expression*

*1p =*· dµ
1
*0*

¶
*1*

µ*p*
*0*

¶

−µd
*p*
*0*

¶
*1*

µ1
*0*

¶¸Á dµ
1
*0*

¶^{2}

*is satisfied approximately. With that we set ˆp*= d*(p/0)/ d(1/0) and ˆ0 = 1/ d(1/0). To finish*
the construction of ˆ*A(q**L**, q**R**), we still need to find the averages of B and χ. Since there*
is no unique way to do so, we might as well compute ˆ*B and ˆa using (28) and set ˆχ =*
[*(2/ ˆ0) − 2 − ˆϕ ˆρ]ˆa ˆρ, where ˆρ = √ρ**L**ρ**R*. Numerical results shown in Section 3.4 indicate
that the set of average states described here is a good one to use for practical multicomponent
problems.

The solution of the linear problem (26) consists of eight discontinuities propagating at
constant speeds (for a system of eight equations). 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*=
X8
*k*=1

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

*where ˆr**k**is the kth eigenvector of ˆA; see (16a) and (16b). The scalar ˆα**k*gives the strength
across the discontinuity that can be determined easily from (29). We find

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

*0*ˆ ˆ

*c*

^{2}

·

−*ˆu*^{2}

2 *1q*^{(1)}*+ ˆu1q*^{(2)}*− 1q*^{(3)}*+ ˆp*¡

*1q*^{(7)}*− 1q** ^{(4)}*¢

*+ 1q*^{(5)}*+ 1q*^{(6)}

¸
*,*
*α*ˆ3= 1

2 ˆ*c*

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

¤*,*

*α*ˆ1*= 1q*^{(1)}*− ˆα*2*− ˆα*3*,* *α*ˆ4*= 1q*^{(4)}*− ˆϕ(1ρ − ˆα*2*),* (30)
*α*ˆ5*= 1q*^{(5)}*− ˆϕ ˆB(1ρ − ˆα*2*), ˆα*6*= 1q*^{(6)}*− ˆχ(1ρ − ˆα*2*),*

*α*ˆ7*= 1q*^{(7)}*,* *α*ˆ8*= 1q*^{(8)}*,*
where ˆ*c*=p

*0[ ˆH − (ˆu*ˆ ^{2}*/2) + ˆp ˆϕ − ˆϕ ˆB − ˆχ] is the speed of sound.*

Note that in this Riemann solution, except the discontinuities for ˆ*λ*1*= ˆu − ˆc, and*
ˆ*λ*3*= ˆu + ˆc, all the other discontinuities (six of them) are propagating at the same speed*
*ˆu. For practical purposes, we may view these discontinuities as a single one with the*
operator *W*2 defined by combining all the jumps across the ˆ*λ*2 wave family; i.e., set
*W*2*= ˆα*2*ˆr*2+P8

*k*=4*α*ˆ*k**ˆr**k*. Clearly, doing so removes the effect of the wave families ˆ*λ*4

to ˆ*λ*8from the solution. With this notation, we also write*W**k**= ˆα**k**ˆr**k*to represent the jump
*across the k-wave for k*= 1 or 3.

*3.3. Wave Propagation Methods*

We use the high-resolution method based on a wave-propagation viewpoint to compute approximate solutions of our multicomponent model introduced in Section 3.1. The method is a variant of the fluctuation-and-signal scheme of Roe [58, 59] in that we solve the Riemann

problems at each cell interface and use the resulting waves (i.e., discontinuities moving at constant speeds) to update the solutions in neighboring grid cells (cf. [34, 39]).

For simlicity, 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 (cf.

*[35, 40]). We use a standard finite-volume formulation in which the value Q*^{n}* _{j}*approximates

*the cell average of the solution over the grid cell [x*

*j*

*, x*

*j*+1

*] at time t*

*n*:

*Q*^{n}* _{j}* ≈ 1

*1x*

Z *x*_{j+1}*x*_{j}

*q(x, t**n**) dx.*

*The time step from the current time t**n**to the next t**n*+1is denoted by*1t.*

*3.3.1. First order method. 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

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

*m*_{w}

X

*k*=1

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

where*λ**k**∈ R and W**k*∈ R^{m}*are solutions of the kth wave family, for k= 1, 2, . . . , m** _{w}*, ob-

*tained from solving the Riemann problems at cell interfaces x*

*j*

*and x*

*j*+1

*, λ*

^{−}

*= min(λ, 0),*and

*λ*

^{+}

*= max(λ, 0). Clearly, the method belongs to a class of upwind schemes (cf.*

[23, 36]), and it will be shown next that the method is quasi-conservative in the sense that when applying the method to (13) not only the conservation equations but also the transport equations are approximated in a consistent manner by the method with the chosen Riemann solver.

To demonstrate that, we first analyze our method for an interface only problem as de-
scribed in Section 3.1. Without loss of generality, we consider a single Riemann problem
*where at cell interface x*_{j}*the initial data consists of uniform pressure p*_{0}and constant par-
*ticle velocity u*0 to the left and right of the interface, but with jumps on the other state
*variables of q. Assuming a positive velocity u*0 *> 0, for example. If the problem is solved*
*by using a shock-only Riemann solver (see Section 3.2.1), from (31) the cell average Q*^{n}* _{j}*
would be updated by

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

*1x(λ*2*W*2*)*^{n}*j**,* (32a)

or equivalently by

*ρ*
*ρu*
*ρE*

*b**ρ*
*γ − 1*
*γ − bρ*

*γ − 1**B*

2*− γ − bρ*
*γ − 1* *aρ*^{2}

1
*γ − 1*

*a*

*n*+1

*j*

=

*ρ*
*ρu*0

*ρE*

*b**ρ*
*γ − 1*
*γ − bρ*

*γ − 1**B*

2*− γ − bρ*
*γ − 1* *aρ*^{2}

1
*γ − 1*

*a*

*n*

*j*

− *1t*
*1xu*_{0}

*1ρ*
*u*0*1ρ*
*1ρE*
*1*_{γ − 1}^{b}^{ρ}*1*^{γ − bρ}_{γ − 1}*B*
*1*^{2}^{− γ − bρ}_{γ − 1}*aρ*^{2}

*1*_{γ − 1}^{1}
*1a*

*n*

*j*

*,* (32b)

when expressing (32a) in terms of the solution states of the problem. Noting that in this
case the difference operator *1 is simply applied to the Riemann data Q*^{n}*j*−1 *and Q*^{n}* _{j}* on