High-Order Accurate Flux-Splitting Scheme for Conservation Laws with Discontinuous Flux Function in Space

Abstract

A new modified Engquist–Osher-type flux-splitting scheme is proposed to approximate the scalar conservation laws with discontinuous flux function in space. The fact that the discontinuity of the fluxes in space results in the jump of the unknown function may be the reason why it is difficult to design a high-order scheme to solve this hyperbolic conservation law. In order to implement the WENO flux reconstruction, we apply the new modified Engquist–Osher-type flux to compensate for the discontinuity of fluxes in space. Together the third-order TVD Runge–Kutta time discretization, we can obtain the high-order accurate scheme, which keeps equilibrium state across the discontinuity in space, to solve the scalar conservation laws with discontinuous flux function. Some examples are given to demonstrate the good performance of the new high-order accurate scheme.

1. Introduction

We focus on the high-order accurate schemes for the conservation laws with discontinuous flux in the form

$$u_t+{\left(F(x,u)\right)}_x=0,$$

(1)

where $u=u(x,t)$ is an unknown function, $(x,t)\in {\mathcal{D}}_{\mathcal{T}}:=\mathbb{R}\times (0,T)$, $F(x,u)=(1-H(x\left)\right)g\left(u\right)+H\left(x\right)f\left(u\right)$ and $H\left(x\right)$ is the Heaviside function. Because of its interesting features and applications in traffic flows, multiple-phase flow in porus media and clarifier-thickener units model and so on, equation of this form (1) has been widely investigated in recent years [1,2,3,4,5,6,7,8].

The discontinuous flux $F(x,u)$ in x results in that the Kružkov-type entropy inequalities [9] do not ensure the uniqueness of solutions to the Cauchy problem of equation of the form (1). In order to overcome this difficulty, different entropy theories have been imposed to obtain distinct choice solutions and the reader is referred to [4,5,6,7,10,11,12,13,14] for details.

Because of the discontinuity of the flux in the spatial variable x, it is very difficult to develop an efficient numerical scheme with general numerical fluxes. This implies that the interface numerical flux plays an important role in designing numerical scheme for conservation laws with discontinuous flux (1). In recent years, some efficient first-order schemes have been proposed to prove the uniqueness of weak solutions to the related model such as conservation laws (1). Based on the Godunov or EO flux, Towers [13,14] studied conservation laws with a discontinuous spatially varying coefficient or spatially varying source terms and established the convergence. Adimurthi et al. [15,16] developed the Godunov-type or DFLU numerical flux to solve scalar conservation laws with discontinuous flux. Then, Sudarshan et al. [17] extended the DFLU interface flux to two-phase multicomponent polymer flooding model and obtained a second-order scheme through the MUSCL-type reconstruction. In addition, Mishra [18] considered mixed fluxes of the convex-concave-type and concave-convex-type in (1) and defined the corresponding interface Godunov-type flux and Engquist–Osher-type flux. Karlsen and Towers [19] used the Lax–Friedrichs scheme to study conservation laws with a discontinuous space-time dependent flux. In particular, Adimurthi et al. [20] claimed that equation of the form (1) had a contractive semigroup solution in $L^1$ for each so-called connection-$(A,B)$. Then, Bürger, Karlsen, Towers [21] put forward the entropy solutions of type $(A,B)$ and utilized the modified EO-type numerical scheme to prove the convergence of this type solution. Recently, Adimurthi et al. [10,22] proposed some numerical flux, including Godunov, Engquist–Osher and modified Lax–Friedrichs-type numerical fluxes to define the interface flux for $(A,B)$-connection. In the aforementioned work, the numerical simulations show that the DFLU scheme [15,16,23] and Engquist–Osher-type scheme [21] work better than others for many models. We also see the modified Local Lax–Friedrichs flux (MLLF) [10,22], Roe-type flux [24] and application of the DFLU flux [23,25] and the references therein.

One purpose of this paper is to propose a simple and an efficient interface numerical flux, which can be easily extended to high-order schemes to approximate the conservation laws with discontinuous flux (1). With the help of the entropy solutions of type $(A,B)$, we propose a new Engquist–Osher-type interface numerical flux, which can be written as the split form.

For the high-order accurate schemes for conservation laws with discontinuous flux (1), it is worth mentioning that Shu and his collaborators [26,27] used the Lax–Friedrichs numerical flux and the WENO reconstruction to approximate the multiclass traffic flow model on homogeneous or inhomogeneous highway. For the details of the WENO reconstruction method, we refer to [28,29,30]. Recently, Bürger, Karlsen and their collaborators [1,2,31] developed a series of second-order accurate schemes to solve the conservation laws with discontinuous flux (1). Then, borrowing the idea of flux-total variation diminishing to conservation laws with discontinuous flux introduced by Towers in [14], Bürger et al. [2] developed a nonlocal limiter algorithm and proposed a second-order flux-total variation diminishing scheme to solve the clarifier-thickener model with discontinuous flux. Adimurthi et al. made a modification of nonlocal limiter algorithm and obtained the second-order flux-TVD scheme for some interface fluxes as defined in [10]. We also refer the readers to the high-order Runge–Kutta discontinuous Galerkin scheme [32] and the semidiscrete central-upwind scheme [33].

However, so far, there are few results on the high-order schemes for (1). To the best of our knowledge, it is not easy to construct high-order schemes for (1) by using the reconstruction of the primitive variables. The reason for this may be that the discontinuity of the flux in the spatial variable x leads to the jump of the unknown function. Fortunately, the flux $F(x,u)=(1-H(x\left)\right)g\left(u\right)+H\left(x\right)f\left(u\right)$ in (1) should satisfy the Rankine–Hugoniot condition

$$g\left(u_{-}\right)=f\left(u_{+}\right)$$

(2)

at the discontinuity $x=0$, where $u_{\pm }$ represent the limits on both sides of the discontinuity $x=0$, respectively. Based on this condition, we could apply the proposed Engquis–Osher type interface flux to connect two fluxes $g\left(u\right)$ and $f\left(u\right)$ and keep the equilibrium to capture the steady-state solution. Then, we can make use of the WENO flux reconstruction method, which introduced in [28,29,30], to develop an efficient, nonoscillatory and high-order accurate scheme for (1). This is another purpose of this paper.

The format of the paper is as follows. In Section 2, we give a new modified Engquist–Osher-type interface numerical flux, which is written into the flux-splitting form. In Section 3, we propose three high-order schemes for (1) by using the modified Engquist–Osher type interface flux or DFLU flux and the WENO reconstruction. In Section 4, some numerical examples are presented. Finally, we make a discussion about three high-order accurate schemes.

2. The Modified Engquist–Osher-Type Interface Numerical Flux

Based on the $(A,B)$-connection [20] and the entropy solutions of type $(A,B)$ [21], we suppose that the flux $F(x,u)=(1-H(x\left)\right)g\left(u\right)+H\left(x\right)f\left(u\right)$ in (1) satisfies:

(A1)$f,g\in Lip\left(\right[0,1\left]\right)$, $f\left(u\right)=g\left(u\right)=0$ for $u=0,1$, f has a single maximum at ${\theta }_f\in [0,1]$ and g has a single maximum at ${\theta }_g\in [0,1]$. The function f (or g) is strictly increasing on $(0,{\theta }_f)$ (or $(0,{\theta }_g)$) and strictly decreasing on $({\theta }_f,1)$ (or $({\theta }_g,1)$).

(A2) both f and g are genuinely nonlinear in the sense that f and g are not linear on any nondegenerate interval.

(A3) there exists at most one point $\widehat{u}\in (0,1)$ such that $f\left(\widehat{u}\right)=g\left(\widehat{u}\right)$.

For convenience, set

$$h\left(u\right)=\left\{\begin{array}{c}g\left(u\right),x<0,\hfill \\ f\left(u\right),x>0.\hfill \end{array}\right.$$

The corresponding Engquist–Osher numerical flux is

$$\widehat{h}(u_{-},u_{+})=\frac{1}{2}(h\left(u_{-}\right)+h\left(u_{+}\right))-\frac{1}{2}{\int }_{u_{-}}^{u_{+}}\left|h^{\prime }\left(\omega \right)\right|d\omega ,$$

(3)

which can be written as the flux-splitting form

$$\widehat{h}(u_{-},u_{+})=h^{+}\left(u_{-}\right)+h^{-}\left(u_{+}\right),$$

(4)

where

$$h^{+}\left(u_{-}\right)=h\left(0\right)+{\int }_0^{u_{-}}max\{h^{\prime }\left(\omega \right),0\}d\omega ,h^{-}\left(u_{+}\right)={\int }_0^{u_{+}}min\{h^{\prime }\left(\omega \right),0\}d\omega .$$

(5)

As stated above, the interface numerical flux plays a key role in designing the numerical methods for (1). So, we modify the Engquist–Osher-type interface numerical flux to approximate (1) with $(A,B)$-connection.

Under the assumptions $\left(A_1\right)$–$\left(A_3\right)$, without the loss of generality, the graphs of functions $g\left(u\right)$, $f\left(u\right)$ and the connection $(A,B)$ are shown in Figure 1. Obviously, there exist $A_g,B_f\in (0,1)$ such that they satisfy $A_g\le {\theta }_g$, $B_f\ge {\theta }_f$ and $g\left(A_g\right)=g\left(A\right)=f\left(B\right)=f\left(B_f\right)$. With the following definitions

$${\tilde{g}}^{\prime }\left(u\right)=\left\{\begin{array}{ll}{g}^{\prime }\left(u\right),\hfill & \mathrm{if}g\left(u\right)<g\left(A\right),\\ 0,\hfill & \mathrm{if}g\left(u\right)\ge g\left(A\right),\end{array}\right.{\tilde{f}}^{\prime }\left(u\right)=\left\{\begin{array}{ll}{f}^{\prime }\left(u\right),\hfill & \mathrm{if}f\left(u\right)<f\left(B\right),\\ 0,\hfill & \mathrm{if}f\left(u\right)\ge f\left(B\right),\end{array}\right.$$

(6)

We can generalize the Engquist–Osher interface numerical flux (4) and (5) in the form

$$\tilde{h}(u_{-},u_{+})={\tilde{g}}^{+}\left(u_{-}\right)+{\tilde{f}}^{-}\left(u_{+}\right)$$

(7)

to capture the entropy solutions of type $(A,B)$, where

$${\tilde{g}}^{+}\left(u_{-}\right)=g\left(0\right)+{\int }_0^{u_{-}}max\{{\tilde{g}}^{\prime }\left(\omega \right),0\}d\omega ,{\tilde{f}}^{-}\left(u_{+}\right)={\int }_0^{u_{+}}min\{{\tilde{f}}^{\prime }\left(\omega \right),0\}d\omega .$$

(8)

In this paper, we call formulas (7) and (8) the modified Engquist–Osher interface flux.

By the discussion case by case, we can obtain that the modified Engquist–Osher interface numerical flux (7) and (8) is equivalent to the interface flux introduced by Bürger, Karlsen, Towers in [21]. Besides this, the new interface flux-splitting form (7) and (8) is simple and can be extended to the high-order approximation by using the WENO reconstruction, see Section 3.

Next, we consider the important properties of the flux $\tilde{h}(x,y)$.

Lemma 1.

Based on the flux function $\tilde{h}(x,y)$, we have:

(1) 

$\tilde{h}(x,y)$ is nondecreasing function in its first variable, nonincreasing function in its second variable, and Lipschitz continuous on both two arguments.

(2) 

$\tilde{h}(x,y)$ is not consistent, but $\tilde{h}(0,0)=0$, $\tilde{h}(1,1)=0$, and $\tilde{h}(A,B)=g\left(A\right)=f\left(B\right)$.

Proof. 

We first demonstrate the statement (1). Based on (6)–(8), we have

$${\displaystyle \frac{\partial \tilde{h}(x,y)}{\partial x}=max\{{\tilde{g}}^{\prime }\left(x\right),0\}\ge 0,\mathrm{and}\frac{\partial \tilde{h}(x,y)}{\partial y}=min\{{\tilde{f}}^{\prime }\left(y\right),0\}\le 0.}$$

For any $(x_1,y),(x_2,y)$, $(x,y_1),(x,y_2)\in [0,1]$, we have

$$\begin{array}{ccc}\hfill \left|\tilde{h}(x_2,y)-\tilde{h}(x_1,y)\right|& =& \left|{\int }_{x_1}^{x_2}max\{{\tilde{g}}^{\prime }\left(\omega \right),0\}d\omega \right|\le \mathcal{M}\left|x_2-x_1\right|,\hfill \\ \hfill \left|\tilde{h}(x,y_2)-\tilde{h}(x,y_1)\right|& =& \left|{\int }_{y_1}^{y_2}min\{{\tilde{f}}^{\prime }\left(\omega \right),0\}d\omega \right|\le \mathcal{M}\left|y_2-y_1\right|,\hfill \end{array}$$

where

$$\mathcal{M}=\left\{|g^{\prime }\left(u\right)|,|f^{\prime }\left(u\right)|\right\},u\in [0,1].$$

(9)

Thus, the statement (1) is true.

The proof of assertion (2) is trivial. So, we proof the Lemma. □

Now, we consider the initial value problem (1) with the initial data

$$u_0\left(x\right)=\left\{\begin{array}{c}{u}_l,x<0,\hfill \\ u_r,x>0,\hfill \end{array}\right.$$

(10)

where $u_l$ and $u_r$ are some constants. For simplicity, we use a uniform space and time grid with space size $\Delta x>0$ and time step $\Delta t>0$ to approximate solution of (1). Define the space grid points $x_{j+1/2}=(j+1/2)\Delta x$ and $x_{-1/2}=x_{1/2}$. Let $I_j=[x_{j-1/2},x_{j+1/2}]$, $\Delta x=x_{j+1/2}-x_{j-1/2}$ for $j\in \mathcal{Z},j\ne 0$. Set $\lambda =\Delta t/\Delta x$ and $t_n=n\Delta t$ for $n=0,1,2,\cdots$. Based on (4)–(9), we construct the numerical scheme for (1) as follows

$$U_j^{n+1}=U_j^n-\lambda \left({\widehat{F}}_{j+\frac{1}{2}}^n-{\widehat{F}}_{j-\frac{1}{2}}^n\right),j\ne 0,$$

(11)

where

$$\begin{array}{c}\hfill {\widehat{F}}_{j+\frac{1}{2}}^n=\widehat{F}(u_j^n,u_{j+1}^n)=\left\{\begin{array}{c}\widehat{g}(u_j^n,u_{j+1}^n),forj\le -2,\hfill \\ \tilde{h}(u_{-1}^n,u_1^n),forj=-1,0,\hfill \\ \widehat{f}(u_j^n,u_{j+1}^n),forj\ge 1\hfill \end{array}\right.\end{array}$$

(12)

In this paper, the first-order accurate scheme (11) and (12) with the CFL condition

$$\lambda \mathcal{M}\le \frac{1}{2}$$

(13)

is referred to as the modified Engquist–Osher-type numerical scheme, denoted by MEO, where the mesh ratio $\lambda =\Delta t/\Delta x$ and the Lipschitz constant $\mathcal{M}$ is defined in (9).

3. High-Order Accurate Schemes for Conservation Laws (1)

In this section, we present two types of high-order accurate schemes to solve (1) by using the idea of WENO reconstruction.

3.1. High-Order WENO Scheme Based on the Modified Interface Engquist–Osher-Type Flux

Although the unknown function u is jump across the discontinuity $x=0$ for (1), the Rankine–Hugoniot condition (2) implies that we could add the interface flux (7) to connect the two fluxes. So, we can take advantage of the idea of WENO reconstruction of fluxes to construct high-order scheme, which is of the form

$${\displaystyle \frac{d}{dt}{u}_j\left(t\right)=\mathcal{L}\left(u_j\right)=-\frac{1}{\Delta x}\left({\widehat{h}}_{j+\frac{1}{2}}-{\widehat{h}}_{j-\frac{1}{2}}\right),}$$

(14)

in which ${\widehat{h}}_{j+\frac{1}{2}}$ is the numerical flux to be given below.

To ensure the numerical stability and correct upwind biasing, we utilize the Engquist–Osher-type flux-splitting form (4) and (7) to construct the numerical flux as follows

$${\widehat{h}}_{j+\frac{1}{2}}={\widehat{h}}_{j+\frac{1}{2}}^{+}+{\widehat{h}}_{j+\frac{1}{2}}^{-}.$$

(15)

Here, we can adopt the fifth-order WENO reconstruction method [28,29,30] to construct ${\widehat{h}}_{j+1/2}^{+}$ and ${\widehat{h}}_{j+1/2}^{-}$. The procedure includes the following two steps.

Step 1. Away from the discontinuity $x=0$.

We give three third-order accurate numerical fluxes for the case $h^{\prime }\left(u\right)\ge 0$ as follows

$$\begin{array}{ccc}\hfill {\widehat{h}}_{j+\frac{1}{2}}^{\left(1\right)}& =& \frac{1}{3}{h}_j^{+}+\frac{5}{6}{h}_{j+1}^{+}-\frac{1}{6}{h}_{j+2}^{+},\hfill \\ \hfill {\widehat{h}}_{j+\frac{1}{2}}^{\left(2\right)}& =& -\frac{1}{6}{h}_{j-1}^{+}+\frac{5}{6}{h}_j^{+}+\frac{1}{3}{h}_{j+1}^{+},\hfill \\ \hfill {\widehat{h}}_{j+\frac{1}{2}}^{\left(3\right)}& =& \frac{1}{3}{h}_{j-2}^{+}-\frac{7}{6}{h}_{j-1}^{+}+\frac{11}{6}{h}_j^{+}.\hfill \end{array}$$

(16)

In this case, we describe the smoothness indicators as

$$\begin{array}{ccc}\hfill {\beta }_1& =& \frac{13}{12}{(h_j^{+}-2h_{j+1}^{+}+h_{j+2}^{+})}^2+\frac{1}{4}{(3h_j^{+}-4h_{j+1}^{+}+h_{j+2}^{+})}^2,\hfill \\ \hfill {\beta }_2& =& \frac{13}{12}{(h_{j-1}^{+}-2h_j^{+}+h_{j+1}^{+})}^2+\frac{1}{4}{(h_{j-1}^{+}-h_{j+1}^{+})}^2,\hfill \\ \hfill {\beta }_3& =& \frac{13}{12}{(h_{j-2}^{+}-2h_{j-1}^{+}+h_j^{+})}^2+\frac{1}{4}{(h_{j-2}^{+}-4h_{j-1}^{+}+3h_j^{+})}^2\hfill \end{array}$$

(17)

and give the linear weights

$$d_1=\frac{3}{10},d_2=\frac{3}{5},d_3=\frac{1}{10}.$$

Then, we use the nonlinear weights

$${\omega }_k=\frac{{\alpha }_k}{{\sum }_{s=1}^3{\alpha }_s}\ge 0,{\alpha }_k=\frac{d_k}{{(\epsilon +{\beta }_k)}^2},k=1,2,3$$

to define

$${\widehat{h}}_{j+\frac{1}{2}}^{+}=\sum _{r=1}^3{\omega }_r{\widehat{h}}_{j+\frac{1}{2}}^{\left(r\right)}.$$

(18)

Here, $\epsilon$ is a small constant and taken as ${10}^{-6}$ to prevent the denominator to be zero. Performing the mirror symmetric with respect to $j+1/2$, we can obtain the reconstruction procedure for the flux ${\widehat{h}}_{j+1/2}^{-}$. We omit the details and refer to [28,29,30] on this issue.

Step 2. Near the discontinuity $x=0$.

Since the flux changes when passing the discontinuity $x=0$, it is not easy for us to reconstruct the ${\widehat{g}}_{-j-1/2}^{\pm }$ and ${\widehat{f}}_{j+1/2}^{\pm }$ for $j=0,1,2$. In order to overcome this difficulty, we could replace $g_{-1}^{+}$ and $f_1^{-}$ with this ${\tilde{g}}_{-1}^{+}$ and ${\tilde{f}}_1^{-}$ respectively, as shown in Figure 2. So we set

$${\widehat{v}}_j=\left\{\begin{array}{c}{g}_j^{+},j\le -2,\hfill \\ {\tilde{g}}_{-1}^{+},j=-1,\hfill \\ f_{j+1}^{+},j\ge 0,\hfill \end{array}\right.$$

(19)

to construct ${\widehat{h}}_{j-1/2}^{+}={\widehat{v}}_{j-1/2}^{+}$ for $j\le 0$ and set

$${\widehat{v}}_j=\left\{\begin{array}{c}{g}_j^{-},j\le -1,\hfill \\ {\tilde{f}}_1^{-},j=0,\hfill \\ f_{j+1}^{-},j\ge 1,\hfill \end{array}\right.$$

(20)

to construct ${\widehat{h}}_{j-1/2}^{-}={\widehat{v}}_{j-1/2}^{-}$ for $j\le -1$ as the procedure in Section 3.1. Similarly, we also reconstruct ${\widehat{h}}_{j+1/2}^{+}$ and ${\widehat{h}}_{j-1/2}^{-}$ for $j\ge 1$.

Now, we apply the third-order TVD Runge–Kutta time discretization introduced by Shu et al. in [34] to (1) and give

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle u^{\left(1\right)}=u^n+\Delta t\mathcal{L}\left(u^n\right),}\\ {\displaystyle u^{\left(2\right)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{\left(1\right)}+\frac{1}{4}\Delta t\mathcal{L}\left(u^{\left(1\right)}\right),}\\ {\displaystyle u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{\left(2\right)}+\frac{2}{3}\Delta t\mathcal{L}\left(u^{\left(2\right)}\right),}\end{array}\right.\end{array}$$

(21)

in which the operator $\mathcal{L}$ is defined in (14). Combining (14) and (21) gives the high-order accurate scheme for (1), which is termed as MEO-WENO5.

3.2. High-Order WENO Schemes Based on the Reconstructions of the Unknown Functions

In this subsection, we adopt the so-called DFLU numerical flux introduced in [15,16] and the WENO reconstruction method to present the other two high-order WENO schemes.

For the fluxes $g\left(u\right)$, $f\left(u\right)$ and the connection $(A,B)$ shown in Figure 1, the DFLU numerical fluxes can be defined as the following form

$$\begin{array}{ccc}\hfill \widehat{h}(x,y)& =& min\left\{h(min(x,{\theta }_h)),h(max(y,{\theta }_h))\right\},h=g\mathrm{or}f,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \tilde{h}(x,y)& =& min\left\{g(min(x,A_g)),f(max(y,B_f))\right\}.\hfill \end{array}$$

(23)

Then, the numerical flux ${\widehat{h}}_{j+1/2}$ in (14) can be written as

$${\widehat{h}}_{j+\frac{1}{2}}=\left\{\begin{array}{c}\widehat{h}(u_{j+\frac{1}{2}}^{n-},u_{j+\frac{1}{2}}^{n+}),\mathrm{for}j\ne -1,0,\hfill \\ {\tilde{h}}_{-\frac{1}{2}}={\tilde{h}}_{\frac{1}{2}}=\tilde{h}(u_{-\frac{1}{2}}^{n-},u_{\frac{1}{2}}^{n+}),\mathrm{for}j=-1,0,\hfill \end{array}\right.$$

(24)

where $u_{j+1/2}^{n\pm }$ are prescribed as follows.

We directly utilize the unknown functions $u_j^n$ ($j\ne 0$) to construct $u_{j+1/2}^{n-}$ and $u_{j-1/2}^{n+}$ for $j\ne 0$. Similar to the procedure in Section 3.1, we set

$${\widehat{v}}_i=\left\{\begin{array}{c}{u}_i^n,i\le -1,\hfill \\ u_{i+1}^n,i\ge 0.\hfill \end{array}\right.$$

(25)

and use five variables ${\widehat{v}}_i$, $j-2\le i\le j+2,i\in {\mathcal{Z}}^{+}$ to reconstruct $u_{j+1/2}^{n-}$ and $u_{j-1/2}^{n+}$ for $j\le -1$ and at the same time set

$${\tilde{v}}_i=\left\{\begin{array}{c}{u}_i^n,i\ge 1,\hfill \\ u_{i-1}^n,i\le 0.\hfill \end{array}\right.$$

(26)

to reconstruct $u_{j+1/2}^{n-}$ and $u_{j-1/2}^{n+}$ for $j\ge 1$. Thus combining (14), (24) and (21) gives the high-order accurate WENO scheme for (1), denoted by DFLU-WENO5.

It is important to capture steady-state solution and prevent spurious oscillations near the discontinuity $x=0$ for (1). Then, we adopt the idea of the well-balanced method [35,36,37] to construct $u_{j+1/2}^{\pm }$. Suppose $(u_A,u_B)$ be a pair of well-balanced solution or steady-state solution to the initial data problem (1) and (10). Different from the above DFLU-WENO5 method, we use

$${\widehat{v}}_i=\left\{\begin{array}{c}{u}_i^n,i\le -1,\hfill \\ u_A,i\ge 0,\hfill \end{array}\right.$$

(27)

satisfying $j-2\le i\le j+2$ and $i\in {\mathcal{Z}}^{+}$, to construct $u_{j+1/2}^{n-}$ and $u_{j-1/2}^{n+}$ for $j\le -1$ and

$${\tilde{v}}_i=\left\{\begin{array}{c}{u}_i^n,i\ge 1,\hfill \\ u_B,i\le 0,\hfill \end{array}\right.$$

(28)

satisfying $j-2\le i\le j+2,i\in {\mathcal{Z}}^{+}$ to construct $u_{j+1/2}^{n-}$ and $u_{j-1/2}^{n+}$ for $j\ge 1$. Once again, we implement (21) to obtain the high-order scheme with the well-balanced method or preserving the steady-state solution, which is termed as DFLU-WENO5B.

4. Numerical Examples

To demonstrate the efficiency of our methods, we here present some numerical results of the new first-order MEO scheme, DFLU scheme consisting of numerical fluxes (22) and (23), and three high-order WENO schemes which are described in Section 2 and Section 3 for comparison. For convenience, the first-order MEO scheme and the first-order DFLU scheme are termed as MEO-1 and DFLU-1, respectively.

Example 1.

We consider the case for traffic flow, where the fluxes can be taken as$g\left(u\right)=u(1-u)$and$f\left(u\right)=1.5u(1-u)$in (1). We take the initial data as follows

$$u_0\left(x\right)=\left\{\begin{array}{c}0.8,x<0,\hfill \\ 0.6,x>0.\hfill \end{array}\right.$$

(29)

For this case, the connect-$(A,B)$can be taken as$(0.5,0.2113)$and the steady-state solution is$(u_A,u_B)=(A,B)$. The numerical parameter$\lambda =\Delta t/\Delta x=0.25$and the compute time$t=1.0$. In this case, the exact Riemann solution can be given in follows

$$u(x,t)=\left\{\begin{array}{c}0.8,x\le -0.6t,\hfill \\ {\displaystyle 0.5(1-\frac{x}{t}),-0.6t<x<0,}\hfill \\ {\displaystyle B,0<x<\frac{f\left(0.6\right)-f\left(B\right)}{0.6-B}t,}\hfill \\ {\displaystyle 0.6,x>\frac{f\left(0.6\right)-f\left(B\right)}{0.6-B}t.}\hfill \end{array}\right.$$

(30)

The $L^1$-errors computed by the above five methods MEO-1, DFLU-1, MEO-WENO5, DFLU-WENO5, and DFLU-WENO5B are given in the

Table 1. $L^1$-errors of two types of MEO and DFLU schemes for Example 1.

$\mathbf{\Delta }\mathit{x}$	MEO-1		DFLU-1		MEO-WENO5		DFLU-WENO5		DFLU-WENO5B
	${\mathit{L}}^1$-Error	Order	${\mathit{L}}^1$-Error	Order	${\mathit{L}}^1$-Error	Order	${\mathit{L}}^1$-error	Order	${\mathit{L}}^1$-Error	Order
1/25	1.88 × 10${}^{-2}$		1.80 × 10${}^{-2}$		8.36 × 10${}^{-3}$		8.70 × 10${}^{-3}$		8.28 × 10${}^{-3}$
1/50	1.14 × 10${}^{-2}$	0.7230	1.11 × 10${}^{-2}$	0.7025	4.57 × 10${}^{-3}$	0.8727	4.57 × 10${}^{-3}$	0.9278	4.57 × 10${}^{-3}$	0.8567
1/100	6.99 × 10${}^{-3}$	0.7049	6.88 × 10${}^{-3}$	0.6846	2.76 × 10${}^{-3}$	0.7270	2.63 × 10${}^{-3}$	0.7983	2.62 × 10${}^{-3}$	0.8023
1/200	4.48 × 10${}^{-3}$	0.6431	4.46 × 10${}^{-3}$	0.6274	1.93 × 10${}^{-3}$	0.5144	1.96 × 10${}^{-3}$	0.4260	1.89 × 10${}^{-3}$	0.4670
1/400	2.29 × 10${}^{-3}$	0.9681	2.25 × 10${}^{-3}$	0.9840	6.13 × 10${}^{-4}$	1.6555	6.11 × 10${}^{-4}$	1.6794	6.11 × 10${}^{-4}$	1.6329

. Although $L^1$-errors of the MEO-1 scheme are slightly larger than those of the DFLU-1 scheme, they have the same numerical results. However, MEO-WENO5 scheme has smaller $L^1$-errors than DFLU-WENO5 starting with $\Delta x=1/25$ and $\Delta x=1/50$, which indicates that the MEO-WENO5 has better performance than the DFLU-WENO5. In addition, we can observe that the DFLU-WENO5B scheme has smaller numerical $L^1$-errors than the MEO-WENO5 and DFLU-WENO5 schemes.

To further analyze these schemes, we present the numerical solutions with space sizes $\Delta x=1/25$ and $\Delta x=1/50$, where the solid line represents exact solutions in the following figures. Figure 3 shows that the MEO-1 and the DFLU-1 perform equally well and they can capture the solution to initial value problem (1) and (29). Figure 4 illustrates that the MEO-WENO5 works better than the DFLU-WENO5 near the discontinuity $x=0$. Moreover, we observe that the numerical solution computed by MEO-WENO5 is in better agreement with the exact solution than the DFLU-WENO5. The main reason for this situation is that the DFLU-WENO5 scheme produces spurious oscillations behind the shock and the discrepancy between the numerical solution and the exact solution can be observed near the discontinuity $x=0$. As the mesh is refined, the above oscillations disappears but the discrepancy still exists, as shown in Figure 5.

We now solve the initial value problem (1) and (29) with space sizes $\Delta x=1/25$ and $\Delta x=1/50$ by MEO-WENO5 and DFLU-WENO5B, respectively. The numerical solutions are shown in Figure 6, from which a good agreement with the exact solution presented by both two schemes is observed. Moreover, the DFLU-WENO5B does not introduce the oscillation behind the shock, but the discrepancy mentioned above on the left of the discontinuity still exists.

Example 2.

To investigate the above two aspects, we take the fluxes$g=u{(1-u)}^2$and$f=u^2(1-u)$coming from [21] and the initial conditions $u_0\left(x\right)=1.0$ for $x<0$ and $u_0\left(x\right)=0.0$ for $x>0$. Set $\lambda =\Delta t/\Delta x=0.15$ and the connection-$(A,B)=(1/3,2/3)$. So we obtain the steady-state solution $(u_A,u_B)=(1/3,2/3)$. For the fluxes $g\left(u\right)$ and $f\left(u\right)$ in this case, there exist $u_1\in [A,1]$ and $u_2\in [0,B]$ such that they satisfy

$${\displaystyle g^{\prime }\left(u_1\right)=\frac{g\left(1\right)-g\left(u_1\right)}{1-u_1}\mathrm{and}{\mathrm{f}}^{\prime }\left({\mathrm{u}}_2\right)=\frac{\mathrm{f}\left({\mathrm{u}}_2\right)-\mathrm{f}\left(0\right)}{{\mathrm{u}}_2-0},}$$

may see Figure 1. Now, we present the exact Riemann solution as

$$u(x,t)=\left\{\begin{array}{c}1,x\le g^{\prime }\left(u_1\right)t,\hfill \\ {\displaystyle {\left(g^{\prime }\right)}^{-1}\left(\frac{x}{t}\right),g^{\prime }\left(u_1\right)t\le x<0,}\hfill \\ {\displaystyle {\left(f^{\prime }\right)}^{-1}\left(\frac{x}{t}\right),0<x\le f^{\prime }\left(u_2\right)t,}\hfill \\ 0,x\ge f^{\prime }\left(u_2\right)t,\hfill \end{array}\right.$$

(31)

which includes a contact discontinuity connecting two states $u_1$ and 1 and another contact discontinuity connecting two states 0 and $u_2$.

We compute the solution up to $t=1.0$ and present the numerical $L^1$-errors of the above five schemes in

Table 2. $L^1$-errors and orders of two types of MEO and DFLU schemes for Example 2.

$\mathbf{\Delta }\mathit{x}$	MEO-1		DFLU-1		MEO-WENO5		DFLU-WENO5		DFLU-WENO5B
	${\mathit{L}}_1$-Error	Order	${\mathit{L}}_1$-Error	Order	${\mathit{L}}_1$-Error	Order	${\mathit{L}}_1$-Error	Order	${\mathit{L}}_1$-Error	Order
1/25	6.52 × 10${}^{-2}$		6.68 × 10${}^{-2}$		2.69 × 10${}^{-2}$		2.49 × 10${}^{-2}$		2.18 × 10${}^{-2}$
1/50	4.55 × 10${}^{-2}$	0.5203	4.66 × 10${}^{-2}$	0.5222	1.82 × 10${}^{-2}$	0.5617	1.73 × 10${}^{-2}$	0.5092	1.59 × 10${}^{-2}$	0.4505
1/100	3.17 × 10${}^{-2}$	0.5225	3.23 × 10${}^{-2}$	0.5253	1.36 × 10${}^{-2}$	0.4280	1.31 × 10${}^{-2}$	0.3933	1.24 × 10${}^{-2}$	0.3626
1/200	1.91 × 10${}^{-2}$	0.7295	1.95 × 10${}^{-2}$	0.7269	7.02 × 10${}^{-3}$	0.9496	6.84 × 10${}^{-3}$	0.9427	6.45 × 10${}^{-3}$	0.9428
1/400	1.13 × 10${}^{-2}$	0.7586	1.16 × 10${}^{-2}$	0.7562	3.58 × 10${}^{-3}$	0.9733	3.48 × 10${}^{-3}$	0.9742	3.29 × 10${}^{-3}$	0.9743

, which shows that the $L^1$-error of MEO-1 is slight smaller than that of DFLU-1, meanwhile both DFLU-WENO5 and DFLU-WENO5B schemes have smaller numerical error than the MEO-WENO5 scheme. Next, we further analyze the numerical solution to explore this situation.

We first compute the initial value problem by the MEO-1 and DFLU-1 schemes with space sizes $\Delta x=1/50$ and $\Delta x=1/100$, respectively. It is easy to see that they have similar behavior, see Figure 7.

Now we present the numerical solutions by the three high-order accurate schemes with space size $\Delta x=1/25$, as illustrated in Figure 8. The solution of two rarefactive waves on both sides of the discontinuity $x=0$ has been well approximated by the MEO-WENO5 scheme. However, the numerical solutions computed by the DFLU-WENO5 and DFLU-WENO5B are slightly different from the exact solution. Especially, the DFLU-WENO5 scheme gives rise to more oscillations.

Finally, the numerical solutions are presented in Figure 9 and Figure 10 as mesh refinement with space size $\Delta x=1/50$. Figure 9 shows that the DFLU-WENO5 introduces the spurious oscillations on both sides of $x=0$, but Figure 10 illustrates that the oscillation disappears for the DFLU-WENO5B. In addition, both the DFLU-WENO5 and DFLU-WENO5B perform better resolution of contact discontinuity away from the discontinuity $x=0$ than the MEO-WENO5. This may be the reason why the $L^1$-errors performed by the DFLU-WENO5 and DFLU-WENO5B are smaller than those performed by the MEO-WENO5.

Example 3

([15]). Considering the two-phase incompressible flow in a porous medium with a rock type changing at $x=0$. The flux functions $g\left(u\right)$ and $f\left(u\right)$ have the following form

$$\begin{array}{ccc}& & {\displaystyle g\left(u\right)=\frac{1}{\varphi }\frac{{\mu }_1}{{\mu }_1+{\mu }_2}[q+(c_1-c_2){\mu }_2]\mathrm{for}x<0,}\hfill \\ & & {\displaystyle f\left(u\right)=\frac{1}{\varphi }\frac{{\lambda }_1}{{\lambda }_1+{\lambda }_2}[q+(c_1-c_2){\lambda }_2]\mathrm{for}x>0.}\hfill \end{array}$$

Under the following data

$$\begin{array}{ccc}& & \varphi =1,q=0,c_1=2,c_2=1,\hfill \\ & & {\mu }_1=50u^2,{\mu }_2=5{(1-u)}^2,{\lambda }_1=10u^2,{\lambda }_2=20{(1-u)}^2.\hfill \end{array}$$

Two fluxes $g\left(u\right)$ and $f\left(u\right)$ are represented in Figure 1. The initial data is $u_0\left(x\right)=1.0$ for $x<0$ and $u_0\left(x\right)=0.0$ for $x>0$ with the connection $(A,B)=(0.317014,0.472372)$. The difference from the previous Example 2 is that the contact discontinuity on the left of the stationary discontinuity $x=0$ is replaced by the constant state B and a composite wave consisting of a rarefaction wave and a contact discontinuity. We take the numerical parameter $\lambda =\Delta t/\Delta x=0.1$ and compute the solution up to time $t=0.5$.

We first make use of the MEO-1 and DFLU-1 schemes to compute this problem with $\Delta x=1/50$ and $\Delta x=1/100$, respectively. The numerical solutions are illustrated in Figure 11. It is apparent that the MEO-1 scheme performs as well as the DFLU-1 scheme.

Now, we utilize three high-order schemes to compute this problem with space size $\Delta x=1/25$ and present the numerical solutions. In Figure 12, we see that the DFLU-WENO5 scheme gives rise to spurious oscillations and has slight deviation from the exact solution on the left of $x=0$. We easily observe from Figure 13 that both MEO-WENO5 and DFLU-WENO5B schemes are in good agreement with the exact solution. However, the discrepancy between the numerical solution computed by the DFLU-WENO5B scheme and the exact solution can still be found in the rarefaction wave region on the left of $x=0$.

5. Conclusions

In this paper, we propose the new first-order modified Engquist–Osher-type scheme (MEO), which performs as well as the DFLU scheme [15,16,25] and the modified EO-scheme [21]. The new modified EO-type interface flux could be written as the split form (7) and (8), which contributes to stability of the scheme and could be taken as a building block for constructing high-order scheme. As mentioned before, the flux $F(x,u)=(1-H(x\left)\right)g\left(u\right)+H\left(x\right)f\left(u\right)$ in (1) should satisfy the Rankine–Hugoniot condition (2), thus the interface numerical flux (7) could be used as a bridge connecting two flux functions $g\left(u\right)$ and $f\left(u\right)$. Therefore, we bring this result to the WENO reconstruction of fluxes and achieve the high-order accurate scheme (MEO-WENO5). Based on the well-balanced method or the steady-state solution, we apply the WENO method to construct another high-order accurate scheme (DFLU-WENO5B). In comparison with the DFLU-WENO5 and DFLU-WENO5B, the MEO-WENO5 makes use of the reconstruction of fluxes, instead of the reconstruction of unknown functions, which shows its conservative to the discontinuity and still keeps nonoscillatory, especially in areas near the discontinuity $x=0$. However, the scheme proposed in this manuscript could not be directly generalized to approximate the system of hyperbolic balance laws involving the discontinuous flux in space, since it is difficult to define the approximate interface flux. However, by virtue of some special methods, such as the idea of discontinuous flux, we could solve some special system of hyperbolic conservation laws involving the discontinuous flux.

In short, the modified Engquist–Osher-type interface flux (7) and (8) and the high-order accurate MEO-WENO5 scheme have the advantages of being simple, efficient, and nonoscillatory. In the future work, we extend this method to solve the corresponding system of hyperbolic equations.