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) = pref(ρ) + 0(ρ)ρ [e − eref(ρ)]. (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 pref, erefare 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 expressions0, pref, and erefis 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(ρu2+ p) = 0, (2)
∂
∂t(ρE) + ∂
∂x(ρEu + pu) = 0,
respectively, where u is the particle velocity, and E= e + u2/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 functions0, pref, and eref. 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 functions0, pref, and eref. 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 ) = eref(V ) + eT(V, T ), (3a) p(V, T ) = pref(V ) +0(V )
V eT(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− eref(V ) = eT(V, T ) = e0(V ) + Z T
T0
CV(V, T0) dT0,
where CV is the specific heat at constant volume: Assume that CV depends only on the specific volume, from the above, we simply get eT = CV(T − T0), yielding T = T0+ (e − eref)/CV. 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 that0 is a function of V only, and takes the form
0(V ) = 00
µV V0
¶α
, (4)
where00= γ0− 1 represents the Gr¨uneisen coefficient at V = V0, γ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 pref and ereflie, the explicit relation between prefand erefwill 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 ( pref, eref) of the model equation of state (1) lies along an isentropic locus from a centering point ( p0, e0),
i.e., the specific entropy, denoted by S, is a constant on the curve. In this case, from the basic thermodynamics relation deref= T dS − prefd V and d S= 0, we obtain easily the condition between prefand erefas pref(V ) = −deref(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 ) = 00
eref(V ) =AV0
R1
exp
µ−R1V V0
¶ +BV0
R2
exp
µ−R2V V0
¶
− e0
(5) pref(V ) = A exp
µ−R1V V0
¶
+ B exp
µ−R2V V0
¶ ,
(ii) The Cochran–Chan (CC) equation of state (for solid explosives [6]), 0(V ) = 00
eref(V ) = −AV0
1− E1
"µ V V0
¶1−E1
− 1
#
+ BV0
1− E2
"µ V V0
¶1−E2
− 1
#
− e0
(6) pref(V ) = A
µV V0
¶−E1
− B µV
V0
¶−E2 .
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 are00, V0, e0, A, B, R1, andR2, while in the latter case, there are00, V0, e0, A, B, E1, and E2. 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.,
σ = c0+ s u. (7)
Hereσ represents the shock velocity, c0is 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 KS= ρ(∂p/∂ρ)|S by(∂ KS/∂p)|S = 4s − 1 (cf. [40]). By virtue of (7), it is easy to deduce that the reference curve for ( pref, eref) is simply a single Hugoniot locus from an initial point ( p0, e0). 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 prefand erefas
pref(V ) = p0+ c20(V0− V ) [V0− s(V0− V )]2
(8) eref(V ) = e0+1
2[ pref(V ) + p0](V0− V );
TABLE I
Typical Material-Dependent Quantities for Three Different Models That Are in the Mie–Gr ¨uneisen Form (1)
JWL EOS ρ0(kg/m3) A (GPa) B (GPa) R1 R2 00 α
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/m3) A (GPa) B (GPa) E1 E2 00 α
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/m3) c0(m/s) s 00 α p0 e0
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 with0 and (pref, eref) 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
c2= µ∂p
∂ρ
¶
S
= µ∂p
∂ρ
¶
e
+ p ρ2
µ∂p
∂e
¶
ρ
= (9) µ
0 + 1 + ρ00 0
¶µp− pref
ρ
¶
+ 0pref
ρ + pref0 − 0ρe0ref
belong to a set of real numbers, where00, p0ref, and e0refare the derivatives of0, pref, and eref
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− pref
0 + ρeref
¶ + u ∂
∂x
µp− pref
0 + ρeref
¶
= 0 (10c)
that is in relation to not only the pressure, but also the material quantities appearing in the functions0, pref, and eref.
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 −(pref/ 0) + ρerefas
∂
∂t µ1
0
¶ + u ∂
∂x µ1
0
¶
= 0, (11a)
∂
∂t µ
−pref
0 + ρeref
¶ + u ∂
∂x µ
−pref
0 + ρeref
¶
= 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 of0, pref, and erefappearing 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
¶
= − µ00
02
¶ µ∂ρ
∂t + u∂ρ
∂x
¶ .
Now by subtracting the term(ρ00/ 02)∂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
¶
− µ00
02
¶ ρ∂u
∂x = − µ00
02
¶ µ∂ρ
∂t + u∂ρ
∂x + ρ∂u
∂x
¶
= 0. (11c)
Analogously, by following the same procedure as for (11a), the modification of (11b) takes the form
∂
∂t µ
−pref 0 + ρeref
¶ + u ∂
∂x µ
−pref 0 + ρeref
¶ +
µ00pref− 0p0ref
02 + eref+ ρe0ref
¶ ρ∂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 of00/ 02is smooth. Recall that 00, p0ref, and e0refare the derivatives of0, pref, and erefwith respect toρ, respectively.
Note that, at each space and time, given the initial conditions for 1/ 0 and −(pref/ 0) + ρeref, to compute the solution of (11c) and (11d) would require five evaluations in total to the mixtures such as00, pref, p0ref, eref, and e0ref, 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 pref and eref, in particular, one simple way to do is to divide (11d) into the following two parts:
∂
∂t µpref
0
¶ + u ∂
∂x µpref
0
¶
− ρ
µ00pref− 0p0ref
02
¶∂u
∂x = 0, (11e)
∂
∂t(ρeref) + u ∂
∂x(ρeref) + ρ¡
eref+ ρeref0
¢∂u
∂x = 0. (11f)
Clearly now instead of a single equation for−(pref/ 0) + ρeref, we have two separate ones for pref/ 0 and ρeref, which together with the solutions of (11c) for0 and the mass conservation equation in (2) forρ are sufficient to determine prefand erefwithout the use of the equation of state. Of course, by doing so we still need to define00, pref0 , and e0refso 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 e0ref= pref/ρ2or e0ref= [pref/ρ + pref0 (ρ/ρ0− 1)]/(2ρ) in a respective manner, provided that the mixture of p0refhas a proper mathematical definition
for numerical purpose also. Thus, in these two situations, to complete the model it is only the mixtures of00and p0refneeded 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 in00and p0ref. Then we set the mixture states of00 and pref0 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, c0, and s, according to an averaging operatorM defined as follows:
M(z) = Xm
i=1
Y(i)z(i),
and compute 00= −α0/ρ, pref0 = c02(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 6im=1Y(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 − 6im=1−1Y(i). In summary, with the Mie–Gr¨uneisen equation of state (1), the mul- ticomponent model we proposed consists of the following system of equations,
∂ρ
∂t + ∂
∂x(ρu) = 0
∂
∂t(ρu) + ∂
∂x(ρu2+ p) = 0
∂
∂t(ρE) + ∂
∂x(ρEu + pu) = 0
∂
∂t µ1
0
¶ + u ∂
∂x µ1
0
¶
− ρ µ00
02
¶∂u
∂x = 0 (12)
∂
∂t µpref
0
¶ + u ∂
∂x µpref
0
¶
− ρ
µ00pref− 0pref0 02
¶∂u
∂x = 0
∂
∂t(ρeref) + u ∂
∂x(ρeref) + ρ(eref+ ρe0ref)∂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ρ +
µpref
0
¶
− (ρeref)
¸Áµ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, pref/ 0, and ρerefin 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, pref
0 , ρeref, Y
¸T
and
A=
0 1 0 0 0 0 0
K− u2 u(2 − 0) 0 −p0 −0 0 0
u(K − H) H − u20 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= (r1, r2, . . . , r7) =
1 1 1 0 0 0 0
u− c u u+ c 0 0 0 0
H− uc u2/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 Ark= λkrk. Here K=0u2/2, H = E + (p/ρ), ϕ = −00/ 02, χ = (0p0ref− 00pref)/
02, andψ = eref+ ρeref0 . 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 qL and qR 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(qL, qR)∂q
∂x = 0, q(x, 0) =
½qL for x < x0
qR for x > x0, (15) where ˆA(qL, qR) 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)= (FR− FL)(i)= [ ˆA(qL, qR)(qR− qL)](i)= [ ˆA(qL, qR)1q](i),
for i = 1, 2, 3, where F ∈ R3 is the usual definition of the fluxes for conservation laws, and1F(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= √ρLzL+ √ρRzR
√ρL+ √ρR
. (16)
Note that in the process of the derivation, as in [44, 45], we have chosen the averages( d1/ 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)/( d1/ 0) and ˆ0 = 1/( d1/ 0). To finish the construction of Aˆ(qL, qR), 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 = qR− qL = X7 k=1
αˆkˆrk, (17)
where ˆrkis the kth eigenvector of ˆA; see (14a) and (14b). The scalar ˆαkgives the strength across the discontinuity that can be determined easily from (17). We find
αˆ2= 1q(1)+ 0ˆ ˆc2
·
−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 W2 defined by combining all the jumps across the ˆλ2wave family, i.e., setW2= ˆα2ˆr2+ P7
k=4αˆkˆrk. With this notation, we also writeWk= ˆαkˆrkto 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 spacing4x, for example. We use a standard finite-volume formulation in which the value Qnj approximates the cell average of the solution over the grid cell [xj, xj+1] at time tn:
Qnj ≈ 1 4x
Z xj+1 xj
q(x, tn) dx.
The time step from the current time tnto the next tn+1is denoted by1t.
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
Qnj+1 = Qnj− 1t 1x
mw
X
k=1
(λ−kWk)nj+1+ (λ+kWk)nj, (19)
whereλk∈ R and Wk∈ Rm are solutions of the kth wave family, for k= 1, 2, . . . , mw, obtained from solving the Riemann problems at cell interfaces xjand xj+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 step1t, and update the cell averages it overlaps. Without going into the detail here
(cf. [24]), with the corrections, (19) is modified by
Qnj+1:= Qnj+1− 1t 21x
mw
X
k=1
·
|λk| µ
1−|λk|1t 1x
¶ Wk
¸n j+1−
·
|λk| µ
1−|λk|1t 1x
¶ Wk
¸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 eachWkin (20) with a limited value Wfkobtained by comparingWkwith the correspondingWkfrom 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),
(W1, W2, W3) = µ
αˆ1ˆr1, ˆα2ˆr2+ X7
k=4
αˆkˆrk, ˆα3ˆr3
¶ ,
by a limited version as
(fW1, fW2, fW3) = µ
α˜1ˆr1, ˜α2ˆr2+ X7 k=4
α˜kˆrk, ˜α3ˆr3
¶ .
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 with1x = 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, e0)L = (1700 kg/m3, 0 m/s, 1012 Pa, 0 kJ/kg), and on the right of the interface, 1/2 m ≤ x ≤ 1 m, we have
(ρ, u, p, e0)R= (1000 kg/m3, 0 m/s, 5 × 1010Pa, 0 kJ/kg).
In this problem, the seven material-dependent quantities:ρ0, A, B, R1, R2, 00, 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 but1x = 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/m3, 7.93 × 109Pa) 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, p0). 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, c0, s, 00,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 p0= 1atm and temperature T0= 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 but1x = 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 size1x = 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 eTwhich can be computed directly from the variables obtained in the algorithm, i.e., eT = (p − pref)/(ρ0).
To go one step further to T by T = eT/CV, we need to do some postprocessing work for the fluid mixture CV. Although there are many ways to get CV, 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 CV = 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, e0)L= (2485.37 kg/m3, 0, 3.7 × 1010Pa, 8149.158 kJ/kg), and on the right when x∈ [0.5, 1]m, we have the copper with data
(ρ, u, p, e0)R = (8900 kg/m3, 0, 105 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) = ˜pref(ρ) + ˜0ρ[e − ˜eref(ρ)] (22) for the copper-explosive mixture, where ˜pref= pref(JWL)+ pref(CC) and ˜eref= e(JWL)ref + e(CC)ref
are defined by simply combining the two different prefand ereffrom (5) and (6) into one, respectively, and ˜0 = 00. 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 eT 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, e0)R = (2260 kg/m3, 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, e0)M= (9961 kg/m3, 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 with1x = 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, e0)L = (11042 kg/m3, 543 m/s, 3 × 1010 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 + XN
j=1
∂
∂xj
(ρuj) = 0
∂
∂t(ρui) + XN
j=1
∂
∂xj
¡ρuiuj+ δi jp¢
= 0 for i = 1, 2, . . . , N
