An Approximate Solver for Multi-medium Riemann Problem with Mie-Grüneisen Equations of State

# An Approximate Solver for Multi-medium Riemann Problem with Mie-Grüneisen Equations of State

Li Chen Ruo Li Chengbao Yao School of Mathematical Sciences, Peking University, Beijing, China HEDPS & CAPT, LMAM & School of Mathematical Sciences, Peking University, Beijing, China Northwest Institute of Nuclear Technology, Xi’an, China
###### Abstract

We propose an approximate solver for multi-medium Riemann problems with materials described by a family of general Mie-Grüneisen equations of state, which are widely used in practical applications. The solver provides the interface pressure and normal velocity by an iterative method. The well-posedness and convergence of the solver is verified with mild assumptions on the equations of state. To validate the solver, it is employed in computing the numerical flux on phase interfaces of a numerical scheme on Eulerian grids that was developed recently for compressible multi-medium flows. Numerical examples are presented for Riemann problems, air blast and underwater explosion applications.

###### keywords:
Multi-medium Riemann problem, approximate Riemann solver, Mie-Grüneisen equation of state, multi-medium flow

## 1 Introduction

Numerical simulations of compressible multi-medium flow are of great interest in practical applications, such as mechanical engineering, chemical industry, and so on. Many conservative Eulerian algorithms perform very well in single-medium compressible flows. However, when such an algorithm is employed to compute multi-medium flows, numerical inaccuracies may occur at the material interfaces Abgrall1996 (); Saurel2007 (); Liu2003 (), due to the great discrepancy of densities and equations of state across the interface. The simulation of compressible multi-medium flow in an Eulerian framework requires special attention in describing the interface connecting distinct fluids, especially for the problems that involve highly nonlinear equations of state. Several techniques have been taken to treat the multi-medium flow interactions. See Liu2003 (); Abgrall2001 (); Karni1994 (); Arienti2004 (); Shyue2001 (); Saurel1999 (); Price2015 () for instance.

A typical procedure of multi-medium compressible flows in Eulerian grids mainly consists of two steps, the interface capture and the interaction between different fluids. There are mainly two different approaches in literatures, the diffuse interface method (DIM) and the sharp interface method (SIM). The former Abgrall1996 (); Abgrall2001 (); Saurel1999 (); Saurel2009 (); Petitpas2009 (); Ansari2013 () assumes the interface as a diffuse zone, and smears out the interface over a certain number of grid cells to avoid pressure oscillations. Diffuse interfaces correspond to artificial mixtures created by numerical diffusion, and the key issue is to fulfill interface conditions within this artificial mixture. The latter assumes the interface to be a sharp contact discontinuity, and different fluids are immiscible. Several approaches such as the volume of fluid (VOF) method Scardovelli1999 (); Noh1976 (), level set method Sethian2001 (); Sussman1994 (), moment of fluid (MOF) method Ahn2007 (); Dyadechko2008 (); Anbarlooei2009 () and front-tracking method Glimm1998 (); Tryggvason2001 () are used extensively to capture the interface. A key element for both diffuse and sharp interface methods, is to determine the interface states. The accurate prediction of the interface states can be used to stabilize the numerical diffusion in diffuse interface methods, and to compute the numerical flux and interface motion in sharp interface methods. One common approach is to solve a multi-medium Riemann problem which contains the fundamentally physical and mathematical properties of the governing equations. Indeed, the Riemann problem plays a key role in understanding the wave structures, since a general fluid flow may be interpreted as a nonlinear superpositions of the Riemann solutions.

The solution of a multi-medium Riemann problem depends not only on the initial states at each side of the interface, but also on the forms of equations of state. For some simple equations of state, such as ideal gas, covolume gas or stiffened gas, the solution of the Riemann problem can be achieved to any desired accuracy with an exact solver. While the Riemann problems for the above equations of state have been fully investigated in Godunov1976 (); Plohr1988 (); Gottlieb1988 (); Toro2008 () for instance, there exist some difficulties in the cases of general nonlinear equations of state due to their high nonlinearity. A variety of methods to solve the corresponding Riemann problems have then been proposed. For example, Larini et al. Larini1992 () developed an exact Riemann solver and applied their methods to a fifth-order virial equation of state (EOS), which is particularly suited to the gaseous detonation products of high explosive compounds. Shyue Shyue2001 () developed a Roe’s approximate Riemann solver for the Mie-Grüneisen EOS with variable Grüneisen coefficient. Quartapelle et al. Quartapelle2003 () proposed an exact Riemann solver by applying the Newton-Raphson iteration to the system of two nonlinear equations imposing the equality of pressure and of velocity, assuming as unknowns the two values of the specific volume at each side of the interface, and implemented it for the van der Waals gas. Arienti et al. Arienti2004 () applied a Roe-Glaster solver to compute the equations combining the Euler equations involving chemical reaction with the Mie-Grüneisen EOS. More recently, Rallu Rallu2009 () and Farhat et al. Farhat2012 () utilized a sparse grid technique to tabulate the solutions of exact multi-medium Riemann problems. Lee et al. Lee2013 () developed an exact Riemann solver for the Mie-Grüneisen EOS with constant Grüneisen coefficient, where the integral terms are evaluated using an iterative Romberg algorithm. Banks Banks2010 () and Kamm Kamm2015 () developed a Riemann solver for the convex Mie-Grüneisen EOS by solving a nonlinear equation for the density increment involved in the numerical integration of rarefaction curves, and chose the JWL (Jones-Wilkins-Lee) EOS as a representative case.

In this paper, we propose an approximate multi-medium Riemann solver for a family of general Mie-Grüneisen equations of state in an iterative manner, which provides a strategy to reproduce the physics of interaction between two mediums across the interface. Several mild conditions on the coefficients of Mie-Grüneisen EOS are assumed to ensure the convexity of equations of state and the well-posedness of our Riemann solver. The algebraic equation of the Riemann problem is solved by an inexact Newton method Dembo1982 (), where the function and its derivative are evaluated approximately depending on the wave configuration. And the convergence of our Riemann solver is analyzed. To validate the proposed solver, we employed it in the computation of two-medium compressible flows with Mie-Grüneisen EOS, as an extension of the numerical scheme that was developed recently for two-medium compressible flows with ideal gases Guo2016 ().

The rest of this paper is arranged as follows. In Section 2, a solution strategy for the multi-medium Riemann problem with arbitrary Mie-Grüneisen equations of state is presented. In Section 3, the procedures of our approximate Riemann solver are outlined, and the well-posedness and convergence are analyzed. In Section 4, the application of our Riemann solver in two-medium compressible flow calculations Guo2016 () is briefly introduced. In Section 5, several classical Riemann problems and applications for air blast and underwater explosions are carried out to validate the accuracy and robustness of our solver. Finally, a short conclusion is drawn in Section 6.

## 2 Multi-medium Riemann Problem

Let us consider the following one-dimensional multi-medium Riemann problem of the compressible Euler equations

 ∂∂τ⎡⎢⎣ρρuE⎤⎥⎦+∂∂ξ⎡⎢⎣ρuρu2+p(E+p)u⎤⎥⎦=0,E=ρe+12ρu2. (1a) Here τ is time and ξ is spatial coordinate, and ρ, u, p, E and e are the density, velocity, pressure, total energy and specific internal energy, respectively. The system has initial values [ρ,u,p]⊤(ξ,τ=0)={[ρl,ul,pl]⊤,ξ<0,⊤,ξ>0. (1b)

Here the equations of state for both mediums under our consideration can be classified into the family known as the Mie-Grüneisen EOS. The Mie-Grüneisen EOS can be used to describe a lot of real materials, for example, the gas, water and gaseous products of high explosives Arienti2004 (); Liu2003 (); Kamm2015 (), which is particularly useful in those practical applications we are studying now. The general form of the Mie-Grüneisen EOS is given by

 p(ρ,e)=Γ(ρ)ρe+h(ρ), (2)

where is the Grüneisen coefficient, and is a reference state associated with the cold contribution resulting from the interactions of atoms at rest Heuze2012 (). Thus the EOS of the multi-medium Riemann problem (1) is given by

 p(ρ,e)={Γl(ρ)ρe+hl(ρ),Γr(ρ)ρe+hr(ρ),

for the medium on the left and the right sides, respectively. For the ease of our analysis, we impose on and the following assumptions:

(C1) ;

(C2) ;

(C3) .

###### Remark 1.

An immediate consequence is that the Grüneisen coefficient must be nonnegative since by the condition (C1).

A lot of equations of state of our interests fulfill these assumptions. Particularly, we collect some equations of state in Appendix A which are used in our numerical tests as examples. These examples include ideal gas EOS, stiffened gas EOS, polynomial EOS, JWL EOS, and Cochran-Chan EOS. The coefficients , and their derivatives for these equations of state are all provided therein.

The Riemann problem for general convex equations of state have been fully analyzed, for example, in Smith1979 (); Menikoff1989 (). Here the problem is more specific, thus we can present the structure of the solution in a straightforward way. The property of EOS is essential on the wave structures in the solution of Riemann problems. It is pointed out that the wave structures are composed solely of elementary waves Menikoff1989 () if the fundamental derivative Thompson1971 ()

 G=1c⋅∂ρc∂ρ∣∣∣s=1+ρ2c2⋅∂c2∂ρ∣∣∣s,

keeps positive 111 When the positivity condition is violated, other configurations of waves may occur, such as composite waves, split waves, expansive shock waves or compressive rarefaction waves Menikoff1989 (). The above anomalous wave structures are common issues in phase transitions. For further discussions on these anomalous wave structures, we refer the readers to Weyl1949 (); Thompson1973 (); Liu1975 (); Liu1976 (); Pego1986 (); Bethe1998 (); Bates2002 (); Voss2005 (); Muller2006 (); Fossati2014 () for instance., where is the specific entropy and is the speed of sound. For the Mie-Grüneisen EOS (2), the speed of sound can be expressed as

 c(p,ρ)=√∂p∂ρ∣∣∣e+pρ2∂p∂e∣∣∣ρ= ⎷(1ρ+Γ′k(ρ)Γk(ρ))(p−hk(ρ))+pρΓk(ρ)+h′k(ρ),

and a tedious calculus gives us that the fundamental derivative is

 G=1c2(12(ρ(ρΓk(ρ))′′+(ρΓk(ρ))′(2+Γk(ρ)))e+12ρh′′k(ρ)+p2ρ(Γ2k(ρ)+2(ρΓk(ρ))′)+12(2+Γk(ρ))h′k(ρ)). (3)

We conclude that

###### Lemma 1.

The solution of the multi-medium Riemann problem (1) consists of only elementary waves if the conditions (C1), (C2) and (C3) are fulfilled.

###### Proof.

With the conditions (C1), (C2) and (C3), a direct check on the terms in (3) gives us that

 G>0,

then the result in Menikoff1989 () is applied to give us the conclusion. ∎

Precisely, the solution of (1) consists of four constant regions connected by a linearly degenerate wave and two genuinely nonlinear waves (either shock wave or rarefaction wave, depending on the initial states), as is shown schematically in Fig. 1. The linearly degenerate wave is actually the phase interface.

Following the convention on notations, we denote the pressure and the velocity by and in the star region, respectively, which have the same value crossing over the phase interface. This allows us to derive a single nonlinear algebraic equation for . Then we will solve this algebraic equation by an iterative method. At the beginning of each step of the iterative method, a provisional value of determines the wave structures from four possible configurations. The wave structures then prescribe the specific formation of the algebraic equation. Based on the residual of the algebraic equation, the value of is updated, which closes a single loop of the iterative method.

Below let us give the details of the plan above. Firstly we need to study the relations of the solution across a nonlinear wave, saying a shock wave or a rarefaction wave. For convenience, we use the subscript or standing for either the left initial state or the right initial state hereafter.

• Solution through a shock wave

If , the corresponding nonlinear wave is a shock wave, and the star region state is connected to the adjacent initial state through a Hugoniot locus. The Rankine-Hugoniot jump conditions Toro2008 () yield

 ∓(u∗−uk)=((p∗−pk)(1ρk−1ρ∗k))1/2, (4)

and

 e(p∗,ρ∗k)−e(pk,ρk)+12(p∗+pk)(1ρ∗k−1ρk)=0, (5)

where

 e(p,ρ)=p−hk(ρ)Γk(ρ)ρ.

Multiplying both sides of the equality (5) by gives rise to

 Γk(ρk)ρk(p∗−hk(ρ∗k))−Γk(ρ∗k)ρ∗k(pk−hk(ρk))−12Γk(ρ∗k)Γk(ρk)(p∗+pk)(ρ∗k−ρk)=0.

For convenience we introduce the Hugoniot function as follows

 Φk(p,ρ):=Γk(ρk)ρk(p−hk(ρ))−Γk(ρ)ρ(pk−hk(ρk))−12Γk(ρk)(p+pk)Γk(ρ)(ρ−ρk), (6)

then the relation (5) boils down to the algebraic equation . The derivative of with respect to the density is found to be

 ∂Φk∂ρ(p,ρ)=−Γk(ρk)ρk(h′k(ρ)−p+pk2Γ′k(ρ))−Γk(ρk)(Γk(ρ)+ρΓ′k(ρ))(ρkek+p+pk2). (7)

As a result, the slope of the Hugoniot locus can be found by the method of implicit differentiation, namely,

 χ(p,ρ):=∂p∂ρ∣∣∣Φk=−2∂Φk(p,ρ)/∂ρΓk(ρk)(2ρk−Γk(ρ)(ρ−ρk)).

Before we discuss the properties of the Hugoniot function , let us introduce the compressive limit of the density such that . By definition solves the algebraic equation . This quantity is uniquely defined since the function

 W(ρ):=(ρρk−1)Γk(ρ)−2,

is monotonically increasing in the interval by the condition (C1), and

 W(ρk)=−2,limρ→+∞W(ρ)≥limρ→+∞(ρρk−1)Γ∞−2=+∞.

We have the following results on the function

###### Lemma 2.

Assume that the functions and satisfy the conditions (C1), (C2) and (C3). Then defined in (6) satisfies the following properties:

• ;

• ;

• ;

• if .

###### Proof.

(1). By definition (6) we have

 Φk(p,ρk)=Γk(ρk)ρk(p−hk(ρk))−Γk(ρk)ρk(pk−hk(ρk))=Γk(ρk)ρk(p−pk)>0.

(2). Since the compressive limit of the density satisfies the relation , we have

 Φk(p,ρmax)=Γk(ρk)ρk(p−hk(ρmax))−Γk(ρmax)ρmax(pk−hk(ρk))−12Γk(ρmax)(ρmax−ρk)Γk(ρk)(p+pk)=−Γk(ρk)ρk(pk+hk(ρmax))−(Γk(ρmax)+2)ρk(pk−hk(ρk)). (8)

Obviously if . On the other hand, suppose that . Rewriting (8) as

 Φk(p,ρmax)=−(Γk(ρk)+Γk(ρmax)+2)ρkpk−ρk(Γk(ρk)hk(ρmax)−(Γk(ρmax)+2)hk(ρk)),

and using the inequality resulting from the conditions (C2) and (C3)

 Γk(ρk)hk(ρmax)≥Γk(ρk)hk(ρk)≥(Γ∞+2)hk(ρk)≥(Γk(ρmax)+2)hk(ρk),

we conclude that .

(3). This is an obvious result from the expression (7).

(4). The second derivative of with respect to the density is

 ∂2Φk∂ρ2(p,ρ)=−Γk(ρk)ρk(h′′k(ρ)−p+pk2Γ′′k(ρ))−Γk(ρk)(2Γ′k(ρ)+ρΓ′′k(ρ))(ρkek+p+pk2),

which is negative if . This completes the proof of the whole theorem. ∎

Instantly by Lemma 2, the density can be uniquely determined from the equation on the interval for any fixed . Also the Hugoniot curve is monotonic due to . Since the equation (5) uniquely defines the interface density for a given value of , the right hand side of (4) can be regarded as a function of the interface pressure alone, formally written as

 fk(p∗)=((p∗−pk)(1ρk−1ρ∗k))1/2,p∗>pk.
• Solution through a rarefaction wave

If, on the other hand, , then the corresponding nonlinear wave is a rarefaction wave, and the interface state is connected to the adjacent initial state through a rarefaction curve. Since the Riemann invariant

 u±∫1ρcdp,

is constant along the right-facing (left-facing) rarefaction curve, we have

 ∓(u∗−uk)=∫p∗pk1ρcdp, (9)

where the density is expressed in terms of by solving the isentropic relation

 ∂p∂ρ=c2(p,ρ). (10)

Similarly, the right hand side of (9) can be expressed as a function of alone. Formally we define

 fk(p∗)=∫p∗pk1ρcdp,p∗≤pk.

Collecting both cases above, we have that

 u∗−ul=−fl(p∗) and u∗−ur=fr(p∗).

Therefore, the interface pressure is the zero of the following pressure function

 f(p):=fl(p)+fr(p)+ur−ul.

And the interface velocity can be determined by

 u∗=12(ul+ur+fr(p∗)−fl(p∗)).

Recall that the formula of the function is given by

 fk(p)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∫ppk1ρcdp,p≤pk,((p−pk)(1ρk−1ρ))1/2,p>pk,

where is determined through either the Hugoniot relation (5) or the isentropic relation (10) for a given . We claim on that

###### Lemma 3.

Assume that the conditions (C1), (C2) and (C3) hold for and , the function is monotonically increasing and concave, i.e.

 f′k(p)>0 and f′′k(p)<0,

if the Hugoniot function is concave with respect to the density, i.e. .

###### Proof.

The first and second derivatives of can be found to be

 f′k(p)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1ρc,p≤pk,12fk(p)(1ρk−1ρ+p−pkρ2χ),p>pk,

and

 f′′k(p)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−Gρ2c3,p≤pk,−14f3k(p)(2(p−pk)2ρ2χ2(1ρk−1ρ)(2ρ+∂χ∂p∣∣∣Φk)+(1ρk−1ρ−p−pkρ2χ)2),p>pk,

where is the fundamental derivative (3) and

 ∂χ∂p∣∣∣Φk=(∂Φk∂ρ(p,ρ))−1(∂2Φk∂ρ2(p,ρ)+χΓk(ρk)(ρkΓ′k(ρ)−(ρΓk(ρ))′)).

The result then follows by direct observation. ∎

The behavior of is related to the existence and uniqueness of the solution of the Riemann problem. The existence and uniqueness of the Riemann solution for gas dynamics under appropriate conditions have been established by Liu Liu1975 () and by Smith Smith1979 (). It is easy to show that the conditions (C1), (C2) and (C3) imply Smith’s “strong” condition . However, for completeness we provide a short proof of the results for the case of Mie-Grüneisen EOS in the following theorem.

###### Theorem 1.

Assume that the Mie-Grüneisen EOS (2) satisfies the conditions (C1), (C2) and (C3). The Riemann problem (1) admits a unique solution (in the class of admissible shocks, interfaces and rarefaction waves separating constant states) if and only if the initial states satisfy the constraint

 ur−ul<∫pl0dpρc+∫pr0dpρc. (11)
###### Proof.

By definition and Lemma 3 we know that the pressure function is monotonically increasing. Next we study the behavior of when tends to infinity. Let represents the density such that . When the pressure , we have , and thus

 f2k(p)=(p−pk)(1ρk−1ρ)>(p−pk)(1ρk−1~ρ).

As a result, tends to positive infinity as and so does .

Based on the behavior of the function , a necessary and sufficient condition for the interface pressure such that to be uniquely defined is given by

 f(0)=fl(0)+fr(0)+ur−ul<0,

or equivalently, the constraint given by (11). This completes the proof. ∎

###### Remark 2.

When the initial states violate the constraint (11), the Riemann problem (1) has no solution in the above sense. One can yet define a solution by introducing a vacuum. However, we are not going to address this issue since it is beyond the scope of our current interests.

## 3 Approximate Riemann Solver

The solution of the Riemann problem (1) is obtained by finding the unique zero of the function . A first try to this problem is to use the Newton-Raphson iteration as

 pn+1=pn−f(pn)f′(pn)=pn−fl(pn)+fr(pn)+ur−ulf′l(pn)+f′r(pn), (12)

with a suitable initial estimate which we may choose as, for example, the acoustic approximation Godunov1976 ()

 p0=ρlclpr+ρrcrpl+ρlclρrcr(ul−ur)ρlcl+ρrcr.

The concavity of the pressure function leads to the following global convergence of the Newton-Raphson iteration

###### Corollary 1.

The Newton-Raphson iteration for (12) converges if .

Unfortunately, there is no closed-form expression for the pressure function or its derivative for equations of state such as polynomial EOS, JWL EOS or Cochran-Chan EOS. Therefore, we have to implement the iteration (12) approximately. Here we adopt the inexact Newton method Dembo1982 () instead, which is formulated as

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩pn+1=pn−FnF′n=pn−Fn,l+Fn,r+ur−ulF′n,l+F′n,r,un=12(ul+ur+Fn,r−Fn,l), (13)

where and approximate and , respectively.

To specify the sequences and , we solve the Hugoniot loci through the numerical iteration and the isentropic curves by the numerical integration. It is natural to expect that the sequences and will tend to and respectively whenever the evaluation errors and are sufficiently small. As a preliminary result, let us introduce the following Lemma 4, which states the local convergence of inexact Newton iterates

###### Lemma 4.

If the initial iterate is sufficiently close to , and the evaluation errors of and satisfy the following constraint

 2|Fn−f(pn)|+|Δpn|⋅|F′n−f′(pn)|≤η|Fn|,

where denotes the step increment and is a fixed constant, then the sequence of inexact Newton iterates defined by

 pn+1=pn−FnF′n,

converges to . Moreover, the convergence is linear in the sense that for .

###### Proof.

Since , there exists sufficiently small that . Choose sufficiently small that

• ;

• ;

• ,

whenever . Now we prove the convergence rate by induction. Let the initial solution satisfy . Suppose that , then

 |f(pn)|=|f(pn)−f(p∗)−f′(p∗)(pn−p∗)+f′(p∗)(pn−p∗)|≤(1+γ)|f′(p∗)||pn−p∗|.

The error of -th iterate can be written as

 pn+1−p∗=f′(pn)−1(rn+(f′(pn)−f′(p∗))(pn−p∗)−(f(pn)−f(p∗)−f′(p∗)(pn−p∗))),

where the residual satisfies

 |rn|=|(f′(pn)−F′n)Δpn+f(pn)−Fn|=|(f′(pn)−F′n)Δpn+(1+η)(f(pn)−Fn)−η(f(pn)−Fn)|≤η|f(pn)|.

Therefore

 |pn+1−p∗|≤(1+γ)|f′(p∗)|−1(η|f(pn)|+γ|f′(p∗)||pn−p∗|+γ|f′(p∗)||pn−p∗|)≤(1+γ)(η(1+γ)+2γ)|pn−p∗|≤ζ|pn−p∗|,

and hence . It follows that converges linearly to . ∎

Recall that and can be decomposed as

 f(p)=fl(p)+fr(p)+ur−ul~{}~{}and~{}~{}f′(p)=f′l(p)+f′r(p).

The following Theorem 2 tells us that if and are evaluated accurately enough, the convergences of iterates and are ensured.

###### Theorem 2.

Suppose that the initial estimate is sufficiently close to . If the evaluation errors of and satisfy

 |Fn,k−fk(pn)|≤16η|Fn|~{}~{}and~{}~{}|F′n,k−f′k(pn)|≤16η|F′n|,

where is a constant. Moreover, assume that the sequence is bounded. Then the sequences of pressure and velocity defined in (13) converge, namely

 pn→p∗~{}~{}and~{}~{}un→u∗.
###### Proof.

Since

 |Fn−f(pn)|≤|Fn,l−fl(pn)|+|Fn,r−fr(pn)|≤13η|Fn|, |F′n−f′(pn)|≤|F′n,l−f′l(pn)|+|F′n,r−f′r(pn)|≤13η|F′n|,

we have

 2|Fn−f(pn)|+|Δpn|⋅|F′n−f′(pn)|≤23η|Fn|+∣∣∣FnF′n∣∣∣⋅13η|F′n|=η|Fn|.

From Lemma 4 we conclude that . Due to the continuity of we have . From and the boundness of we know that . It follows that

 |Fn,k−fk(p∗)|≤|Fn,k−fk(pn)|+|fk(pn)−fk(p∗)|→0.

Finally, by the definition of and we have

 |un−u∗|=12|(Fn,r−fr(p∗))−(Fn,l−fl(p∗))|→0,

or . ∎

By Theorem 2, the convergence is guaranteed by a posteriori control on the evaluation errors of and . The evaluation errors of and depend on the residual of the algebraic equation as well as the truncation error of the ordinary differential equation. Therefore, to achieve better accuracy one may effectively reduce the residual term and the truncation error. To achieve a trade-off between the accuracy and efficiency in a practical implementation, we apply the Newton-Raphson method to solve the Hugoniot loci, and the Runge-Kutta method to solve the isentropic curves.

Precisely, for the given -th iterator , if , then we solve the following algebraic equation

 Φk(pn,~ρn,k)=0, (14)

to obtain to a prescribed tolerance with the Newton-Raphson method

 ρn,k,m+1=ρn,k,m−Φk(pn,ρn,k,m)∂Φk(pn,ρn,k,m)/∂ρ.

By Lemma 2, we arrive at the following convergence result

###### Corollary 2.

The Newton-Raphson iteration for (14) must converge if .

The values of and for the shock branch are thus taken as

 Fn,k =((pn−pk)(1ρk−1~ρn,k))1/2, (15) F′n,k =12Fn,k(1ρk−1~ρn,k+pn−pkρ2n,kχ(pn,~ρn,k)). (16)

If, on the other hand, , then we solve the following initial value problem

 dfkdp=1ρc, fk|p=pk=0,

coupled with the initial value problem of the isentropic relations

 dρdp=1c2, ρ|p=pk=ρk,

backwards until using fourth-order Runge-Kutta method. The values of and are then taken as

 Fn,k =−16(pk−pn)(1Z1+2Z2+2Z3+1Z4), (17) F′n,k =1~ρn,kc(pn,~ρn,k), (18)

where

 ~ρn,k=ρk−16(pk−pn)(1c21+2c22+2c23+1c24),

and the intermediate states are given by

 c1 =c(pk,ρk), Z1 =ρkc1, c2 =c(pk+pn2,ρk−pk−pn2c21), Z2 =(ρk−p