A Riemann problem solution methodology for a class of evolutionary mixture equations with an arbitrary number of components

Applied Numerical Mathematics 76 (2014) 145–165

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A Riemann problem solution methodology for a class of evolutionary mixture equations with an arbitrary number of components M.W. Crochet 1 , K.A. Gonthier ∗,2 Department of Mechanical and Industrial Engineering, Louisiana State University, Baton Rouge, LA 70803, USA

a r t i c l e

i n f o

Article history: Received 29 April 2013 Received in revised form 28 August 2013 Accepted 13 September 2013 Available online 19 September 2013 Keywords: Multiphase flows Granular mixtures Riemann solver Hyperbolic systems Nonlinear algebraic equations Nonconservative sources

a b s t r a c t The solution of the two-phase Riemann problem is a critical component of upwind finitevolume numerical schemes used to solve systems of evolutionary equations, which are routinely used to model compaction and combustion phenomena in gas–granular explosive mixtures. Extensions of a common two-phase model are currently being used to analyze the thermomechanics and combustion of explosive mixtures consisting of N components. Although a solution to the two-phase Riemann problem has been formulated, there is currently no available analogue for the N-phase system in the literature, due to the inherent difficulty of determining the correct wave ordering within the Riemann solver. The development of a solution for these systems is therefore an important step in the formulation of numerical schemes applied to N-phase mixtures. Here, an extension of the exact two-phase solution methodology is proposed for the N-phase case, which may be utilized in the construction of finite-volume schemes for multiphase systems, and can be used with general, convex equations of state. Finally, example problems for three-phase mixtures are considered to illustrate the accuracy of the solution compared to the results of a centered numerical scheme. These solutions also demonstrate the complexity of the possible wave configurations that arise when multiple solid phases are present, as well as the algorithmic challenges which must be addressed to provide a robust implementation. © 2013 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Systems of hyperbolic partial differential equations are frequently used to analyze wave phenomena in two-phase flows [3,2,23,24,12]. In particular, the two-phase Baer–Nunziato (BN) model [4] was formulated to study deflagration-to-detonation transition (DDT) in granular energetic materials, where a mixture of reactant granular solid and product gas simultaneously exist. The system of equations consists of mass, momentum, and energy balance equations for each phase, as well as a volume fraction evolution equation required for closure. Volumetric source terms in these equations account for mass, momentum, and energy interactions between phases due to combustion, drag, heat transfer, and compaction. In this study we ignore phase interaction processes and focus on terms that account for nonlinear wave formation and propagation, which can be computationally challenging to resolve. The effects of volumetric source terms can be computationally examined separately based on standard operator splitting techniques [27].

* 1 2

Corresponding author. Tel.: +1 225 578 5915. E-mail addresses: [email protected] (M.W. Crochet), [email protected] (K.A. Gonthier). Graduate Research Assistant. Associate Professor. Tel.: +1 225 578 5915.

0168-9274/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2013.09.004

146

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

The one-dimensional form of the two-phase equations that accounts for nonlinear convection is given by:

∂φs ∂ q ∂ f(q) + = g(q) , ∂t ∂x ∂x

(1)

where

q = [φs ρs , φ g ρ g , φs ρs u s , φ g ρ g u g , φs ρs E s , φ g ρ g E g , φs ] ,



f(q) = φs ρs u s , φ g ρ g u g , φs ρs u 2s + φs P s , φ g ρ g u 2g + φ g P g , φs ρs u s ( E s + P s /ρs ), φ g ρ g u g ( E g + P g /ρ g ), 0



(2)

,



g(q) = [0, 0, P g , − P g , P g u s , − P g u s , −u s ] .

(3) (4)

E i = e i + u 2i /2

Here, t is time and x is spatial position. For i = s, g, φi is volume fraction; ρi is density; u i is velocity; is the total specific energy, where e i is the specific internal energy, and P i is pressure. The subscripts s and g denote the solid and gas phases, respectively. This system is closed by equations of state (EOS’s) e i = e i (ρi , P i ) for both the gas and solid phases, and the saturation constraint φs + φ g = 1. The compaction equation convects φs at the local velocity u s . Some multiphase models convect φs at a different velocity U [24,18], such as the mass-weighted velocity of the phases. The Riemann solver presented here is generally not applicable in such cases. The sources g(q)∂φs /∂ x are referred to as nonconservative products in the literature, and are required to satisfy the strong form of the entropy inequality for the mixture. Formulating accurate discretizations of these products within finite-volume schemes has proven to be particularly challenging [17,6], and their proper mathematical treatment remains an active area of research [16,11,9,22,20,21]. The two-phase BN model has been extended to include mixtures containing an arbitrary number of condensed phases to examine the thermomechanics and combustion of multi-component energetic solids [14,28,5]. Here, the solid subscript s is replaced by the index i, for i = 1, 2, . . . , M, where M is the total number of solid phases. The gas–solid mixture consists of N = M + 1 components, and the one-dimensional system of equations is given by:

∂φi ∂ q ∂ f(q)  + = gi (q) , ∂t ∂x ∂x M

(5)

i =1

where



q = {φi ρi }iM=1 , φ g ρ g , {φi ρi u i }iM=1 , φ g ρ g u g , {φi ρi E i }iM=1 , φ g ρ g E g , {φi }iM=1



f(q) = {φi ρ

M i u i } i =1 , φ g



2 i ui

M



,

(6)

2 gug

ρ g u g , φ i ρ + φ i P i i =1 , φ g ρ + φg P g ,   M φi ρi u i ( E i + P i /ρi ) i =1 , φ g ρ g u g ( E g + P g /ρ g ), {0}iM=1 ,   gi (q) = {0}iM=1 , 0, P g ai , − P g , P g u i ai , − P g u i , −u i ai .

(7) (8)

The notation {·}iM=1 indicates a sequence of elements for i = 1, 2, . . . , M, and the ai are vectors of length M such that the j-th component ai j is given by:



ai j =

0 if j = i , 1 if j = i .

(9)

Finite-volume numerical methods are commonly used to solve systems of hyperbolic partial differential equations such as those given by Eqs. (1) and (5). Many of these methods are based on upwind Godunov methods, which approximate the solution as piecewise-continuous polynomials within each computational cell. Discontinuities at cell interfaces constitute initial states for Riemann problems that are locally solved to obtain interface fluxes needed to evolve the solutions in space and time. Embid and Baer first analyzed the eigenstructure of the BN equations [10], and the solution of the associated two-phase Riemann problem has been investigated by Andrianov and Warnecke [1] and Deledicque and Papalexandris [8], with Schwendeman et al. [25] providing the first direct solution for ideal and stiffened equations of state. This solution is also utilized in [25] to evaluate the nonconservative products analytically. However, there is currently no available solution to the Riemann problem for Eq. (5) for M > 1, primarily due to the complexity of possible wave configurations and difficulties that arise in determining their spatial order during the solution procedure. An analysis of this solution is a critical component in the construction of upwind finite-volume numerical schemes needed to predict flows governed by Eq. (5). In this paper, a direct solution methodology for the Riemann problem for an arbitrary number of solid phases is formulated that can be used with any convex EOS. This solution can be used in the verification of numerical schemes applied to multiphase flow problems and may provide a framework for the implementation of exact or approximate Riemann solvers for Eq. (5) within existing upwind numerical schemes. The solution approach is based on the methodology used in [25], and utilizes an inner Newton method to solve for the phase shock speeds and an outer Newton iteration to obtain the specific volumes. All other flow properties may be calculated from these quantities. While it is possible to eliminate the shock speeds and solve the system of nonlinear equations with a single Newton iteration for the phase pressures, the resulting relations include square-root terms which can make the system discontinuous in large sub-domains of R W , where W is

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

147

the total number of nonlinear algebraic equations. This may negatively affect the global convergence of Newton’s method, particularly for large W . In addition, non-ideal EOS’s (e.g., virial, Mie–Grüneisen) that are commonly applied to gas and solid phases may render the use of phase pressures impractical as independent unknown variables, because other flow properties are difficult to express as functions of the pressures in closed form. For these reasons the solution procedure presented here may generally provide a more robust multiphase Riemann solver than alternative methods. Because the wave configurations of the N-phase Riemann problem are complex, a system of checks must also be utilized in the solution algorithm to ensure that physically admissible configurations are selected. The solution must satisfy the solid contact wave jump conditions defined by the Riemann invariants [10,25], in addition to the classical relations for shock and expansion waves. Consequently, the specific set of nonlinear equations to be solved depends on the wave configuration, which must be detected adaptively within the Newton iteration. The practical implementation of wave ordering detection criteria is a key consideration in developing a robust solver. The plan of the paper is as follows. The Riemann solver is derived in Section 2. The two-phase solution is first considered to demonstrate the procedure used to solve for properties within both the nonlinear fields and the additional intermediate state introduced by the second phase. This problem also illustrates the dependence of the governing set of algebraic equations on the wave ordering, and motivates the use of adaptive detection criteria. The full solution for an arbitrary number of phases is then formulated on this foundation. Several test problems are then presented in Section 3 to demonstrate the accuracy of the Riemann problem solution compared to the results of a numerical scheme, and conclusions are discussed in Section 4. 2. Riemann problem solution procedure The two-phase Riemann problem solution for the BN system is first considered in this section, followed by the corresponding N-phase solution procedure. The strategy used here deviates somewhat from that presented in [25]. The first distinction concerns the choice of parameterization for relations obtained from the Riemann invariants and Rankine– Hugoniot (RH) relations across expansion waves and shocks, respectively. The two-phase solver devised by Schwendeman et al. in [25] is based in part on the single-phase solver formulated by Toro [27]. The single-phase Riemann problem solution to the Euler equations is a canonical problem in gas dynamics, and a detailed description of the solution procedure is omitted here for brevity. In this solver the local velocity u and density ρ are expressed as functions of pressure P for specific equations of state. Because the ideal EOS was considered in [27], and both the ideal and stiffened EOS’s were used in [25], closed-form expressions for u ( P ) and ρ ( P ) may be obtained in a relatively straightforward manner due to the simple functional form of these equations of state. The single-phase solution is therefore obtained completely by solving the velocity Riemann invariant for the Euler equations across the contact wave for P . However, for arbitrary convex equations of state, it is generally not possible to obtain such closed-form expressions; therefore, the use of an iterative technique is required to compute ρ and u. As a result, additional unknown variables and equations must be introduced to obtain the single-phase solution. For general EOS’s the single-phase solution may then be computed in a number of ways. In one approach outlined by Moiseev and Mukhamadieva in [19], the two Riemann problem intermediate state densities ρβ , where β = 1, 2, together with the single intermediate pressure P ∗ , may be obtained simultaneously by solving the velocity Riemann invariant across the contact wave together with the energy RH relation and/or an isentrope ordinary differential equation (ODE). The system of equations F(x) = 0 may then be solved for x using Newton’s method. However, this strategy increases the number of equations to three; in addition, the RH momentum equation expressed in terms of ρ and P introduces square-root terms which render F(x) discontinuous in large subsets of x ∈ R3 . These factors may degrade the global convergence behavior of Newton’s method for poor initial estimates of x and increase computational expense; the solution procedure will fail altogether if Newton’s method produces negative corrective specific volume iterations. To mitigate these potential issues, the left and right shock speeds S L and S R are used here to parameterize u and P when shocks are present. This leads to the second major departure from the solution strategy used in [27,25], which involves the use of an inner Newton iteration to solve the energy RH relation for S α , α = L , R, using fixed values of the specific volumes v β = ρβ−1 . A system of two equations, composed of the velocity and pressure Riemann invariants across the contact wave, is then solved in an outer Newton iteration for v β , β = 1, 2. The outer system, parameterized by S α and v β , is continuous on R2 except on manifolds where v β = 0, and numerical experiments indicate that the use of the additional inner iteration does not significantly increase computational expense. The two-phase solution procedure of Schwendeman et al. extends the concepts developed by Toro by parameterizing all phase quantities using the gas pressures P β, g , solid pressures P β,s , and a separate intermediate gas density ρ0 . However, ideal and stiffened equations of state are again considered in [25] exclusively, and a system of nonlinear equations extended to higher state space dimensions in a similar manner to the single-phase case would be necessary for general equations of state. It should be noted that there is some flexibility in the choice of parameterization for more complex equations of state, such as the virial EOS [26], but the isentrope ODE’s and/or energy RH relations must admit closed-form expressions for the phase variables in terms of the pressures to retain the advantage of this approach. In the present work separate inner Newton iterations are used to solve for shock speeds in the gas and solid phases, which are used together with the specific volumes in each phase to parameterize a system of equations solved in an outer Newton iteration. This solution procedure is discussed in detail in the following subsection; the remaining differences in methodology between the present strategy

148

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

and that in [25] are minor, and arise from the order in which the governing algebraic–differential equations are solved. Due to the convergence considerations discussed for the single-phase case, the solution approach presented here is considered to be a prudent choice for arbitrary equations of state applied to the multiphase case. The basic algorithmic approach for the two-phase system is extended in a relatively straightforward manner to the Riemann problem associated with the general N-phase system described by Eq. (5). It is also noted that expansion waves may require numerical integration of the isentrope ODE’s depending on the EOS used; however, conventional integration schemes may be used to obtain high accuracy in these cases. 2.1. Two-phase problem In the Riemann problem associated with the BN equations, the solid-phase wave configurations are identical to those from the single-phase theory, with two nonlinear waves (shocks or expansions) and a contact wave. The solid and gas are coupled through the compaction equation; therefore, gas-phase properties experience a discontinuous jump across the compaction (i.e., solid contact) wave. These conditions are given by:

[u s ]+ − = 0, [η g ]+ − = 0, +  φ g (u s − u g ) vg



−

φ g P g + φs P s +  eg + P g v g +

1 2

= 0, φ g (u s − u g )2 vg

(u s − u g )2

+ −

+ −

= 0,

= 0,

(10)

where η is entropy; v = ρ −1 , and [·]+ − indicates a jump across the solid contact discontinuity from the left “−” state to the right “+” state. Unlike the single-phase solution, there are several possible wave configurations for the gas-phase properties, as shown in Fig. 1, determined by the relative positions of the solid contact and nonlinear gas waves. Here, configurations for which |u s − u g | < c g are referred to as “subsonic” in the literature, and those corresponding to |u s − u g | > c g are denoted as “supersonic”. The special cases of coincident waves corresponding to either u s = u g , or a parabolic degeneracy where |u s − u g | = c g , have been analyzed in the literature [8,13], but are rare events and are not considered here. In the two-phase problem, each state is defined by v s , u s , P s , v g , u g , P g , and φs . Since φs and φ g change only across the solid contact, we note that only the “L” and “R” subscripts associated with the initial states are necessary. For all other quantities, the states in Fig. 1 are denoted by the same convention used for the single-phase case, with an additional intermediate “0” gas state. Here, cases for which the solid contact is positioned to the left of the gas contact are denoted by Configuration A, as shown in Figs. 1(a) and 1(c). Cases where the solid contact is located to the right of the gas contact are referred to as Configuration B, which are indicated by Figs. 1(b) and 1(d). As noted in [1,25], the Riemann solver must include a method to determine which of the possible wave configurations is valid for given initial left and right states. However, although the subsonic configurations are more physically relevant in actual multiphase flows due to large interphase drag, it is generally not sufficient to consider only subsonic cases during the solution iteration [7]. The states adjacent to the nonlinear waves in the solid phase are solved in an identical manner to the single-phase case, where the RH relations and Riemann invariants remain valid across shocks and expansions, respectively:



P β,s =

u β,s =



P α ,s + 1 −

v β,s (u α ,s − S α ,s )2 , v α ,s v α ,s

P η,s ( v β,s ; v α ,s ), v β,s (u α , s − S α , s ) + v α ,s

v c (v ) u α ,s ∓ v β,s − s v s α ,s s

v β,s < v α ,s , v β,s  v α ,s ,

S α ,s ,

v β,s < v α ,s ,

dv s ,

v β,s  v α ,s .

(11)

(12)

For all equations presented in this work, β = 1 when the left initial state α = L is used; β = 2 when the right initial state α = R is used, and S denotes the shock speed. Here, the function P η,s ( v ) is the solution of the isentrope ODE for the solid phase:

dP dv

2  c(v , P ) =− ,

(13)

v

where the sound speed c is determined using the EOS:

 c2 (v , P ) = v 2

P + ∂∂ ve | P ∂e ∂ P |v

 ;

(14)

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

149

Fig. 1. Possible gas-phase wave ordering for the two-phase Riemann problem for (a)–(b) subsonic and (c)–(d) supersonic configurations.

a similar function P η, g is applied to the gas phase. Calculation of the intermediate “0” gas state depends on the wave configurations shown in Fig. 1. Therefore, Eq. (10) gives:

u0 =

⎧ ⎪ ⎪ u 2, s + ⎪ ⎪ ⎪ ⎪ ⎨u + 2, s

⎪ ⎪ u 1, s + ⎪ ⎪ ⎪ ⎪ ⎩u +

φL, g φR,g φL, g φR,g φR,g φL, g φR,g φL, g

v0 v L, g

(u L , g − u 1,s ),

Supersonic A,

v0 v 1, g

(u 1, g − u 1,s ),

Subsonic A,

v0 v R,g

(u R , g − u 2,s ), Supersonic B,

(15)

v0 v 2, g

(u 2, g − u 2,s ), Subsonic B, ⎧  φ L , g (u L , g −u 1,s )2 φ ( u − u )2  1 ⎪ − φ R ,s P 2,s − R , g v0 0 2,s , ⎪ v L, g ⎪ φ R , g φ L , g P L , g + φ L , s P 1, s + ⎪ ⎪ ⎪ 1  φ L , g (u 1, g −u 1,s )2 φ ( u − u )2  ⎪ ⎨ − φ R ,s P 2,s − R , g v0 0 2,s , φ R , g φ L , g P 1, g + φ L , s P 1, s + v 1, g P0 =  φ R , g (u R , g −u 2,s )2 φ ( u − u )2  1 ⎪ ⎪ − φ L ,s P 1,s − L, g 0v 0 1,s , ⎪ φ L , g φ R , g P R , g + φ R , s P 2, s + v R,g ⎪ ⎪ ⎪ ⎪ ⎩ 1 φ P + φ P + φ R , g (u 2, g −u 2,s )2 − φ P − φL, g (u 0 −u 1,s )2 , 1, s

φL, g

R,g

2, g

R ,s

2, s

v 2, g

L ,s

1, s

⎧ e L , g + P L , g v L , g + 12 (u L , g − u 1,s )2 − P 0 v 0 − 12 (u 0 − u 2,s )2 , ⎪ ⎪ ⎪ ⎪ ⎨ e + P v + 1 (u − u )2 − P v − 1 (u − u )2 , 1, g 1, g 1, g 1, g 1, s 0 0 0 2, s 2 2 e0 = 1 1 2 ⎪ e R , g + P R , g v R , g + 2 (u R , g − u 2,s ) − P 0 v 0 − 2 (u 0 − u 1,s )2 , ⎪ ⎪ ⎪ ⎩ e 2, g + P 2, g v 2, g + 12 (u 2, g − u 2,s )2 − P 0 v 0 − 12 (u 0 − u 1,s )2 ,

v0

Supersonic A, Subsonic A,

(16)

Supersonic B, Subsonic B,

Supersonic A, Subsonic A, Supersonic B,

(17)

Subsonic B.

Thus, the gas-phase solution across the nonlinear waves also depends on the particular wave configuration. The supersonic cases require special treatment of quantities adjacent to the “0” state; all other quantities for both subsonic and supersonic configurations are determined in the same manner as the single-phase case:

150

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Supersonic A:

P 1, g =

⎧ ⎨



P0 + 1 −

v 1, g v0

( u L , g − S L , g )2

P η, g ( v 1, g ; v 0 ),

v0

, v 1, g < v 0 , v 1, g  v 0 ,

(18)

v 1, g v0

( u 0 − S L , g ) + S L , g , v 1, g < v 0 , u 1, g =

⎩ u 0 − v g ,1 − c g ( v g ) dv g , v 1, g  v 0 , v0 vg

 v 2, g (u 0 − S R , g )2 P0 + 1 − v , v 2, g < v 0 , v0 0 Supersonic B: P 2, g = P η, g ( v 2, g ; v 0 ), v 2, g  v 0 , ⎧v ⎨ v2, g (u 0 − S R , g ) + S R , g , v 2, g < v 0 , 0 u 2, g =

⎩ u 0 + v 2, g − c g ( v g ) dv g , v 2, g  v 0 , v0 vg

 v ( u α − S α , g )2 P α , g + 1 − v β, g , v β, g < v α , g , vα α, g Else: P β, g = P ( v ; v ), v β, g  v α , g , ⎧ v η, g β, g α , g ⎨ v β, g (u α , g − S α , g ) + S α , g , v β, g < v α , g , α, g u β, g =

⎩ u α , g ∓ v β, g − c g ( v g ) dv g , v β, g  v α , g . v α, g vg

(19)

(20)

Thus, all quantities are expressed in terms of the specific volumes and shock speed(s). If shocks are present, the energy RH equation is used to express the shock speeds in terms of the specific volumes:

e β,s ( v β,s , P β,s ) = e α ,s ( v α ,s , P α ,s ) + P α ,s ( v α ,s − v β,s ) +

1



2

v β,s v α ,s

2 (u α ,s − S α ,s )2 ,

−1

v β,s < v α ,s .

(21)

The gas-phase relations again depend on the configuration considered. Special treatment is required for regions adjacent to the “0” state in supersonic cases; all other flow states for both subsonic and supersonic configurations are determined in a manner identical to Eq. (21):

Supersonic A: e 1, g ( v 1, g , P 1, g ) = e 0 + P 0 ( v 0 − v 1, g ) + Supersonic B: e 2, g ( v 2, g , P 2, g ) = e 0 + P 0 ( v 0 − v 2, g ) +

1 2 1 2

 

v 1, g v0 v 2, g v0

2 −1

v 1, g < v 0 ,

(22)

(u 0 − S R , g )2 ,

v 2, g < v 0 ,

(23)

2 −1

Else: e β, g ( v β, g , P β, g ) = e α , g ( v α , g , P α , g ) + P α , g ( v α , g − v β, g ) + v β, g < v α , g .

(u 0 − S L , g )2 ,

1 2



v β, g v α, g

2 −1

(u α , g − S α , g )2 , (24)

Because all flow states are expressed solely in terms of the specific volumes, the system of equations to be solved consists of five relations in the variables v 1,s , v 2,s , v 1, g , v 2, g , and v 0 . Two equations are the constant-pressure and constant-velocity conditions across the gas contact wave obtained from Eq. (10); one is a constant-velocity condition across the solid contact; one is an EOS equality across the solid contact, and the final equation enforces constant gas entropy across the solid contact.

⎧ ⎨ P 2, g = P 0 , P 1, g = P 0 , ⎩P = P , 1, g 2, g ⎧ ⎨ u 2, g = u 0 , u 1, g = u 0 , ⎩u = u , 1, g 2, g

Subsonic A, Subsonic B, Supersonic A or B,

(25)

Subsonic A, Subsonic B, Supersonic A and B,

(26)

u 1, s = u 2, s ,

(27)

e( v 0 , P 0 ) = e0 ,

(28)

⎧ P η, g ( v 0 ; v L ) = P 0 , ⎪ ⎨ P η, g ( v 0 ; v R ) = P 0 , ⎪ ⎩ P η, g ( v 0 ; v 1 ) = P 0 , P η, g ( v 0 ; v 2 ) = P 0 ,

Supersonic A, Supersonic B, Subsonic A, Subsonic B.

(29)

The final expression in Eqs. (25)–(29) is obtained by exploiting Eq. (13) to ensure isentropic gas flow, rather than attempting an explicit calculation of the gas-phase entropy.

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

151

The following procedure is then used to obtain the solution: (0)

(0)

(0)

(0)

(0)

1. Initial estimates for v 1,s , v 2,s , v 1, g , v 2, g , and v 0 are made, which are used to compute velocities and sound speeds in regions “1” and “2”. Schwendeman et al. [25] used the single-phase Riemann problem solutions as initial guesses, since these provide the solution to the two-phase problem in the limit as |φ L ,s − φ R ,s | → 0 where the phases are decoupled. (0)

(0)

(0)

(0)

(0)

A separate value of v 0 must also be provided. Initial values of the shock speeds, denoted by S L , g , S R , g , S L ,s , and S R ,s are then obtained, if necessary. These guesses are used to estimate the initial wave configuration: • Supersonic A:



( 0)

( 0)

( 0)

u 1, s < S L , g

if v 1, g < v L , g ,

( 0)

( 0)

u 1, s < u L , g − c L , g

if v 1, g  v L , g ,

• Supersonic B:



( 0)

( 0)

( 0)

( 0)

( 0)

u 2, s > S R , g

if v 2, g < v R , g , ( 0)

( 0)

u 2, s > u R , g + c R , g

if v 2, g > v R , g ,

(0)

(0)

(0)

(0)

(0)

(0) < u 2, g

(0)

• Subsonic A: u 1, g − c 1, g < u s < u g , • Subsonic B: u g < u s

(0) + c 2, g .

These adaptive wave configuration detection criteria were not considered by Schwendeman et al. in [25], since either a subsonic or supersonic configuration was assumed a priori. (m) (m) (m) (m) 2. After the m-th iteration, the specific volumes v 1, g , v 2, g , v 1,s , and v 2,s are known. For m = 0, the following wave configuration detection procedure may be ignored. Otherwise, the solid velocity, gas velocity, and gas sound speeds from the (m − 1)-th iteration are used to determine the wave configuration adaptively using criteria similar to that given in step 1: • Supersonic A:



(m−1)

u 1, s

(m−1) u 1, s

(m−1)

(m−1)

< S L, g

(m−1)

< u0

if v 1, g (m−1)

− c0

(m−1) if v 1, g

(m−1)

< v0

(m−1)

 v0

, ,

• Supersonic B:



(m−1)

> S (Rm, g−1)

(m−1)

> u (0m−1) + c 0(m−1) if v (2m, g−1) > v (0m−1) ,

u 2, s u 2, s

(m−1)

• Subsonic A: u 1, g

(m−1)

(m−1)

− c 1, g

(m−1)

if v 2, g

(m−1)

< us

< v (0m−1) ,

(m−1)

< ug

,

(m−1)

(m−1) (m−1) • Subsonic B: u g < us < u 2, g + c 2, g . (m) (m) 3. The solid shock speeds S L ,s and S R ,s are computed, if necessary, from Eq. (21) using an inner Bisection method or (m) (m) (m) (m) Newton iteration. The solid-phase quantities u 1,s , P 1,s , u 2,s , and P 2,s are then calculated from Eqs. (11) and (12).

4. For supersonic states: (m) (m) (m) (a) The “0” state properties u 0 , P 0 , and e 0 are first computed using Eqs. (15)–(17). (b) If necessary, the shock speeds in both phases are computed using separate inner Bisection method or Newton method iterations for Eqs. (22) and (24). (m) (m) (m) (m) (c) The remaining gas-phase quantities u 1, g , P 1, g , u 2, g , and P 2, g are then determined from Eqs. (18)–(20). For subsonic states: (m) (m) (a) If necessary, the shock speeds S L , g and S L ,s are computed using separate inner Bisection method and Newton method iterations for Eqs. (22) and (24). (m) (m) (m) (m) (b) The gas-phase quantities u 1,s , P 1,s , u 2,s , and P 2,s are then determined from Eqs. (18)–(20). (m)

(m)

(m)

(c) The remaining “0” state quantities u 0 , P 0 , and e 0

are computed. (m+1)

5. The system of Eqs. (25)–(29) is solved using a Newton iteration to obtain the updated quantities v 1,s (m+1)

(m+1)

(m+1)

, v 2, s

(m+1)

, v 1, g

,

. v 2, g , and v 0 6. Steps 2–5 are repeated until convergence. Some additional comments on the selection of initial estimates are in order, since this choice can significantly affect the convergence behavior of the iterative method. The single-phase solutions obtained by using the appropriate initial data (0) (0) for each phase and ignoring the volume fraction, as outlined in [25], are a logical choice for the selection of v 1,s , v 2,s ,

152

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Fig. 2. General gas-phase wave diagram for a system with M solid phases.

Table 1 Two-phase wave configuration nomenclature.

Supersonic A Subsonic A Supersonic B Subsonic B

(0)

l

b

r

1 0 0 0

0 1 0 0

2 2 1 2

(0)

(0)

v 1, g , and v 2, g . The difficulty arises in choosing a suitably accurate estimate for v 0 . A heuristic selection is the “nearest neighbor” gas state immediately adjacent to the solid contact wave bordering the “0” region. A linearization of the system equations based on a Taylor expansion about this “nearest neighbor” state provides a more accurate estimate [26], but is dependent on the equation of state, where an explicit closed-form solution may not be generally available. Similar issues arise in choosing initial estimates for the general N-phase system; therefore, the formulation of alternative selection criteria suitable for general equations of state remains an ongoing topic of research. 2.2. Solution for multiple solid phases The procedure for constructing a Riemann solver for M solid phases is effectively an extension of the two-phase methodology. For each solid phase the wave configurations remain unchanged from the two-phase case; however, the gas phase now includes M intermediate “0” states. The wave ordering is now significantly more complex, since there may exist a combination of subsonic and supersonic wave configurations, as shown in Fig. 2. Here, three special intermediate regions are identified, because they allow for a systematic classification of wave configurations in constructing the solution. The supersonic intermediate states immediately adjacent to the left and right nonlinear waves are designated by the subscripts l and r, respectively; and the intermediate subsonic field located to the immediate left of the gas contact is denoted by the subscript b. The intermediate states are ordered from left to right, with subscripts increasing from 1 to M. For the special case where a left supersonic configuration does not exist, l = 0, which corresponds to the left initial state L. Similarly, when a right supersonic wave configuration does not exist, r = M + 1, corresponding to the right initial state R. The subscript convention used for phase properties in this problem is now defined. Quantities with the subscript 0, j refer to gas-phase properties in the j-th intermediate “0” region, for j = 1, 2, . . . , M, as shown in Fig. 2. The remaining gas fields corresponding to states “1” and “2” in Fig. 2 are denoted by the subscripts 1, g and 2, g, respectively. For the j-th solid phase only the single-phase “1” and “2” states exist, and the associated quantities are denoted by the subscripts 1, j and 2, j, respectively, for j = 1, 2, . . . , M. To illustrate the use of this notation, Table 1 gives the values of l, r, and b corresponding to the configurations shown in Fig. 1. Here, the solution procedure consists of a series of recursive calculations from the left and right initial states to obtain relations in the intermediate regions adjacent to the gas contact wave, denoted by b and b + 1 in Fig. 2. From the numbering convention defined here, the gas wave solution consists of l supersonic states adjacent to the left nonlinear gas wave, and M − (r − 1) supersonic states adjacent to the right nonlinear gas wave. Since there are a total of M intermediate “0” states, there are r − l − 1 subsonic states, in addition to the “1” and “2” regions. Thus, the speeds of the solid contacts must be computed and arranged such that u 1,1 < u 1,2 < · · · < u 1, M . Since the solid volume fractions are initially specified, the gas volume fractions φ0, j are also known, and are given by:

φ 0, j =

1−

j

φ k=1 R ,k  j −1 1 − k=1 φ R ,k

− −

M

k= j +1 φ L ,k ,

M

k= j

φ L ,k ,

j = 1, 2, . . . , b , j = b + 1, b + 2, . . . , M ,

(30)

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

153

Fig. 3. Type I gas-phase wave configuration cases for (a) l = 0, (b) r = M + 1, and (c) l = r − 1 > 0.

where φ L ,k and φ R ,k are the initial left and right solid volume fractions of the k-th solid phase. As in the two-phase scenario, linear degenerate cases are not considered here for simplicity. 2.2.1. Special cases Although Fig. 2 illustrates the general combinations of gas-phase wave configurations, several specific cases must also be considered in constructing the Riemann solver to satisfy jump conditions across the gas contact wave. Type I is a scenario in which only supersonic intermediate regions are present, where l = r − 1. This scenario is shown in Fig. 3, and can be summarized as:



Type I: l = r − 1

⇒

l = M, l = 0, 0 < l < M,

all left supersonic waves, all right supersonic waves, left and right supersonic waves.

(31)

Type II corresponds to the case where b = l, and no subsonic solid contact waves are present between the left nonlinear gas wave and gas contact wave. Similarly, Type III configurations arise when b = r − 1, and no subsonic solid contacts exist between the gas contact and the right nonlinear gas wave. Schematics of the gas-phase wave configurations for Types II and III are shown in Figs. 4 and 5, respectively. Here, it should be noted that these are not necessarily mutually exclusive designations, and configurations belonging to Types II or III may also be classified as Type I for the case where b = l = r − 1. 2.2.2. Riemann solver procedure (0) (0) 1. For j = 1, 2, . . . , M, the j-th solid specific volume initial values v 1, j and v 2, j are obtained using initial estimates, such as those obtained from the corresponding single-phase solutions for each phase [25], and the initial shock speed values (0) (0) S L , j and S R , j , are computed if necessary. Initial values for the intermediate gas specific volumes v 0, j must also be (0)

(0)

(0)

(0)

provided. The gas specific volumes v 1, g , v 2, g and shock speeds S L , g , and S R , g , are also computed from the single-phase analysis. The configuration parameters l, b and r are initially determined using the gas shock speed, sound speed, and local velocity estimates: • Region l:



( 0)

( 0)

( 0)

( 0)

( 0)

u 0,l−1 < u 0,l < S L , g u 0,l−1 < u 0,l < u L , g − c L , g

( 0)

if v 1, g < v L , g , ( 0)

if v 1, g  v L , g ,

154

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Fig. 4. Type II gas-phase wave configuration (b = l).

Fig. 5. Type III gas-phase wave configuration (b = r − 1).

• Region r:



( 0)

( 0)

( 0)

( 0)

( 0)

( 0)

u 0,r +1 > u 0,r > S R , g

if v 2, g < v R , g ,

u 0,r +1 > u 0,r > u R , g + c R , g (0)

(0)

( 0)

if v 2, g  v R , g ,

(0)

• Region b: u 0,b < u 1, g < u 0,b+1 .

(m)

(m)

(m)

(m)

(m)

2. After the m-th iteration, the specific volumes v 1, g , v 2, g , v 1, j , v 2, j , and v 0, j are known. For m = 0, the following wave configuration detection procedure may be ignored. Otherwise, the solid contact waves are arranged such that (m−1) (m−1) (m−1) u 1, 1 < u 1,2 < · · · < u 1, M . Using the gas velocity and gas sound speeds obtained from the (m − 1)-th iteration, the wave configuration for the m-th iteration is determined adaptively by updating the values of l, b, and r, using criteria similar to that outlined in step 1: • Region l:



(m−1)

(m−1)

(m−1) u 0,l−1

(m−1) < u 0,l

u 0,l−1 < u 0,l

(m−1)

(m−1)

< S L, g

(m−1) < u 0,l

if v 1, g (m−1) − c 0,l

(m−1) if v 1, g

(m−1)

< v 0,l 

,

(m−1) v 0,l ,

• Region r:



(m−1)

(m−1)

(m−1) u 0,r +1

(m−1) > u 0,r

u 0,r +1 > u 0,r

(m−1)

• Region b: u 0,b

(m−1)

(m−1)

> S R,g

(m−1) > u 0,r

(m−1)

< u 1, g

if v 2, g (m−1) + c 0,r

(m−1) if v 2, g

(m−1)

< v 0,r 

,

(m−1) v 0,r ,

(m−1)

< u 0,b+1 .

3. The m-th iterations of the remaining flow properties are now computed. The supersonic wave states are calculated first. After dropping the (m) superscript for convenience, applying the solid contact jump conditions to these regions yields:

u 0, j =

⎧ ⎨ u 2, j + ⎩ u 1, j +

φ0, j−1 φ0, j φ0, j+1 φ0, j

v 0, j ( u 0, j −1 v 0, j −1 v 0, j ( u 0, j +1 v 0, j +1

− u 1, j ),

j = 1, 2, . . . , l ,

− u 2, j ),

j = M , M − 1, . . . , r ,

(32)

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

⎧ ⎪ ⎪ ⎪ ⎪ ⎨

P 0, j =

⎪ ⎪ ⎪ ⎪ ⎩

1

φ0, j

 φ 0, j −1 P 0, j −1 + φ L , j P 1 , j +

j = 1, 2, . . . , l , 1

φ0, j

φ0, j−1 v 0, j −1

 φ 0, j +1 P 0, j +1 + φ R , j P 2 , j +

φ0, j+1 v 0, j +1

(u 0, j −1 − u 1, j )2 − φ R , j P 2, j −

φ0, j

(u 0, j +1 − u 2, j )2 − φ L , j P 1, j −

φ0, j

v 0, j

v 0, j

155

 (u 0, j − u 2, j )2 ,  (u 0, j − u 1, j )2 ,

j = M , M − 1, . . . , r , ⎧ e 0, j −1 + P 0, j −1 v 0, j −1 + 12 (u 0, j −1 − u 1, j )2 − P 0, j v 0, j − 12 (u 0, j − u 2, j )2 , ⎪ ⎪ ⎪ ⎨ j = 1, 2, . . . , l , e 0, j = ⎪ e 0, j +1 + P 0, j +1 v 0, j +1 + 12 (u 0, j +1 − u 2, j )2 − P 0, j v 0, j − 12 (u 0, j − u 1, j )2 , ⎪ ⎪ ⎩ j = M , M − 1, . . . , r .

(33)

(34)

4. The nonlinear wave fields are computed:

P 1, g =

P 2, g =



P 0,l + 1 −

v g ,1 (u 0,l − S L , g )2 , v 0,l v 0,l

v 1, g < v 0,l ,

P η, g ( v 1, g ; v 0,l ),



(u 0,r − S R , g )2

v 2, g v 0,r

P 0,r + 1 −

(35)

v 1, g  v 0,l , v 0,r

, v 2, g < v 0,r ,

P η ( v 2, g ; v 0,r ),

(36)

v 2, g  v 0,r .

Here, u 1, g = u 0,l − f L , g and u 2, g = u 0,r + f R , g , where v 1, g (u 0,l v

v 1, g 0,cl − v dv , v 0,l

 f L, g =

1−



− S L , g ), v 1, g < v 0,l ,

v 2, g (u 0,r v

v 2 0c,r − v dv , v 0,r

f R,g =

(37)

v 1, g  v 0,l ,

1−

− S R , g ), v 2, g < v 0,r ,

(38)

v 2, g  v 0,r .

Similar expressions for the j-th solid phase are also obtained for j = 1, 2, . . . , M:

P 1, j =

P 2, j =



P L, j + 1 −

v 1, j v L, j

( u L , j − S L , j )2 v L, j

, v 1, j < v L , j ,

P η, j ( v 1, j ; v L , j ),



P R, j + 1 −

v 2, j v R, j

(39)

v 1, j  v L , j ,

( u R , j − S R , j )2 v R, j

, v 2, j < v R , j ,

P η, j ( v 2, j ; v R , j ),

(40)

v 2, j  v R , j .

The velocities u 1, j = u L , j − f L , j and u 2, j = u R , j + f R , j , and

 f L, j =

v 1, j (u L , j v

v 1, j L,cj − v dv , v L, j

1−

 f R, j =

v 2, j v

v 2, j R ,cj −v v R, j

1−



− S L , j ), v 1, j < v L , j ,

(41)

v 1, j  v L , j ,

(u R , j − S R , j ), v 2, j < v R , j ,

dv ,

(42)

v 2, j  v R , j . (m)

(m)

(m)

(m)

(m)

(m)

(m)

(m)

(m)

5. Taking the m-th outer Newton iterations v 1, g , v 2, g , v 0, j , v 1, j , v 2, j , the shock speeds S L , g , S R , g , S L , j , and S R , j are obtained from an inner Newton loop if necessary. Dropping the (m) superscript for convenience:

e g ( v 1, g , P 1, g ) = e g ( v 0,l , P 0,l ) + P 0,l ( v 0,l − v 1, g ) +

1



2

e g ( v 2, g , P 2, g ) = e g ( v 0,r , P 0,r ) + P 0,r ( v 0,r − v 2, g ) +

v 0,l

1



2

e j ( v β, j , P β, j ) = e j ( v α , j , P α , j ) + P α , j ( v α , j − v β, j ) +

v 1, g

1 2

2 −1

v 2, g v 0,r



(u 0,l − S L , g ),

(43)

2 −1

v β, j v α, j

(u 0,r − S R , g ),

(44)

2 −1

(u α , j − S α , j ).

(45)

6. The subsonic wave field properties are now computed using the solid contact jump relations. For the gas state adjacent to the “1” region,

156

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

φ0,l v 0,l+1 (u 1, g − u 1,l+1 ), φ0,l+1 v 1, g

u 0,l+1 = u 2,l+1 + P 0,l+1 =



1

φ0,l+1 −

φ0,l P 1, g + φ L ,l+1 P 1,l+1 +

φ0,l+1 v 0,l+1

e 1, g + P 1, g v 1, g +

(u 0,l+1 − u 2,l+1 )2 , 1 2

φ0,l v 1, g

(46)

(u 1, g − u 1,l+1 )2 − φ R ,l+1 P 2,l+1 (47) 1

(u 1, g − u 1,l+1 )2 − P 0,l+1 v 0,l+1 − (u 0,l+1 − u 2,l+1 )2 ;

(48)

2

and for the state adjacent to the “2” state,

u 0,r −1 = u 1,r −1 + P 0,r −1 =



1

φ0,r P 2, g + φ R ,r −1 P 2,r −1 +

φ0,r −1 −

φ0,r v 0,r −1 (u 2, g − u 2,r −1 ), φ0,r −1 v 2, g

φ0,r −1 v 0,r −1

(u 0,r −1 − u 1,r −1 )2 ,

e 0,r −1 = e 2, g + P 2, g v 2, g +

1 2

φ0,r v 2, g

(49)

(u 2, g − u 2,r −1 )2 − φ L ,r −1 P 1,r −1 (50) 1

(u 2, g − u 2,r −1 )2 − P 0,r −1 v 0,r −1 − (u 0,r −1 − u 1,r −1 )2 . 2

(51)

The remaining intermediate region properties are now computed from the solid contact jump relations:

u 0, j =

⎧ ⎨ u 2, j + ⎩ u 1, j +

φ0, j−1 φ0, j φ0, j+1 φ0, j

v 0, j ( u 0, j −1 v 0, j −1 v 0, j v 0, j +1

− u 1, j ),

(u 0, j +1 − u 2, j ),

j = l + 2, l + 3, . . . , b ,

(52)

j = r − 2, r − 3, . . . , b + 1,

⎧   φ0, j−1 φ0, j 1 2 2 ⎪ ⎪ φ0, j φ0, j −1 P 0, j −1 + φ L , j P 1, j + v 0, j−1 (u 0, j −1 − u 1, j ) − φ R , j P 2, j − v 0, j (u 0, j − u 2, j ) , ⎪ ⎪ ⎨ j = l + 2, l + 3, . . . , b , P 0, j =   φ φ 1 ⎪ ⎪ φ P + φ R , j P 2, j + v 00,, jj++11 (u 0, j +1 − u 2, j )2 − φ L , j P 1, j − v 00,, jj (u 0, j − u 1, j )2 , ⎪ ⎪ φ0, j 0, j +1 0, j +1 ⎩ j = r − 2, r − 3, . . . , b + 1, ⎧ e 0, j −1 + P 0, j −1 v 0, j −1 + 12 (u 0, j −1 − u 1, j )2 − P 0, j v 0, j − 12 (u 0, j − u 2, j )2 , ⎪ ⎪ ⎨ j = l + 2, l + 3, . . . , b , e 0, j = 1 1 2 2 ⎪ ⎪ ⎩ e 0, j +1 + P 0, j +1 v 0, j +1 + 2 ( u 0, j +1 − u 2 , j ) − P 0, j v 0, j − 2 ( u 0, j − u 1 , j ) , j = r − 2, r − 3, . . . , b + 1. (m+1)

(m+1)

(m+1)

(m+1)

(53)

(54)

(m+1)

7. The outer Newton iteration is used to solve for v 1 , v2 , v 1, j , v 2, j , and v 0, j . The system of 2 + 3M equations is given by the EOS and Riemann invariants across the gas and solid contact waves. Because gas velocity and pressure must be preserved across the gas contact wave, the special cases denoted by Types I–III must also be considered here:

⎧ P 1, g = P 2, g , ⎪ ⎨ P 1, g = P 0,b+1 , ⎪ ⎩ P 0,b = P 2, g , P 0,b = P 0,b+1 , ⎧ u 1, g = u 2, g , ⎪ ⎨ u 1, g = u 0,b+1 , ⎪ ⎩ u 0,b = u 2, g , u 0,b = u 0,b+1 , u 1, j = u 2, j ,

Type I, Type II, Type III, otherwise, Type I, Type II, Type III, otherwise,

j = 1, 2, . . . , M ,

e ( v 0, j , P 0, j ) = e 0, j ,



P 0, j =

j = 1, 2, . . . , M ,

P η, g ( v 0, j ; v 0, j −1 ), P η, g ( v 0, j ; v 0, j +1 ),

j = 1, 2, . . . , b , j = M , M − 1, . . . , b + 1.

8. Steps 2–7 are then repeated until convergence.

(55)

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

157

For analytically simple EOS’s, such as the ideal or stiffened EOS considered in [25], closed-form solutions for the shock speeds can be obtained from the RH energy conservation equation. For these special cases the inner Newton iterations are not required, and the shock speeds may be expressed in terms of the unknown specific volumes solved in the outer iteration. Similarly, the isentrope ODE (13) can be solved analytically to obtain explicit expressions for the specific volumes in the intermediate “0” gas states. This eliminates the need to include the gas entropy Riemann invariant in Eq. (10) within the set of nonlinear equations solved in the outer iteration. Furthermore, the velocity integrals appearing in the Riemann invariant expressions may be evaluated analytically, and do not require the additional computational expense of a numerical quadrature. (0) As in the two-phase case, suitably accurate selections of the initial estimates v 0, j for j = 1, 2, . . . , M are challenging to formulate for complex equations of state. A further complication is the linearization of the N-phase system, where it is unclear whether simple closed-form expressions for the specific volume perturbations can be obtained in an analogous manner to two-phase systems. Due to these difficulties, the simple “nearest neighbor” method described for the two-phase solution is applied to the example problems discussed in the following section. 3. Example problems To demonstrate the algorithmic approach, solutions for Riemann problems are obtained for mixtures containing three phases (M = 2). This allows for an examination of the salient features of multiple wave configurations, without introducing unnecessary wave ordering complexity associated with larger M. These are compared to numerical results obtained from a centered finite-volume method. This scheme is an extension of the method formulated by Kurganov et al. [15], which is modified to formally discretize the nonconservative sources in Eq. (5) [6]. The principal advantage of this centered scheme is that the use of a Riemann solver is not necessary. Only knowledge of the left and right initial velocities and sound speeds in each phase at each cell boundary is required. However, as discussed in [6], accuracy of the numerical scheme is reduced compared to upwind techniques. This is due to errors arising from the approximate discretization of the nonconservative sources, which are not present when Riemann solvers are utilized, as well as the increased numerical diffusion in the vicinity of contact waves common to centered schemes in general. For the cases presented here, both subsonic and supersonic wave configurations are considered. For simplicity in specifying predetermined wave configurations, the first two problems utilize the ideal gas EOS for each phase. For these problems the discussion is focused on numerical and algorithmic issues inherent to the Riemann solver, such as convergence properties and suitable initial estimates, rather than the physical implications of the selected EOS. Because the solver has been formulated for a general convex EOS, the resulting wave configurations are independent of the chosen constitutive relations. To examine the behavior of gas–granular solid mixtures using representative material properties, the third problem applies a Mie–Grüneisen EOS to both of the solid phases. For each problem a resolution of 800 computational cells is used in the numerical scheme, and the open-source Fortran program nleq1.f is used to solve the system of nonlinear algebraic equations that comprise the Riemann solver with Newton’s method. 3.1. Problem 1 Here, both the gas and “solid” phases are governed by the ideal gas EOS:

e( v , P ) =

Pv

γ −1

,

(56)

where γ is the specific heat ratio, and γ1 = γ2 = γ g = 1.4; the subscript g denotes the gas phase, and the subscripts 1 and 2 denote the corresponding “solid” phases, respectively. Initial guesses for the specific volumes in regions “1” and “2” for each phase are obtained from the single-phase theory, in an analogous manner to the two-phase case discussed in [25]. Initial estimates for the intermediate gas specific volumes v 0, j within the “0” regions are then taken to be equal to the values in the adjacent regions:

⎧ v , j  l, L, g ⎪ ⎪ ⎪ ( 0) ⎨ v ( 0) 1, g , l < j  b , v 0, j = ⎪ v ( 0) , b < j < r , ⎪ ⎪ ⎩ 2, g v R,g , j  r,

(57)

for j = 1, 2. The initial left and right states represent a multiphase shock tube problem, and are provided in Table 2. The gas- and “solid”-phase solutions for this problem are shown in Figs. 6 and 7, respectively. There are no supersonic contact wave configurations, and u 1, g < u 1,1 < u 1,2 . Here, l = 0, b = 0, and r = 3. Good agreement between the iterative and numerical solutions exists. Unlike the single-phase solution to the classical Euler equations, the “solid” velocities are not equal across its contact wave, as shown in Fig. 7.

158

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Table 2 Initial data for problem 1.

Left: Right:

ρ g (kg m−3 )

u g (m/s)

P g (MPa)

φ g (–)

3.0 1.23

0 0

0.3 0.101

0.3 0.75

ρ1 (kg m−3 )

u 1 (m/s)

P 1 (MPa)

φ1 (–)

Left: Right:

4.0 1.23

0 0

0.6 0.101

0.4 0.1

ρ2 (kg m−3 )

u 2 (m/s)

P 2 (MPa)

φ2 (–)

Left: Right:

5.0 1.23

0 0

1.0 0.101

0.3 0.15

Fig. 6. Gas-phase solution profile for problem 1 for (a) density, (b) pressure, and (c) velocity at t = 0.7 ms.

3.2. Problem 2 All phases are described by an ideal EOS, with γ1 = γ2 = γ g = 1.4. Here, the initial left and right shock tube states given in Table 3 are prescribed to produce a mixed subsonic–supersonic wave configuration. For this problem l = 0, b = 0, and r = 2. However, the Newton method failed to converge when the single-phase theory was used to obtain the initial states, as described in step 1 in Section 2.2.2. The results of the unsteady numerical method were used to select alternate initial

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

159

Fig. 7. “Solid”-phase solution profile for problem 1 for (a) phase 1 density, (b) phase 1 pressure, (c) phase 2 density, and (d) phase 2 pressure at t = 0.7 ms.

Table 3 Initial data for problem 2.

ρ g (kg m−3 )

u g (m/s)

P g (MPa)

φ g (–)

Left: Right:

3.0 1.23

0 0

0.3 0.101

0.3 0.75

ρ1 (kg m−3 )

u 1 (m/s)

P 1 (MPa)

φ1 (–)

Left: Right:

8.0 1.23

0 0

0.6 0.101

0.4 0.1

Left: Right:

ρ2 (kg m−3 )

u 2 (m/s)

P 2 (MPa)

φ2 (–)

0.5 1.23

0 0

1.0 0.101

0.3 0.15

guesses near the solution to ensure rapid convergence. These results are shown in Figs. 8 and 9, and again show proper agreement, although more accurate initial estimates were required. The sensitivity of solution approach to the initial guesses was further examined when the estimates used to achieve convergence were perturbed (<3%). The iterative method then converged to an incorrect, nonphysical solution, which is shown in Fig. 10 for the gas phase. This was due to the inability of the adaptive technique to determine the correct wave configuration during the Newton iteration. When the single-phase theory was used to obtain the initial guesses according to the first step of the solution procedure, this caused the wave ordering detection method to enter an infinite loop, which oscillated

160

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Fig. 8. Gas-phase solution profile for problem 2 for (a) density, (b) pressure, and (c) velocity at t = 0.6 ms.

between configurations where r = 3 and r = 2 and resulted in nonconvergence. However, prescribing the perturbed initial estimates caused the wave configuration to initially predict the nonphysical state where r = 3, and the adaptive detection technique was unable to produce stable corrective iterations, resulting in convergence to the incorrect solution. These difficulties were likely exacerbated by the close proximity of the second “solid” contact and gas shock waves shown in Fig. 10. This nearly results in a linear degeneracy and degrades stability within the wave ordering detection algorithm. It should be noted that the existence of non-unique Riemann problem solutions for two-phase flows is well-documented [1,8,25], and methods that can be used to obtain physical solutions within the Newton iteration remains a significant challenge. The development of a robust, adaptive wave configuration technique is therefore a critical element for successful implementation of the iterative solution presented here. 3.3. Problem 3 In this problem a granular solid–gas mixture representing a metalized explosive and air is considered. The solid species are the high explosive HMX (C4 H8 N8 O8 ) and aluminum, which can each be described by a Mie–Grüneisen EOS, given by:

ei ( v i , P i ) =

v ∗i 

Γi



P i − f i (v i ) ,

⎧ ⎨ P H ( v i ) − vΓ∗i e H ( v i ) if y iy−1 v ∗i < v i  v ∗i , i i f i (v i ) =  ⎩ ω2 1 − 1∗ if v > v ∗i , i i vi v i

(58)

(59)

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

161

Fig. 9. “Solid”-phase solution profile for problem 2 for (a) phase 1 density and (b) phase 1 pressure at t = 0.6 ms; and (c) phase 2 density and (d) phase 2 pressure at t = 0.26 ms.

where

 P H (v i ) =

ω v ∗i − y i ( v ∗i − v i )

2





v ∗i − v i ,

e H (v i ) =

1



ωi ( v ∗i − v i )

2 v ∗i − y i ( v ∗i − v i )

2 ;

(60)

here, Γi , y i , v ∗i , and ωi are material-dependent parameters, and i = 1, 2 for HMX and aluminum, respectively. These properties are provided in Table 4. The gas phase is taken to be air, and obeys the ideal EOS with γ g = 1.4. The initial Riemann problem states are provided in Table 5. Spatial profiles for the phase properties are shown in Figs. 11 and 12. This case corresponds to l = 0, b = 2, and r = 3, where only subsonic configurations are present, and the Al contact wave precedes the HMX contact due to the lower inertia of the more dense Al phase. Initial estimates were provided by the single-phase theory, and good agreement with the numerical solution is obtained. However, when considering stiff EOS’s for the solid phases, convergence of this method is achieved only for small jumps in solid volume fraction across the initial discontinuity in each phase (<0.2). A possible explanation is the stiffness of the Mie–Grüneisen EOS, where small specific volume perturbations result in large fluctuations in the phase pressure which magnifies numerical errors produced within the corrective Newton iterations. These can be mitigated by obtaining more accurate initial guesses than those provided by the single-phase solutions, and motivates the need to identify a method capable of producing improved estimates to enhance the robustness of the N-phase Riemann solver.

162

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Fig. 10. Nonphysical gas-phase solution profile for problem 2 for (a) density, (b) pressure, and (c) velocity at t = 0.7 ms.

Table 4 Mie–Grüneisen solid EOS parameter data [14].

HMX Al

Γ (–)

y (–)

( v ∗s )−1 (kg/m3 )

ω (m/s)

3.0 2.0

2.6 1.4

1900 2700

2740 5380

Table 5 Initial data for problem 3.

Left: Right:

ρ g (kg m−3 )

u g (m/s)

P g (MPa)

φ g (–)

5.0 1.23

0 0

1.0 0.101

0.55 0.75

ρ1 (kg m−3 )

u 1 (m/s)

P 1 (MPa)

φ1 (–)

Left: Right:

1980 1900

0 0

2500 0.101

0.2 0.1

ρ2 (kg m−3 )

u 2 (m/s)

P 2 (MPa)

φ2 (–)

Left: Right:

2750 2700

0 0

5000 0.101

0.25 0.15

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

163

Fig. 11. Gas-phase solution profile for problem 3 for (a) density, (b) pressure, and (c) velocity at t = 0.7 ms.

4. Conclusions Upwind finite-volume numerical methods are frequently used to solve systems of hyperbolic equations governing multiphase flows. These schemes are based on the solutions of Riemann problems constructed at the computational cell interfaces. In this work a multiphase Riemann solver has been formulated for an extension of the Baer–Nunziato equations, which governs flows in mixtures containing an arbitrary number of components. Solutions were obtained for several test problems, and were verified using the results of a finite-volume numerical scheme to determine the accuracy of the Riemann solver, showing proper agreement. These cases also demonstrated that the robustness of the solver is limited by both the quality of the initial estimates used in the Newton method that serves as the basis of the Riemann solver, as well as the quality of the adaptive technique used to determine the wave configuration. One possible strategy to improve the initial guesses is a linearization of the system of algebraic equations defined in this work, in a manner similar to that proposed in [25]. Although it is possible to obtain a closed-form linearization for more complex EOS’s in two-phase flows ( M = 1) [26], it is unclear whether a similar technique is feasible for N = M + 1 > 2. For the N-phase case it may be necessary to implement a linear solver, e.g., LU decomposition, to calculate the first-order perturbations required to improve the accuracy of the initial estimates. However, this approach may introduce significant extra computational expense. Another potential avenue is the use of a formal asymptotics-based approach with more general functions used in the expansion to obtain more accurate initial estimates compared to the linear approximation. Both strategies may also form the basis for an approximate Riemann solver, accurate in the limit of small |φ R , j − φ L , j | for j = 1, 2, . . . , M, which can significantly improve computational efficiency in a manner similar to the implementation discussed by Schwendeman et al. [25]. These techniques appear to be a logical strategy to provide more accurate initial estimates for the N-phase case, and may be explored in a future work. Despite the drawbacks of the current solver methodology, the present work represents the first

164

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

Fig. 12. Solid-phase solution profile for problem 3 for (a) HMX density, (b) HMX pressure, (c) Al density, and (d) Al pressure at t = 0.05 ms.

proposed N-phase analogue of the two-phase solver introduced in [25], and therefore may serve as a framework for exact or approximate solvers utilized in finite-volume schemes. Acknowledgements This work was supported by NSF-IGERT on Computational Fluid Dynamics at Louisiana State University, grant number DGE-0504507. References [1] N. Andrianov, G. Warnecke, The Riemann problem for the Baer–Nunziato two-phase flow model, J. Comput. Phys. 195 (2004) 434–464. [2] M.R. Baer, Numerical studies of dynamic compaction of inert and energetic granular materials, J. Appl. Mech. 55 (1988) 36–43. [3] M.R. Baer, R.J. Gross, J.W. Nunziato, An experimental and theoretical study of deflagration-to-detonation transition (DDT) in the granular explosive, CP, Combust. Flame 65 (1986) 15–30. [4] M.R. Baer, J.W. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials, Int. J. Multiph. Flow 12 (1986) 861–889. [5] M.W. Crochet, Modeling, numerical analysis, and predictions for the detonation of multi-component energetic solids, Ph.D. thesis, Louisiana State University, 2013, in preparation. [6] M.W. Crochet, K.A. Gonthier, Numerical investigation of a modified family of centered schemes applied to multiphase equations with nonconservative sources, J. Comput. Phys. (2013), http://dx.doi.org/10.1016/j.jcp.2013.08.010. [7] M.W. Crochet, K.A. Gonthier, J.E. Tohline, Application of a generalized multiphase Riemann solver to a finite-volume method with nozzling sources, in: Proc. APS SCCM, 2011, pp. 1471–1474. [8] V. Deledicque, M.V. Papalexandris, An exact Riemann solver for compressible two-phase flow models containing non-conservative products, J. Comput. Phys. 222 (2007) 217–245.

M.W. Crochet, K.A. Gonthier / Applied Numerical Mathematics 76 (2014) 145–165

165

[9] M. Dumbser, E.F. Toro, A simple extension of the Osher Riemann solver to non-conservative hyperbolic systems, J. Sci. Comput. 48 (2011) 70–88. [10] P. Embid, M. Baer, Mathematical analysis of a two-phase continuum mixture theory, Contin. Mech. Thermodyn. 4 (1992) 279–312. [11] P. Goatin, P.G. LeFloch, The Riemann problem for a class of resonant hyperbolic systems of balance laws, Ann. Inst. Henri Poincaré, Anal. Non Linéaire 21 (2004) 881–902. [12] K.A. Gonthier, C.F. Cox, Predictions for the impact of a granular solid having spatially non-uniform bulk porosity, Comput. Mech. 39 (2007) 419–437. [13] K.A. Gonthier, J.M. Powers, A high-resolution numerical method for a two-phase model of deflagration-to-detonation transition, J. Comput. Phys. 163 (2000) 376–433. [14] K.A. Gonthier, C.G. Rumchik, ME-TS1-09, Technical Report, Louisiana State University, 2009. [15] A. Kurganov, S. Noelle, G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton–Jacobi equations, SIAM J. Sci. Comput. 23 (2001) 707–740. [16] P. LeFloch, T.P. Liu, Existence theory for nonlinear hyperbolic systems in nonconservative form, Forum Math. 5 (1993) 261–280. [17] C.A. Lowe, Two-phase shock tube problems and numerical methods of solution, J. Comput. Phys. 204 (2005) 598–632. [18] T.P. McGrath II, Numerical modeling of multiphase explosions, Ph.D. thesis, University of Maryland, 2008. [19] N.Y. Moiseev, T.A. Mukhamadieva, Newton’s method as applied to the Riemann problem for media with general equations of state, Comput. Math. Math. Phys. 48 (2008) 1102–1110. [20] M.L. Muñoz-Ruiz, C. Parés, Godunov method for nonconservative hyperbolic systems, ESAIM, Math. Model. Numer. Anal. 41 (2007) 169–185. [21] C. Parés, M. Castro, On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems, ESAIM, Math. Model. Numer. Anal. 38 (2004) 821–852. [22] C. Parés, M.L. Muñoz-Ruiz, On some difficulties of the numerical approximation of nonconservative hyperbolic systems, Bol. Soc. Esp. Mat. Apl. 47 (2009) 19–48. [23] J.M. Powers, D.S. Stewart, H. Krier, Theory of two-phase detonation, Combust. Flame 80 (1990) 264–303. [24] R. Saurel, R. Abgrall, A multiphase Godunov method for compressible multifluid and multiphase flows, J. Comput. Phys. 150 (1999) 425–467. [25] D.W. Schwendeman, C.W. Wahle, A.K. Kapila, The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow, J. Comput. Phys. 212 (2006) 490–526. [26] D.W. Schwendeman, C.W. Wahle, A.K. Kapila, A study of detonation evolution and structure for a model of compressible two-phase reactive flow, Combust. Theory Model. 12 (2008) 159–204. [27] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 2nd edition, Springer-Verlag, 1999. [28] A.W. Vreman, Macroscopic theory of multicomponent flows: irreversibility and well-posed equations, Physica D 225 (2007) 94–111.