Stiffened gas approximation and GRP resolution for fluid flows of real materials

The equation of state (EOS) embodies thermodynamic properties of compressible fluid materials and usually has very complicated forms in real engineering applications, subject to the physical requirements of thermodynamics. The complexity of EOS in form gives rise to the difficulty in analyzing relevant wave patterns. Concerning the design of numerical algorithms, the complex EOS causes the inefficiency of Riemann solvers and even the loss of robustness, which hampers the development of Godunov-type numerical schemes. In this paper, a strategy of local stiffened gas approximation is proposed for real materials. The stiffened gas EOS is used to approximate general EOS locally at each interface of computational control volumes so that the Riemann solver can be significantly simplified. In the meantime, the generalized Riemann problem (GRP) solver is adopted not only for high resolution purpose but effective reflection of the local thermodynamics as well. The resulting scheme is demonstrated to be efficient and robust and numerical examples display the excellent performance of such an approximation.

Comments

There are no comments yet.

Authors

• 185 publications
• 7 publications
06/11/2019

A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows

This work demonstrates a computational framework for simulating vaporizi...
08/29/2017

On the Complexity of Instationary Gas Flows

We study a simplistic model of instationary gas flows consisting of a se...
07/13/2021

One-sided GRP Solver and Numerical Boundary Conditions for compressible fluid flows

In the computation of compressible fluid flows, numerical boundary condi...
12/03/2019

Mathematical aspects relative to the fluid statics of a self-gravitating perfect-gas isothermal sphere

In the present paper we analyze and discuss some mathematical aspects of...
11/02/2020

Second-order accurate BGK schemes for the special relativistic hydrodynamics with the Synge equation of state

This paper extends the second-order accurate BGK finite volume schemes f...
09/29/2021

Pigment Adsorption Optimization in Various Low Cost Adsorbents

Abstract: The applications of adsorption are very important. The followi...
03/30/2021

State-of-the-art SPH solver DualSPHysics: from fluid dynamics to multiphysics problems

DualSPHysics is a weakly compressible smoothed particle hydrodynamics (S...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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) TdS=de+pdτ

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) p=p(ρ,S).

In practice, the EOS is often expressed in terms of another pair of independent variables

 (1.3) p=p(ρ,e)

such as the Mie-Grüneisen EOS [14]. Our familiar examples are the EOS for polytropic gases

 (1.4) p=(γ−1)ρe,   γ>1,

and stiffened gases

 (1.5) p=(γ−1)ρe−γp∞,

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) p=p(ρ,e)=κ(ρ)e+χ(ρ),

where and are two given functions of density . The EOS describing the fluid is convex [14], i.e., the fundamental derivative

 (1.7) G:=−12τ[(∂2p/∂τ2)S/(∂p/∂τ)S]

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 Jones-Wilkins-Lee (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 Godunov-type 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 Godunov-type 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) p=~p(ρ,e)=(γ(ρ0)−1)ρe−γ(ρ0)p∞(ρ0), (1.9) γ(ρ0)=1+κ(ρ0)/ρ0, p∞(ρ0)=−χ(ρ0)/γ(ρ0).

As such an approximation is applied to the Riemann problem for compressible fluid flows of real materials, a new two-material Riemann problem is formulated. As far as this strategy is applied for Godunov-type numerical schemes, the stiffened gas approximation is made over each control volume and the two-material 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 non-conservative equation in a two-phase 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 Godunov-type 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.

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

2. 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.

3. 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 two-material Riemann problem with stiffened gases instead of solving a single-material 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) ∂U∂t+∂F(U)∂x=0,U=(ρ,ρu,ρE)⊤,   F(U)=(ρu,ρu2+p,u(ρE+p))⊤,

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)

 p=p(ρ,e)=κ(ρ)e+χ(ρ).

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) U(x,0)={UL,x<0,UR,x>0.

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:

1. numerical integration and interpolation for the rarefaction waves;

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

3. reasonable “starting” estimates of the star-state pressure;

4. 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 two-material stiffened gas Riemann problem for (2.1) and (2.2). Specifically, we have two cases.

1. Case 1: and . There is no need making any approximation. The Riemann solution is

 (2.3) (ρ,u,p)(x,t)={(ρL,uL,pL),xuLt=uRt.

It turns out that the internal energy is defined through the EOS

 (2.4) pL=κ(ρL)e(x,t)+χ(ρL),xuLt=uRt.
2. Case 2: or . The stiffened gas approximation is applied for the EOS locally. Then a two-material Riemann problem is formulated. The EOS of the fluid in the left-hand side of the material interface is approximated as

 (2.5) p(ρ,e)=(γ(ρL)−1)ρe−γ(ρL)p∞(ρL);

while the EOS in the right-hand side is

 (2.6) p(ρ,e)=(γ(ρR)−1)ρe−γ(ρR)p∞(ρR).

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 single-material 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 self-similar, consisting of a contact discontinuity and a combination of (expansive) rarefaction waves and (compressive) shocks. This approximate two-material 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 Rarefaction-Contact-Shock 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) uL+2cLγ(ρL)−1=u∗+2c∗Lγ(ρL)−1,(pL+p∞(ρL))ρ−γ(ρL)L=(p∗+p∞(ρL))ρ−γ(ρL)∗L,

where is sound speed, given by . The Rankine-Hugoniot (RH) relations provide the basis for resolution of the shock wave in the right,

 (2.8) [ρu]=σ[ρ], [ρu2+p]=σ[ρu], [u(ρE+p)]=σ[ρE],

where and is the shock speed. We combine (2.7) and (2.8) to obtain

 u∗=uL−2cLγ(ρL)−1[(p∗+p∞(ρL)pL+p∞(ρL))γ(ρL)−12γ(ρL)−1]=:uL−fL(p∗,ρL,pL), u∗=uR+(p∗−pR)⎡⎢ ⎢⎣2γ(ρR)+11ρRp∗+γ(ρR)−1γ(ρR)+1pR+2γ(ρR)γ(ρR)+1p∞(ρR)⎤⎥ ⎥⎦12=:uR+fR(p∗,ρR,pR).

The equality of the pressure and velocity across the material interface provides the bridge to solve the Riemann problem exactly

 Δu+fR(p∗,ρR,pR)+fL(p∗,ρL,pL)=0, Δu=uR−uL.

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.

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

2. 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.

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

4. 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) U(x,0)={PL(x),x<0,PR(x),x>0,

where and are smooth functions, but generally speaking. The GRP solver serves to compute instantaneous values

 (3.2) U∗=limt→0+U(0,t),   (∂U∂t)∗=limt→0+∂U∂t(0,t).

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) U(x,0)={UL+U′Lx,x<0,UR+U′Rx,x>0,

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) limt→0+U(λt,t)=RA(λ;UL,UR),    λ>0.

In Figure 3.1, we display a wave configuration of rarefaction-contact-shock 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 two-material 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 two-material 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 left-going rarefaction wave and the right-going 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 Godunov-type (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) dp=∂p(ρ,e)∂ρdρ+∂p(ρ,e)∂ede.

The Gibbs relation implies

 (3.6) ∂p(ρ,e)∂eTdS=dp−(∂p(ρ,e)∂ep/ρ2+∂p(ρ,e)∂ρ)dρ.

Specified to the EOS (1.6), we have

 (3.7) dp=(κ′(ρ)e+χ′(ρ))dρ+κ(ρ)de.

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) κ(ρ)TdS=dp−(κ′(ρ)e+χ′(ρ)+κ(ρ)p/ρ2)dρ.

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) κ(ρi)TiS′i=p′i−(κ′(ρi)ei+χ′(ρi)+κ(ρi)pi/ρ2i)ρ′i,   i=L,R,

and by using (1.6). Hence at the GRP level, there is no need making any approximation about the entropy variation.

As the left-going 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) T(∂S/∂x)∗TLS′L=θ1μ2L+1∗,  θ∗=c∗LcL,  μ2L=μ2(ρL)=γ(ρL)−1γ(ρL)+1,

where is the local sound speed for approximate stiffened gases,

 (3.11) c2=γ(p+p∞)ρ.

The ratio (3.10) shows the entropy evolution across the (curved in GRP) non-isentropic 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 so-called “Riemann invariants” ,

 (3.12) ∂ϕ∂t+(u−c)∂ϕ∂x=T∂S∂x,∂ψ∂t+(u+c)∂ψ∂x=T∂S∂x,∂S∂t+u∂S∂x=0,

where the Riemann invariants are

 ϕ=u−∫ρc(ω,S)ωdω,  ψ=u+∫ρc(ω,S)ωdω,

with differential relations

 dϕ=du−1ρcdp−1cTdS,   dψ=du+1ρcdp+1cTdS.

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) (DuDt)∗+1ρ∗Lc∗L(DpDt)∗=dL,  DDt=∂∂t+u∗∂∂x,

where has an explicit expression

 dL=[1+μ2L1+2μ2Lθ1/(2μ2L)∗+μ2L1+2μ2Lθ(1+μ2L)/μ2L∗]TLS′L−cLθ1/(2μ2L)∗ψ′L.

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 right-moving shock wave associated with , the Rankine-Hugoniot relations (2.8) hold and consist of two parts: The Hugoniot condition

 (3.14) e∗R−eR+(τ∗R−τR)p∗+pR2=0,

describes the relation of thermodynamic variables across the shock; and the kinematic-thermodynamic jump relation

 (3.15) [ρu]2=[ρ][ρu2+p],

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) aR(DuDt)∗+bR(DpDt)∗=dR,

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) dp=κ(ρ)de+(κ′(ρ)e+χ′(ρ))dρ.

The energy equation in (2.1) can be written as

 (3.18) ∂e∂t+u∂e∂x+pρ∂u∂x=0.

Therefore, as the Lax-Wendroff 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.

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

2. Using the Rankine-Hugoniot 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

 Unj=1Δxj∫IjU(x,tn)dt.

A typical Godunov type scheme for (2.1) assumes a piecewise polynomial approximation at each time step over control volume ,

 (4.1) U(x,tn)=pnj(x)∈Pk(Ij),   x∈Ij,

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.

1. Average advancing. We advance the solution average of (2.1) and (4.1) to the next time step according to the formula

 (4.2) Un+1j=Unj−ΔtΔx[F∗j+12−F∗j−12],

where is the numerical flux with desired accuracy at the cell boundary .

2. Projection of data. We project the solution of (2.1)-(4.1) to the space of piecewise polynomials .

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) U(x,0)={PL(x):=pnj(x),x<0,PR(x):=pnj+1(x),x>0.

The stiffened gas approximation in the last sections can be applied.

Specified to the second order GRP scheme, the data (4.1) takes

 (4.4) pnj(x)=Unj+(Ux)nj(x−xj),   x∈Ij,

where and are constants. Reformulating this data to the form (3.3), we have

 UL=Unj+Δx2(Ux)nj,      U′L=(Ux)nj, UR=Unj+1−Δx2(Ux)nj+1, U′R=(Ux)nj+1.

Then the GRP numerical fluxes are computed as

 (4.5) F∗j+12=F(Un+12j+12), Un+12j+12:=Unj+12+Δt2(∂U∂t)nj+12,

where and are the associated Riemann solution in Section 2 and in Section 3,

 (4.6) Unj+12=RA(0;UL,UR),   (∂U∂t)nj+12=(∂U∂t)∗.

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 one-dimensional Riemann problems, a two-dimensional Riemann problem, and a shock-bubble 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 “Godunov-exact EOS”, “Godunov-approximate EOS” and “GRP-approximate 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. 1-D examples

5.1.1. Contact problem with the Cochran-Chan (C-C) EOS

This Riemann problem adopted from [17] shows the advection of a density discontinuity in the flow of uniform velocity and pressure with the C-C EOS,

 κ(ρ)=Γρ, χ(ρ)=−Γρeref(ρ)+pref(ρ), eref(ρ)=−A(1−ϵ1)ρ0[(ρ0ρ)1−ϵ1−1]+B(1−ϵ2)ρ0[(ρ0ρ)1−ϵ2−1]−e0, pref(ρ)=A(ρ0ρ)−ϵ1−B(ρ0ρ)−ϵ2.

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 non-physical pressure and velocity profiles (abeit small deviation) across the contact discontinuity. As discussed in an earlier paper [12], these non-physical 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.

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) κ(ρ)=Γρ, χ(ρ)=−Γρeref(ρ)+pref(ρ), (5.2) eref(ρ)=AR1ρ0exp(−R1ρ0ρ)+BR2ρ0exp(−R2ρ0ρ)−e0, (5.3) pref(ρ)=Aexp(−R1ρ0ρ)+Bexp(−R2ρ0ρ),

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.

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 LX-17. 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].

5.2. 2-D examples

5.2.1. 2-D Riemann problem

The two-dimensional 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 non-dimensional 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 PBX-9502 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.

5.2.2. 2-D shock-bubble interaction with JWL EOS

This last numerical experiment shows the interaction between a right shock and a bubble. The initial condition used for the shock-bubble 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).

The numerical solutions of the shock-bubble 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.

1. 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.

2. Following the same analogy, the GRP solver can be extend to multi-material fluid flows and the arbitrary Lagrangian-Euler (ALE) framework. For a specific EOS, the Gibbs relation has the corresponding implication which will be reflected to the corresponding numerical fluxes.

3. For now, we only consider the EOS with the form (1.6). It would be significant for engineering application to use the similar idea for solving fluid with more general EOSs, such as for the virial EOS and the van der Waals gas [14] or the tabular forms.

Acknowledgment

Jiequan Li is supported by NSFC (Nos. 11771054, 12072042,91852207), the Sino-German 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 (Rarefaction-Contact-Shock) and extract the detailed formulation of the two-material 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:

1. 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

 dxdt=u+c,   dxdt=u−c.

Denote . Then we have

 Θ(β)Θ(βL)=(c(0,β)cL)12μ2L,   μ2L=μ2(ρL)=γ(ρL)−1γ(ρL)+1.

where . when .

2. 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,

 T∂S/∂x(0,β)TLS′L=(c(0,β)cL)1μ2L+1.

Given the initial state in (4.1), the initial variation of the entropy is known thanks to the Gibbs relation (1.1),

 TLS′L=e′L−pLρ2Lρ′L.
3. 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

 aLDuDt(0,β)+bLDpDt(0,β)=dL(0,β),

where the total (material) derivative , and

 aL=1,   bL=1ρ(0,β)c(0,β),   ψ′L=u′L+1ρLcLp′L+1cLTLS′L, dL=[1+μ2L1+2μ2L(θ(β))1/(2μ2L)+μ2L1+2μ2L(θ(β))(1+μ2L)/μ2L]TLS′L−cL(θ(β))1/(2μ2L