doi:10.1006/jcph.2001.6801, available online at http://www.idealibrary.com on

### A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with Mie–Gr ¨uneisen 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 October 10, 2000; revised April 6, 2001

A simple interface-capturing approach proposed previously by the author for effi-
cient numerical resolution of multicomponent problems with a van der Waals fluid
*[J. Comput. Phys., 156 (1999), pp. 43–88] is extended to a more general case with*
real materials characterized by a Mie–Gr¨uneisen equation of state. As before, the
flow regime of interests is assumed to be homogeneous with no jumps in the pres-
sure and velocity (the normal component of it) across the interfaces that separate
two regions of different fluid components. The algorithm uses a mixture type of
the model system that is formed by combining the Euler equations of gas dyna-
mics for the basic conserved variables and an additional set of effective equations
for the problem-dependent material quantities. In this approach, the latter equations
are introduced in the algorithm primarily for an easy computation of the pressure
from the equation of state, and are derived so as to ensure a consistent modeling
of the energy equation near the interfaces where two or more fluid components are
present in a grid cell, and also the fulfillment of the mass equation in the other single
component regions. A standard high-resolution wave propagation method designed
originally for single component flows is generalized to solve the proposed system for
multicomponent flows, giving an efficient implementation of the algorithm. Several
numerical results are presented in both one and two space dimensions that show
the feasibility of the method with the Roe Riemann solver as applied to a reason-
able class of practical problems without introducing any spurious oscillations in the
pressure near the interfaces. This includes results obtained using a multicomponent
version of the AMRCLAW software package of Berger and LeVeque for the simu-
lation of the impact of an underwater aluminum plate to a copper plate in two space
dimensions. ° 2001 Academic Pressc

*Key Words: Godunov-type scheme; impact problems; Mie–Gr¨uneisen equation of*
state; multicomponent flows; Roe approximate Riemann solver.

678 0021-9991/01 $35.00

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

**1. INTRODUCTION**

In this paper, we describe extensions of a fluid-mixture type algorithm proposed previ- ously by the author for efficient numerical resolution of multicomponent problems with a van der Waals gas (cf. [45]) to a more general case with materials characterized by a Mie–Gr¨uneisen equation of state of the form

*p(ρ, e) = p*ref*(ρ) + 0(ρ)ρ [e − e*ref*(ρ)].* (1)
*Here p, ρ, and e denote the pressure, density, and specific internal energy of the flow,*
respectively;*0 = (1/ρ)(∂p/∂e)|**ρ**is the Gr¨uneisen coefficient, and p*ref*, e*refare the properly
chosen states of the pressure and internal energy along some reference curve (e.g., along
an isentrope, a single shock Hugoniot, or the other empirically fitting curves) in order to
match the experimental data of the material being examined. Note that, for simplicity, each
of the expressions*0, p*ref*, and e*refis taken as a function of the density only. Even with this
simplification, the analytical form of the equation of state (1) is an adequate approximation
to a wide variety of materials of interest. This includes some gaseous or solid explosives
and solid metals under high pressure; see Section 2 for the details.

It is known that for a general multicomponent flow system (compressible or not), de- pending specifically on conditions such as the topological structure of the interfaces and jumps of fluid properties across them, one can distinguish various type of flow regimes of practical importance, e.g., annular flow, slug flow, bubbly flow, and so on (cf. [10, 49, 51, 53]). Among them, in this work (cf. [44, 45, 46] also), we are interested in problems arising from a so-called homogeneous flow in which there is typically a strong coupling between the motion of each fluid component, and assumes a simple flow condition with no jumps in the pressure and velocity (the normal component of it) across interfaces that separate two different fluid components. Consider a one-dimensional inviscid compressible flow, for example. The basic conservation laws for the fluid mixtures of mass, momentum, and energy are

*∂ρ*

*∂t* + *∂*

*∂x(ρu) = 0,*

*∂*

*∂t(ρu) +* *∂*

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

*∂*

*∂t(ρE) +* *∂*

*∂x(ρEu + pu) = 0,*

*respectively, where u is the particle velocity, and E= e + u*^{2}*/2 is the specific total energy.*

Clearly, (2) takes the same form as the standard Euler equations of gas dynamics for a single component flow, and has been used quite extensively in modeling the behavior of a homogeneous flow (cf. [44, 45, 46]). Note that, in contrast to the case mentioned above, the use of a separate set of equations for each fluid component is often preferred for nonhomogeneous multicomponent problems; see [2, 12, 39] for an example.

To solve a compressible multicomponent problem with a general. Mie–Gr¨uneisen equa- tion of state (1), we want to use an Eulerian formulation of the equations as in the form described in (2), and to employ a state-of-the-art shock capturing method on a uniform rect- angular grid for numerical approximation. Aside from the basic properties that a numerical method should follow in regions where the solutions contain only a single component

(cf. [7]), one major problem in the method development of a multicomponent solver is the
need to devise a proper model and treatment of the numerical mixing between more than
one fluid component within a grid cell. For the homogeneous flow problems considered
here, in particular, it is imperative to construct the method so that both the pressure and
velocity remain in equilibrium without introducing any spurious oscillations for these mix-
ture cells. With applications to materials modeled by (1), some representative methods of
the previous efforts in this direction are the volume-of-fluid approach of Miller and Puckett
[31], the two-phase flow approach of Saurel and Abgrall [41], and the ghost-fluid approach
*of Fedkiw et al. [13]; see also [47] for an Lagrangian–Eulerian approach and [18, 42] for*
other up-to-date multicomponent algorithms.

With (1), our approach to model grid cells that contain more than one fluid component follows essentially the same idea as developed in [44, 45] for stiffened and van der Waals gases, and is a further generalization of the quasi-conservative method of Abgrall [1] for ideal gases. That is to say, we begin by considering an interface-only problem in one dimension where both the pressure and velocity are constants in the domain, while there are jumps in the other material-dependent variables across some interfaces. Then, from the energy equation, we derive a set of effective equations for the mixtures of the problem- dependent material quantities near the interfaces, see (11a) and (11b), so as to ensure the pressure remains in equilibrium for this problem. As in the previous work [45], in order to keep the material quantities unchanged as it should be in a single component region for a more general problem with shock and rarefaction waves as well, we proceed to modify these equations and obtain (11c) and (11d).

Note that here because of the strongly nonlinear coupling between many of the material
quantities in the Mie–Gr¨uneisen equation of state (1), see Section 2, it is not possible to
manipulate those equations further to find out a suitable effective equation for each of
the material quantities as we have hope for in the van der Waals gas case [45], which
yields the calculation of some of the material quantities from the equation of state an
explicit step out of the question. To remedy this situation, through a process of splitting
from the equation for the internal energy, we come up with a set of three equations, i.e.,
Eqs. (11c), (11e), and (11f ), which together with a local model based on the volume fraction
of fluid components within a grid cell can be used to the determination of the material-
dependent functions*0, p*ref*, and e*ref. Therefore, we are able to compute the pressure from
the equation of state in an easy manner with a reasonable amount of cost. A combination
of the Euler’s equation (2) with this set of three equations and the evolution equations
for volume fractions gives a complete model system that is a viable one to use in our
algorithm for numerical approximation of multicomponent problems. This will be discussed
further in Section 3 for the one-dimensional case, and Section 6 for the multidimensional
extension.

It should be mentioned that the multicomponent model we have derived, i.e., Eq. (12) or (23), is not written in the full conservation form, but is rather a quasi-conservative system of equations. Nevertheless, as in the case for single component flows, this model is a hyperbolic system when each physically relevant value of the state variables of the flow is defined in the region of thermodynamic stability; see Sections 2 and 3. As before (cf. [44, 45]), here we use the high-resolution method based on the wave-propagation viewpoint to compute approximate solution of the problem, giving an efficient implementation of the algorithm and also very accurate results for a variety of one- and two-dimensional problems; see Sections 5 and 6.1 for the details.

This paper is organized as follows. In Section 2, we discuss two important types of the curves (i.e., the isentropic and Hugoniot loci) for the reference states in the Mie–

Gr¨uneisen equation of state, and give some examples of interests for an explicit expression
of the material-dependent functions*0, p*ref*, and e*_{ref}. In Section 3, we describe in detail
the construction of our fluid-mixture type multicomponent model in one dimension. The
numerical method used to find approximate solution of the model system is briefly reviewed
in Section 4. This includes some discussion of the approximate Riemann solver of Roe. One-
dimensional results obtained using our multicomponent algorithm are shown in Section 5.

In Section 6, we extend the one-dimensional algorithm to multiple space dimensions, and show some numerical results in two dimensions.

**2. EQUATIONS OF STATE**

We are interested in a model for real materials (cf. [56]) where the thermodynamic behavior, such as the specific internal energy and the pressure of the material, can be characterized by the following two-terms relations

*e(V, T ) = e*ref*(V ) + e**T**(V, T ),* (3a)
*p(V, T ) = p*ref*(V ) +0(V )*

*V* *e**T**(V, T ).* (3b)

*Here V* *= 1/ρ denotes the specific volume, T denotes the temperature, and the subscripts*
*ref of ( p, e) and T of e refer to the “reference” and “thermal” states of the variables,*
*respectively. Note that to determine the value of T from those of V and e, we use the*
well-known relation in thermodynamics,

*e− e*ref*(V ) = e**T**(V, T ) = e*0*(V ) +*
Z *T*

*T*_{0}

*C**V**(V, T*^{0}*) dT*^{0}*,*

*where C**V* *is the specific heat at constant volume: Assume that C**V* depends only on the
*specific volume, from the above, we simply get e**T* *= C**V**(T − T*0*), yielding T = T*0*+ (e −*
*e*ref*)/C**V*. Clearly when we choose*ρ and e as our nature state-variables, from (3a) and (3b),*
we simply get the Mie–Gr¨uneisen equation of state (1).

Here, for simplicity, we assume that*0 is a function of V only, and takes the form*

*0(V ) = 0*0

µ*V*
*V*0

¶_{α}

*,* (4)

where*0*0*= γ*0*− 1 represents the Gr¨uneisen coefficient at V = V*0*, γ*0*> 1 is the usual def-*
inition of the ratio of specific heats, and*α ∈ [0, 1] is a dimensionless parameter. Depending*
*on the specific reference curve on which the states of the functions p*_{ref} *and e*_{ref}lie, the
*explicit relation between p*ref*and e*refwill be different. We next discuss two typical cases
of practical importance; see [26, 55] also for some other possible instances.

*2.1. Reference State along a Isentropic Locus*

*We begin by considering a class of materials where the thermodynamic state ( p*ref*, e*ref) of
*the model equation of state (1) lies along an isentropic locus from a centering point ( p*0*, e*0),

*i.e., the specific entropy, denoted by S, is a constant on the curve. In this case, from the*
*basic thermodynamics relation de*_{ref}*= T dS − p*ref*d V and d S*= 0, we obtain easily the
*condition between p*ref*and e*ref*as p*ref*(V ) = −de*ref*(V )/dV . Among many materials that*
belong to this type, in this paper we are mainly concerned with the following two sample
examples which have been used quite extensively for modeling the behavior of explosives
and other materials (cf. [35] for an example of solids).

(i) The Jones–Wilkins–Lee (JWL) equation of state (for gaseous explosives [9, 54]),
*0(V ) = 0*0

*e*ref*(V ) =AV*0

*R*1

exp

µ*−R*1*V*
*V*0

¶
+*BV*0

*R*2

exp

µ*−R*2*V*
*V*0

¶

*− e*0

(5)
*p*ref*(V ) = A exp*

µ*−R*1*V*
*V*0

¶

*+ B exp*

µ*−R*2*V*
*V*0

¶
*,*

(ii) The Cochran–Chan (CC) equation of state (for solid explosives [6]),
*0(V ) = 0*0

*e*ref*(V ) =* *−AV*0

1*− E*1

"µ
*V*
*V*_{0}

¶1*−E*1

− 1

#

+ *BV*0

1*− E*2

"µ
*V*
*V*_{0}

¶1*−E*2

− 1

#

*− e*0

(6)
*p*_{ref}*(V ) = A*

µ*V*
*V*0

¶_{−E}_{1}

*− B*
µ*V*

*V*0

¶_{−E}_{2}
*.*

Note that in each of these cases we have a total of seven material-dependent quantities in the
description of the material property, i.e., in the former case, there are*0*0*, V*0*, e*0*, A, B, R*1*,*
and*R*2, while in the latter case, there are*0*0*, V*0*, e*0*, A, B, E*1*, and E*2. Table I shows typical
set of numerical values for some sample materials of interest.

*2.2. Reference State along a Hugoniot Locus*

Our next example is concerned with a popular model for solid media such as metals. In this instance, in the absence of pronounced dynamic yielding effects or phase transitions, the hydrostatic pressure is commonly expressed by the Mie–Gr¨uneisen equation of state (1) together with a linear fit assumption for the shock velocity as a function of the particle velocity, i.e.,

*σ = c*0*+ s u.* (7)

Here*σ represents the shock velocity, c*0is the zero-pressure isentropic speed of sound, and
*s is a dimensionless parameter which is related to the pressure derivative of the isentropic*
*bulk modulus K**S**= ρ(∂p/∂ρ)|**S* by*(∂ K**S**/∂p)|**S* *= 4s − 1 (cf. [40]). By virtue of (7), it*
*is easy to deduce that the reference curve for ( p*_{ref}*, e*ref) is simply a single Hugoniot locus
*from an initial point ( p*0*, e*0). With this in mind, using the standard Rankine–Hugoniot jump
conditions for the Euler equations (2), after some simple algebraic manipulations, we find
*the explicit expression for p*_{ref}*and e*_{ref}as

*p*ref*(V ) = p*0+ *c*^{2}_{0}*(V*0*− V )*
*[V*0*− s(V*0*− V )]*^{2}

(8)
*e*ref*(V ) = e*0+1

2*[ p*ref*(V ) + p*0]*(V*0*− V );*

**TABLE I**

**Typical Material-Dependent Quantities for Three Different Models That**
**Are in the Mie–Gr ¨uneisen Form (1)**

JWL EOS *ρ*0*(kg/m*^{3}*)* *A (GPa)* *B (GPa)* *R*1 *R*2 *0*0 *α*

TNT 1840 854.5 20.5 4.6 1.35 0.25 0

Water 1004 1582 −4.67 8.94 1.45 1.17 0

CC EOS *ρ*0*(kg/m*^{3}*)* *A (GPa)* *B (GPa)* *E*1 *E*2 *0*0 *α*

Copper 8900 145.67 147.75 2.99 1.99 2 0

TNT 1840 12.87 13.42 4.1 3.1 0.93 0

Shock EOS *ρ*0*(kg/m*^{3}*)* *c*0(m/s) *s* *0*0 *α* *p*0 *e*0

Aluminum 2785 5328 1.338 2.0 1 0 0

Copper 8924 3910 1.51 1.96 1 0 0

Molybdenum 9961 4770 1.43 2.56 1 0 0

MORB 2660 2100 1.68 1.18 1 0 0

Water 1000 1483 2.0 2.0 10^{−4} 0 0

*Note. Data adapted from [26, 27, 55].*

see [28] for the details. Note that with*0 and (p*ref*, e*ref) defined by (4) and (8), respectively,
the resulting form of the Mie–Gr¨uneisen equation of state is often called the shock wave or
HOM equation of state [17, 26].

It had been discussed in detail (cf. [29]) that this shock wave equation of state has certain limitations. Nevertheless, it is observed experimentally that the model considered here is an adequate approximation for many metals, when the pressure is up to several megabars. A typical set of parameter values for metals, such as aluminum and copper, is given in Table I for the reference (cf. [27]). See [17, 40] for a more general discussion of the equation of state when (7) is replaced by a higher-order polynomial in the particle velocity.

It should be mentioned that to fulfill the conditions for the thermodynamic stability of
the materials of interests, we assume that for each given physical state the speed of sound
*c defined by*

*c*^{2}=
µ*∂p*

*∂ρ*

¶

*S*

=
µ*∂p*

*∂ρ*

¶

*e*

+ *p*
*ρ*^{2}

µ*∂p*

*∂e*

¶

*ρ*

= (9) µ

*0 + 1 + ρ0*^{0}
*0*

¶µ*p− p*ref

*ρ*

¶

*+ 0p*_{ref}

*ρ* *+ p*ref^{0} *− 0ρe*^{0}ref

belong to a set of real numbers, where*0*^{0}*, p*^{0}ref*, and e*^{0}_{ref}are the derivatives of*0, p*ref*, and e*ref

with respect to*ρ, respectively. Of course, it is both interesting and important to include the*
cavitation and spallation effects to materials modeled by (1) in a region where the pressure
drops to a critical value. But this subject matter is a very difficult one, and is beyond the
scope of this paper.

**3. EQUATIONS OF MOTION**

The basic governing equations in our multicomponent model consist of two parts. We
use (2) as a model system that describes the motion of the fluid mixtures of the conserved
variables*ρ, ρu, and ρE in a multicomponent grid cell. Assume a homogeneous flow with*
a single velocity and pressure on grid cells that contain more than one fluid components.

From the basic physical principles of mass and energy conservations, we derive a set of effective equations for the problem-dependent material functions in those cells (see below) that can be used easily to the determination of the pressure from the equation of state.

Combining these two set of the equations together with the equation of state constitutes a complete model system that is fundamental in our algorithm for numerical approximation of multicomponent problems.

To find out the aforementioned effective equations for the mixture of material quantities in a general Mie–Gr¨uneisen equation of state (1), similar to the previous work (cf. [44, 45]), we begin by considering an interface only problem where both the pressure and particle velocity are constants in the domain, while the other variables such as the density and the material quantities are having jumps across some interfaces. In this case, from the Euler Eqs. (2), it is easy to obtain equations for the time-dependent behavior of the density and total internal energy as

*∂ρ*

*∂t* *+ u∂ρ*

*∂x* *= 0,* (10a)

*∂*

*∂t(ρe) + u* *∂*

*∂x(ρe) = 0,* (10b)

in a respective manner. By inserting the Mie–Gr¨uneisen equation of state (1) into (10b), we have an equation of the form

*∂*

*∂t*

µ*p− p*ref

*0* *+ ρe*ref

¶
*+ u* *∂*

*∂x*

µ*p− p*ref

*0* *+ ρe*ref

¶

= 0 (10c)

that is in relation to not only the pressure, but also the material quantities appearing in the
functions*0, p*ref*, and e*ref.

In our algorithm, to maintain the pressure in equilibrium as it should be for our model
interface only problem, we split (10c) into the following two equations for the fluid mixtures
of 1*/ 0 and −(p*ref*/ 0) + ρe*refas

*∂*

*∂t*
µ1

*0*

¶
*+ u* *∂*

*∂x*
µ1

*0*

¶

*= 0,* (11a)

*∂*

*∂t*
µ

−*p*ref

*0* *+ ρe*ref

¶
*+ u* *∂*

*∂x*
µ

−*p*ref

*0* *+ ρe*ref

¶

*= 0,* (11b)

respectively. We emphasize that in order to have the correct pressure equilibrium in (10c)
near the interfaces, these are the two key equations that should be satisfied and approximated
consistently (when the problem is solved numerically) for any given expressions of*0, p*ref,
*and e*refappearing in the equation of state. As before (cf. [45]), because the solution of (11a)
and (11b) would depend on not only the material quantities, but also the density, to be able
to handle more general problems with shock and rarefaction waves, we need to modify each
of them so that the mass-conserving behavior of the solution in the single component region
can be obtained as well.

To accomplish this, consider the simpler case with (11a) as an example. Our basic ap- proach begins with a proper smoothness assumption of the density (such as in the case of rarefaction waves), and so we may apply the chain rule from differential calculus to the

partial derivatives in (11a), yielding easily the equivalent relation

*∂*

*∂t*
µ1

*0*

¶
*+ u* *∂*

*∂x*
µ1

*0*

¶

= −
µ*0*^{0}

*0*^{2}

¶ µ*∂ρ*

*∂t* *+ u∂ρ*

*∂x*

¶
*.*

Now by subtracting the term*(ρ0*^{0}*/ 0*^{2}*)∂u/∂x from the above relation on the both sides,*
and using the mass conservation on the right, we arrive at an equation of the form

*∂*

*∂t*
µ1

*0*

¶
*+ u* *∂*

*∂x*
µ1

*0*

¶

−
µ*0*^{0}

*0*^{2}

¶
*ρ∂u*

*∂x* = −
µ*0*^{0}

*0*^{2}

¶ µ*∂ρ*

*∂t* *+ u∂ρ*

*∂x* *+ ρ∂u*

*∂x*

¶

*= 0. (11c)*

Analogously, by following the same procedure as for (11a), the modification of (11b) takes the form

*∂*

*∂t*
µ

−*p*_{ref}
*0* *+ ρe*ref

¶
*+ u* *∂*

*∂x*
µ

−*p*_{ref}
*0* *+ ρe*ref

¶ +

µ*0*^{0}*p*_{ref}*− 0p*^{0}ref

*0*^{2} *+ e*ref*+ ρe*^{0}ref

¶
*ρ∂u*

*∂x* *= 0.*

(11d)
Clearly, (11c) and (11d) reduce to (11a) and (11b), respectively, for a solution near the
interfaces where*∂u/∂x = 0, and to the same mass conservation equation of the fluid mixture*
for a solution near rarefaction waves where the variation of*0*^{0}*/ 0*^{2}is smooth. Recall that
*0*^{0}*, p*^{0}ref*, and e*^{0}_{ref}are the derivatives of*0, p*ref*, and e*refwith respect to*ρ, respectively.*

Note that, at each space and time, given the initial conditions for 1*/ 0 and −(p*ref*/ 0) +*
*ρe*ref, to compute the solution of (11c) and (11d) would require five evaluations in total to the
mixtures such as*0*^{0}*, p*ref*, p*^{0}ref*, e*ref*, and e*^{0}_{ref}, from the equation of state. Here because of the
strongly nonlinear coupling between the material quantities (see Section 2 for an example),
from (11c) and (11d), it is not possible to come up with additional conditions for the further
details of the related parameters that makes the evaluation of any of the aforementioned
quantities in an explicit step. This is in contrast to the van der Waals case considered in
[45], and poses some difficulties in the realization of our multicomponent algorithm for
materials modeled by (1).

*To get by the problems involving the extra evaluations of the terms p*ref *and e*ref, in
particular, one simple way to do is to divide (11d) into the following two parts:

*∂*

*∂t*
µ*p*_{ref}

*0*

¶
*+ u* *∂*

*∂x*
µ*p*_{ref}

*0*

¶

*− ρ*

µ*0*^{0}*p*_{ref}*− 0p*^{0}ref

*0*^{2}

¶*∂u*

*∂x* *= 0,* (11e)

*∂*

*∂t(ρe*ref*) + u* *∂*

*∂x(ρe*ref*) + ρ*¡

*e*ref*+ ρe*ref^{0}

*¢∂u*

*∂x* *= 0.* (11f)

Clearly now instead of a single equation for*−(p*ref*/ 0) + ρe*ref, we have two separate
*ones for p*ref*/ 0 and ρe*ref, which together with the solutions of (11c) for*0 and the mass*
conservation equation in (2) for*ρ are sufficient to determine p*ref*and e*refwithout the use
of the equation of state. Of course, by doing so we still need to define*0*^{0}*, p*_{ref}^{0} *, and e*^{0}_{ref}so
as to have a working model system.

Note that if the reference state of the Mie–Gr¨uneisen equation of state (1) lies either along
an isentropic or a shock Hugoniot locus, from the basic thermodynamic relations described in
*Section 2, it is an easy matter to set e*^{0}_{ref}*= p*ref*/ρ*^{2}*or e*^{0}_{ref}*= [p*ref*/ρ + p*ref^{0} *(ρ/ρ*0*− 1)]/(2ρ)*
*in a respective manner, provided that the mixture of p*^{0}_{ref}has a proper mathematical definition

for numerical purpose also. Thus, in these two situations, to complete the model it is only
the mixtures of*0*^{0}*and p*^{0}_{ref}needed to be defined. (In fact, provided that some modification
of the equations is made, this is also true in a more general case where the reference state
lies partially on an isentrope and partially on a Hugoniot. But we will not discuss that case
here.) Although there may be other better ways, encouraged by the simplicity and also the
success of the previous work for a van der Waals gas case [45], we first introduce a local
model based on the volume-fraction formulation to the computation of all the remaining
undefined material quantities appearing in*0*^{0}*and p*^{0}_{ref}. Then we set the mixture states of*0*^{0}
*and p*_{ref}^{0} from the equation of state as in the single component case of the problem.

*To be more specific, consider an m-component flow problem with materials modeled*
by the shock wave equation of state (8), for example, we assign the material-dependent
mixtures:*α, ρ*0*, c*0*, and s, according to an averaging operatorM defined as follows:*

*M(z) =*
X*m*

*i*=1

*Y*^{(i)}*z*^{(i)}*,*

and compute *0*^{0}*= −α0/ρ, p*ref^{0} *= c*0^{2}*(1 − η)*^{2}*(1 + sη)/(1 − sη)*^{3} in an explicit manner.

*Here Y*^{(i)}*∈ [0, 1] is the volume-fraction function of the ith fluid component with a property*
*6**i** ^{m}*=1

*Y*

^{(i)}*= 1, z*

^{(i)}*is a material quantity belonging to the i th component, and*

*η = 1 −*

*(ρ*0

*/ρ). We use the evolution equation of the form*

*∂Y*^{(i)}

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

*∂x* *= 0,*

*for the motion of Y*^{(i)}*(see [31] for the other possibility in choosing the equation), i* =
1*, 2, . . . , m − 1, where u is the underlying particle velocity of the fluid mixture, and set*
*Y*^{(m)}*= 1 − 6**i** ^{m}*=1

^{−1}

*Y*

*. In summary, with the Mie–Gr¨uneisen equation of state (1), the mul- ticomponent model we proposed consists of the following system of equations,*

^{(i)}*∂ρ*

*∂t* + *∂*

*∂x(ρu) = 0*

*∂*

*∂t(ρu) +* *∂*

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

*∂*

*∂t(ρE) +* *∂*

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

*∂*

*∂t*
µ1

*0*

¶
*+ u* *∂*

*∂x*
µ1

*0*

¶

*− ρ*
µ*0*^{0}

*0*^{2}

¶*∂u*

*∂x* = 0 (12)

*∂*

*∂t*
µ*p*_{ref}

*0*

¶
*+ u* *∂*

*∂x*
µ*p*_{ref}

*0*

¶

*− ρ*

µ*0*^{0}*p*_{ref}*− 0p*_{ref}^{0}
*0*^{2}

¶*∂u*

*∂x* = 0

*∂*

*∂t(ρe*ref*) + u* *∂*

*∂x(ρe*ref*) + ρ(e*ref*+ ρe*^{0}_{ref}*)∂u*

*∂x* = 0

*∂Y*^{(i)}

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

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

*This gives us a system of m*+ 5 equations in total that is independent of the number of
material quantities involved in the equation of state (e.g., there are seven of them in (5) or
*(6)), for an m-component flow problem; m*≥ 1. It is easy to see 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 ones are the effective equations that*
are introduced to ensure the correct mixing of the problem-dependent material variables
near the interfaces. With a system written in this way, there is no problem to compute the
pressure from the equation of state

*p*=

·

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

µ*p*ref

*0*

¶

*− (ρe*ref*)*

¸Áµ1
*0*

¶
*.*

The initialization of the state variables in (12) for fluid-mixture cells can be made in a standard way as described in [45] for numerical simulation.

*Note that, when m*= 1 (single component flow), the effect to the introduction of the
equations for 1*/ 0, p*ref*/ 0, and ρe*refin the model is to reduce extra equation-of-state com-
putations in a numerical method to the least possible amount. It is easy to see that our
multicomponent model is a hyperbolic system by first writing (12) in a quasi-linear system
of equations

*∂q*

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

*∂x* *= 0.* (13)

*Here, for simplicity, in a two-component version of the model, we have the state vector q*
*and the matrix A defined by*

*q*=

·

*ρ, ρu, ρE,* 1
*0,* *p*ref

*0* *, ρe*ref*, Y*

¸*T*

and

*A*=

0 1 0 0 0 0 0

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

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

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

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

*−ψu* *ψ* 0 0 0 *u* 0

0 0 0 0 0 0 *u*

*.*

*We then compute the eigen-structure of the matrix A. It is a straightforward to show that*
*for each variables q defined in the region of thermodynamic stability the eigen-structure of*
*the matrix A possesses real eigenvalues*

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

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

1 1 1 0 0 0 0

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

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

*ϕ* 0 *ϕ* 1 0 0 0

*χ* 0 *χ* 0 1 0 0

*ψ* 0 *ψ* 0 0 1 0

0 0 0 0 0 0 1

(14b)

*with Ar**k**= λ**k**r**k**. Here K=0u*^{2}*/2, H = E + (p/ρ), ϕ = −0*^{0}*/ 0*^{2}*, χ = (0p*^{0}ref*− 0*^{0}*p*ref*)/*

*0*^{2}, and*ψ = e*ref*+ ρe*ref^{0} . Regarding discontinuous solutions of the system, such as shock
waves or contact discontinuities, it is not difficult to show that (12) has the usual form of
the Rankine–Hugoniot jump conditions across the waves; see Section 4.1 for more details.

**4. NUMERICAL METHODS**

To find approximate solutions of our model system (12) for multicomponent problems, we use a high-resolution wave propagation method developed by LeVeque [20, 23] for general hyperbolic systems of partial differential equations. This method is a variant of the fluctuation-and-signal scheme of Roe [37, 38] 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. To achieve high resolution, we introduce slopes and limiters to the method as in many other high resolution schemes for conservation laws [21, 50].

*4.1. Roe Riemann Solver*

Clearly, one of the major steps in our multicomponent algorithm is the numerical res-
olution of the Riemann problem at each cell interface. Here, with materials characterized
by the Mie–Gr¨uneisen equation of state (1), this amounts to solving the nonlinear system
*(12) with piecewise constant data q*_{L}*and q** _{R}* to the left and right of the interface. It is
well-known that, except under certain extreme conditions (cf. [31, 34, 55, 56]), the solution
of this Riemann problem would consist of two genuinely nonlinear waves, such as shock
and rarefaction, and a linearly degenerate wave (contact discontinuity); this is just like the
Riemann problem for a perfect gas (cf. [48]). In Fig. 1, we plot a typical solution structure
and the variables involved in the Riemann problem considered here. Because in general it

**FIG. 1.** Typical solution structure of the Riemann problem for our multicomponent model discussed in
Section 3. The key step in obtaining this solution is to find the midstate*(u*∗*, p*∗*) in the u − p phase plane. In*
general, it is a difficult task to do both exactly and efficiently.

is too complicated to solve the problem exactly, even in the single component case for real materials (cf. [8, 36, 43]), we discuss an approximate Riemann solver of Roe; see [31, 36]

for another approach based on the two-shock approximation.

In a Roe’s approximate Riemann solver, we replace the nonlinear Riemann problem mentioned above by a linear problem as

*∂q*

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

*∂x* *= 0,* *q(x, 0) =*

½*q**L* *for x* *< x*0

*q**R* *for x* *> x*0*,* (15)
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 in (13) about an average state. To find that matrix, as it is often done in*
many other Roe solvers (cf. [5, 14, 15]), we want to seek an average state that the difference
of the fluxes in the conservation part of (12) (i.e., the first three equations of the system) 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]*^{(i)}*,*

*for i* *= 1, 2, 3, where F ∈ R*^{3} 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 and ˆH 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*

*.* (16)

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

*1p =*

"µd
1
*0*

¶
*1*

µ*p*
*0*

¶

−µd
*p*
*0*

¶
*1*

µ1
*0*

¶# Á dµ
1
*0*

¶^{2}

is satisfied approximately (cf. [33] for an review of the other up-to-date approaches for real
*gases). With that we set ˆp= ( dp/ 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 ϕ, χ, and ψ. Since there is no unique way to*
do so, we might as well compute them using the Roe-average (16) also. It is our experience
that the set of average states described here is a reasonable one to use for many practical
multicomponent problems (see numerical results present in Section 5) as long as the flow
condition is not too extreme (i.e., with very large density and pressure ratios) across the
interfaces, (cf. [11, 45] for more discussions and the possible cures for that matters).

In contrast to the solution structure for a nonlinear Riemann problem (see Fig. 1), the
solution of the linear problem (15) consists of seven discontinuities propagating at constant
speeds (for a two-component system of seven equations). The jump across each discon-
tinuity 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* =
X7
*k*=1

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

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

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

*0*ˆ

*ˆc*

^{2}

·

−*u*ˆ^{2}

2 *1q*^{(1)}*+ ˆu1q*^{(2)}*− 1q*^{(3)}*+ ˆp1q*^{(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*),* *α*ˆ5*= 1q*^{(5)}*− ˆχ(1ρ − ˆα*2*),*
*α*ˆ6*= 1q** ^{(6)}*− ˆ

*ψ(1ρ − ˆα*2

*),*

*α*ˆ7

*= 1q*

^{(7)}*,*

(18)

*where ˆc*=p

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

Notice that in this Riemann solution, except the discontinuities for ˆ*λ*1*= ˆu − ˆc and ˆλ*3=
ˆ

*u+ ˆc, all the other discontinuities (five 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 ˆ*λ*2wave family, i.e., set*W*2*= ˆα*2*ˆr*2+
P7

*k*=4*α*ˆ*k**ˆr**k*. With this notation, we also write*W**k**= ˆα**k**ˆr**k*to represent the jump across the
*k-wave for k*= 1 or 3. Thus, without causing any confusion, we may assume that the wave
family in total is 3 for the solution of this Riemann problem.

*4.2. High-Resolution Wave Propagation Scheme*

Consider a uniform grid with fixed mesh spacing*4x, for example. 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

*4x*

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.*

In this numerical discretization 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**,* (19)

where*λ**k**∈ R and W**k*∈ R^{m}*are solutions of the kth wave family, for k= 1, 2, . . . , m**w*,
*obtained from solving the Riemann problems at cell interfaces x**j**and x**j*+1; see Section 4.1.

As usual, we define*λ*^{−}*= min(λ, 0) and λ*^{+}*= max(λ, 0). Clearly, the method belongs to a*
class of upwind schemes (cf. [15, 21]), and by following the same procedure as described
in [45], it is quasi-conservative in the sense that when applying the method to (12) 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 achieve high resolution in this method, we begin by introducing correction waves
in a piecewise-linear form with zero mean value. We then propagate each wave over the
time step*1t, and update the cell averages it overlaps. Without going into the detail here*

(cf. [24]), with the corrections, (19) is modified by

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

*m*_{w}

X

*k*=1

·

*|λ**k*|
µ

1*−|λ**k*|*1t*
*1x*

¶
*W**k*

¸*n*
*j*+1−

·

*|λ**k*|
µ

1*−|λ**k*|*1t*
*1x*

¶
*W**k*

¸*n*
*j*

*. (20)*

It is important to mention that, in practice, the jump of each wave in the above formula
should be limited by using a “slope-limiter” (cf. [21]) to avoid unnecessary fluctuations
near discontinuities. We want to do this by replacing each*W**k*in (20) with a limited value
*W*f*k*obtained by comparing*W**k*with the corresponding*W**k*from the neighboring Riemann
problem to the left (if*λ**k**> 0) or to the right (if λ**k**< 0).*

Now with the use of the Roe solver to the computations, it is quite common to limit over
each strength of the wave ˆ*α**k j* *via a limiter functionφ (e.g., by using the minmod function*
*φ(θ) = max(0, min(1, θ)) or some others as discussed in [50]), and set*

*α*˜*k j* *= φ(θ**k j**) ˆα**k j* with *θ**k j* =*α*ˆ*k J*

*α*ˆ*k j*

*, J =*

(*j− 1 if ˆλ**k j* ≥ 0

*j+ 1 if ˆλ**k j**< 0,* (21)
*for k= 1, 2, . . . , 7 (cf. [15, 22, 23]). In this approach, we then replace the waves in (20),*

*(W*1*, W*2*, W*3*) =*
µ

*α*ˆ1*ˆr*1*, ˆα*2*ˆr*2+
X7

*k*=4

*α*ˆ*k**ˆr**k**, ˆα*3*ˆr*3

¶
*,*

by a limited version as

*(fW*1*, fW*2*, fW*3*) =*
µ

*α*˜1*ˆr*1*, ˜α*2*ˆr*2+
X7
*k*=4

*α*˜*k**ˆr**k**, ˜α*3*ˆr*3

¶
*.*

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 (20). Moreover, we obtain a better resolution of the result as compared to the first- order result. Concerning stability of the method, it is observed numerically that the method is stable under the usual CFL (Courant–Friedrichs–Lewy) condition for hyperbolic systems of conservation laws; see Section 5 for an example.

**5. NUMERICAL RESULTS IN ONE DIMENSION**

We now present some sample numerical results obtained using our multicomponent algorithm with the Roe solver described in Section 4.

*5.1. Single-Component Case*

As a preliminary, we begin by showing results for problems with only a single fluid component presence in the problem formulation.

EXAMPLE5.1. Our first test problem is a Riemann problem in a shock tube with the material inside the tube modeled by the Jones–Wilkins–Lee equation of state (5). For comparison purposes, we take the similar initial data as studied by Rider [36], where on the

**FIG. 2.** High-resolution results for a single component Riemann problem with the gaseous explosives at time
*t**= 12 µs. The solid line is the fine grid solution computed by 1x = 1/2000, and the points show the solution*
with*1x = 1/100. The dashed line in each subplot is the initial condition at time t = 0. The gaseous explosive is*
modeled by the Jones–Wilkins–Lee equation of state (5).

left of the interface, 0*≤ x < 1/2 m, we have*

*(ρ, u, p, e*0*)**L* *= (1700 kg/m*^{3}*, 0 m/s, 10*^{12} Pa*, 0 kJ/kg),*
and on the right of the interface, 1*/2 m ≤ x ≤ 1 m, we have*

*(ρ, u, p, e*0*)**R**= (1000 kg/m*^{3}*, 0 m/s, 5 × 10*^{10}Pa*, 0 kJ/kg).*

In this problem, the seven material-dependent quantities:*ρ*0*, A, B, R*1*, R*2*, 0*0, and*α, have*
been chosen for the product gases of the explosive TNT as given in Table I.

In Fig. 2, we show results for the density, velocity, pressure, and the speed of sound at
*time t= 12 µs, where the test has been carried out by using the high-resolution method*
with theMINMODlimiter, the Courant number*µ = 0.9, and the mesh size 1x = 1/100. By*
comparing the computed solution with the fine grid solution obtained using the same method
but*1x = 1/2000, we observe good agreement in the region of rarefaction wave where the*
flow is smooth, and reasonable resolution in the region of shock and contact discontinuity
where the flow is not smooth ( judging from the approximate location and the monotonicity
of the solution profile for the discontinuity). In addition, it is easy to make comparisons and
see that our solution agrees quite well with the result present in [36] where a MUSCL-type
scheme with an approximate Riemann solver based on the two-shock approximation was
used in the computation.

EXAMPLE5.2. We are next concerned with an impact problem in which a precompressed
semi-infinite aluminum slab at rest with*(ρ, p) = (4000 kg/m*^{3}*, 7.93 × 10*^{9}Pa) is being hit
by an ambient aluminum slab traveling at the speed 2 km/s from the right to the left with
the reference state*(ρ, p) = (ρ*0*, p*0*). As in [29, 31, 36] and references therein, we use the*
popular shock wave equation of state (8) to model the thermodynamical behavior of the
aluminum; see Table I for the numerical values of the material constants:*ρ*0*, c*0*, s, 0*0,and*α.*

In this setup, it is not difficult to show that the exact solution of this problem would consist of a leftward going shock wave to the stationary aluminum, a material interface,

**FIG. 3.** High-resolution results for a single component impact problem with two aluminum slabs at time
*t**= 50 µs. The aluminum is modeled by the shock wave equation of state (8). The graphs of the solutions are*
displayed in the same manner as in Fig. 2.

and a rightward going shock wave to the moving aluminum. Figure 3 shows the numerical
*result for this problem at time t= 50 µs. As compared to the fine grid solution, which is a*
good approximation to the exact solution, it is clear that our result gives the correct solution
behavior of this problem; see [36] also for a similar calculation. Here the computation was
performed in the same manner as in Example 5.1, where the initial point of the projectile
impact was set at the center of a meter-wide computational domain.

*5.2. Multicomponent Case*

We now show results for examples with more than one fluid component in the problem formulation.

EXAMPLE5.3. To begin, we are interested in a two-component impact problem of Saurel
and Abgrall [41]. Initially, under the atmospheric condition (i.e., with uniform pressure
*p*_{0}*= 1atm and temperature T*0= 300 K throughout the domain), there is a rightward going
*copper plate with the speed u*= 1500 m/s interacting with a solid explosive (considered as
an inert material) at rest on the right of the plate. In this problem, to model the material
properties of the copper and (solid) explosive, we use the same Cochran–Chan equation
of state (6), but with a different set of material-dependent quantities for each of them, see
numerical values given in Table I.

As in Example 5.2, the exact solution of this impact problem is composed of a leftward-
going shock wave to the copper, a rightward-going shock waves to the inert explosive, and a
material interface lying in between that separates these two different materials. We run this
problem using exactly the same method as performed in the previous examples for single
*component flow, and show the resulting solution in Fig. 4 at time t* *= 85 µs for the density,*
velocity, pressure, and the thermal internal energy. By comparing the computed solution with
the fine grid one obtained using the same method but*1x = 1/2000, we observe reasonable*
behavior of the solution with the correct shock speeds and free of spurious oscillations in
the pressure near the interface. Checking our result with the displayed solution appearing

**FIG. 4.** *High-resolution results for a two-component (solid explosive-copper) impact problem at time t*=
85*µs. The solid line is the fine grid solution computed by 1x = 1/2000, and the points show the solution with*
*1x = 1/100. The dashed line in each subplot is the initial condition at time t = 0. Both the solid explosive and*
copper are modeled by the Cochran–Chan equation of state (6), but with a different set of material quantities for
each of them.

in [41] with the same mesh size*1x = 1/100, we find excellent agreement in the density,*
pressure, and velocity. Clearly, for detonation problems, it is often necessary to report the
*solution of the temperature T as well. As we have seen in the figure (see Fig. 5 also), the*
*algorithm did quite a good job to the resolution of thermal internal energy e** _{T}*which can be

*computed directly from the variables obtained in the algorithm, i.e., e*

*T*

*= (p − p*ref

*)/(ρ0).*

*To go one step further to T by T* *= e**T**/C**V*, we need to do some postprocessing work for the
*fluid mixture C**V**. Although there are many ways to get C**V*, say by using the volume-fraction

**FIG. 5.** High-resolution results for a two-component (gaseous explosive–copper) Riemann problem at time
*t**= 73 µs. The gaseous explosive is modeled by the Jones–Wilkins–Lee equation of state (5), while the copper is*
modeled by the Cochran–Chan equation of state (6). The graphs of the solutions are displayed in the same manner
as in Fig. 4.

function, for example, this is not really in the heart of the whole algorithm, and so the plot
*of the temperature is not shown here. Note that, we have C*_{V}*= 393 and 1087 J/(kg · K) for*
copper and explosive, respectively.

EXAMPLE5.4. Our next example concerns a two-component Riemann problem of Saurel and Abgrall [41] that involves the interaction of gaseous detonation products with a copper plate. In this test, as in Example 5.3, copper is modeled by the Cochran–Chan equation of state (6), while the detonation products are modeled by the widely employed Jones–

*Wilkins–Lee equation of state (5). Initially, on the left when x∈ [0, 0.5)m, we have the*
detonation product with the data

*(ρ, u, p, e*0*)**L**= (2485.37 kg/m*^{3}*, 0, 3.7 × 10*^{10}Pa*, 8149.158 kJ/kg),*
*and on the right when x∈ [0.5, 1]m, we have the copper with data*

*(ρ, u, p, e*0*)**R* *= (8900 kg/m*^{3}*, 0, 10*^{5} Pa*, 117.9 kJ/kg).*

We note that the data on the left is at the Chapman–Jouget state (see [41] for the details), while the data on the right is at the usual atmospheric conditions. In Table I, we list the material quantities of these two substances to this run. For this problem, it is known that the exact solution consists of a shock wave moving to the right in the copper and a rarefaction wave propagating to the left in the explosive; see [41].

To solve this problem numerically, we need to define a hybrid version of the equation of state that is necessary in the algorithm for the numerical mixing between these two different materials. This can be done by following the same approach as described in [45]

for a case with the mixing between stiffened and van der Waals gases, yielding easily a Mie–Gr¨uneisen equation of state of the form

*p(ρ, e) = ˜p*ref*(ρ) + ˜0ρ[e − ˜e*ref*(ρ)]* (22)
*for the copper-explosive mixture, where ˜p*ref*= p*ref^{(JWL)}*+ p*ref^{(CC)}*and ˜e*ref*= e** ^{(JWL)}*ref

*+ e*

*ref*

^{(CC)}*are defined by simply combining the two different p*ref*and e*reffrom (5) and (6) into one,
respectively, and ˜*0 = 0*0. Here the computation was performed in the same way as before,
*and the results are shown in Fig. 5 at time t= 73 µs for the variables ρ, u, p, and e**T* also.

Comparing our solution with the one shown in [41] using a two-phase flow solver, we again observe good agreement for this problem.

EXAMPLE5.5. To end this section, we test our algorithm for a model shock-contact prob-
lem that involves the interaction of a shock wave in molybdenum and an encapsulated MORB
(Mid-Ocean Ridge Basalt) liquid (this problem is motivated by a two-dimensional test of
Miller and Puckett [31]). The initial condition is composed of a stationary (molybdenum-
*MORB) interface at x* *= 0.6 m and a rightward going Mach 1.163 shock wave in molybde-*
*num at x= 0.4 m traveling from left to right in a shock tube of unit length. The material on*
the right of the interface is a MORB liquid modeled by the shock wave equation of state (8)
with the data

*(ρ, u, p, e*0*)**R* *= (2260 kg/m*^{3}*, 0 m/s, 0 Pa, 0 kJ/kg),*

and the material on the left of the interface (i.e., on the middle and the preshock state) is molybdenum modeled by the shock wave equation of state also with data

*(ρ, u, p, e*0*)**M**= (9961 kg/m*^{3}*, 0 m/s, 0 Pa, 0 kJ/kg).*

**FIG. 6.** High-resolution results for a shock wave in molybdenum interacting with an encapsulated MORB
liquid at time 120*µs. The solid line is the fine grid solution computed by 1x = 1/2000, and the points show*
the solution with*1x = 1/100. The dashed line in each subplot is the initial condition at time t = 0. Both the*
molybdenum and MORB are modeled by the shock wave equation of state (8), but with different material constants
for each of them.

The state behind the shock in the molybdenum is

*(ρ, u, p, e*0*)**L* *= (11042 kg/m*^{3}*, 543 m/s, 3 × 10*^{10} Pa*, 0 kJ/kg),*

see the dashed line shown in Fig. 6 for illustration. We note that this gives us one example in which the (molybdenum-MORB) interface is accelerated by a shock wave coming from the heavy-fluid to the light-fluid region, and it is known that the resulting wave pattern after the interaction would consist of a transmitted shock wave, an interface, and a reflected rarefaction wave (cf. [4, 16]).

*Numerical results for this problem are shown in Fig. 6 at time t= 120 µs for the states*
*ρ, u, p, and 0. Clearly, we observe sensible resolution and convergence of the solution*
structure as the mesh is refined. Note that because of the passage of the transmitted
shock wave, the MORB liquid is compressed, yielding the increase of the density, ve-
locity, and pressure. A two-dimensional version of this problem will be considered in
Section 6.1.

**6. EXTENSION TO MULTIPLE DIMENSIONS**

The multidimensional version of our model system (12) for compressible multicomponent problems with the Mie–Gr¨uneisen equation of state (1) takes the form

*∂ρ*

*∂t* +
X*N*

*j*=1

*∂*

*∂x**j*

*(ρu**j**) = 0*

*∂*

*∂t(ρu**i**) +*
X*N*

*j*=1

*∂*

*∂x**j*

¡*ρu**i**u*_{j}*+ δ*^{i j}*p*¢

*= 0 for i = 1, 2, . . . , N*