A new smoothness indicator for improving the ... - Semantic Scholar
Journal of Computational Physics 269 (2014) 329–354
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A new smoothness indicator for improving the weighted essentially non-oscillatory scheme Ping Fan a,∗ , Yiqing Shen b , Baolin Tian c , Chao Yang a a b c
Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China Institute of Applied Physics and Computational Mathematics, Beijing 100088, China
a r t i c l e
i n f o
Article history: Received 16 September 2013 Received in revised form 10 March 2014 Accepted 14 March 2014 Available online 29 March 2014 Keywords: WENO scheme Smoothness indicator Hyperbolic conservation law Euler equation
a b s t r a c t In this work, a new smoothness indicator that measures the local smoothness of a function in a stencil is introduced. The new local smoothness indicator is defined based on the Lagrangian interpolation polynomial and has a more succinct form compared with the classical one proposed by Jiang and Shu [12]. Furthermore, several global smoothness indicators with truncation errors of up to 8th-order are devised. With the new local and global smoothness indicators, the corresponding weighted essentially non-oscillatory (WENO) scheme can present the fifth order convergence in smooth regions, especially at critical points where the first and second derivatives vanish (but the third derivatives are not zero). Also, the use of higher order global smoothness indicators incurs less dissipation near the discontinuities of the solution. Numerical experiments are conducted to demonstrate the performance of the proposed scheme. © 2014 Elsevier Inc. All rights reserved.
1. Introduction High-order accurate weighted essentially non-oscillatory (WENO) schemes have recently been developed to solve the hyperbolic conservation law
u t + f (u ) x = 0
(1)
Based on the successful essentially non-oscillatory (ENO) scheme in [10], WENO schemes use the idea of adaptive stencils in the reconstruction procedure to automatically achieve high-order accuracy and non-oscillatory property near discontinuities [22,23]. The first WENO scheme was developed by Liu, Osher and Chan [14] for a third-order of accuracy and the classic fifth-order WENO scheme with a general framework for the design of the smoothness indicators and nonlinear weights was constructed in [12], which is referred to hereafter as WENO-JS. Later, very high order WENO schemes were developed in [1] and [8]. WENO schemes use a nonlinear convex combination of all the ENO candidate sub-stencils and assign each sub-stencil a weight between 0 and 1 based on local smoothness indicators. The basic weighting strategy is to use optimal weights for each of the lower order polynomials to combine an upwind scheme of maximum order in smooth regions of the solution, while assign small weights to those lower order polynomials whose underlying stencils contain discontinuities so that
*
Corresponding author. E-mail addresses: [email protected], [email protected] (P. Fan), [email protected] (Y. Shen), [email protected] (B. Tian), [email protected] (C. Yang). http://dx.doi.org/10.1016/j.jcp.2014.03.032 0021-9991/© 2014 Elsevier Inc. All rights reserved.
330
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
the ENO property is achieved. Although being successful in a wide number of applications, the WENO methodology is still in development to improve its convergence rate in smooth regions and decrease the dissipation near the discontinuities. Henrick, Aslam and Powers pointed out that the actual convergence rate of the fifth-order WENO-JS was less than 5th for many problems [11]. They demonstrated that the smoothness indicator of WENO-JS failed to recover the maximum order of the scheme at critical points where the first or higher order derivatives of the function vanished. In the same paper, they proposed a new WENO-M scheme to disentangle this problem by modifying the smoothness indicator of WENO-JS with a mapping procedure to satisfy the sufficient criteria for fifth-order convergence. Compared to the WENO-JS scheme, the WENO-M scheme can achieve the optimal convergence order at critical points of smooth parts of the solution and reduce the numerical dissipation near discontinuities. With a different weighting formulation, Borges, Carmona, Costa, and Don [5] introduced another version of the fifth-order WENO scheme (called WENO-Z) which uses a global higher order reference value for the smoothness indicators and drives the WENO weights to the optimal values faster than the WENO-M scheme. The WENO-Z presents less dissipation than WENO-JS but at the first order critical points the convergence order is 4th and will degrade to 2nd when higher order critical points are encountered. In [6], a closed-form formula for the WENO-Z scheme of higher than fifth-order accuracy was derived. A new smoothness indicator for the WENO-Z scheme was recently developed in [9] and a new 6th-order global smoothness indicator was devised. The convergence rate of the resulted WENO scheme can recover the optimal fifth-order even at the first order critical point, while, it drops to a same lower rate as WENO-Z when a second order critical point is met. Recently, developments have also been made in designing WENO schemes on unstructured meshes in more space dimensions. For example, a third-order finite volume WENO scheme was constructed on three-dimensional tetrahedral meshes [27] and an arbitrary high order finite volume WENO scheme on unstructured grids in two and three space dimensions was developed by using the ADER (ADER stands for Arbitrary Derivative Riemann problem) approach [7]. In [13] a third-order compact central WENO scheme was presented for multidimensional conservation laws. Excepting applications for the hydrodynamics problems, the WENO schemes have also been applied for solving the magnetohydrodynamics (MHD) problems. For example, the divergence-free WENO schemes for simulating the MHD flows were designed by Balsara [2] and it was found that for MHD the ADER-WENO schemes are better than the Runge– Kutta ( R K )-WENO schemes. In [3,4] efficient techniques for WENO interpolation have been proposed in developing the ADER-WENO schemes for the divergence-free MHD. They presented very elegant and compact formulations of WENO reconstruction as well as their implementation details by expressing the interpolating functions in modal space of Hermite polynomial. In this article, we introduce a new smoothness indicator η for evaluating the local smoothness of the numerical solution in a stencil. The new local smoothness indicator is calculated based on the Lagrangian interpolation polynomial and has a more succinct form compared with the classical one of the WENO-JS scheme. The WENO scheme that equipped with the new smoothness indicator η is designated as the WENO-η scheme. Observations from the given numerical experiments manifest that the WENO-η scheme has an equivalent performance as the WENO-JS scheme. Furthermore, several higher order (up to 8th-order) global smoothness indicators τ are constructed, with which the associated WENO schemes (hereafter, denoted by WENO-Zη ) can achieve fifth convergence order in smooth regions, even at the second order critical points where the first and second derivatives vanish. Also, it is found that the use of higher order global smoothness indicators incurs less dissipation near the discontinuities of the solution. Numerical experiments are presented to show that the WENO-Zη schemes provide at least similar or improved behaviors over the WENO-Z scheme. The rest of this paper is organized as follows. In Section 2, a brief review of the WENO schemes for one-dimensional scalar conservation laws is provided. In Section 3, we present a new local smoothness indicator and several higher order global smoothness indicators for improving the WENO scheme. In Section 4, some numerical results are presented to show the capacity of the proposed WENO scheme. Section 5 concludes this paper. 2. Review of the WENO schemes For the hyperbolic conservation law (1), the flux function f (u ) is split into two parts as f (u ) = f + (u ) + f − (u ) with 0 and df − (u )/du 0. The semi-discretization form of (1) can be written as
df + (u )/du
du i (t ) dt
=−
1
x
(hi + 1 − hi − 1 ), 2
(2)
2
− + where the numerical flux is h i +1/2 = h+ i +1/2 + h i +1/2 . Hereinafter, only the positive part h i +1/2 is described and the superscript “+” is dropped for brevity. The flux of the classical fifth-order WENO-JS scheme is built through the convex combination of interpolated values R i −2+l,3 (xi +1/2 ) (l = 0, 1, 2)
hi+ 1 = 2
2 l =0
ωl R i−2+l,3 (xi+ 1 ) + O x5 , 2
in which R i −2+l,3 is the reconstruction polynomial on stencil S i −2+l,3 = { I i −2+l , I i −1+l , I i +l } ( I i ≡ {xi −1/2 , xi +1/2 })
(3)
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
R l,3 (xi + 1 ) = 2
2
αl
m =0
(4)
m =0
ωl are defined as
The non-oscillatory weights
ωl = 2
i = 0, . . . , N
clm f i +l−2+m ,
331
αm
αl =
,
dl
(ε + βl ) p
(5)
,
where dl are called the ideal weights since they generate the central upstream fifth-order scheme for the 5-points stencil S i −2,5 = { I i −2 , . . . , I i +2 } and βl are the local smoothness indicators. The general idea of the weights definition (5) is that if the function f (x) is smooth in all of the candidate stencils the smoothness indicators βl are all small and about the same size, generating weights ωl that approximate the ideal weights dl with third-order of accuracy. Then, the convex combination of all R i −2+l,3 guarantees fifth-order of accuracy to approximate the flux at the cell boundary. On the other hand, if f (x) has a discontinuity in the sub-stencil S i −2+l,3 but is smooth in at least one of the other stencils, then βl is O (1) and the corresponding weight ωl is small relatively to the other weights. The scheme is then third-order accurate and emulation of ENO is achieved. Henrick, Aslam and Powers [11] indicated that to ensure the overall scheme retaining fifth-order accuracy, it is sufficient to require
ωl± − dl = O x3
(6)
To meet the above requirement they proposed a mapping function to make ωl approximating the ideal weights dl with increased accuracy. The mapping function used in the WENO-M scheme is defined as
ω(dl + dl2 − 3dl ω + ω2 ) dl2 + ω(1 − 2dl )
gl (ω) =
(7)
Borges, Carmona, Costa, and Don [5] introduced the absolute difference between β0 and β2 to devise a new set of WENO weights
ωlz = 2
αlz
m =0
αmz
αlz = dl 1 +
,
τ5
p ,
ε + βl
where the global higher order smoothness indicator
l = 0, 1 , 2 ,
(8)
τ5 is defined as
τ5 = |β0 − β2 |
(9)
With this implementation the convex combination of the WENO-Z weights becomes closer to the central scheme than that of the WENO-JS and hence generates less dissipation. 3. New smoothness indicators 3.1. Classic smoothness indicators The smoothness indicator that evaluates the local smoothness of a function in a stencil was defined originally in [14], which can be written as
I Sl =
r −n r −1 ( f [i + l + m − r , n])2 , r −n
(10)
n =1 m =1
where [· , ·] is the nth undivided difference operator. With this smoothness indicator, a (r + 1)th-order WENO scheme is converted from a rth-order ENO scheme. To improve accuracy, Jiang and Shu [12] designed a modified smoothness indicator as x
βl =
2 m =1
x2m−1 x
1
i+ 2
dm
2
R i −r +1+l,r (x) dxm
dx
i− 1 2
For r = 3, the explicit formulae for βl in terms of the cell averaged value of f are
(11)
332
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
⎧ ⎪ ⎪ β0 = ⎪ ⎪ ⎪ ⎨ β1 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ β2 =
13 12 13 12 13 12
1
( f i −2 − 2 f i −1 + f i )2 + ( f i −2 − 4 f i −1 + 3 f i )2 ; 2
4 1
( f i −1 − 2 f i + f i +1 ) + ( f i −1 − f i +1 )2 ;
(12)
4 1
( f i − 2 f i +1 + f i +2 )2 + (3 f i − 4 f i +1 + f i +2 )2 , 4
and their Taylor series expansions at xi are
⎧ 2 13 2 ⎪ ⎪ = f x + − f β ⎪ 0 i ⎪ ⎪ 12 i ⎪ ⎪ ⎨ 2 13 2 β1 = f i x + + f ⎪ 12 i ⎪ ⎪ ⎪ ⎪ 2 13 2 ⎪ ⎪ ⎩ β2 = f x + − f i
13 1 (4) x5 + O x6 ; fi fi − fi fi 6 2
1 x4 + O x6 ; f f 3 i i
2 x4 + f f 3 i i
i
12
2 x4 − f f 3 i i
(13)
13 1 (4) x5 + O x6 fi fi − fi fi 6 2
The smoothness indicator defined by (11) has currently been widely used in the WENO schemes, with which a rth-order ENO scheme can be converted to a (2r −1)th-order WENO scheme. However, to evaluate the explicit form of βl is essentially a very hard work, since it involves an integral calculation on the square of derivatives of the reconstruction polynomials. Although the explicit formula of βl for r = 4, 5, 6 and r = 7, 8, 9 have been presented by Balsara and Shu [1] and Gerolymos, Schal and Vallet [8], respectively, the works for evaluating the explicit form of βl are still tremendously burdensome especially for large r. Apart from the above mentioned works on the explicit formulae for the smoothness indicator, it is worthy of a particular mention that new expressions that are more compact than those in Balsara and Shu [1] and Jiang and Shu [12] have also been obtained recently in [3] and [4]. By casting the problem in a modal space of Hermite polynomials
q0 (ξ ) = 1,
q1 (ξ ) = ξ,
q2 (ξ ) = ξ 2 −
1 12
,
q3 (ξ ) = ξ 3 −
3 20
ξ,
(14)
Balsara et al. [3,4] have succeeded in deriving a series of more elegant and compact formulations for the smoothness indicators. For r = 3, these smoothness indicators can be expressed explicitly as
ISi = ( f xi )2 + where
13 3
( f xxi )2 (i = 0, 1, 2),
(15)
⎧ 1 1 ⎪ ⎪ f x0 = ( f i −2 − 4 f i −1 + 3 f i ); f xx0 = ( f i −2 − 2 f i −1 + f i ) ⎪ ⎪ 2 2 ⎪ ⎨ 1 1 f x1 = ( f i +1 − f i −1 ); f xx1 = ( f i −1 − 2 f i + f i +1 ) ⎪ 2 2 ⎪ ⎪ ⎪ ⎪ ⎩ f x2 = − 1 (3 f i − 4 f i +1 + f i +2 ); f xx2 = 1 ( f i − 2 f i +1 + f i +2 ) 2
(16)
2
Although, it can be seen that the above results have an equivalent form to the classic one suggested by Jiang and Shu [12], however, for r larger than 3, the explicit formulae presented in [3] and [4] are much more compact than the classic ones by Balsara and Shu [1]. For the sake of brevity, expressions of the mentioned smoothness indicator for r = 4 and 5 will not be presented in this paper and the interested readers who want to acquire the exact form are referred to Refs. [3] and [4]. 3.2. The new local smoothness indicator In this paper, we present a new succinct formula for the local smoothness indicators defined by
ηl =
r −1
(m) 2
Pl
=
m =1
r −1
m =1
2 (m) xm P i −r +1+l,r (xi ) ,
0 l r − 1,
(17)
(m)
where P l is a polynomial that approximates the mth-order derivative of the flux function f at point xi . P i ,r (x) is the Lagrangian interpolation polynomial for approximating the value of the flux function f at any point x based upon a r-point stencil S i ,r = {xi , . . . , xi +r −1 }
P i ,r (x) =
r −1 m =0
f i +m
r −1 l=0,l=m
x − xi +l xi +m − xi +l
,
(18)
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
333
(m)
and P i ,r (x) is the m-th derivative of P i ,r (x). (m) can be written explicitly as For r = 3, the polynomials P l
⎧ 1 3 (1 ) (1 ) ⎪ ⎪ P 0 = xP i −2,3 (xi ) = f i −2 − 2 f i −1 + f i , ⎪ ⎪ 2 2 ⎪ ⎨ 1 1 (1 ) (1 ) P 1 = xP i −1,3 (xi ) = − f i −1 + f i +1 , ⎪ 2 2 ⎪ ⎪ ⎪ (1 ) ⎪ ⎩ P = xP (1) (xi ) = − 3 f i + 2 f i +1 − 1 f i +2 2 i ,3 2 2 ⎧ (2 ) 2 (2 ) P = x P i −2,3 (xi ) = f i −2 − 2 f i −1 + f i , ⎪ ⎪ ⎨ 0 (2 ) (2 ) P 1 = x2 P i −1,3 (xi ) = f i −1 − 2 f i + f i +1 , ⎪ ⎪ ⎩ (2 ) (2 ) P 2 = x2 P i ,3 (xi ) = f i − 2 f i +1 + f i +2
Accordingly, the new smoothness indicators
⎧ ⎪ ⎪ η0 = ⎪ ⎪ ⎪ ⎨ η1 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ η2 =
1 4 1 4 1 4
(19)
(20)
ηl in terms of the cell averaged value of f can be expressed explicitly as
( f i −2 − 4 f i −1 + 3 f i )2 + ( f i −2 − 2 f i −1 + f i )2 , ( f i −1 − f i +1 )2 + ( f i −1 − 2 f i + f i +1 )2 ,
(21)
(3 f i − 4 f i +1 + f i +2 )2 + ( f i − 2 f i +1 + f i +2 )2
Using the Taylor series expansion, we find
⎧ ⎪ ⎪ (1 ) ⎪ P0 = ⎪ ⎪ ⎪ ⎨ (1 ) P1 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ P (1 ) = 2 ⎧ (2 ) ⎪ ⎪ ⎨ P0 = (2 ) P1 = ⎪ ⎪ ⎩ (2 ) P2 =
(3 )
(1 )
fi
x −
fi
3
x3 +
1 (4 ) 4 f x + O x5 , 4 i
1 (3 ) 3 f i x + O x5 , 6 1 (3 ) 1 (4 ) (1 ) f i x − f i x3 − f i x4 + O x5 3 4 (1 )
fi
(22)
x +
x2 − f i(3) x3 + O x4 , (2 ) f i x2 + O x4 , (2 ) (3 ) f i x2 + f i x3 + O x4 (2 )
fi
Thus, the Taylor series expansions of
⎧ ( 1 ) 2 ⎪ (2 )2 ⎪ η = f x + fi − ⎪ 0 i ⎪ ⎪ ⎪ ⎪ ⎨ 2 η1 = f i(1) x + f i(2)2 + ⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ ⎪ ⎩ η2 = f ( 1 ) x + f ( 2 ) 2 − i
i
ηl read
2 (1 ) (3 ) x4 + fi f 3 i
(23)
1 (1 ) (4 ) (2 ) (3 ) x5 + O x6 fi fi − 2 fi fi 2
1 (1 ) (3 ) x4 + O x6 fi f 3 i
2 (1 ) (3 ) x4 − fi f 3 i
ηl are 1 + O x2 ,
(24)
1 (1 ) (4 ) (2 ) (3 ) x5 + O x6 fi fi − 2 fi fi 2
Obviously, the truncation errors of
2
ηl = f i(1) x
(25)
showing a close resemblance to the asymptotic features of βl . Nevertheless, the computation for ηl is much simplified since it needs only a simple computation for the Lagrangian interpolation polynomials. Furthermore, higher order global smoothness indicators can be easily devised with ηl , making the associated WENO scheme less dissipative than the classic ones. 3.3. New higher order global smoothness indicators (1)
τ5
The novel idea of the WENO-Z scheme introduced in [5] is the construction of the higher order global smoothness indicator τ5 from the classic local smoothness indicator βl . In this paper, we find that a new 5th-order global smoothness
334
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
indicator can be acquired from a similar operation on the new local smoothness indicator smoothness indicator can be calculated by
ηl . As in (9), the new global
τ5 = |η0 − η2 |
(26)
It is straightforward to see from (24) that the truncation error of the new
τ5 is
(1 ) (4 ) f f − 4 f (2) f (3) x5 + O x6 , i
i
i
(27)
i
which is equivalent to that introduced in [5]. In the latter part of this paper, we will show by numerical tests that the results obtained by the WENO-Z scheme with βl and ηl and the corresponding τ5 have little differences. (2)
τ6
As argued previously, it is found in this paper that the order of τ can in fact be raised to higher than 5th-order without using information outside the given large 5-point stencil S i −2,5 = {xi −2 , . . . , xi +2 }. We will not present the method for devising the 6th-order τ from βl , since it is analogous to what will be introduced in the following by using the new local smoothness indicator ηl . To build this higher order τ , we should first introduce a reference value for evaluating the smoothness of the large 5-point stencil S i −2,5 , which is written by
η5 =
2
2 (m) xm P i −2,5 (xi )
m =1
=
1
144
( f i −2 − 8 f i −1 + 8 f i +1 − f i +2 )2 + ( f i −2 − 16 f i −1 + 30 f i − 16 f i +1 + f i +2 )2
From using the Taylor series expansion, we see that
2
η5 = f i(1) x + f i(2) x2
2
(1 ) (5 )
−
fi
fi
15
6
(28)
η5 can be expanded as
x6 + O x7
(29)
On the other hand, one can note from (24) that the following combination of
4 η1 + η0 + η2
ηl can be expanded as
2 2 = f i(1) x + f i(2) x2 + O x6
Therefore, a new 6th-order global smoothness indicator
4 η1 + η0 + η2 = O x6 τ6 = η5 − 6
(30)
τ6 can be acquired by the following ways (31)
Comparing with τ5 , the rate of convergence of WENO-Z with τ6 can recover 5th-order at a first order critical point with no need of increasing the power p in (8) from 1 to 2. Also, it will be shown that results of τ6 are less dissipative than those of τ5 . (3)
τ8
Other than the above introduced τ5 and τ6 , this paper finds that the truncation error of τ can be raised up to 8 th-order by a quite different approach. The new 8th-order global smoothness indicator τ8−1 is defined by
τ8−1 = P 0(1) − P 2(1) P 0(2) + P 2(2) − 2P 1(2)
(32)
From the Taylor series expansion, one can find that the truncation error of
1 τ8−1 = f i(4) x4 + O x5 O x4 = O x8
τ8−1 is
(33)
2
It is worthy of note that there is another way to construct the 8th-order
2 2 τ8−2 = P (1) − P (1) + P (2) + P (2) − 2P (2) 0
2
0
2
1
τ , as defined by (34)
It is easy to verify that τ8−2 is also of 8th-order, while, in contrast to τ8−1 , the results of τ8−2 are relatively more dissipative than those of τ8−1 . This can be shown in the latter numerical experiments. Furthermore, comparing with τ5 and τ6 , the benefit of τ8 is that the rate of convergence of WENO-Z with τ8 can recover that of the formal 5th-order, even at a second order critical point. This will be demonstrated in the following section.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
335
Table 1a Rate of convergence at a first order critical point (k = 2) for WENO-JS and WENO-Z schemes. N
WENO-Z (p = 1)
WENO-JS
23 24 25 26 27 28 29
WENO-Z (p = 2)
Error
Order
Error
Order
Error
Order
5.8085E–02 6.1571E–03 4.8105E–04 3.2076E–05 1.1798E–06 1.2947E–08 1.1654E–10
−
1.6448E–02 9.9119E–04 5.5215E–05 3.1169E–06 1.8038E–07 1.0713E–08 6.4960E–10
– 4.0526 4.1660 4.1469 4.1110 4.0736 4.0437
5.4640E–02 4.0202E–03 1.3079E–04 3.0220E–06 6.5728E–08 1.5353E–09 3.9560E–11
– 3.7646 4.9419 5.4356 5.5229 5.4199 5.2783
3.2378 3.6780 3.9066 4.7649 6.5098 6.7956
Table 1b Rate of convergence at a first order critical point (k = 2) for WENO-Zη(τ5 ) schemes. WENO-Zη(τ5 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ5 ) ( p = 2)
Error
Order
Error
Order
1.6237E–02 9.8866E–04 5.5315E–05 3.1259E–06 1.8086E–07 1.0737E–08 6.5090E–10
– 4.0377 4.1597 4.1453 4.1113 4.0742 4.0440
5.4380E–02 4.0045E–03 1.3009E–04 2.9932E–06 6.4699E–08 1.5011E–09 3.8461E–11
– 3.7634 4.9440 5.4417 5.5318 5.4297 5.2865
Table 1c Rate of convergence at a first order critical point (k = 2) for WENO-Zη(τ6 ) schemes. WENO-Zη(τ6 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ6 ) ( p = 2)
Error
Order
Error
Order
6.1221E–03 1.4197E–04 3.2037E–06 7.9428E–08 2.1549E–09 6.2224E–11 1.9593E–12
– 5.4304 5.4697 5.3339 5.2040 5.1140 4.9891
1.7041E–02 7.2599E–05 1.7949E–06 5.7051E–08 1.7984E–09 5.6529E–11 1.8119E–12
– 7.8748 4.9532 4.9755 4.9875 4.9916 4.9634
3.4. Convergence at critical points Critical points have become a focus of discussion since in [11] it was shown that the standard value of the parameter
ε = 10−6 in fact hided the loss of accuracy of WENO-JS which achieves only third order at critical points of smooth solutions.
WENO-M was presented as a fix to this situation, for it used a mapping that corrected the weights of WENO-JS and recovered the formal fifth-order at critical points even when using very small values of ε . Furthermore, it was demonstrated in [5] that the WENO-Z scheme is fourth-order accurate at first order critical points of a smooth solution and the fifth-order convergence can be recovered by increasing parameter p in (8) from p = 1 to p = 2. However, this modification increases the dissipation of the WENO-Z scheme. We now perform a numerical experiment to compare the behavior of the classic WENO-JS, WENO-Z and the new WENO-Zη in the presence of smooth solutions containing critical points. Unless indicated otherwise, the values of parameter ε are taken as ε = 10−6 for the WENO-JS scheme and ε = 10−40 for the others. Consider the following test function
f (x) = xk exp(x),
x ∈ [−1, 1]
(35)
in which its first k − 1 derivatives f ( j ) (0) = 0, j = 0, 1, . . . , k − 1. That is, this function has a critical point of order ncp = k − 1 at x = 0. Tables 1 and 2 show the L ∞ convergence behaviors for these schemes at the first and second order critical points, i.e., k = 2 and 3, respectively. We use the fixed value of p = 2 for WENO-JS since its rate of convergence does not change when p varies. It should be mentioned particularly that in this subsection only the convergence of WENO-Zη(τ8−2 ) is considered and referred to as WENO-Zη(τ8 ). As shown in Table 1a and Table 2a, increasing value of p from 1 to 2 allows a speed up of convergence rate of WENO-Z at a first order critical point, but this implementation loses its effect at a second order critical point. What in Tables 1b and 2b shows that WENO-Zη(τ5 ) takes an equivalent performance as WENO-Z. In Tables 1c and 2c, it can be seen that WENO-Zη(τ6 ) recovers the fifth-order convergence rate at a first critical point, regardless of value of p. However, at a second order critical point, the rate of convergence of WENO-Zη(τ6 ) degrades to
336
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Table 1d Rate of convergence at a first order critical point (k = 2) for WENO-Zη(τ8 ) schemes. WENO-Zη(τ8 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ8 ) ( p = 2)
Error
Order
Error
Order
1.7427E–03 5.5593E–05 1.7949E–06 5.7051E–08 1.7984E–09 5.6529E–11 1.8598E–12
– 4.9703 4.9529 4.9755 4.9875 4.9916 4.9258
1.6775E–03 5.5604E–05 1.7949E–06 5.7051E–08 1.7984E–09 5.6529E–11 1.8119E–12
– 4.9150 4.9532 4.9755 4.9875 4.9916 4.9634
Table 2a Rate of convergence at a second order critical point (k = 3) for WENO-JS and WENO-Z schemes. N
23 24 25 26 27 28 29
WENO-Z (p = 1)
WENO-JS
WENO-Z (p = 2)
Error
Order
Error
Order
Error
Order
1.6430E–01 3.5553E–02 6.3056E–03 5.0357E–05 8.2116E–07 1.4425E–08 2.7234E–10
−
1.4191E–01 3.0455E–02 7.0249E–03 1.6684E–03 3.9888E–04 9.4831E–05 2.2365E–05
– 2.2202 2.1161 2.0740 2.0644 2.0725 2.0841
1.6429E–01 3.5573E–02 8.3065E–03 2.0061E–03 4.9003E–04 1.1848E–04 2.7570E–05
– 2.2074 2.0985 2.0498 2.0335 2.0482 2.1035
2.2083 2.4953 6.9683 5.9384 5.8310 5.7270
Table 2b Rate of convergence at a second order critical point (k = 3) for WENO-Zη(τ5 ) schemes. N
23 24 25 26 27 28 29
WENO-Zη(τ5 ) ( p = 1)
WENO-Z
η(τ5 ) ( p = 2)
Error
Order
Error
Order
1.4080E–01 3.0176E–02 6.9518E–03 1.6488E–03 3.9351E–04 9.3377E–05 2.1994E–05
– 2.2222 2.1179 2.0760 2.0669 2.0753 2.0860
1.6423E–01 3.5562E–02 8.3032E–03 2.0048E–03 4.8918E–04 1.1790E–04 2.7301E–05
– 2.2073 2.0986 2.0502 2.0350 2.0528 2.1105
Table 2c Rate of convergence at a second order critical point (k = 3) for WENO-Zη(τ6 ) schemes. N 23 24 25 26 27 28 29
WENO-Zη(τ6 ) ( p = 1) Error
Order
WENO-Z Error
η(τ6 ) ( p = 2)
1.2804E–01 2.5667E–02 5.6184E–03 1.2991E–03 3.1093E–04 7.5945E–05 1.8759E–05
– 2.3186 2.1917 2.1126 2.0629 2.0336 2.0174
1.6377E–01 3.5277E–02 8.1766E–03 1.9626E–03 4.7987E–04 1.1856E–04 2.9457E–05
Order – 2.2149 2.1092 2.0587 2.0321 2.0170 2.0089
second order, as do WENO-Z and WENO-Zη(τ5 ). Tables 1d and 2d illustrate that WENO-Zη(τ8 ) recovers the fifth-order convergence rate at a first critical point. At a second critical point, rate of convergence of WENO-Zη(τ8 ) can achieve 5th-order when value of p is increased from 1 to 2. 4. Numerical examples In this section, we compare the numerical performance of WENO-η and WENO-Zη with the classic WENO-JS and WENOZ. The numerical presentation of this section starts with the solution of the scalar advection equation and followed by the solution of the one-dimensional Euler equations with Riemann initial value problems such as the Lax, Sod, shock-density wave interaction and the interacting blast problems. Finally, the two-dimensional problems on the Richtmyer–Meshkov instability, shock-vortex interaction, Rayleigh–Taylor instability, double-Mach shock reflection and the forward facing step problems are simulated. The time dependent problems are all solved via the 3rd-order Runge–Kutta TVD method [20] with CFL = 0.5.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
337
Table 2d Rate of convergence at a second order critical point (k = 3) for WENO-Zη(τ8 ) schemes. WENO-Zη(τ8 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ8 ) ( p = 2)
Error
Order
Error
Order
9.3881E–02 9.5484E–03 1.4256E–03 1.1070E–04 7.4017E–06 4.6962E–07 2.9410E–08
– 3.2975 2.7437 3.6868 3.9027 3.9783 3.9971
1.5496E–01 1.3660E–02 1.6515E–03 3.5043E–05 5.5463E–07 8.6262E–09 1.3503E–10
– 3.5039 3.0481 5.5585 5.9815 6.0067 5.9974
Fig. 1. Numerical solutions of the linear advection equation (36) with N = 27 at t = 6: (a) presents the comparison between the exact solution, the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
4.1. 1D scalar test problem Consider the one-dimensional linear advection equation
∂C ∂C + = 0, ∂t ∂x
x ∈ [−1, 1]
with an initial condition given by
(36)
338
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 2. Numerical solutions of the Lax problem with N = 200 at t = 1.3: (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
C (x, 0) =
1,
x ∈ [−0.5, 0.5),
0,
otherwise
(37)
The advection equation (36) was solved until the final time t = 6 on the interval [−1, 1] with 27 grid points and the periodic boundary conditions. Figs. 1(a) and (a-1) compare the results between the exact solution, the WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ) schemes. It is seen that the new WENO-η and WENO-Zη(τ5 ) display nearly equal performance to the WENO-JS and WENO-Z, respectively. Furthermore, the results of WENO-Zη with different global smoothness indicators τ are compared and presented in Figs. 1(b) and (b-1). It can be observed that among these schemes the WENO-Zη(τ5 ) and WENO-Zη(τ8−1 ) achieve, respectively, the largest and least dissipations near the discontinuities. 4.2. 1D Euler systems In this section, we present numerical experiment with the one-dimensional Euler equations for gas dynamics in strong conservation form
∂U ∂F + = 0, ∂t ∂x where
(38)
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
339
Fig. 3. Numerical solutions of the Sod problem with N = 200 at t = 0.2: (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
U = (ρ , ρ u , E ) T , F=
ρ u , ρ u 2 + P , ( E + P )u
T
(39)
The equation of state is given by
P = (γ − 1) E −
1 2
ρ u 2 , γ = 1.4,
(40)
where ρ , u, P and E are the density, velocity, pressure and total energy, respectively. These equations are solved by following the general WENO methodology [12]. Specifically, the Roe approximation is used for the characteristic decomposition at the cell faces and the Lax–Friedrichs formulation is used for the numerical fluxes. The numerical experiments are conducted using 200 grid points. (I). Riemann initial value problem of Lax For the Lax problem, the initial condition is
(ρ , U , p ) =
(0.445, 0.698, 0.3528), −5 x < 0, 0 x 5, (0.5, 0, 0.571),
(41)
340
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 4. Numerical solutions of the shock-density interaction with N = 200 at t = 1.8: (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
and the final time is 1.3. The simulated density of the Lax problem is shown in Figs. 2. We refer to the solution computed by the WENO-JS with N = 4000 grid points as the ‘exact’ solution. From Figs. 2(a) and (a-1) we see that the new WENO-η and WENO-Zη(τ5 ) schemes show equivalent performance as those of WENO-JS and WENO-Z, respectively. Figs. 2(b) and (b-1) show that the WENO-Zη(τ8−1 ) presents the least dissipation near the discontinuities compared to WENO-Zη(τ5 ), WENO-Zη(τ6 ) and WENO-Zη(τ8−2 ). (II). Riemann initial value problem of Sod For the Sod problem, the initial condition is
(ρ , U , p ) =
x 0, (1.0, 0, 1.0), (0.125, 0, 0.1), x > 0,
(42)
and the final time is 0.2. The simulated density of the Sod problem is shown in Figs. 3. As before, the performance of WENO-η and WENO-Z η(τ5 ) are nearly the same as those of WENO-JS and WENO-Z, respectively, while, the WENO-Zη(τ8−1 ) shows the least dissipation than WENO-Zη(τ5 ), WENO-Zη(τ6 ) and WENO-Zη(τ8−2 ).
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
341
Fig. 5. Numerical solutions of the interacting blast waves problem with N = 200 at t = 0.038; (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) gives the local enlarged plot of (a).
(III). Shock-density wave interaction The initial condition for the shock-density wave interaction problem [21] is given by
(ρ , U , p ) =
(3.857143, 2.629369, 10.3333), x < −4, x −4 (1 + 0.2 sin(5x), 0, 1.0),
(43)
This case represents a Mach 3 shock wave interacting a sine entropy wave. The results at time t = 1.8 with mesh size of N = 200 are plotted in Figs. 4. As before, WENO-η and WENO-Zη(τ5 ) show good agreement with WENO-JS and WENO-Z, respectively. Moreover, the results of WENO-Zη(τ6 ) and WENO-Zη(τ8−1 ) are much less dissipative than those of the other schemes. (IV). Interacting blast wave problem The initial condition for the interacting blast wave example [28] is
(ρ , U , p ) =
(1, 0, 1000), 0 x < 0.1, (1, 0, 0.01), 0.1 x < 0.9, (1, 0, 100), 0.9 x 1,
(44)
342
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 6. Initial location of the interface in the RM instability.
and the final time is 0.038. Since this problem is much more stringent than those previous ones, only results from WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ) are presented. Simulations under the use of WENO-Zη schemes equipped with higher order global smoothness indicators are found to be unstable, which may be largely ascribed to the less dissipation of these schemes. From Figs. 5(a) and (b) it can be seen that the new WENO-Zη(τ5 ) scheme is as good as WENO-Z, both of which are less dissipative than WENO-JS and WENO-η . 4.3. 2D Euler systems In this subsection, we extend the new WENO-η and WENO-Zη schemes to solve the 2D compressible system of the form:
∂U ∂F ∂G + + = 0, ∂t ∂x ∂x
(45)
where
U = (ρ , ρ u , ρ v , E ) T , F= G=
T
ρ u , ρ u 2 + P , ρ uv , ( E + P )u ,
ρ v , ρ uv , ρ v 2 + P , ( E + P ) v
T
(46)
The equation of state is given by
1 2 2 , P = (γ − 1) E − ρ u + v 2
γ = 1.4
(47)
Here ρ , u, v, P and E are the density, components of velocity in the x and y coordinate directions, pressure and total energy, respectively. The Steger–Warming flux vector splitting method is used for the 2D problems in this paper. (I) Richtmyer–Meshkov Instability (RMI) In this test, we consider the standard RMI unstable interface problem [15,17,24]. The initial condition in the tube is composed of an interface separating two fluids of different densities and a shock wave approaching the interface. The shock tube is [0, 1.5] × [0, 1] in the 2D domain. It is well known that this interface becomes unstable after the passage of a shock wave, irrespective of the side of the fluid that the shock is incident upon. Here we consider a single mode perturbation of an interface between two fluids, and the initial location of the interface (as shown in Fig. 6) is represented by x = x0 + a0 cos(2π y ), where x0 = 0.2 is the location of the unperturbed interface, a0 = 0.05 is the amplitude of the perturbation. We use the ideal gas equation of state for air and SF6 with γ = 1.4 and γ = 1.093, respectively. We choose the density ratio equal to 6 with a corresponding Atwood number
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 7. Density profiles of the RMI driven by a planar shock: (a) WENO-JS, (b) WENO-Z, (c) WENO-Zη(τ6 ), (d) WENO-Zη(τ8−1 ), (e) WENO-Zη(τ8−2 ).
343
344
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 8. The pressure contour of the shock/vortex interaction by WENO-Zη(τ6 ).
Fig. 9(a). Pressure distributions of the shock/vortex interaction on central line y = 0.5.
0.71. To trigger the instability, a planar Mach 1.2 shock wave propagates from the left to the interface. This gives us one example that the interface is accelerated by a shock wave coming from the light-fluid to the heavy-fluid region, and the resulting wave pattern after the interaction would consist of a transmitting shock wave, an interface, and a reflected rarefaction wave. In Figs. 7 we compare the numerical performances of the WENO-JS, WENO-Z, WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ) at t = 6 with a mesh system 241 × 161 is used. An examine of these results reveals that WENO-Zη(τ6 ) and WENO-Zη(τ8−1 ) display better resolution of the fine structure near the interface than the other WENO schemes. (II) Two-dimensional shock vortex interaction In this subsection the two-dimensional shock vortex interaction problem taken from [12] is solved. This problem describes the interaction between a stationary shock and a vortex. The computational domain is taken to be [0, 2] × [0, 1]. A stationary Mach 1.1 shock is positioned at x = 0.5 and normal to the x-axis. Its left state is (ρ , u , v , p ) = (1, 1.1γ 1/2 , 0, 1). A small vortex is superimposed to the flow on the left of the shock and is centered at (xc , y c ) = (0.25, 0.5). The vortex is described as a perturbation to the velocity (u , v ), temperature T = p /ρ , and entropy S = ln( p /ρ γ ) of the mean flow and the parameters are taken the same as in [12,19]. The time step is taken as follows [16]
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
345
Fig. 9(b). Locally enlarged plot of Fig. 9(a).
Fig. 10. Density contours for the Rayleigh–Taylor instability: (a) WENO-JS; (b) WENO-Z; (c) WENO-Zη(τ6 ); (d) WENO-Zη(τ8−1 ); (e) WENO-Zη(τ8−2 ).
t x t y , t x + t y x t x = , maxi , j (|u i , j | + c i , j )
t = σ
t y =
y , maxi , j (| v i , j | + c i , j )
(48)
where c is the speed of sound and σ = 0.5 is the CFL number. Fig. 8 is the pressure contours at t = 0.6 calculated by the WENO-Zη(τ6 ) scheme. Figs. 9(a) and 9(b) are the comparisons of the pressure along the center line of y = 0.5. In order to show the accuracy of the new scheme, the results obtained by the WENO-Z scheme with a refined mesh of 2001 × 401 is given as the “exact” solution. The results with the coarse mesh of 251 × 101 are compared. Figs. 9(a) and 9(b) show that the WENO-Zη(τ6 ) has an equivalent effect to the WENO-Zη(τ8−1 ), both of which have a sharper shock profile than WENO-JS, WENO-Z and WENO-Zη(τ8−2 ) due to lower numerical dissipations.
346
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 11. Density contours of the double-Mach shock reflection as computed by the WENO-η scheme at time t = 0.2. The mesh resolution is 800 × 200 zone and thirty contours were fit between a range of 1.3 and 22.
(III) Two-dimensional Rayleigh–Taylor instability
This problem has been simulated by many authors [18,25,26] to test the numerical dissipation of various high order shock capturing schemes. The computational domain is [0, 0.25] × [0, 1], and the initial conditions are
⎧ ⎪ ⎨ (2, 0, −0.025 53 ρp cos(8π x), 2 y + 1), if 0 y < 12 (ρ , u , v , p ) = ⎪ ⎩ (1, 0, −0.025 5 p cos(8π x), y + 3 ), if 1 y < 1 3ρ 2 2
(49)
This problem represents the interface instability between fluids with different densities when an acceleration is directed from heavy fluid to light one. The gravitational effect is introduced by adding ρ and ρ v to the right hand side of the y-momentum and the energy equations, respectively. Reflective boundary conditions are imposed for the left and right boundaries, and (ρ , u , v , p ) = (1, 0, 0, 2.5) and (ρ , u , v , p ) = (2, 0, 0, 1) are set as top and bottom boundary conditions, respectively. The simulation time is t = 1.95. As analyzed in [18,25], because the inviscid Euler equation is solved, the details of the complex structures due to the physical instability are related to the specific form of numerical viscosity of the scheme. From the density contours with meshes of 121× 481 plotted in Fig. 10, it can be seen that the WENO-Zη(τ8−1 ) scheme obtains more complex structures than the other schemes.
(IV) Double Mach reflection of a strong shock
This problem has been carefully examined by Woodward and Colella [28]. We use the same setup as the authors did to test the problem where a Mach 10 shock hits a reflecting wall at the bottom of the domain given by [0, 4] × [0, 1]. The angle between the shock and the horizontal axis is 60◦ . The equations are the two-dimensional Euler equations (γ = 1.4) and initial conditions are
(ρ , u , v , p ) =
(8, 8.25 cos π6 , −8.25 sin π6 , 116.5), x < x0 + (1.4, 0, 0, 1.0),
√y
3
y x x0 + √
(50)
3
with x0 = 1/6. Boundary conditions at x = 0 are inflow, with post-shock values as above, and the right boundary is set to be outflow boundary. At y = 0, reflecting boundary conditions are applied to the interval [x0 , 4]. The top boundary condition imposes the exact motion of a Mach 10 shock in the flow variables. Fig. 11 shows the density variable at t = 0.2 in [0, 3] × [0, 1] with a mesh resolution of 800 × 200 zones as computed by the WENO-η scheme. In Figs. 12 we display the region around the double Mach stems in order to observe more clearly the
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
347
Fig. 12. Density contours of the double-Mach shock reflection as computed by the WENO-JS (a); WENO-η (b); WENO-Z (c); WENO-Zη(τ5 ) (d); WENO-Zη(τ6 ) (e); WENO-Zη(τ8−1 ) (f); and WENO-Zη(τ8−2 ) schemes (g). The mesh resolution is 800× 200 zone and thirty contours were fit between a range of 1.3 and 22.
numerical solutions of the different WENO schemes. In general, the global structure of the solution is very similar across different schemes. However, the dissipation of the various schemes can be distinguished by the number of small vortices generated along the slip line and the wall jet behind the lower half of the right moving shock. As is shown clearly in Figs. 12, the WENO-Zη(τ6 ) and WENO-Zη(τ8−1 ) exhibit the least dissipation, while, the WENO-JS and WENO-η are the most dissipative. Fig. 13 shows the density variable under different mesh resolutions which are 400 × 100, 800 × 200 and 1600 × 400, respectively. It is clear that the small vortices along the slip line are captured and become gradually obvious with the increase of the mesh resolutions. (V) Forward facing step problem This test problem was first presented by Woodward and Colella [28]. Recently, Balsara et al. [3] carried out a resolution study using the fourth order ADER-WENO scheme and showed that increasing the resolution enabled them to capture important details such as the roll-up of the vortex sheet via Kelvin–Helmholtz instability. In this paper, we also successfully captured the roll-up of the vortex sheet not only from increasing the grid resolution but also from using the new less
348
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 12. (continued)
dissipative WENO-Zη schemes. Furthermore, the robustness of the new WENO schemes gets substantially demonstrated from the success of simulating this severely stringent problem. The setup of the problem is the following: the wind tunnel spans a domain of [0, 3] × [0, 1]. A forward-facing step is located at (0.6, 0.2). The problem is initialized by a right-going Mach 3 flow with a density of 1. 4 and a pressure of unity. Reflective boundary conditions are applied along the walls of the tunnel. Inflow and outflow boundary conditions are applied, respectively, at the left and right boundaries. The simulation was run until a time of 4 and the ratio of specific heats is given by 1.4. Figs. 14 and 15 show the density at the final time at mesh resolutions of 600 × 200 and 1200× 400 zones, respectively. In all of these simulations, five 5th-order WENO schemes are used, which, in turns, are (a) WENO-JS; (b) WENO-η ; (c) WENOZ; (d) WENO-Zη(τ5 ) and (e) WENO-Zη(τ8−2 ). All the shocks are found to have sharp profiles and captured properly on the computing grid. From comparing Figs. 14 and Figs. 15, we see that with the increase of the mesh resolution from 600 × 200 to 1200 × 400, the vortex sheets are observed more evidently that emanates from the Mach stem. Moreover, the roll-up of the vortex sheet becomes the most intensely when the WENO-Zη(τ8−2 ) scheme is used. This presents a powerful evidence for the advantage of the new WENO-Zη scheme that has the less dissipation and hence has a better resolution in capturing details of the flowing structures.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
349
Fig. 13. Density contours of the double-Mach shock reflection as computed by the WENO-JS and WENO-η schemes. The mesh resolutions are 400 × 100 in (a) and (a-1); 800 × 200 in (b) and (b-1); and 1600 × 400 in (c) and (c-1). Thirty equally spaced contours are shown ranging from 1.3 to 22.
350
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 14. Density contours of the forward facing step problem as computed by the WENO-JS (a); WENO-η (b); WENO-Z (c); WENO-Zη(τ5 ) (d); and WENO-Zη(τ8−2 ) schemes (e). The resolution is 600 × 200 zones and time t = 4. Thirty equally spaced contours are shown ranging from 0.1 to 6.2.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
351
Fig. 14. (continued)
5. Conclusions
In this paper, we introduce a new local smoothness indicator for the WENO scheme which is defined based on the Lagrangian interpolation polynomial and is more compactable than the classical ones. The WENO scheme using the new local smoothness indicator, named as WENO-η , has an equal ENO property as the classic WENO-JS scheme in dealing with the shock-capturing problem. Furthermore, several higher order global smoothness indicators are devised, ensuring the resulted new WENO-Zη scheme can provide the fifth convergence order in smooth regions, especially at the first and second order critical points. Numerical experiments are run to show that the WENO-Zη scheme with 5th-order global smoothness indicator τ5 achieves a similar behavior as WENO-Z and can induce less dissipation near the discontinuities of the solution if higher order global smoothness indicators τ6 , τ8−1 and τ8−2 are used.
Acknowledgements
This work is supported by the NSFC (11172299), 973 Program (2012CB224806), CAS Program for Cross & Cooperative Team of the Science & Technology Innovation, and the National Natural Science Foundation of China (21376243). The helpful comments from the anonymous reviewers are also appreciated.
352
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 15. Density contours of the forward facing step problem as computed by the WENO-JS (a); WENO-η (b); WENO-Z (c); WENO-Zη(τ5 ) (d); and WENO-Zη(τ8−2 ) schemes (e). The resolution is 1200 × 400 zones and time t = 4. Thirty contours were fit between a range of 0.1 and 6.2.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
353
Fig. 15. (continued)
References [1] D. Balsara, C.W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy, J. Comput. Phys. 160 (2000) 405–452. [2] D.S. Balsara, Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics, J. Comput. Phys. 228 (2009) 5040–5056. [3] D.S. Balsara, T. Rumpf, M. Dumbser, C.-D. Munz, Efficient, high-accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics, J. Comput. Phys. 228 (2009) 2480–2516. [4] D.S. Balsara, M. Dumbser, C. Meyer, H. Du, Z. Xu, Efficient implementation of ADER schemes for Euler and magnetohydrodynamic flow on structured meshes – comparison with Runge–Kutta methods, J. Comput. Phys. 235 (2013) 934–969. [5] R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (2008) 3101–3211. [6] M. Castro, B. Costa, W.S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys. 230 (2011) 1766–1792. [7] M. Dumbser, M. Käser, Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems, J. Comput. Phys. 221 (2007) 693–723. [8] G.A. Gerolymos, D. Schal, I. Vallet, Very high order WENO schemes, J. Comput. Phys. 228 (2009) 8481–8524. [9] Y. Ha, C.H. Kim, Y.J. Lee, J. Yoon, An improved weighted essentially non-oscillatory scheme with a new smoothness indicator, J. Comput. Phys. 232 (2013) 68–86. [10] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order essentially non-oscillatory schemes, III, J. Comput. Phys. 71 (1987) 231–303. [11] A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points, J. Comput. Phys. 207 (2005) 542–567. [12] G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [13] D. Levy, G. Puppo, G. Russo, Compact central WENO schemes for multidimensional conservation laws, SIAM J. Sci. Comput. 22 (2000) 656–672. [14] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [15] K.O. Mikaelian, Rayleigh–Taylor and Richtmyer–Meshkov instabilities and mixing in stratified cylindrical shells, Phys. Fluids 17 (2005) 094105. [16] S. Pirozzoli, Conservative hybrid compact-WENO schemes for shock–turbulence interaction, J. Comput. Phys. 178 (2002) 81–117.
354
[17] [18] [19] [20] [21] [22] [23]
[24] [25] [26] [27] [28]
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
R.D. Richtmyer, Taylor instability in shock acceleration of compressible fluids, Commun. Pure Appl. Math. 8 (1960) 297–319. J. Shi, Y.T. Zhang, C.W. Shu, Resolution of high order WENO schemes for complicated flow structures, J. Comput. Phys. 186 (2003) 690–696. Y.Q. Shen, G.C. Zha, Generalized finite compact difference scheme for shock/complex flowfield interaction, J. Comput. Phys. 230 (2011) 4419–4436. C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78. C.W. Shu, High order weighted essentially non-oscillatory schemes for convection dominated problems, SIAM Rev. 51 (2009) 82–126. C.W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: B. Cockburn, C. Johnson, C.-W. Shu, E. Tadmor (Eds.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, in: Lecture Notes in Mathematics, Springer-Verlag, Berlin/New York, 1998, pp. 325–432. B.L. Tian, W.D. Shen, S. Jiang, S.H. Wang, Y. Liu, A global arbitrary Lagrangian–Eulerian method for stratified Richtmyer–Meshkov instability, Comput. Fluids 46 (2011) 113–121. Z.F. Xu, C.-W. Shu, Anti-diffusive flux corrections for high order finite difference WENO schemes, J. Comput. Phys. 205 (2005) 458–485. Y.N. Young, H. Tufo, A. Dubey, R. Rosner, On the miscible Rayleigh–Taylor instability: two and three dimensions, J. Fluid Mech. 447 (2001) 377–408. Y.-T. Zhang, C.-W. Shu, Third order WENO scheme on three dimensional tetrahedral meshes, Commun. Comput. Phys. 5 (2009) 836–848. P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115–173.
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A new smoothness indicator for improving the weighted essentially non-oscillatory scheme Ping Fan a,∗ , Yiqing Shen b , Baolin Tian c , Chao Yang a a b c
Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China Institute of Applied Physics and Computational Mathematics, Beijing 100088, China
a r t i c l e
i n f o
Article history: Received 16 September 2013 Received in revised form 10 March 2014 Accepted 14 March 2014 Available online 29 March 2014 Keywords: WENO scheme Smoothness indicator Hyperbolic conservation law Euler equation
a b s t r a c t In this work, a new smoothness indicator that measures the local smoothness of a function in a stencil is introduced. The new local smoothness indicator is defined based on the Lagrangian interpolation polynomial and has a more succinct form compared with the classical one proposed by Jiang and Shu [12]. Furthermore, several global smoothness indicators with truncation errors of up to 8th-order are devised. With the new local and global smoothness indicators, the corresponding weighted essentially non-oscillatory (WENO) scheme can present the fifth order convergence in smooth regions, especially at critical points where the first and second derivatives vanish (but the third derivatives are not zero). Also, the use of higher order global smoothness indicators incurs less dissipation near the discontinuities of the solution. Numerical experiments are conducted to demonstrate the performance of the proposed scheme. © 2014 Elsevier Inc. All rights reserved.
1. Introduction High-order accurate weighted essentially non-oscillatory (WENO) schemes have recently been developed to solve the hyperbolic conservation law
u t + f (u ) x = 0
(1)
Based on the successful essentially non-oscillatory (ENO) scheme in [10], WENO schemes use the idea of adaptive stencils in the reconstruction procedure to automatically achieve high-order accuracy and non-oscillatory property near discontinuities [22,23]. The first WENO scheme was developed by Liu, Osher and Chan [14] for a third-order of accuracy and the classic fifth-order WENO scheme with a general framework for the design of the smoothness indicators and nonlinear weights was constructed in [12], which is referred to hereafter as WENO-JS. Later, very high order WENO schemes were developed in [1] and [8]. WENO schemes use a nonlinear convex combination of all the ENO candidate sub-stencils and assign each sub-stencil a weight between 0 and 1 based on local smoothness indicators. The basic weighting strategy is to use optimal weights for each of the lower order polynomials to combine an upwind scheme of maximum order in smooth regions of the solution, while assign small weights to those lower order polynomials whose underlying stencils contain discontinuities so that
*
Corresponding author. E-mail addresses: [email protected], [email protected] (P. Fan), [email protected] (Y. Shen), [email protected] (B. Tian), [email protected] (C. Yang). http://dx.doi.org/10.1016/j.jcp.2014.03.032 0021-9991/© 2014 Elsevier Inc. All rights reserved.
330
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
the ENO property is achieved. Although being successful in a wide number of applications, the WENO methodology is still in development to improve its convergence rate in smooth regions and decrease the dissipation near the discontinuities. Henrick, Aslam and Powers pointed out that the actual convergence rate of the fifth-order WENO-JS was less than 5th for many problems [11]. They demonstrated that the smoothness indicator of WENO-JS failed to recover the maximum order of the scheme at critical points where the first or higher order derivatives of the function vanished. In the same paper, they proposed a new WENO-M scheme to disentangle this problem by modifying the smoothness indicator of WENO-JS with a mapping procedure to satisfy the sufficient criteria for fifth-order convergence. Compared to the WENO-JS scheme, the WENO-M scheme can achieve the optimal convergence order at critical points of smooth parts of the solution and reduce the numerical dissipation near discontinuities. With a different weighting formulation, Borges, Carmona, Costa, and Don [5] introduced another version of the fifth-order WENO scheme (called WENO-Z) which uses a global higher order reference value for the smoothness indicators and drives the WENO weights to the optimal values faster than the WENO-M scheme. The WENO-Z presents less dissipation than WENO-JS but at the first order critical points the convergence order is 4th and will degrade to 2nd when higher order critical points are encountered. In [6], a closed-form formula for the WENO-Z scheme of higher than fifth-order accuracy was derived. A new smoothness indicator for the WENO-Z scheme was recently developed in [9] and a new 6th-order global smoothness indicator was devised. The convergence rate of the resulted WENO scheme can recover the optimal fifth-order even at the first order critical point, while, it drops to a same lower rate as WENO-Z when a second order critical point is met. Recently, developments have also been made in designing WENO schemes on unstructured meshes in more space dimensions. For example, a third-order finite volume WENO scheme was constructed on three-dimensional tetrahedral meshes [27] and an arbitrary high order finite volume WENO scheme on unstructured grids in two and three space dimensions was developed by using the ADER (ADER stands for Arbitrary Derivative Riemann problem) approach [7]. In [13] a third-order compact central WENO scheme was presented for multidimensional conservation laws. Excepting applications for the hydrodynamics problems, the WENO schemes have also been applied for solving the magnetohydrodynamics (MHD) problems. For example, the divergence-free WENO schemes for simulating the MHD flows were designed by Balsara [2] and it was found that for MHD the ADER-WENO schemes are better than the Runge– Kutta ( R K )-WENO schemes. In [3,4] efficient techniques for WENO interpolation have been proposed in developing the ADER-WENO schemes for the divergence-free MHD. They presented very elegant and compact formulations of WENO reconstruction as well as their implementation details by expressing the interpolating functions in modal space of Hermite polynomial. In this article, we introduce a new smoothness indicator η for evaluating the local smoothness of the numerical solution in a stencil. The new local smoothness indicator is calculated based on the Lagrangian interpolation polynomial and has a more succinct form compared with the classical one of the WENO-JS scheme. The WENO scheme that equipped with the new smoothness indicator η is designated as the WENO-η scheme. Observations from the given numerical experiments manifest that the WENO-η scheme has an equivalent performance as the WENO-JS scheme. Furthermore, several higher order (up to 8th-order) global smoothness indicators τ are constructed, with which the associated WENO schemes (hereafter, denoted by WENO-Zη ) can achieve fifth convergence order in smooth regions, even at the second order critical points where the first and second derivatives vanish. Also, it is found that the use of higher order global smoothness indicators incurs less dissipation near the discontinuities of the solution. Numerical experiments are presented to show that the WENO-Zη schemes provide at least similar or improved behaviors over the WENO-Z scheme. The rest of this paper is organized as follows. In Section 2, a brief review of the WENO schemes for one-dimensional scalar conservation laws is provided. In Section 3, we present a new local smoothness indicator and several higher order global smoothness indicators for improving the WENO scheme. In Section 4, some numerical results are presented to show the capacity of the proposed WENO scheme. Section 5 concludes this paper. 2. Review of the WENO schemes For the hyperbolic conservation law (1), the flux function f (u ) is split into two parts as f (u ) = f + (u ) + f − (u ) with 0 and df − (u )/du 0. The semi-discretization form of (1) can be written as
df + (u )/du
du i (t ) dt
=−
1
x
(hi + 1 − hi − 1 ), 2
(2)
2
− + where the numerical flux is h i +1/2 = h+ i +1/2 + h i +1/2 . Hereinafter, only the positive part h i +1/2 is described and the superscript “+” is dropped for brevity. The flux of the classical fifth-order WENO-JS scheme is built through the convex combination of interpolated values R i −2+l,3 (xi +1/2 ) (l = 0, 1, 2)
hi+ 1 = 2
2 l =0
ωl R i−2+l,3 (xi+ 1 ) + O x5 , 2
in which R i −2+l,3 is the reconstruction polynomial on stencil S i −2+l,3 = { I i −2+l , I i −1+l , I i +l } ( I i ≡ {xi −1/2 , xi +1/2 })
(3)
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
R l,3 (xi + 1 ) = 2
2
αl
m =0
(4)
m =0
ωl are defined as
The non-oscillatory weights
ωl = 2
i = 0, . . . , N
clm f i +l−2+m ,
331
αm
αl =
,
dl
(ε + βl ) p
(5)
,
where dl are called the ideal weights since they generate the central upstream fifth-order scheme for the 5-points stencil S i −2,5 = { I i −2 , . . . , I i +2 } and βl are the local smoothness indicators. The general idea of the weights definition (5) is that if the function f (x) is smooth in all of the candidate stencils the smoothness indicators βl are all small and about the same size, generating weights ωl that approximate the ideal weights dl with third-order of accuracy. Then, the convex combination of all R i −2+l,3 guarantees fifth-order of accuracy to approximate the flux at the cell boundary. On the other hand, if f (x) has a discontinuity in the sub-stencil S i −2+l,3 but is smooth in at least one of the other stencils, then βl is O (1) and the corresponding weight ωl is small relatively to the other weights. The scheme is then third-order accurate and emulation of ENO is achieved. Henrick, Aslam and Powers [11] indicated that to ensure the overall scheme retaining fifth-order accuracy, it is sufficient to require
ωl± − dl = O x3
(6)
To meet the above requirement they proposed a mapping function to make ωl approximating the ideal weights dl with increased accuracy. The mapping function used in the WENO-M scheme is defined as
ω(dl + dl2 − 3dl ω + ω2 ) dl2 + ω(1 − 2dl )
gl (ω) =
(7)
Borges, Carmona, Costa, and Don [5] introduced the absolute difference between β0 and β2 to devise a new set of WENO weights
ωlz = 2
αlz
m =0
αmz
αlz = dl 1 +
,
τ5
p ,
ε + βl
where the global higher order smoothness indicator
l = 0, 1 , 2 ,
(8)
τ5 is defined as
τ5 = |β0 − β2 |
(9)
With this implementation the convex combination of the WENO-Z weights becomes closer to the central scheme than that of the WENO-JS and hence generates less dissipation. 3. New smoothness indicators 3.1. Classic smoothness indicators The smoothness indicator that evaluates the local smoothness of a function in a stencil was defined originally in [14], which can be written as
I Sl =
r −n r −1 ( f [i + l + m − r , n])2 , r −n
(10)
n =1 m =1
where [· , ·] is the nth undivided difference operator. With this smoothness indicator, a (r + 1)th-order WENO scheme is converted from a rth-order ENO scheme. To improve accuracy, Jiang and Shu [12] designed a modified smoothness indicator as x
βl =
2 m =1
x2m−1 x
1
i+ 2
dm
2
R i −r +1+l,r (x) dxm
dx
i− 1 2
For r = 3, the explicit formulae for βl in terms of the cell averaged value of f are
(11)
332
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
⎧ ⎪ ⎪ β0 = ⎪ ⎪ ⎪ ⎨ β1 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ β2 =
13 12 13 12 13 12
1
( f i −2 − 2 f i −1 + f i )2 + ( f i −2 − 4 f i −1 + 3 f i )2 ; 2
4 1
( f i −1 − 2 f i + f i +1 ) + ( f i −1 − f i +1 )2 ;
(12)
4 1
( f i − 2 f i +1 + f i +2 )2 + (3 f i − 4 f i +1 + f i +2 )2 , 4
and their Taylor series expansions at xi are
⎧ 2 13 2 ⎪ ⎪ = f x + − f β ⎪ 0 i ⎪ ⎪ 12 i ⎪ ⎪ ⎨ 2 13 2 β1 = f i x + + f ⎪ 12 i ⎪ ⎪ ⎪ ⎪ 2 13 2 ⎪ ⎪ ⎩ β2 = f x + − f i
13 1 (4) x5 + O x6 ; fi fi − fi fi 6 2
1 x4 + O x6 ; f f 3 i i
2 x4 + f f 3 i i
i
12
2 x4 − f f 3 i i
(13)
13 1 (4) x5 + O x6 fi fi − fi fi 6 2
The smoothness indicator defined by (11) has currently been widely used in the WENO schemes, with which a rth-order ENO scheme can be converted to a (2r −1)th-order WENO scheme. However, to evaluate the explicit form of βl is essentially a very hard work, since it involves an integral calculation on the square of derivatives of the reconstruction polynomials. Although the explicit formula of βl for r = 4, 5, 6 and r = 7, 8, 9 have been presented by Balsara and Shu [1] and Gerolymos, Schal and Vallet [8], respectively, the works for evaluating the explicit form of βl are still tremendously burdensome especially for large r. Apart from the above mentioned works on the explicit formulae for the smoothness indicator, it is worthy of a particular mention that new expressions that are more compact than those in Balsara and Shu [1] and Jiang and Shu [12] have also been obtained recently in [3] and [4]. By casting the problem in a modal space of Hermite polynomials
q0 (ξ ) = 1,
q1 (ξ ) = ξ,
q2 (ξ ) = ξ 2 −
1 12
,
q3 (ξ ) = ξ 3 −
3 20
ξ,
(14)
Balsara et al. [3,4] have succeeded in deriving a series of more elegant and compact formulations for the smoothness indicators. For r = 3, these smoothness indicators can be expressed explicitly as
ISi = ( f xi )2 + where
13 3
( f xxi )2 (i = 0, 1, 2),
(15)
⎧ 1 1 ⎪ ⎪ f x0 = ( f i −2 − 4 f i −1 + 3 f i ); f xx0 = ( f i −2 − 2 f i −1 + f i ) ⎪ ⎪ 2 2 ⎪ ⎨ 1 1 f x1 = ( f i +1 − f i −1 ); f xx1 = ( f i −1 − 2 f i + f i +1 ) ⎪ 2 2 ⎪ ⎪ ⎪ ⎪ ⎩ f x2 = − 1 (3 f i − 4 f i +1 + f i +2 ); f xx2 = 1 ( f i − 2 f i +1 + f i +2 ) 2
(16)
2
Although, it can be seen that the above results have an equivalent form to the classic one suggested by Jiang and Shu [12], however, for r larger than 3, the explicit formulae presented in [3] and [4] are much more compact than the classic ones by Balsara and Shu [1]. For the sake of brevity, expressions of the mentioned smoothness indicator for r = 4 and 5 will not be presented in this paper and the interested readers who want to acquire the exact form are referred to Refs. [3] and [4]. 3.2. The new local smoothness indicator In this paper, we present a new succinct formula for the local smoothness indicators defined by
ηl =
r −1
(m) 2
Pl
=
m =1
r −1
m =1
2 (m) xm P i −r +1+l,r (xi ) ,
0 l r − 1,
(17)
(m)
where P l is a polynomial that approximates the mth-order derivative of the flux function f at point xi . P i ,r (x) is the Lagrangian interpolation polynomial for approximating the value of the flux function f at any point x based upon a r-point stencil S i ,r = {xi , . . . , xi +r −1 }
P i ,r (x) =
r −1 m =0
f i +m
r −1 l=0,l=m
x − xi +l xi +m − xi +l
,
(18)
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
333
(m)
and P i ,r (x) is the m-th derivative of P i ,r (x). (m) can be written explicitly as For r = 3, the polynomials P l
⎧ 1 3 (1 ) (1 ) ⎪ ⎪ P 0 = xP i −2,3 (xi ) = f i −2 − 2 f i −1 + f i , ⎪ ⎪ 2 2 ⎪ ⎨ 1 1 (1 ) (1 ) P 1 = xP i −1,3 (xi ) = − f i −1 + f i +1 , ⎪ 2 2 ⎪ ⎪ ⎪ (1 ) ⎪ ⎩ P = xP (1) (xi ) = − 3 f i + 2 f i +1 − 1 f i +2 2 i ,3 2 2 ⎧ (2 ) 2 (2 ) P = x P i −2,3 (xi ) = f i −2 − 2 f i −1 + f i , ⎪ ⎪ ⎨ 0 (2 ) (2 ) P 1 = x2 P i −1,3 (xi ) = f i −1 − 2 f i + f i +1 , ⎪ ⎪ ⎩ (2 ) (2 ) P 2 = x2 P i ,3 (xi ) = f i − 2 f i +1 + f i +2
Accordingly, the new smoothness indicators
⎧ ⎪ ⎪ η0 = ⎪ ⎪ ⎪ ⎨ η1 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ η2 =
1 4 1 4 1 4
(19)
(20)
ηl in terms of the cell averaged value of f can be expressed explicitly as
( f i −2 − 4 f i −1 + 3 f i )2 + ( f i −2 − 2 f i −1 + f i )2 , ( f i −1 − f i +1 )2 + ( f i −1 − 2 f i + f i +1 )2 ,
(21)
(3 f i − 4 f i +1 + f i +2 )2 + ( f i − 2 f i +1 + f i +2 )2
Using the Taylor series expansion, we find
⎧ ⎪ ⎪ (1 ) ⎪ P0 = ⎪ ⎪ ⎪ ⎨ (1 ) P1 = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ P (1 ) = 2 ⎧ (2 ) ⎪ ⎪ ⎨ P0 = (2 ) P1 = ⎪ ⎪ ⎩ (2 ) P2 =
(3 )
(1 )
fi
x −
fi
3
x3 +
1 (4 ) 4 f x + O x5 , 4 i
1 (3 ) 3 f i x + O x5 , 6 1 (3 ) 1 (4 ) (1 ) f i x − f i x3 − f i x4 + O x5 3 4 (1 )
fi
(22)
x +
x2 − f i(3) x3 + O x4 , (2 ) f i x2 + O x4 , (2 ) (3 ) f i x2 + f i x3 + O x4 (2 )
fi
Thus, the Taylor series expansions of
⎧ ( 1 ) 2 ⎪ (2 )2 ⎪ η = f x + fi − ⎪ 0 i ⎪ ⎪ ⎪ ⎪ ⎨ 2 η1 = f i(1) x + f i(2)2 + ⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ ⎪ ⎩ η2 = f ( 1 ) x + f ( 2 ) 2 − i
i
ηl read
2 (1 ) (3 ) x4 + fi f 3 i
(23)
1 (1 ) (4 ) (2 ) (3 ) x5 + O x6 fi fi − 2 fi fi 2
1 (1 ) (3 ) x4 + O x6 fi f 3 i
2 (1 ) (3 ) x4 − fi f 3 i
ηl are 1 + O x2 ,
(24)
1 (1 ) (4 ) (2 ) (3 ) x5 + O x6 fi fi − 2 fi fi 2
Obviously, the truncation errors of
2
ηl = f i(1) x
(25)
showing a close resemblance to the asymptotic features of βl . Nevertheless, the computation for ηl is much simplified since it needs only a simple computation for the Lagrangian interpolation polynomials. Furthermore, higher order global smoothness indicators can be easily devised with ηl , making the associated WENO scheme less dissipative than the classic ones. 3.3. New higher order global smoothness indicators (1)
τ5
The novel idea of the WENO-Z scheme introduced in [5] is the construction of the higher order global smoothness indicator τ5 from the classic local smoothness indicator βl . In this paper, we find that a new 5th-order global smoothness
334
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
indicator can be acquired from a similar operation on the new local smoothness indicator smoothness indicator can be calculated by
ηl . As in (9), the new global
τ5 = |η0 − η2 |
(26)
It is straightforward to see from (24) that the truncation error of the new
τ5 is
(1 ) (4 ) f f − 4 f (2) f (3) x5 + O x6 , i
i
i
(27)
i
which is equivalent to that introduced in [5]. In the latter part of this paper, we will show by numerical tests that the results obtained by the WENO-Z scheme with βl and ηl and the corresponding τ5 have little differences. (2)
τ6
As argued previously, it is found in this paper that the order of τ can in fact be raised to higher than 5th-order without using information outside the given large 5-point stencil S i −2,5 = {xi −2 , . . . , xi +2 }. We will not present the method for devising the 6th-order τ from βl , since it is analogous to what will be introduced in the following by using the new local smoothness indicator ηl . To build this higher order τ , we should first introduce a reference value for evaluating the smoothness of the large 5-point stencil S i −2,5 , which is written by
η5 =
2
2 (m) xm P i −2,5 (xi )
m =1
=
1
144
( f i −2 − 8 f i −1 + 8 f i +1 − f i +2 )2 + ( f i −2 − 16 f i −1 + 30 f i − 16 f i +1 + f i +2 )2
From using the Taylor series expansion, we see that
2
η5 = f i(1) x + f i(2) x2
2
(1 ) (5 )
−
fi
fi
15
6
(28)
η5 can be expanded as
x6 + O x7
(29)
On the other hand, one can note from (24) that the following combination of
4 η1 + η0 + η2
ηl can be expanded as
2 2 = f i(1) x + f i(2) x2 + O x6
Therefore, a new 6th-order global smoothness indicator
4 η1 + η0 + η2 = O x6 τ6 = η5 − 6
(30)
τ6 can be acquired by the following ways (31)
Comparing with τ5 , the rate of convergence of WENO-Z with τ6 can recover 5th-order at a first order critical point with no need of increasing the power p in (8) from 1 to 2. Also, it will be shown that results of τ6 are less dissipative than those of τ5 . (3)
τ8
Other than the above introduced τ5 and τ6 , this paper finds that the truncation error of τ can be raised up to 8 th-order by a quite different approach. The new 8th-order global smoothness indicator τ8−1 is defined by
τ8−1 = P 0(1) − P 2(1) P 0(2) + P 2(2) − 2P 1(2)
(32)
From the Taylor series expansion, one can find that the truncation error of
1 τ8−1 = f i(4) x4 + O x5 O x4 = O x8
τ8−1 is
(33)
2
It is worthy of note that there is another way to construct the 8th-order
2 2 τ8−2 = P (1) − P (1) + P (2) + P (2) − 2P (2) 0
2
0
2
1
τ , as defined by (34)
It is easy to verify that τ8−2 is also of 8th-order, while, in contrast to τ8−1 , the results of τ8−2 are relatively more dissipative than those of τ8−1 . This can be shown in the latter numerical experiments. Furthermore, comparing with τ5 and τ6 , the benefit of τ8 is that the rate of convergence of WENO-Z with τ8 can recover that of the formal 5th-order, even at a second order critical point. This will be demonstrated in the following section.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
335
Table 1a Rate of convergence at a first order critical point (k = 2) for WENO-JS and WENO-Z schemes. N
WENO-Z (p = 1)
WENO-JS
23 24 25 26 27 28 29
WENO-Z (p = 2)
Error
Order
Error
Order
Error
Order
5.8085E–02 6.1571E–03 4.8105E–04 3.2076E–05 1.1798E–06 1.2947E–08 1.1654E–10
−
1.6448E–02 9.9119E–04 5.5215E–05 3.1169E–06 1.8038E–07 1.0713E–08 6.4960E–10
– 4.0526 4.1660 4.1469 4.1110 4.0736 4.0437
5.4640E–02 4.0202E–03 1.3079E–04 3.0220E–06 6.5728E–08 1.5353E–09 3.9560E–11
– 3.7646 4.9419 5.4356 5.5229 5.4199 5.2783
3.2378 3.6780 3.9066 4.7649 6.5098 6.7956
Table 1b Rate of convergence at a first order critical point (k = 2) for WENO-Zη(τ5 ) schemes. WENO-Zη(τ5 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ5 ) ( p = 2)
Error
Order
Error
Order
1.6237E–02 9.8866E–04 5.5315E–05 3.1259E–06 1.8086E–07 1.0737E–08 6.5090E–10
– 4.0377 4.1597 4.1453 4.1113 4.0742 4.0440
5.4380E–02 4.0045E–03 1.3009E–04 2.9932E–06 6.4699E–08 1.5011E–09 3.8461E–11
– 3.7634 4.9440 5.4417 5.5318 5.4297 5.2865
Table 1c Rate of convergence at a first order critical point (k = 2) for WENO-Zη(τ6 ) schemes. WENO-Zη(τ6 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ6 ) ( p = 2)
Error
Order
Error
Order
6.1221E–03 1.4197E–04 3.2037E–06 7.9428E–08 2.1549E–09 6.2224E–11 1.9593E–12
– 5.4304 5.4697 5.3339 5.2040 5.1140 4.9891
1.7041E–02 7.2599E–05 1.7949E–06 5.7051E–08 1.7984E–09 5.6529E–11 1.8119E–12
– 7.8748 4.9532 4.9755 4.9875 4.9916 4.9634
3.4. Convergence at critical points Critical points have become a focus of discussion since in [11] it was shown that the standard value of the parameter
ε = 10−6 in fact hided the loss of accuracy of WENO-JS which achieves only third order at critical points of smooth solutions.
WENO-M was presented as a fix to this situation, for it used a mapping that corrected the weights of WENO-JS and recovered the formal fifth-order at critical points even when using very small values of ε . Furthermore, it was demonstrated in [5] that the WENO-Z scheme is fourth-order accurate at first order critical points of a smooth solution and the fifth-order convergence can be recovered by increasing parameter p in (8) from p = 1 to p = 2. However, this modification increases the dissipation of the WENO-Z scheme. We now perform a numerical experiment to compare the behavior of the classic WENO-JS, WENO-Z and the new WENO-Zη in the presence of smooth solutions containing critical points. Unless indicated otherwise, the values of parameter ε are taken as ε = 10−6 for the WENO-JS scheme and ε = 10−40 for the others. Consider the following test function
f (x) = xk exp(x),
x ∈ [−1, 1]
(35)
in which its first k − 1 derivatives f ( j ) (0) = 0, j = 0, 1, . . . , k − 1. That is, this function has a critical point of order ncp = k − 1 at x = 0. Tables 1 and 2 show the L ∞ convergence behaviors for these schemes at the first and second order critical points, i.e., k = 2 and 3, respectively. We use the fixed value of p = 2 for WENO-JS since its rate of convergence does not change when p varies. It should be mentioned particularly that in this subsection only the convergence of WENO-Zη(τ8−2 ) is considered and referred to as WENO-Zη(τ8 ). As shown in Table 1a and Table 2a, increasing value of p from 1 to 2 allows a speed up of convergence rate of WENO-Z at a first order critical point, but this implementation loses its effect at a second order critical point. What in Tables 1b and 2b shows that WENO-Zη(τ5 ) takes an equivalent performance as WENO-Z. In Tables 1c and 2c, it can be seen that WENO-Zη(τ6 ) recovers the fifth-order convergence rate at a first critical point, regardless of value of p. However, at a second order critical point, the rate of convergence of WENO-Zη(τ6 ) degrades to
336
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Table 1d Rate of convergence at a first order critical point (k = 2) for WENO-Zη(τ8 ) schemes. WENO-Zη(τ8 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ8 ) ( p = 2)
Error
Order
Error
Order
1.7427E–03 5.5593E–05 1.7949E–06 5.7051E–08 1.7984E–09 5.6529E–11 1.8598E–12
– 4.9703 4.9529 4.9755 4.9875 4.9916 4.9258
1.6775E–03 5.5604E–05 1.7949E–06 5.7051E–08 1.7984E–09 5.6529E–11 1.8119E–12
– 4.9150 4.9532 4.9755 4.9875 4.9916 4.9634
Table 2a Rate of convergence at a second order critical point (k = 3) for WENO-JS and WENO-Z schemes. N
23 24 25 26 27 28 29
WENO-Z (p = 1)
WENO-JS
WENO-Z (p = 2)
Error
Order
Error
Order
Error
Order
1.6430E–01 3.5553E–02 6.3056E–03 5.0357E–05 8.2116E–07 1.4425E–08 2.7234E–10
−
1.4191E–01 3.0455E–02 7.0249E–03 1.6684E–03 3.9888E–04 9.4831E–05 2.2365E–05
– 2.2202 2.1161 2.0740 2.0644 2.0725 2.0841
1.6429E–01 3.5573E–02 8.3065E–03 2.0061E–03 4.9003E–04 1.1848E–04 2.7570E–05
– 2.2074 2.0985 2.0498 2.0335 2.0482 2.1035
2.2083 2.4953 6.9683 5.9384 5.8310 5.7270
Table 2b Rate of convergence at a second order critical point (k = 3) for WENO-Zη(τ5 ) schemes. N
23 24 25 26 27 28 29
WENO-Zη(τ5 ) ( p = 1)
WENO-Z
η(τ5 ) ( p = 2)
Error
Order
Error
Order
1.4080E–01 3.0176E–02 6.9518E–03 1.6488E–03 3.9351E–04 9.3377E–05 2.1994E–05
– 2.2222 2.1179 2.0760 2.0669 2.0753 2.0860
1.6423E–01 3.5562E–02 8.3032E–03 2.0048E–03 4.8918E–04 1.1790E–04 2.7301E–05
– 2.2073 2.0986 2.0502 2.0350 2.0528 2.1105
Table 2c Rate of convergence at a second order critical point (k = 3) for WENO-Zη(τ6 ) schemes. N 23 24 25 26 27 28 29
WENO-Zη(τ6 ) ( p = 1) Error
Order
WENO-Z Error
η(τ6 ) ( p = 2)
1.2804E–01 2.5667E–02 5.6184E–03 1.2991E–03 3.1093E–04 7.5945E–05 1.8759E–05
– 2.3186 2.1917 2.1126 2.0629 2.0336 2.0174
1.6377E–01 3.5277E–02 8.1766E–03 1.9626E–03 4.7987E–04 1.1856E–04 2.9457E–05
Order – 2.2149 2.1092 2.0587 2.0321 2.0170 2.0089
second order, as do WENO-Z and WENO-Zη(τ5 ). Tables 1d and 2d illustrate that WENO-Zη(τ8 ) recovers the fifth-order convergence rate at a first critical point. At a second critical point, rate of convergence of WENO-Zη(τ8 ) can achieve 5th-order when value of p is increased from 1 to 2. 4. Numerical examples In this section, we compare the numerical performance of WENO-η and WENO-Zη with the classic WENO-JS and WENOZ. The numerical presentation of this section starts with the solution of the scalar advection equation and followed by the solution of the one-dimensional Euler equations with Riemann initial value problems such as the Lax, Sod, shock-density wave interaction and the interacting blast problems. Finally, the two-dimensional problems on the Richtmyer–Meshkov instability, shock-vortex interaction, Rayleigh–Taylor instability, double-Mach shock reflection and the forward facing step problems are simulated. The time dependent problems are all solved via the 3rd-order Runge–Kutta TVD method [20] with CFL = 0.5.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
337
Table 2d Rate of convergence at a second order critical point (k = 3) for WENO-Zη(τ8 ) schemes. WENO-Zη(τ8 ) ( p = 1)
N
23 24 25 26 27 28 29
WENO-Zη(τ8 ) ( p = 2)
Error
Order
Error
Order
9.3881E–02 9.5484E–03 1.4256E–03 1.1070E–04 7.4017E–06 4.6962E–07 2.9410E–08
– 3.2975 2.7437 3.6868 3.9027 3.9783 3.9971
1.5496E–01 1.3660E–02 1.6515E–03 3.5043E–05 5.5463E–07 8.6262E–09 1.3503E–10
– 3.5039 3.0481 5.5585 5.9815 6.0067 5.9974
Fig. 1. Numerical solutions of the linear advection equation (36) with N = 27 at t = 6: (a) presents the comparison between the exact solution, the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
4.1. 1D scalar test problem Consider the one-dimensional linear advection equation
∂C ∂C + = 0, ∂t ∂x
x ∈ [−1, 1]
with an initial condition given by
(36)
338
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 2. Numerical solutions of the Lax problem with N = 200 at t = 1.3: (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
C (x, 0) =
1,
x ∈ [−0.5, 0.5),
0,
otherwise
(37)
The advection equation (36) was solved until the final time t = 6 on the interval [−1, 1] with 27 grid points and the periodic boundary conditions. Figs. 1(a) and (a-1) compare the results between the exact solution, the WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ) schemes. It is seen that the new WENO-η and WENO-Zη(τ5 ) display nearly equal performance to the WENO-JS and WENO-Z, respectively. Furthermore, the results of WENO-Zη with different global smoothness indicators τ are compared and presented in Figs. 1(b) and (b-1). It can be observed that among these schemes the WENO-Zη(τ5 ) and WENO-Zη(τ8−1 ) achieve, respectively, the largest and least dissipations near the discontinuities. 4.2. 1D Euler systems In this section, we present numerical experiment with the one-dimensional Euler equations for gas dynamics in strong conservation form
∂U ∂F + = 0, ∂t ∂x where
(38)
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
339
Fig. 3. Numerical solutions of the Sod problem with N = 200 at t = 0.2: (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
U = (ρ , ρ u , E ) T , F=
ρ u , ρ u 2 + P , ( E + P )u
T
(39)
The equation of state is given by
P = (γ − 1) E −
1 2
ρ u 2 , γ = 1.4,
(40)
where ρ , u, P and E are the density, velocity, pressure and total energy, respectively. These equations are solved by following the general WENO methodology [12]. Specifically, the Roe approximation is used for the characteristic decomposition at the cell faces and the Lax–Friedrichs formulation is used for the numerical fluxes. The numerical experiments are conducted using 200 grid points. (I). Riemann initial value problem of Lax For the Lax problem, the initial condition is
(ρ , U , p ) =
(0.445, 0.698, 0.3528), −5 x < 0, 0 x 5, (0.5, 0, 0.571),
(41)
340
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 4. Numerical solutions of the shock-density interaction with N = 200 at t = 1.8: (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) presents the comparison between the exact solution and the numerical results of WENO-Zη(τ5 ), WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ); (a-1) and (b-1) are the local enlarged plots of (a) and (b), respectively.
and the final time is 1.3. The simulated density of the Lax problem is shown in Figs. 2. We refer to the solution computed by the WENO-JS with N = 4000 grid points as the ‘exact’ solution. From Figs. 2(a) and (a-1) we see that the new WENO-η and WENO-Zη(τ5 ) schemes show equivalent performance as those of WENO-JS and WENO-Z, respectively. Figs. 2(b) and (b-1) show that the WENO-Zη(τ8−1 ) presents the least dissipation near the discontinuities compared to WENO-Zη(τ5 ), WENO-Zη(τ6 ) and WENO-Zη(τ8−2 ). (II). Riemann initial value problem of Sod For the Sod problem, the initial condition is
(ρ , U , p ) =
x 0, (1.0, 0, 1.0), (0.125, 0, 0.1), x > 0,
(42)
and the final time is 0.2. The simulated density of the Sod problem is shown in Figs. 3. As before, the performance of WENO-η and WENO-Z η(τ5 ) are nearly the same as those of WENO-JS and WENO-Z, respectively, while, the WENO-Zη(τ8−1 ) shows the least dissipation than WENO-Zη(τ5 ), WENO-Zη(τ6 ) and WENO-Zη(τ8−2 ).
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
341
Fig. 5. Numerical solutions of the interacting blast waves problem with N = 200 at t = 0.038; (a) presents the comparison between the exact solution and the numerical results of WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ); (b) gives the local enlarged plot of (a).
(III). Shock-density wave interaction The initial condition for the shock-density wave interaction problem [21] is given by
(ρ , U , p ) =
(3.857143, 2.629369, 10.3333), x < −4, x −4 (1 + 0.2 sin(5x), 0, 1.0),
(43)
This case represents a Mach 3 shock wave interacting a sine entropy wave. The results at time t = 1.8 with mesh size of N = 200 are plotted in Figs. 4. As before, WENO-η and WENO-Zη(τ5 ) show good agreement with WENO-JS and WENO-Z, respectively. Moreover, the results of WENO-Zη(τ6 ) and WENO-Zη(τ8−1 ) are much less dissipative than those of the other schemes. (IV). Interacting blast wave problem The initial condition for the interacting blast wave example [28] is
(ρ , U , p ) =
(1, 0, 1000), 0 x < 0.1, (1, 0, 0.01), 0.1 x < 0.9, (1, 0, 100), 0.9 x 1,
(44)
342
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 6. Initial location of the interface in the RM instability.
and the final time is 0.038. Since this problem is much more stringent than those previous ones, only results from WENO-JS, WENO-η , WENO-Z and WENO-Zη(τ5 ) are presented. Simulations under the use of WENO-Zη schemes equipped with higher order global smoothness indicators are found to be unstable, which may be largely ascribed to the less dissipation of these schemes. From Figs. 5(a) and (b) it can be seen that the new WENO-Zη(τ5 ) scheme is as good as WENO-Z, both of which are less dissipative than WENO-JS and WENO-η . 4.3. 2D Euler systems In this subsection, we extend the new WENO-η and WENO-Zη schemes to solve the 2D compressible system of the form:
∂U ∂F ∂G + + = 0, ∂t ∂x ∂x
(45)
where
U = (ρ , ρ u , ρ v , E ) T , F= G=
T
ρ u , ρ u 2 + P , ρ uv , ( E + P )u ,
ρ v , ρ uv , ρ v 2 + P , ( E + P ) v
T
(46)
The equation of state is given by
1 2 2 , P = (γ − 1) E − ρ u + v 2
γ = 1.4
(47)
Here ρ , u, v, P and E are the density, components of velocity in the x and y coordinate directions, pressure and total energy, respectively. The Steger–Warming flux vector splitting method is used for the 2D problems in this paper. (I) Richtmyer–Meshkov Instability (RMI) In this test, we consider the standard RMI unstable interface problem [15,17,24]. The initial condition in the tube is composed of an interface separating two fluids of different densities and a shock wave approaching the interface. The shock tube is [0, 1.5] × [0, 1] in the 2D domain. It is well known that this interface becomes unstable after the passage of a shock wave, irrespective of the side of the fluid that the shock is incident upon. Here we consider a single mode perturbation of an interface between two fluids, and the initial location of the interface (as shown in Fig. 6) is represented by x = x0 + a0 cos(2π y ), where x0 = 0.2 is the location of the unperturbed interface, a0 = 0.05 is the amplitude of the perturbation. We use the ideal gas equation of state for air and SF6 with γ = 1.4 and γ = 1.093, respectively. We choose the density ratio equal to 6 with a corresponding Atwood number
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 7. Density profiles of the RMI driven by a planar shock: (a) WENO-JS, (b) WENO-Z, (c) WENO-Zη(τ6 ), (d) WENO-Zη(τ8−1 ), (e) WENO-Zη(τ8−2 ).
343
344
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 8. The pressure contour of the shock/vortex interaction by WENO-Zη(τ6 ).
Fig. 9(a). Pressure distributions of the shock/vortex interaction on central line y = 0.5.
0.71. To trigger the instability, a planar Mach 1.2 shock wave propagates from the left to the interface. This gives us one example that the interface is accelerated by a shock wave coming from the light-fluid to the heavy-fluid region, and the resulting wave pattern after the interaction would consist of a transmitting shock wave, an interface, and a reflected rarefaction wave. In Figs. 7 we compare the numerical performances of the WENO-JS, WENO-Z, WENO-Zη(τ6 ), WENO-Zη(τ8−1 ) and WENO-Zη(τ8−2 ) at t = 6 with a mesh system 241 × 161 is used. An examine of these results reveals that WENO-Zη(τ6 ) and WENO-Zη(τ8−1 ) display better resolution of the fine structure near the interface than the other WENO schemes. (II) Two-dimensional shock vortex interaction In this subsection the two-dimensional shock vortex interaction problem taken from [12] is solved. This problem describes the interaction between a stationary shock and a vortex. The computational domain is taken to be [0, 2] × [0, 1]. A stationary Mach 1.1 shock is positioned at x = 0.5 and normal to the x-axis. Its left state is (ρ , u , v , p ) = (1, 1.1γ 1/2 , 0, 1). A small vortex is superimposed to the flow on the left of the shock and is centered at (xc , y c ) = (0.25, 0.5). The vortex is described as a perturbation to the velocity (u , v ), temperature T = p /ρ , and entropy S = ln( p /ρ γ ) of the mean flow and the parameters are taken the same as in [12,19]. The time step is taken as follows [16]
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
345
Fig. 9(b). Locally enlarged plot of Fig. 9(a).
Fig. 10. Density contours for the Rayleigh–Taylor instability: (a) WENO-JS; (b) WENO-Z; (c) WENO-Zη(τ6 ); (d) WENO-Zη(τ8−1 ); (e) WENO-Zη(τ8−2 ).
t x t y , t x + t y x t x = , maxi , j (|u i , j | + c i , j )
t = σ
t y =
y , maxi , j (| v i , j | + c i , j )
(48)
where c is the speed of sound and σ = 0.5 is the CFL number. Fig. 8 is the pressure contours at t = 0.6 calculated by the WENO-Zη(τ6 ) scheme. Figs. 9(a) and 9(b) are the comparisons of the pressure along the center line of y = 0.5. In order to show the accuracy of the new scheme, the results obtained by the WENO-Z scheme with a refined mesh of 2001 × 401 is given as the “exact” solution. The results with the coarse mesh of 251 × 101 are compared. Figs. 9(a) and 9(b) show that the WENO-Zη(τ6 ) has an equivalent effect to the WENO-Zη(τ8−1 ), both of which have a sharper shock profile than WENO-JS, WENO-Z and WENO-Zη(τ8−2 ) due to lower numerical dissipations.
346
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 11. Density contours of the double-Mach shock reflection as computed by the WENO-η scheme at time t = 0.2. The mesh resolution is 800 × 200 zone and thirty contours were fit between a range of 1.3 and 22.
(III) Two-dimensional Rayleigh–Taylor instability
This problem has been simulated by many authors [18,25,26] to test the numerical dissipation of various high order shock capturing schemes. The computational domain is [0, 0.25] × [0, 1], and the initial conditions are
⎧ ⎪ ⎨ (2, 0, −0.025 53 ρp cos(8π x), 2 y + 1), if 0 y < 12 (ρ , u , v , p ) = ⎪ ⎩ (1, 0, −0.025 5 p cos(8π x), y + 3 ), if 1 y < 1 3ρ 2 2
(49)
This problem represents the interface instability between fluids with different densities when an acceleration is directed from heavy fluid to light one. The gravitational effect is introduced by adding ρ and ρ v to the right hand side of the y-momentum and the energy equations, respectively. Reflective boundary conditions are imposed for the left and right boundaries, and (ρ , u , v , p ) = (1, 0, 0, 2.5) and (ρ , u , v , p ) = (2, 0, 0, 1) are set as top and bottom boundary conditions, respectively. The simulation time is t = 1.95. As analyzed in [18,25], because the inviscid Euler equation is solved, the details of the complex structures due to the physical instability are related to the specific form of numerical viscosity of the scheme. From the density contours with meshes of 121× 481 plotted in Fig. 10, it can be seen that the WENO-Zη(τ8−1 ) scheme obtains more complex structures than the other schemes.
(IV) Double Mach reflection of a strong shock
This problem has been carefully examined by Woodward and Colella [28]. We use the same setup as the authors did to test the problem where a Mach 10 shock hits a reflecting wall at the bottom of the domain given by [0, 4] × [0, 1]. The angle between the shock and the horizontal axis is 60◦ . The equations are the two-dimensional Euler equations (γ = 1.4) and initial conditions are
(ρ , u , v , p ) =
(8, 8.25 cos π6 , −8.25 sin π6 , 116.5), x < x0 + (1.4, 0, 0, 1.0),
√y
3
y x x0 + √
(50)
3
with x0 = 1/6. Boundary conditions at x = 0 are inflow, with post-shock values as above, and the right boundary is set to be outflow boundary. At y = 0, reflecting boundary conditions are applied to the interval [x0 , 4]. The top boundary condition imposes the exact motion of a Mach 10 shock in the flow variables. Fig. 11 shows the density variable at t = 0.2 in [0, 3] × [0, 1] with a mesh resolution of 800 × 200 zones as computed by the WENO-η scheme. In Figs. 12 we display the region around the double Mach stems in order to observe more clearly the
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
347
Fig. 12. Density contours of the double-Mach shock reflection as computed by the WENO-JS (a); WENO-η (b); WENO-Z (c); WENO-Zη(τ5 ) (d); WENO-Zη(τ6 ) (e); WENO-Zη(τ8−1 ) (f); and WENO-Zη(τ8−2 ) schemes (g). The mesh resolution is 800× 200 zone and thirty contours were fit between a range of 1.3 and 22.
numerical solutions of the different WENO schemes. In general, the global structure of the solution is very similar across different schemes. However, the dissipation of the various schemes can be distinguished by the number of small vortices generated along the slip line and the wall jet behind the lower half of the right moving shock. As is shown clearly in Figs. 12, the WENO-Zη(τ6 ) and WENO-Zη(τ8−1 ) exhibit the least dissipation, while, the WENO-JS and WENO-η are the most dissipative. Fig. 13 shows the density variable under different mesh resolutions which are 400 × 100, 800 × 200 and 1600 × 400, respectively. It is clear that the small vortices along the slip line are captured and become gradually obvious with the increase of the mesh resolutions. (V) Forward facing step problem This test problem was first presented by Woodward and Colella [28]. Recently, Balsara et al. [3] carried out a resolution study using the fourth order ADER-WENO scheme and showed that increasing the resolution enabled them to capture important details such as the roll-up of the vortex sheet via Kelvin–Helmholtz instability. In this paper, we also successfully captured the roll-up of the vortex sheet not only from increasing the grid resolution but also from using the new less
348
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 12. (continued)
dissipative WENO-Zη schemes. Furthermore, the robustness of the new WENO schemes gets substantially demonstrated from the success of simulating this severely stringent problem. The setup of the problem is the following: the wind tunnel spans a domain of [0, 3] × [0, 1]. A forward-facing step is located at (0.6, 0.2). The problem is initialized by a right-going Mach 3 flow with a density of 1. 4 and a pressure of unity. Reflective boundary conditions are applied along the walls of the tunnel. Inflow and outflow boundary conditions are applied, respectively, at the left and right boundaries. The simulation was run until a time of 4 and the ratio of specific heats is given by 1.4. Figs. 14 and 15 show the density at the final time at mesh resolutions of 600 × 200 and 1200× 400 zones, respectively. In all of these simulations, five 5th-order WENO schemes are used, which, in turns, are (a) WENO-JS; (b) WENO-η ; (c) WENOZ; (d) WENO-Zη(τ5 ) and (e) WENO-Zη(τ8−2 ). All the shocks are found to have sharp profiles and captured properly on the computing grid. From comparing Figs. 14 and Figs. 15, we see that with the increase of the mesh resolution from 600 × 200 to 1200 × 400, the vortex sheets are observed more evidently that emanates from the Mach stem. Moreover, the roll-up of the vortex sheet becomes the most intensely when the WENO-Zη(τ8−2 ) scheme is used. This presents a powerful evidence for the advantage of the new WENO-Zη scheme that has the less dissipation and hence has a better resolution in capturing details of the flowing structures.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
349
Fig. 13. Density contours of the double-Mach shock reflection as computed by the WENO-JS and WENO-η schemes. The mesh resolutions are 400 × 100 in (a) and (a-1); 800 × 200 in (b) and (b-1); and 1600 × 400 in (c) and (c-1). Thirty equally spaced contours are shown ranging from 1.3 to 22.
350
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 14. Density contours of the forward facing step problem as computed by the WENO-JS (a); WENO-η (b); WENO-Z (c); WENO-Zη(τ5 ) (d); and WENO-Zη(τ8−2 ) schemes (e). The resolution is 600 × 200 zones and time t = 4. Thirty equally spaced contours are shown ranging from 0.1 to 6.2.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
351
Fig. 14. (continued)
5. Conclusions
In this paper, we introduce a new local smoothness indicator for the WENO scheme which is defined based on the Lagrangian interpolation polynomial and is more compactable than the classical ones. The WENO scheme using the new local smoothness indicator, named as WENO-η , has an equal ENO property as the classic WENO-JS scheme in dealing with the shock-capturing problem. Furthermore, several higher order global smoothness indicators are devised, ensuring the resulted new WENO-Zη scheme can provide the fifth convergence order in smooth regions, especially at the first and second order critical points. Numerical experiments are run to show that the WENO-Zη scheme with 5th-order global smoothness indicator τ5 achieves a similar behavior as WENO-Z and can induce less dissipation near the discontinuities of the solution if higher order global smoothness indicators τ6 , τ8−1 and τ8−2 are used.
Acknowledgements
This work is supported by the NSFC (11172299), 973 Program (2012CB224806), CAS Program for Cross & Cooperative Team of the Science & Technology Innovation, and the National Natural Science Foundation of China (21376243). The helpful comments from the anonymous reviewers are also appreciated.
352
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
Fig. 15. Density contours of the forward facing step problem as computed by the WENO-JS (a); WENO-η (b); WENO-Z (c); WENO-Zη(τ5 ) (d); and WENO-Zη(τ8−2 ) schemes (e). The resolution is 1200 × 400 zones and time t = 4. Thirty contours were fit between a range of 0.1 and 6.2.
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
353
Fig. 15. (continued)
References [1] D. Balsara, C.W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy, J. Comput. Phys. 160 (2000) 405–452. [2] D.S. Balsara, Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics, J. Comput. Phys. 228 (2009) 5040–5056. [3] D.S. Balsara, T. Rumpf, M. Dumbser, C.-D. Munz, Efficient, high-accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics, J. Comput. Phys. 228 (2009) 2480–2516. [4] D.S. Balsara, M. Dumbser, C. Meyer, H. Du, Z. Xu, Efficient implementation of ADER schemes for Euler and magnetohydrodynamic flow on structured meshes – comparison with Runge–Kutta methods, J. Comput. Phys. 235 (2013) 934–969. [5] R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (2008) 3101–3211. [6] M. Castro, B. Costa, W.S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys. 230 (2011) 1766–1792. [7] M. Dumbser, M. Käser, Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems, J. Comput. Phys. 221 (2007) 693–723. [8] G.A. Gerolymos, D. Schal, I. Vallet, Very high order WENO schemes, J. Comput. Phys. 228 (2009) 8481–8524. [9] Y. Ha, C.H. Kim, Y.J. Lee, J. Yoon, An improved weighted essentially non-oscillatory scheme with a new smoothness indicator, J. Comput. Phys. 232 (2013) 68–86. [10] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order essentially non-oscillatory schemes, III, J. Comput. Phys. 71 (1987) 231–303. [11] A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points, J. Comput. Phys. 207 (2005) 542–567. [12] G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [13] D. Levy, G. Puppo, G. Russo, Compact central WENO schemes for multidimensional conservation laws, SIAM J. Sci. Comput. 22 (2000) 656–672. [14] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [15] K.O. Mikaelian, Rayleigh–Taylor and Richtmyer–Meshkov instabilities and mixing in stratified cylindrical shells, Phys. Fluids 17 (2005) 094105. [16] S. Pirozzoli, Conservative hybrid compact-WENO schemes for shock–turbulence interaction, J. Comput. Phys. 178 (2002) 81–117.
354
[17] [18] [19] [20] [21] [22] [23]
[24] [25] [26] [27] [28]
P. Fan et al. / Journal of Computational Physics 269 (2014) 329–354
R.D. Richtmyer, Taylor instability in shock acceleration of compressible fluids, Commun. Pure Appl. Math. 8 (1960) 297–319. J. Shi, Y.T. Zhang, C.W. Shu, Resolution of high order WENO schemes for complicated flow structures, J. Comput. Phys. 186 (2003) 690–696. Y.Q. Shen, G.C. Zha, Generalized finite compact difference scheme for shock/complex flowfield interaction, J. Comput. Phys. 230 (2011) 4419–4436. C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78. C.W. Shu, High order weighted essentially non-oscillatory schemes for convection dominated problems, SIAM Rev. 51 (2009) 82–126. C.W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: B. Cockburn, C. Johnson, C.-W. Shu, E. Tadmor (Eds.), Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, in: Lecture Notes in Mathematics, Springer-Verlag, Berlin/New York, 1998, pp. 325–432. B.L. Tian, W.D. Shen, S. Jiang, S.H. Wang, Y. Liu, A global arbitrary Lagrangian–Eulerian method for stratified Richtmyer–Meshkov instability, Comput. Fluids 46 (2011) 113–121. Z.F. Xu, C.-W. Shu, Anti-diffusive flux corrections for high order finite difference WENO schemes, J. Comput. Phys. 205 (2005) 458–485. Y.N. Young, H. Tufo, A. Dubey, R. Rosner, On the miscible Rayleigh–Taylor instability: two and three dimensions, J. Fluid Mech. 447 (2001) 377–408. Y.-T. Zhang, C.-W. Shu, Third order WENO scheme on three dimensional tetrahedral meshes, Commun. Comput. Phys. 5 (2009) 836–848. P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115–173.
