Improved Weighted Essentially Non-Oscillatory Schemes by Mixed Stencils

Shaogang Guo^*, Wang Li

Beijing Institute of Control Engineering, Beijing, China.

DOI: 10.4236/ojfd.2021.113009                

Abstract

Improved Weighted Essentially Non-oscillatory Scheme is a high order finite volume method. The mixed stencils can be obtained by a combination of r + 1 order and r order stencils. We improve the weights by the mapping method. The restriction that conventional ENO or WENO schemes only use r order stencils, is removed. Higher resolution can be achieved by introducing the r + 1 order stencils. This method is verified by three cases, i.e. the interaction of a moving shock with a density wave problem, the interacting blast wave problem and the double mach reflection problem. The numerical results show that the Improved Weighted Essential Non-oscillatory method is a stable, accurate high-resolution finite volume scheme.

Keywords

Finite Volume Method, WENO Schemes, Weights, Shock, Stencils, High-Resolution

1. Introduction

The ENO idea proposed in [1] seems to be the first successful attempt to obtain a self-similar, uniformly high order accurate, yet essentially non-oscillatory interpolation for piecewise smooth functions. The generic solution for hyperbolic conservation laws is in the class of piecewise smooth functions. The reconstruction in [1] was a natural extension of an earlier second order version of Harten and Osher [2]. In [1], Harten, Engquist, Osher and Chakravarthy investigated different ways of measuring local smoothness to determine the local stencil, and developed a hierarchy that begins with one or two cells, then adds one cell at a time to the stencil from the two candidates on the left and right, based on the size of the two relevant Newton divided differences. Although there are other reasonable strategies to choose the stencil based on local smoothness, such as comparing the magnitudes of the highest degree divided differences among all candidate stencils and picking the one with the least absolute value, experience shows that the hierarchy proposed in [1] is the most robust for a wide range of grid sizes, both before and inside the asymptotic regime. Especially, in 1994, Liu, Osher and Chan [3] introduced Weighted ENO schemes which use a convex combination of all candidates to achieve the essentially non-oscillatory property, Weighted ENO improves robustness of the essentially non-oscillatory schemes greatly. In [4], Jiang and Shu proposed a new way of measuring the smoothness of a numerical solution, and used a sum of the L² norms of all the derivatives of the interpolation polynomial over the corresponding interval as a new smoothness measurement, which improve the capacity in resolving shock and complicated flow structures greatly.

Since the WENO schemes were introduced, it has been improved greatly. Recently, Bryson and Doron [5] used a suitable high-order WENO-type reconstruction to develop a high-order semi-discrete central-upwind scheme. Kim, Jang [6] introduced a hybrid methodology to construct a hybrid central-WENO scheme. At the same time, the application of WENO schemes was expanded greatly, Johnsen and Colonius applied WENO schemes in compressible multicomponent flow problems [7], Caleffi, Valiani and Bernin [8] treated with fourth-order balanced source term in central WENO schemes for shallow water equations.

But these schemes only used rth-order stencils ignore larger than rth-order stencils without reference to ENO or WENO. The resolution can improve by using a larger stencil in theory. So it is possible to construct new WENO schemes by introducing larger than rth-order stencils. Despite the (2r − 1)-order convergence behavior often exhibited by WENO schemes, its actual rate of convergence is less than (2r − 1) order for many problems. In fact, the analysis done in [4] does not identify all the properties of the weights which are necessary for WENO schemes to converge at fifth-order. Furthermore, the incomplete set of properties given in [4] is not satisfied by the WENO weights when the first derivative of the solution vanishes [9] [10] [11] [12]. In particular, at critical points where the third derivative does not simultaneously vanish with the first, WENO schemes suffer a loss in accuracy.

In this paper, we present the improved WENO schemes by introducing a large mixed stencil, which is constructed by combining rth-order stencil and (r+ 1)th-order stencil. Simultaneity, we improve the weights by the mapping method [9] [12]. The schemes are tested in three examples, respectively, the interaction of a moving shock with a density wave, the interacting blast waves problem and double mach reflection problem. The numerical results show improved WENO schemes are the high resolution difference methods.

2. WENO Schemes

Firstly, we present the conventional essentially non-oscillatory (ENO) [2] and weighted essentially non-oscillatory (WENO) finite volume methodology [3] [13] [14], for solving nonlinear hyperbolic conservation laws

$u_t+f{\left(u\right)}_x=0$ (2.1)

Given cells $I_j=\left[x_{i-1/2},x_{i+1/2}\right]$ (intervals in one dimension), the partial discretization is written as

$\frac{\text{d}}{\text{d}t}{u}_t+\frac{1}{\Delta x}\left({\hat{f}}_{j+1/2}-{\hat{f}}_{j-1/2}\right)=0$ (2.2)

where ${\hat{f}}_{j+1/2}$ is a numerical approximation of $f\left(u\left(x_{j+1/2}\right)\right)$. Once the right-hand side of this expression has been evaluated, numerical techniques for solving ordinary differential equations, such as Runge–Kutta methods may be employed to advance the solution in time. To ensure stability, $f\left(u\right)$ is generally split into $f^{+}\left(u\right)$, which has a strictly non-negative derivative, and $f^{-}\left(u\right)$, which has a strictly non-positive one [4].

WENO schemes compute ${\hat{f}}_{j+1/2}^{+}$ through interpolating polynomials on a number of overlapping candidate stencils, each containing r grid points. In Jiang and Shu’s method, there are r candidate stencils. The one most upwinded candidate stencil ranges over mesh point indices $j-r+1$ toj, the fully downwinded candidate stencil ranges over j to $j+r-1$, and the other candidate stencils fall in between. Figure 1 provides a schematic of this arrangement.

If the flux approximation on stencil k, which contains r grid points, is designated $q_k^r$ which is the interpolation polynomial on stencil $S_k$ and the weight assigned to that stencil is ${\omega }_k$, the final numerical approximation becomes

${\tilde{f}}_{j+1/2}={\displaystyle \underset{k=0}{\overset{r-1}{\sum }}{\omega }_k{q}_k^r}$ (2.3)

The weights ${\omega }_k$ must adapt to the relative smoothness of each candidate stencil, the smoother the stencil is, the larger the weights are given. A discontinuous stencil is assigned a zero weight (or close to zero), so the (2r − 1)th WENO schemes decrease to rth-order near shocks. At the same time, in [4], optimum weights $C_k^r$ is introduced in the definition of ${\omega }_k$ to make the WENO schemes achieve (2r − 1)th order accuracy in smooth domain. The optimum weights $C_k^r$ is calculated as

$q_{r-1}^{2r-1}={\displaystyle \underset{k=0}{\overset{r-1}{\sum }}{C}_k^r{q}_k^r}$ (2.4)

the left term is a (2r − 1)th-order approximation of $h_{j+1/2}$.

Figure 1. The assembly method of the 3th-order stencils in Jiang and Shu’s method.

${\omega }_k=\frac{{\alpha }_k}{{\alpha }_0+\cdots +{\alpha }_{r-1}}$, (2.5)

where

${\alpha }_k=\frac{C_k^r}{{\left(\epsilon +I{S}_k\right)}^p},k=0,1,\cdots ,r-1$ (2.6)

The constant $\epsilon$ prevents division by zero, and $p=2$ is chosen to increase or decrease WENO adaptation sensitivity. The smoothness measurement $I{S}_k$ becomes large when discontinuities are present within stencil k and remains relatively small otherwise. The strict definition of $I{S}_k$ as presented by Jiang and Shu [4] may be expressed in the form.

$I{S}_k={\displaystyle \underset{l=1}{\overset{k-1}{\sum }}{\displaystyle \underset{x_{j-1/2}}{\overset{x_{j+1/2}}{\int }}\Delta x^{2l-1}{\left(\frac{{\partial }^l{q}_k^r\left(x\right)}{\partial x^l}\right)}^2\text{d}x}}$ (2.7)

The corresponding stencil diagram for ${\tilde{f}}_{j+1/2}^{-}$ is simply a mirror image of Figure 1 [10]. It is not overstated here.

3. New WENO Scheme

But if we take too high order stencils, it is probable that some grid points of different stencils overlap, so it is impossible to obtain high resolution. Of course, we can reduce the order of the stencils according to [15] [16], but a considerable amount of calculation will be needed. What’s more, this iterative calculation is often unsuccessful, because the iteration always begins with the (2r − 1)th-order stencil regardless of the structure of stencils (even obvious repetitive stencils). As for different equations or different grid points, it is usually impossible that reasonable amount of stencils than rth-order stencils are selected. Now, we only select (r+ 1)th-order to combine with rth-order to construct a new WENO schemes.

Now, we introduce a large mixed stencil T which is combined by using all the rth and (r + 1)th order stencils:

$\begin{array}{l}{S}_k^r=\left(I_{j+k-r+1},I_{j+k-r+2},\cdots ,I_{j+k}\right),k=0,1,\cdots ,r-1\\ S_k^{r+1}=\left(I_{j+k-r+1+1},I_{j+k-r+1+2},\cdots ,I_{j+k}\right),k=1,2,\cdots ,r-1\end{array}$ (3.1)

$T={\displaystyle {\cup }_{k=0}^{k=r-1}{S}_k^r}\cup {\displaystyle {\cup }_{k=1}^{k=r-1}{S}_k^{r+1}}$ (3.2)

The assembly method of the rth and (r+ 1)th order stencils (take r= 3 for example) is shown in Figure 1 and Figure 2.

The rth-order stencils: (Figure 2).

The (r+ 1)th-order stencils: (Figure 3).

We can construct (2r − 2)th order polynomial $P\left(x\right)$ in T, Based on the idea of the polynomial construction in [10]. Simple algebra gives the coefficients $C_k^m$ such that

$P\left(x_{j+1/2}\right)={\displaystyle \underset{k=0}{\overset{2r-1}{\sum }}{C}_k^m{q}_k^m\left(x_{j+1/2}\right)}$ (3.3)

Figure 2. The assembly method of the rth-order stencils.

Figure 3. The assembly method of the (r + 1)th-order stencils.

here, $q_k^m\left(x\right)$ is the corresponding polynomial of the kth small stencil, for example, when the stencil is (r + 1)th-order, $q_k^m\left(x\right)$ is the (r + 1)th-order polynomial. Obviously, the methodology to obtain the coefficients $C_k^m$ is similar to Weighted ENO schemes in [4].

When r= 3, there are the results of calculation in [17], obviously, there are five (2r − 1) small stencil $S_k^m\left(m=r,r+1\right)$, and then we rearrange these small stencils and assemble a large stencil,

$\begin{array}{l}{S}_0=S_1^4=\left(I_{j-2},I_{j-1},I_j,I_{j+1}\right),S_1=S_0^3=\left(I_{j-2},I_{j-1},I_j\right),\\ S_2=S_1^3=\left(I_{j-1},I_j,I_{j+1}\right),S_3=S_2^3=\left(I_j,I_{j+1},I_{j+2}\right),\\ S_4=S_2^4=\left(I_{j-1},I_{j-1},I_j,I_{j+2}\right),T={\displaystyle {\cap }_{k=0}^4{S}_k}\end{array}$ (3.4)

The procedure of calculating numerical fluxes is similar to the WENO scheme proposed by Jiang and Shu [4], and the interpolation to obtain the cell-face fluxes ${\hat{f}}_{j+1/2}^{\pm }$ is described in [4].

For brevity, we define a general flux ${\hat{f}}_{j+1/2}$ that denote the positive flux $f^{+}$ or the negative flux $f^{-}$, which is split by the Rusanov-type flux splitting method [18] [19] [20] as

$\left(\begin{array}{l}{g}_{j-2}=f_{j-2}^{+}\\ g_{j-21}=f_{j-1}^{+}\\ g_j=f_j^{+}\\ g_{j+1}=f_{j+1}^{+}\\ g_{j+2}=f_{j+1}^{+}\end{array}\right),\left(\begin{array}{l}{g}_{j-1}=f_{j+3}^{-}\\ g_j=f_{j+2}^{-}\\ g_{j+1}=f_{j+1}^{-}\\ g_{j+2}=f_j^{-}\\ g_{j+3}=f_{j-1}^{-}\end{array}\right)$ (3.5)

${\hat{f}}_{j+1/2}={\hat{f}}_{j+1/2}^{+}+{\hat{f}}_{j+1/2}^{-}$ (3.6)

We take ${\hat{f}}_{j+1/2}^{+}$ for example, ${\hat{f}}_{j+1/2}^{+}$ is calculated as

$\begin{array}{l}{\hat{f}}_{j+1/2}^{+}={\omega }_0\left(\frac{1}{12}{g}_{j-2}-\frac{5}{12}{g}_{j-1}+\frac{13}{12}{g}_j+\frac{1}{4}{g}_{j+1}\right)+{\omega }_1\left(\frac{1}{3}{g}_{j-2}-\frac{7}{6}{g}_{j-1}+\frac{11}{6}{g}_j\right)\\ +{\omega }_2\left(-\frac{1}{6}{g}_{j-2}+\frac{5}{6}{g}_{j-1}+\frac{1}{3}{g}_j\right)+{\omega }_3\left(\frac{1}{3}{g}_{j-2}+\frac{5}{6}{g}_{j-1}-\frac{1}{6}{g}_j\right)\\ +{\omega }_4\left(\frac{1}{4}{g}_{j-1}+\frac{13}{12}{g}_j-\frac{5}{12}{g}_{j+1}+\frac{1}{12}{g}_{j+2}\right)\end{array}$ (3.7)

In order to obtain the smallest truncation error, we write a Taylor series expansion for $g_{(•)}$ around $x_{i+1/2}$ :

$g_j=g_{j+1/2}-{g^{\prime }}_{j+1/2}\left(h/2\right)+{g^{″}}_{j+1/2}{\left(h/2\right)}^2/2!+\cdots +{\left(-1\right)}^n{g}_{j+1/2}^{\left(n\right)}{\left(h/2\right)}^n/n!+\cdots$ (3.8)

So, we can write the Taylor series expansion as the sum of all the term in (3.9):

$\begin{array}{l}\frac{1}{12}{g}_{j-2}-\frac{5}{12}{g}_{j-1}+\frac{13}{12}{g}_j+\frac{1}{4}{g}_{j+1}\\ =g_{j+1/2}-0.0416667h^2{g^{″}}_{j+1/2}+0.0512153h^4{g}_{j+1/2}^{\left(4\right)}-0.0416667h^5{g}_{j+1/2}^{\left(5\right)}\\ +0.0216942h^6{g}_{j+1/2}^{\left(6\right)}+\cdots +\frac{1}{3}{g}_{j-2}-\frac{7}{6}{g}_{j-1}+\frac{11}{6}{g}_j\end{array}$

$\begin{array}{l}=g_{j+1/2}-0.0416667h^2{g^{″}}_{j+1/2}-0.25h^3{g^{‴}}_{j+1/2}+0.301215h^4{g}_{j+1/2}^{\left(4\right)}\\ -0.197917h^5{g}_{j+1/2}^{\left(5\right)}+0.0946108h^6{g}_{j+1/2}^{\left(6\right)}+\cdots -\frac{1}{6}{g}_{j-2}+\frac{5}{6}{g}_{j-1}+\frac{1}{3}{g}_j\\ =g_{j+1/2}-0.0416667h^2{g^{″}}_{j+1/2}+0.0833333h^3{g^{‴}}_{j+1/2}-0.0321181h^4{g}_{j+1/2}^{\left(4\right)}\\ +0.0104167h^5{g}_{j+1/2}^{\left(5\right)}-0.0026114h^6{g}_{j+1/2}^{\left(6\right)}+\cdots +\frac{1}{3}{g}_{j-2}+\frac{5}{6}{g}_{j-1}-\frac{1}{6}{g}_j\end{array}$

$\begin{array}{l}=g_{j+1/2}-0.0416667h^2{g^{″}}_{j+1/2}-0.0833333h^3{g^{‴}}_{j+1/2}-0.0321181h^4{g}_{j+1/2}^{\left(4\right)}\\ -0.0104167h^5{g}_{j+1/2}^{\left(5\right)}-0.0026114h^6{g}_{j+1/2}^{\left(6\right)}+\cdots \\ +\frac{1}{4}{g}_{j-1}+\frac{13}{12}{g}_j-\frac{5}{12}{g}_{j+1}+\frac{1}{12}{g}_{j+2}\\ =g_{j+1/2}-0.0416667h^2{g^{″}}_{j+1/2}-0.0321181h^4{g}_{j+1/2}^{\left(4\right)}-0.0026114h^6{g}_{j+1/2}^{\left(6\right)}+\cdots \end{array}$ (3.9)

Firstly, Supposed that the coefficients of ${\hat{f}}_{j+1/2}^{+}$ is A,

$A=\left[\begin{array}{ccccc}\frac{1}{12}& \frac{-5}{12}& \frac{13}{12}& \frac{1}{4}& 0\\ \frac{1}{3}& \frac{-7}{6}& \frac{11}{6}& 0& 0\\ 0& \frac{-1}{6}& \frac{5}{6}& \frac{1}{3}& 0\\ 0& 0& \frac{1}{3}& \frac{5}{6}& \frac{-1}{6}\\ 0& \frac{-1}{12}& \frac{7}{12}& \frac{7}{12}& \frac{-1}{12}\end{array}\right],Rank\left(A\right)=Rank\left(\left[\begin{array}{ccc}\frac{1}{3}& \frac{-7}{6}& \frac{11}{6}\\ \frac{-1}{6}& \frac{5}{6}& \frac{1}{3}\\ \frac{1}{3}& \frac{5}{6}& \frac{-1}{6}\end{array}\right]\right)$ (3.10)

Here, the rank of matrix A is 3, so it is so impossible to get the solution of $C_j$ like WENO5 [4]. But we can obtain the optimal solution by complementary condition. The primary objective is to make the truncation error the least, so we make the term including ${g^{″}}_{j+1/2}$ of (3.9) equal to zero [17].

Now, we can calculate the coefficients $C_k^m$, and the results is written as

$C_0=C_1^4=0.2,C_1=C_0^3=0.05,C_2=C_1^3=0.3,C_3=C_2^3=0.15,C_4=C_2^4=0.3$ (3.11)

Because of ${\omega }_i=C_i+O\left(h\right)$, so the optimal coefficients can be tested by replacing ${\omega }_j$ by $C_j$.

$\begin{array}{l}{\hat{f}}_{j+1/2}^{+}\approx 0.2\left(\frac{1}{12}{g}_{j-2}-\frac{5}{12}{g}_{j-1}+\frac{13}{12}{g}_j+\frac{1}{4}{g}_{+1}\right)+0.05\left(\frac{1}{3}{g}_{j-2}-\frac{7}{6}{g}_{j-1}+\frac{11}{6}{g}_j\right)\\ +0.3\left(-\frac{1}{6}{g}_{j-2}+\frac{5}{6}{g}_{j-1}+\frac{1}{3}{g}_j\right)+0.15\left(\frac{1}{3}{g}_{j-2}+\frac{5}{6}{g}_{j-1}-\frac{1}{6}{g}_j\right)\\ +0.3\left(\frac{1}{4}{g}_{j-1}+\frac{13}{12}{g}_j-\frac{5}{12}{g}_{j+1}+\frac{1}{12}{g}_{j+2}\right)\end{array}$ (3.12)

Corresponding to $C_k^m$, the smoothness coefficients $I{S}_k^m$ is calculated by (2.7):

$I{S}_0=I{S}_1^4,I{S}_1=I{S}_0^3,I{S}_2=I{S}_1^3,I{S}_3=I{S}_2^3,I{S}_4=I{S}_2^4$ (3.13)

Here,

$\begin{array}{l}I{S}_0=(g_{j-2}\left(267g_{j-2}-1642g_{j-1}+1602g_j-494g_{j+1}\right)\\ +g_{j-1}\left(2843g_{j-1}-5966g_j+1922g_{j+1}\right)\\ +g_j\left(3443g_j-2522g_{j+1}\right)+547g_{j+1}^2)/240\\ I{S}_1=\frac{13}{12}{\left(g_j-2g_{j-1}+g_{j-2}\right)}^2+\frac{1}{4}{\left(3g_j-4g_{j-1}+g_{j-2}\right)}^2\\ I{S}_2=\frac{13}{12}{\left(g_{j+1}-2g_j+g_{j-1}\right)}^2+\frac{1}{4}{\left(g_{j+1}-g_{j-1}\right)}^2\end{array}$

$\begin{array}{l}I{S}_3=\frac{13}{12}{\left(g_{j+2}-2g_{j+1}+g_j\right)}^2+\frac{1}{4}{\left(3g_j-4g_{j+1}+g_{j+2}\right)}^2\\ I{S}_4=(g_{j-1}\left(547g_{j-1}-2522g_j+1922g_{j+1}-494g_{j+2}\right)\\ +g_j\left(3443g_j-5966g_{j+1}+1602g_{j+2}\right)\\ +g_{j+1}\left(2843g_{j+1}-1642g_{j+2}\right)+267g_{j+2}^3)/240\end{array}$ (3.14)

Here, the weights ${\omega }_k^m$ :

${\omega }_k={\alpha }_k/{\displaystyle \underset{k=0}{\overset{4}{\sum }}{\alpha }_k},k=0,1,2,3,4$ (3.15)

${\alpha }_k=C_k/\left(\epsilon +I{S}_k\right),k=0,1,2,3,4$ (3.16)

To increase the accuracy of these weights, consider the functions [9].

$g_k\left(\omega \right)=\frac{\omega \left({\bar{\omega }}_k+{\bar{\omega }}_k^2-3{\bar{\omega }}_k\omega +{\omega }^2\right)}{{\bar{\omega }}_k^2+\omega \left(1-2{\bar{\omega }}_k\right)},{\bar{\omega }}_k\in \left(0,1\right)$ (3.17)

for $k=0,1,2,3,4$. All of these functions have the following features: they are monotonically increasing with finite slope, $g_k\left(0\right)=0$, $g_k\left(1\right)=1$, $g_k\left({\bar{\omega }}_k\right)={\bar{\omega }}_k$, ${g^{\prime }}_k\left({\bar{\omega }}_k\right)=0$, ${g^{″}}_k\left({\bar{\omega }}_k\right)=0$.

A more accurate approximation of the weights is given by

${\alpha }_k^{\ast }=g_k\left(w_k\right)$ (3.18)

The modified weights are defined according to

$w_k^{\ast }=\frac{{\alpha }_k^{*}}{{\displaystyle \underset{i=0}{\overset{4}{\sum }}{\alpha }_i^{*}}}$ (3.19)

which satisfies ${\displaystyle \underset{i=0}{\overset{4}{\sum }}{w}_k^{\ast }}=1$, Here a superscript ( $\ast$ ) has been added to signify the improved weights, which is called MWENO.

4. Numerical Results

Calculations are carried out for three different examples of compressible flows. Example 1 is the interaction of a moving shock with a density wave, example 2 is interacting blast waves, example 3 is double mach reflection problem.

Example 1: Interaction of a moving shock with a density wave.

We applied the developed improved Weighted essentially non-oscillatory scheme to the Shu and Osher’s 1D shock/turbulence interaction model problem [21]. This problem describes the interaction between a moving Mach 3 shock wave and a fluctuating density since wave [22]. The initial flow condition is given by:

${\left(\rho ,u,p\right)}^{\text{T}}=\{\begin{array}{l}\left(1+0.2\sin \left(5x\right),0,1\right)x\ge -4.0\\ {\left(3.857143,2.629369,10.33333\right)}^{\text{T}}x<-4.0\end{array}$.

The computational domain, $-5\le x\le 15$, is covered with 400 nodes. We obtain a reference solution with 3200 grid points by WENO5 schemes. The Roe scheme is used for calculating the fluxes in the system of hyperbolic conservation laws and the third-order Runge–Kutta scheme is used to advance in time, and the CFL number is 0.2. The solutions at t = 4.5 are given in Figure 4 and Figure 5 for the WENO5 and the present WENO5 schemes, respectively. It is clear that the present scheme resolves the salient features of the flow with higher fidelity.

Figure 4. Density distributions of the 1D shock/turbulence interaction at t = 4.5 (WENO5).

Figure 5. Density distributions of the 1D shock/turbulence interaction at t = 4.5 (Present).

Example 2: Interacting blast waves.

The interacting blast wave example [23] is another difficult test of shock capturing schemes. The Euler equation is solved with initial conditions with reflection boundary conditions

${\left(\rho ,u,p\right)}^{\text{T}}=\{\begin{array}{l}\left(1,0,1000\right)0\le x<0.1\\ \left(1,0,0.01\right)0.1\le x<0.9\\ {\left(1,0,100\right)}^{\text{T}}0.9\le x\le 1\end{array}$, at $x=0$ and $x=1$.

Again, the Roe scheme and third-order Runge–Kutta are too used in constructing the numerical solution. The solutions at t = 0.038 are given in Figure 6 and Figure 7 for the WENO5 and the present schemes, respectively. Each computational cells is N = 400. Careful examination reveals that the present scheme has better agreement with the exact solution in resolving the salient features of the flow. The three contact waves, near x = 0.595, 0.765, 0.799, are better resolved as are the various peaks and valleys.

Example 3: Double mach reflection problem of a shock wave.

This problem was initially proposed and studied in detail by Woodward and Colella [23]. It has been used extensively in the literature as a test case for high resolution schemes. In this case, a strong shock wave impinges on a ramp with a 30˚ angle. The inlet flow is strong shock wave whose mach number is 10. In order to simplify the problem, the inlet shock wave is inclined with a 60˚ angle with a reflecting wall, and the computational domain is beginning with x = 1/6. The computational domain for this problem is chosen to be (0, 4) × (0, 1). For the top boundary of the computational domain, the solution is set to describe the exact motion of the Mach 10 shock. The left boundary is set as the exact post-shock condition, while the right boundary is set as outflow boundary. The initial condition

Figure 6. Interacting blast waves. Density plot of WENO5 solution at t = 0.038 with N = 400, and WENO5 solution with N = 6400, solid line.

Figure 7. Interacting blast waves. Density plot of the present solution at t = 0.038 with N = 400, and WENO5 solution with N = 6400, solid line.

is set: $u\left(x,y,0\right)=\{\begin{array}{l}{u}_L,y\ge h\left(x,0\right)\\ u_R,y<h\left(x,0\right)\end{array}$, here, left and right state values and shock wave’s height:

$u_L=\left(8,57.1597,-33.0012,563.544\right)$, $u_R=\left(1.4,0,0,2.5\right)$, $h\left(x,t\right)=\sqrt{3}\left(x-\frac{1}{6}\right)-20t$,

and computational grid number is 960 × 240. All data are obtained at t = 0.2 and the CFL number is set to 0.5.

Figure 8 shows density contours using the WENO5 scheme, which shows the result of using the WENO5 scheme with a 960 × 240 uniform grid system. Figure 9 shows the results of using the present scheme with the same uniform grid systems. The complicated flow features in the ‘‘blow-up’’ region around the double Mach stems are resolved better by the present scheme. The present scheme can capture the vortical flow structures along the slip line very well. The close-up view around the ‘‘blow-up’’ region is shown in Figure 10(a) and Figure 10(b). It is seen that the present scheme gives greatly improved resolution compared with the WENO5 scheme.

Figure 8. Double mach reflection problem (WENO5), t = 0.2, density contours (960 × 240).

Figure 9. Double mach reflection problem (the present scheme),t = 0.2, density contours (960 × 240).

(a) (b)

Figure 10. Close-up view of Figure 8 and Figure 9, WENO5 (Left), the present scheme (right).

5. Conclusion

A new WENO method is developed by combined (r + 1)th-order stencil and rth-order stencil to construct multiple stencils. The method is tested by three customary cases: interaction of a moving shock with a density wave, the interacting blast waves problem, and the double mach reflection problem. The numerical results have shown the method is the stable, accurate difference scheme. Furthermore, in the high-velocity area and complicated boundary layer, the scheme can distinguish and capture discontinuous physical phenomena like shock waves, and can obtain higher resolution and restrain numerical oscillation effectively. Especially, the improved WENO schemes make full use of the selected stencils over the general WENO schemes in the smooth domain, and the schemes improve the capacity of capturing discontinuities by introducing the higher order stencils.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Harten, A. (1986) On High-Order Accurate Interpolation for Non-Oscillatory Shock Capturing Schemes. In: Dafermos, C., et al., Eds., Oscillation Theory. Computation and Methods of Compensated Compactness, Springer-Verlag, New York, 71-105. https://doi.org/10.1007/978-1-4613-8689-6_4
[2]	Harten, A. and Osher, S. (1987) Uniformly High Order Accurate Essential Non-Oscillatory Schemes, I. SIAM Journal on Numerical Analysis, 24, 279-309. https://doi.org/10.1137/0724022
[3]	Liu, X., Osher, S. and Chan, T. (1994) Weighted Essentially Non-Oscillatory Schemes. Journal of Computational Physics, 115, 200-212. https://doi.org/10.1006/jcph.1994.1187
[4]	Jiang, G.-S. and Shu, C.-W. (1996) Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics, 126, 202-228. https://doi.org/10.1006/jcph.1996.0130
[5]	Bryson, S., and Levy, D. (2003) High-Order Semi-Discrete Central-Upwind Schemes for Multi-Dimensional Hamilton-Jacobi Equations. Journal of Computational Physics, 189, 63-87. https://doi.org/10.1016/S0021-9991(03)00201-8
[6]	Kim, D. and Kwon, J.H. (2005) A High-Order Accurate Hybrid Scheme Using a Central Flux Scheme and a WENO Scheme for Compressible Flowfield Analysis. Journal of Computational Physics, 210, 554-583. https://doi.org/10.1016/j.jcp.2005.04.023
[7]	Johnsen, E. and Colonius, T. (2006) Implementation of WENO Schemes in Compressible Multicomponent Flow Problems. Journal of Computational Physics, 219, 715-732. https://doi.org/10.1016/j.jcp.2006.04.018
[8]	Caleffi, V., Valiani, A. and Bernin, A. (2006) Fourth-Order Balanced Source Term Treatment in Central WENO Schemes for Shallow Water Equations. Journal of Computational Physics, 218, 228-245. https://doi.org/10.1016/j.jcp.2006.02.001
[9]	Henrick, A.K., Aslam, T.D. and Powers, J.M. (2005) Mapped Weighted Essentially Non-Oscillatory Schemes: Achieving Optimal Order near Critical Points. Journal of Computational Physics, 207, 542-567. https://doi.org/10.1016/j.jcp.2005.01.023
[10]	Martin, M.P., Taylor, E.M., Wu, M. and Weirs, V.G. (2006) A Bandwidth-Optimized WENO Scheme for the Effective Direct Numerical Simulation of Compressible Turbulence. Journal of Computational Physics, 220, 270-289. https://doi.org/10.1016/j.jcp.2006.05.009
[11]	Zheng, F., Shu, C.W. and Qiu, J. (2019) High Order Finite Difference Hermite WENO Schemes for the Hamilton-Jacobi Equations on Unstructured Meshes. Computers & Fluids, 183, 53-65. https://doi.org/10.1016/j.compfluid.2019.02.010
[12]	Wu, L., Zhang, Y.T., Zhang, S., et al. (2016) High Order Fixed-Point Sweeping WENO Methods for Steady State of Hyperbolic Conservation Laws and Its Convergence Study. Communications in Computational Physics, 20, 835-869. https://doi.org/10.4208/cicp.130715.010216a
[13]	Shu, C.-W. (1997) Essential Non-Oscillatory and Weighted Essential Non-Oscillatory Schemes for Hyperbolic Conservation Laws. NASA CR-97-206253. https://doi.org/10.1007/BFb0096355
[14]	Shu, C.-W. (1997) Preface to the Republication of “Uniformly High Order Essentially Non-Oscillatory Schemes, III,” by Harten, Engquist, Osher, and Chakravarthy. Journal of Computational Physics, 131, 1-2. https://doi.org/10.1006/jcph.1996.5630
[15]	Zong, W.G., Deng, X.G. and Zhang, H.X. (2003) Double Weighted Essential Non-Oscillatory Shock-Capturing Schemes. Acta Aerodynamica Sinica, 21, 218-224.
[16]	Zong, W.G., Deng, X.G. and Zhang, H.X. (2003) Implementation of DWENO Schemes and Its Application. Acta Aerodynamica Sinica, 21, 399-407.
[17]	Jun, Z. and Ning, Z. (2005) A Kind of MWENO Schemes and Its Application. Acta Aerodynamica Sinica, 23, 330-336.
[18]	Rusanov, V.V. (1961) Calculation of Interaction of Non-Steady Shock Waves with Obstacles. National Research Council of Canada, Ottawa.
[19]	Levy, D., Puppo, G. and Russo, G. (1999) Central WENO Schemes for Hyperbolic Systems of Conservation Laws. Mathematical Modelling and Numerical Analysis, 33, 547-571. https://doi.org/10.1051/m2an:1999152
[20]	Levy, D., Puppo, G. and Russo, G. (2002) A Fourth-Order Central WENO Scheme for Multidimensional Hyperbolic Systems of Conservation Laws. SIAM Journal on Scientific Computing, 24, 480-506. https://doi.org/10.1137/S1064827501385852
[21]	Shu, C. (1998) Essentialy Non-Oscilatory and Weighted Essential Non-Oscilatory Schemes for Hyperbolic Conservation Laws. In: Ockburn, B., Jo Hnsonc, S., Hu, C. and Tadmor, E., Eds., Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Lecture Notes in Mathematics, Vol. 1697, Springer-Verlag, Berlin, 325-432. https://doi.org/10.1007/BFb0096355
[22]	Shu, C.W. and Osher, S. (1998) Efficient Implementation of Essentially Non-Oscillatory Shock Capturing Schemes II. Journal of Computational Physics, 83, 32-78. https://doi.org/10.1016/0021-9991(89)90222-2
[23]	Woodward, P. and Colella, P. (1984) The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks. Journal of Computational Physics, 54, 115-173. https://doi.org/10.1016/0021-9991(84)90142-6