1. Introduction
The dynamical evolution of compressible fluid flows is determined by the conservation laws of mass, momentum, and energy. To obtain a complete mathematical description of fluids, the conservation laws must be supplemented with the equation of state (EOS) besides constitutive relations characterizing the material properties. In this context, the fundamental Gibbs relation
(1.1) 
connects all thermodynamic variables, the temperature , the entropy , the internal energy , the pressure , the density and the specific volume . This relation also shows that just two of them are independent. We choose and as the independent variables and express the EOS in a complete form, e.g.,
(1.2) 
In practice, the EOS is often expressed in terms of another pair of independent variables
(1.3) 
such as the MieGrüneisen EOS [14]. Our familiar examples are the EOS for polytropic gases
(1.4) 
and stiffened gases
(1.5) 
where is a constant. Note that the formula of EOS can be expressed in terms of by using the Gibbs relation (1.1) [14]. In this paper we consider the following family of EOSs.
Assumption. The pressure has a linear relation with the internal energy ,
(1.6) 
where and are two given functions of density . The EOS describing the fluid is convex [14], i.e., the fundamental derivative
(1.7) 
is strictly positive.
The family of EOS (1.6) has extensive applications as shown in numerical examples in Section 5. With this EOS, the Euler system of conservation laws consisting of mass, momentum and energy is closed and can be used to describe the dynamics of compressible fluid flows, for which the Riemann problem plays a fundamental role not only in understanding the wave patterns but also in designing numerical algorithms. There are abundant literatures about how to solve the Riemann problem [3, 8, 14, 15, 20, 21]. Many of the associated relations can be evaluated in closed form for the polytropic gases. However, several of these relations must be solved numerically for general EOSs. Kamm [11] summarizes an algorithm for solving the Riemann problem of compressible flow with general convex EOSs. The approach for the JonesWilkinsLee (JWL) EOS is coded in the ExactPack package to derive exact solutions for code verification. Since there involve an appreciable number of EOS calls in the procedure of finding the zeroes of the nonlinear equations and numerically solving ODEs for rarefaction wavces, this approach is prohibitively inefficient as the “exact” Riemann solver for a general compressible flow algorithm. As far as approximated Riemann solvers are concerned, the oversimplification may destroy the thermodynamic property that the original EOS embodies. Hence, the complexity of EOS hampers the practical application of associated Riemann solvers and the development of corresponding Godunovtype numerical schemes.
The EOS embodies the subtle thermodynamic properties of compressible fluids. The interaction of kinematics and thermodynamics determines the whole dynamical process. The more significant the thermodynamic effect is in the fluid dynamics, the more important role the EOS plays. For numerical purpose in practice, proper approximation has to be made so that numerical fluxes can be constructed effectively. However, we have to bear in mind that any approximation should inherit and preserve the inherent thermodynamic property of the EOS. To this end we propose a new strategy, a Stiffened Gas Approximation of EOS for real materials, in order to suit for the design of high resolution Godunovtype schemes. Considering EOS as a part of the dynamical system of fluids, the stiffened gas approximation means that at each local background state , the EOS (1.6) is approximated as
(1.8)  
(1.9) 
As such an approximation is applied to the Riemann problem for compressible fluid flows of real materials, a new twomaterial Riemann problem is formulated. As far as this strategy is applied for Godunovtype numerical schemes, the stiffened gas approximation is made over each control volume and the twomaterial Riemann problem is solved at each cell boundary with approximate stiffened gas EOSs. It turns out that the local Riemann solver at each cell boundary has the same simplicity as that for polytropic gases. It is noted that a similar idea can be found in [16] to derive the discretization of the nonconservative equation in a twophase flow model. However, significance of the stiffened gas approximation and its application on numerical fluxes have not been discussed yet.
In order to manifest the thermodynamic properties numerically that EOS embodies, we adopt the generalized Riemann problem (GRP) solver locally [4, 5, 6]. As we discussed in [13], the GRP solver not only plays the role of accuracy in the design of high resolution Godunovtype schemes, but also provides an approach including the thermodynamic effect in the resulting scheme. The passage is taken through precisely measuring the entropy variation (by the Gibbs relation) and its interaction with kinematic and other thermodynamic variables. Such consideration provides an evident and immediate performance when thermodynamic effect is severe such as in the LeBlanc problem (the problem for large ratios of density or pressure [19, 13]). As an immediate product of the stiffened gas approximation to the Riemann problem and the generalized Riemann problem (GRP), a new type of approximate Godunov type scheme is designed for general EOSs (1.6). The main advantages using the stiffened gas approximation are listed as follows.

Difficulty of the exact Riemann solver with general EOS is avoided, thus the resulting scheme can be more efficient and robust.

The thermodynamic relations are derived explicitly and manifested in the design of the approximate GRP solver to guarantee resulting numerical solutions obey the same or at least an approximate analogue.

The stiffened gas approximation is made just locally at each interface of control volumes. Therefore the property of EOS is maintained as much as possible (even in form).
This paper is organized by six sections. Besides the introduction here, we will show the stiffened gas approximation of the EOS and formulate an approximate twomaterial Riemann problem with stiffened gases instead of solving a singlematerial Riemann problem with general EOS in Section 2. The stiffened gas approximation of the generalized Riemann problem is presented in Section 3 and the issue of thermodynamics is discussed during the approximation. Numerical schemes are summarized in Section 4 and numerical examples are presented in Section 5 to demonstrate the performance of the resulting scheme. Concluding remarks are given in Section 6.
2. stiffened gas approximation of the Riemann problem
The Euler equations for compressible fluids are written in the form
(2.1) 
where the total energy consists of the internal energy and the kinematic energy . The system is closed by an equation of state using two independent thermodynamical variables. Here we consider the convex equation of state in the form of (1.6)
For instance, the stiffened gas takes and where and are constant parameters [9]. It can be reduced to the polytropic gas with .
The Riemann problem for (2.1) is subject to the initial data
(2.2) 
As we solve the Riemann problem (2.2) numerically with the original complex EOS (1.6), computational complexity is unavoidable, as mentioned in [11] for the exact Riemann solver:

numerical integration and interpolation for the rarefaction waves;

nonlinear equations with the intricate relations of shocks and rarefaction waves;

reasonable “starting” estimates of the starstate pressure;

many EOS calls during the procedure for solving nonlinear equations.
Due to the above reasons, the exact Riemann solver for complex EOSs is prohibitively inefficient for a general compressible flow algorithm.
In order to make the Riemann solver more efficient and robust as well as maintain details of wave interactions of the Riemann solution, we make the following stiffened gas approximation and formulate an approximate twomaterial stiffened gas Riemann problem for (2.1) and (2.2). Specifically, we have two cases.

Case 1: and . There is no need making any approximation. The Riemann solution is
(2.3) It turns out that the internal energy is defined through the EOS
(2.4) 
Case 2: or . The stiffened gas approximation is applied for the EOS locally. Then a twomaterial Riemann problem is formulated. The EOS of the fluid in the lefthand side of the material interface is approximated as
(2.5) while the EOS in the righthand side is
(2.6) Here . Note that the material interface is determined together with the solution and not known a priori, the same as the contact discontinuity for a singlematerial Riemann problem.
The outcome of the stiffened gas approximation is evident. Since the stiffened gas is a typical type of the general EOS (1.6), the Riemann solution of (2.1) and (2.2) is selfsimilar, consisting of a contact discontinuity and a combination of (expansive) rarefaction waves and (compressive) shocks. This approximate twomaterial stiffened gas Riemann problem does not give rise to any new difficulty. The corresponding Riemann solver is analogous to that for polytropic gases [5, 8, 9, 10, 21]. Let us take a RarefactionContactShock structure as an example. Then Riemann invariants and a isentropic relation provide a key for the expansion of the rarefaction wave in the left
(2.7) 
where is sound speed, given by . The RankineHugoniot (RH) relations provide the basis for resolution of the shock wave in the right,
(2.8) 
where and is the shock speed. We combine (2.7) and (2.8) to obtain
The equality of the pressure and velocity across the material interface provides the bridge to solve the Riemann problem exactly
So any nonlinear iteration scheme for this algebraic equation is the same as that for polytropic gases in the Riemann solver [21].
We close this section by concluding the gains from the stiffened gas approximation for the Riemann problem.

For the Riemann problem with general EOS, a secant method is needed since the formulations can not be written explicitly. The NewtonRaphson iterative procedure can be applied for the stiffened gas approximation so that the iteration will be more efficient.

For the exact Riemann solver of complex EOS, except the iteration for the star state and , nonlinear equations are also solved inside the rarefaction waves and shocks. However, only one nonlinear equation is needed for the approximation case. As a result, the EOS calls are substantially reduced.

Since there is the explicit relation for rarefaction waves, numerical integration and interpolation are avoided.

The Riemann problem is still “precisely” solved since all of the information about the wave configuration is known.
3. Stiffened gas approximation of the generalized Riemann problem and thermodynamic effects
In this section, we will discuss the generalized Riemann problem (GRP) for (2.1) subject to the discontinuous initial data
(3.1) 
where and are smooth functions, but generally speaking. The GRP solver serves to compute instantaneous values
(3.2) 
Denote , , and ,. Our previous studies [4, 5, 6] show that the instantaneous values (3.2) only depend on the limiting values and . Therefore we might as well discuss the “linear GRP” for (2.1) directly with the piecewise linear distribution
(3.3) 
and compute the corresponding values (3.2). Such a discussion is not limited to the theoretical significance. It is the basis of high resolution numerical schemes [4, 5, 6, 7]. More importantly, the thermodynamic effect is included in numerical fluxes through the GRP solver [13], which is a main issue that we focus on now.
As justified in [5, 7], the resolution of GRP is based on the connection with the associated Riemann problem, i.e., the Riemann problem for (2.1) subject to the initial data . This associated Riemann solution is denoted as , and the solution of GRP (2.1)(3.1) is . Then we have
(3.4) 
In Figure 3.1, we display a wave configuration of rarefactioncontactshock for GRP and its connection with the associated Riemann solution. With such a connection, we see that the local wave configuration is determined by the associated Riemann solution. Hence the local twomaterial Riemann problem, i.e., the stiffened gas approximation with the background states and , is formulated to approximate the local wave configuration. The value is computed using the twomaterial Riemann solver that we proposed in the last section.
Based on this, the GRP can be used to compute the instantaneous derivative value in (3.2) and resolved for the approximate stiffened gas generalized Riemann problem. We adopt the direct Eulerian version in [6]. With reference to Figure 3.1, the leftgoing rarefaction wave and the rightgoing shock should be resolved, respectively. Details are put in Appendix. As we discussed in [13], the GRP solver is not only used to increase accuracy of order for high resolution Godunovtype (e.g. GRP) schemes, but also reflect the thermodynamic effect. Let’s comment on this issue from the entropy variation in terms of the EOS.
In the case of EOS (1.6), we have
(3.5) 
The Gibbs relation implies
(3.6) 
Specified to the EOS (1.6), we have
(3.7) 
Hereafter the prime refers to the corresponding derivative with respect to its independent variable, and the explanation is omitted when no confusion is caused. We continue to use the Gibbs relation (1.1) and obtain all differential relations, e.g.,
(3.8) 
Hence the entropy variation can be expressed in terms of the variations of primitive variables . In particular, the initial entropy variations are measured as
(3.9) 
and by using (1.6). Hence at the GRP level, there is no need making any approximation about the entropy variation.
As the leftgoing rarefaction wave is considered, the initial entropy variation gives rise to the future variation. As in [13] (the notations are adopted there and also in Appendix), we find that there holds
(3.10) 
where is the local sound speed for approximate stiffened gases,
(3.11) 
The ratio (3.10) shows the entropy evolution across the (curved in GRP) nonisentropic rarefaction wave and reflects the thermodynamic effect, which further affects the variation of other (kinematic and thermodynamic) variables.
To see this point, we write the Euler equations in terms of socalled “Riemann invariants” ,
(3.12) 
where the Riemann invariants are
with differential relations
From (3.12), it is observed that unless the flow is isentropic or the initial distribution of data is uniform, the entropy variation () plays an important role in the dynamics. We still exemplify the wave configuration in Figure 3.1 and obtain
(3.13) 
where has an explicit expression
From this expression, we see immediately how the initial entropy variation (the variation of initial distribution) takes effect on the variation of velocity and the pressure quantitatively, i.e., the influence of thermodynamics on kinematics.
For the rightmoving shock wave associated with , the RankineHugoniot relations (2.8) hold and consist of two parts: The Hugoniot condition
(3.14) 
describes the relation of thermodynamic variables across the shock; and the kinematicthermodynamic jump relation
(3.15) 
describes how the thermodynamic change affects the kinematic jump, and vice versa. Here notations refer to (2.8) and Appendix. We track the shock trajectory to obtain
(3.16) 
where , and are put in Appendix. We remark that there is no need approximating the EOS for the derivation of (3.16). Indeed, the EOS (1.6) provides,
(3.17) 
The energy equation in (2.1) can be written as
(3.18) 
Therefore, as the LaxWendroff approach (the exchange of temporal and spatial derivatives) is used for (3.16), the stiffened gas approximation is not necessary for the derivation of (3.16). Of course, the adoption of stiffened gas approximation can simplify the calculation.
Thanks to the fact that the continuity of and across the material interface implies the continuity of the corresponding material derivatives, we combine (3.13) and (3.16) to obtain , , and subsequently , . The value follows from the EOS. Thus the instantaneous derivative value in (3.2) is available.
As a summary, we comment on the thermodynamic effects by two points.

The Gibbs relation brings the initial thermodynamics into the entropy variation across the rarefaction waves. The “kinematicthermodynamic” variables and provide a bridge between kinematical and thermodynamic variables.

Using the RankineHugoniot relations to track the singularity, the relations between kinematics and thermodynamics are built automatically across for shock waves.
4. Godunov type schemes with stiffened gas approximation
Once the stiffened gas approximation for the general EOSs (1.6) is available, the Godunov type scheme for (2.1) is immediate. We just state the new scheme in one space dimension, and there is no difference at the point of the approximation of EOS for the extension in multiple space dimensions.
Assume that the spatial computational domain is covered completely by pairwise disjoint spatial elements , with and the cell average of within is defined at time as
A typical Godunov type scheme for (2.1) assumes a piecewise polynomial approximation at each time step over control volume ,
(4.1) 
where is the space of polynomials of degree on the cell . Let be the time step size and satisfy the usual CFL constraint. The numerical solution is updated in two steps.
The Godunov type schemes with the stiffened gas approximation focuses on the computation of flux functions . Locally, the generalized Riemann problem for (2.1) is solved at each singularity point , which is shifted to , subject to the initial data
(4.3) 
The stiffened gas approximation in the last sections can be applied.
5. Numerical examples
In this section we present several test problems for the Euler equations with general nonlinear equation of states to display the performance of the resulting Godunov and GRP schemes. They are three onedimensional Riemann problems, a twodimensional Riemann problem, and a shockbubble problem. The numerical solutions by the Godunov method with exact EOSs and the Godunov/GRP methods with approximate stiffened gas EOS are denoted in plots as “Godunovexact EOS”, “Godunovapproximate EOS” and “GRPapproximate EOS”, respectively. For one dimensional Riemann problems, the exact solution obtained from the exact Riemann solver is shown as a reference solution, denoted as “exact”. The CFL number is chosen as for all cases.
5.1. 1D examples
5.1.1. Contact problem with the CochranChan (CC) EOS
This Riemann problem adopted from [17] shows the advection of a density discontinuity in the flow of uniform velocity and pressure with the CC EOS,
The material parameters for Nitromethane are given in Table 1. The initial conditions are provided in Table 2. The computational domain is for . Figure 5.1 shows the numerical solutions and the exact solution for the contact problem using cells. The numerical solutions are computed using the Godunov method with exact EOSs and the Godunov/GRP methods with approximate stiffened gas EOS, respectively. In Figure 5.1, we can observe the contact is precisely captured by the current schemes. The resolutions are almost the same by the exact Riemann solver and the approximate stiffened gas Riemann solver. The approximate GRP method have better resolution thanks to second order accuracy and the inclusion of thermodynamic effect. We also observe nonphysical pressure and velocity profiles (abeit small deviation) across the contact discontinuity. As discussed in an earlier paper [12], these nonphysical phenomena are produced by conservative schemes and arise from the discrepancy of the physical properties and the nonlinearity of the EOS. However the GRP scheme evidently improves the result compared with the first order Godunov scheme.
CC  
EOS  Mbar  Mbar  
1.134  0.0  1.19  8192  1508  4.53  1.42 
CC  
EOS  
50  40  1.134  2e4  0.1  0.5  2e4  0.1 
5.1.2. Shyue Shock Tube
This problem, proposed by Shyue [18] as a test problem for algorithm development, uses the JWL EOS
(5.1)  
(5.2)  
(5.3) 
on either side of a massless discontinuity. The JWL parameters, listed in Table 3, are an approximation for the EOS of TNT. This problem is set over the interval , with the initial conditions in Table 4.
Results for this problem at are given in Figure 5.2 with computational cells, which contains snapshot plots of the density, velocity, pressure, and internal energy. The red line represents results obtained with the general EOS exact solution code, the green circles are obtained by the Godunov scheme with the exact Riemann solver, and the purple and blue lines are obtained by the approximate Godunov and GRP schemes using the proposed stiffened gas approximation, respectively. The numerical results computed by the Godunov scheme with approximate stiffened gas EOS show the same resolution compared with the one with exact JWL EOS. The approximate GRP scheme shows better resolution with second order accuracy.
JWL  

EOS  Mbar  Mbar  
1.84  0.0  0.25  8.545  0.205  4.6  1.35 
JWL  
EOS  
50  12  1.7  10.0  0.0  1.0  0.5  0.0 
5.1.3. Lee Shock Tube
Lee et al. [12] consider another JWL 1D Riemann problem. As in the previous problem, the same JWL EOS is used on either side of the initial discontinuity. The JWL parameters, listed in Table 5, are an approximation for the EOS of LX17. Like the previous example, this problem is also run on the interval , but with the different initial conditions, as given in Table 6.
Results for this problem are given in Figure 5.3, which contains snapshot plots of the density, velocity, pressure, and internal energy. Numerical results computed by the Godunov and GRP schemes with approximate stiffened gas EOS show reasonable resolution. Spurious oscillation in pressure and velocity appears due to the sufficiently nonlinear EOS [3, 12, 17].
JWL  

EOS  Mbar  Mbar  
1.905  0.0  0.8938  632.1  0.04472  11.3  1.13 
JWL  
EOS  
50  20  0.9525  1.0  0.0  3.810  2.0  0.0 
5.2. 2D examples
5.2.1. 2D Riemann problem
The twodimensional test case consists of four initial states divided by the horizontal and vertical diaphragms, as shown in Figure 5.4 [2]. Here, just for convenience, we use nondimensional variables denoted with . Variables are normalised by a reference pressure of , a reference density chosen to be in Table 7, a reference length of , and a reference time of . For the domain of in space, the intersection of two diaphragms is located at . The initial conditions for each region are summarised in Table. 7. The JWL EOS for the detonation product of PBX9502 is used. The numerical solutions of the Riemann problem at using the Godunov and GRP schemes with approximated stiffened gas EOS are shown in Figure 5.5 in the density contours with equidistant computational cells of (upper) and (lower). We can observe their respective performances clearly.
JWL  HE  
EOS  Mbar  Mbar  
PBX9502  1.895  0  0.5  1.362e6  7.199e4  6.2  2.2 
JWL EOS  Region I  Region II  Region III  Region IV 
1  0.5  0.125  0.5  
0  0  1  1  
0  1  1  0  
1  0.25  0.025  0.25 
5.2.2. 2D shockbubble interaction with JWL EOS
This last numerical experiment shows the interaction between a right shock and a bubble. The initial condition used for the shockbubble interaction case is summarised in Table 9. For the domain of , the shocked state is set for . The radius of the bubble is and the center of it is located at . The JWL EOS with parameters for TNT is applied for this case (Table 3).
states  shocked  unshocked  bubble 

4.545  1  2  
1.737  0  0  
0  0  0  
4.369  0.5  0.5 
The numerical solutions of the shockbubble interaction problem at and using the Godunov scheme the GRP scheme, respectively, with approximated stiffened gas EOS are shown in Figure 5.7 in the density contours equidistant computational cells of (upper) and (lower). It can be seen that the second order GRP scheme with approximate EOS shows better resolution and small structures compared with the Godunov scheme.
6. More remarks
In this paper we present Godunov type schemes for general EOS by using stiffened gas approximation and highlight the thermodynamic effects of the GRP scheme. Although this work is done just for a single material, the idea is of general significance. As for the use of the approximate GRP solver and its extension, some remarks are in order.

The original nonlinear GRP solver is useful for strong nonlinear waves, since thermodynamic effect becomes significant. The GRP solver with stiffened gas EOS approximation brings a new idea for engineering applications in order to make the scheme more robust and efficient while the thermodynamic effects are also embodied into the scheme.

Following the same analogy, the GRP solver can be extend to multimaterial fluid flows and the arbitrary LagrangianEuler (ALE) framework. For a specific EOS, the Gibbs relation has the corresponding implication which will be reflected to the corresponding numerical fluxes.
Acknowledgment
Jiequan Li is supported by NSFC (Nos. 11771054, 12072042,91852207), the SinoGerman Research Group Project (No. GZ1465) and Foundation of LCP.
Appendix A The details of the stiffened gas approximate GRP solver
In this appendix we assume a typical configuration (RarefactionContactShock) and extract the detailed formulation of the twomaterial stiffened gas GRP solver. As convention [4, 5, 6], the GRP solver consists four parts, consistent with the associated Riemann solver: (A) Resolution of rarefaction waves; (B) Tracking of shocks; (C) Bridge by contact discontinuity; and (D) Compute the time derivative of density. More details can be found in [13].
(A) Resolution of rarefaction waves.
The key idea is to resolve the singularity at . It consists of three steps:

Measurement of the expansion of the rarefaction wave. Characteristic coordinates are applied to characterize the expansion of such a rarefaction wave. Let and be the integral curves, respectively, of
Denote . Then we have
where . when .

The rate of entropy variation across curved rarefaction waves. Across the curved rarefaction wave associated with , the entropy variation in the neighborhood of the singularity point has the change rate,
Given the initial state in (4.1), the initial variation of the entropy is known thanks to the Gibbs relation (1.1),

The interaction of kinematics and thermodynamics. The entropy variation strongly affects the dynamics of kinematic variables. Apply (3.12) and return to the frame, then we have
where the total (material) derivative , and
Comments
There are no comments yet.