^{1}

^{*}

^{1}

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.

The ENO idea proposed in [^{2} 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 [

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 [

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 [

Firstly, we present the conventional essentially non-oscillatory (ENO) [

u t + f ( u ) x = 0 (2.1)

Given cells I j = [ x i − 1 / 2 , x i + 1 / 2 ] (intervals in one dimension), the partial discretization is written as

d d t u t + 1 Δ x ( f ^ j + 1 / 2 − f ^ j − 1 / 2 ) = 0 (2.2)

where f ^ j + 1 / 2 is a numerical approximation of f ( u ( x j + 1 / 2 ) ) . 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 ( u ) is generally split into f + ( u ) , which has a strictly non-negative derivative, and f − ( u ) , which has a strictly non-positive one [

WENO schemes compute 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.

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 ω k , the final numerical approximation becomes

f ˜ j + 1 / 2 = ∑ k = 0 r − 1 ω k q k r (2.3)

The weights ω 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 [

q r − 1 2 r − 1 = ∑ k = 0 r − 1 C k r q k r (2.4)

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

ω k = α k α 0 + ⋯ + α r − 1 , (2.5)

where

α k = C k r ( ε + I S k ) p , k = 0 , 1 , ⋯ , r − 1 (2.6)

The constant ε 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 [

I S k = ∑ l = 1 k − 1 ∫ x j − 1 / 2 x j + 1 / 2 Δ x 2 l − 1 ( ∂ l q k r ( x ) ∂ x l ) 2 d x (2.7)

The corresponding stencil diagram for f ˜ j + 1 / 2 − is simply a mirror image of

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 [

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

S k r = ( I j + k − r + 1 , I j + k − r + 2 , ⋯ , I j + k ) , k = 0 , 1 , ⋯ , r − 1 S k r + 1 = ( I j + k − r + 1 + 1 , I j + k − r + 1 + 2 , ⋯ , I j + k ) , k = 1 , 2 , ⋯ , r − 1 (3.1)

T = ∪ k = 0 k = r − 1 S k r ∪ ∪ 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

The rth-order stencils: (

The (r+ 1)th-order stencils: (

We can construct (2r − 2)th order polynomial P ( x ) in T, Based on the idea of the polynomial construction in [

P ( x j + 1 / 2 ) = ∑ k = 0 2 r − 1 C k m q k m ( x j + 1 / 2 ) (3.3)

here, q k m ( x ) is the corresponding polynomial of the kth small stencil, for example, when the stencil is (r + 1)th-order, q k m ( x ) is the (r + 1)th-order polynomial. Obviously, the methodology to obtain the coefficients C k m is similar to Weighted ENO schemes in [

When r= 3, there are the results of calculation in [

S 0 = S 1 4 = ( I j − 2 , I j − 1 , I j , I j + 1 ) , S 1 = S 0 3 = ( I j − 2 , I j − 1 , I j ) , S 2 = S 1 3 = ( I j − 1 , I j , I j + 1 ) , S 3 = S 2 3 = ( I j , I j + 1 , I j + 2 ) , S 4 = S 2 4 = ( I j − 1 , I j − 1 , I j , I j + 2 ) , T = ∩ k = 0 4 S k (3.4)

The procedure of calculating numerical fluxes is similar to the WENO scheme proposed by Jiang and Shu [

For brevity, we define a general flux 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 [

( 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 + ) , ( 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 − ) (3.5)

f ^ j + 1 / 2 = f ^ j + 1 / 2 + + f ^ j + 1 / 2 − (3.6)

We take f ^ j + 1 / 2 + for example, f ^ j + 1 / 2 + is calculated as

f ^ j + 1 / 2 + = ω 0 ( 1 12 g j − 2 − 5 12 g j − 1 + 13 12 g j + 1 4 g j + 1 ) + ω 1 ( 1 3 g j − 2 − 7 6 g j − 1 + 11 6 g j ) + ω 2 ( − 1 6 g j − 2 + 5 6 g j − 1 + 1 3 g j ) + ω 3 ( 1 3 g j − 2 + 5 6 g j − 1 − 1 6 g j ) + ω 4 ( 1 4 g j − 1 + 13 12 g j − 5 12 g j + 1 + 1 12 g j + 2 ) (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 ′ j + 1 / 2 ( h / 2 ) + g ″ j + 1 / 2 ( h / 2 ) 2 / 2 ! + ⋯ + ( − 1 ) n g j + 1 / 2 ( n ) ( h / 2 ) n / n ! + ⋯ (3.8)

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

1 12 g j − 2 − 5 12 g j − 1 + 13 12 g j + 1 4 g j + 1 = g j + 1 / 2 − 0.0416667 h 2 g ″ j + 1 / 2 + 0.0512153 h 4 g j + 1 / 2 ( 4 ) − 0.0416667 h 5 g j + 1 / 2 ( 5 ) + 0.0216942 h 6 g j + 1 / 2 ( 6 ) + ⋯ + 1 3 g j − 2 − 7 6 g j − 1 + 11 6 g j

= g j + 1 / 2 − 0.0416667 h 2 g ″ j + 1 / 2 − 0.25 h 3 g ‴ j + 1 / 2 + 0.301215 h 4 g j + 1 / 2 ( 4 ) − 0.197917 h 5 g j + 1 / 2 ( 5 ) + 0.0946108 h 6 g j + 1 / 2 ( 6 ) + ⋯ − 1 6 g j − 2 + 5 6 g j − 1 + 1 3 g j = g j + 1 / 2 − 0.0416667 h 2 g ″ j + 1 / 2 + 0.0833333 h 3 g ‴ j + 1 / 2 − 0.0321181 h 4 g j + 1 / 2 ( 4 ) + 0.0104167 h 5 g j + 1 / 2 ( 5 ) − 0.0026114 h 6 g j + 1 / 2 ( 6 ) + ⋯ + 1 3 g j − 2 + 5 6 g j − 1 − 1 6 g j

= g j + 1 / 2 − 0.0416667 h 2 g ″ j + 1 / 2 − 0.0833333 h 3 g ‴ j + 1 / 2 − 0.0321181 h 4 g j + 1 / 2 ( 4 ) − 0.0104167 h 5 g j + 1 / 2 ( 5 ) − 0.0026114 h 6 g j + 1 / 2 ( 6 ) + ⋯ + 1 4 g j − 1 + 13 12 g j − 5 12 g j + 1 + 1 12 g j + 2 = g j + 1 / 2 − 0.0416667 h 2 g ″ j + 1 / 2 − 0.0321181 h 4 g j + 1 / 2 ( 4 ) − 0.0026114 h 6 g j + 1 / 2 ( 6 ) + ⋯ (3.9)

Firstly, Supposed that the coefficients of f ^ j + 1 / 2 + is A,

A = [ 1 12 − 5 12 13 12 1 4 0 1 3 − 7 6 11 6 0 0 0 − 1 6 5 6 1 3 0 0 0 1 3 5 6 − 1 6 0 − 1 12 7 12 7 12 − 1 12 ] , R a n k ( A ) = R a n k ( [ 1 3 − 7 6 11 6 − 1 6 5 6 1 3 1 3 5 6 − 1 6 ] ) (3.10)

Here, the rank of matrix A is 3, so it is so impossible to get the solution of C j like WENO5 [

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 ω i = C i + O ( h ) , so the optimal coefficients can be tested by replacing ω j by C j .

f ^ j + 1 / 2 + ≈ 0.2 ( 1 12 g j − 2 − 5 12 g j − 1 + 13 12 g j + 1 4 g + 1 ) + 0.05 ( 1 3 g j − 2 − 7 6 g j − 1 + 11 6 g j ) + 0.3 ( − 1 6 g j − 2 + 5 6 g j − 1 + 1 3 g j ) + 0.15 ( 1 3 g j − 2 + 5 6 g j − 1 − 1 6 g j ) + 0.3 ( 1 4 g j − 1 + 13 12 g j − 5 12 g j + 1 + 1 12 g j + 2 ) (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,

I S 0 = ( g j − 2 ( 267 g j − 2 − 1642 g j − 1 + 1602 g j − 494 g j + 1 ) + g j − 1 ( 2843 g j − 1 − 5966 g j + 1922 g j + 1 ) + g j ( 3443 g j − 2522 g j + 1 ) + 547 g j + 1 2 ) / 240 I S 1 = 13 12 ( g j − 2 g j − 1 + g j − 2 ) 2 + 1 4 ( 3 g j − 4 g j − 1 + g j − 2 ) 2 I S 2 = 13 12 ( g j + 1 − 2 g j + g j − 1 ) 2 + 1 4 ( g j + 1 − g j − 1 ) 2

I S 3 = 13 12 ( g j + 2 − 2 g j + 1 + g j ) 2 + 1 4 ( 3 g j − 4 g j + 1 + g j + 2 ) 2 I S 4 = ( g j − 1 ( 547 g j − 1 − 2522 g j + 1922 g j + 1 − 494 g j + 2 ) + g j ( 3443 g j − 5966 g j + 1 + 1602 g j + 2 ) + g j + 1 ( 2843 g j + 1 − 1642 g j + 2 ) + 267 g j + 2 3 ) / 240 (3.14)

Here, the weights ω k m :

ω k = α k / ∑ k = 0 4 α k , k = 0 , 1 , 2 , 3 , 4 (3.15)

α k = C k / ( ε + I S k ) , k = 0 , 1 , 2 , 3 , 4 (3.16)

To increase the accuracy of these weights, consider the functions [

g k ( ω ) = ω ( ω ¯ k + ω ¯ k 2 − 3 ω ¯ k ω + ω 2 ) ω ¯ k 2 + ω ( 1 − 2 ω ¯ k ) , ω ¯ k ∈ ( 0 , 1 ) (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 ( 0 ) = 0 , g k ( 1 ) = 1 , g k ( ω ¯ k ) = ω ¯ k , g ′ k ( ω ¯ k ) = 0 , g ″ k ( ω ¯ k ) = 0 .

A more accurate approximation of the weights is given by

α k ∗ = g k ( w k ) (3.18)

The modified weights are defined according to

w k ∗ = α k * ∑ i = 0 4 α i * (3.19)

which satisfies ∑ i = 0 4 w k ∗ = 1 , Here a superscript ( ∗ ) has been added to signify the improved weights, which is called MWENO.

Calculations are carried out for three different examples of compressible ﬂows. 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 [

( ρ , u , p ) T = { ( 1 + 0.2 sin ( 5 x ) , 0 , 1 ) x ≥ − 4.0 ( 3.857143 , 2.629369 , 10.33333 ) T x < − 4.0 .

The computational domain, − 5 ≤ x ≤ 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

Example 2: Interacting blast waves.

The interacting blast wave example [

( ρ , u , p ) T = { ( 1 , 0 , 1000 ) 0 ≤ x < 0.1 ( 1 , 0 , 0.01 ) 0.1 ≤ x < 0.9 ( 1 , 0 , 100 ) T 0.9 ≤ x ≤ 1 , 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

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

This problem was initially proposed and studied in detail by Woodward and Colella [

is set: u ( x , y , 0 ) = { u L , y ≥ h ( x , 0 ) u R , y < h ( x , 0 ) , here, left and right state values and shock wave’s height:

u L = ( 8 , 57.1597 , − 33.0012 , 563.544 ) , u R = ( 1.4 , 0 , 0 , 2.5 ) , h ( x , t ) = 3 ( x − 1 6 ) − 20 t ,

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

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.

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

Guo, S.G. and Li, W. (2021) Improved Weighted Essentially Non-Oscillatory Schemes by Mixed Stencils. Open Journal of Fluid Dynamics, 11, 153-165. https://doi.org/10.4236/ojfd.2021.113009