A direct Eulerian GRP scheme for radiation hydrodynamical equations in diffusion limit

TL;DR

This paper introduces a second-order accurate Eulerian GRP scheme for radiation hydrodynamical equations in the diffusion limit, effectively handling nonlinear waves and discontinuities without explicit flux expressions.

Contribution

It develops a novel direct Eulerian GRP scheme that analytically resolves nonlinear waves using generalized Riemann invariants and jump conditions, improving accuracy in radiation hydrodynamics.

Findings

Achieves second-order accuracy in numerical tests

Effectively captures strong discontinuities

Demonstrates high resolution of nonlinear waves

Abstract

The paper proposes a second-order accurate direct Eulerian generalized Riemann problem (GRP) scheme for the radiation hydrodynamical equations (RHE) in the zero diffusion limit. The difficulty comes from no explicit expression of the flux in terms of the conservative vector. The characteristic fields and the relations between the left and right states across the elementary-waves are first studied, and then the solution of the one-dimensional Riemann problem is analyzed and given. Based on those, the direct Eulerian GRP scheme is derived by directly using the generalized Riemann invariants and the Runkine-Hugoniot jump conditions to analytically resolve the left and right nonlinear waves of the local GRP in the Eulerian formulation. Several numerical examples show that the GRP scheme can achieve second-order accuracy and high resolution of strong discontinuity.

Figures

Tables

Table 1. Table 4.1 : Example 4.1 : The l p superscript 𝑙 𝑝 l^{p} errors in density at t = 0.5 𝑡 0.5 t=0.5 and corresponding convergence rates of the one-dimensional GRP scheme, where p = 1 , 2 , ∞ 𝑝 1 2 p=1,2,\infty .

N	$l^1$-error	$l^1$-order	$l^2$-error	$l^2$-order	$l^{\infty }$-error	$l^{\infty }$-order
10	7.91e-4	$-$	2.85e-3	$-$	1.64e-2	$-$
20	2.23e-4	1.83	8.83e-4	1.69	6.01e-3	1.45
40	5.93e-5	1.91	2.72e-4	1.70	2.42e-3	1.31
80	1.40e-5	2.09	8.16e-5	1.73	9.57e-4	1.34
160	3.37e-6	2.05	2.43e-5	1.75	3.70e-4	1.37
320	8.14e-7	2.05	7.18e-6	1.76	1.41e-4	1.39

Table 2. Table 4.2 : Example 4.6 : The l p superscript 𝑙 𝑝 l^{p} errors in density at t = 0.5 𝑡 0.5 t=0.5 and convergence rates of the two-dimensional GRP scheme, where p = 1 , 2 , ∞ 𝑝 1 2 p=1,2,\infty .

N	$l^1$-error	$l^1$-order	$l^2$-error	$l^2$-order	$l^{\infty }$-error	$l^{\infty }$-order
10	1.09e-2	$-$	1.35e-2	$-$	2.29e-2	$-$
20	3.68e-3	1.57	4.74e-3	1.51	1.02e-2	1.16
40	1.08e-3	1.77	1.46e-3	1.70	3.93e-3	1.38
80	2.67e-4	2.01	4.38e-4	1.74	1.50e-3	1.39
160	6.35e-5	2.07	1.30e-4	1.75	5.67e-4	1.41
320	1.54e-5	2.04	3.82e-5	1.76	2.12e-4	1.42

Equations

$$\frac{\partial\mbox{\boldmath\small$U$}}{\partial t}+\frac{\partial\mbox{\boldmath\small$F(\mbox{\boldmath\small$U$})$}}{\partial x}=0,$$

$$\mbox{\boldmath\small$U$}=(\rho,m,E)^{T},\quad\mbox{\boldmath\small$F(\mbox{\boldmath\small$U$}$})=(\rho u,\rho u^{2}+p_{tot},u(E+p_{tot}))^{T}.$$

$$\mbox{\boldmath\small$E$}=\frac{1}{2}\rho^{-1}m^{2}+\rho e+\frac{\hat{a}_{R}e^{4}}{C_{v}^{4}},$$

$$\frac{\partial\mbox{\boldmath\small$W$}}{\partial t}+\mbox{\boldmath\small$\tilde{A}$}(\mbox{\boldmath\small$W$})\frac{\partial\mbox{\boldmath\small$W$}}{\partial x}=0,$$

$$\mbox{\boldmath\small$\tilde{A}$}(\mbox{\boldmath\small$W$})=\left(\begin{array}[]{ccc}u&\rho&0\\ 0&u&\frac{1}{\rho}\\ 0&\rho c^{2}&u\\ \end{array}\right),$$

$$c:=\sqrt{{\left(\frac{\partial p_{tot}}{\partial \rho }\right)}_S}=\sqrt{\frac{p}{\rho }+\frac{{\hat{K}}^2}{\rho e\tilde{K}}},\hat{K}:=p+\frac{4{\hat{a}}_R{e}^4}{3C_v^4},\tilde{K}:=\rho +\frac{4{\hat{a}}_R{e}^3}{C_v^4}.$$

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

$$S=C_v\ln \frac{p}{{\rho }^{\gamma }}+\frac{4{\hat{a}}_R{e}^3}{3\rho C_v^3}.$$

$${\lambda }_1=u-c,{\lambda }_2=u,{\lambda }_3=u+c,$$

$$\mbox{\boldmath\small$\tilde{R}$}_{1}=\left(\begin{array}[]{c}1\\ -\frac{c}{\rho}\\ c^{2}\\ \end{array}\right),\ \mbox{\boldmath\small$\tilde{R}$}_{2}=\left(\begin{array}[]{c}1\\ 0\\ 0\\ \end{array}\right),\ \mbox{\boldmath\small$\tilde{R}$}_{3}=\left(\begin{array}[]{c}1\\ \frac{c}{\rho}\\ c^{2}\\ \end{array}\right).$$

$$\nabla_{\mbox{\boldmath\small$W$}}\lambda_{2}(\mbox{\boldmath\small$W$})\cdot\mbox{\boldmath\small$\tilde{R}$}_{2}(\mbox{\boldmath\small$W$})=0,$$

$$\begin{split}\nabla_{\mbox{\boldmath\small$W$}}\lambda_{1}(\mbox{\boldmath\small$W$})\cdot\mbox{\boldmath\small$\tilde{R}$}_{1}(\mbox{\boldmath\small$W$})=-\frac{c}{\rho}-\frac{\partial c}{\partial\rho}-c^{2}\frac{\partial c}{\partial p_{tot}}=-\frac{\partial\rho c}{\rho\partial\rho}-c^{2}\frac{\partial c}{\partial p_{tot}},\end{split}$$

$$\rho c=\frac{1}{3}\rho \sqrt{e}\sqrt{\frac{(9{\gamma }_1^2+60{\gamma }_1{r}_e+9{\gamma }_1+16{\gamma }_1^2)}{4r_e+1}}.$$

$$\frac{\partial \rho c}{\partial \rho }=\frac{(540{\gamma }_1^3{r}_e+1152{\gamma }_1^2{r}_e^2+1344{\gamma }_1{r}_e^3+256r_e^4+27{\gamma }_1^3+72{\gamma }_1^2{r}_e+72{\gamma }_1{r}_e)\sqrt{e}}{6(4r_e+1{)}^{\frac{3}{2}}(3{\gamma }_1+4r_e)\sqrt{9{\gamma }_1^2+60{\gamma }_1{r}_e+16r_e^2+9{\gamma }_1}},$$

$$\frac{\partial c}{\partial p_{tot}}=\frac{256r_e^3+240{\gamma }_1{r}_e^2+112r_e^2+168{\gamma }_1{r}_e-72{\gamma }_1^2{r}_e+9{\gamma }_1^2+9{\gamma }_1}{2\rho (4r_e+1{)}^{\frac{3}{2}}(3{\gamma }_1+4r_e)\sqrt{e}\sqrt{9{\gamma }_1^2+60{\gamma }_1{r}_e+16r_e^2+9{\gamma }_1}}.$$

$$\displaystyle\nabla_{\mbox{\boldmath\small$W$}}\lambda_{1}(\mbox{\boldmath\small$W$})$$

$$f(r_e,{\gamma }_1):=$$

$$+(7488r_e^3+4464r_e^2+756r_e+54){\gamma }_1+1792r_e^4+640r_e^3.$$

$$\frac{\partial f}{\partial {\gamma }_1}=81(-8r_e+1){\gamma }_1^2+18(192r_e^2+120r_e+9){\gamma }_1+7488r_e^3+4464r_e^2+756r_e+54.$$

$$\frac{\partial f}{\partial {\gamma }_1}(r_e,0)>0,\frac{{\partial }^{2}f}{\partial {\gamma }_1^2}=162(-8r_e+1){\gamma }_1+54(64r_e^2+40r_e+3),$$

$$\mathcal{R}:=\{{\gamma }_1|f(r_e,{\gamma }_1)>0,\forall r_e\ge 0\}=\{{\gamma }_1|0<{\gamma }_1<{\hat{\gamma }}_1\}\subset [0,\infty ),$$

$${\lambda }_1=u-c:$$

$${\lambda }_2=u:$$

$${\lambda }_3=u+c:$$

$${\rho }_L{\hat{u}}_L={\rho }_R{\hat{u}}_R,$$

$${\rho }_L{\hat{u}}_L^2+p_{tot,L}={\rho }_R{\hat{u}}_R^2+p_{tot,R},$$

$${\hat{u}}_L({\hat{E}}_L+p_{tot,L})={\hat{u}}_R({\hat{E}}_R+p_{tot,R}),$$

$$\begin{array}[]{cc}p_{tot,L}<p_{tot,R},\quad u_{L}>u_{R},&\mbox{ for shock},\\ p_{tot,L}>p_{tot,R},\quad u_{L}<u_{R},&\mbox{ for rarefaction}.\end{array}$$

$$p_{tot,L}=p_{tot,R},u_L=u_R.$$

$$\begin{array}[]{cc}p_{tot,L}>p_{tot,R},\quad u_{L}>u_{R},&\mbox{ for shock},\\ p_{tot,L}<p_{tot,R},\quad u_{L}<u_{R},&\mbox{ for rarefaction}.\end{array}$$

$$u_R-c_R>u_L-c_L,$$

$$u_L-u_R+{\int }_{{\rho }_R}^{{\rho }_L}\frac{c}{\rho }d\rho =0.$$

$$c_L-c_R+{\int }_{{\rho }_R}^{{\rho }_L}\frac{c}{\rho }d\rho >0.$$

Full text

1 Introduction

The radiation hydrodynamical equations (RHE) frequently appear in astrophysics [14]. In the zero diffusion limit, they can be written as a hyperbolic system of conservation laws, but are different from the Euler equations of gas dynamics. One of major difficulties of standard numerical methods for the RHEs is to resolve and keep track of strong shock waves, while another major difficulty comes from no explicit expression of the flux in terms of the conservative vector due to the nonlinearity in the radiation pressure and energy. Dai and Woodward [6] proposed the Godunov type schemes including linear and nonlinear Riemann solvers for the RHEs. The numerical results showed that such method kept the principle advantages of the Godunov scheme, but was much time-consuming. The gas kinetic schemes are also extended to the RHEs in the zero diffusion limit. Tang and Wu first studied the kinetic flux vector splitting schemes [18], and then Jiang and Sun [9] extended it to the second-order BGK schemes.

The GRP scheme, as an analytic high-order accurate extension of the Godunov method, was originally devised for compressible fluid dynamics [1] by utilizing a piecewise linear function to approximate the “initial” data and then analytically resolving a local GRP at each interface to yield numerical flux, see the comprehensive description in [2]. There exist two versions of the original GRP scheme: the Lagrangian and Eulerian. The Eulerian version is always derived by using the Lagrangian framework with a transformation, which is quite delicate, particularly for the sonic and multi-dimensional cases. To avoid those difficulties, second-order accurate direct Eulerian GRP schemes were respectively developed for the shallow water equations [10], the Euler equations [4], and a more general weakly coupled system [3] by directly resolving the local GRPs in the Eulerian formulation via the generalized Riemann invariants and Rankine-Hugoniot jump conditions. A comparison of the GRP scheme with the gas-kinetic scheme showed the good performance of the GRP solver for some inviscid flow simulations [11]. Combined the GRP scheme with the moving mesh method [19], the adaptive direct Eulerian GRP scheme was well developed in [7] with improved resolution as well as accuracy. The accuracy and performance of the adaptive GRP scheme were further studied in [8] in simulating two-dimensional complex wave configurations formulated with the two-dimensional Riemann problems of Euler equations. The adaptive GRP scheme was also extended to unstructured triangular meshes [12]. The second order and third-order accurate extensions of the direct Eulerian GRP scheme were recently studied for one- and two-dimensional non-relativistic and relativistic Euler equations in [24, 25, 22, 23, 21] and the general one-dimensional hyperbolic balance laws in [15]. Moreover, a two-stage fourth order time-accurate GRP scheme was proposed for hyperbolic conservation laws in [13].

The aim of the paper is to develop a second-order accurate direct Eulerian GRP scheme for the RHEs. Due to the nonlinearity in the radiation pressure and energy, there is no explicit expression of the flux in terms of the conservative vector and developing the GRP schemes for the RHEs becomes not trivial and much more technical than for the Euler equations of gas dynamics. The paper is organized as follows. Section 2 introduces the one-dimensional RHEs in the zero diffusion limit, gives corresponding eigenvalues and eigenvectors, Rankine-Hugoniot jump conditions, and generalized Riemann invariants as well as their total differentials, and discusses the characteristic fields. The second-order accurate direct Eulerian GRP schemes are developed in Section 3. Section 3.1 gives the outline of the one-dimensional scheme. Section 3.2 resolves the local classical and generalized Riemann probems analytically, The two-dimensional extension of the GRP scheme is given in Section 3.3. Several numerical experiments are conducted in Section 4 to demonstrate the performance and accuracy of the proposed GRP schemes. Section 5 concludes the paper.

2 Preliminaries and notations

This section introduces the one-dimensional radiation hydrodynamical equations (RHEs) in the zero diffusion limit, gives corresponding eigenvalues and eigenvectors, Rankine-Hugoniot jump conditions, and generalized Riemann invariants as well as their total differentials, and discusses the characteristic fields.

2.1 One-dimensional RHEs

If assuming that the radiative temperature is equal to the fluid temperature, and the gas is (extremely) radiative opaque, that is, the equilibrium diffusion is dealt with and the mean-free-path of photons is much smaller than the typical length of flows, then the one-dimensional RHEs without the radiative heat diffusivity can be described as follows [14]

[TABLE]

where

[TABLE]

Here the total pressure $p_{tot}:=p+\frac{1}{3}\hat{a}_{R}T^{4}$, the momentum $m=\rho u$, $\rho,u,E,p$ and $T$ denote the density, velocity, total energy, pressure and temperature, respectively, $\hat{a}_{R}$ is the radiation coefficient and assumed to be a nonnegative constant. The total energy $E$ is described as $E=\frac{1}{2}\rho u^{2}+\rho e_{tot}$, where the total specific internal energy $e_{tot}:=e+\hat{a}_{R}\rho^{-1}T^{4}$, the specific internal energy $e$ satisfies $e=C_{v}T$, and $C_{v}$ is the specific heat at constant volume and assumed to be a positive constant. Without loss of generality, $C_{v}$ is assumed to be 1 [6, 18, 9]. The total specific enthalpy $h_{tot}$ can be written into $h_{tot}:=e_{tot}+\rho^{-1}p_{tot}$. The terms $\frac{1}{3}\hat{a}_{R}T^{4}$ and $\hat{a}_{R}T^{4}$ denote the radiation pressure and energy, respectively. To close the above system, the equation of state (EOS) is needed, $p=p(\rho,e)$. The paper only considers the EOS for the perfect gas, that is, $p=(\gamma-1)\rho e$, where $\gamma$ denotes adiabatic index. If the radiation coefficient $\hat{a}_{R}$ is zero, then the RHEs (2.1) reduces to the one-dimensional Euler equations in gas dynamics. The readers are referred to [14] for the details of the RHEs.

For the purpose of numerically solving the RHEs (2.1), the “primitive” variable vector \mbox{\boldmath\smallW}=(\rho,u,p_{tot})^{T} can be first calculated from the known conservative vector $U$ at each time step in order to calculate the values of flux \mbox{\boldmath\smallF(\mbox{\boldmath\small$U$}}) and corresponding eigenvalues etc. It is trivial for the Euler equations in gas dynamics, but becomes nontrivial for (2.1) via (numerically) solving a quartic equation in terms of the specific internal energy

[TABLE]

by any standard root-finding algorthm, e.g. Newton’s method, to get the specific internal energy $e$, where the values of $U$ are given. Fortunately, the right hand side of the equation is convex so that the Newton’s method converges with any positive initial guess. As soon as the value of $e$ is gotten, the value of the total pressure $p_{tot}$ can be calculated.

When the solution of the RHEs (2.1) is smooth, it can be cast in the quasi-linear form

[TABLE]

with

[TABLE]

where $c$ is the speed of sound given by

[TABLE]

Here $S$ denotes the entropy related to the other thermodynamical variables via the thermodynamic relation

[TABLE]

and the definition and meaning of the entropy are similar to those in the Euler equations, see [14], so that one has

[TABLE]

Three real eigenvalues and corresponding right eigenvectors of matrix \mbox{\boldmath\smallA}(\mbox{\boldmath\smallW}) can be derived as follows

[TABLE]

and

[TABLE]

2.2 Characteristic fields, and relations between the left and right states across

the elementary-waves

For the one-dimensional RHE (2.1), three characteristic fields defined by $\lambda_{i}$ [20, Section 2.4.3] satisfy the following properties, $i=1,2,3$.

Lemma 2.1

(i) The second characteristic field is linearly degenerate. (ii) the first and third characteristic fields are genuinely nonlinear if $1<\gamma<15.95$ but not genuinely nonlinear if $\gamma>15.96$.

• Proof

• (i)

From (2.5) and (2.6), it is easy to verify that the eigenvalue $\lambda_{2}$ and corresponding eigenvector \mbox{\boldmath\small\tilde{R}}_{2} satisfy

[TABLE]

where \nabla_{\mbox{\boldmath\smallW}}\lambda_{2}(\mbox{\boldmath\smallW})=(\partial\lambda_{2}/\partial\rho,\partial\lambda_{2}/\partial u,\partial\lambda_{2}/\partial p_{tot}). It tells us that the second characteristic field is linearly degenerate.

• (ii)

For the first characteristic field, one has

[TABLE]

where $c$ has been considered as a function of density and total pressure, i.e. $c=c(\rho,p_{tot})$.

For the sake of convenience, denote $\gamma_{1}:=\gamma-1$, and $r_{e}:={\hat{a}_{R}T^{4}}/{\rho e},$ which denotes the ratio of radiation energy to internal energy. Thus the product $\rho c$ can be written as follows

[TABLE]

Taking the first order partial derivative of $\rho c$ with respect to $\rho$ gives

[TABLE]

and taking the first order partial derivative of $c$ with respect to $p_{tot}$ yields

[TABLE]

Substituting them into (2.7) gives

[TABLE]

where

[TABLE]

whose first order partial derivative with respect to $\gamma_{1}$ is equal to

[TABLE]

In the following, the zero of $f(r_{e},\gamma_{1})$ is discussed case by case.

Case (i): $0\leq r_{e}\leq\frac{1}{8}$. It is obvious that $\frac{\partial f}{\partial\gamma_{1}}$ is always bigger than [math] for $\gamma_{1}\in[0,\infty)$ so that $f(r_{e},\gamma_{1})>f(r_{e},0)>0$.

Case (ii): $r_{e}>\frac{1}{8}$. Because

[TABLE]

the function $\frac{\partial f}{\partial\gamma_{1}}$ is first increasing and then decreasing monotonically as $\gamma_{1}$ increases, and has only one positive zero with respect to $\gamma_{1}\in[0,\infty)$. Thus the function $f(r_{e},\gamma_{1})$ is also first increasing and then decreasing monotonically in the interval $[0,\infty)$ as $\gamma_{1}$ increases. On the other hand, because $f(r_{e},0)>0$, the function $f(r_{e},\gamma_{1})$ has only one positive zero with respect to $\gamma_{1}$. Hence, if there exist some $r_{e}\geq 0$ and $\tilde{\gamma}_{1}$ such that $f(r_{e},\tilde{\gamma}_{1})=0$, then the inequality $f(r_{e},\gamma_{1})<0$ holds for such $r_{e}$ and all $\gamma_{1}$ bigger than $\tilde{\gamma}_{1}$. It implies that

[TABLE]

where $\hat{\gamma}_{1}:={\inf}\{\gamma_{1}|\exists r_{e}\geq 0,~{}\mbox{\rm s.t. }f(r_{e},\gamma_{1})<0\}$. So one may observe the maximum real zero of $f(r_{e},\gamma_{1})$ with respect to the independent variable $r_{e}$, and denote it by $r_{e}^{\max}$. If $r_{e}^{\max}\geq 0$, then there exist some $r_{e}=r_{e}^{\max}\geq 0$ and $\tilde{\gamma}_{1}$ such that $f(r_{e}^{\max},\tilde{\gamma}_{1})=0$ so that $\gamma_{1}\notin\mathcal{R}$. Fig. 2.1 shows the maximum real zero of $f(r_{e},\gamma_{1})$ with respect to $r_{e}$. It is seen that $\hat{\gamma}_{1}\in(14.95,14.96)$.

Combing the above two cases gives that if $1<\gamma<15.95$, then $f(r_{e},\gamma_{1})>0$ for all $r_{e}$ bigger than zero so that \nabla_{\mbox{\boldmath\smallW}}\lambda_{1}(\mbox{\boldmath\smallW})\cdot\mbox{\boldmath\small\tilde{R}}_{1}(\mbox{\boldmath\smallW})<0 for any physically admissible state $W$. Thus the first characteristic field is genuinely nonlinear when $1<\gamma<15.95$, but not genuinely nonlinear when $\gamma>15.96$.

• (iii)

The conclusion on the third characteristic field can be similarly verified. The proof is completed.

Let us study the relations between the left and right states across the elementary-waves, which need the generalized Riemann invariants and the Rankine-Hugoniot relations.

The generalized Riemann invariants provide a powerful tool of analysis of hyperbolic conservation laws. For the RHEs (2.1), the generalized Riemann invariants corresponding to three characteristic fields can be calculated as follows

[TABLE]

whose total differentials will be presented later. For a discontinuous wave solution of speed $s$, the Rankine-Hugoniot relations of the RHEs (2.1) state

[TABLE]

where $\hat{u}_{L}:=u_{L}-s$, $\hat{u}_{R}:=u_{R}-s$, $\hat{E}_{L}:=\rho_{L}e_{tot,L}+\frac{1}{2}\rho_{L}\hat{u}^{2}_{L}$, $\hat{E}_{R}:=\rho_{R}e_{tot,R}+\frac{1}{2}\rho_{R}\hat{u}^{2}_{R}$, and the subscript $L$ and $R$ denote the left and right (limited) states of the physical variables across the discontinuity.

Based on the above generalized Riemann invariants and Rankine-Hugoniot relations, the following theorem can be established.

Theorem 2.1

If $\gamma\in{(1,14.39]}$, then the following conclusions hold:

(i)* The left and right states across the 1-wave satisfy*

[TABLE]

(ii)* For the contact discontinuity, one has*

[TABLE]

(iii)* The left and right states across the 3-wave satisfy*

[TABLE]

• Proof

(i) For the 2-wave, i.e. contact discontinuity, because $p_{tot}$ and $u$ are two $2$-Riemann invariants, see (2.8), two identities in (2.11) hold.

(ii) If assuming that the 1-wave is the rarefaction wave, then the values of $\lambda_{1}$ at the left- and right-hand sides of the rarefaction wave satisfy

[TABLE]

and the second 1-Riemann invariant in (2.8) satisfies

[TABLE]

Combing them gives

[TABLE]

On the other hand, within the rarefaction wave, one has $S=S_{L}$ due to the first 1-Riemann invariant $S$ in (2.8), so that the total pressure $p_{tot}$ can be considered as a function of a single variable $\rho$ within this rarefaction wave, i.e. $p_{tot}=p_{tot}(\rho)$. Thus within the first rarefaction wave, $c(\rho,p_{tot})$ can also be expressed as a function of a single variable $\rho$, i.e. $c=c(\rho)$, and one has

[TABLE]

Substituting the above identity into (2.14) gives

[TABLE]

Because of the proof of Lemma 2.1, the inequality

[TABLE]

holds for any $\rho>0$ and $\gamma\in(1,14.39]\subset(1,15.95)$. Combing the above inequality with (2.15) gives

[TABLE]

Due to (2.3), $p_{tot}^{\prime}(\rho)>0$, which implies

[TABLE]

Moreover, the inequality $\rho_{L}>\rho_{R}$ and (2.13) imply

[TABLE]

(iii) For the $1$-shock wave, using the Rankine-Hugoniot conditions (2.9) gives

[TABLE]

and

[TABLE]

On the other hand, the Lax shock wave conditions

[TABLE]

gives that $\hat{u}_{L}>c_{L}>0$ and $0<\hat{u}_{R}<c_{R}$.

The proof by contradiction is considered here. Assume $p_{tot,L}\geq p_{tot,R}$. Using such assumption and the first and second identities in (2.9) gives

[TABLE]

Combing it with (2.18) yields

[TABLE]

Again using the first identity in (2.9) gives

[TABLE]

Substituting it into (2.17) derives

[TABLE]

The following attempts to find the contradiction to (2.20) in two cases.

Case (i): $\gamma_{1}>\frac{1}{3}$. Because $0<c_{L}<\hat{u}_{L}\leq\hat{u}_{R}<c_{R}$, and (2.16), one has

[TABLE]

Consider $\frac{c^{2}}{2}+h$ as a function of $\rho$ and $e_{tot}$, and then one has

[TABLE]

[TABLE]

where

[TABLE]

If $r_{e}<\frac{1}{8}$, then $g(r_{e},\gamma_{1})$ is always positive; otherwise, because

[TABLE]

the function $g(r_{e},\gamma_{1})$ is first increasing and then decreasing monotonically as $\gamma_{1}$ increases. Thanks to the fact that $g(r_{e},0)>0$, $g(r_{e},\gamma_{1})$ has only one positive zero with respect to $\gamma_{1}$. If there exists $\tilde{\gamma}_{1}$ and $r_{e}$ such that $g(r_{e},\tilde{\gamma}_{1})=0$, then for all $\gamma_{1}>\tilde{\gamma}_{1}$ and such $r_{e}$, $g(r_{e},\gamma_{1})<0$. Hence, one has

[TABLE]

Let us observe the maximum real zero of $g(r_{e},\gamma_{1})$ with respect to the independent variable $r_{e}$, and denote it by $r_{e}^{\max}$. If $r_{e}^{\max}\geq 0$, then there exist some $r_{e}=r_{e}^{\max}\geq 0$ and $\tilde{\gamma}_{1}$ such that $g(r_{e},\tilde{\gamma}_{1})=0$. It implies $\gamma_{1}\notin\mathcal{\tilde{R}}$. Fig. 2.2 shows the maximum real zero of $g(r_{e},\gamma_{1})$ with respect to $r_{e}$. It is seen that if $\bar{\gamma}_{1}\in(13.39,13.4)$, then $\frac{c^{2}}{2}+h$ is monotone increasing with respect to $\rho$ and $e_{tot}$ when $\frac{1}{3}<\gamma_{1}<\bar{\gamma}_{1}$, thus one has

[TABLE]

which conflicts with (2.19) or (2.20).

Case (ii): $0<\gamma_{1}\leq\frac{1}{3}$. The first equation of (2.9) and Lax conditions (2.18) gives

[TABLE]

If $\rho c$ is considered as a function of $\rho$ and $e_{tot}$, then one has

[TABLE]

and

[TABLE]

where

[TABLE]

Similarly, if $r_{e}<\frac{1}{8}$, then $\tilde{g}(r_{e},\gamma_{1})$ is positive; otherwise, because

[TABLE]

$\tilde{g}(r_{e},\gamma_{1})$ is first increasing and then decreasing monotonically with respect to $\gamma_{1}$. On the other hand, because

[TABLE]

$\tilde{g}(r_{e},\gamma_{1})>0$ for $0<\gamma_{1}\leq\frac{1}{3}$, and $\rho c$ is monotone increasing with respect to $\rho$ and $e_{tot}$. Thus one has

[TABLE]

which conflicts with (2.19) or (2.20).

In summary, one has

[TABLE]

Combing it with the first and second equations in (2.9) yields $u_{L}>u_{R}$.

(iv) The conclusion on the third characteristic field can be similarly verified. The proof is completed.

Before ending this section, the total differentials of the generalized Riemann invariants are given here. For smooth flows, the RHEs (2.1) can be equivalently written in the following form

[TABLE]

where $\frac{D}{Dt}=\frac{\partial}{\partial t}+u\frac{\partial}{\partial x}$ is the material derivative operator.

Because the third equation in (2.21) shows that the entropy $S$ is constant along the trajectory of each fluid particle, and the entropy $S$ is one of the i-Riemann invariants associated with the characteristic field $\lambda_{i}$, $i=1$ or $3$, see (2.8), all thermodynamic variables can be considered as functions of $\rho$ and $S$. With the help of such consideration, the total differentials of the i-Riemann invariants $\psi_{i}=\psi_{i}(u,{\rho},S)$ are

[TABLE]

where $\frac{\partial\psi_{i}}{\partial S}=\int^{\rho}\frac{1}{\omega}\frac{\partial c(\omega,S)}{\partial S}d\omega$, and

[TABLE]

with

[TABLE]

and

[TABLE]

which is obtained by (2.4), here the function ${\rm LambertW}(x)$ is the inverse of $x\exp(x)$, so-called the Lambert W function [5].

Remark 2.1

The integral $\int^{\rho}_{0}\frac{1}{\omega}\frac{\partial c(\omega,S)}{\partial S}d\omega$ in (2.23) can not be exactly evaluated, so that some numerical quadrature has to be considered. Our work uses the QAGS adaptive integration with singularities in the GNU Scientific Library.

3 Numerical schemes

This section gives the direct Eulerian GRP schemes for the one- and two-dimensional RHEs.

3.1 The outline of the one-dimensional GRP scheme

This section gives the outline of the one-dimensional direct Eulerian GRP scheme. For the sake of simplicity, the domain $\Omega$ is divided into a uniform mesh $\{x_{j}=j\Delta x,j\in\mathbb{Z}\}$ and $I_{j+\frac{1}{2}}=[x_{j},x_{j+1}]$ denotes the $(j+\frac{1}{2})$th cell. A (nonuniform) partition of the time interval $[0,T]$ is also given as $\{t_{0}=0,t_{n+1}=t_{n}+\Delta t_{n},n\geq 0\}$, where the time step size $\Delta t_{n}$ is determined by

[TABLE]

here \mbox{\boldmath\smallU}_{j+\frac{1}{2}}^{n} approximates the cell average value of of \mbox{\boldmath\smallU}(x,t_{n}) over the cell $I_{j+\frac{1}{2}}$ at time $t_{n}$, and $C_{cfl}$ denotes the CFL number.

If assuming that the (initial) data at time $t_{n}$ are piecewise linear functions with slopes (\mbox{\boldmath\smallU^{\prime}})_{j+\frac{1}{2}}^{n}, i.e.

[TABLE]

then the solution vector $U$ at time $t_{n+1}$ of (2.1) can be approximately evolved by using a second order accurate Godunov-type scheme

[TABLE]

where \mbox{\boldmath\smallU}_{j}^{n+1/2} denotes a second-order approximation of \mbox{\boldmath\smallU}(x_{j},t_{n}+\frac{1}{2}\Delta t_{n}), and is analytically derived by the second order accurate resolution of the local generalized Riemann problem (GRP) at each point $(x_{j},t_{n})$, i.e.

[TABLE]

More specifically, \mbox{\boldmath\smallU}_{j}^{n+1/2} is usually calculated by

[TABLE]

where the calculations of \left(\frac{\partial\mbox{\boldmath\smallU}}{\partial t}\right)_{j}^{n} is one of the key elements in the GRP scheme, and will be completed by resolving the problem (3.2), see Sections 3.2.2-3.2.5, and \mbox{\boldmath\smallU}_{j}^{RP,n}:=\mbox{\boldmath\small\omega}(0;\mbox{\boldmath\smallU}_{j,L}^{n},\mbox{\boldmath\smallU}_{j,R}^{n}), here \mbox{\boldmath\small\omega}(\frac{x-x_{j}}{t-t_{n}};\mbox{\boldmath\smallU}_{j,L}^{n},\mbox{\boldmath\smallU}_{j,R}^{n}) is the exact solution of the (classical) Riemann problem for (2.1) centered at $(x_{j},t_{n})$, i.e. the Cauchy problem

[TABLE]

and \mbox{\boldmath\smallU}_{j,L}^{n} and \mbox{\boldmath\smallU}_{j,R}^{n} are the left and right limiting values of \mbox{\boldmath\smallU}_{h}(x) at $(x_{j},t_{n})$. The exact solution of Riemann problem (3.4) will be discussed in Section 3.2.1.

After getting \mbox{\boldmath\smallU}_{j+\frac{1}{2}}^{n+1}, the approximate slope \mbox{\boldmath\smallU}^{\prime} at time $t_{n+1}$ is evolved by

[TABLE]

where \mbox{\boldmath\smallR}_{j+\frac{1}{2}}:=\mbox{\boldmath\smallR}(\mbox{\boldmath\smallU}_{j+\frac{1}{2}}^{n+1}), the parameter $\theta\in[1,2)$, and

[TABLE]

Here the limiter function is used to suppress possible numerical oscillations near the discontinuities, and \mbox{\boldmath\smallR}(\mbox{\boldmath\smallU}) is the right eigenvector matrix of the Jacobian matrix \frac{\partial\mbox{\boldmath\smallF}(\mbox{\boldmath\smallU})}{\partial\mbox{\boldmath\smallU}} as follows

[TABLE]

where $H:=\frac{E+p_{tot}}{\rho}$ denotes the total enthalpy.

3.2 Resolution of classical and generalized Riemann problems

This section derives \mbox{\boldmath\smallU}_{j}^{RP,n} and \left(\frac{\partial\mbox{\boldmath\smallU}}{\partial t}\right)_{j}^{n} in (3.3) by respectively resolving the local classical Riemann problem (3.4) and GRP (3.2), which can be changed into

[TABLE]

and

[TABLE]

by a coordinate transformation and ignoring the subscript $j$ in the initial data, where \mbox{\boldmath\smallU}_{L},\mbox{\boldmath\smallU}_{R},\mbox{\boldmath\smallU}^{\prime}_{L} and \mbox{\boldmath\smallU}^{\prime}_{R} are corresponding constant vectors.

The initial structure of the solution \mbox{\boldmath\smallU}^{GRP}(x,t) of (3.7) is determined by the solution \mbox{\boldmath\small\omega}(x/t;\mbox{\boldmath\smallU}_{L},\mbox{\boldmath\smallU}_{R}) of the Riemann problem (3.6), and

[TABLE]

The local wave configuration around the singularity point $(x,t)=(0,0)$ of the GRP (3.7) is usually piecewise smooth and consists of the elementary waves such as rarefaction wave, the shock wave and the contact discontinuity, as the schematic descriptions shown in Figs. 3.1 and 3.2.

The flow is isentropic for the rarefaction waves as a part of the solution \mbox{\boldmath\small\omega}(x/t;\mbox{\boldmath\smallU}_{L},\mbox{\boldmath\smallU}_{R}) of the Riemann problem (3.6), so the generalized Riemann invariants $\psi_{i}$ and $S$ are constant and their derivatives vanish within the $i$th-rarefaction wave, $i=1,3$. Unfortunately, those properties do not hold for the GRP (3.7) generally, because the general (curved) rarefaction waves have to be considered in the solution of the GRP (3.7). However, the solution of (3.7) can be regarded as a perturbation of the solution of (3.6) when $0\leq t\ll 1$, so that it can still be expected that $\psi_{i}$ and $S$ are regular within the $i$th-rarefaction wave as a part of the solution \mbox{\boldmath\smallU}^{GRP}(x,t) of (3.7) around the singularity point $(x,t)=(0,0)$, $i=1,3$, and the generalized Riemann invariants are used to resolve the rarefaction waves around the singularity point $(x,t)=(0,0)$.

From now on, our attention is restricted to the local wave configuration displayed in Figs. 3.1 and 3.2, where a rarefaction wave moves to the left and a shock moves to the right, the intermediate region is separated by a contact discontinuity, and the intermediate states in the two subregions are denoted by \mbox{\boldmath\smallU}_{1} and \mbox{\boldmath\smallU}_{2}, respectively. Other local wave configurations can be similarly considered. Finally, we denote by \mbox{\boldmath\smallU}_{*}, \left(\frac{\partial\mbox{\boldmath\smallU}}{\partial t}\right)_{*}, and \left(\frac{D\mbox{\boldmath\smallU}}{Dt}\right)_{*} the limiting states at $x=0$, as $t\rightarrow 0+$. The major feature of the direct Eulerian GRP scheme is to form a linear algebraic system

[TABLE]

by resolving respectively the left and right waves in Fig. 3.1. Solving this system gives the values of the quantities $\left(\frac{Du}{Dt}\right)_{*}$ and $\left(\frac{Dp_{tot}}{Dt}\right)_{*}$, and then gets $\left(\frac{\partial\rho}{\partial t}\right)_{*}$, $\left(\frac{\partial u}{\partial t}\right)_{*}$, $\left(\frac{\partial p_{tot}}{\partial t}\right)_{*}$ and closes the calculation in (3.3) and the scheme (3.1).

3.2.1 Exact solution of classical Riemann problem

This section discusses the exact solution of the Riemann problem (3.6) in order to get \mbox{\boldmath\smallU}_{*}. For the sake of convenience, $(\rho_{1*},u_{*},p_{tot,*})$ and $(\rho_{2*},u_{*},p_{tot,*})$ are used to denote the left and right states of the contact discontinuity in Fig. 3.2.

For the rarefaction and shock waves in Fig. 3.2, $u_{*}$ can be represented by $p_{tot,*}$ as follows

[TABLE]

and

[TABLE]

which are deduced from (2.8) and (2.9). The quantity $\rho_{1*}$ in the first equation of (3.9) can be represented in terms of $p_{tot,*}$ through the second equation in (3.9). Similarly, $\rho_{2*}$ in the first equation of (3.10) can also be represented in terms of $p_{tot,*}$ through the second equation in (3.10). In practice, both the second equations with a fixed $p_{tot,*}$ in (3.9) and (3.10) are numerically solved by the bisection method with a proper interval.

Remark 3.1

If the left wave is the shock wave, then (3.9) is replaced with

[TABLE]

If the right wave is the rarefaction wave, then (3.10) is replaced with

[TABLE]

Combining (3.9) with (3.10) gives a nonlinear equation with respect to $p_{tot,*}$ as follows

[TABLE]

The following theorem tells us when (3.11) has a solution.

Theorem 3.1

If $\gamma\in{(1,14.39]}$ and the initial data in (3.6) satisfy the positive pressure condition

[TABLE]

then (3.11) has an unique solution.

• Proof

Taking the partial derivative of (3.9) and (3.10) with respect to $p_{tot,*}$ gives

[TABLE]

and

[TABLE]

Thanks to Theorem 2.1, $\frac{d\rho_{*}(p_{tot,*})}{dp_{tot,*}}>0$ and $\rho_{2*}>\rho_{R}$ for the right shock wave, thus it holds

[TABLE]

For the right rarefaction wave and the left shock wave, the similar conclusions can be obtained. Hence $u_{L}(p_{tot,*})-u_{R}(p_{tot,*})$ is monotone decreasing.

It is not difficult to show that

[TABLE]

and the inequality

[TABLE]

holds for $p_{tot}>p_{tot,\max}:=\max\{p_{tot,L},p_{tot,R}\}$. The last inequality implies

[TABLE]

Thus, if $u_{L}(p_{tot,\min})-u_{R}(p_{tot,\min})>0$ and $u_{L}(p_{tot,\max})-u_{R}(p_{tot,\max})<0$, where $p_{tot,\min}:=\min\{p_{tot,L},p_{tot,R}\}$, then there exists an unique solution $p_{tot,*}\in(p_{tot,\min},p_{tot,\max})$; if $u_{L}(p_{tot,\max})-u_{R}(p_{tot,\max})>0$, then there exists an unique solution $p_{tot,*}\in(p_{tot,\max},+\infty)$ due to (3.14); if $u_{L}(p_{tot,\min})-u_{R}(p_{tot,\min})<0$, then there exists an unique solution $p_{tot,*}\in(0,p_{tot,\min})$ due to (3.13). The proof is completed.

In practice, (3.11) is solved by using the Newton’s iteration, and then $\rho_{1*}$, $\rho_{2*}$ and $u_{*}$ can be gotten with the help of the generalized Riemann invariants and the Rankine-Hugoniot conditions. Finally, \mbox{\boldmath\smallU}_{*}=\mbox{\boldmath\smallU}_{1} if $u_{*}>0$; otherwise \mbox{\boldmath\smallU}_{*}=\mbox{\boldmath\smallU}_{2}.

For the sonic case, i.e. the $t$-axis is within the left rarefaction region. As $t\rightarrow 0+$ along $t$-axis, $u_{*}=c_{*}$. The limiting value \mbox{\boldmath\smallU}_{*} can be simply gotten as follows. The equation

[TABLE]

is first solved to obtain $\rho_{*}$ by the bisection method with the interval $[\rho_{1*},\rho_{L}]$ and then calculate $u_{*}=c(\rho_{*},S_{L})$ and $p_{tot,*}=p_{tot}(\rho_{*},S_{L})$.

Based on the above discussion, the exact solution of the Riemann problem (3.6) can be completely obtained by mimicking the procedure for the exact solution of the Riemann problem of the Euler equations in gas dynamics in [20, Chapter 4].

3.2.2 Resolution of rarefaction waves

This section resolves the left rarefaction wave in the solution of local GRP (3.7) as shown in Fig. 3.1 and gives the first equation in (3.8). Before that, some similar notations are introduced to those in [4, 21, 22, 23, 24, 25].

The region of the left rarefaction wave for (3.7) can be described by the set $\mathcal{R}=\{(\alpha(x,t),\beta(x,t))|\beta\in[\beta_{L},\beta_{*}],-\infty<\alpha<0\}$, where \beta_{L}=\lambda_{1}(\mbox{\boldmath\smallU}_{L}), \beta_{*}=\lambda_{1}(\mbox{\boldmath\smallU}_{*}), and $\beta=\beta(x,t)$ and $\alpha=\alpha(x,t)$ are the integral curves of the following equations

[TABLE]

respectively. Here $\beta$ and $\alpha$ have been denoted as follows: $\beta$ is the initial value of the slope $\lambda_{1}$ at the singularity point $(0,0)$, and $\alpha$ for the transversal characteristic curves $\frac{dx}{dt}=\lambda_{3}$ is the $x$-coordinate of the intersection point with the leading $\beta$-curve, see Fig. 3.1. Similarly, a local coordinate transformation is also introduced within the region of the left rarefaction wave in Fig. 3.2. To avoid confusion, it is denoted as $x_{\rm ass}=x({\alpha,\beta})$ and $t_{\rm ass}=t({\alpha,\beta})$.

The major result of this section is given in the following lemma.

Lemma 3.1

The limiting values of $\frac{Du}{Dt}$ and $\frac{Dp_{tot}}{Dt}$ satisfy

[TABLE]

with

[TABLE]

[TABLE]

where the expressions of $A_{L}(0,\beta)$, $\frac{\partial t}{\partial\alpha}(0,\beta)$, and $B_{L}(0,\beta)$ are given as follows

[TABLE]

[TABLE]

and

[TABLE]

here $S^{\prime}_{L}=\left(\frac{\partial S}{\partial x}\right)_{L}$, $\psi_{1,L}^{\prime}=-\frac{1}{\rho_{L}c_{L}}\left(\frac{\partial p_{tot}}{\partial x}\right)_{L}+\left(\frac{\partial u}{\partial x}\right)_{L}+K(\rho_{L},S_{L})S^{\prime}_{L}$, and $\beta=\lambda_{1}$.

• Proof

Following the proof of Lemma 3.1 in [4], one can obtain (3.18), (3.19), and

[TABLE]

[TABLE]

[TABLE]

Using

[TABLE]

gives (3.22) and

[TABLE]

and then derives (3.20) and (3.21). The proof is completed.

Remark 3.2

For the right rarefaction wave, corresponding coefficients $a_{R}$, $b_{R}$, and $d_{R}$ in (3.8) can be similarly given by

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

here $S^{\prime}_{R}=\left(\frac{\partial S}{\partial x}\right)_{R}$, $\psi_{3,R}^{\prime}=-\frac{1}{\rho_{R}c_{R}}\left(\frac{\partial p_{tot}}{\partial x}\right)_{R}+\left(\frac{\partial u}{\partial x}\right)_{R}-K(\rho_{R},S_{R})S^{\prime}_{R}$, and $\beta=\lambda_{3}$.

Remark 3.3

The integral in (3.22) or (3.27) can be exactly evaluated. In our computations, the three point Gauss-Lobatto quadrature is used only with an additional calculation the physical states at the internal point $\frac{1}{2}(c_{L}+c_{1,*})$ of the integral interval. Moreover, if $\hat{a}_{R}=0$, the above corresponding coefficients $a_{L}$, $b_{L}$, and $d_{L}$ (or $a_{R}$, $b_{R}$, and $d_{R}$) are the same as those in [4].

3.2.3 Resolution of shock waves

This section resolves the right shock wave in the solution for the local GRP (3.7) as shown in Fig. 3.1 and gives the second equation in (3.8) through differentiating the shock relation along the shock trajectory.

Let $x=x(t)$ be the shock trajectory associated with the third characteristic field and assume that it propagates with the speed $\sigma:=x^{\prime}(t)>0$ to the right, see Fig. 3.1. Denote the left and right states of the shock wave by \mbox{\boldmath\smallU}(t) and \mbox{\boldmath\small\bar{U}}(t), respectively, i.e. \mbox{\boldmath\smallU}(t)=\mbox{\boldmath\smallU}(x(t)-0,t) and \mbox{\boldmath\small\bar{U}}(t)=\mbox{\boldmath\smallU}(x(t)+0,t). Across the shock wave, using the Rankine-Hugoniot jump conditions gives

[TABLE]

where

[TABLE]

and the density $\rho$ is implicitly defined by the equation

[TABLE]

that is to say, $\rho$ can be regarded as a function of $p_{tot},\bar{p}_{tot},\bar{\rho}$, denoted by $\rho=H(p_{tot},\bar{p}_{tot},\bar{\rho})$.

Along the shock trajectory, one has

[TABLE]

where $D_{\sigma}:=\frac{\partial}{\partial t}+\sigma\frac{\partial}{\partial x}$ denotes the directional derivative operator along the shock trajectory.

The major result in this section is given as follows.

Lemma 3.2

The limiting values of $\frac{Du}{Dt}$ and $\frac{Dp_{tot}}{Dt}$ satisfy

[TABLE]

with

[TABLE]

where

[TABLE]

and $\sigma_{R}=\frac{\rho_{2*}u_{*}-\rho_{R}u_{R}}{\rho_{2*}-\rho_{R}}$. Here $\Phi_{i}$, $i=1,2,3$, are defined by

[TABLE]

and

[TABLE]

[TABLE]

where

[TABLE]

• Proof

It is similar to the proof of Lemma 4.1 of [4].

Remark 3.4

For the left shock wave, corresponding coefficients $a_{L}$, $b_{L}$, $d_{L}$ in (3.8) can be similarly derived as follows

[TABLE]

with

[TABLE]

where $\sigma_{L}=\frac{\rho_{1*}u_{*}-\rho_{L}u_{L}}{\rho_{1*}-\rho_{L}}$.

3.2.4 Time derivatives of solutions at singularity point

This section solves the $2\times 2$ linear system (3.8) formed by (3.17) and (3.32) to give the total derivatives $\left(\frac{Du}{Dt}\right)_{*}$ and $\left(\frac{Dp_{tot}}{Dt}\right)_{*}$, and further derives the time derivatives $\left(\frac{\partial u}{\partial t}\right)_{*}$, $\left(\frac{\partial p_{tot}}{\partial t}\right)_{*}$ and $\left(\frac{\partial\rho}{\partial t}\right)_{*}$ at the singularity point $(x,t)=(0,0)$.

Case 1: Nonsonic case.

It is necessary to check whether such linear system has an unique solution.

Theorem 3.2

The linear system (3.17) and (3.32) has an unique solution.

• Proof

In order to prove the conclusion, it needs to prove

[TABLE]

For the left (rarefaction) wave, (3.18) implies

[TABLE]

For the right (shock) wave, because $\rho_{2*}>\rho_{R}$, and $p_{tot,*}>p_{tot,R}$, one has

[TABLE]

where

[TABLE]

[TABLE]

Using $\sigma_{R}>u_{*}$ and (3.33) gives

[TABLE]

In summary, one has

[TABLE]

so that the linear system (3.17) and (3.32) has an unique solution.

The time derivatives of solutions at singularity point can be calculated as follows.

Theorem 3.3

Assuming that the $t$ axis is not included in the rarefaction wave, then the limiting values of time derivatives \left(\frac{\partial\mbox{\boldmath\smallW}}{\partial t}\right)_{*} can be calculated by

[TABLE]

and

[TABLE]

or

[TABLE]

where $A_{L}(0,\beta_{*})$ and $K(\rho_{*},S_{*})$ are given in (3.20) and (2.23), respectively, and $g_{u}^{R},g_{p_{tot}}^{R}$, and $g^{R}_{\rho}$ are constant, depending on \mbox{\boldmath\smallU}_{R} and \mbox{\boldmath\smallU}_{*}, and their expressions are given by

[TABLE]

Here $H_{i}$ is defined in (3.38), $i=1,2,3$.

Case 2:: Sonic case.

When the $t$ axis ($x=0$) is located in the rarefaction fan, the sonic case happens and Theorems 3.2 and 3.3 are not available. However, since one of the characteristic curves is tangential to the $t$-axis, the sonic case becomes much simpler.

Theorem 3.4

If the $t$ axis is within the left rarefaction wave, then one has

[TABLE]

[TABLE]

where $A_{L}(0,\beta)$, $d_{L}(\beta)$ and $K(\rho_{*},S_{*})$ are given in (3.20), (3.19) and (2.23), respectively.

Remark 3.5

The proofs of Theorems 3.3 and 3.4 are similar to those of several theorems in [4, Section 5].

3.2.5 Acoustic case

This section turns to the acoustic case of the GRP (3.7), i.e. \mbox{\boldmath\smallU}_{L}=\mbox{\boldmath\smallU}_{R} and \mbox{\boldmath\smallU}^{\prime}_{L}\neq\mbox{\boldmath\smallU}^{\prime}_{R}. In this case, \mbox{\boldmath\smallU}_{L}=\mbox{\boldmath\smallU}_{R}=\mbox{\boldmath\smallU}_{*} and only linear waves emanate from the origin (0, 0) so that the resolution of the GRP becomes simpler than the general case discussed before. For the sake of simplicity, in the following we only reserve the subscripts related to the state \mbox{\boldmath\smallU}_{L} or \mbox{\boldmath\smallU}_{*} or \mbox{\boldmath\smallU}_{R} for the derivatives of all physical variables.

Theorem 3.5

In the acoustic case of that \mbox{\boldmath\smallU}_{L}=\mbox{\boldmath\smallU}_{R} and \mbox{\boldmath\smallU}^{\prime}_{L}\neq\mbox{\boldmath\smallU}^{\prime}_{R}, if $u_{*}-c_{*}<0$ and $u_{*}+c_{*}>0$, then $\left(\frac{\partial u}{\partial t}\right)_{*},\left(\frac{\partial p_{tot}}{\partial t}\right)_{*}$ and $\left(\frac{\partial\rho}{\partial t}\right)_{*}$ can be obtained by

[TABLE]

and

[TABLE]

• Proof

Mimicking the proof of Theorem 6.1 of [4] can complete the proof.

3.3 Two-dimensional extension

This section extends the previous GRP scheme to the two-dimensional RHEs

[TABLE]

with

[TABLE]

where the physical meanings of $\rho,e_{tot}$, and $p_{tot}$ are the same as those in the one-dimensional case, $(u,v)$ denotes the velocity vector, and the total energy $E$ is given by $E=\frac{1}{2}\rho(u^{2}+v^{2})+\rho e_{tot}$.

The extension will be with the help of the Strang splitting technique, see [16, 17] and [2, Chapter 7]. First, the (3.49) is split into two subsystems such as

[TABLE]

Then the second-order accurate Strang splitting method can be given by

[TABLE]

where $\mathcal{L}_{x}(\Delta t)$ and $\mathcal{L}_{y}(\Delta t)$ denote the one-dimensional evolution operators for one time step $\Delta t$ of the above subsystems in (3.50), respectively.

The rest of the tasks only replace $\mathcal{L}_{x}(\Delta t)$ and $\mathcal{L}_{y}(\Delta t)$ with corresponding one-dimensional GRP evolution operators. Without loss of generality, take the $x$-split system in (3.50) as an example. In our GRP scheme, the term \mathcal{L}_{x}(\Delta t)\mbox{\boldmath\smallU}^{n} is almost the right-hand side of the scheme (3.1). In fact, in the $x$-split system, the third equation can be simplified as follows

[TABLE]

and decoupled from the $x$-split system, and the other three equations are the same as those in the one-dimensional RHEs (2.1). Thus, both the 1- and 4-waves of $x$-split system are nonlinear, and the 2- and 3-waves are contact discontinuity and shear wave with the speed of $u$, respectively. The previous resolution of one-dimensional GRP can be used to derive the values of $\rho_{*}$, $u_{*}$, and $e_{*}$ as well as the limiting values of their first order derivatives at $x=0$, as $t\rightarrow 0^{+}$, while the value of $v_{*}$ is calculated as follows

[TABLE]

and the calculation of $\left(\frac{\partial v}{\partial t}\right)_{*}$ can be similarly done in the way used in [4]. For the case that the 1- and 4-waves are the left-moving rarefaction wave and right-moving shock wave, respectively, and the $t$-axis is located in the intermediate region, see Fig. 3.1, it is done as follows.

Theorem 3.6

For the $x$-split system of (3.49) in the case that the 1- and 4-waves are the left-moving rarefaction wave and right-moving shock wave, respectively, and the $t$-axis is located in the intermediate region, one has:

(i) For the case of $u_{*}>0$, $(\partial v/\partial t)_{*}$ is calculated as follows

[TABLE]

(ii) For the case of $u_{*}<0$, $(\partial v/\partial t)_{*}$ is calculated as follows

[TABLE]

• Proof

Using Theorem 7.1 of [4] can complete the proof.

It is worth noting that the previous results are still available if the $t$ axis is within the rarefaction wave. Specially, if the $t$ axis is within the left rarefaction wave, then $\left(\frac{\partial v}{\partial t}\right)_{*}$ can be obtained by (3.52).

4 Numerical experiments

This section will solve several initial value or initial-boundary-value problems of the one- and two-dimensional RHEs (2.1) and (3.49) to verify the accuracy and the discontinuity resolving capability of the proposed second order accurate GRP schemes. Unless specifically stated, the adiabatic index $\gamma$ and radiation coefficient $\hat{a}_{R}$ are chosen as $5/3$ and $1$, respectively, and the CFL number $C_{cfl}$ and the limiter parameter $\theta$ in (3.5) are taken as $0.45$ and $1.5$, respectively.

4.1 One-dimensional case

Five one-dimensional examples are considered here. The first is used to check the accuracy of the one-dimensional GRP scheme and the others are the Riemann problems considered in [6, 18] and used to validate the performance of the GRP scheme in resolving the discontinuities. The numerical solutions obtained by the GRP scheme will be drawn in the symbols “$\circ$” and compared to the solutions of the MUSCL-Hancock scheme [20, Chapter 14.4] (in the symbol “+”) and the exact or reference solutions in the solid lines. The computational domain $\Omega$ is chosen as $[0,1]$.

Example 4.1** (Accuracy test)**

The first example is used to check the accuracy of the one-dimensional GRP scheme for the smooth solution of RHEs

[TABLE]

which describes a sine wave propagating periodically within the domain $\Omega=[0,1]$.

The computational domain $\Omega$ is divided into $N$ uniform cells and the periodic boundary conditions are specified at $x=0,1$. Table 4.1 gives the $l^{p}$ errors in the density at $t=0.5$ and corresponding convergence rates of the GRP scheme, where $p=1,2,\infty$. The data show that the convergence rates of second-order can be almost obtained in $l^{1}$ and $l^{2}$ norms, but the $l^{\infty}$ convergence rates is lower than 1.5 due to the nonlinear limiter (3.5).

Example 4.2** (Riemann problem 1)**

The initial data of this problem are given as follows

[TABLE]

The exact solutions at time $t>0$ involve two shock waves and a contact discontinuity. Fig. 4.1 gives the numerical solutions at $t=0.04$ obtained by using the GRP and MUSCL-Hancock schemes with 200 uniform cells with $\theta=1$. It is seen that the numerical solutions are in good agreement with the exact and the GRP scheme captures the shock waves and contact discontinuity better than the MUSCL-Hancock scheme.

Example 4.3** (Riemann problem 2)**

It is similar to the first Riemann problem, but with much stronger shock waves and larger initial velocity jump. The initial data are given by

[TABLE]

Fig. 4.2 gives the numerical solutions at $t=0.018$ obtained by using the GRP and MUSCL-Hancock schemes with 200 uniform cells with $\theta=1$. It is seen that the numerical solutions are in good agreement with the exact solutions, and the GRP scheme does work well for very strong shock waves and captures them and the contact discontinuity better than the MUSCL-Hancock scheme.

Example 4.4** (Riemann problem 3)**

The initial data of the third one-dimensional Riemann problem are taken as

[TABLE]

Fig. 4.3 gives the numerical solutions at $t=0.2$ obtained by using the GRP and MUSCL-Hancock schemes. It is seen that the solutions contain the left- and right-moving rarefaction waves, there is a undershoot (resp. overshoot) “bubble” in the computed density (resp. temperature) due to the numerical wall-heating, and the GRP scheme resolve those rarefaction waves much better than the MUSCL-Hancock scheme and produces the less serious wall-heating phenomenon at $x=0.5$.

Example 4.5** (Riemann problem 4)**

The initial data of the last one-dimensional problem are

[TABLE]

Fig. 4.4 displays the solutions at $t=0.01$ obtained by using the second order GRP and MUSCL-Hancock schemes with 200 uniform cells. It is shown that the solutions involve a left-moving rarefaction wave, a contact discontinuity, and a right-moving shock wave, the GRP scheme resolves the narrow region between the contact discontinuity and shock wave much better than the MUSCL-Hancock scheme, the obvious overshoot is observed near the tail of the numerical rarefaction waves.

4.2 Two-dimensional case

Three two-dimensional examples are considered here. The first is used to check the accuracy of the two-dimensional GRP scheme and the others are the shock and cloud interaction problem and the wind and cloud interaction problem considered in [6, 18] and used to validate the performance of the GRP scheme in resolving the discontinuities and complex wave structure.

Example 4.6** (Accuracy test)**

It is used to check the accuracy of our GRP scheme for the smooth solution of two-dimensional RHEs (3.49)

[TABLE]

which describes a sine wave propagating periodically at an angle of $\frac{\pi}{4}$ relative to the $x$-axis in the domain $\Omega=[0,1]\times[0,1]$ with the periodic boundary conditions.

Table 4.2 lists the $l^{p}$ errors in the density $\rho$ at $t=0.5$ and convergence rates of the GRP scheme, where $p=1,2,\infty$, and $\Omega$ is divided into $N\times N$ uniform rectangular cells. It is shown that the $l^{1}$ and $l^{2}$ convergence rates of GRP scheme are close to $2$, but the $l^{\infty}$ convergence rates is lower than 1.5 due to the nonlinear limiter (3.5).

Example 4.7** (Shock and cloud interaction)**

This problem is about the interaction between a shock wave and a denser cloud in the domain $\Omega=[0,2]\times[0,1]$ with the boundary conditions of zero gradient for the flow variables. Initially, there is a strong planer shock wave with the Mach number $256.7$ at $x_{0}=1.64$, which propagates in the $x$-direction. Its pre- and post-states are given by

[TABLE]

The circular cloud is 100 denser than the pre-shock state, and its radius $R$ and center are $0.15$ and $(x_{0}+0.18,0.5)$, respectively. Moreover, $\hat{a}_{R}=0.01$.

Figs. 4.5 and 4.6 display the contour plots of density and temperature at $t=0.07$ obtained by using the present GRP scheme with $256\times 128$ uniform cells, and Fig. 4.7 draws the density and temperature at $y=0.5$, in comparison with those obtained by using the MUSCL-Hancock scheme. It is seen that they are obviously better than those obtained by using the MUSCL-Hancock scheme and the second-order accurate high resolution KFVS scheme [18] and second order BGK scheme [9].

Example 4.8** (Wind and cloud interaction)**

It is about the interaction between a wind and a denser cloud within the computational domain $\Omega=[0,2]\times[0,1]$. Initially, there is a circular cloud of radius $R=0.15$ in $\Omega$, whose center is located at $(x_{0},y_{0})=(0.3,0.5)$. The cloud is 25 times denser than the ambient gas with the state $(\rho,T,u,v)=(1,0.09,0,0)$. The wind is introduced through the left boundary, at which the state (\rho,T,u,v)(0,y,t)=\big{(}1,0.09,6(1-e^{-10t}),0\big{)} is always assigned. The right, upper and lower boundary conditions are zero gradient for the flow variables. Specially, the initial data in $\Omega$ are given by

[TABLE]

Figs. 4.8 and 4.9 show the contour plots of density and temperature obtained by using the GRP scheme with $256\times 128$ uniform rectangular cells at $t=0.6$, respectively, and Fig. 4.10 plots the density and temperature at $y=0.5$, in comparison with those obtained by using the MUSCL-Hancock scheme. The interaction between the wind and denser cloud is captured much better than those resolved by the MUSCL-Hancock scheme and the high resolution KFVS and second order BGK schemes in the literature.

5 Conclusions

The paper developed the second-order accurate direct Eulerian generalized Riemann problem (GRP) scheme for the radiation hydrodynamical equations (RHE) in the zero diffusion limit. The main difficulty came from no explicit expression of the flux in terms of the conservative vector due to the nonlinearity in the radiation pressure and energy. The characteristic fields and the relations between the left and right states across the elementary-waves were first studied, and then the exact solution of the one-dimensional Riemann problem was analyzed and given. The initial data reconstruction in space were done for the characteristic variables. Based on those, the direct Eulerian GRP scheme was derived by directly using two main ingredients, the generalized Riemann invariants and the Runkine-Hugoniot jump conditions, to analytically resolve the left and right nonlinear waves of the local GRP in the Eulerian formulation and to obtain the limiting values of the time derivatives of the conservative variables along the cell interface and the numerical flux for the GRP scheme. Several one- and two-dimensional numerical experiments were conducted to demonstrate the accuracy and resolution of the proposed GRP schemes, in comparison with the second-order accurate MUSCL-Hancock scheme.

Acknowledgements

This work was partially supported by the Special Project on High-performance Computing under the National Key R&D Program (No. 2016YFB0200603), Science Challenge Project (No. JCKY2016212A502), and the National Natural Science Foundation of China (Nos. 91330205, 91630310, 11421101). The authors would also like to thank Mr. Pu Chen for his discussion.