An Approximate Solver for Multimedium Riemann Problem with MieGrüneisen Equations of State
Abstract
We propose an approximate solver for multimedium Riemann problems with materials described by a family of general MieGrü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 wellposedness 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 multimedium flows. Numerical examples are presented for Riemann problems, air blast and underwater explosion applications.
keywords:
Multimedium Riemann problem, approximate Riemann solver, MieGrüneisen equation of state, multimedium flow1 Introduction
Numerical simulations of compressible multimedium 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 singlemedium compressible flows. However, when such an algorithm is employed to compute multimedium 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 multimedium 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 multimedium flow interactions. See Liu2003 (); Abgrall2001 (); Karni1994 (); Arienti2004 (); Shyue2001 (); Saurel1999 (); Price2015 () for instance.
A typical procedure of multimedium 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 fronttracking 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 multimedium 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 multimedium 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 fifthorder 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 MieGrüneisen EOS with variable Grüneisen coefficient. Quartapelle et al. Quartapelle2003 () proposed an exact Riemann solver by applying the NewtonRaphson 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 RoeGlaster solver to compute the equations combining the Euler equations involving chemical reaction with the MieGrüneisen EOS. More recently, Rallu Rallu2009 () and Farhat et al. Farhat2012 () utilized a sparse grid technique to tabulate the solutions of exact multimedium Riemann problems. Lee et al. Lee2013 () developed an exact Riemann solver for the MieGrü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 MieGrüneisen EOS by solving a nonlinear equation for the density increment involved in the numerical integration of rarefaction curves, and chose the JWL (JonesWilkinsLee) EOS as a representative case.
In this paper, we propose an approximate multimedium Riemann solver for a family of general MieGrü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 MieGrüneisen EOS are assumed to ensure the convexity of equations of state and the wellposedness 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 twomedium compressible flows with MieGrüneisen EOS, as an extension of the numerical scheme that was developed recently for twomedium compressible flows with ideal gases Guo2016 ().
The rest of this paper is arranged as follows. In Section 2, a solution strategy for the multimedium Riemann problem with arbitrary MieGrüneisen equations of state is presented. In Section 3, the procedures of our approximate Riemann solver are outlined, and the wellposedness and convergence are analyzed. In Section 4, the application of our Riemann solver in twomedium 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 Multimedium Riemann Problem
Let us consider the following onedimensional multimedium Riemann problem of the compressible Euler equations
(1a)  
Here is time and is spatial coordinate, and , , , and are the density, velocity, pressure, total energy and specific internal energy, respectively. The system has initial values  
(1b) 
Here the equations of state for both mediums under our consideration can be classified into the family known as the MieGrüneisen EOS. The MieGrü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 MieGrüneisen EOS is given by
(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 multimedium Riemann problem (1) is given by
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 CochranChan 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 ()
keeps positive ^{1}^{1}1 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 MieGrüneisen EOS (2), the speed of sound can be expressed as
and a tedious calculus gives us that the fundamental derivative is
(3) 
We conclude that
Lemma 1.
The solution of the multimedium 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
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 RankineHugoniot jump conditions Toro2008 () yield
(4) and
(5) where
Multiplying both sides of the equality (5) by gives rise to
For convenience we introduce the Hugoniot function as follows
(6) then the relation (5) boils down to the algebraic equation . The derivative of with respect to the density is found to be
(7) As a result, the slope of the Hugoniot locus can be found by the method of implicit differentiation, namely,
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
is monotonically increasing in the interval by the condition (C1), and
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
(2). Since the compressive limit of the density satisfies the relation , we have
(8) Obviously if . On the other hand, suppose that . Rewriting (8) as
and using the inequality resulting from the conditions (C2) and (C3)
we conclude that .
(3). This is an obvious result from the expression (7).
(4). The second derivative of with respect to the density is
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


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
is constant along the rightfacing (leftfacing) rarefaction curve, we have
(9) where the density is expressed in terms of by solving the isentropic relation
(10) Similarly, the right hand side of (9) can be expressed as a function of alone. Formally we define
Collecting both cases above, we have that
Therefore, the interface pressure is the zero of the following pressure function
And the interface velocity can be determined by
Recall that the formula of the function is given by
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.
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
and
where is the fundamental derivative (3) and
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 MieGrüneisen EOS in the following theorem.
Theorem 1.
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
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
or equivalently, the constraint given by (11). This completes the proof. ∎
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 NewtonRaphson iteration as
(12) 
with a suitable initial estimate which we may choose as, for example, the acoustic approximation Godunov1976 ()
The concavity of the pressure function leads to the following global convergence of the NewtonRaphson iteration
Corollary 1.
The NewtonRaphson iteration for (12) converges if .
Unfortunately, there is no closedform expression for the pressure function or its derivative for equations of state such as polynomial EOS, JWL EOS or CochranChan EOS. Therefore, we have to implement the iteration (12) approximately. Here we adopt the inexact Newton method Dembo1982 () instead, which is formulated as
(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
where denotes the step increment and is a fixed constant, then the sequence of inexact Newton iterates defined by
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
The error of th iterate can be written as
where the residual satisfies
Therefore
and hence . It follows that converges linearly to . ∎
Recall that and can be decomposed as
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
where is a constant. Moreover, assume that the sequence is bounded. Then the sequences of pressure and velocity defined in (13) converge, namely
Proof.
Since
we have
From Lemma 4 we conclude that . Due to the continuity of we have . From and the boundness of we know that . It follows that
Finally, by the definition of and we have
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 tradeoff between the accuracy and efficiency in a practical implementation, we apply the NewtonRaphson method to solve the Hugoniot loci, and the RungeKutta method to solve the isentropic curves.
Precisely, for the given th iterator , if , then we solve the following algebraic equation
(14) 
to obtain to a prescribed tolerance with the NewtonRaphson method
By Lemma 2, we arrive at the following convergence result
Corollary 2.
The NewtonRaphson iteration for (14) must converge if .
The values of and for the shock branch are thus taken as
(15)  
(16) 
If, on the other hand, , then we solve the following initial value problem
coupled with the initial value problem of the isentropic relations
backwards until using fourthorder RungeKutta method. The values of and are then taken as
(17)  
(18) 
where
and the intermediate states are given by