- Email: [email protected]
The 6th-order weighted ENO schemes for hyperbolic conservation laws
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
Computers and Fluids 0 0 0 (2018) 1–12
Contents lists available at ScienceDirect
Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid
The 6th-order weighted ENO schemes for hyperbolic conservation laws Fuxing Hu Department of Mathematics, Huizhou University, Huizhou, Guangdong 516007, PR China
a r t i c l e
i n f o
Article history: Received 21 November 2017 Revised 14 March 2018 Accepted 13 July 2018 Available online xxx MSC: 65M06 65M20
a b s t r a c t In this work, a 6th-order central-upwind weighted essentially non-oscillatory scheme is devised by additionally introducing a large stencil. The constrained condition of choosing this large stencil is that it should contain upwind information to raise the numerical stability of schemes. The implementation of scheme presented here is as compact as the classical weighted essentially non-oscillatory scheme. Several numerical examples show the robust and high-resolution performances of the scheme which is comparable with the recently developed 6th-order weighted essentially non-oscillatory scheme. © 2018 Elsevier Ltd. All rights reserved.
Keywords: WENO scheme Spurious oscillation Hyperbolic conservation laws
1. Introduction The distinguishing characteristic of nonlinear hyperbolic conservation laws is that their solutions may contain discontinuities such as shock waves even for smooth initial conditions. Due to this characteristic, the high-order numerical methods constructed by traditional ways could generate spurious oscillations around the discontinuous regions. Thus, to cure these oscillations, some extra techniques are needed to devise the high-order numerical methods. Several high-order numerical methods, e.g. essentially nonoscillatory schemes (ENO) [1,2], weighted ENO schemes (WENO) [3,4], discontinuous Galerkin finite element methods (DG FEM) [5,6] and so on, are capable of producing satisfactory numerical results. The classical WENO schemes [4], as a kind of the most popular schemes, have attracted a lot of attention due to the feature of easy implementation and high resolution. Several improvements have been developed to decrease the numerical dissipation of the classical WENO schemes. One kind of improvement, e.g. [7,8] and so on, is to modify the weights and distribute a little more weights to the less smooth stencils. Another improvement, e.g. [9–13], is to couple the low-dissipation scheme which is applied on the smooth region and the WENO scheme which is applied on the non-smooth region. Recently, based on the often used 5th-order WENO scheme, several 6th-order central-upwind WENO schemes [14–16] have been
developed by adding one more stencil. The idea to devise the 6thorder WENO originates from [14] in which the schemes are designed to maximize order of accuracy and bandwidth, and meanwhile minimize dissipation. By including an additional downwind stencil, the overall central schemes WENO-SYMOO and WENOSYMBO in [14] are constructed on the basis of optimal order of accuracy and optimal bandwidth efficiency, respectively. However, the order of WENO-SYMBO scheme degenerates even in smooth regions, and the WENO-SYMOO scheme introduces numerical instabilities near the contact surfaces even when only moderate discontinuities are involved. In [15], with the same stencils as [14], the authors define the smoothness indicator of downwind stencil using all six grid values. To achieve the optimal order of accuracy, a reference smoothness indicator τ is devised which is similar to that first presented in [8]. In addition, it is found that this scheme needs extra artificial dissipation to maintain the numerical stability. In [16], another 6th-order scheme is developed which is analogous with that in [15] but with different smoothness indicator for the newly-introduced downwind stencil and no additional artificial dissipation. We will choose this scheme as comparison since it is recently developed and easily implemented. In numerical simulations, it is found that this scheme cannot tolerate large CFL number. That is, with large CFL number, the scheme generates visible oscillations around discontinuities. Instead of disappearing with increasing cell number, the oscillations show the apparent trend of growth.
E-mail address: [email protected] https://doi.org/10.1016/j.compfluid.2018.07.008 0045-7930/© 2018 Elsevier Ltd. All rights reserved.
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 2
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
To alleviate (or eliminate) the effect of large CFL number on the stability of scheme, we introduce a different stencil into the classical 5th-order WENO scheme [4]. Since the introduction of downwind stencil may lead to certain numerical instability especially for large CFL number, the stencil introduced should contain upwind information. Using the smoothness indicator consisted of all six grid values to measure the smoothness of the downwind stencil consisted of three grid values, as did in [14–16], is rather abrupt. The smoothness indicator of newly-introduced stencil here is directly measured by the formulation presented by Jiang and Shu in [4]. In general, to make the scheme arrive the 6th-order of accuracy and each stencil contain upwind information, this new stencil needs to be larger than original stencils in theory. The nonlinear weights we get cannot satisfy the sufficient condition of retaining the optimal order, by Taylor analysis, and in fact only the 5thorder of accuracy can be obtained. But, this 5th-order scheme still has the smaller errors and higher resolution when compared with the classical 5th-order WENO scheme [4]. In order to cure the degeneration in accuracy, there are two reliable methods generally: embodying a mapping function [7,17,18] into the construction of weights or devising a reference smoothness indicator [8]. In addition, the implementation of schemes presented here is compact like the schemes in [7] or [8] and no artificial parameter is used to control the resolution of solutions. The organization of this paper is as follows. In Section 2, the recently developed 6th-order WENO scheme [16] is reviewed and we also investigate the numerical oscillations generated by this scheme under large CFL number. In Section 3, we present the new 6th-order WENO schemes. The mapping function in [18] is used in order to achieve the optimal accuracy. Also, as another way, a reference smoothness indicator is also devised to recover the optimal accuracy. Finally, in Section 4, several common examples are used to demonstrate the numerical performances of the new 6th-order WENO schemes. 2. The 6th-order WENO-CU6 scheme
In this section, we briefly review the recently developed 6thorder central-upwind WENO scheme [16] for solving hyperbolic systems of conservation laws
( u1 , · · ·
a ≤ x ≤ b,
t ≥ 0,
, um )T
(1)
where u = and g(u ) = system (1) is hyperbolic if Jacobian matrix eigenvalues,
(g1 (u ), · · · gm (u ))T . dg du
We say the has real and distinct
k = 1, · · · , m
(3)
where the positive flux f jk,+ = f jk,−
=
1 k 2(fj
− α k vkj ). The scalar
1 k 2(fj α k is
+ α k vkj ) and negative flux
for all values of u. First, the computational domain [a, b] is divided into N cells a = x1/2 < x3/2 < · · · < xN−1/2 < xN+1/2 = b. Denote the interval [xi−1/2 , xi+1/2 ] by Ii and the length of the cell Ii by xi . For simplicity, we will generally assume a uniform grid, but this is not required. Solving the conservation laws (1) with a conservative difference approximation to the spatial derivative
(2)
where ui (t) is the numerical approximation to the column vector u(xi , t). what to do in the next is how to compute the numerical flux gˆi±1/2 properly. Following the procedure 2.10 in [19], we use the way of characteristic-wise decomposition to project the relevant vectors uj and g(uj ) into the local characteristic fields, where j is in a neighborhood of i. Denote by vj and fj the obtained vectors in the local characteristic fields. To avoid the presentation of
taken as α = max λk (ui ) 1≤i≤N
over the global range of u. The WENO-CU6 reconstruction is used ± to computed each component fˆik, at xi+1/2 and then we obtain +1/2 k, + k, − k fˆ = fˆ + fˆ . Finally, transforming the fˆi+1/2 in characteri+1/2
i+1/2
i+1/2
istic fields back to the physical space, we obtain the numerical flux gˆi+1/2 . In general, the difference between the different WENO schemes lies on the reconstruction of fˆk,± . Here we only give the i+1/2
+ − procedure of obtaining the positive numerical flux fˆik, , for fˆik, +1/2 +1/2 a similar way can be used. For simplicity, we will omit the super± script k in fˆik, in the following procedure. +1/2 In [16], the positive numerical flux fˆi++1/2 can be achieved by the following three steps: 1. Compute the smoothness indicators of four substencils S0 = {Ii−2 , Ii−1 , Ii }, S1 = {Ii−1 , Ii , Ii+1 }, S2 = {Ii , Ii+1 , Ii+2 }, S3 = {Ii+1 , Ii+2 , Ii+3 } by
βr =
2 l=1
xi+1/2 xi−1/2
x2l−1 ( p(rl ) (x ))2 dx,
r = 0, 1, 2
(4)
and
β3 =
5 l=1
xi+1/2 xi−1/2
x2l−1 ( p(3l ) (x ))2 dx
(5)
where pr (x), r = 0, 1, 2, are the 2nd-degree interpolation polynomials over the interval [xi+r−5/2 , xi+r+1/2 ] and p3 (x) is the 5th-degree interpolation polynomial over the whole interval [xi−5/2 , xi+7/2 ]. 2. Compute the nonlinear weights wr by the following formulas
ar , ar = dr C + a0 + a1 + a2 + a3 1 = β3 − (β0 + 4β1 + β2 ) 6
τ
τ
βr + ε
, (6)
and where d0 = 1/20, d1 = 9/20, d2 = 9/20 and d3 = 1/20 are the optimal weights and ε is a parameter used to avoid a division by zero in the denominator and set to be ε = 10−40 in all numerical examples. 3. Compute the convex combination of four substencils,
fˆi++1/2 = w0 fˆi++1,0/2 + w1 fˆi++1,1/2 + w2 fˆi++1,2/2 + w3 fˆi++1,3/2
λ1 (u ) < λ2 (u ) < · · · < λm (u ),
1 dui (t ) =− gˆi+1/2 − gˆi−1/2 , dt x
f jk = f jk,+ + f jk,− ,
wr =
2.1. Review the construction of the WENO-CU6 scheme
ut + g(u )x = 0,
entropy violating solutions, for the each component of characteristic variables, the simple scalar Lax–Friedrichs splitting is used to divide the each component f jk of fj into
(7)
,r fˆi++1 /2
where are the 2nd-order numerical fluxes obtained on stencil Sr . Substituting fˆi++1/2 and fˆi−+1/2 obtained by the above procedure into (3), we can get the numerical flux fˆi+1/2 at the boundary of cell xi+1/2 . Similarly, we can get the numerical flux fˆi−1/2 at xi−1/2 and finally form the scheme (2). The discretization of the ODEs (2) in time obtained from spatial discretization using the methods of lines approach are solved by the third-order TVD Runge–Kutta method [20,21]
u(1) = un + tL(un ), 3 1 u ( 2 ) = un + u ( 1 ) + 4 4 1 n 2 (2 ) n+1 u = u + u + 3 3
1 tL u(1) , 4 2 tL u(2) , 3
where L is the spatial discretization operator and at time tn .
(8) un
is the solution
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
3
Fig. 1. The numerical solutions are computed by the WENO-CU6 for the initial value problem (27) and (29) with CFL = 0.6 and C = 20. N denotes the cell number.
2.2. The problem appeared in WENO-CU6 scheme The WENO-CU6 scheme is developed by adding a downwind stencil into the classical 5th-order WENO scheme [4]. To alleviate the instability brought by downwind information, the smoothness indicator of this stencil is defined on whole stencil S = {Ii−2 , Ii−1 , Ii , Ii+1 , Ii+2 , Ii+3 }. However, the scheme only achieves the 4th-order accuracy if we directly define the nonlinear weights as that in [4]
wr =
ar , a0 + a1 + a2 + a3
ar =
dr , (βr + ε )2
ously, in the case of C = 1, the scheme admits CFL number up to 0.5 which is relaxed when compared with the case of C = 20. In addition, we find that from Fig. 2, with the same CFL number, e.g. CFL = 0.6, the solution with larger parameter C has larger oscillation. For much clearer to see this trend, the Fig. 3 shows the oscillations behind the discontinuity with the fixed CFL=0.6 but varying the parameter C. And now, we can easily find that the oscillations grow larger with the choice of larger parameters C. 3. The construction of new 6th-order WENO schemes
r = 0, 1, 2, 3.
(9)
This claim can be easily confirmed by symbolic operations or numerical tests. It is noted that, although the scheme with nonlinear weights (9) cannot achieves optimal order, it can tolerate a reasonable CFL number as that of the classical 5th-order WENO scheme [4]. To cure the degeneration of order of accuracy, the authors apply the technique proposed in [8] to define the weights as in (6). And a new parameter C = 20 is introduced in (6) for increasing the contribution of optimal weights when the smoothness indicators have comparable magnitudes. But, for the WENO-CU6 scheme, the construction of weight in (6) and the introduction of parameter C = 20 both result in the weak tolerance of large CFL number. That is, the scheme will generate oscillations around the discontinuities if one chooses a large CFL number. Considering the simple example (27) and (29), we first set C = 20 and t = CFLx where CFL = 0.6. From the Fig. 1, one can find there are obvious oscillations behind the discontinuities at x = −0.4 and x = −0.2. The close-up view around x = −0.2 is shown in right subgraph of Fig. 1. In addition, we find the oscillations do not yet disappear as increasing cell number to N = 3200, but present the trend of growth. Next, we show the influence of different CFL numbers on the WENO-CU6 scheme with C = 20 and N = 3200. The left subgraph in Fig. 2 shows the numerical discontinuity around x = −0.2 with CFL =0.3, 0.4, 0.5, 0.6, 0.7. The oscillations appear visibly when CFL = 0.5 and increase gradually with increasing CFL number. Hence, at least for this example, a safe choice of CFL number is 0.4 (or even smaller 0.3). As for the parameter C, it also affects the choice of CFL number. Although the larger parameter C generally produces the slightly higher-resolution numerical solutions, it often imposes more requirement on CFL number. In contrast, the lesser parameter C makes the choice of CFL number more relaxed. For example, decreasing the parameter C to 1, it will be found that the scheme can tolerate the larger CFL number. The right subgraph in Fig. 2 are computed with C = 1 and N = 3200 which has same cell number with left figure and the only difference is the choice of C. Obvi-
In this section, we present the procedure of constructing the new 6th-order WENO scheme. As other 6th-order WENO schemes shown in [14–16], we also add an extra stencil S3 shown in left subgraph of Fig. 4 (or S0 in right subgraph). The difference is that the new stencil introduced here contains upwind information which may aid in numerical stability. From Fig. 4, to construct the positive numerical flux fˆi++1/2 and the negative one fˆ− with the 6th order of accuracy, the steni+1/2
cils
S0+ = {Ii−2 , Ii−1 , Ii }, S1+ = {Ii−1 , Ii , Ii+1 }, S2+ = {Ii , Ii+1 , Ii+2 }, S3+ = {Ii , Ii+1 , Ii+2 , Ii+3 }. and
S0− = {Ii−2 , Ii−1 , Ii , Ii+1 }, S1− = {Ii−1 , Ii , Ii+1 }, S2− = {Ii , Ii+1 , Ii+2 }, S3− = {Ii+1 , Ii+2 , Ii+3 }.
are used, respectively. Using each stencil Sr± , r = 0, · · · , 3, one gets four low-order numerical fluxes
7 1 11 + fˆi++1,0/2 = fi+−2 − fi+−1 + f , 3 6 6 i 5 − 1 − 13 − 1 − 0 fˆi−, = f − f + f + fi+1 , +1/2 12 i−2 12 i−1 12 i 4 ˆf +,1 = − 1 f + + 5 f + + 1 f + , i+1/2 6 i−1 6 i 3 i+1 1 5 1 −, 1 fˆi+1/2 = − fi−−1 + fi− + fi−+1 , 6 6 3 ˆf +,2 = 1 f + + 5 f + − 1 f + , i+1/2 3 i 6 i+1 6 i+2 1 1 5 2 fˆi−, = fi− + fi−+1 − fi−+2 , +1/2 3 6 6 5 + 1 13 + 1 + fˆi++1,3/2 = fi+ + f − f + f , 4 12 i+1 12 i+2 12 i+3 7 11 − 1 3 fˆi−, = f − f− + f− , +1/2 6 i+1 6 i+2 3 i+3
(10)
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 4
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
Fig. 2. The left is computed with C = 20 and right with C = 1.
Inserting (10) into (11), through simple computing one gets the linear optimal weights
d0+ =
1 9 3 1 , d1+ = , d2+ = , d3+ = , 20 20 10 5
and
d0− =
1 3 9 1 , d1− = , d2− = , d3− = . 5 10 20 20
Remark 1. In above procedure, taking the positive flux fˆi++1/2 for
example, we choose S3+ = {Ii , Ii+1 , Ii+2 , Ii+3 } as the fourth stencil. Actually, the stencil S3+ = {Ii−1 , Ii , Ii+1 , Ii+2 , Ii+3 } consisted of five points can also be chosen as the candidate. Under this choice, the overall scheme is also stable numerically. However, through comparison, it is found that the first choice has higher resolution and lower computational cost. So we only discuss the scheme formed by the first choice in this paper.
Fig. 3. The numerical solutions are computed by WENO-CU6 scheme with N = 1600. We fix CFL = 0.6 but vary parameter C.
Following the idea of the classical WENO scheme, a nonlinear convex combination of low-order numerical flux fˆi±,r similar to +1/2 (7) is needed to control the spurious oscillations around the discontinuities ± ˆ±,1 ± ˆ±,2 ± ˆ±,3 ˆ±,0 fˆi±+1/2 = w± 0 f i+1/2 + w1 f i+1/2 + w2 f i+1/2 + w3 f i+1/2 ,
,3 0 and the numerical flux fˆi++1 computed on S3+ and fˆi−, computed /2 +1/2
on S0− have 4th-order accuracy. Through the linear convex combination of four low-order stencils S0 , S1 , S2 and S3 , one wishes to get the 6th-order central interpolation on 6-point stencil S
1 + fi−2 − 8 fi+−1 + 37 fi+ + 37 fi++1 − 8 fi++2 + fi++3 , 60 d− fˆ−,0 + d− fˆ−,1 + d− fˆ−,2 + d− fˆ−,3 =
1
i+1/2
2
i+1/2
3
a± 0
+
a± 1
a± r , + a± + a± 2 3
a± r =
(β
dr± , + ε )2
± r
r = 0, 1, 2, 3. (13)
The smoothness indicators are computed by (2.61) in [19]
β2+
i+1/2
1 − fi−2 − 8 fi−−1 + 37 fi− + 37 fi−+1 − 8 fi−+2 + fi−+3 . 60
w± r =
β1+
=
i+1/2
where the nonlinear weights w± r can be obtained as that in [4]
2 1 + 2 13 + f − 2 fi+−1 + fi+ + f − 4 fi+−1 + 3 fi+ , 12 i−2 4 i−2 2 1 + 2 13 + = fi−1 − 2 fi+ + fi++1 + fi−1 − fi++1 , 12 4 2 13 + 1 + + + 2 = f − 2 fi+1 + fi+2 + 3 fi − 4 fi++1 + fi++2 , 12 i 4 2 1 = 11 fi+ − 18 fi++1 + 9 fi++2 − 2 fi++3 36
β0+ =
d0+ fˆi++1,0/2 + d1+ fˆi++1,1/2 + d2+ fˆi++1,2/2 + d3+ fˆi++1,3/2
0
(12)
(11)
β3+
Fig. 4. Candidate stencils Sr for the numerical flux fˆi++1/2 (left) and fˆi−+1/2 (right).
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
2 13 + 2 fi − 5 fi++1 + 4 fi++2 − fi++3 12 2 781 + + fi − 3 fi++1 + 3 fi++2 − fi++3 , 720 +
(14)
β0− =
β
− 1
β2− β3−
w r − d r = O x 4 , r = 0 , 1 , 2 , w 3 − d 3 = O x 3 .
(16)
In order to more easily implement, the sufficient conditions (16) can be extended to the more compact formulas
w r − d r = O x 4 , r = 0 , · · · , 3 .
(15)
where, for β3+ and β0− , we have used the more compact explicit formulas with the sum of squares in [18]. Finally, by similar way we can obtain the negative numerical flux fˆi−+1/2 and finally the numerical flux fˆi+1/2 at xi+1/2 is
fˆi+1/2 = fˆi++1/2 + fˆi−+1/2 . Remark 2. Despite the fact that the formulas of linear weights dr , low-order numerical fluxes fi±,r and smoothness indicators are all +1/2 not symmetrical, the overall scheme to compute the numerical flux fˆi+1/2 at cell boundary xi+1/2 is symmetrical at smooth region. At non-smooth region, the scheme is upwind. Hence this kind of 6thorder scheme is also called central-upwind scheme. Up to present, we have presented the numerical flux fˆi+1/2 at xi+1/2 . One can find that the construction of this scheme is as compact as the classical WENO scheme in [4,22]. It is supposed to obtain the 6th-order approximation
fˆi+1/2 − fˆi−1/2 = fx (u(xi , t )) + O (x6 ) x since we have constructed fˆi+1/2 and fˆi−1/2 with 7 flux values
f (ui−3 ), f (ui−2 ), f (ui−1 ), f (ui ), f (ui+1 ), f (ui+2 ), f (ui+3 ). Unfortunately, the mixed smoothness indicators of three 2nddegree and one 3rd-degree polynomials can only guarantee
generates much higher-resolution solutions and presents lower numerical errors when compared with the classical 5th-order WENO scheme. As in [7], we use similar way to derive the sufficient conditions for making the scheme achieve the 6th order of accuracy
and
2 1 − 2 fi−2 − 9 fi−−1 + 18 fi− − 11 fi−+1 36 2 13 − + f − 4 fi−−1 + 5 fi− − 2 fi−+1 12 i−2 2 781 − + fi−2 − 3 fi−−1 + 3 fi− − fi−+1 , 720 2 1 − 2 13 − = f − 2 fi− + fi−+1 + f − 4 fi− + 3 fi−+1 , 12 i−1 4 i−1 2 1 − 2 13 − = fi − 2 fi−+1 + fi−+2 + fi − fi−+2 , 12 4 2 13 − 1 − − − 2 = f − 2 fi+2 + fi+3 + 3 fi+1 − 4 fi−+2 + fi−+3 , 12 i+1 4
5
2
w r = d r + O x , r = 0 , · · · 3 , which make the scheme achieve the 5th-order of accuracy at most when the solution is smooth and meanwhile no critical point appeared. Detailedly, by Taylor expansion, we indeed only get
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 5 . x At critical points, the order of accuracy will degenerate seriously. For example, for f = 0, we get
wr = dr + O (x ), r = 0, · · · 3, and
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 3 . x
(17)
So, to recover the optimal order of accuracy, we only need to make the nonlinear weights approximate the linear weights with O (x4 ). In following parts of this section, we will present two methods to recover the optimal order of accuracy of the scheme: one is to apply a mapping function, as in [7], on the nonlinear weights computed above, and the other is to construct a reference smoothness indicator as in [8]. The common purpose of these two remedy methods is to make the nonlinear weights converge to the linear weights fast enough as x → 0. 3.1. The mapping function The mapping function devised in [7]
g¯ (wr ) = dr +
dr2
( wr − dr )3 + ( 1 − 2dr )wr
(18)
is enough to make the nonlinear weights satisfy the sufficient condition (17). At the smooth region with no critical point, since the nonlinear weights achieve wr − dr = O (x2 ), r = 0, · · · , 3, applying the mapping function (18) on nonlinear weights gives
w r − d r = O x 6 , r = 0 , · · · , 3 .
(19)
However, when at critical point f = 0 which leads to wr = dr + O (x ), this mapping function only ensures wr = dr + O (x3 ) which is not enough to guarantee the optimal order. Actually, in this case, the scheme only achieves 5th-order accuracy with the mapping function (18). Hence, the generalized mapping function in [18] is preferred to be used
g( w r ) = d r +
(wr − dr )7 C0 + C1 wr + C2 w2r + C3 w3r
,
(20)
where C0 = dr6 , C1 = −7dr5 , C2 = 21dr4 , C3 = (1 − dr )6 − (C0 + C1 + C2 ). This function apparently recovers the optimal order of 6thorder scheme even at critical points f = 0. Besides, when compared with the mapped WENO in [7], it allocates much smaller weights to discontinuous stencils which could further compress the spurious oscillations. After obtaining the mapping function (20), the improved nonlinear weights which satisfy the sufficient condition (17) are computed by
aM r wM r = k=3 k=0
aM k
,
aM r = g( w r ) ,
r = 0, · · · , 3.
Hence, for the numerical tests in Section 4, we always choose the mapping function (20) unless otherwise stated.
That is, this “6th-order” scheme only achieves 3rd-order of accuracy around critical points with f = 0.
3.2. The reference smoothness indicator
Remark 3. Although this “6th-order” scheme cannot yet achieve the optimal accuracy, by numerical experiment, we find that it still
In this section, we plan to present the other method appeared in [8] to recover the optimal order of 6th-order scheme. The idea
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 6
of this correction is to construct a high-order reference smoothness indicator. Expanding the smoothness indicators (14), we get
13 β0 = fi2 x2 + f 2 − 12 i x 5 + O x 6 , 13 β1 = fi2 x2 + f 2 + 12 i 13 β2 = fi2 x2 + f 2 − 12 i x 5 + O x 6 , β3
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
2 f f x 4 + 3 i i
1 2
fi fi(4) −
13 f f 6 i i
1 f i f i x 4 + O x 6 , 3 1 2 13 fi fi x4 + − fi fi(4) + fi fi 3
2
6
(21)
13 2 4 1 (4) 5 = fi2 x2 + f x + f i f i x + O x 6 . 12 i 2
(21)
Since the reference smoothness indicator here is preferred to contain all 6 points, the highest-order one we can obtain is the following linear combination
1 τ = β3 − (β0 + 4β1 + β2 ) = O x5 .
(22)
6
Remark 4. Actually, there exists another combination τ = ((β1 − β0 ) − 32 (β3 − β2 ) which also arrives O (x5 ). However, the numerical results with this reference smoothness indicator are worse than that with (22). So, in this work, we only consider the reference smoothness indicator (22). Finally, the nonlinear weights can be computed similarly by
aZr
wZr = k=3 k=0
aZk
,
aZr = dr 1 +
p
τ
,
βr +
r = 0, · · · , 3, (23)
where the parameter p is used to control the order of accuracy and dissipation. For the case of p = 1, through simple computation, it seems to get that
wZr = dr + O x3 .
(24)
Obviously, this “doesn’t” satisfy the sufficient condition (17). But, it was not as straightforward as we were led to believe. By detailed Taylor expansion, one can get
aZr = dr 1 + Ax3 + Bx4 + O x5 ,
Remark 5. It is noted that this 6th-order scheme has less computational cost compared with the scheme in [16]. The difference of computational cost between these two schemes mainly lies on the computation of smoothness indicator of new introduced stencil. Comparing the smoothness indicator (14) (or (15)) and the formula (25) in [16] (which we used is the more compact the sum of squared), the latter has two sum of squared more than the former. Up to now, two methods have been used to recover the optimal order of accuracy. The use of mapping function can achieve the optimal order efficiently even around the critical points. The flaw of this method is the raise of computational cost compared to the classical WENO. The use of reference smoothness indicator with the parameter p = 1 also achieves the optimal order efficiently but not including the critical points. A compromised way to recover the 6th-order accuracy at critical points is to enlarge the parameter p. Nevertheless, this will increase the dissipation of the scheme and the scheme presents the low resolution around the discontinuities which is confirmed in [8]. 4. Numerical results In this section, several 1D and 2D numerical examples are used to test the 6th-order WENO schemes. The 6th-order WENO scheme with the mapping function is denoted by WENO-M6 and with reference smoothness indicator is denoted by WENO-Z6. For comparison, we choose the recently developed 6th-order scheme in [16] denoted by WENO-CU6 and the parameter C = 20. Although the similar parameter C can also be introduced in the WENO-Z6 scheme to increase the resolution, we abandon it for the compaction of scheme and expect to make the numerical solutions look monotone. Overall, the WENO-CU6 and WENO-Z6 present similar errors and accuracy if the same parameters C and p are set in schemes. 4.1. The linear advection equation
f (4 ) 2f
In this section, we consider the simple case
ut + ux = 0, x ∈ [−1, 1], t ≥ 0,
and
u(x, 0 ) = sin
Finally, the nonlinear weights in (23) arrive
d r 1 + A x 3 + B x 4 + O x 5 1 + A
x3
+ B
x4
+O
x5
πx −
sinπ x
(28)
π
and
5
= d r + O x .
(26)
Hence, the nonlinear weights constructed in (23) actually approach the linear weights with O (x5 ) and conform to the requirement (17) of retaining the 6th-order accuracy at smooth region with no critical point. For the 1st-order critical point, by symbolic operation the nonlinear weights achieve
wZr = dr + O x3
(27)
with the following two initial conditions to test the schemes,
2 4 (5) 91 3 (4) 481 3 (3) 2 f f − f f f + f f 36 720 B= 3 . 5 f
wZr =
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 6 . x
(25)
where
A=
To alleviate the degeneration of accuracy at 1st-order critical point, the favourite way is to enlarge the parameter p in (23). For instance, in the case of p = 2, around the smooth region with the 1st-order critical point, the scheme recovers the optimal order
⎧ ⎪ ⎪(G(x, β , z − δ ) + G(x, β , z + δ ) + 4G(x, β , z ) )/6 ⎪ −0.8 ≤ x ≤ −0.6, ⎪ ⎪ ⎪ 1 − 0.4 ≤ x ≤ −0.2, ⎨ 1 − |10x − 1|0.0 ≤ x ≤ 0.2, u(x, 0 ) = ⎪ ⎪ (F (x, α , a − δ ) + F (x, α , a + δ ) + 4F (x, α , a ) )/6 ⎪ ⎪ ⎪ ⎪ ⎩ 0.4 ≤ x ≤ 0.6, 0
otherwise,
and the order of accuracy of scheme is
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 5 . x
(29) where
G(x, β , z ) = exp(−β (x − z )2 ),
F (x, α , a ) =
max(1 − (α (x − a ))2 , 0 ), and the parameters are given by
a = 0.5, z = −0.7,
δ = 0.005, α = 10, β =
ln 2 . 36δ 2
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
7
Fig. 5. The l1 and l∞ errors are computed by the 6th-order linear, WENO-CU6, WENO-M6 and WENO-Z6 schemes, respectively, for the initial value problem (27) and (28). Table 1 The l1 errors and rate of convergence are shown around the square wave in the problem (27) and (29) with CFL = 0.3 and the output time t = 2.
x 2.00e−02 1.00e−02 5.00e−03 2.50e−03
WENO-CU6
WENO-M6
WENO-Z6
l1 error
r
l1 error
r
l1 error
r
4.07342e−02 2.21530e−02 1.23156e−02 6.85770e−03
0.8787 0.8470 0.8447
4.83906e−02 2.54545e−02 1.39010e−02 7.58246e−03
0.9268 0.8727 0.8745
4.53699e−02 2.48965e−02 1.37859e−02 7.62225e−03
0.8658 0.8527 0.8549
The first condition with smooth solution is used to test the order of accuracy and the second with discontinuous and high-frequency waves is used to test the essentially non-oscillatory property of schemes. All three methods are used to compute the initial problem (27) and (28). The time step is t = x2 and ε = 10−40 . The l1 and l∞ errors of schemes are shown in Fig. 5 and all the WENO schemes present the 6th order of accuracy and show the almost same errors as the linear central 6th-order scheme. Here, the parameter p in (23) is chosen to be 2 for recovering the order at critical points (for all other examples p = 1). It is noted that the error of WENO-Z6 will be enlarged slightly if p = 1 for this smooth example. Of course, the introduction of C = 20 in WENO-Z6 will decrease the error to that of WENO-CU6 scheme even p = 1. Besides the accuracy of 6th-order schemes for the smooth solution, the rate of convergence around the discontinuity is also tested here. For nonlinear shock waves, the high-resolution schemes generally produce a 1st-order of accuracy. However, we are interested in the case of linear (or linearly degenerate) discontinuities which will dominate the global rate of convergence. As the analysis in [23], the 6th-order WENO schemes here should produce the expected rates of convergence O (x6/7 ). We choose the problem (27) and (29) which contains Gaussian, square, triangle and semiellipse waves to test the rate of convergence around the square (which is linearly discontinuous). To eliminate the influence of other three waves, we only present the numerical errors around [−0.5, −0.1]. Table 1 shows the rate of convergence of the 6thorder schemes. For the problem (27) and (29), the initial condition contains Gaussian, square, triangle and semi-ellipse waves. Fig. 6 shows the solutions computed by three methods with cell number N = 200 and time step t = 0.3x. They all present high-resolution solutions around Gaussian, triangle and semi-ellipse waves except for the square. But, from the right amplified subgraph in Fig. 6, the WENO-CU6 produces much better than other two schemes. This
owes to the introduction of parameter C = 20 into the WENO-CU6 scheme and the WENO-Z6 will generate the similar solution with WENO-CU6 if we introduce the parameter C = 20 into WENO-Z6. But, as shown in Section 2.2, it has been found that, in the WENOCU6 scheme, the introduction of C = 20 partly leads to the poor tolerance of large CFL number. That is the WENO-CU6 scheme easily produces visible oscillations behind discontinuity if we increase the CFL number. And the magnitude of oscillation do not disappear (or even alleviate) with increasing the cell number. Through this example, we wonder how the CFL number affects the numerical stability of WENO-M6 and WENO-Z6 schemes. Fig. 7 shows the results of these two schemes at x = −0.2 computed with N = 3200 and several different CFL numbers. In the left subgraph, the WENO-M6 scheme is almost unaffected by CFL number and does not present any oscillation even for CFL = 0.8. However, from the right subgraph, we can find that no oscillation is observed with CFL number up to 0.6 and there are clear oscillation when CFL ≥ 0.7 (here we do not consider the situation when CFL number between 0.6 and 0.7). It is noted that, under the case of C = 1, the WENO-CU6 conserves no oscillation with CFL number up to 0.5. Hence the WENO-Z6 has stronger tolerance on the CFL number than WENO-CU6 scheme at least for this example. And in Section 4.2, we will find that similar conclusion is still hold even for nonlinear Euler systems. In addition, comparing the results between the right subgraph of Fig. 2 (where C = 1 in WENOCU6) and the right subgraph of Fig. 7, the magnitude of oscillation generated by WENO-CU6 is slight larger than WENO-Z6 under the same CFL number.
4.2. 1D Euler equations In this example, we consider 1D Euler equations since one of the main application areas of weighted essentially non-oscillatory
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 8
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
Fig. 6. The numerical solutions are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes, respectively, for the initial value problem (27) and (29) with N = 200 and CFL = 0.3.
Fig. 7. The numerical solutions are computed by the WENO-M6 (left) and WENO-Z6 (right) for the initial value problem (27) and (29) with N = 3200.
schemes is compressible gas dynamics,
ρ ρu
+
E
t
ρu ρ u2 + p ( E + p )u
= 0,
(30)
x
where ρ , u, p, E are density, velocity, pressure and total energy, respectively. The system of equations is closed by the equation of state for an ideal polytropic gas:
E=
γ
p 1 + ρ u2 , −1 2
U (l, i ) = U (l, N ), l = 1, 2, 3; i = N + 1, N + 2, N + 3.
1. (Sod problem)
(ρ , u, p) =
( 1, 0, 1 ), (0.125, 0, 0.1 ),
if − 5 ≤ x < 0, if 0 ≤ x ≤ 5,
(31)
2. (Shock entropy wave interaction problem)
(ρ , u, p) =
⎧
√ ⎨ 27 4 35 31 ⎩
, , , 7 9 3 (1 + 0.2 sin(5x ), 0, 1 ),
if − 5 ≤ x < −4, if − 4 ≤ x ≤ 5, (32)
3. (Blast wave interaction problem)
(1, 0, 10 0 0 ), (ρ , u, p) = (1, 0, 0.01 ), (1, 0, 100 ),
if − 5 ≤ x < −4, if − 4 ≤ x < 4, if 4 ≤ x ≤ 5.
U (l, i ) = U (l, 1 ), l = 1, 2, 3; i = −2, −1, 0, and
where the ratio of specific heats γ = 1.4. The following three initial conditions are considered
The shock tube problem (30) and (31) is called Sod problem [24]. This problem contains a rarefaction wave, a contact discontinuity and a shock. At the output time t = 2, Since waves created by the discontinuity at x = 0 do not reach the boundaries, the zero gradient boundary conditions are applied on the left and right boundaries. Specifically, denote U (l, i ) = (ρ (i ), ρ (i )u(i ), E (i ))T , l = 1, 2, 3, the practical implementations of boundary condition are
(33)
We first test the rate of convergence around the contact discontinuity. Also, to avoid the influence of other waves, we choose to only compute the l1 errors in interval [1,2.75]. Table 2 shows the numerical errors of 6th-order schemes and all the rates of convergence are approximately equal to the expected order O (x6/7 ). In addition, Fig. 8 shows the density distributions of the 6thorder schemes and they all achieve high resolution. These numerical density is calculated with moderate CFL = 0.6. However, the overshoot will appear clearly around discontinuities if CFL number is increased to 0.8. Fig. 9 presents the results of these methods under CFL = 0.8 with cell numbers N = 20 0, 40 0, 80 0, 160 0. For the WENO-CU6 scheme, there exist obvious oscillations behind the contact discontinuity and the magnitude of oscillations has no trend of disappearance with increasing the cell number. Another observation from Fig. 9 is that there are little overshoot behind
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
9
Table 2 The l1 errors and rate of convergence are shown around the square wave in the problem (30) and (31) with CFL = 0.3 and the output time t = 2.
x 1.00e−02 5.00e−03 2.50e−03 1.25e−03
WENO-CU6
WENO-M6
WENO-Z6
l1 error
r
l1 error
r
l1 error
r
1.81614e−03 1.03725e−03 5.73400e−04 3.19969e−04
0.8081 0.8551 0.8416
2.06224e−03 1.17450e−03 6.35782e−04 3.45954e−04
0.8122 0.8854 0.8631
2.03386e−03 1.16065e−03 6.35977e−04 3.48802e−04
0.8096 0.8679 0.8666
Fig. 8. The density profiles of problem (30) and (31) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 200 and CFL = 0.6.
Fig. 9. The density profiles of problem (30) and (31) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 200, 400, 800, 1600 and CFL = 0.8.
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 10
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
Fig. 10. The density profiles of problem (30) and (32) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 400 and CFL = 0.6.
Fig. 11. The density profiles of problem (30) and (33) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 400 and CFL = 0.6.
the shock and they disappear with increasing the cell number. For the WENO-Z6 scheme, there exist only little overshoot behind the shock with low-resolution cell number which disappears as the WENO-CU6 scheme with increasing the cell number. The WENOM6 scheme again behaves robustly and produces look-monotone solutions. Hence, as in the last example (27) and (29), the ability of tolerating the large CFL number of WENO-CU6 is less than that of WENO-M6 and WENO-Z6. The shock entropy-wave interaction problem (30) and (32) in [25] simulates the interaction of a right-moving Mach 3 shock with a wavelike perturbed density whose magnitude is much smaller than the shock. Also, the zero gradient boundary conditions are applied on the left and right boundaries. As a result, a flow field of compressed and amplified wave trails is created right behind the shock. The “exact” solution is computed by the 5th-order WENOJS scheme with N = 6400. From Fig. 10, we can find that all the schemes generate high-resolution results. By amplifying the part of entropy waves, the WENO-CU6 scheme shows slightly higher resolution than the WENO-M6 and WENO-Z6 schemes. The last example is the problem (30) and (33) of two interactive blast waves [26] with reflection boundary conditions. the practical implementations of numerical boundary conditions are
U (l, i ) = (−1 )
l+1
U (l, 1 − i ), l = 1, 2, 3; i = −2, −1, 0
for left boundary and
U (l, i ) = (−1 )l+1U (l, 2N + 1 − i ), l = 1, 2, 3; i = N + 1, N + 2, N + 3
for right boundary. The “exact” solution is also computed by the 5th-order WENO-JS scheme with N = 6400. Again, all the schemes perform with high accuracy and the WENO-CU6 scheme shows slightly higher resolution than the WENO-M6 and WENO-Z6 schemes around x = 2.5 in Fig. 11. Remark 6. Through simulating above several linear scalar equations and nonlinear systems, one can find that these three schemes all generate high-accuracy solution. In general, the WENO-CU6 scheme has slightly better numerical solutions than WENO-M6 and WENO-Z6 schemes. The expense of WENO-CU6 is the spurious oscillations behind the linear discontinuity in some case and the oscillations do not seem to disappear with increasing cell number. This weak tolerance is largely due to the introduction of parameter C. 4.3. 2D Euler equations
⎡
In this section, we consider 2D Euler equations,
⎤ ⎡ ⎤ ⎡ ⎤ ρu ρv ρ 2 ⎢ρ u⎥ ⎢ ρ u + p ⎥ ⎢ ρ uv ⎥ ⎣ρv⎦ + ⎣ ρ uv ⎦ + ⎣ ρv2 + p ⎦ = 0, E t ( E + p )u x ( E + p )v y
(34)
where the equation of state is p = (γ − 1 )(E − 12 ρ (u2 + v2 )). Here we apply the 6th-order schemes to the 2D double-Mach shock reflection problem which contains strong shock waves and
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
JID: CAF
ARTICLE IN PRESS
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
11
Fig. 12. The density profiles are computed respectively by the WENO-CU6, WENO-M6 and WENO-Z6 schemes. 30 equally spaced contours for the density from 1.731 to 20.9705 with N = 800 × 200 and CFL = 0.6.
contact discontinuity. This problem was first proposed by Woodward and Colella [26] and had been taken extensively as a test example for high-order schemes. The computational domain is chosen to be [0, 4] × [0, 1] and the reflective wall lies on the bottom of the computational domain for 16 ≤ x ≤ 4. In the beginning, a right-moving Mach 10 shock is located at x = 16 , y = 0 and makes an angle 60° with the x-axis. For the boundary conditions, the exact postshock condition is imposed for bottom boundary from x = 0 to x = 16 , and the reflective boundary condition is imposed for the rest; the flows are imposed on the top boundary such that there is no interaction with the Mach 10 shock; inflow and outflow boundary conditions are set for the left and right boundaries respectively. The unshocked fluid has a density of 1.4, a pressure of 1. Specifically, the boundary conditions can be implemented numerically as follows. Denote by Nx the
grid number in x-axis direct, Ny in y-axis direct and U (l, i, j ) = (ρ (i, j ), ρ (i, j )u(i, j ), ρ (i, j )v(i, j ), E (i, j ))T , where l = 1, · · · , 4, i = 1, · · · , Nx and j = 1, · · · , Ny . At the left boundary,
U (l, i, j ) = U (l, 1, j ), i = −2, −1, 0; at the right boundary,
U (l, i, j ) = U (l, Nx , j ), i = Nx + 1, Nx + 2, Nx + 3; at the bottom boundary, if the cell center x(i) ≤ 1/6 then
√ U (1, i, j ) = 8, U (2, i, j ) = 33 3, U (3, i, j ) = −33, U (4, i, j ) = 563.5, or else
U (1, i, j ) = U (1, i, 1 − j ), U (2, i, j ) = U (2, i, 1 − j ),
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
JID: CAF 12
ARTICLE IN PRESS
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
U (3, i, j ) = −U (3, i, 1 − j ), U (4, i, j ) = U (4, i, 1 − j ), where j = −2, −1 √, 0; at the top boundary, if the cell center x(i ) ≤ 1/6 + (1 + 20t )/ 3 then
√ U (1, i, j ) = 8, U (2, i, j ) = 33 3, U (3, i, j ) = −33, U (4, i, j ) = 563.5, or else
U (1, i, j ) = 1.4, U (2, i, j ) = 0, U (3, i, j ) = 0, U (4, i, j ) = 2.5, where j = Ny + 1, Ny + 2, Ny + 3. The Fig. 12 shows the density profiles computed by the WENOCU6, WENO-M6 and WENO-Z6 schemes with cell number N = 800 × 200 and output time t = 0.2. The WENO-CU6 scheme gives the better solution than that of WENO-M6 and WENO-Z6 schemes, and the difference between the WENO-CU6 and WENO-M6 is subtle. 5. Conclusions In this work, we present new 6th-order WENO schemes for solving hyperbolic conservation laws. Similar to other 6th-order WENO, we also introduce an additional stencil into the classical 5th-order WENO. The difference is that the new-introduced stencil here contains upwind information and is not completely down-wind stencil. This choice looks more reasonable to control the numerical stability of scheme. A mapping function and reference smoothness indicator are used to make the scheme achieve optimal order, respectively. Then two slightly different schemes, WENO-M6 and WENO-Z6, are developed. We choose the recently developed 6th-order WENO-CU6 scheme as the comparison with the new schemes in this paper. These three schemes all generate the high-resolution solutions and the WENO-CU6 is better relatively under close-up view. Through numerical simulations, we find that the WENO-CU6 scheme has weak tolerance of large CFL number. This weak tolerance is partly due to the introduction of parameter C. The WENO-Z6 also cannot tolerate the large CFL number, but is more robust than the WENO-CU6. The WENO-M6 shows the strong-robust characteristic and is almost unaffected by CFL number. Acknowledgments This work was supported by the Technology Foundation of Huizhou (No. 2017C0403019) and Huizhou University (No. HZUX1201620).
[2] Harten A, Engquist B, Osher S, Chakravarthy S. Uniformly high order essentially non-oscillatory schemes III. J Comput Phys 1987;71:231–303. [3] Liu XD, Osher S, Chan T. Weighted essentially non-oscillatory schemes. J Comput Phys 1994;115:200–12. [4] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phys 1996;126:202–28. [5] Cockburn B, Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Math Comp 1989;52:411–35. [6] Cockburn B, Hou S, Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case. Math Comp 1990;54:545–81. [7] Henrick AK, Aslam TD, Powers JM. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points. J Comput Phys 2005;207:542–67. [8] Borges R, Carmona M, Costa B, Don WS. An improved weighted essentially non-oscillation scheme for hypebolic conservation laws. J Comput Phys 2008;227:3191–211. [9] Adams NA, Shariff K. A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems. J Comput Phys 1996;127:27–51. [10] Pirozzoli S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction. J Comput Phys 2002;178:81–117. [11] Ren YX, Liu M, Zhang H. A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws. J Comput Phys 2003;192:365–86. [12] Kim D, Kwon HK. A high-order accurate hybrid scheme using a central flux scheme and a WENO scheme for compressible flowfield analysis. J Comput Phys 2005;210:554–83. [13] Wang QJ, Ren YX, Sun ZS, Sun YT. Low dispersion finite volume scheme based on reconstruction with minimized dispersion and controllable dissipation. Sci China Phys Mech Astron 2013;56:423–31. [14] Martín MP, Taylor EM, Wu M, Weirs VG. A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence. J Comput Phys 2006;220:270–89. [15] Yamaleev NK, Carpenter MH. A systematic methodology for constructing high-order energy stable WENO schemes. J Comput Phys 2009;228:4248–72. [16] Hu XY, Wang Q, Adams NA. An adaptive central-upwind weighted essentially non-oscillatory scheme. J Comput Phys 2010;229:8952–65. [17] Feng H, Hu F, Wang R. A new mapped weighted essentially non-oscillatory scheme. J Sci Comput 2012;51:449–73. [18] Wang R, Feng H, Huang C. A new mapped weighted essentially non-oscillatory method using rational mapping function. J Sci Comput 2015;67:540–80. [19] Shu CW. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In: Advanced numerical approximation of nonlinear hyperbolic equations. In: Lecture notes in mathematics, 1697. Berlin: Springer-Verlag; 1998. p. 325–432. [20] Shu CW. Total-variation-diminishing time discretizations. SIAM J Sci Statist Comput 1988;9(6):1073–84. [21] Gottlieb S, Shu CW. Total variation diminishing Runge-Kutta schemes. Math Comp 1998;67:73–85. [22] Balsara DS, Shu CW. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. J Comput Phys 20 0 0;160:405–52. [23] Banks JW, Aslam T, Rider WJ. On sub-linear convergence for linearly degenerate waves in capturing schemes. J Comput Phys 2008;227:6985–7002. [24] Sod GA. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J Comput Phys 1978;107:1–31. [25] Shu CW, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes II. J Comput Phys 1989;83:32–78. [26] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comput Phys 1984;54:115–73.
References [1] Harten A, Osher S. Uniformly high order essentially non-oscillatory schemes I. SIAM J Numer Anal 1987;24:279–309.
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
JID: CAF
[m5G;July 19, 2018;13:27]
Computers and Fluids 0 0 0 (2018) 1–12
Contents lists available at ScienceDirect
Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid
The 6th-order weighted ENO schemes for hyperbolic conservation laws Fuxing Hu Department of Mathematics, Huizhou University, Huizhou, Guangdong 516007, PR China
a r t i c l e
i n f o
Article history: Received 21 November 2017 Revised 14 March 2018 Accepted 13 July 2018 Available online xxx MSC: 65M06 65M20
a b s t r a c t In this work, a 6th-order central-upwind weighted essentially non-oscillatory scheme is devised by additionally introducing a large stencil. The constrained condition of choosing this large stencil is that it should contain upwind information to raise the numerical stability of schemes. The implementation of scheme presented here is as compact as the classical weighted essentially non-oscillatory scheme. Several numerical examples show the robust and high-resolution performances of the scheme which is comparable with the recently developed 6th-order weighted essentially non-oscillatory scheme. © 2018 Elsevier Ltd. All rights reserved.
Keywords: WENO scheme Spurious oscillation Hyperbolic conservation laws
1. Introduction The distinguishing characteristic of nonlinear hyperbolic conservation laws is that their solutions may contain discontinuities such as shock waves even for smooth initial conditions. Due to this characteristic, the high-order numerical methods constructed by traditional ways could generate spurious oscillations around the discontinuous regions. Thus, to cure these oscillations, some extra techniques are needed to devise the high-order numerical methods. Several high-order numerical methods, e.g. essentially nonoscillatory schemes (ENO) [1,2], weighted ENO schemes (WENO) [3,4], discontinuous Galerkin finite element methods (DG FEM) [5,6] and so on, are capable of producing satisfactory numerical results. The classical WENO schemes [4], as a kind of the most popular schemes, have attracted a lot of attention due to the feature of easy implementation and high resolution. Several improvements have been developed to decrease the numerical dissipation of the classical WENO schemes. One kind of improvement, e.g. [7,8] and so on, is to modify the weights and distribute a little more weights to the less smooth stencils. Another improvement, e.g. [9–13], is to couple the low-dissipation scheme which is applied on the smooth region and the WENO scheme which is applied on the non-smooth region. Recently, based on the often used 5th-order WENO scheme, several 6th-order central-upwind WENO schemes [14–16] have been
developed by adding one more stencil. The idea to devise the 6thorder WENO originates from [14] in which the schemes are designed to maximize order of accuracy and bandwidth, and meanwhile minimize dissipation. By including an additional downwind stencil, the overall central schemes WENO-SYMOO and WENOSYMBO in [14] are constructed on the basis of optimal order of accuracy and optimal bandwidth efficiency, respectively. However, the order of WENO-SYMBO scheme degenerates even in smooth regions, and the WENO-SYMOO scheme introduces numerical instabilities near the contact surfaces even when only moderate discontinuities are involved. In [15], with the same stencils as [14], the authors define the smoothness indicator of downwind stencil using all six grid values. To achieve the optimal order of accuracy, a reference smoothness indicator τ is devised which is similar to that first presented in [8]. In addition, it is found that this scheme needs extra artificial dissipation to maintain the numerical stability. In [16], another 6th-order scheme is developed which is analogous with that in [15] but with different smoothness indicator for the newly-introduced downwind stencil and no additional artificial dissipation. We will choose this scheme as comparison since it is recently developed and easily implemented. In numerical simulations, it is found that this scheme cannot tolerate large CFL number. That is, with large CFL number, the scheme generates visible oscillations around discontinuities. Instead of disappearing with increasing cell number, the oscillations show the apparent trend of growth.
E-mail address: [email protected] https://doi.org/10.1016/j.compfluid.2018.07.008 0045-7930/© 2018 Elsevier Ltd. All rights reserved.
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 2
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
To alleviate (or eliminate) the effect of large CFL number on the stability of scheme, we introduce a different stencil into the classical 5th-order WENO scheme [4]. Since the introduction of downwind stencil may lead to certain numerical instability especially for large CFL number, the stencil introduced should contain upwind information. Using the smoothness indicator consisted of all six grid values to measure the smoothness of the downwind stencil consisted of three grid values, as did in [14–16], is rather abrupt. The smoothness indicator of newly-introduced stencil here is directly measured by the formulation presented by Jiang and Shu in [4]. In general, to make the scheme arrive the 6th-order of accuracy and each stencil contain upwind information, this new stencil needs to be larger than original stencils in theory. The nonlinear weights we get cannot satisfy the sufficient condition of retaining the optimal order, by Taylor analysis, and in fact only the 5thorder of accuracy can be obtained. But, this 5th-order scheme still has the smaller errors and higher resolution when compared with the classical 5th-order WENO scheme [4]. In order to cure the degeneration in accuracy, there are two reliable methods generally: embodying a mapping function [7,17,18] into the construction of weights or devising a reference smoothness indicator [8]. In addition, the implementation of schemes presented here is compact like the schemes in [7] or [8] and no artificial parameter is used to control the resolution of solutions. The organization of this paper is as follows. In Section 2, the recently developed 6th-order WENO scheme [16] is reviewed and we also investigate the numerical oscillations generated by this scheme under large CFL number. In Section 3, we present the new 6th-order WENO schemes. The mapping function in [18] is used in order to achieve the optimal accuracy. Also, as another way, a reference smoothness indicator is also devised to recover the optimal accuracy. Finally, in Section 4, several common examples are used to demonstrate the numerical performances of the new 6th-order WENO schemes. 2. The 6th-order WENO-CU6 scheme
In this section, we briefly review the recently developed 6thorder central-upwind WENO scheme [16] for solving hyperbolic systems of conservation laws
( u1 , · · ·
a ≤ x ≤ b,
t ≥ 0,
, um )T
(1)
where u = and g(u ) = system (1) is hyperbolic if Jacobian matrix eigenvalues,
(g1 (u ), · · · gm (u ))T . dg du
We say the has real and distinct
k = 1, · · · , m
(3)
where the positive flux f jk,+ = f jk,−
=
1 k 2(fj
− α k vkj ). The scalar
1 k 2(fj α k is
+ α k vkj ) and negative flux
for all values of u. First, the computational domain [a, b] is divided into N cells a = x1/2 < x3/2 < · · · < xN−1/2 < xN+1/2 = b. Denote the interval [xi−1/2 , xi+1/2 ] by Ii and the length of the cell Ii by xi . For simplicity, we will generally assume a uniform grid, but this is not required. Solving the conservation laws (1) with a conservative difference approximation to the spatial derivative
(2)
where ui (t) is the numerical approximation to the column vector u(xi , t). what to do in the next is how to compute the numerical flux gˆi±1/2 properly. Following the procedure 2.10 in [19], we use the way of characteristic-wise decomposition to project the relevant vectors uj and g(uj ) into the local characteristic fields, where j is in a neighborhood of i. Denote by vj and fj the obtained vectors in the local characteristic fields. To avoid the presentation of
taken as α = max λk (ui ) 1≤i≤N
over the global range of u. The WENO-CU6 reconstruction is used ± to computed each component fˆik, at xi+1/2 and then we obtain +1/2 k, + k, − k fˆ = fˆ + fˆ . Finally, transforming the fˆi+1/2 in characteri+1/2
i+1/2
i+1/2
istic fields back to the physical space, we obtain the numerical flux gˆi+1/2 . In general, the difference between the different WENO schemes lies on the reconstruction of fˆk,± . Here we only give the i+1/2
+ − procedure of obtaining the positive numerical flux fˆik, , for fˆik, +1/2 +1/2 a similar way can be used. For simplicity, we will omit the super± script k in fˆik, in the following procedure. +1/2 In [16], the positive numerical flux fˆi++1/2 can be achieved by the following three steps: 1. Compute the smoothness indicators of four substencils S0 = {Ii−2 , Ii−1 , Ii }, S1 = {Ii−1 , Ii , Ii+1 }, S2 = {Ii , Ii+1 , Ii+2 }, S3 = {Ii+1 , Ii+2 , Ii+3 } by
βr =
2 l=1
xi+1/2 xi−1/2
x2l−1 ( p(rl ) (x ))2 dx,
r = 0, 1, 2
(4)
and
β3 =
5 l=1
xi+1/2 xi−1/2
x2l−1 ( p(3l ) (x ))2 dx
(5)
where pr (x), r = 0, 1, 2, are the 2nd-degree interpolation polynomials over the interval [xi+r−5/2 , xi+r+1/2 ] and p3 (x) is the 5th-degree interpolation polynomial over the whole interval [xi−5/2 , xi+7/2 ]. 2. Compute the nonlinear weights wr by the following formulas
ar , ar = dr C + a0 + a1 + a2 + a3 1 = β3 − (β0 + 4β1 + β2 ) 6
τ
τ
βr + ε
, (6)
and where d0 = 1/20, d1 = 9/20, d2 = 9/20 and d3 = 1/20 are the optimal weights and ε is a parameter used to avoid a division by zero in the denominator and set to be ε = 10−40 in all numerical examples. 3. Compute the convex combination of four substencils,
fˆi++1/2 = w0 fˆi++1,0/2 + w1 fˆi++1,1/2 + w2 fˆi++1,2/2 + w3 fˆi++1,3/2
λ1 (u ) < λ2 (u ) < · · · < λm (u ),
1 dui (t ) =− gˆi+1/2 − gˆi−1/2 , dt x
f jk = f jk,+ + f jk,− ,
wr =
2.1. Review the construction of the WENO-CU6 scheme
ut + g(u )x = 0,
entropy violating solutions, for the each component of characteristic variables, the simple scalar Lax–Friedrichs splitting is used to divide the each component f jk of fj into
(7)
,r fˆi++1 /2
where are the 2nd-order numerical fluxes obtained on stencil Sr . Substituting fˆi++1/2 and fˆi−+1/2 obtained by the above procedure into (3), we can get the numerical flux fˆi+1/2 at the boundary of cell xi+1/2 . Similarly, we can get the numerical flux fˆi−1/2 at xi−1/2 and finally form the scheme (2). The discretization of the ODEs (2) in time obtained from spatial discretization using the methods of lines approach are solved by the third-order TVD Runge–Kutta method [20,21]
u(1) = un + tL(un ), 3 1 u ( 2 ) = un + u ( 1 ) + 4 4 1 n 2 (2 ) n+1 u = u + u + 3 3
1 tL u(1) , 4 2 tL u(2) , 3
where L is the spatial discretization operator and at time tn .
(8) un
is the solution
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
3
Fig. 1. The numerical solutions are computed by the WENO-CU6 for the initial value problem (27) and (29) with CFL = 0.6 and C = 20. N denotes the cell number.
2.2. The problem appeared in WENO-CU6 scheme The WENO-CU6 scheme is developed by adding a downwind stencil into the classical 5th-order WENO scheme [4]. To alleviate the instability brought by downwind information, the smoothness indicator of this stencil is defined on whole stencil S = {Ii−2 , Ii−1 , Ii , Ii+1 , Ii+2 , Ii+3 }. However, the scheme only achieves the 4th-order accuracy if we directly define the nonlinear weights as that in [4]
wr =
ar , a0 + a1 + a2 + a3
ar =
dr , (βr + ε )2
ously, in the case of C = 1, the scheme admits CFL number up to 0.5 which is relaxed when compared with the case of C = 20. In addition, we find that from Fig. 2, with the same CFL number, e.g. CFL = 0.6, the solution with larger parameter C has larger oscillation. For much clearer to see this trend, the Fig. 3 shows the oscillations behind the discontinuity with the fixed CFL=0.6 but varying the parameter C. And now, we can easily find that the oscillations grow larger with the choice of larger parameters C. 3. The construction of new 6th-order WENO schemes
r = 0, 1, 2, 3.
(9)
This claim can be easily confirmed by symbolic operations or numerical tests. It is noted that, although the scheme with nonlinear weights (9) cannot achieves optimal order, it can tolerate a reasonable CFL number as that of the classical 5th-order WENO scheme [4]. To cure the degeneration of order of accuracy, the authors apply the technique proposed in [8] to define the weights as in (6). And a new parameter C = 20 is introduced in (6) for increasing the contribution of optimal weights when the smoothness indicators have comparable magnitudes. But, for the WENO-CU6 scheme, the construction of weight in (6) and the introduction of parameter C = 20 both result in the weak tolerance of large CFL number. That is, the scheme will generate oscillations around the discontinuities if one chooses a large CFL number. Considering the simple example (27) and (29), we first set C = 20 and t = CFLx where CFL = 0.6. From the Fig. 1, one can find there are obvious oscillations behind the discontinuities at x = −0.4 and x = −0.2. The close-up view around x = −0.2 is shown in right subgraph of Fig. 1. In addition, we find the oscillations do not yet disappear as increasing cell number to N = 3200, but present the trend of growth. Next, we show the influence of different CFL numbers on the WENO-CU6 scheme with C = 20 and N = 3200. The left subgraph in Fig. 2 shows the numerical discontinuity around x = −0.2 with CFL =0.3, 0.4, 0.5, 0.6, 0.7. The oscillations appear visibly when CFL = 0.5 and increase gradually with increasing CFL number. Hence, at least for this example, a safe choice of CFL number is 0.4 (or even smaller 0.3). As for the parameter C, it also affects the choice of CFL number. Although the larger parameter C generally produces the slightly higher-resolution numerical solutions, it often imposes more requirement on CFL number. In contrast, the lesser parameter C makes the choice of CFL number more relaxed. For example, decreasing the parameter C to 1, it will be found that the scheme can tolerate the larger CFL number. The right subgraph in Fig. 2 are computed with C = 1 and N = 3200 which has same cell number with left figure and the only difference is the choice of C. Obvi-
In this section, we present the procedure of constructing the new 6th-order WENO scheme. As other 6th-order WENO schemes shown in [14–16], we also add an extra stencil S3 shown in left subgraph of Fig. 4 (or S0 in right subgraph). The difference is that the new stencil introduced here contains upwind information which may aid in numerical stability. From Fig. 4, to construct the positive numerical flux fˆi++1/2 and the negative one fˆ− with the 6th order of accuracy, the steni+1/2
cils
S0+ = {Ii−2 , Ii−1 , Ii }, S1+ = {Ii−1 , Ii , Ii+1 }, S2+ = {Ii , Ii+1 , Ii+2 }, S3+ = {Ii , Ii+1 , Ii+2 , Ii+3 }. and
S0− = {Ii−2 , Ii−1 , Ii , Ii+1 }, S1− = {Ii−1 , Ii , Ii+1 }, S2− = {Ii , Ii+1 , Ii+2 }, S3− = {Ii+1 , Ii+2 , Ii+3 }.
are used, respectively. Using each stencil Sr± , r = 0, · · · , 3, one gets four low-order numerical fluxes
7 1 11 + fˆi++1,0/2 = fi+−2 − fi+−1 + f , 3 6 6 i 5 − 1 − 13 − 1 − 0 fˆi−, = f − f + f + fi+1 , +1/2 12 i−2 12 i−1 12 i 4 ˆf +,1 = − 1 f + + 5 f + + 1 f + , i+1/2 6 i−1 6 i 3 i+1 1 5 1 −, 1 fˆi+1/2 = − fi−−1 + fi− + fi−+1 , 6 6 3 ˆf +,2 = 1 f + + 5 f + − 1 f + , i+1/2 3 i 6 i+1 6 i+2 1 1 5 2 fˆi−, = fi− + fi−+1 − fi−+2 , +1/2 3 6 6 5 + 1 13 + 1 + fˆi++1,3/2 = fi+ + f − f + f , 4 12 i+1 12 i+2 12 i+3 7 11 − 1 3 fˆi−, = f − f− + f− , +1/2 6 i+1 6 i+2 3 i+3
(10)
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 4
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
Fig. 2. The left is computed with C = 20 and right with C = 1.
Inserting (10) into (11), through simple computing one gets the linear optimal weights
d0+ =
1 9 3 1 , d1+ = , d2+ = , d3+ = , 20 20 10 5
and
d0− =
1 3 9 1 , d1− = , d2− = , d3− = . 5 10 20 20
Remark 1. In above procedure, taking the positive flux fˆi++1/2 for
example, we choose S3+ = {Ii , Ii+1 , Ii+2 , Ii+3 } as the fourth stencil. Actually, the stencil S3+ = {Ii−1 , Ii , Ii+1 , Ii+2 , Ii+3 } consisted of five points can also be chosen as the candidate. Under this choice, the overall scheme is also stable numerically. However, through comparison, it is found that the first choice has higher resolution and lower computational cost. So we only discuss the scheme formed by the first choice in this paper.
Fig. 3. The numerical solutions are computed by WENO-CU6 scheme with N = 1600. We fix CFL = 0.6 but vary parameter C.
Following the idea of the classical WENO scheme, a nonlinear convex combination of low-order numerical flux fˆi±,r similar to +1/2 (7) is needed to control the spurious oscillations around the discontinuities ± ˆ±,1 ± ˆ±,2 ± ˆ±,3 ˆ±,0 fˆi±+1/2 = w± 0 f i+1/2 + w1 f i+1/2 + w2 f i+1/2 + w3 f i+1/2 ,
,3 0 and the numerical flux fˆi++1 computed on S3+ and fˆi−, computed /2 +1/2
on S0− have 4th-order accuracy. Through the linear convex combination of four low-order stencils S0 , S1 , S2 and S3 , one wishes to get the 6th-order central interpolation on 6-point stencil S
1 + fi−2 − 8 fi+−1 + 37 fi+ + 37 fi++1 − 8 fi++2 + fi++3 , 60 d− fˆ−,0 + d− fˆ−,1 + d− fˆ−,2 + d− fˆ−,3 =
1
i+1/2
2
i+1/2
3
a± 0
+
a± 1
a± r , + a± + a± 2 3
a± r =
(β
dr± , + ε )2
± r
r = 0, 1, 2, 3. (13)
The smoothness indicators are computed by (2.61) in [19]
β2+
i+1/2
1 − fi−2 − 8 fi−−1 + 37 fi− + 37 fi−+1 − 8 fi−+2 + fi−+3 . 60
w± r =
β1+
=
i+1/2
where the nonlinear weights w± r can be obtained as that in [4]
2 1 + 2 13 + f − 2 fi+−1 + fi+ + f − 4 fi+−1 + 3 fi+ , 12 i−2 4 i−2 2 1 + 2 13 + = fi−1 − 2 fi+ + fi++1 + fi−1 − fi++1 , 12 4 2 13 + 1 + + + 2 = f − 2 fi+1 + fi+2 + 3 fi − 4 fi++1 + fi++2 , 12 i 4 2 1 = 11 fi+ − 18 fi++1 + 9 fi++2 − 2 fi++3 36
β0+ =
d0+ fˆi++1,0/2 + d1+ fˆi++1,1/2 + d2+ fˆi++1,2/2 + d3+ fˆi++1,3/2
0
(12)
(11)
β3+
Fig. 4. Candidate stencils Sr for the numerical flux fˆi++1/2 (left) and fˆi−+1/2 (right).
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
2 13 + 2 fi − 5 fi++1 + 4 fi++2 − fi++3 12 2 781 + + fi − 3 fi++1 + 3 fi++2 − fi++3 , 720 +
(14)
β0− =
β
− 1
β2− β3−
w r − d r = O x 4 , r = 0 , 1 , 2 , w 3 − d 3 = O x 3 .
(16)
In order to more easily implement, the sufficient conditions (16) can be extended to the more compact formulas
w r − d r = O x 4 , r = 0 , · · · , 3 .
(15)
where, for β3+ and β0− , we have used the more compact explicit formulas with the sum of squares in [18]. Finally, by similar way we can obtain the negative numerical flux fˆi−+1/2 and finally the numerical flux fˆi+1/2 at xi+1/2 is
fˆi+1/2 = fˆi++1/2 + fˆi−+1/2 . Remark 2. Despite the fact that the formulas of linear weights dr , low-order numerical fluxes fi±,r and smoothness indicators are all +1/2 not symmetrical, the overall scheme to compute the numerical flux fˆi+1/2 at cell boundary xi+1/2 is symmetrical at smooth region. At non-smooth region, the scheme is upwind. Hence this kind of 6thorder scheme is also called central-upwind scheme. Up to present, we have presented the numerical flux fˆi+1/2 at xi+1/2 . One can find that the construction of this scheme is as compact as the classical WENO scheme in [4,22]. It is supposed to obtain the 6th-order approximation
fˆi+1/2 − fˆi−1/2 = fx (u(xi , t )) + O (x6 ) x since we have constructed fˆi+1/2 and fˆi−1/2 with 7 flux values
f (ui−3 ), f (ui−2 ), f (ui−1 ), f (ui ), f (ui+1 ), f (ui+2 ), f (ui+3 ). Unfortunately, the mixed smoothness indicators of three 2nddegree and one 3rd-degree polynomials can only guarantee
generates much higher-resolution solutions and presents lower numerical errors when compared with the classical 5th-order WENO scheme. As in [7], we use similar way to derive the sufficient conditions for making the scheme achieve the 6th order of accuracy
and
2 1 − 2 fi−2 − 9 fi−−1 + 18 fi− − 11 fi−+1 36 2 13 − + f − 4 fi−−1 + 5 fi− − 2 fi−+1 12 i−2 2 781 − + fi−2 − 3 fi−−1 + 3 fi− − fi−+1 , 720 2 1 − 2 13 − = f − 2 fi− + fi−+1 + f − 4 fi− + 3 fi−+1 , 12 i−1 4 i−1 2 1 − 2 13 − = fi − 2 fi−+1 + fi−+2 + fi − fi−+2 , 12 4 2 13 − 1 − − − 2 = f − 2 fi+2 + fi+3 + 3 fi+1 − 4 fi−+2 + fi−+3 , 12 i+1 4
5
2
w r = d r + O x , r = 0 , · · · 3 , which make the scheme achieve the 5th-order of accuracy at most when the solution is smooth and meanwhile no critical point appeared. Detailedly, by Taylor expansion, we indeed only get
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 5 . x At critical points, the order of accuracy will degenerate seriously. For example, for f = 0, we get
wr = dr + O (x ), r = 0, · · · 3, and
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 3 . x
(17)
So, to recover the optimal order of accuracy, we only need to make the nonlinear weights approximate the linear weights with O (x4 ). In following parts of this section, we will present two methods to recover the optimal order of accuracy of the scheme: one is to apply a mapping function, as in [7], on the nonlinear weights computed above, and the other is to construct a reference smoothness indicator as in [8]. The common purpose of these two remedy methods is to make the nonlinear weights converge to the linear weights fast enough as x → 0. 3.1. The mapping function The mapping function devised in [7]
g¯ (wr ) = dr +
dr2
( wr − dr )3 + ( 1 − 2dr )wr
(18)
is enough to make the nonlinear weights satisfy the sufficient condition (17). At the smooth region with no critical point, since the nonlinear weights achieve wr − dr = O (x2 ), r = 0, · · · , 3, applying the mapping function (18) on nonlinear weights gives
w r − d r = O x 6 , r = 0 , · · · , 3 .
(19)
However, when at critical point f = 0 which leads to wr = dr + O (x ), this mapping function only ensures wr = dr + O (x3 ) which is not enough to guarantee the optimal order. Actually, in this case, the scheme only achieves 5th-order accuracy with the mapping function (18). Hence, the generalized mapping function in [18] is preferred to be used
g( w r ) = d r +
(wr − dr )7 C0 + C1 wr + C2 w2r + C3 w3r
,
(20)
where C0 = dr6 , C1 = −7dr5 , C2 = 21dr4 , C3 = (1 − dr )6 − (C0 + C1 + C2 ). This function apparently recovers the optimal order of 6thorder scheme even at critical points f = 0. Besides, when compared with the mapped WENO in [7], it allocates much smaller weights to discontinuous stencils which could further compress the spurious oscillations. After obtaining the mapping function (20), the improved nonlinear weights which satisfy the sufficient condition (17) are computed by
aM r wM r = k=3 k=0
aM k
,
aM r = g( w r ) ,
r = 0, · · · , 3.
Hence, for the numerical tests in Section 4, we always choose the mapping function (20) unless otherwise stated.
That is, this “6th-order” scheme only achieves 3rd-order of accuracy around critical points with f = 0.
3.2. The reference smoothness indicator
Remark 3. Although this “6th-order” scheme cannot yet achieve the optimal accuracy, by numerical experiment, we find that it still
In this section, we plan to present the other method appeared in [8] to recover the optimal order of 6th-order scheme. The idea
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 6
of this correction is to construct a high-order reference smoothness indicator. Expanding the smoothness indicators (14), we get
13 β0 = fi2 x2 + f 2 − 12 i x 5 + O x 6 , 13 β1 = fi2 x2 + f 2 + 12 i 13 β2 = fi2 x2 + f 2 − 12 i x 5 + O x 6 , β3
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
2 f f x 4 + 3 i i
1 2
fi fi(4) −
13 f f 6 i i
1 f i f i x 4 + O x 6 , 3 1 2 13 fi fi x4 + − fi fi(4) + fi fi 3
2
6
(21)
13 2 4 1 (4) 5 = fi2 x2 + f x + f i f i x + O x 6 . 12 i 2
(21)
Since the reference smoothness indicator here is preferred to contain all 6 points, the highest-order one we can obtain is the following linear combination
1 τ = β3 − (β0 + 4β1 + β2 ) = O x5 .
(22)
6
Remark 4. Actually, there exists another combination τ = ((β1 − β0 ) − 32 (β3 − β2 ) which also arrives O (x5 ). However, the numerical results with this reference smoothness indicator are worse than that with (22). So, in this work, we only consider the reference smoothness indicator (22). Finally, the nonlinear weights can be computed similarly by
aZr
wZr = k=3 k=0
aZk
,
aZr = dr 1 +
p
τ
,
βr +
r = 0, · · · , 3, (23)
where the parameter p is used to control the order of accuracy and dissipation. For the case of p = 1, through simple computation, it seems to get that
wZr = dr + O x3 .
(24)
Obviously, this “doesn’t” satisfy the sufficient condition (17). But, it was not as straightforward as we were led to believe. By detailed Taylor expansion, one can get
aZr = dr 1 + Ax3 + Bx4 + O x5 ,
Remark 5. It is noted that this 6th-order scheme has less computational cost compared with the scheme in [16]. The difference of computational cost between these two schemes mainly lies on the computation of smoothness indicator of new introduced stencil. Comparing the smoothness indicator (14) (or (15)) and the formula (25) in [16] (which we used is the more compact the sum of squared), the latter has two sum of squared more than the former. Up to now, two methods have been used to recover the optimal order of accuracy. The use of mapping function can achieve the optimal order efficiently even around the critical points. The flaw of this method is the raise of computational cost compared to the classical WENO. The use of reference smoothness indicator with the parameter p = 1 also achieves the optimal order efficiently but not including the critical points. A compromised way to recover the 6th-order accuracy at critical points is to enlarge the parameter p. Nevertheless, this will increase the dissipation of the scheme and the scheme presents the low resolution around the discontinuities which is confirmed in [8]. 4. Numerical results In this section, several 1D and 2D numerical examples are used to test the 6th-order WENO schemes. The 6th-order WENO scheme with the mapping function is denoted by WENO-M6 and with reference smoothness indicator is denoted by WENO-Z6. For comparison, we choose the recently developed 6th-order scheme in [16] denoted by WENO-CU6 and the parameter C = 20. Although the similar parameter C can also be introduced in the WENO-Z6 scheme to increase the resolution, we abandon it for the compaction of scheme and expect to make the numerical solutions look monotone. Overall, the WENO-CU6 and WENO-Z6 present similar errors and accuracy if the same parameters C and p are set in schemes. 4.1. The linear advection equation
f (4 ) 2f
In this section, we consider the simple case
ut + ux = 0, x ∈ [−1, 1], t ≥ 0,
and
u(x, 0 ) = sin
Finally, the nonlinear weights in (23) arrive
d r 1 + A x 3 + B x 4 + O x 5 1 + A
x3
+ B
x4
+O
x5
πx −
sinπ x
(28)
π
and
5
= d r + O x .
(26)
Hence, the nonlinear weights constructed in (23) actually approach the linear weights with O (x5 ) and conform to the requirement (17) of retaining the 6th-order accuracy at smooth region with no critical point. For the 1st-order critical point, by symbolic operation the nonlinear weights achieve
wZr = dr + O x3
(27)
with the following two initial conditions to test the schemes,
2 4 (5) 91 3 (4) 481 3 (3) 2 f f − f f f + f f 36 720 B= 3 . 5 f
wZr =
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 6 . x
(25)
where
A=
To alleviate the degeneration of accuracy at 1st-order critical point, the favourite way is to enlarge the parameter p in (23). For instance, in the case of p = 2, around the smooth region with the 1st-order critical point, the scheme recovers the optimal order
⎧ ⎪ ⎪(G(x, β , z − δ ) + G(x, β , z + δ ) + 4G(x, β , z ) )/6 ⎪ −0.8 ≤ x ≤ −0.6, ⎪ ⎪ ⎪ 1 − 0.4 ≤ x ≤ −0.2, ⎨ 1 − |10x − 1|0.0 ≤ x ≤ 0.2, u(x, 0 ) = ⎪ ⎪ (F (x, α , a − δ ) + F (x, α , a + δ ) + 4F (x, α , a ) )/6 ⎪ ⎪ ⎪ ⎪ ⎩ 0.4 ≤ x ≤ 0.6, 0
otherwise,
and the order of accuracy of scheme is
fˆi+1/2 − fˆi−1/2 = f x ( u ( x i , t ) ) + O x 5 . x
(29) where
G(x, β , z ) = exp(−β (x − z )2 ),
F (x, α , a ) =
max(1 − (α (x − a ))2 , 0 ), and the parameters are given by
a = 0.5, z = −0.7,
δ = 0.005, α = 10, β =
ln 2 . 36δ 2
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
7
Fig. 5. The l1 and l∞ errors are computed by the 6th-order linear, WENO-CU6, WENO-M6 and WENO-Z6 schemes, respectively, for the initial value problem (27) and (28). Table 1 The l1 errors and rate of convergence are shown around the square wave in the problem (27) and (29) with CFL = 0.3 and the output time t = 2.
x 2.00e−02 1.00e−02 5.00e−03 2.50e−03
WENO-CU6
WENO-M6
WENO-Z6
l1 error
r
l1 error
r
l1 error
r
4.07342e−02 2.21530e−02 1.23156e−02 6.85770e−03
0.8787 0.8470 0.8447
4.83906e−02 2.54545e−02 1.39010e−02 7.58246e−03
0.9268 0.8727 0.8745
4.53699e−02 2.48965e−02 1.37859e−02 7.62225e−03
0.8658 0.8527 0.8549
The first condition with smooth solution is used to test the order of accuracy and the second with discontinuous and high-frequency waves is used to test the essentially non-oscillatory property of schemes. All three methods are used to compute the initial problem (27) and (28). The time step is t = x2 and ε = 10−40 . The l1 and l∞ errors of schemes are shown in Fig. 5 and all the WENO schemes present the 6th order of accuracy and show the almost same errors as the linear central 6th-order scheme. Here, the parameter p in (23) is chosen to be 2 for recovering the order at critical points (for all other examples p = 1). It is noted that the error of WENO-Z6 will be enlarged slightly if p = 1 for this smooth example. Of course, the introduction of C = 20 in WENO-Z6 will decrease the error to that of WENO-CU6 scheme even p = 1. Besides the accuracy of 6th-order schemes for the smooth solution, the rate of convergence around the discontinuity is also tested here. For nonlinear shock waves, the high-resolution schemes generally produce a 1st-order of accuracy. However, we are interested in the case of linear (or linearly degenerate) discontinuities which will dominate the global rate of convergence. As the analysis in [23], the 6th-order WENO schemes here should produce the expected rates of convergence O (x6/7 ). We choose the problem (27) and (29) which contains Gaussian, square, triangle and semiellipse waves to test the rate of convergence around the square (which is linearly discontinuous). To eliminate the influence of other three waves, we only present the numerical errors around [−0.5, −0.1]. Table 1 shows the rate of convergence of the 6thorder schemes. For the problem (27) and (29), the initial condition contains Gaussian, square, triangle and semi-ellipse waves. Fig. 6 shows the solutions computed by three methods with cell number N = 200 and time step t = 0.3x. They all present high-resolution solutions around Gaussian, triangle and semi-ellipse waves except for the square. But, from the right amplified subgraph in Fig. 6, the WENO-CU6 produces much better than other two schemes. This
owes to the introduction of parameter C = 20 into the WENO-CU6 scheme and the WENO-Z6 will generate the similar solution with WENO-CU6 if we introduce the parameter C = 20 into WENO-Z6. But, as shown in Section 2.2, it has been found that, in the WENOCU6 scheme, the introduction of C = 20 partly leads to the poor tolerance of large CFL number. That is the WENO-CU6 scheme easily produces visible oscillations behind discontinuity if we increase the CFL number. And the magnitude of oscillation do not disappear (or even alleviate) with increasing the cell number. Through this example, we wonder how the CFL number affects the numerical stability of WENO-M6 and WENO-Z6 schemes. Fig. 7 shows the results of these two schemes at x = −0.2 computed with N = 3200 and several different CFL numbers. In the left subgraph, the WENO-M6 scheme is almost unaffected by CFL number and does not present any oscillation even for CFL = 0.8. However, from the right subgraph, we can find that no oscillation is observed with CFL number up to 0.6 and there are clear oscillation when CFL ≥ 0.7 (here we do not consider the situation when CFL number between 0.6 and 0.7). It is noted that, under the case of C = 1, the WENO-CU6 conserves no oscillation with CFL number up to 0.5. Hence the WENO-Z6 has stronger tolerance on the CFL number than WENO-CU6 scheme at least for this example. And in Section 4.2, we will find that similar conclusion is still hold even for nonlinear Euler systems. In addition, comparing the results between the right subgraph of Fig. 2 (where C = 1 in WENOCU6) and the right subgraph of Fig. 7, the magnitude of oscillation generated by WENO-CU6 is slight larger than WENO-Z6 under the same CFL number.
4.2. 1D Euler equations In this example, we consider 1D Euler equations since one of the main application areas of weighted essentially non-oscillatory
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 8
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
Fig. 6. The numerical solutions are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes, respectively, for the initial value problem (27) and (29) with N = 200 and CFL = 0.3.
Fig. 7. The numerical solutions are computed by the WENO-M6 (left) and WENO-Z6 (right) for the initial value problem (27) and (29) with N = 3200.
schemes is compressible gas dynamics,
ρ ρu
+
E
t
ρu ρ u2 + p ( E + p )u
= 0,
(30)
x
where ρ , u, p, E are density, velocity, pressure and total energy, respectively. The system of equations is closed by the equation of state for an ideal polytropic gas:
E=
γ
p 1 + ρ u2 , −1 2
U (l, i ) = U (l, N ), l = 1, 2, 3; i = N + 1, N + 2, N + 3.
1. (Sod problem)
(ρ , u, p) =
( 1, 0, 1 ), (0.125, 0, 0.1 ),
if − 5 ≤ x < 0, if 0 ≤ x ≤ 5,
(31)
2. (Shock entropy wave interaction problem)
(ρ , u, p) =
⎧
√ ⎨ 27 4 35 31 ⎩
, , , 7 9 3 (1 + 0.2 sin(5x ), 0, 1 ),
if − 5 ≤ x < −4, if − 4 ≤ x ≤ 5, (32)
3. (Blast wave interaction problem)
(1, 0, 10 0 0 ), (ρ , u, p) = (1, 0, 0.01 ), (1, 0, 100 ),
if − 5 ≤ x < −4, if − 4 ≤ x < 4, if 4 ≤ x ≤ 5.
U (l, i ) = U (l, 1 ), l = 1, 2, 3; i = −2, −1, 0, and
where the ratio of specific heats γ = 1.4. The following three initial conditions are considered
The shock tube problem (30) and (31) is called Sod problem [24]. This problem contains a rarefaction wave, a contact discontinuity and a shock. At the output time t = 2, Since waves created by the discontinuity at x = 0 do not reach the boundaries, the zero gradient boundary conditions are applied on the left and right boundaries. Specifically, denote U (l, i ) = (ρ (i ), ρ (i )u(i ), E (i ))T , l = 1, 2, 3, the practical implementations of boundary condition are
(33)
We first test the rate of convergence around the contact discontinuity. Also, to avoid the influence of other waves, we choose to only compute the l1 errors in interval [1,2.75]. Table 2 shows the numerical errors of 6th-order schemes and all the rates of convergence are approximately equal to the expected order O (x6/7 ). In addition, Fig. 8 shows the density distributions of the 6thorder schemes and they all achieve high resolution. These numerical density is calculated with moderate CFL = 0.6. However, the overshoot will appear clearly around discontinuities if CFL number is increased to 0.8. Fig. 9 presents the results of these methods under CFL = 0.8 with cell numbers N = 20 0, 40 0, 80 0, 160 0. For the WENO-CU6 scheme, there exist obvious oscillations behind the contact discontinuity and the magnitude of oscillations has no trend of disappearance with increasing the cell number. Another observation from Fig. 9 is that there are little overshoot behind
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
9
Table 2 The l1 errors and rate of convergence are shown around the square wave in the problem (30) and (31) with CFL = 0.3 and the output time t = 2.
x 1.00e−02 5.00e−03 2.50e−03 1.25e−03
WENO-CU6
WENO-M6
WENO-Z6
l1 error
r
l1 error
r
l1 error
r
1.81614e−03 1.03725e−03 5.73400e−04 3.19969e−04
0.8081 0.8551 0.8416
2.06224e−03 1.17450e−03 6.35782e−04 3.45954e−04
0.8122 0.8854 0.8631
2.03386e−03 1.16065e−03 6.35977e−04 3.48802e−04
0.8096 0.8679 0.8666
Fig. 8. The density profiles of problem (30) and (31) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 200 and CFL = 0.6.
Fig. 9. The density profiles of problem (30) and (31) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 200, 400, 800, 1600 and CFL = 0.8.
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
ARTICLE IN PRESS
JID: CAF 10
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
Fig. 10. The density profiles of problem (30) and (32) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 400 and CFL = 0.6.
Fig. 11. The density profiles of problem (30) and (33) are computed by the WENO-CU6, WENO-M6 and WENO-Z6 schemes with N = 400 and CFL = 0.6.
the shock and they disappear with increasing the cell number. For the WENO-Z6 scheme, there exist only little overshoot behind the shock with low-resolution cell number which disappears as the WENO-CU6 scheme with increasing the cell number. The WENOM6 scheme again behaves robustly and produces look-monotone solutions. Hence, as in the last example (27) and (29), the ability of tolerating the large CFL number of WENO-CU6 is less than that of WENO-M6 and WENO-Z6. The shock entropy-wave interaction problem (30) and (32) in [25] simulates the interaction of a right-moving Mach 3 shock with a wavelike perturbed density whose magnitude is much smaller than the shock. Also, the zero gradient boundary conditions are applied on the left and right boundaries. As a result, a flow field of compressed and amplified wave trails is created right behind the shock. The “exact” solution is computed by the 5th-order WENOJS scheme with N = 6400. From Fig. 10, we can find that all the schemes generate high-resolution results. By amplifying the part of entropy waves, the WENO-CU6 scheme shows slightly higher resolution than the WENO-M6 and WENO-Z6 schemes. The last example is the problem (30) and (33) of two interactive blast waves [26] with reflection boundary conditions. the practical implementations of numerical boundary conditions are
U (l, i ) = (−1 )
l+1
U (l, 1 − i ), l = 1, 2, 3; i = −2, −1, 0
for left boundary and
U (l, i ) = (−1 )l+1U (l, 2N + 1 − i ), l = 1, 2, 3; i = N + 1, N + 2, N + 3
for right boundary. The “exact” solution is also computed by the 5th-order WENO-JS scheme with N = 6400. Again, all the schemes perform with high accuracy and the WENO-CU6 scheme shows slightly higher resolution than the WENO-M6 and WENO-Z6 schemes around x = 2.5 in Fig. 11. Remark 6. Through simulating above several linear scalar equations and nonlinear systems, one can find that these three schemes all generate high-accuracy solution. In general, the WENO-CU6 scheme has slightly better numerical solutions than WENO-M6 and WENO-Z6 schemes. The expense of WENO-CU6 is the spurious oscillations behind the linear discontinuity in some case and the oscillations do not seem to disappear with increasing cell number. This weak tolerance is largely due to the introduction of parameter C. 4.3. 2D Euler equations
⎡
In this section, we consider 2D Euler equations,
⎤ ⎡ ⎤ ⎡ ⎤ ρu ρv ρ 2 ⎢ρ u⎥ ⎢ ρ u + p ⎥ ⎢ ρ uv ⎥ ⎣ρv⎦ + ⎣ ρ uv ⎦ + ⎣ ρv2 + p ⎦ = 0, E t ( E + p )u x ( E + p )v y
(34)
where the equation of state is p = (γ − 1 )(E − 12 ρ (u2 + v2 )). Here we apply the 6th-order schemes to the 2D double-Mach shock reflection problem which contains strong shock waves and
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
JID: CAF
ARTICLE IN PRESS
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
11
Fig. 12. The density profiles are computed respectively by the WENO-CU6, WENO-M6 and WENO-Z6 schemes. 30 equally spaced contours for the density from 1.731 to 20.9705 with N = 800 × 200 and CFL = 0.6.
contact discontinuity. This problem was first proposed by Woodward and Colella [26] and had been taken extensively as a test example for high-order schemes. The computational domain is chosen to be [0, 4] × [0, 1] and the reflective wall lies on the bottom of the computational domain for 16 ≤ x ≤ 4. In the beginning, a right-moving Mach 10 shock is located at x = 16 , y = 0 and makes an angle 60° with the x-axis. For the boundary conditions, the exact postshock condition is imposed for bottom boundary from x = 0 to x = 16 , and the reflective boundary condition is imposed for the rest; the flows are imposed on the top boundary such that there is no interaction with the Mach 10 shock; inflow and outflow boundary conditions are set for the left and right boundaries respectively. The unshocked fluid has a density of 1.4, a pressure of 1. Specifically, the boundary conditions can be implemented numerically as follows. Denote by Nx the
grid number in x-axis direct, Ny in y-axis direct and U (l, i, j ) = (ρ (i, j ), ρ (i, j )u(i, j ), ρ (i, j )v(i, j ), E (i, j ))T , where l = 1, · · · , 4, i = 1, · · · , Nx and j = 1, · · · , Ny . At the left boundary,
U (l, i, j ) = U (l, 1, j ), i = −2, −1, 0; at the right boundary,
U (l, i, j ) = U (l, Nx , j ), i = Nx + 1, Nx + 2, Nx + 3; at the bottom boundary, if the cell center x(i) ≤ 1/6 then
√ U (1, i, j ) = 8, U (2, i, j ) = 33 3, U (3, i, j ) = −33, U (4, i, j ) = 563.5, or else
U (1, i, j ) = U (1, i, 1 − j ), U (2, i, j ) = U (2, i, 1 − j ),
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008
JID: CAF 12
ARTICLE IN PRESS
[m5G;July 19, 2018;13:27]
F. Hu / Computers and Fluids 000 (2018) 1–12
U (3, i, j ) = −U (3, i, 1 − j ), U (4, i, j ) = U (4, i, 1 − j ), where j = −2, −1 √, 0; at the top boundary, if the cell center x(i ) ≤ 1/6 + (1 + 20t )/ 3 then
√ U (1, i, j ) = 8, U (2, i, j ) = 33 3, U (3, i, j ) = −33, U (4, i, j ) = 563.5, or else
U (1, i, j ) = 1.4, U (2, i, j ) = 0, U (3, i, j ) = 0, U (4, i, j ) = 2.5, where j = Ny + 1, Ny + 2, Ny + 3. The Fig. 12 shows the density profiles computed by the WENOCU6, WENO-M6 and WENO-Z6 schemes with cell number N = 800 × 200 and output time t = 0.2. The WENO-CU6 scheme gives the better solution than that of WENO-M6 and WENO-Z6 schemes, and the difference between the WENO-CU6 and WENO-M6 is subtle. 5. Conclusions In this work, we present new 6th-order WENO schemes for solving hyperbolic conservation laws. Similar to other 6th-order WENO, we also introduce an additional stencil into the classical 5th-order WENO. The difference is that the new-introduced stencil here contains upwind information and is not completely down-wind stencil. This choice looks more reasonable to control the numerical stability of scheme. A mapping function and reference smoothness indicator are used to make the scheme achieve optimal order, respectively. Then two slightly different schemes, WENO-M6 and WENO-Z6, are developed. We choose the recently developed 6th-order WENO-CU6 scheme as the comparison with the new schemes in this paper. These three schemes all generate the high-resolution solutions and the WENO-CU6 is better relatively under close-up view. Through numerical simulations, we find that the WENO-CU6 scheme has weak tolerance of large CFL number. This weak tolerance is partly due to the introduction of parameter C. The WENO-Z6 also cannot tolerate the large CFL number, but is more robust than the WENO-CU6. The WENO-M6 shows the strong-robust characteristic and is almost unaffected by CFL number. Acknowledgments This work was supported by the Technology Foundation of Huizhou (No. 2017C0403019) and Huizhou University (No. HZUX1201620).
[2] Harten A, Engquist B, Osher S, Chakravarthy S. Uniformly high order essentially non-oscillatory schemes III. J Comput Phys 1987;71:231–303. [3] Liu XD, Osher S, Chan T. Weighted essentially non-oscillatory schemes. J Comput Phys 1994;115:200–12. [4] Jiang GS, Shu CW. Efficient implementation of weighted ENO schemes. J Comput Phys 1996;126:202–28. [5] Cockburn B, Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Math Comp 1989;52:411–35. [6] Cockburn B, Hou S, Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case. Math Comp 1990;54:545–81. [7] Henrick AK, Aslam TD, Powers JM. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points. J Comput Phys 2005;207:542–67. [8] Borges R, Carmona M, Costa B, Don WS. An improved weighted essentially non-oscillation scheme for hypebolic conservation laws. J Comput Phys 2008;227:3191–211. [9] Adams NA, Shariff K. A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems. J Comput Phys 1996;127:27–51. [10] Pirozzoli S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction. J Comput Phys 2002;178:81–117. [11] Ren YX, Liu M, Zhang H. A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws. J Comput Phys 2003;192:365–86. [12] Kim D, Kwon HK. A high-order accurate hybrid scheme using a central flux scheme and a WENO scheme for compressible flowfield analysis. J Comput Phys 2005;210:554–83. [13] Wang QJ, Ren YX, Sun ZS, Sun YT. Low dispersion finite volume scheme based on reconstruction with minimized dispersion and controllable dissipation. Sci China Phys Mech Astron 2013;56:423–31. [14] Martín MP, Taylor EM, Wu M, Weirs VG. A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence. J Comput Phys 2006;220:270–89. [15] Yamaleev NK, Carpenter MH. A systematic methodology for constructing high-order energy stable WENO schemes. J Comput Phys 2009;228:4248–72. [16] Hu XY, Wang Q, Adams NA. An adaptive central-upwind weighted essentially non-oscillatory scheme. J Comput Phys 2010;229:8952–65. [17] Feng H, Hu F, Wang R. A new mapped weighted essentially non-oscillatory scheme. J Sci Comput 2012;51:449–73. [18] Wang R, Feng H, Huang C. A new mapped weighted essentially non-oscillatory method using rational mapping function. J Sci Comput 2015;67:540–80. [19] Shu CW. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In: Advanced numerical approximation of nonlinear hyperbolic equations. In: Lecture notes in mathematics, 1697. Berlin: Springer-Verlag; 1998. p. 325–432. [20] Shu CW. Total-variation-diminishing time discretizations. SIAM J Sci Statist Comput 1988;9(6):1073–84. [21] Gottlieb S, Shu CW. Total variation diminishing Runge-Kutta schemes. Math Comp 1998;67:73–85. [22] Balsara DS, Shu CW. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. J Comput Phys 20 0 0;160:405–52. [23] Banks JW, Aslam T, Rider WJ. On sub-linear convergence for linearly degenerate waves in capturing schemes. J Comput Phys 2008;227:6985–7002. [24] Sod GA. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J Comput Phys 1978;107:1–31. [25] Shu CW, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes II. J Comput Phys 1989;83:32–78. [26] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comput Phys 1984;54:115–73.
References [1] Harten A, Osher S. Uniformly high order essentially non-oscillatory schemes I. SIAM J Numer Anal 1987;24:279–309.
Please cite this article as: F. Hu, The 6th-order weighted ENO schemes for hyperbolic conservation laws, Computers and Fluids (2018), https://doi.org/10.1016/j.compfluid.2018.07.008