A Fourth Order Entropy Stable Scheme for Hyperbolic Conservation Laws

Abstract

This paper develops a fourth order entropy stable scheme to approximate the entropy solution of one-dimensional hyperbolic conservation laws. The scheme is constructed by employing a high order entropy conservative flux of order four in conjunction with a suitable numerical diffusion operator that based on a fourth order non-oscillatory reconstruction which satisfies the sign property. The constructed scheme possesses two features: (1) it achieves fourth order accuracy in the smooth area while keeping high resolution with sharp discontinuity transitions in the nonsmooth area; (2) it is entropy stable. Some typical numerical experiments are performed to illustrate the capability of the new entropy stable scheme.

1. Introduction

During the past few decades, abundant numerical methods for solving hyperbolic conservation laws have been designed; one can consult the review papers [1,2] and the references therein. Among the various methods, high order schemes, such as total variation diminishing (TVD) schemes, weighted essentially non-oscillatory (WENO) schemes and discontinuous Galerkin (DG) schemes, have achieved great success. However, it seems that there are few entropy stability results for high order numerical schemes for hyperbolic conservation laws, especially for the nonlinear hyperbolic conservation laws. In view of the above discussion, we limit our research on high order entropy stable methods for solving hyperbolic conservation laws.

In accordance with the idea of Tadmor, entropy schemes are often constructed by utilizing entropy conservative flux in conjunction with suitable numerical diffusion [3]. For example, physical viscosity was discretized and used as numerical diffusion in [4]. Roe-type numerical diffusion was selected for various systems like Euler equations, shallow water equations and ideal magnetohydrodynamic equations [5,6,7]. Based on the limiter mechanism, Liu et al. constructed a family of entropy consistent schemes with high resolution [8]. Recently, Dubey discussed the amount of suitable diffusion for the sake of devising non-oscillatory entropy stable schemes in the TVD sense [9]. To obtain high order entropy stable schemes, Fjordholm et al. firstly proved that ENO reconstruction satisfies a sign property and presented entropy stable schemes based on ENO reconstruction of arbitrary order accuracy [10,11]. To overcome the drawbacks of ENO reconstruction, we presented a third order reconstruction which is non-oscillatory and satisfies the sign property to construct a third order entropy stable scheme [12]. WENO reconstruction was modified to satisfy the sign property and also used to construct an entropy stable scheme in [13,14], although they are limited to third order accuracy.

In this paper, we are aiming at presenting a new entropy stable scheme of fourth order accuracy to solve hyperbolic conservation laws in one dimension. First, we employ a fourth order entropy conservative flux which is a linear combination of two-point entropy conservative fluxes. Then, we present a non-oscillatory reconstruction of fourth order accuracy which possesses the sign property to obtain a fourth order accurate numerical diffusion operator. By adding this numerical diffusion operator to the entropy conservative flux, the resulting flux is entropy stable.

The remainder of this paper is given as follows. Section 2 describes the procedure for building our fourth order entropy scheme. After that, we present some typical numerical experiments to show the effectiveness of the newly developed scheme. Finally, concluding remarks are given in Section 4.

2. Numerical Method

Consider the scalar conservation law

$$\frac{\partial }{\partial t}u+\frac{\partial }{\partial x}f\left(u\right)=0.$$

(1)

The weak solutions for Equation (1) are not unique and we are interested in the so-called entropy solution, which satisfies the entropy condition

$$\frac{\partial }{\partial t}\eta \left(u\right)+\frac{\partial }{\partial x}q\left(u\right)\le 0.$$

(2)

Here, the entropy function $\eta \left(u\right)$ is convex and the entropy flux function $q\left(u\right)$ satisfies $\frac{\partial q}{\partial u}=v\frac{\partial f}{\partial u}$ with the entropy variable $v=\frac{\partial \eta }{\partial u}$. To solve Equation (1) with the conservative difference method, we have the semi-discrete scheme

$$\frac{d{u}_i}{dt}=-\frac{f_{i+1/2}-f_{i-1/2}}{\Delta x},$$

(3)

where $u_i$ represents the point valve at the node $x_i$ on a uniform Cartesian mesh with the mesh size $x_{i+1}-x_i=\Delta x$. The scheme (3) is said to be entropy stable if it satisfies a discrete version of the entropy condition (2), namely,

$$\frac{\partial }{\partial t}\eta \left(u_i\right)+\frac{q_{i+1/2}-q_{i-1/2}}{\Delta x}\le 0;$$

(4)

the scheme is said to be entropy conservative if it satisfies the discrete entropy equality

$$\frac{\partial }{\partial t}\eta \left(u_i\right)+\frac{q_{i+1/2}-q_{i-1/2}}{\Delta x}=0.$$

(5)

In this work, we utilize the familiar Runge–Kutta scheme of order four to solve the ordinary Equation (3). Now, we focus on how to build the entropy stable numerical flux $f_{i+1/2}$.

First, the numerical flux is split into an entropy conservative part and a numerical diffusion part

$$f_{i+1/2}=f_{i+1/2}^{EC}-\frac{1}{2}{\delta }_{i+1/2}\left(v_{i+1/2}^{+}-v_{i+1/2}^{-}\right).$$

(6)

Here, $f_{i+1/2}^{EC}$ is a fourth order entropy conservative flux [15] expressed as

$$f_{i+1/2}^{EC}=\frac{4}{3}\tilde{f}\left(u_i,u_{i+1}\right)-\frac{1}{6}\left(\tilde{f}\left(u_{i-1},u_{i+1}\right)-\tilde{f}\left(u_i,u_{i+2}\right)\right)$$

(7)

with $\tilde{f}(\alpha ,\beta )$ being the second order entropy conservative flux defined by Tadmor [3]. Generally, ${\delta }_{i+1/2}$ in the numerical diffusion part is chosen to be $\left|f^{\prime }\left(u_{i+1/2}\right)\right|$. If the reconstructed values $v_{i+1/2}^{\pm }$ at the interfaces from the entropy variable $v_i$ satisfy

$$sign(v_{i+1/2}^{+}-v_{i+1/2}^{-})=sign(v_{i+1}-v_i),$$

(8)

the numerical flux (6) achieves entropy stability. Property (8) in which the jump of the reconstructed point values at each cell interface has the same sign as the jump of the underlying point values across that interface is called the sign property. We present the fourth order reconstruction satisfying the sign property as follows to obtain $v_{i+1/2}^{\pm }$.

Consider a polynomial reconstruction of the form

$$p_i\left(x\right)=v_i+d_i\left(\frac{x-x_i}{\Delta x}\right)+\left(\frac{v_{i-1}-2v_i+v_{i+1}}{2}\right){\left(\frac{x-x_i}{\Delta x}\right)}^2+\left(\frac{-v_{i-1}+v_{i+1}-2d_i}{2}\right){\left(\frac{x-x_i}{\Delta x}\right)}^3$$

(9)

with $d_i$ being a slope function. For convenience, we introduce the following notation

$$\begin{array}{c}d{s}_i=\frac{2}{3}{v}_{i+1}-\frac{2}{3}{v}_{i-1}-\frac{1}{12}{v}_{i+2}+\frac{1}{12}{v}_{i-2},\hfill \\ W{C}_i=v_{i+1}-v_{i-1},W{R}_i=v_{i+1}-v_i,WC{2}_i=v_{i+2}-v_{i-2},\hfill \\ ds{1}_i=0,ds{2}_i=\frac{1}{2}\left(W{C}_i-8W{R}_i\right),ds{3}_i=\frac{1}{2}\left(8W{R}_i-7W{C}_i\right),\hfill \\ S_i=sign\left(W{C}_i\right),\hfill \\ C_1=\frac{\sqrt{3}}{6},C_2=\frac{6}{12+\sqrt{3}},\hfill \end{array}$$

(10)

and define the slope $d_i$ in the following way:

• If $S_i=0$, then $d_i=0$.
• If $S_i\ne 0$ and $\left(2S_i·W{C}_i\ge S_i·WC{2}_i\right)$, then $d_i=d{s}_i$.
• If $S_i\ne 0$ and $\left(2S_i·W{C}_i<S_i·WC{2}_i\right)$, then the following hold:

  (a)

  If $v_i=\frac{v_{i+1}+v_{i-1}}{2}$, we define

  $$d_i=\left\{\begin{array}{cc}max(ds{1}_i,d{s}_i),\hfill & \hfill \mathrm{if}{S}_i>0,\\ min(ds{1}_i,d{s}_i),\hfill & \hfill \mathrm{if}{S}_i<0.\end{array}\right.$$

  (b)

  If $v_i\ne \frac{v_{i+1}+v_{i-1}}{2}$, then the following hold:

  i.

  If $\left|W{R}_i-\frac{1}{2}W{C}_i\right|\ge \frac{1}{8}\left|WC{2}_i-2W{C}_i\right|$, then

  $$d_i=\left\{\begin{array}{cc}max(ds{2}_i,ds3,d{s}_i),\hfill & \hfill \mathrm{if}{S}_i>0,\\ min(ds{2}_i,ds3,d{s}_i),\hfill & \hfill \mathrm{if}{S}_i<0.\end{array}\right.$$

  ii.

  If $\left|W{R}_i-\frac{1}{2}W{C}_i\right|<\frac{1}{8}\left|WC{2}_i-2W{C}_i\right|$, then

  $$d_i=\left\{\begin{array}{cc}\frac{W{C}_i}{2}-S_i·C1·\left|2W{R}_i-W{C}_i\right|,\hfill & \hfill \mathrm{if}\left|\frac{W{R}_i}{W{C}_i}-\frac{1}{2}\right|\le C2,\\ \frac{W{C}_i}{2},\hfill & \hfill \mathrm{if}\left|\frac{W{R}_i}{W{C}_i}-\frac{1}{2}\right|>C2.\end{array}\right.$$

Theorem 1.

With this definition of the slope, the polynomial reconstruction (9) fulfills the shape preserving property, namely,

• $p_i\left(x\right)$ is monotone in $\left[x_{i-1/2},x_{i+1/2}\right]$ iff the point values $\{v_{i-1},v_i,v_{i+1}\}$ are;
• $p_i\left(x\right)$ generates an extremum in the interior of $\left[x_{i-1/2},x_{i+1/2}\right]$ iff $v_i$ is an extremum value.

The proof of this theorem can be found in [16]. According to this property, the polynomial $p_i\left(x\right)$ does not create new extrema in the interior of $\left[x_{i-1/2},x_{i+1/2}\right]$. In other words, the polynomial $p_i\left(x\right)$ may create new extrema at the interfaces $\left\{x_{i\pm 1/2}\right\}$. To keep away from such spurious extrema, we can employ a similar strategy as in [17]. Let us modify the polynomial reconstruction to be the following form:

$${\phi }_i\left(x\right)=(1-{\theta }_i)v_i+{\theta }_i{p}_i\left(x\right),0<{\theta }_i<1.$$

(11)

The limiter ${\theta }_i$ is determined by

$${\theta }_i=\left\{\begin{array}{cc}\min \left\{\frac{M_{i+1/2}-v_i}{M_i-v_i},\frac{m_{i-1/2}-v_i}{m_i-v_i},1\right\},\hfill & \mathrm{if}{v}_{i-1}<v_i<v_{i+1},\hfill \\ \min \left\{\frac{M_{i-1/2}-v_i}{M_i-v_i},\frac{m_{i+1/2}-v_i}{m_i-v_i},1\right\},\hfill & \mathrm{if}{v}_{i-1}>v_i>v_{i+1},\hfill \\ 1,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(12)

with

$$M_i=\underset{x\in [x_{i-1/2},x_{i+1/2}]}{\max }{p}_i\left(x\right),m_i=\underset{x\in [x_{i-1/2},x_{i+1/2}]}{\min }{p}_i\left(x\right),$$

and

$$\left\{\begin{array}{c}{M}_{i\pm 1/2}=\max \left\{\frac{1}{2}(v_i+v_{i\pm 1}),p_{i\pm 1}\left(x_{i\pm 1/2}\right)\right\},\hfill \\ m_{i\pm 1/2}=\min \left\{\frac{1}{2}(v_i+v_{i\pm 1}),p_{i\pm 1}\left(x_{i\pm 1/2}\right)\right\}.\hfill \end{array}\right.$$

Theorem 2.

The polynomial reconstruction (11) is fourth order accurate and fulfills the sign property (8).

Theorem 3.

With the definitions of the reconstructed values$v_{i+1/2}^{+}={\phi }_{i+1}\left(x_{i+1/2}\right)$and$v_{i+1/2}^{-}={\phi }_i\left(x_{i+1/2}\right)$, the numerical flux (6) is entropy stable and fourth order accurate.

Remark 1.

The proofs for Theorems 2 and 3 can be carried out very similarly as in [12,17]. We omit the simple but trivial procedure here.

Remark 2.

For hyperbolic systems, the reconstruction procedure should be implemented in the local characteristic directions for the purpose of achieving entropy stability. The detailed steps can be found in our previous paper [12].

3. Numerical Examples

In this section, we illustrate the effectiveness of the presented scheme which is abbreviated by ES4 by means of three typical examples. Numerical results include the convergence order and the capacity of dealing with discontinuous problems.

Example 1.

Consider the linear advection equation$u_t+u_x=0$on the domain$[-1,1]$with two initial data$u(x,0)=sin\left(\pi x\right)$and$u(x,0)=si{n}^4\left(\pi x\right)$. For a comparison, a third order entropy stable scheme (ES3) [12] is also implemented for this example. The numerical errors ($L^1$error and$L^{\infty }$error) and convergence orders are displayed in

Table 1. The numerical errors and convergence orders for $u_t+u_x=0$ with $u(x,0)=sin\left(\pi x\right)$ at $t=8$.

Method	N	${\mathit{L}}^1$ Error	${\mathit{L}}^1$ Order	${\mathit{L}}^{\mathit{\infty }}$ Error	${\mathit{L}}^{\mathit{\infty }}$ Order
	40	7.8941 $\times {10}^{-3}$	2.9838	6.1563 $\times {10}^{-3}$	2.9528
	80	9.8686 $\times {10}^{-4}$	2.9999	7.7371 $\times {10}^{-4}$	2.9922
ES3	160	1.2332 $\times {10}^{-4}$	3.0004	9.6812 $\times {10}^{-5}$	2.9985
	320	1.5413 $\times {10}^{-5}$	3.0002	1.2104 $\times {10}^{-5}$	2.9997
	640	1.9266 $\times {10}^{-6}$	3.0000	1.5131 $\times {10}^{-6}$	2.9999
	40	7.2962 $\times {10}^{-4}$		5.5350 $\times {10}^{-4}$
	80	4.4746 $\times {10}^{-5}$	4.0273	3.4513 $\times {10}^{-5}$	4.0034
ES4	160	2.7650 $\times {10}^{-6}$	4.0164	2.1523 $\times {10}^{-6}$	4.0032
	320	1.7178 $\times {10}^{-7}$	4.0086	1.3430 $\times {10}^{-7}$	4.0023
	640	1.0702 $\times {10}^{-8}$	4.0046	8.3862 $\times {10}^{-9}$	4.0013

and

Table 2. The numerical errors and convergence orders for $u_t+u_x=0$ with $u(x,0)=si{n}^4\left(\pi x\right)$ at $t=1$.

Method	N	${\mathit{L}}^1$ Error	${\mathit{L}}^1$ Order	${\mathit{L}}^{\mathit{\infty }}$ Error	${\mathit{L}}^{\mathit{\infty }}$ Order
	80	3.8419 $\times {10}^{-3}$		3.7835 $\times {10}^{-3}$
	160	4.9986 $\times {10}^{-4}$	2.9422	5.8676 $\times {10}^{-4}$	2.6889
ES3	320	6.7841 $\times {10}^{-5}$	2.8813	8.6901 $\times {10}^{-5}$	2.7553
	640	1.0357 $\times {10}^{-5}$	2.7115	1.5669 $\times {10}^{-5}$	2.4715
	1280	1.4533 $\times {10}^{-6}$	2.5625	4.9399 $\times {10}^{-6}$	1.6654
	80	8.1572 $\times {10}^{-4}$		1.2001 $\times {10}^{-3}$
	160	5.9041 $\times {10}^{-5}$	3.7883	1.4743 $\times {10}^{-4}$	3.0249
ES4	320	4.9006 $\times {10}^{-6}$	3.5907	2.7023 $\times {10}^{-5}$	1.3202
	640	1.6184 $\times {10}^{-7}$	4.9203	1.4716 $\times {10}^{-7}$	7.5207
	1280	1.0013 $\times {10}^{-8}$	4.0146	8.6289 $\times {10}^{-9}$	4.0921

. We can observe that the fourth order convergence of the proposed scheme is confirmed and ES4 performs better than ES3.

Example 2.

Consider the Burgers equation

$$u_t+{(u^2/2)}_x=0$$

(13)

subjected to the initial data

$$u(x,0)=\left\{\begin{array}{cc}1,\hfill & \mathrm{for}\left|x\right|\le 1/3,\hfill \\ -1,\hfill & \mathrm{for}1/3<\left|x\right|\le 1.\hfill \end{array}\right.$$

For this problem, we can deduce the analytical solution that evolves a rarefaction fan and a stationary shock on the left-hand and right-hand side, respectively. Figure 1 presents the numerical result at time$t=0.3$on a mesh of 100 grids. Our scheme resolves the shock wave and the rarefaction wave very well.

Example 3.

Consider the Euler equations from aerodynamics

$$\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ \rho \mu \\ E\end{array}\right)+\frac{\partial }{\partial x}\left(\begin{array}{c}\rho \mu \\ \rho {\mu }^2+p\\ \mu (E+p)\end{array}\right)=0$$

(14)

with $\rho ,\mu ,p$ and E being the density, velocity, pressure and total energy, respectively. For an idea gas, the total energy E is given by the relation

$$E=\frac{p}{\gamma -1}+\frac{1}{2}\rho {\mu }^2$$

(15)

with the specific heats ratio$\gamma =1.4$. Three Riemann problems are tested by the presented scheme.

Case 1: Sod’s shock tube problem. The initial data is given as

$$(\rho ,\mu ,p)=\left\{\begin{array}{cc}(1,0,1),\hfill & \mathrm{for}x<0,\hfill \\ (0.125,0,0.1),\hfill & \mathrm{for}x>0.\hfill \end{array}\right.$$

The numerical simulation is carried out on a mesh of 200 grids on $[-0.5,0.5]$ up to time $t=0.16$. The computed density is plotted in Figure 2. We can see that the ES4 scheme performs well by capturing the shock, the contact discontinuity and the rarefaction wave accurately.

Case 2: Toro’s 123 problem. The initial data is given as

$$(\rho ,\mu ,p)=\left\{\begin{array}{cc}(1,-2,0.4),\hfill & \mathrm{for}x<0,\hfill \\ (1,2,0.4),\hfill & \mathrm{for}x>0.\hfill \end{array}\right.$$

The difficulty for simulating this problem lies in the fact that the pressure between the evolved rarefactions is very small (near vacuum) and may bring about the blow-ups of the code if the numerical method is not robust. The numerical simulation is carried out on a mesh of 200 grids on $[-0.5,0.5]$ up to time $t=0.1$. Figure 3 displays the computed density. It can be observed that the computed result by ES4 compares well with the reference solution.

Case 3: Shu–Osher problem. The initial data is given as

$$(\rho ,\mu ,p)=\left\{\begin{array}{cc}(3.857,2.629,10.333),\hfill & \mathrm{for}x<-4,\hfill \\ (1+0.2sin(5x),0,1),\hfill & \mathrm{for}x>-4.\hfill \end{array}\right.$$

This problem, also called the shock density-wave interaction problem, describes a moving Mach 3 shock interacting with sine waves in density. The numerical simulation is carried out on a mesh of 500 grids on $[-5,5]$ up to time $t=1.8$. We present the results of density in Figure 4. It can be seen clearly that the ES4 scheme produces accurate results and captures the sine wave well.

4. Conclusions

This paper presents a fourth order entropy stable scheme for solving one-dimensional hyperbolic conservation laws. Along the lines of [10,12], our scheme is also obtained by utilizing entropy conservative flux in conjunction with suitable numerical diffusion. We first select the existing fourth order entropy conservative scheme based on the combination of the two-point entropy conservative flux. The main novelty lies in the construction of numerical diffusion by presenting a fourth order non-oscillatory reconstruction possessing the sign property. Compared to other high order schemes, the main advantage of our scheme is the entropy stability. Some numerical results are displayed to show the accuracy and shock capturing capacity of our scheme. Ongoing work involves generalizing the idea of this paper to multidimensional cases and other hyperbolic systems such as shallow water equations and magnetohydrodynamic equations.