Next Article in Journal
Cluster-Based Aircraft Fuel Estimation Model for Effective and Efficient Fuel Budgeting on New Routes
Next Article in Special Issue
Trajectory Optimization with Complex Obstacle Avoidance Constraints via Homotopy Network Sequential Convex Programming
Previous Article in Journal
Research on Trajectory Prediction of a High-Altitude Zero-Pressure Balloon System to Assist Rapid Recovery
Previous Article in Special Issue
Framework for Estimating Performance and Associated Uncertainty for Modified Aircraft Configurations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A New Mapped WENO Method for Hyperbolic Problems

1
Department of Mechanical, Materials and Aerospace Engineering, University of Liverpool, Liverpool L69 3BX, UK
2
Department of Aerospace Engineering, University of Bristol, Bristol BS8 1TR, UK
3
School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore
*
Author to whom correspondence should be addressed.
Aerospace 2022, 9(10), 623; https://doi.org/10.3390/aerospace9100623
Submission received: 11 July 2022 / Revised: 9 October 2022 / Accepted: 11 October 2022 / Published: 19 October 2022
(This article belongs to the Special Issue Advances in Aerospace Sciences and Technology III)

Abstract

:
In this study, a new family of rational mapping functions gRM(ω;k,m,s) is introduced for seventh-order WENO schemes. gRM is a more general family of mapping functions, which includes other mapping functions such as gM and gIM as special cases. The mapped WENO scheme WENO-IM(2,0.1), which uses gIM, performs excellently at fifth order but rather poorly at seventh order. The reason for this loss of accuracy was found to be the over-amplification of very small weights by the mapping process, which can be traced back to the large slope of gIM at ω = 0. For m > 1, gRM can be designed to have a unit slope at ω = 0, which will preserve small weights with little to no amplification. It has been demonstrated through several one-dimensional linear advection test cases that the mapped WENO scheme WENO-RM(6,3,2 × 103), which uses the mapping function gRM(ω;6,3,2 × 103), outperforms both WENO-M and WENO-IM(2,0.1) at seventh order. The proposed scheme also performs better at a number of one-dimensional inviscid gas flow problems compared to other popular WENO schemes such as the WENO-Z scheme.

1. Introduction

Hyperbolic problems arise in many scientific and engineering applications. The Euler equations governing the dynamics of an inviscid fluid is a notable example of a non-linear hyperbolic problem which admits discontinuities to develop and persist in the solutions. These discontinuities, such as shock waves and contact surfaces, are key features in high-speed flows and are regularly encountered in aerospace engineering applications [1]. Needless to say, it is imperative to capture these discontinuities accurately. In addition to capturing discontinuities, satisfactory accuracy remains a requirement for smooth flow regions to resolve the finer-scale turbulence as well as pressure and velocity fluctuations, which are essential for the prediction of unsteady aerodynamic behaviour and aerodynamically generated noise [2]. As a result, the development of high-order schemes has attracted significant attention from both aerospace research and engineering communities over the past decade. Conventional high-order schemes, which perform excellently in the smooth regions of the solution, tend to produce spurious oscillations in the vicinity of discontinuities. Harten et al. [3] designed the essentially non-oscillatory (ENO) schemes to overcome this problem by choosing the smoothest stencil to reconstruct the solution, thereby effectively avoiding reconstruction across a discontinuity. Later, Liu et al. [4] proposed the weighted ENO (WENO) scheme in which a weighted reconstruction is obtained from all available stencils, with the weight assigned to each stencil being inversely proportional to the smoothness of the reconstruction within that stencil. In addition to remaining (essentially) non-oscillatory near discontinuities, the WENO method on a (2r − 1)-cell stencil was able to attain the (r +1)th order of convergence in smooth regions, one order higher than the approximations from the constituent r-cell stencils, each of which converge at rth order. Subsequently, Jiang and Shu [5] improved upon the WENO scheme by proposing a new smoothness indicator which allowed the WENO method to achieve the optimal (2r − 1)th order of convergence in smooth regions. This improved scheme is commonly referred to as WENO-JS.
The WENO-JS method has since become a popular choice for solving hyperbolic problems numerically due to its robustness and accuracy. However, a crucial detail remained unnoticed until the work of Henrick et al. [6], which showed that the fifth-order WENO-JS scheme fails to achieve the optimal fifth-order convergence near critical points in the solution where the first derivative of the solution vanished. They explained that the loss of accuracy was masked when a relatively large value such as 10−6 was used for ϵ , a parameter that was thought to be inconsequential and introduced solely to prevent division by zero (see Equation (10)). In addition, they proposed a mapping function gM to modify the WENO-JS stencil weights to recover the optimal order of accuracy near critical points. This scheme is referred to as WENO-M. More recently, it was shown that gM actually belongs to a more general family of functions gIM [7]. Using a different set of parameters, Feng et al. [7] showed that vastly superior results could be obtained for fifth-order WENO. Interestingly, this improved scheme, referred to as WENO-IM(2,0.1), was able to capture discontinuities with very little dissipation compared to WENO-M and WENO-JS. However, when applied to a seventh-order WENO scheme, it was reported to perform poorly [7]. There have been attempts to improve the performance of a WENO-JS scheme that involves modifying the linear weights. Martín et al. [8], for instance, proposed a bandwidth-optimised WENO scheme with a different set of linear weights for improved spectral properties. However, it is likely that such methods also suffer from the same pitfall as the WENO-JS scheme near critical points since the non-linear stencil weights are computed in an identical manner. Mapped WENO methods can likely be applied in conjunction with such methods as well to provide further improvements.
Indeed, high-order WENO schemes are considered appropriate for numerical simulations of high-speed flow scenarios, which require both the ability to capture discontinuities and resolve fine-scale turbulence [9]. Shen and Zha [10] developed a robust seventh-order WENO for a transonic convergent–divergent nozzle and RAE aerofoils. Min [11] implemented a hybrid numerical scheme based on seventh-order WENO for improving the numerical accuracy over a range of 2D and 3D unsteady aerodynamic problems. Sun et al. [4] studied the blade–vortex interaction noise using a seventh-order WENO-Z scheme and concluded that the higher-order scheme was computationally cheaperthan the lower-order schemes with comparable solution accuracy. More recently, Han et al. [12] resolved the helicopter wake flow during hover flight using seventh-order WENO-K scheme and obtained excellent agreement with the measured data at both subsonic and transonic flow regimes. Clearly, high-order WENO schemes are constantly adopted and implemented for improved solution accuracies in various flow problems relevant to aerospace engineering applications.
In this study, the performances of the mapped methods near discontinuities are analyzed for fifth- and seventh-order WENO schemes. In particular, the cause behind the poor performance of the seventh-order WENO-IM(2,0.1) scheme is studied in detail. Based on the insights gleaned, a new family of rational mapping functions is introduced, which, incidentally, includes the mapping functions gIM as a special case.
This paper is organised as follows: In Section 2, a brief review of WENO-JS and the mapped WENO schemes is provided. Next, in Section 3, the new family of rational mapping functions gRM(ω;k,m,s) is introduced and a parametric study is performed to determine suitable coefficients k, m and s. In Section 4, the proposed mapped WENO method, WENO-RM(k,m,s), is compared with existing WENO schemes for several linear advection and Euler cases to demonstrate the superior performance of the new method. Finally, the concluding remarks are provided in Section 5.

2. Numerical Methods

2.1. WENO-JS Scheme

To begin with, consider a one-dimensional domain x [ a , b ] discretized into N cells of width Δ x = ( b a ) / N . The center and interface locations of the ith cell are denoted by xi and x i ± 1 / 2 = x i ± Δ x / 2 , respectively. In the finite volume methodology, the discrete solution to a problem is determined in terms of cell averages. Given a function u(x), the ith cell average is denoted by u i ¯ . In this section, the left-biased fifth- and seventh-order WENO-JS approximations to the point value u i + 1 / 2 L = lim x x i + 1 / 2 u ( x ) are briefly outlined. The stencils for the fifth- and seventh-order reconstructions are shown in Figure 1.
The WENO-JS procedure first computes the rth order approximations to u i + 1 / 2 L as:
u i + 1 / 2 L = ( u ^ i + 1 / 2 L , r ) j + O ( Δ x r ) j [ 0 , r 1 ] .
For the fifth-order WENO-JS scheme (r = 3), three third order approximations are obtained using the expressions:
( u ^ i + 1 / 2 L , 3 ) 0 = 1 6 ( 2 u i 2 ¯ 7 u i 1 ¯ + 11 u i ¯ ) , ( u ^ i + 1 / 2 L , 3 ) 1 = 1 6 ( u i 1 ¯ + 5 u i ¯ + 2 u i + 1 ¯ ) , ( u ^ i + 1 / 2 L , 3 ) 2 = 1 6 ( 2 u i ¯ + 5 u i + 1 ¯ u i + 2 ¯ ) .
Similarly, for the seventh-order WENO-JS scheme (r = 4), four fourth-order approximations are obtained using the expressions:
( u ^ i + 1 / 2 L , 4 ) 0 = 1 12 ( 3 u i 3 ¯ + 13 u i 2 ¯ 23 u i 1 ¯ + 25 u i ¯ ) , ( u ^ i + 1 / 2 L , 4 ) 1 = 1 12 ( u i 2 ¯ 5 u i 1 ¯ + 13 u i ¯ + 3 u i + 1 ¯ ) , ( u ^ i + 1 / 2 L , 4 ) 2 = 1 12 ( u i 1 ¯ + 7 u i ¯ + 7 u i + 1 ¯ u i + 2 ¯ ) , ( u ^ i + 1 / 2 L , 4 ) 3 = 1 12 ( 3 u i 1 ¯ + 13 u i ¯ 5 u i + 1 ¯ + u i + 2 ¯ ) .
By performing a weighted average of these rth order approximations using the optimal weights d j ( r ) , it is possible to obtain the (2r – 1)th order upstream central approximation u i + 1 / 2 L , C ( 2 r 1 ) as:
u i + 1 / 2 L , C ( 2 r 1 ) = j = 0 r 1 d j ( r ) ( u ^ i + 1 / 2 L , r ) j .
The optimal weights for fifth- and seventh-order WENO-JS schemes are given in Equations (5) and (6), respectively.
d 0 ( 3 ) = 1 10 , d 1 ( 3 ) = 6 10 , d 2 ( 3 ) = 3 10
d 0 ( 4 ) = 1 35 , d 1 ( 4 ) = 12 35 , d 2 ( 4 ) = 18 35 , d 3 ( 4 ) = 4 35
For each rth order approximation, a corresponding smoothness indicator is computed. Following Ref. [5], the rth order smoothness indicator I S j ( r ) is defined as:
I S j ( r ) = l = 1 r 1 Δ x 2 l 1 x i 1 / 2 x i + 1 / 2 [ d l p j ( x ) d x l ] 2 d x ,
where p j ( x ) refers to the rth order polynomial reconstructed on the jth stencil, which extends from cell (i + jr +1) to cell (i + j). The explicit forms of I S j ( r ) are given for r = 3 and r = 4 in Equations (8) and (9), respectively.
I S 0 ( 3 ) = 1 4 ( u i 2 ¯ 4 u i 1 ¯ + 3 u i ¯ ) 2 + 13 12 ( u i 2 ¯ 2 u i 1 ¯ + u i ¯ ) 2 I S 1 ( 3 ) = 1 4 ( u i 1 ¯ u i + 1 ¯ ) 2 + 13 12 ( u i 1 ¯ 2 u i ¯ + u i + 1 ¯ ) 2 I S 2 ( 3 ) = 1 4 ( 3 u i ¯ 4 u i + 1 ¯ + u i + 2 ¯ ) 2 + 13 12 ( u i ¯ 2 u i + 1 ¯ + u i + 2 ¯ ) 2
I S 0 ( 4 ) = ( 1 3 u i 3 ¯ 3 2 u i 2 ¯ + 3 u i 1 ¯ 11 6 u i ¯ ) 2 + 13 12 ( u i 3 ¯ 4 u i 2 ¯ + 5 u i 1 ¯ 2 u i ¯ ) 2 + 781 720 ( u i 3 ¯ 3 u i 2 ¯ + 3 u i 1 ¯ u i ¯ ) 2 I S 1 ( 4 ) = ( 1 6 u i 2 ¯ u i 1 ¯ + 1 2 u i ¯ + 1 3 u i + 1 ¯ ) 2 + 13 12 ( u i 1 ¯ 2 u i ¯ + u i + 1 ¯ ) 2 + 781 720 ( u i 2 ¯ 3 u i 1 ¯ + 3 u i ¯ u i + 1 ¯ ) 2 I S 2 ( 4 ) = ( 1 3 u i 1 ¯ + 1 2 u i ¯ u i + 1 ¯ + 1 6 u i + 2 ¯ ) 2 + 13 12 ( u i 1 ¯ 2 u i ¯ + u i + 1 ¯ ) 2 + 781 720 ( u i 1 ¯ 3 u i ¯ + 3 u i + 1 ¯ u i + 2 ¯ ) 2 I S 3 ( 4 ) = ( 11 6 u i ¯ 3 u i + 1 ¯ + 3 2 u i + 2 ¯ 1 3 u i + 3 ¯ ) 2 + 13 12 ( 2 u i ¯ 5 u i + 1 ¯ + 4 u i + 2 ¯ u i + 3 ¯ ) 2 + 781 720 ( u i ¯ 3 u i + 1 ¯ + 3 u i + 2 ¯ u i + 3 ¯ ) 2
With the smoothness indicators computed, the non-linear WENO-JS stencil weights ωj can be computed as
α j = d j ( I S j + ϵ ) 2 , ω j = α j j = 0 r 1 α j .
The superscript (r) will be dropped for the sake of conciseness. The value of ϵ is taken to be 10−6 for WENO-JS scheme and 10−40 for all the mapped WENO schemes, as recommended in Ref. [6]. Finally, the (2r −1)th order WENO-JS approximation u i + 1 / 2 L , W is obtained by taking a weighted average of the rth order approximations using ωj, as shown in Equation (11).
u i + 1 / 2 L , W = j = 0 r 1 ω j ( u ^ i + 1 / 2 L ) j

2.2. Mapped WENO Schemes

In order for the WENO approximation in Equation (11) to converge at the optimal (2r −1)th order in smooth regions, ωj has to satisfy the criterion
ω j d j = O ( Δ x r 1 ) .
This holds true for all values of r. However, the WENO-JS scheme does not satisfy this criterion near critical points. A critical point is defined as a point where one or more derivatives of the function u(x) vanish. At a critical point x = xc of order ncp (i.e., a critical point of nth order), the first ncp derivatives of the function u(x) vanish, i.e., u ( x c ) = = u ( n c p ) ( x c ) = 0 . It can be shown that the WENO-JS approximation exhibits the convergence behaviour
ω j d j = O ( Δ x r 1 n c p ) ,
for n c p = 1 , , r 1 instead of that in Equation (12) [7]. Thus, when no derivatives vanish (ncp = 0), WENO-JS approximation converges at optimal order, whereas the order of convergence deteriorates as ncp increases.
In order to restore the order of convergence, and, hence, the accuracy near critical points, Henrick et al. [6] suggested the use of a mapping function g(ω) with the following properties:
(a)
g ( 0 ) = 0 , g ( d ) = d and g ( 1 ) = 1 ;
(b)
Monotone increasing with finite slopes for ω [ 0 , 1 ] ;
(c)
g ( d ) = = g ( k ) ( d ) = 0 g ( k + 1 ) ( d ) , i.e., the function flattens near ω = d .
The mapping function g(ω) can be used to calculate the mapped weights ω ˜ j from the WENO-JS weights ωj obtained from Equation (10), as shown in Equation (14). The mapped WENO approximation of u i + 1 / 2 L is obtained by replacing ωj in Equation (11) with the mapped weights
ω ˜ j = g ( ω j ) j = 0 r 1 g ( ω j ) .
Following their analyses, the mapping function g(ω) can be expanded using Taylor series about ωj = dj to demonstrate that the mapping process is indeed able to recover the optimal order of accuracy near critical points.
g ( ω j ) = g ( d j ) + l = 1 k 1 l ! g ( l ) ( d j ) ( ω j d j ) l + 1 ( k + 1 ) ! g ( k + 1 ) ( d j ) ( ω j d j ) k + 1 +
Equation (15) can be simplified using the properties of the mapping function mentioned earlier. From property (a), g(dj) = dj, and all the terms inside the summation vanish due to property (c). Finally, substituting Equation (13) into the result and rearranging yields
g ( ω j ) = d j + 1 ( k + 1 ) ! g ( k + 1 ) ( d j ) ( ω j d j ) k + 1 + = d j + O ( Δ x ( k + 1 ) ( r 1 n c p ) ) .
Substituting Equation (16) into Equation (14) and using the fact that dj sum to unity, results in
ω ˜ j d j = O ( Δ x ( k + 1 ) ( r 1 n c p ) ) .
Comparing Equation (17) with Equation (12), it can be concluded that as long as ( k + 1 ) ( r 1 n c p ) r 1 , the mapped weights ω ˜ j satisfy the required criterion. For example, for fifth-order convergence at a critical point of order ncp = 1, k must be equal to or greater than 1.
The mapping functions put forward by Henrick et al. [6] and Feng et al. [7] denoted by gM(ω) and gIM(ω;k,A), respectively, are given in Equations (18) and (19) as a reference prior to introducing the proposed mapping function in the present study.
g M ( ω ) = ω ( d + d 2 3 d ω + ω 2 ) d 2 + ( 1 2 d ) ω
g IM ( ω ; k , A ) = d + A ( ω d ) k + 1 A ( ω d ) k + ω ( 1 ω )
In the function gIM(ω;k,A), k is a positive even integer and A is a positive real number. It can be shown directly that g M ( ω ) = g IM ( ω ; 2 , 1 ) , i.e., gM belongs to the gIM family of functions. Moreover, the gIM family of functions satisfies all the properties (a)–(c) mentioned earlier.
Feng et al. [7] showed that g IM ( ω ; 2 , 0.1 ) vastly outperforms gM(ω) (i.e., g IM ( ω ; 2 , 1 ) ) for fifth-order WENO schemes. Both the mapping functions are shown in Figure 2 for the smallest and largest values of the optimal weights d j ( r ) for r = 3 and r = 4. Notice that g IM ( ω ; 2 , 0.1 ) is much flatter than gM(ω) near ωj = dj, signifying that g IM ( ω ; 2 , 0.1 ) attempts to map a wider range of weights in the neighbourhood of ωj = dj to the optimal weight dj. As a result, g IM ( ω ; 2 , 0.1 ) is quicker to bring the mapped weights ω ˜ j close to the optimal weights dj compared to gM(ω).

3. A New Family of Mapping Functions

3.1. Motivation

Consider the one-dimensional linear advection problem
u t + u x = 0 ,
which represents one of the most straightforward hyperbolic problems. Given the initial condition u ( x , t = 0 ) = u 0 ( x ) , its exact solution is simply u ( x , t ) = u 0 ( x t ) , which makes it a suitable starting point for the assessment of numerical schemes. It must be remarked that the Euler equations reduce to linear advection of density in the absence of velocity and pressure gradients.
Averaging Equation (20) over the ith cell results in the ordinary differential equation (ODE) for the cell average u ¯ i
d u ¯ i d t = 1 Δ x ( u i + 1 / 2 L u i 1 / 2 L ) .
Notice that the right-hand side of Equation (21) involves the left-biased approximations of u(x) at the cell interfaces, which can be obtained through the WENO schemes described in the previous sections. In the present study, this ODE was solved in time using the explicit third-order TVD Runge–Kutta time-marching scheme [13].
With the numerical schemes defined, consider the initial condition of two constant states separated by discontinuities at x = 0 and x = ±0.5
C a s e   A 1   u ¯ i = { + 1 , 0.5 x i < 0 1 , 0 x i 0.5 .
This problem was solved on the periodic domain x [ 0.5 , 0.5 ] discretized into N = 100 uniform cells using fifth- and seventh-order WENO schemes until t = 100 s (100 cycles) using CFL = 0.1. Note that the mapped WENO methods using gM(ω) and g IM ( ω ; 2 , 0.1 ) are denoted by WENO-M and WENO-IM(2,0.1), respectively. The final profiles are shown in Figure 3.
It can be observed from Figure 3a that at fifth order, the mapped WENO schemes WENO-M and WENO-IM(2,0.1) are able to capture the discontinuities with lesser numerical diffusion compared to the WENO-JS scheme. Recall that mapped WENO methods were developed to achieve optimal accuracy near smooth critical points. Yet, they seem to be able to capture sharper discontinuities than the original unmapped WENO-JS scheme. Between the mapped WENO schemes, WENO-IM(2,0.1) performs better as it captures the flat regions on both sides of the discontinuities more accurately. In contrast, WENO-M suffers a loss of accuracy to the right of the discontinuities and performs similar to WENO-JS scheme. On the other hand, the results become different at seventh order, as seen from Figure 3b. Both the mapped WENO schemes deliver a similar, and at some regions even worse, performance compared to WENO-JS scheme. This can be determinantal when simulating the propagation of discontinuities and sharp gradients over long distances.
WENO-M and WENO-IM(2,0.1) both suffer a particularly significant loss in accuracy at the flat regions slightly to the left of the discontinuities. To better understand this loss in accuracy, the mapped weights were analysed in detail. The mapped weights ω ˜ j of the leftmost and rightmost stencils, i.e., j = 0 and j = 2, at the end of the first cycle (t = 1 s) are plotted for fifth-order WENO-M and WENO-IM(2,0.1) schemes in Figure 4. Comparing Figure 4a,b, a common trend can be discerned—the leftmost stencil weight ω ˜ 0 is larger than the rightmost stencil weight ω ˜ 2 in the region to the left of the discontinuity and vice versa to the right. Although the regions on either side of the discontinuity appear smooth and uniform, they are in fact highly non-smooth in the sense that the “flatness” of the profile changes vastly between one stencil to the next. This results in the smoothness indicators ISj of adjacent stencils being an order of magnitude apart or more. Since a single WENO reconstruction involves multiple stencils (three for fifth order, four for seventh order), the resulting minimum and maximum values of ISj are several orders of magnitudes apart. This explains the behaviour of the weights mentioned earlier. Interestingly, across the discontinuity itself, the weights remain close to the optimal weights dj since the discontinuity has been sufficiently “smoothened” by numerical diffusion.
A close examination of the weights near discontinuities suggests that the fifth-order WENO-IM(2,0.1) scheme allocates smaller weights to the non-smooth stencils (j = 2 and j = 0 to the left and right of the discontinuity, respectively) than the fifth-order WENO-M scheme at the end of first cycle, as shown in Figure 4a,b. This is quite unexpected since WENO-IM(2,0.1) scheme generally allocates greater weights to non-smooth stencils than WENO-M scheme, as evident from the shapes of their mapping functions shown in Figure 2. The only conclusion is that WENO-M scheme produces a smoother, more diffused profile close to and across the discontinuity compared to WENO-IM(2,0.1) scheme at the end of the first cycle. This poses the question as to why the fifth-order WENO-IM(2,0.1) scheme captures a much sharper discontinuity compared to the fifth-order WENO-M scheme. It could be related to the behaviour of the mapping process across the discontinuity itself. Closer examination of the mapped weights across the discontinuity revealed that the mapped weights from the fifth-order WENO-IM(2,0.1) remained very close to their optimal values while those from the fifth-order WENO-M scheme deviated slightly. This is due to g IM ( ω ; 2 , 0.1 ) having a wider flat portion in the neighbourhood of ω j = d j than g M ( ω ) (see Figure 2). As a result, g IM ( ω ; 2 , 0.1 ) attempts to map a wider range of weights to the optimal weights. Therefore, it is postulated that the fifth-order WENO-IM(2,0.1) scheme captures a sharper discontinuity because it is quicker to revert to the upstream central scheme at the smeared discontinuity and because it amplifies small weights ( ω j < d j ) to a greater extent in general.
Moving on to seventh-order WENO schemes, the mapped WENO weights of the leftmost (j = 0) and the rightmost (j = 3) stencils WENO-M and WENO-IM(2,0.1) schemes at the end of the first cycle (t = 1 s) are similarly shown in Figure 5. The behaviour of the seventh-order mapped weights is similar to the fifth-order ones across the discontinuity and to its right. To the left of the discontinuity, however, the mapped weights ω ˜ 3 are relatively larger; for WENO-IM(2,0.1), they are even comparable to ω ˜ 0 at some locations. This is a crucial difference between the fifth- and seventh-order mapped WENO schemes. It is believed that this difference is due to the small value of d 0 ( 4 ) = 1 35 . Referring to the close-up view in Figure 2c, it can be seen that both mapping functions g M ( ω ) and g IM ( ω ; 2 , 0.1 ) result in particularly severe amplifications to the small weights compared to the identity map (WENO-JS). In fact, the slope of g IM ( ω ; k , A ) at ω = 0
g IM ( 0 ; k , A ) = 1 + 1 A d k 1 ,
is strongly proportional to the inverse of d. Feng et al. [7] have pointed this out and recommended the use of the WENO-M scheme instead of the WENO-IM(2,0.1) scheme for seventh order since g M ( 0 ) < g IM ( 0 ; 2 , 0.1 ) . However, it is clear from Figure 3b and Figure 5 that the WENO-M scheme, while somewhat better than the WENO-IM(2,0.1) scheme, is still inadequate in overcoming the amplification problem for seventh order.

3.2. Proposed Solution

Based on the discussion in the previous section, it appears that limiting the amplification of small weights is important for seventh-order mapped WENO schemes. This idea was pursued in a couple of earlier studies [14,15] by ensuring that the slope of the mapping function vanishes at ω = 0. This actively suppresses near-zero weights. In the present paper, a different, rather simpler approach is suggested to overcome the amplification problem. Consider a general family of rational mapping functions g RM given as
g RM ( ω ; k , m , s ) = d + ( ω d ) k + 1 ( ω d ) k + s [ ω ( 1 ω ) ] m ,
where k is a positive even integer, m is a positive integer and s is a positive scaling factor. The corresponding mapped WENO scheme will be referred to as WENO-RM(k, m, s). gRM includes gIM as a special case for m = 1 and s = A 1 . Since ω is a positive real value less than or equal to 1, ω(1 − ω) is a positive quantity with an upper bound of 1 4 . Therefore, increasing the value of m allows the term [ω(1 − ω)]m in the denominator to diminish and, therefore, causes the mapping function to move closer to the identity map g(ω) = ω. In fact, for m > 1, g RM ( 0 ; k , m , s ) = 1 , i.e., gRM functions are able to closely follow the identity map near ω = 0 as shown in the inset plots in Figure 6. This means that, in contrast to previous studies in which small weights were actively suppressed, small weights are preserved with little or no amplification upon mapping with gRM functions. Therefore, for m > 1, the proposed mapping function gRM is more suitable for mapped WENO schemes at seventh order and above.
Increasing the values of k and s have the opposite effect of increasing m. Increasing the value of k causes more derivatives of gRM to vanish at ω = d. Similarly, increasing the value of s would increase the denominator of the rational expression causing the constant term d to dominate. In both cases, the flat portion of the mapping function in the neighbourhood of ω = d is widened, thereby extending the range over which weights get mapped to optimal weights. Feng et al. [7] did not investigate g IM ( ω ; k , A ) mapping functions beyond k = 2 because increasing k severely amplifies the small weights (see Equation (23)). However, the introduction of parameter m overcomes this restriction and, therefore, a greater range of k values can be explored for mapping functions, which couple better with higher-order WENO methods with m > 1.
Before assessing the performance of any mapping function, it is crucial to verify that it satisfies all the three properties (a)–(c) mentioned earlier. It can be easily verified that gRM satisfies properties (a) and (c). However, it does not satisfy property (b) unconditionally. gRM can be non-monotone and determining the exact conditions for monotonicity is not very straightforward. Nevertheless, it can be shown that km − 1 is a sufficient condition for monotonicity regardless of the value of s and d (refer to Appendix A). Therefore, the present discussion will be restricted to km − 1. It is useful to mention that such restriction is considered reasonable because there are valid values of k that satisfy the monotonicity condition for every value of m.

3.3. Parametric Study

A parametric study was conducted to determine good combinations of parameters (k,m,s) using a simple one-dimensional linear advection problem [Equation (20)] with the initial condition given as
C a s e   A 2   u ¯ i = { 1 6 [ G ( x i , + δ ) + G ( x i , δ ) + 4 G ( x i , 0 ) ] 0.8 x i 0.6 1 0.4 x i 0.2 1 | 10 ( x i 0.1 ) | 0 x i 0.2 1 6 [ F ( x i , + δ ) + F ( x i , δ ) + 4 F ( x i , 0 ) ] 0.4 x i 0.6 0 otherwise ,
where G ( x , z ) = exp [ log 2 36 δ 2 ( x + 0.7 z ) 2 ] , F ( x , a ) = max [ 0 , 1 100 ( x 0.5 a ) 2 ] and δ = 0.005. From left to right, the initial profile consists of a Gaussian, a square wave, a triangle wave and half-ellipse. This initial condition was first used by Jiang and Shu [5] in their classic paper on WENO schemes. Henceforth, it has become a popular choice for benchmarking numerical methods for hyperbolic problems.
The problem was solved on the domain x [ 1 , 1 ] at two resolutions, N = 400 and 800 cells. The simulation was run until t = 200 s (100 cycles) using CFL = 0.1. The parameter combinations tested are listed in Table 1. The corresponding L1 error norms are plotted in Figure 7 and Figure 8 for N = 400 and 800 cells, respectively.
A consistent trend can be observed in all cases. For any given combination of k and m, the errors tend to be large for small values of s. Then, the errors reduce with increasing s before increasing slightly at very large values of s. The explanation for this trend is quite intuitive—increasing the value of s extends the range of weights which are mapped to their optimal values. Too small a value of s would not map the weights to optimal weights ‘rapidly’ enough resulting in greater numerical diffusion (see discussion in Section 3.1). On the other hand, too large a value would cause the weights to be mapped to optimal weights over-aggressively producing spurious oscillations. This effect is illustrated in Figure 9 for k = 4 and m = 2 at three different values of s. When s = 2, all four profiles suffer distortions. Upon increasing s to 200, the distortions are nearly completely eliminated, and all four profiles are captured well. A small spurious oscillation occurs at the foot of the half-ellipse at x ≈ 0.6. When increased further to s = 2 × 104, prominent spurious oscillations form at the square wave and at the foot of the triangle wave.
Another notable trend is that increasing the value of m almost always results in larger errors at lower values of s. Increasing m causes the mapping function to follow the identity map more closely and decreases the range over which weights are mapped to optimal weights. As such, a large value of s is required to counteract its effect and widen the range before the errors start to reduce. Finally, increasing the value of k tends to have a favourable effect, suggesting that a wide flat region about ωj = dj is crucial with the only exception taking place when k is increased from 4 to 6 with m = 2.
Among the various combinations of parameters k, m and s tested, k = 6 and m = 3 consistently results in very low errors over the entire range of s values at both resolutions. However, a suitable value for s appears to be dependent on the problem and resolution. For instance, among the values of s tested, s = 2 × 103 results in the smallest L1 error norm when using N = 400 cells. It reduces to s = 2 when using N = 800 cells. However, the choice of the value of s may not be too critical since it is evident from Figure 7 and Figure 8 that the error norms remain small over a wide range of values of s when k = 6 and m = 3.
To conclude this section, let us return to linear advection of sharp discontinuities, which motivated the development of this new family of mapping functions. Figure 10a,b show the mapped weights for this problem at t = 1 s (1 cycle) for two different values of s with k = 6 and m = 3. Notice that in both cases, the mapped weights behave in the desired manner without any severe amplification in the flat regions on either side of the smeared discontinuity. It can also be seen that the mapped weights for the non-smooth stencils (j = 0 to the left of the discontinuity and j = 3 to its right) reduce when s is increased. This indicates that the discontinuity is captured with less dissipation when using s = 2 × 103. The final profiles at t = 100 s (100 cycles) using both values of s are shown in Figure 10c. In contrast to the results shown in Figure 3b, the discontinuities are captured well without any distortions for both values of s even at relatively long output times. The close-up view shows that using s = 2 × 103 indeed results in a slightly sharper discontinuity compared to s = 2. The value of s was chosen to be 2 × 103 for all the test cases reported in the following section with the exception of the double-Mach reflection Euler problem for which a value of s = 2 was chosen to ensure a stable run.

4. Results

In this section, the proposed schemes are compared with other WENO schemes such as WENO-JS, WENO-M, WENO-IM(2,0.1), WENO-RM(260) [15] and WENO-Z [16]. The value of ϵ was set to 10−40 for all schemes except for WENO-JS and WENO-RM(260) for which it was set to 10−6 and 10−99 [15], respectively. The values of parameters k and m were fixed at 6 and 3, respectively, for all cases. Parameter s was set to 2 × 103 or 2 corresponding to the WENO-RM(6,3,2 × 103) or WENO-RM(6,3,2) schemes, respectively. The CFL number was fixed at 0.1 and 0.5 when solving the advection equation and the Euler equations, respectively.

4.1. Approximate Dispersion Relation

The dispersion relation of a numerical scheme can be used to assess the scheme’s ability to resolve Fourier modes of different wavelengths on a uniform grid. Given a Fourier mode with a wavelength l, its associated wavenumber on a grid with cell spacing Δ x is φ = 2 π Δ x / l . Numerical schemes with finite stencils invariably introduce errors in resolving the mode, which results in a modified wavenumber number Φ . However, it is not possible to derive the exact dispersion relations of non-linear schemes. Instead, Pirozzoli [17] proposed an empirical method to compute approximate dispersion relations of non-linear shock-capturing schemes. The approximate dispersion relations can be used to assess the resolution of the reconstruction scheme in wavenumber space. A uniform N-point mesh supports wavenumbers φ k = 2 π k / N where k [ 0 , N / 2 ] . For each φ k , the one-dimensional linear advection equation Equation (20) is initialized with the sinusoidal mode
u j ( 0 ) = cos ( φ k j Δ x ) ,
for j [ 0 , N 1 ] and advanced up to a small-time step Δt using the first-order Euler scheme. Next, Fourier transforms are applied to the initial and advanced profiles, and the complex amplitudes u ^ k ( 0 ) and u ^ k ( Δ t ) associated with the kth wavenumber φ k are used to compute the kth modified wavenumber Φ k as:
Φ k = Δ x i Δ t log ( u ^ k ( Δ t ) u ^ k ( 0 ) ) .
The real and imaginary parts of Φ k pertain to the dispersion and diffusion properties of the numerical scheme, respectively. The approximate dispersion relations for the various seventh-order schemes are plotted in Figure 11. The seventh-order upstream central scheme is the linear scheme whereby one simply takes the sub-stencil weights to be the optimal weights, i.e., ωj = dj. It represents the optimal resolution that one could achieve from a WENO scheme that relies on the same optimal weights. Hence, it is more meaningful to compare the WENO schemes with the upstream central scheme rather than with spectral methods. The results show that the unmapped WENO-JS scheme has the most deviation from the upstream central scheme and, therefore, the poorest resolution among all schemes. While all mapped WENO schemes and WENO-Z scheme bring about significant improvement over the WENO-JS scheme, the WENO-RM(6,3,2 × 103) scheme comes the closest to the upstream central scheme in terms of both dispersion and diffusion properties, making it ideal for resolving small-scale flow structures, their convective effects and the induced wave propagation in a high-fidelity simulation of compressible flows where shock waves may also be present. Though the WENO7-IM(2,0.1) scheme follows close behind the WENO7-RM(6,3,2 × 103) scheme, one must bear in mind that it performs poorly for cases with discontinuities, as demonstrated earlier in Section 3.1.

4.2. Linear Advection

The one-dimensional linear advection problem [Equation (20)] was solved on the periodic domain x [ 0.5 , 0.5 ] with the smooth initial condition
C a s e   A 3   u ¯ i = sin 9 ( 2 π x i ) ,
at three different resolutions N = 100, 200 and 400 cells. The L1 error norms at t = 100 s (100 cycles) are listed in Table 2.
It can be seen that WENO-JS scheme performs better than both WENO-M and WENO-IM(2,0.1) schemes at all three resolutions. At N = 400 cells, it even surpasses WENO-Z. This is mainly because the relatively large value of ϵ = 10 6 overpowers the much smaller values of the smoothness indicators ISj in the nearly flat regions. As a result, the WENO-JS scheme behaves essentially like the upstream central scheme (ωjdj). Both WENO-M and WENO-IM(2,0.1) schemes perform poorly for this problem, which shows that they are not suitable at seventh order and higher. The proposed scheme WENO-RM(6,3,2 × 103) delivers a better performance than all other schemes except WENO-RM(260). Nevertheless, its performance improves with increasing resolution, such that at the finest resolution of N = 400 cells, it performs nearly identical to WENO-RM(260).
For the next advection problem, the initial condition
C a s e   A 4   u ¯ i = { [ 1 sin ( 4 π x i ) ] / 2 1 / 8 < x i 1 / 2 1 / 2 1 / 2 < x i 7 / 8 0 otherwise ,  
was obtained from Ref. [18]. It consists of discontinuities in u at x = 7 / 8 , u at x = 1 / 2 and u at x = 1 / 8 . It also possesses a smooth maximum at x = 3 / 8 . Once again, the problem was solved on the periodic domain x [ 0 , 1 ] at three different resolutions N = 100, 200 and 400 cells. The L1 error norms at t = 100 s (100 cycles) are listed in Table 3. The final profiles for N = 100 cells are plotted in Figure 12.
For this case, WENO-RM(6,3,2 × 103) scheme consistently delivers the best performance at all three resolutions. Once again, both mapped WENO schemes WENO-M and WENO-IM(2,0.1) deliver poor performances especially near the discontinuity in u at x = 7/8 and the discontinuity in u at x = 1/8, as shown in Figure 12. The WENO-JS scheme is unable to capture the smooth maximum at the coarse resolution of N = 100 cells, but it is able to do so at the finer resolutions. WENO-Z and WENO-RM(260) schemes deliver similar performances except at the coarsest resolution of N = 100 cell, whereby the WENO-RM(260) scheme results in a small overshoot on the left side of the discontinuity at x = 7/8. This is quite surprising since WENO-RM(260) is a mapped WENO scheme that actively suppresses small weights. The overshoot does not occur at the finer resolutions though.

4.3. Euler Equations

The system of one-dimensional Euler equations for an ideal gas is given as
t ( ρ ρ u ρ e ) + x ( ρ u ρ u 2 + p ρ h u ) = 0 ,
where ρ, u and p refer to the gas density, velocity and pressure, respectively. The specific total energy is given by e = p / ( γ 1 ) ρ + u 2 / 2 and the specific total enthalpy by h = e + p / ρ . γ is the ratio of specific heats with a value of 1.4. The system was solved using a hybrid method coupling a MUSCL-type flux approach and HLL-based characteristic flux approach using a shock sensor [19].
The first case Is a shock tube problem proposed by Lax [20]. The initial condition
C a s e   E 1   ( ρ , u , p ) = { ( 0.445 , 0.698 , 3.528 ) x < 0 ( 0.500 , 0.000 , 0.571 ) x 0 ,
consists of two constant states separated by a discontinuity at x = 0. The problem was solved on the domain x [ 0.5 , 0.5 ] until t = 0.13 s. The final density profiles are shown in Figure 13.
For this case, all the improved WENO schemes capture a sharper contact discontinuity at x = 0.2 and a sharper shock at x = 0.32 compared to the WENO-JS scheme. WENO-RM(6,3,2×103) scheme captures the sharpest discontinuities among all schemes. All the schemes, including the WENO-JS scheme, result in spurious oscillations to the right of the contact discontinuity to some extent, but WENO-Z is the only scheme to cause a large overshoot and a small undershoot to left and right of the shock, respectively. It must be mentioned that the WENO-RM(260) scheme results in the smallest spurious oscillation among the improved WENO schemes.
The next two test cases are shock-entropy wave interaction problems proposed by Shu and Osher [21] and Titarev and Toro [22]. The corresponding initial conditions are given in Equations (32) and (33), respectively.
C a s e   E 2   ( ρ , u , p ) = { ( 3.857143 , 2.629369 , 10.33333 ) x < 4 ( 1 + 0.2 sin ( 5 x ) , 0 , 1 ) x 4
C a s e   E 3   ( ρ , u , p ) 0 = { ( 1.515695 , 0.523346 , 1.80500 ) x < 4.5 ( 1 + 0.1 sin ( 20 x ) , 0 , 1 ) x 4.5  
The interaction of the shock with the entropy wave results in short-wavelength structures whose amplitudes are a good measure of the numerical viscosity of the numerical scheme. These structures are representative of the large-scale eddies encountered, for instance, in the wake of an airfoil undergoing flow separation. Case E3 is a more severe version of case E2 as it results in structures of much shorter-wavelength structures compared to the latter case, potentially more representative of the finer-scale turbulent structures and their fluctuations in the flow. Case E2 was solved on the domain x [ 5 , 5 ] discretized into N = 200 uniform cells until t = 1.8 s. Case E3 was solved on the same domain discretized into N = 1000 cells until t = 5 s. The reference solutions for both cases were obtained on a fine grid with ten times the resolution. The final density profiles are shown in Figure 14 and Figure 15.
The performance trends are similar for both cases. Once again, WENO-JS results in the smallest amplitudes for the short-wavelength structures among all schemes. The best performance is delivered by the WENO-RM(6,3,2×103) scheme, which is closely followed by WENO-IM(2,0.1) and WENO-RM(260). Surprisingly, WENO-Z does not perform as well as in the advection cases; its performance is only slightly better than that of WENO-M.
For the last case, we consider the two-dimensional double Mach reflection (DMR) problem proposed by Woodward and Colella [23], in which a strong Mach 10 shock wave impinges obliquely onto a ramp and undergoes a complex reflection producing intricate flow features. The WENO schemes were extended to higher dimensions by applying the one-dimensional schemes along the local face-normal direction. This approach is significantly cheaper and more robust compared to a fully multi-dimensional formulation. The DMR problem was computed using one of the alternative setups suggested by U S Vevek et al. [24], whereby the shock wave makes an angle of θ = tan 1 ( 7 4 ) 60.26 ° with the ramp. This enables exact tracking of the shock wave along the top boundary on a uniform structured mesh. The problem was computed on the domain x × y [ 0 , 4 ] × [ 0 , 1 ] discretized into the uniform cells of width Δ x = Δ y = 1 240 . The ramp extends along the bottom boundary from x = 1 to the right. Reflective boundary conditions are applied along the ramp and Neumann boundary conditions are applied along the remainder of the bottom boundary from x = 0. The exact pre- and post-shock conditions
C a s e   E 4   ( ρ , u , v , p ) p r e = ( 1.4 , 0 , 0 , 1 ) , ( ρ , u , v , p ) p o s t = ( 8 , 8.25 cos θ , 8.25 θ , 116.5 ) ,
were prescribed at the right and left boundaries, respectively. The terms u and v denote the horizontal and vertical velocity components, respectively. The shock wave was initialized to intersect the bottom boundary at x = 0.1 and the simulation was computed until t = 0.27 s with a re-initialization at an intermediate time t = 0.08 s to prevent undesirable numerical artefacts. The final results are plotted in Figure 16, using 41 equally spaced density contours in the range ρ [ 1.4 , 21.4 ] .
The WENO-IM(2,0.1), WENO-Z and WENO-RM(6,3,2 × 103) schemes became unstable for this problem within the first few time steps. These schemes are most likely too aggressive near the shock wave in that they do not introduce sufficient numerical diffusion. It must be recalled from the previous section that both s = 2 and s = 2 × 103 led to distortion-free results for the linear advection of the discontinuities problem (see Figure 10), and the latter was chosen on the grounds that it produced a slightly sharper discontinuity. However, since this leads to an unstable scheme for the DMR problem, the value of parameter s was chosen to be 2 for this problem alone for a stable computation. It can be seen from Figure 16 that there are not many notable differences in the results computed using the different schemes except some minor differences along the slip line. The lack of roll-ups along the slip line may be surprising to some, but it has been demonstrated in [19] that the choice of reconstructed variables and flux scheme has an effect on the extent of the details captured along the slip lines and that the usual notion of more roll-ups being indicative of a better, less diffusive scheme need not necessarily hold true in all cases.
The results obtained using the WENO-RM(6,3,2) scheme for the other Euler cases were similar to those obtained using the WENO-RM(260) scheme, but the former scheme was computationally more efficient. For instance, for the DMR problem, the WENO-RM(6,3,2) scheme was 16% slower than the unmapped WENO-JS scheme, whereas the WENO-RM(260) scheme was 83% slower.

5. Conclusions

This paper was motivated by the poor performance of the mapped WENO schemes WENO-M and WENO-IM(2,0.1) in capturing discontinuities accurately over long output times. While WENO-M performed poorly at both fifth and seventh orders, WENO-IM(2,0.1) performed poorly only at seventh order. By observing the mapped stencil weights ω ˜ j near the discontinuities, it was argued that the superior performance of the fifth-order WENO-IM(2,0.1) scheme was attributed to the wider flat portion of g IM ( ω ; 2 , 0.1 ) in the neighbourhood ωj = dj. The fifth-order WENO-M scheme, on the other hand, had a narrower flat portion at about ωj = dj, which resulted in increased numerical diffusion since it did not bring the weights close to their optimal values as quickly.
At seventh order, both WENO-M and WENO-IM(2,0.1) schemes performed poorly due to severe amplification of small weights. While such amplification also occurred at fifth order, it was not as severe due to the relatively larger numerical values of dj at fifth order compared to the seventh order. A simple solution was suggested in the form of a more general family of mapping functions g RM ( ω ; k , m , s ) which includes g IM ( ω ; k , A ) as a special case. For m > 1, these mapping functions have a unit slope at ω = 0 regardless of the values of dj. This prevents the over-amplification of small weights, making it particularly suitable for seventh order and above. The values of parameters k, m and s were determined to be 6, 3 and 2 × 103, respectively, based on a parametric study. The resulting mapped WENO scheme is referred to as WENO-RM(6,3,2 × 103).
It was shown through a number of one-dimensional linear advection problems that the WENO-RM(6,3,2 × 103) scheme outperformed both the WENO-M and WENO-IM(2,0.1) schemes. It was able to deliver similar or even superior performances compared to other popular WENO schemes such as WENO-RM(260) and WENO-Z. It was also shown that the proposed scheme performed the best among all the schemes tested, for a number of one-dimensional inviscid gas flow problems. It was able to capture the sharpest discontinuities in the shock tube problem and greatest amplitudes for the short wavelength structured in shock-entropy wave interaction problems. It should also be quickly remarked that these conclusions are equally applicable to the finite difference method as the finite difference and finite volume methods are equivalent in one dimension. The WENO-RM(6,3,2) scheme was used for the two-dimensional DMR problem alone because the WENO-RM(6,3,2 × 103) scheme was unstable for this problem. The key lesson from the DMR problem is that the value of parameter s is case dependent. In fact, even for a single case, it would be better to allow parameter s to vary locally based on the solution.
The insights gleaned from this study were the basis for the successful development of the adaptive mapped WENO scheme [25]. The adaptive scheme relies on the gRM family of mapping functions introduced in the present paper but, instead of using a constant value for s, it uses an adaptive value calculated based on the local smoothness of the solution.

Author Contributions

Conceptualization, T.H.N.; Formal analysis, U.S.V. and B.Z.; Funding acquisition, T.H.N.; Investigation, U.S.V. and B.Z.; Methodology, U.S.V.; Project administration, T.H.N.; Supervision, T.H.N.; Validation, B.Z.; Writing—original draft, U.S.V.; Writing—review & editing, B.Z. and T.H.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Singapore Ministry of Education AcRF Tier-2 Grant (grant number MOE2014-T2-1-002) and MAE Seed Project (M4082026.050) from the School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The mapping function g RM ( ω ; k , m , s ) [Equation (24)] would be monotone increasing if its first derivative g RM 0 . g RM can be written as:
g RM ( ω ; k , m , s ) = ( ω d ) k h ( ω ) { ( ω d ) k + s [ ω ( 1 ω ) ] m } 2 ,
where h(ω) is given by
h ( ω ) = ( ω d ) k I + s [ ω ( 1 ω ) ] m II [ ( 2 m k 1 ) ω 2 + ( k + 1 m 2 m d ) ω + m d ] III .
Since ω [ 0 , 1 ] and k is an even integer, all the terms in Equation (A1) are positive except for h(ω), which may become negative. Deriving the conditions for the positivity of h(ω) is non-trivial. Therefore, it is broken into three parts, as shown in Equation (A2). Terms I and II are always positive, but term III can be non-positive. Let term III be denoted by f(ω). Deriving the conditions for the positivity of f(ω) is simpler and will be pursued next. If f(ω) is positive, it directly follows that g RM ( ω ; k , m , s ) is monotone.
f(ω) is a quadratic function in ω passing through f ( 0 ) = m d > 0 and f ( 1 ) = m ( 1 d ) > 0 . There are three cases to consider:
  • Case 1: k > 2 m 1
In this case, the coefficient of ω2 in f(ω) is negative and f(ω) is a parabola with a maximum point. A parabola with a maximum point has the property that, given any interval [a,b], the function within that interval always lies above the minimum of the endpoints of the interval, i.e., f ( ω ) min [ f ( a ) , f ( b ) ] . Hence, f ( ω ) min [ m d , m ( 1 d ) ] 0 in the interval ω [ 0 , 1 ] .
  • Case 2: k = 2 m 1
In this case, the coefficient of ω2 vanishes and f(ω) becomes a linear function, which is always positive. However, this case does not occur since k is an even integer and so, k 2 m 1 .
  • Case 3: k < 2 m 1
In this case, the coefficient of ω2 in f(ω) is positive. Hence, f(ω) is a parabola with a minimum point, which may dip below zero. The coordinates of the minimum point, found by solving f ( ω ) = 0 , are:
ω m i n = m ( 1 + 2 d ) k 1 2 ( 2 m k 1 ) , f m i n = f ( ω min ) = m d [ m ( 1 + 2 d ) k 1 ] 2 4 ( 2 m k 1 ) .
The expression for fmin given in Equation (A3) is a quadratic polynomial in d with a negative coefficient for d2. Therefore, it is a parabola with a maximum point and, thus, f m i n ( d ) min [ f m i n ( a ) , f m i n ( b ) ] .
If ωmin lies outside the interval [0,1], f(ω) is monotone within the interval. Since the endpoints of f(ω) lie above zero, the function also lies above zero in the interval [0,1]. However, if ωmin lies within the interval [0,1], one must ensure that fmin ≥ 0. Substituting the expression for ωmin from Equation (A3) into the inequality 0 ≤ ωmin ≤ 1 results in the following inequality for d.
1 2 + k + 1 2 m d 3 2 k + 1 2 m
Substituting the maximum and minimum values of d from Equation (A4) into the expression for fmin in Equation (A3) yields the following results.
d = 1 2 + k + 1 2 m f m i n = k + 1 m 2 d = 3 2 k + 1 2 m f m i n = k + 1 m 2
Interestingly, both values of d lead to the same result for fmin. So, if the minimum point were to occur within the interval ω [ 0 , 1 ] , for fmin ≥ 0, km − 1. Therefore, for km − 1, f(ω) ≥ 0.

References

  1. Zhang, S.; Zhu, J.; Shu, C.-W. A brief review on the convergence to steady state solutions of Euler equations with high-order WENO schemes. Adv. Aerodyn. 2019, 1, 16. [Google Scholar] [CrossRef] [Green Version]
  2. Sun, Y.; Shi, Y.; Xu, G. Application of high-order WENO scheme in the CFD/FW-H method to predict helicopter rotor blade-vortex interaction tonal noise. Aerospace 2022, 9, 196. [Google Scholar] [CrossRef]
  3. Harten, A.; Engquist, B.; Osher, S.; Chakravarthy, S.R. Uniformly high order accurate essentially non-oscillatory schemes, III. J. Comput. Phys. 1987, 71, 231–303. [Google Scholar] [CrossRef]
  4. Liu, X.-D.; Osher, S.; Chan, T. Weighted essentially non-oscillatory schemes. J. Comput. Phys. 1994, 115, 200–212. [Google Scholar] [CrossRef] [Green Version]
  5. Jiang, G.S.; Shu, C.W. Efficient implementation of weighted ENO schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef] [Green Version]
  6. Henrick, A.K.; Aslam, T.D.; Powers, J.M. Mapped weighted essentially non-oscillatory schemes: Achieving optimal order near critical points. J. Comput. Phys. 2005, 207, 542–567. [Google Scholar] [CrossRef]
  7. Feng, H.; Huang, C.; Wang, R. An improved mapped weighted essentially non-oscillatory scheme. Appl. Math. Comput. 2014, 232, 453–468. [Google Scholar] [CrossRef]
  8. Martín, M.P.; Taylor, E.M.; Wu, M.; Weirs, V.G. A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence. J. Comput. Phys. 2006, 220, 270–289. [Google Scholar] [CrossRef]
  9. Klioutchnikov, I.; Ballmann, J.; Oliver, H.; Hermes, V.; Alshabu, A. Unsteady transonic airfoil flow simulations using high-order WENO schemes. In Hyperbolic Problems: Theory, Numerics, Applications; Springer: Berlin/Heidelberg, Germany, 2008; pp. 617–624. [Google Scholar]
  10. Shen, Y.; Zha, G. A robust seventh-order WENO scheme and its applications. In Proceedings of the 46th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 7–10 January 2008; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2008. [Google Scholar] [CrossRef] [Green Version]
  11. Min, B.-Y. A physics based investigation of Gurney flaps for enhancement of rotorcraft flight characteristics. Ph.D. Thesis, Georgia Institute of Technology, Savannah, GA, USA, 2010. [Google Scholar]
  12. Han, S.; Song, W.; Han, Z. A novel high-order scheme for numerical simulation of wake flow over helicopter rotors in hover. Chin. J. Aeronaut. 2022, 35, 260–274. [Google Scholar] [CrossRef]
  13. Gottlieb, S.; Shu, C.W. Total variation diminishing Runge-Kutta schemes. Math. Comput. 1998, 67, 73–85. [Google Scholar] [CrossRef] [Green Version]
  14. Feng, H.; Hu, F.; Wang, R. A new mapped weighted essentially non-oscillatory scheme. J. Sci. Comput. 2012, 51, 449–473. [Google Scholar] [CrossRef]
  15. Wang, R.; Feng, H.; Huang, C. A new mapped weighted essentially non-oscillatory method using rational mapping function. J. Sci. Comput. 2016, 67, 540–580. [Google Scholar] [CrossRef]
  16. Borges, R.; Carmona, M.; Costa, B.; Don, W.S. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. J. Comput. Phys. 2008, 227, 3191–3211. [Google Scholar] [CrossRef]
  17. Pirozzoli, S. On the spectral properties of shock-capturing schemes. J. Comput. Phys. 2006, 219, 489–497. [Google Scholar] [CrossRef]
  18. Abgrall, R.; De Santis, D. High order residual distribution scheme for Navier-Stokes equations. In Proceedings of the 20th AIAA Computational Fluid Dynamics Conference, Honolulu, HI, USA, 27–30 June 2011; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2011. [Google Scholar] [CrossRef] [Green Version]
  19. Vevek, U.S.; Zang, B.; New, T.H. An efficient hybrid method for solving Euler equations. J. Sci. Comput. 2019, 81, 732–762. [Google Scholar] [CrossRef]
  20. Lax, P.D. Weak solutions of nonlinear hyperbolic equations and their numerical computation. Commun. Pure Appl. Math. 1954, 7, 159–193. [Google Scholar] [CrossRef]
  21. Shu, C.-W.; Osher, S. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II. In Upwind and High-Resolution Schemes; Springer: Berlin/Heidelberg, Germany, 1989; pp. 328–374. [Google Scholar]
  22. Titarev, V.A.; Toro, E.F. Finite-volume WENO schemes for three-dimensional conservation laws. J. Comput. Phys. 2004, 201, 238–260. [Google Scholar] [CrossRef]
  23. Woodward, P.; Colella, P. The numerical simulation of two-dimensional fluid flow with strong shocks. J. Comput. Phys. 1984, 54, 115–173. [Google Scholar] [CrossRef]
  24. Vevek, U.S.; Zang, B.; New, T.H. On alternative setups of the double Mach reflection problem. J. Sci. Comput. 2019, 78, 1291–1303. [Google Scholar] [CrossRef]
  25. Vevek, U.S.; Zang, B.; New, T.H. Adaptive mapping for high order WENO methods. J. Comput. Phys. 2019, 381, 162–188. [Google Scholar] [CrossRef]
Figure 1. Schematic of stencils for fifth (r = 3) and seventh (r = 4) order WENO reconstructions.
Figure 1. Schematic of stencils for fifth (r = 3) and seventh (r = 4) order WENO reconstructions.
Aerospace 09 00623 g001
Figure 2. Mapping functions gM(ω) and g IM ( ω ; 2 , 0.1 ) for the smallest (left) and largest (right) values of d j ( 3 ) (top) and d j ( 4 ) (right).
Figure 2. Mapping functions gM(ω) and g IM ( ω ; 2 , 0.1 ) for the smallest (left) and largest (right) values of d j ( 3 ) (top) and d j ( 4 ) (right).
Aerospace 09 00623 g002
Figure 3. Performance of (a) fifth- and (b) seventh-order WENO schemes for linear advection case A1 (N = 100 cells, t = 100 s).
Figure 3. Performance of (a) fifth- and (b) seventh-order WENO schemes for linear advection case A1 (N = 100 cells, t = 100 s).
Aerospace 09 00623 g003
Figure 4. Mapped weights ω ˜ j for leftmost and rightmost stencils for fifth-order mapped WENO schemes at the end of first cycle (t = 1 s).
Figure 4. Mapped weights ω ˜ j for leftmost and rightmost stencils for fifth-order mapped WENO schemes at the end of first cycle (t = 1 s).
Aerospace 09 00623 g004
Figure 5. Mapped weights ω ˜ j for leftmost and rightmost stencils for seventh-order mapped WENO schemes at the end of first cycle (t = 1 s).
Figure 5. Mapped weights ω ˜ j for leftmost and rightmost stencils for seventh-order mapped WENO schemes at the end of first cycle (t = 1 s).
Aerospace 09 00623 g005
Figure 6. Mapping functions g RM ( ω ; 4 , m , 10 ) for different values of m for the smallest and largest values of d j ( 4 ) . The identity map g ( ω ) = ω is shown in grey.
Figure 6. Mapping functions g RM ( ω ; 4 , m , 10 ) for different values of m for the smallest and largest values of d j ( 4 ) . The identity map g ( ω ) = ω is shown in grey.
Aerospace 09 00623 g006
Figure 7. Norms for linear advection case A2 using N = 400 cells at t = 200 s.
Figure 7. Norms for linear advection case A2 using N = 400 cells at t = 200 s.
Aerospace 09 00623 g007
Figure 8. Norms for linear advection case A2 using N = 800 cells at t = 200 s.
Figure 8. Norms for linear advection case A2 using N = 800 cells at t = 200 s.
Aerospace 09 00623 g008
Figure 9. Effect of increasing s for linear advection case A2 using N = 400 cells at t = 200 s.
Figure 9. Effect of increasing s for linear advection case A2 using N = 400 cells at t = 200 s.
Aerospace 09 00623 g009
Figure 10. Mapped weights ω ˜ j for leftmost and rightmost stencils for seventh-order WENO-RM(6,3,s) at t = 1 s for (a) s = 2 and (b) s = 2 × 103. (c) Performance of both schemes for linear advection case A1 using N = 100 cells at t = 100 s.
Figure 10. Mapped weights ω ˜ j for leftmost and rightmost stencils for seventh-order WENO-RM(6,3,s) at t = 1 s for (a) s = 2 and (b) s = 2 × 103. (c) Performance of both schemes for linear advection case A1 using N = 100 cells at t = 100 s.
Aerospace 09 00623 g010aAerospace 09 00623 g010b
Figure 11. Approximate dispersion relations for various seventh-order schemes.
Figure 11. Approximate dispersion relations for various seventh-order schemes.
Aerospace 09 00623 g011
Figure 12. Performance of seventh-order WENO schemes for linear advection case A4 using N = 100 cells at t = 100 s (100 cycles).
Figure 12. Performance of seventh-order WENO schemes for linear advection case A4 using N = 100 cells at t = 100 s (100 cycles).
Aerospace 09 00623 g012
Figure 13. Performance of seventh-order WENO schemes for Lax shock tube problem (case E1) using N = 100 cells at t = 0.13 s.
Figure 13. Performance of seventh-order WENO schemes for Lax shock tube problem (case E1) using N = 100 cells at t = 0.13 s.
Aerospace 09 00623 g013
Figure 14. Performance of seventh-order WENO schemes for Shu–Osher problem (case E2) using N = 200 cells at t = 1.8 s.
Figure 14. Performance of seventh-order WENO schemes for Shu–Osher problem (case E2) using N = 200 cells at t = 1.8 s.
Aerospace 09 00623 g014
Figure 15. Performance of seventh-order WENO schemes for Titarev–Toro problem (case E3) using N = 1000 cells at t = 5 s.
Figure 15. Performance of seventh-order WENO schemes for Titarev–Toro problem (case E3) using N = 1000 cells at t = 5 s.
Aerospace 09 00623 g015
Figure 16. Performance of seventh-order WENO schemes for the DMR problem (case E4).
Figure 16. Performance of seventh-order WENO schemes for the DMR problem (case E4).
Aerospace 09 00623 g016
Table 1. Parameter combinations used in parametric study.
Table 1. Parameter combinations used in parametric study.
Kms
22, 31, 2, 5, 10, 20, 50, 100, 200, 500, 1 × 103, 2 × 103, 5 × 103, 1 × 104, 2 × 104, 5 × 104, 1 × 105
42, 3, 4, 5
62, 3, 4, 5, 6
Table 2. Error norms for linear advection case A3 at t = 100 s (100 cycles).
Table 2. Error norms for linear advection case A3 at t = 100 s (100 cycles).
SchemesN
100200400
WENO-JS5.1352 × 10−37.1533 × 10−56.4335 × 10−6
WENO-M7.7547 × 10−32.6823 × 10−47.9258 × 10−6
WENO-IM(2,0.1)9.7562 × 10−38.9693 × 10−47.8301 × 10−4
WENO-Z1.6348 × 10−36.0435 × 10−56.4376 × 10−6
WENO-RM(260)1.4095 × 10−35.7145 × 10−56.4227 × 10−6
WENO-RM(6,3,2 × 103)1.5083 × 10−35.7983 × 10−56.4225 × 10−6
Table 3. Error norms for linear advection case A4 at t = 100 s (100 cycles).
Table 3. Error norms for linear advection case A4 at t = 100 s (100 cycles).
SchemesN
100200400
WENO-JS1.6364 × 10−27.7120 × 10−35.4452 × 10−3
WENO-M1.5900 × 10−21.0101 × 10−25.2535 × 10−3
WENO-IM(2,0.1)1.5687 × 10−29.4957 × 10−34.1429 × 10−3
WENO-Z1.0286 × 10−25.1822 × 10−32.7195 × 10−3
WENO-RM(260)1.0898 × 10−25.1936 × 10−32.7181 × 10−3
WENO-RM(6,3,2 × 103)1.0106 × 10−25.0941 × 10−32.6748 × 10−3
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Vevek, U.S.; Zang, B.; New, T.H. A New Mapped WENO Method for Hyperbolic Problems. Aerospace 2022, 9, 623. https://doi.org/10.3390/aerospace9100623

AMA Style

Vevek US, Zang B, New TH. A New Mapped WENO Method for Hyperbolic Problems. Aerospace. 2022; 9(10):623. https://doi.org/10.3390/aerospace9100623

Chicago/Turabian Style

Vevek, U S, Bin Zang, and T. H. New. 2022. "A New Mapped WENO Method for Hyperbolic Problems" Aerospace 9, no. 10: 623. https://doi.org/10.3390/aerospace9100623

APA Style

Vevek, U. S., Zang, B., & New, T. H. (2022). A New Mapped WENO Method for Hyperbolic Problems. Aerospace, 9(10), 623. https://doi.org/10.3390/aerospace9100623

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop