Thermodynamical Effects and High Resolution Methods for Compressible Fluid Flows

TL;DR

This paper refines the GRP solver to incorporate thermodynamical effects into high resolution methods for compressible fluid flows, emphasizing the importance of thermodynamics in accurately resolving strong waves and demonstrating solver versatility.

Contribution

The paper introduces a refined GRP solver that integrates thermodynamic variations and works with generic equations of state, enhancing high resolution numerical methods for compressible flows.

Findings

Thermodynamic effects significantly influence wave resolution.

The refined GRP solver effectively handles strong waves.

The solver is versatile across different equations of state.

Abstract

One of the fundamental differences of compressible fluid flows from incompressible fluid flows is the involvement of thermodynamics. This difference should be manifested in the design of numerical methods and seems often be neglected in addition that the entropy inequality, as a conceptual derivative, is taken into account to reflect irreversible processes and verified for some first order schemes. In this paper, we refine the GRP solver to illustrate how the thermodynamical variation is integrated into the design of high resolution methods for compressible fluid flows and demonstrate numerically the importance of thermodynamic effect in the resolution of strong waves. As a by-product, we show that the GRP solver works for generic equations of state, and is independent of technical arguments.

Equations

$$Tds=de-\frac{p}{{\rho }^2}d\rho ,$$

$${\mathbf{u}}_t+\mathbf{f}(\mathbf{u}{)}_x=\mathbf{G}(x,\mathbf{u}),t>0,$$

$$\mathbf{u}=(\rho ,\rho u,\rho E{)}^{\top },\mathbf{f}(\mathbf{u})=(\rho u,\rho u^2+p,u(\rho E+p){)}^{\top },$$

$$\mathbf{u}(x,t_n)={\mathbf{P}}_j^{k,n}(x),x\in I_j.$$

$$\begin{array}[]{l}\mathbf{u}_{j}^{n+1}=\mathbf{u}_{j}^{n}-\displaystyle{\frac{\Delta t}{\Delta x}}[\mathbf{F}_{j+\frac{1}{2}}^{*}-\mathbf{F}_{j-\frac{1}{2}}^{*}]+\Delta t\mathbf{G}_{j}^{*},\\[8.53581pt] \mathbf{u}_{j}^{n}=\displaystyle{\frac{1}{\Delta x}}\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\mathbf{u}(x,t_{n})dx,\end{array}$$

$${\mathbf{F}}_{j+\frac{1}{2}}^{\ast }\approx \frac{1}{\Delta t}{\int }_{t_n}^{t_{n+1}}\mathbf{f}(\mathbf{u}(x_{j+\frac{1}{2}},t))dt.$$

$$\begin{array}[]{l}\mathbf{u}_{t}+\mathbf{f}(\mathbf{u})_{x}=\mathbf{G}(x,\mathbf{u}),\ \ \ x\in(x_{j-\frac{1}{2}},x_{j+\frac{3}{2}}),\ \ \ t_{n}<t<t_{n+1},\\[8.53581pt] \mathbf{u}(x,t_{n})=\left\{\begin{array}[]{ll}\mathbf{P}_{j}^{k,n}(x),&x\in I_{j},\\ \mathbf{P}_{j+1}^{k,n}(x),&x\in I_{j+1}.\end{array}\right.\end{array}$$

$$\begin{array}[]{c}\mathbf{u}_{L}=\mathbf{P}_{j}^{k,n}(0_{-}),\ \ \ \ \mathbf{u}_{R}=\mathbf{P}_{j+1}^{k,n}(0_{+}),\\[8.53581pt] \mathbf{u}_{L}^{\prime}=\displaystyle{\frac{d}{dx}}\mathbf{P}_{j}^{k,n}(0_{-}),\ \ \ \mathbf{u}_{L}^{\prime}=\displaystyle{\frac{d}{dx}}\mathbf{P}_{j+1}^{k,n}(0_{+}).\end{array}$$

$${\mathbf{u}}_{\ast }=\underset{t\to 0+}{\lim }\mathbf{u}(0,t),{\left(\frac{\partial \mathbf{u}}{\partial t}\right)}_{\ast }=\underset{t\to 0+}{\lim }\frac{\partial \mathbf{u}}{\partial t}(0,t).$$

$$\mathbf{u}(x_{j+\frac{1}{2}},t_n+\frac{\Delta t}{2})={\mathbf{u}}_{\ast }+\frac{\Delta t}{2}{\left(\frac{\partial \mathbf{u}}{\partial t}\right)}_{\ast }+\mathcal{O}(\Delta t^2),$$

$${\mathbf{F}}_{j+\frac{1}{2}}^{\ast }:=\mathbf{f}(\mathbf{u}(x_{j+\frac{1}{2}},t_n+\frac{\Delta t}{2})),$$

$${\mathbf{F}}_{j+\frac{1}{2}}^{\ast }:=\mathbf{f}({\mathbf{u}}_{\ast })+\frac{\Delta t}{2}\frac{\partial \mathbf{f}({\mathbf{u}}_{\ast })}{\partial \mathbf{u}}{\left(\frac{\partial \mathbf{u}}{\partial t}\right)}_{\ast },$$

$$\begin{array}[]{l}\mathbf{u}^{A}_{t}+\mathbf{f}(\mathbf{u}^{A})_{x}=0,\ \ \ x\in\mathbb{R},\ \ \ t>0,\\[8.53581pt] \mathbf{u}^{A}(x,0)=\left\{\begin{array}[]{ll}\mathbf{u}_{L},&x<0,\\[8.53581pt] \mathbf{u}_{R},&x>0.\end{array}\right.\end{array}$$

$$\underset{t\to 0+}{\lim }\mathbf{u}(\alpha t,t)=R^A(\alpha ;{\mathbf{u}}_L,{\mathbf{u}}_R).$$

$${\mathbf{u}}_{\ast }=R^A(0;{\mathbf{u}}_L,{\mathbf{u}}_R).$$

$$\frac{dx}{dt}=u+c,\frac{dx}{dt}=u-c.$$

$$\frac{\partial \alpha }{\partial x}dx+\frac{\partial \alpha }{\partial t}dt=0,$$

$$\frac{\partial \beta }{\partial x}dx+\frac{\partial \beta }{\partial t}dt=0.$$

$$\frac{\partial x}{\partial \alpha }=(u-c)\frac{\partial t}{\partial \alpha },\frac{\partial x}{\partial \beta }=(u+c)\frac{\partial t}{\partial \beta }.$$

$$\frac{\Theta (\beta )}{\Theta ({\beta }_L)}=\exp \left[{\int }_{{\beta }_L}^{\beta }\frac{1}{2c(0,\xi )}d\xi \right]=:\Pi (c;\beta ,{\beta }_L).$$

$$\frac{\Theta (\beta )}{\Theta ({\beta }_L)}={\left(\frac{c(0,\beta )}{c_L}\right)}^{\frac{1}{2{\mu }^2}},{\mu }^2=\frac{\gamma -1}{\gamma +1}.$$

$$\begin{array}[]{l}\displaystyle{\frac{\partial^{2}x}{\partial\alpha\partial\beta}}=\displaystyle{\frac{\partial(u-c)}{\partial\beta}}\displaystyle{\frac{\partial t}{\partial\alpha}}+(u-c)\displaystyle{\frac{\partial^{2}t}{\partial\alpha\partial\beta}},\\[8.53581pt] \displaystyle{\frac{\partial^{2}x}{\partial\alpha\partial\beta}}=\displaystyle{\frac{\partial(u+c)}{\partial\alpha}}\displaystyle{\frac{\partial t}{\partial\beta}}+(u+c)\displaystyle{\frac{\partial^{2}t}{\partial\alpha\partial\beta}}.\end{array}$$

$$2c\frac{{\partial }^{2}t}{\partial \alpha \partial \beta }=\frac{\partial (u-c)}{\partial \beta }\frac{\partial t}{\partial \alpha }-\frac{\partial (u+c)}{\partial \alpha }\frac{\partial t}{\partial \beta }.$$

$$\frac{\partial }{\partial \beta }\Theta (\beta )=\frac{1}{2c(0,\beta )}\Theta (\beta ).$$

$$\frac{\partial s}{\partial t}+u\frac{\partial s}{\partial x}=0.$$

$$T_L{s}_L^{\prime }=e_L^{\prime }-\frac{p_L}{{\rho }_L^2}{\rho }_L^{\prime },$$

$$\frac{s_x(0,\beta )}{s_L^{\prime }}=\frac{c_L}{c(0,\beta )}\exp \left[{\int }_{{\beta }_L}^{\beta }-\frac{1}{c(0,\xi )}d\xi \right]=\frac{c_L}{c(0,\beta )}\cdot {\Pi }^{-2}(c;\beta ,{\beta }_L).$$

$$\frac{T{s}_x(0,\beta )}{T_L{s}_L^{\prime }}={\left(\frac{c(0,\beta )}{c_L}\right)}^{\frac{1}{{\mu }^2}+1}.$$

$$\frac{\partial s}{\partial t}+(u+c)\frac{\partial s}{\partial x}=c\frac{\partial s}{\partial x},$$

$$\frac{\partial s}{\partial t}+(u-c)\frac{\partial s}{\partial x}=-c\frac{\partial s}{\partial x}.$$

$$\frac{\partial s}{\partial \beta }=\frac{\partial t}{\partial \beta }\frac{\partial s}{\partial t}+\frac{\partial x}{\partial \beta }\frac{\partial s}{\partial x}=\frac{\partial t}{\partial \beta }\left[\frac{\partial s}{\partial t}+(u+c)\frac{\partial s}{\partial x}\right]=\frac{\partial t}{\partial \beta }\cdot \left(c\frac{\partial s}{\partial x}\right),$$

Full text

Key words: GRP solver, Thermodynamical effects, Entropy, Riemann invariants, Kinematic-thermodynamic variables, Nonlinear geometrical optics

1. Introduction

In the study of compressible fluid flows, the thermodynamical (Gibbs) relation

[TABLE]

has always the fundamental importance, where $T$ is the temperature, $s$ is the entropy, $e$ is the internal energy, $p$ is the pressure and $\rho$ is the density. In general, the internal energy includes the static energy, chemical energy in the field of combustion, and stress tensor in elasticity etc. This relation distinguishes compressible fluid flows from incompressible fluid ones. The more compressible the flows are, the more dominant role it needs to play. This feature should be manifested in the design of numerical methods in order to guarantee resulting numerical solutions obey the same or at least an approximate analogue.

Let us recall the (first order) Godunov scheme [8] for inviscid compressible Euler equations to roughly inspect how the thermodynamics works for numerical methods since it is the reference of any first order numerical schemes and has become the foundation of modern CFD. The Godunov scheme assumes, as most first order finite volume schemes do, that the physical state is uniform over each computational cell at each time step and the flow variation is described through the jumps of states across neighboring cell boundaries. The cellwise uniformity of the initial data implies that the resulting rarefaction waves emanating from the jumps are always isentropic. Moreover, the local self-similarity of the solution implies that the entropy is constant along each cell boundary. It turns out that no thermodynamical process is included in numerical fluxes. It does not matter if the thermodynamical effect is weak and the dynamical process is not severe, just as exhibited in many popular numerical examples in literature, e.g. [14]. However, once the thermodynamical process becomes significant, an effective numerical scheme has to include the thermodynamical variation. Otherwise, it might have some defects, as shown in [13] for the problem of large ratio of density or pressure. This is insurmountable in the framework of first order schemes or higher order accurate schemes when first order numerical flux (Riemann, approximate Riemann) solvers are adopted.

In this paper, we will demonstrate how the thermodynamical effect is integrated into high order accurate numerical schemes through the study of the generalized Riemann problem (GRP). Given piecewise polynomials of high degree as initial data, all waves from the initial jumps at cell boundaries are curved and in particular rarefaction waves are no longer isentropic. The initial entropy variation activates the interaction between the kinematical and thermodynamical quantities. As the waves are sufficiently strong, the interaction cannot be neglected and the entropy variation should be plugged into numerical fluxes in order to precisely characterize the thermodynamical process, as described in the GRP approach in [4]. We recognize by careful inspection on the resolution of the GRP that the thermodynamical quantities play a fundamental role in the following sense:

1. (i)

The expansion of waves can be characterized in terms of local sound speed alone; 2. (ii)

The entropy variation rate across rarefaction waves only depends on the local sound speed; 3. (iii)

The interaction of kinematic and thermodynamical quantities is strongly influenced by the entropy variation.

These explain why GRP-based numerical schemes work well for problems under extreme conditions (e.g. high temperature and high pressure etc.) and the thermodynamics plays an important role in the design of high resolution schemes.

As another purpose, this paper aims to refine the GRP solvers that were derived before, such as the original GRP solver in [2, 3], the Eulerian version [4, 5] and high order extension [11, 12, 17]. In the development of the GRP solver, an important progress was made in [9] where the Riemann invariants were introduced so that the GRP solver can be extended to general hyperbolic balance laws. The current contribution emphasizes that the GRP solver is independent of technical analyses although the tricky technique of ”nonlinear geometric optics” is applied to the singularity point. It turns out that the technicality has nothing to do with any specific equation of state and readers are relieved of the tedious description of local characteristic coordinates in the previous versions of GRP solver in [3, 4, 5] even though all ingredients are already cooked there. Thus the conclusion is applicable beyond gas dynamics.

In order to keep the clarity of our presentation, we confine our discussion in the finite volume framework, following van Leer’s philosophy [16]. In Section 2, we simply summarize numerical methods for compressible fluid flows and relate them to the generalized Riemann problem (GRP). In Section 3, we refine the arguments on the resolution of rarefaction waves and highlight the role of thermodynamics. The technique of “nonlinear geometric optics” is applied to measure how rarefaction waves expand and how the thermodynamics takes effect on the kinematic quantities. In Section 4, we resolve shocks via the singularity tracking. The GRP solver is summarized in Section 5. In Section 6, a numerical demonstration is given to emphasize the importance of thermodynamic effect included in numerical fluxes. In Section 7, more remarks are presented about the GRP solver. All notations we use are put in Table I in Appendix.

2. Set-up of high order schemes and the generalized Riemann problem

We write the flow equations in general form,

[TABLE]

where $\mathbf{f}(\mathbf{u})$ is the flux function of physical vector $\mathbf{u}$, $\mathbf{G}(x,\mathbf{u})$ is the source term, $x$ is the spatial variable, $t$ is the temporal variable. The prototype of (2.1) is the compressible Euler equations with external forces or geometrical effects,

[TABLE]

where the primitive variables $\rho$, $u$, $p$ are density, velocity and pressure, respectively; the total energy consists of the kinematic energy $u^{2}/2$ and the internal energy $e$, $E=\frac{u^{2}}{2}+e$, the internal energy $e$ is defined through the Gibbs relation (1.1).

A high order finite volume scheme for (2.1) assumes piecewise polynomials of degree $k$ at each time step $t=t_{n}$ over each computational cell $I_{j}=(x_{j-\frac{1}{2}},x_{j+\frac{1}{2}})$,

[TABLE]

where $x_{j}=j\Delta x$, $x_{j+\frac{1}{2}}=\frac{1}{2}(x_{j}+x_{j+1})$, and $\Delta x$ is the spatial mesh size. Let $\Delta t$ be the time step size and satisfy the usual CFL constraint. The numerical solution is updated in two steps.

1. (i)

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

[TABLE]

where $\mathbf{F}_{j+\frac{1}{2}}^{*}$ is the numerical flux at the cell boundary $x=x_{j+\frac{1}{2}}$ and $\mathbf{G}_{j}^{*}$ is the proper evaluation of the source term over the control volume $I_{j}\times(t_{n},t_{n+1})$, $t_{n+1}=t_{n}+\Delta t$. 2. (ii)

Projection of data. We project the solution $\mathbf{u}(x,t_{n+1}-)$ to the space of piecewise polynomials $\mathbf{u}(x,t_{n+1})=\mathbf{P}_{j}^{k,n+1}(x)$, $x\in I_{j}$.

This paper focuses on the approximation of flux function,

[TABLE]

Such an approximation depends on the resolution of the generalized Riemann problem (GRP) at each singularity point $(x_{j+\frac{1}{2}},t_{n})$,

[TABLE]

We shift $(x_{j+\frac{1}{2}},t_{n})=(0,0)$ and denote

[TABLE]

The task of the GRP solver aims to evaluating the value $\mathbf{u}(0,t)$ along the cell boundary $x=0$,

[TABLE]

Once these are available, we evaluate

[TABLE]

and the flux across each interface $x=x_{j+\frac{1}{2}}$ is approximated as

[TABLE]

or

[TABLE]

within second order accuracy. We remark that we are satisfied with the second order accurate GRP solver because it is sufficient to be a building block so as to construct multi-stage higher order schemes adopting the strategy in [10].

Certainly, there are many approaches approximating the numerical flux through the technique of interpolation. The numerical example in Section 6 shows that any interpolation should reflect the thermodynamical effect precisely once the effect becomes prominent. This is not a trivial task. Instead of plausible interpolations, we want to see the exact information along the cell boundary so that the numerical flux can be properly constructed.

The resolution of the GRP is closely related to the associated Riemann problem with the asymptotic property, as sketched in the following proposition.

Proposition 2.1**.**

Let $\mathbf{u}(x,t)$ be the solution of (2.6), and $\mathbf{u}^{A}(x,t)=R^{A}(x/t;\mathbf{u}_{L},\mathbf{u}_{R})$ be the solution to the associated Riemann problem,

[TABLE]

Then at the singularity point $(0,0)$, there holds

[TABLE]

for any $\alpha\in\mathbb{R}$.

With this proposition, we have

[TABLE]

Although the Riemann solution $\mathbf{u}^{A}(x,t)$ is self-similar, all waves in the solution of GRP are curved. In particular, the curved rarefaction wave is no longer isentropic, and the initial variation of entropy should be carefully quantified, which is the key to plug the thermodynamical effect into the numerical flux.

Without loss of generality, we assume that a rarefaction wave emanating from $(0,0)$ moves to the left, a shock moves to the right and the $t$-axis is located inside the intermediate region, as shown in Figure 2.1 111If all waves move to one side of $x=0$, the GRP is solved upwind.. The rarefaction wave is associated with the characteristic field $u-c$, and the transversal characteristic field is $u+c$, where $u$ is the flow velocity and $c$ is the local sound speed, $c^{2}=\frac{\partial p(\rho,s)}{\partial\rho}$.

3. Resolution of rarefaction waves

In this section we will refine the arguments on the resolution of rarefaction waves in [4] although all ingredients are already cooked there. The point here is to emphasize that all arguments are independent of the specific equation of state and highlight the thermodynamic effect.

3.1. The thermodynamical effect on the expansion of rarefaction waves

In order to understand how the rarefaction wave expands. we resolve the singularity at the origin using a technique of “nonlinear geometrical optics” and come to conclusion that the local sound speed determines the expansion. The involvement of characteristic coordinates is crucial, but the conclusion is independent of technical arguments.

Let $\alpha(x,t)=C_{1}$ and $\beta(x,t)=C_{2}$ be the integral curves, respectively, of

[TABLE]

Away from the vacuum and the singularity point $(0,0)$, there is a one-to-one correspondence $(x,t)\rightarrow(\alpha,\beta)$ such that,

[TABLE]

Regarding the inverse $(\alpha,\beta)\rightarrow(x,t)$, we have

[TABLE]

In the following discussion, we often use both pairs of independent variables $(x,t)$ or $(\alpha,\beta)$. For example, as we use $(\alpha,\beta)$ as independent variables, we allude to the fact that $(x,t)=(x(\alpha,\beta),t(\alpha,\beta))$, if no confusion is caused.

Thanks to the asymptotics of the GRP to the associated Riemann problem at the singularity point, we denote by $\beta_{L}$ the slope of the wave head, by $\beta$ the speed of the wave speed inside the rarefaction wave, and by $\beta_{*}$ the speed of wave tail. Then we have the following fact.

Proposition 3.1**.**

Consider the curved rarefaction wave associated with $u-c$ and denote by $\Theta(\beta):=\frac{\partial t}{\partial\alpha}(0,\beta)$. Then we have

[TABLE]

In particular, for polytropic gases, we have

[TABLE]

Proof.

From (3.2), we obtain

[TABLE]

The subtraction of one from the other yields

[TABLE]

We can set $t(0,\beta)=0$ and $\frac{\partial(u-c)}{\partial\beta}(0,\beta)=1$. By noting that $\frac{\partial t}{\partial\beta}(0,\beta)\equiv 0$ and $\frac{\partial(u+c)}{\partial\alpha}$ is bounded, we obtain

[TABLE]

Integrating from $\beta_{L}$ to $\beta$ yields (3.3). For the polytropic gases, we have (3.4). ∎

Remark 3.2*.*

This proposition characterizes how the rarefaction wave expands near the singularity point in terms of the characteristic coordinate $\alpha$. The local sound speed $c$ solely determines the degree of expansion.

3.2. The rate of entropy variation across curved rarefaction waves

Now we want to investigate how the entropy varies across the curved rarefaction wave associated with $u-c$. The entropy function $s$ just keeps constant along any specific particle trajectory

[TABLE]

But it does not mean that the entropy is uniform in the whole regime of the curved rarefaction wave fan. Given the initial state in (2.3), the initial variation of the entropy is known thanks to the Gibbs relation (1.1),

[TABLE]

where the superscript “prime” represents the derivative of corresponding variables with respect to $x$, the subscript “L” denotes the limit from the left. Then the entropy variation across the rarefaction wave has a rate only depending on the local sound speed $c$.

Proposition 3.3**.**

Across the curved rarefaction wave associated with $u-c$, the entropy variation $s_{x}(0,\beta)$ in the neighborhood of the singularity point has the change rate,

[TABLE]

For polytropic gases, it becomes

[TABLE]

Proof.

We rewrite the equation (3.8) as,

[TABLE]

In terms of the characteristic coordinate $(\alpha,\beta)$, they become

[TABLE]

In order to derive (3.10), it suffices to measure $\displaystyle{\frac{\partial s}{\partial\alpha}}(0,\beta)$. For this purpose, we differentiate (3.12) with respect to $\alpha$,

[TABLE]

Inserting (3.6) and noting $\frac{\partial t}{\partial\beta}(0,\beta)\equiv 0$, we obtain

[TABLE]

Substituting (3.13) into this equation yields

[TABLE]

This is an ODE for $\frac{\partial s}{\partial\alpha}(0,\beta)$. We integrate it from $\beta_{L}$ to $\beta$ to get

[TABLE]

We go back to the frame of $(x,t)$-coordinate, by using (3.13), to obtain

[TABLE]

This is just (3.10). Specified to the polytropic gases, we obtain (3.11). ∎

Corollary 3.4**.**

The instantaneous change of the entropy along the interface $x=0$ is

[TABLE]

if the interface is located inside the intermediate regime.

This corollary describes the entropy change along the cell interface and reflects the thermodynamical process. The rate is independent of the specific form of the equation of state and even has nothing to do with geometric effects possibly involved. However, if certain more physical factors are included in the Gibbs relation (1.1), the variation of these factors is added into the initial entropy variation through the internal energy variation, in view of (3.9), and propagates into the intermediate region. In the next section, we will see how the entropy variation takes effect on numerical fluxes.

3.3. The interaction of kinematics and thermodynamics

In the design of numerical methods for compressible fluid flows, we often use characteristic quantities or Riemann invariants. We introduce

[TABLE]

This quantity is named as the “Riemann invariant” associated with $u-c$ when the flow is isentropic, and describes the inherent relation between the kinematic quantity $u$ and thermodynamical quantities

[TABLE]

For ideal gases, $K(\rho,s)=T/c$. However, we would like to rename it as a “ kinematic-thermodynamic” variable for non-isentropic case because it is not invariant over the whole rarefaction wave region (equivalently not invariant along the vector field defined by the eigenvector associated with $u-c$). In general, the quantity $\psi$ satisfies the equation of form,

[TABLE]

where $G(x,t)$ results from external forces or geometrical effects. Here we consider the case that $G(x,t)\equiv 0$ in order to see how the pure entropy variation takes effect. The analysis is done in the same matter as that for the entropy variation.

Since the right hand side of (3.22) is already known from (3.10), we denote it by $H(x,t)$ as a given function,

[TABLE]

We write (3.22) in terms of characteristic coordinates $(\alpha,\beta)$,

[TABLE]

Differentiating with respect to $\alpha$ yields

[TABLE]

Again, we use (3.6) and the fact that $\frac{\partial t}{\partial\beta}(0,\beta)$=0 to obtain

[TABLE]

which provides by integrating from $\beta_{L}$ to $\beta$

[TABLE]

Then there remain two issues unanswered: (i) one is about the initial value $\frac{\partial\psi}{\partial\alpha}(0,\beta_{L})$; (ii) the other is the variation of $\psi$ in terms of the physical independent variables $(x,t)$.

(i) The initial value $\frac{\partial\psi}{\partial\alpha}(0,\beta_{L})$. Note that

[TABLE]

and plug (3.22) into this identity. By setting $\beta=\beta_{L}$ we have

[TABLE]

(ii) The return to the $(x,t)$-frame. The combination of (3.28) and (3.22) for $\psi$ gives

[TABLE]

We collect (3.27), (3.29) and (3.30) to obtain

[TABLE]

We continue using (3.22) to obtain at $(0,\beta)$

[TABLE]

Denote the total (material) derivative by $D_{0}/Dt=\partial/\partial t+u\partial/\partial x$. Then we have

[TABLE]

where $d_{L}$ takes ,

[TABLE]

The integral in $d_{L}$ can be either approximated numerically for very general cases (no explicit equation of state), or integrated out explicitly. For polytropic gases $d_{L}$ is (by denoting $\theta(\beta)=c(0,\beta)/c_{L}$),

[TABLE]

Proposition 3.5**.**

Consider the curved rarefaction wave associated with $u-c$. The interaction of kinematic-thermodynamic variables can be described as

[TABLE]

where $d_{L}$ is given in (3.34) for a generic equation of state or (3.35) for polytropic gases.

4. Resolution of shocks

Let $x=x(t)$ be a shock with speed $\sigma=x^{\prime}(t)$ and separate two states $\mathbf{u}(x,t)$ in the wave front and $\bar{\mathbf{u}}(x,t)$ in the wave back. This shock is defined by the Rankine-Hugoniot relations,

[TABLE]

The second identity, the kinematic-thermodynamic relation, can be written as

[TABLE]

The last one is named the Hugoniot relation, defining the jump of thermodynamical quantities alone. With the condition that

[TABLE]

we express $\rho$ in terms of $p$,

[TABLE]

We substitute this into (4.2) and obtain

[TABLE]

Therefore the Rankine-Hugoniot relations comprise of a kinematic-thermodynamic relation, a (pure thermodynamic) Hugoniot relation and an identification of the propagation speed,

[TABLE]

The signs “$\pm$” correspond to $u\pm c$”, respectively. All the details can be found in [6].

As a key part of the GRP solver, we need to track the singularity. Inherently, we make differentiation along the shock trajectory $x=x(t)$. Denote

[TABLE]

We specify to the shock associated with $u+c$ and take the plus sign in (4.6). Then we have

[TABLE]

Note that in smooth regions there hold

[TABLE]

The Lax-Wendroff methodology is adopted to replace the spatial derivatives of solutions in the wave front by the corresponding temporal derivatives, $\frac{\partial\mathbf{u}}{\partial x}\rightarrow\frac{\partial\mathbf{u}}{\partial t}$; and replace the temporal derivatives in the wave back by the corresponding spatial derivatives $\frac{\partial\bar{\mathbf{u}}}{\partial t}\rightarrow\frac{\partial\bar{\mathbf{u}}}{\partial x}$. Taking the limit along the shock trajectory $x=x(t)$ with

[TABLE]

we can express the temporal variation of solution in the intermediate region in terms of the local Riemann solution and the spatial variation of initial data from the right, as stated in the following proposition.

Proposition 4.1**.**

Consider the shock associated with the characteristic field $u+c$. The temporal variation of solution in the intermediate region is described as

[TABLE]

where the coefficients $a_{R}$, $b_{R}$ and $d_{R}$ are given in terms of the intermediate state $\mathbf{u}_{*}$ and the initial data from the right,

[TABLE]

and

[TABLE]

Once $(D_{0}u/Dt)_{*}$ and $(D_{0}p/Dt)_{*}$ are available, we can derive the temporal variation of density in the intermediate region $(\partial\rho/\partial t)_{*}$ by tracking the shock trajectory $x=x(t)$.

Proposition 4.2**.**

If the intermediate state is located between a contact discontinuity and the right-moving shock associated with $u+c$, then we have

[TABLE]

where $g_{\rho}^{R}$, $g_{p}^{R}$, $g_{u}^{R}$ and $f_{R}$ are constant, depending on the initial data (2.3) in the right hand side and the Riemann solution $R^{A}(0;\mathbf{u}_{L},\mathbf{u}_{R})$. They are expressed in the following,

[TABLE]

and $H_{i}$, $i=1,2,3$, are

[TABLE]

We emphasize that (4.6) is independent of the specific equation of state. Specified to the polytropic cases, we refer to [4] for details.

5. The GRP solver

The GRP solver has two main versions: Acoustic version and nonlinear version. The acoustic version applies to the case that the waves from each cell boundary are weak and so to the regions in which the flow is smooth. Of course, if one would like to use a simplified version, it is a choice. The second order ADER method is the acoustic version of GRP solver. If the waves from a cell boundary is very strong, the nonlinear GRP has to be used because the thermodynamic effect becomes significant. We just provide the GRP solver when the cell interface $x=0$ is located inside the intermediate region or the sonic case. If all waves from the singularity point $(0,0)$ move to one side of the cell boundary $x=0$, $\mathbf{u}_{*}$ and $(\partial\mathbf{u}/\partial t)_{*}$ are just taken upwind.

5.1. Acoustic GRP

As $\|\mathbf{u}_{L}-\mathbf{u}_{R}\|\ll 1$, the acoustic GRP solver is adopted. Then $\mathbf{u}_{*}=\mathbf{u}_{L}=\mathbf{u}_{R}$, and $(\partial\mathbf{u}/\partial t)_{*}$ is obtained by solving the linearized system

[TABLE]

Specified to the compressible Euler equations, we have the proposition.

Proposition 5.1** (Acoustic case).**

As $\mathbf{u}_{L}=\mathbf{u}_{*}=\mathbf{u}_{R}$ and $\mathbf{u}_{L}^{\prime}\neq\mathbf{u}_{R}^{\prime}$, we have the acoustic case. If $u_{*}-c_{*}<0$ and $u_{*}+c_{*}>0$, then $(\partial u/\partial t)_{*}$ and $(\partial p/\partial t)_{*}$ can be solved as

[TABLE]

Then the quantity $(\partial\rho/\partial t)_{*}$ is calculated from the equation of state $p=p(\rho,s)$,

[TABLE]

5.2. Nonlinear GRP

As $\|\mathbf{u}_{L}-\mathbf{u}_{R}\|\gg 1$, the nonlinear GRP solver has to be applied. The numerical example in the present paper shows the importance of the nonlinear GRP solver. The nonlinear GRP solver for compressible fluid flows comprises of two parts:

1. (i)

*Kinematic-thermodynamic relation. * The material derivatives of $u$ and $p$ are obtained by solving the linear algebraic system

[TABLE]

where $(a_{L},b_{L},d_{L})$ and $(a_{R},b_{R},d_{R})$ are given in terms of local Riemann solution $R^{A}(x/t;\mathbf{u}_{L},\mathbf{u}_{R})$ and the corresponding initial data. 2. (ii)

The temporal variation of pure thermodynamical quantities. Once $(D_{0}u/Dt)_{*}$ and $(D_{0}p/Dt)_{*}$ are available, we can determine other thermodynamic quantities such as the density. If $x=0$ is located between the rarefaction wave and the contact discontinuity, the Gibbs relation or the equation of state, $p=p(\rho,s)$, is used directly,

[TABLE]

Otherwise, if $x=0$ is located between the shock and the contact discontinuity, (4.14) is used to obtain $(\partial\rho/\partial t)_{*}$.

From the derivation of (5.4), the thermodynamics is included in the coefficients $d_{L}$ and $d_{R}$ to exert on $(D_{0}u/Dt)_{*}$ and $(D_{0}p/Dt)_{*}$ and then on $(\partial\rho/\partial t)_{*}$. In the end the thermodynamic effect is exhibited in the numerical fluxes to influence numerical solutions.

6. Numerical demonstration of thermodynamical effect

In order to show the importance of thermodynamic effects, we choose the example from [13], which is equivalent to the Leblanc problem. The initial data is of Riemann-type,

[TABLE]

The polytropic index $\gamma=1.4$. The solution consists of a very strong left-going rarefaction wave, a relatively weak right-moving shock, and a contact discontinuity in the middle, besides constant states. The rarefaction wave spans the density from $10^{4}$ to $80$, while the density just jumps from around $1$ to $8$ across the shock. If the rarefaction wave is displayed with 100 grid points, the jump of density values at neighboring grid points is about $100$. Hence, the difference of density (or pressure) values at neighboring grid points inside the rarefaction wave is much larger than that near the shock. Hence any possible poor resolution of the rarefaction wave may cause the numerical inaccurate capturing of the shock (including the wave speed and the strength). The following numerical demonstration is carried out with CFL number $0.5$ for first order schemes and $0.32$ for second order schemes. All figures are displayed with $66$ points.

We start the numerical demonstration with first order methods, as shown in Figures 6.1. The first order solvers we choose are the popularly-used exact Riemann solver (Godunov scheme), HLLC solver and Roe solver with entropy fix. We carry out the computations using different grid points: $200$ points, $1000$ points and $10^{4}$ grid points, respectively. It is observed that there are large disparities of numerical solutions from the exact one even with $10^{4}$ grid points. The reason, we think, is the following. Almost all first order schemes assume that the flow is uniform at every time step inside each computational cell, which results in the fact that the associated rarefaction wave is isentropic. The uniformity of entropy leads to the fact that the entropy variation is difficult to be included in numerical fluxes unless it can be made up through the jumps from cell boundaries.

We proceed to look into the situation that schemes have accuracy of second order in space and first order in time, shown in Figure 6.2. At each time step, the MUSCL type data reconstruction is adopted, but only first order flux solvers are used. It is observed that there is no significant improvement at all and even oscillations are present near the contact discontinuities. The large disparities of numerical solutions from the exact solution result from the inaccurate resolution of rarefaction waves and the presence oscillations are due to the inconsistency of spatial and temporal accuracy.

In Figure 6.3, we keep the MUSCL type data reconstruction at each time step, but increase the temporal accuracy to second order accuracy. We first use the two-stage Rung-Kutta temporal iteration with first order flux solvers as building blocks. The result is improved a little bit but not prominently. Hence we make a trial to execute a one-stage method with acoustic GRP solver (consistent with ADER [15]). The result is displayed in Figure 6.4. The result is not satisfactory until the grid points are taken to be $10^{4}$, and even the oscillations near the edge of the rarefaction wave are more severe.

Finally, we use the nonlinear GRP solver in the strong wave regions to simulate this problem with $100$, $200$ and $300$ grid points, respectively. The result is displayed in Figure 6.5. We can see that the computation with $300$ grid points can effectively cope with the resolution of strong waves.

7. More remarks

In this paper, we refine the GRP solver and highlight the thermodynamic effect on the simulation of strong waves. Although this work is done just for compressible Euler equations, the conclusion is of general significance. As for the use of GRP solver and its extension, some remarks are in order.

1. (i)

There are two versions of the GRP solver: Acoustic and nonlinear. In smooth regions of flows or weak waves, just the acoustic solver is needed. However, once nonlinear waves become strong or thermodynamic effect becomes significant, the nonlinear GRP solver has to be used. It turns out that in our simulations, the acoustic solver is used for most of the computational time, and the nonlinear GRP solver works just near strong waves. 2. (ii)

We just refine the 1-D GRP solver with a single physical effect in the present paper. As more physical factors are included, the GRP solver can be derived similarly. To be more precise, these factors take effect on the kinematical-thermodynamic variables and the equation (3.22) is replaced by

[TABLE]

where $G_{1}$, $G_{2}$, $\cdots$, represent the aforementioned physical effects. It turns out that in (3.33) $d_{L}$ is replaced by $\tilde{d}_{L}$,

[TABLE]

The integrals for $d_{iL}$ can be evaluated either precisely (if possible) or numerically. In parallel, $d_{R}$ should include other physical effects through the Lax-Wendroff methodology. 3. (iii)

For a specific equation of state, the Gibbs relation has the corresponding implication. For example, in chemical thermodynamics [7], the internal energy comprises of a static energy $e_{0}$, the chemical potential $e_{i}$

[TABLE]

Then the Gibbs relation becomes

[TABLE]

Recall that in the computation of entropy variation, any initial variation in (3.9) is reflected in the corresponding numerical fluxes, and so does chemical potential $de_{i}$. This is an example how the GRP solver includes additional effects in practical applications. 4. (iv)

The role of GRP solver for second order schemes is just like that of the exact Riemann solver for first order schemes, but it has its own significance. We can regard the schemes based on the GRP solver as discontinuous Lax-Wendroff schemes. The spatial-temporal coupling can effectively absorb the full (physical) information provided by the governing equations, which is the spirit of the Cauchy-Kowalevski methodology.

Acknowledgement

Jiequan Li is supported by by NSFC (No. 11371063, 91130021), the doctoral program from the Education Ministry of China (No. 20130003110004). Yue Wang is supported by NSFC (No. 11501040).

TABLE I: Basic notations

[TABLE]