A highly efficient low-dissipation hybrid weighted essentially non-oscillatory scheme

A highly efficient low-dissipation hybrid weighted essentially non-oscillatory scheme

X. Y. Hu B. Wang N. A. Adams Lehrstuhl für Aerodynamik, Technische Universität München
85748 Garching, Germany
School of Aerospace, Tsinghua University
100084 Beijing, China
Abstract

In this paper, we propose a simple hybrid WENO scheme to increase computational efficiency and decrease numerical dissipation. Based on the characteristic-decomposition approach, the scheme switches between the numerical fluxes with characteristic-wise WENO reconstruction and with the component-wise corresponding optimal linear reconstruction according to a new discontinuity detector. This non-dimensional discontinuity detector measures the resolution limit of the linear scheme and does not degenerate at critical points. A number of numerical examples on inviscid and viscid flow problems computed with 5th-order WENO scheme suggest that, while achieving small numerical dissipation and good robustness, the hybrid scheme generally has the total amount calls for the WENO scheme less than 2%.

keywords:
compressible flow, high-resolution method, high-order scheme

, and

1 Introduction

High-order weighted essentially non-oscillatory (WENO) schemes [4] are non-linear schemes which are generally suitable for the simulation of shock-turbulence interaction due to their high-resolution properties. However, there are still two important issues hindering their extensive application. One is that the cost of computing nonlinear weights and implementing the local characteristic decomposition is very high. The other is that WENO schemes are more dissipative than many low-dissipative linear schemes, particularly in regions without strong shock wave but with large density variation or shear rates.

One promising approach to overcome these drawbacks is using hybrid methods which switch or blend between the numerically dissipative but stable non-linear WENO scheme and a less dissipative but more accurate linear scheme according to a shock sensor or discontinuity detector [1, 5, 12, 10]. Generally, there are two types of hybrid methods based on the trouble-cell hybridization [1] and the characteristic-wise hybridization [13]. In the trouble-cell hybridization, the trouble or discontinuity-crossing cells are first detected and tagged by measuring the flow properties such as density, velocity, pressure and entropy. Then the numerical fluxes at the cell-faces of trouble cells and their neighboring cells are computed by the WENO scheme. One issue of this hybridization is the lack of smooth transition between the linear and WENO schemes. This may be due to the suddenly introduced excessive numerical dissipation in the neighbors of trouble cells, or the very different numerical properties between the WENO scheme and the chosen linear scheme [9]. Another issue is that the discontinuity detector based on flow properties may not be able to detect all discontinuities or multiple discontinuities closely located and, hence, lead to numerical instabilities [6]. Furthermore, the performance, i.e. the overall numerical dissipation, numerical stability and computational efficiency, of the trouble-cell hybridization is strongly dependent on problem, grid resolution and discontinuity-detector parameter [12, 10]. Compared to the trouble-cell hybridization, since the characteristic-wise hybridization switches between linear and WENO schemes at each cell face without taking account its neighboring cells and on each characteristic variable, it provides a much more smoother transition and is able to detect all possible discontinuities. Also, the discontinuity detector used usually measures the less problem-dependent non-dimensional quantities corresponding to the ratio between the combinations of the approximated low-order derivatives [13, 7], or the ratio between the approximated high- and low-order derivatives [17]. One issue of the characteristic-wise hybridization is that, when the solution is near critical points with zero low-order spatial derivatives, even the solution is still smooth, the discontinuity detector degenerates and switches on the WENO scheme. To overcome such difficulty, very often multiple dimensional parameters are still tuned to obtain more accurate or more stable solution for different problems. Another issue is that, up to now, the computational efficiency of the characteristic-wise hybridization usually is much less than that of the trouble-cell hybridization due to the large computational effort for the characteristic decomposition [13, 7, 17].

In this paper, we propose a simple, highly efficient hybrid WENO scheme with low dissipation based on the characteristic-wise hybridization. The scheme switches the numerical fluxes of each characteristic variables between those of the WENO scheme and its corresponding optimal linear scheme without computing the WENO weights. Different from previous characteristic-wise hybridization [13, 7, 17], the proposed linear scheme is specially designed for much higher efficiency by applying component-wise reconstruction and omitting most of the characteristic-decomposition operations. The scheme employs a new non-dimensional discontinuity detector, which measures the resolvability of the linear scheme. Compared with previous discontinuity detectors [13, 7, 17], since the new detector does not degenerate at critical points and it decreases the numerical dissipation and improves considerably the accuracy of the solution in the smooth region. Furthermore, we show by extensive numerical examples that, by choosing a broadly effective switch threshold, the hybrid scheme has achieved the total amount calls for the WENO scheme generally less than 2% without compensating numerical robustness.

2 Method

We assume that the fluid is inviscid and compressible, described by the one-dimensional Euler equations as

(1)

where , and . This set of equations describes the conservation laws for mass density , momentum density and total energy density , where is the specific internal energy. To close this set of equations, the ideal-gas equation of state with a constant is used.

2.1 Characteristic-wise WENO scheme

For completeness, we briefly recall the classical 5th-order WENO scheme [4] for solving Eq. (1). The discretization is within the spatial domain such that , , where is the spatial step, the semi-discretized form by the method of lines yields a system of ordinary differential equations

(2)

where are the numerical flux at , respectively. Once the right-hand side of this expression has been evaluated, a TVD Runge-Kutta method is employed to advance the solution in time. In the typical characteristic-wise finite-difference WENO scheme, the are usually reconstructed within the local characteristic fields. Let us take the matrix to be the Roe-average Jacobian at . We denote by , (column vector) and (row vector) the th eigenvalue, and the right and left eigenvectors of , respectively. The physical fluxes and conservative variables on the respective reconstruction stencil are mapped to the characteristic field by the characteristic-projection step

(3)

where and . For each component of the characteristic variables, the corresponding split numerical fluxes are constructed by

(4)

where for a Roe flux (RF). Alternatively one can use , where represents the entire computational domain for a Lax-Friedrichs flux (LF) or the neighborhood of for a local Lax-Friedrichs flux (LLF). The WENO reconstruction gives

(5)

where are 3rd-order reconstructions from upwind 3-point stencils, and are WENO weights defined by

(6)

where , and are optimal weights. These optimal weights generate the 5th-order upwind scheme, by which the numerical flux is reconstructed from a 5-point stencil. prevents division by zero, or is chosen to adjust the distinct weights, and are the smoothness indicators. We refer to Jiang and Shu [4] for the details of WENO reconstruction. The numerical flux in each characteristic field is then computed by

(7)

At last, this numerical flux is projected back to the physical space by

(8)

It can be found that the major part of floating-point operations in the characteristic-wise WENO scheme are due to the characteristic-projection step (matrix-vector product) of Eq. (3), and the computation of the WENO weights in Eq. (6).

2.2 Hybrid WENO scheme

In order to decrease the number of floating-point operations, we propose to hybridize the WENO scheme with its optimal linear scheme by the characteristic-decomposition approach. A straightforward implementation is to switch between Eq. (7) and the corresponding flux of the optimal linear scheme based on the characteristic-variable projection. Here, however, we consider a component-wise reconstruction for the optimal linear scheme to further save the computational effort for the characteristic-projection step.

2.2.1 Optimal linear scheme with component-wise reconstruction

In this formulation, the approximated physical fluxes and the differences of approximated conservative variables at are computed component-by-component as

(9)

With the optimal linear scheme of the 5th-order WENO scheme, one has

(10)
(11)

Note that Taylor-series expansion of Eq. (11) suggests

(12)

Then and are projected onto the characteristic field by

(13)

For each component of the characteristic variables, the corresponding numerical flux is constructed by

(14)

where . Note that such choice of RF fluxes leads to less numerical dissipation compared to that of LLF or LF fluxes. At last, the numerical flux obtained in each characteristic field is projected back to the component space by Eq. (8). Note that, compared to the characteristic-wise WENO scheme, the linear scheme omits the computation of the WENO weights in Eq. (6), and decreases by the characteristic-projection operations of Eq. (3). Also note that such properties of the present single-component characteristic-projection allows the hybridization applied within the characteristic field, as shown in the next section, as in Refs. [13, 7, 17] with much less computational effort.

2.2.2 Hybridization and discontinuity detector

Since the linear scheme is numerically unstable for solutions with discontinuities, it should be switched on only in smooth regions of the solution according to a discontinuity detector. A less problem-independent discontinuity detector than those in the trouble-cell hybridization [1, 11, 10] are usually based on non-dimensional measurement on the characteristic variables. A usual way to achieve this is using the ratio between low order undivided differences [13, 7], or the ratio between high and low order undivided differences [10]. Since undivided difference corresponds to approximation of a spatial derivative, this type of discontinuity detectors degenerate near critical points with zero low order derivatives as the denominators approach zero and switches on the WENO scheme. To decrease the over-dissipation associated with the degeneration, dimensional parameters must be introduced and tuned for different problems.

In this paper, we achieve a general effective measure in a different way. By noticing that the characteristic variables have the dimension of density we define a non-dimensional discontinuity detector by

(15)

where is the Roe-average density of . Note that there may be different forms of left and right eigenvectors which lead to different dimension of the characteristic variables. In this case, the expected dimension can be achieved by multiplying or diving a factor to corresponding eigenvectors. From Eq. (12) it can be found that also is related to the 5th-order derivatives of the characteristic variables due to the linearized characteristic projection. Since the optimal 5th-order linear scheme reconstructs with 4th-degree polynomials and is not able to resolve 5th or higher order derivatives, is actually a measure on the resolvability of the linear scheme. If there is no vacuum in the flow (which is almost always the case), the denominator of does not degenerate due to the finite value of .

With a given threshold expression

(16)

where is a constant, is the characteristic length scale of the problem and is an integer, the hybrid scheme can be constructed as follows: for each component of the characteristic variables the numerical flux is obtained by linear flux of Eq. (14) if ; otherwise it is obtained by WENO flux of Eq. (7). Note that neither Eq. (15) or Eq. (16) introduces any dimensional parameter. As will be shown in the next section also, is introduced to achieve convergence for capturing discontinuities. Assume Eq. (15) measures on the cell-face at next to a discontinuity at , the hybrid scheme tolerates, i.e. by applying the optimal linear scheme, a small overshoot or undershoot scaling with . Due to the self-similar property of the discontinuity, is unchanged with increasing resolution. In this case, convergence requires that decreases, where determines the rate, with increasing resolution.

Alo note that the above approach can be directly applied to multiple dimensions by a dimension-by-dimension approach, and can be very easily extended to higher-order or modified WENO schemes. The only modification for other WENO schemes is that Eq. (9) should be computed with the corresponding optimal linear schemes.

3 Numerical examples

The following numerical examples are provided to illustrate the potential of the proposed hybrid WENO scheme. The governing equations are the one- and two-dimensional compressible Euler or Navier-Stokes equations. While the original 5th-order WENO is denoted as WENO-5, the present hybrid scheme is denoted as H-WENO-5. The 3rd-order TVD Runge-Kutta scheme is used for time integration [14]. In order to achieve an objective measurement on the computational cost, the local efficiency on each cell-face through the entire computation time is defined by , where is the number of operations with WENO flux and is the number of operations with linear flux. The overall computational efficiency is obtained by averaging over the entire computational domain. In order to show the general effective discontinuity detector of Eq. (15), We set the parameters and in Eq. (16) for all numerical examples (see also a study of the sensitivity of in Sec. 3.3). If not mentioned otherwise, all the computations are carried out with a CFL number of 0.6.

3.1 Propagation of broadband sound waves

This problem, taken from Sun et al. [16], corresponds to the propagation of a sound wave packet which contains acoustic turbulent structure with various length scales. The initial condition is

where is a random number between 0 and 1 with uniform distribution, , , is the speed of sound and

is the energy spectrum which reaches its maximum at . A periodic boundary condition is applied at and . Computations have been carried out on a -point grid using a CFL number of 0.2 for one period of time.

Figure 1: Propagation of broadband sound waves computed on a 128 points grid: (a) pressure and (b) local computational efficiency distribution.

The numerical results in Fig. 1a show that the present methods gives considerably better resolved sound waves than WENO-5, especially for the regions near the critical points, at which previous discontinuity detectors, such as that used in Sun et al. [16], can switch on WENO scheme if no dimensional parameter is introduced. As shown in Fig. 1b, no WENO flux is switched on during the entire computation, indicating that the present hybrid scheme recovers the optimal linear scheme, therefore is able to achieve the formal 5th-order accuracy even near the critical points, where WENO-5 degenerates the order of accuracy and introduces extra numerical dissipation.

3.2 Shock-tube problems

Here, we show that the proposed scheme H-WENO-5 passes the shock-tube test problems: the Sod problem (Sod 1978), the Lax problem (Lax 1954) and the 1-2-3 problem (Einfeldt et al. 1991). For the Sod problem, the initial condition is

and the final time is . For the Lax problem, the initial condition is

and the final time is . For the 123 problem, the initial condition is

and the final time is . We examine the numerical solution with difference resolutions to study the performance of the present method with increasing resolution.

Figure 2: Shock-tube problems: (a, b) Sod problem; (c, d) Lax problem; (e, f) 1-2-3 problem.

Figures 2a, c and e gives the density distributions computed on a 200-point grid and local computational efficiency for three grid resolutions. It can be observed that good agreement between the present solution and that of WENO-5. As shown in Fig. 2b, d and f, the computational efficiency increases with resolution. It can also be observed that the WENO-5 is essentially not switched on in the smooth regions, even those with 1st or high-order singularities near the fringe of expansion waves. For all three problems, the overall computational efficiency achieves 99.1% for the lowest resolution and 99.7% for the highest resolution. Note that, the computational efficiency achieved here is even much higher than those obtained by the trouble-cell hybridization using several different discontinuity detectors [10]. Also note that the hybrid scheme introduces very small overshoots near the discontinuities, as shown in Fig. 2a and c. If the threshold is chosen in the form of Eq. (16), as shown by Fig. 3,

Figure 3: Sod problem computed with (a) and (b) for different grid resolutions.

these overshoots vanish with increasing resolution and the captured discontinuities converge to their exact solution. However, if a constant threshold is used, these overshoots remain at approximately the same magnitude therefore no convergence is achieved.

3.3 Interacting blast waves

We consider a two-blast-wave interaction problem, which is taken from Woodward and Colella [2] . The initial condition is

and the final time is . The reflective boundary condition is applied at both and . We examine the numerical solution on a 400-point grid with the parameters , and . The reference ”exact” solution is a high-resolution solution on a -point grid computed by WENO-5.

Figure 4: Interacting blast waves: (a) density profile with a closer view; (b) local computational efficiency.

Figure 4 give the profiles of density and local computational efficiency. Again, good agreement with the reference solutions is observed and the computational efficiency is much higher than those obtained by the trouble-cell hybridization [10]. Due to smaller numerical dissipation H-WENO-5 shows an slightly improved solution compared to WENO-5. As shown in Fig. 4b, the computational efficiency is only weakly dependent on the choice of . Actually, by increasing 100 times, the extra gain of overall efficiency is less than 0.7% (from 98.3% to 99.0%). Note that, since the computation (results not shown here) is still stable even with , the numerical robustness of the present method is also weakly depends on the choice of .

3.4 Shock-density wave interaction

We consider a shock density-wave interaction problem [14]. The initial condition is set by a Mach 3 shock interacting with a perturbed density field

and the final time is . A zero-gradient boundary condition is applied at and . We examine the numerical solution on 200-point and 400-point grids, and the reference ”exact” solution is a high-resolution solution on a -point grid computed by WENO-5.

Figure 5: Shock-density-wave interaction: density and velocity profiles on (a) a 200-point grid and (b) a 400-point grid.

Fig. 5a show the calculated density and velocity profile. It can be observed that, when the grid resolution is low, H-WENO-5 gives considerably better resolved density waves behind the shock than WENO-5, which is too dissipative. By considering all cases discussed in previous sections, it can be find that the present method gains more significant improvement in smoothed region than near the discontinuity.

3.5 Double Mach reflection a strong shock

We consider the problem from Woodward and Colella [2] on the double Mach reflection of a strong shock. A Mach 10 shock in air is reflected from the wall with incident angle of . For this case, the initial condition is

and the final time is . The computational domain of this problem is . Initially, the shock extends from the point at the bottom to the top of the computational domain. Along the bottom boundary, at , the region from to is always assigned post-shock conditions, whereas a reflecting wall condition is set from to . Inflow and outflow boundary conditions are applied at the left and right ends of the domain, respectively. The values at the top boundary are set to describe the exact motion of a Mach 10 shock. We examine the numerical solution with three resolutions on grids of , and points.

Figure 6 shows contours of density and local computational efficiency computed on a grid.

Figure 6: Double-Mach reflection of a Mach 10 shock wave at computed on grid: (a) 40 density contours from 1.88783 to 20.9144, (b) 40 contours of local computational efficiency from 0.5 to 1. Note that in the most of the computational domain the calls for WENO-5 scheme are very few, i.e. less than 1.3% which gives the first contour in (b).

It is observed a good agreement with the results of Kim and Kwon [7] (their Figs. 12 and 13) computed with fine-tuned parameters at the same resolution. Compared to the results in Zhou et al. [17] computed with fine-tuned parameters and higher grid resolution, much more fine-scale structures are obtained by the present scheme. Note that, as shown in Fig. 6b, the WENO flux is seldom calculated in the entire computation domain, except near the end of the main reflection wave. The overall computational efficiencies respectively for three resolutions are , and , which decreases with increasing grid resolution. Such computational efficiencies indicate that, the computational cost of the WENO scheme is essentially negligible.

3.6 Viscous shock tube problem

We consider the two-dimensional viscous flow problem in a square shock tube with unit height and insulated walls [3]. In this problem, the propagation of the incident shock wave and contact discontinuity lead to a thin boundary layer. After its reflection on the right wall, the shock wave interacts with this boundary layer and results a separation region and the formation of a typical ”-shape like shock pattern”. The initial condition is

The fluid is assumed as ideal gas with and constant dynamics viscosity and Prandtl number and satisfying the Stokes assumption. If the reference values are chosen as the initial speed of sound, unit density and unit length, the Reynolds number is 200. By applying the symmetry condition at the upper boundary, only the lower half domain is actually computed. For other boundaries, the no-slip and adiabatic wall conditions are applied. The viscous and heat transfer is calculated by 6th-order accuracy, Here, we examine the numerical solution with two resolutions on grids of 250125 and 500 points.

As shown by the density contours on the coarse grid in Fig. 7a, the computed result, especially with respect to the shape and the tilting angle of the primary vortex, is already considerably better than those obtained by the WENO-5, the 6th-order artificial compression method (ACM) with wavelet based filter scheme [15] (their Fig. 4) and the 5th-order multi-dimensional limiting process (MLP) cocombined with M-AUSMPW+ scheme [8] (their Fig. 41) with the same grid resolution.

Figure 7: Viscous shock tube problem at : 20 density contours from 15 to 125 computed on (a) a grid and (b) a grid.

Grid convergence is indicated for the results computed on the fine grid because no notable differences, except the near shock region, can be identified upon comparison with the grid-converged density contours of Sjögreen and Yee [15] (their Fig. 2) obtained with much higher grid resolutions. Note that there are small disturbances near the shock wave in the low resolution results (see. Fig. 7a). These small disturbances are numerical errors rather than instabilities, since they do not propagate to other regions and their magnitudes decreases by increasing resolution (see. Fig. 7b). Also note that the obtained overall computational efficiency for the two resolutions are , , which again imply that the WENO flux is actually hardly utilized throughout the entire computation.

4 Concluding remarks

We propose a simple hybrid WENO scheme to increase computational efficiency and decrease numerical dissipation. Based on characteristic-wise hybridization, the scheme switches the numerical flux of each characteristic variables between that of the WENO scheme and its optimal linear scheme according to a discontinuity detector measuring the resolvability of the linear scheme. As shown by a number of numerical examples the overall computational efficiency, measured by the fraction of linear flux used, is always higher than , shows that the hybridization increases the computational efficiency greatly. Also, by choosing a general effective threshold for the discontinuity detector, compared to the original WENO scheme, considerable lower numerical dissipation is achieved in smooth region without compensating the robustness for capturing discontinuities. Note that the method can be applied for multiple-species flows with even higher computational efficiency due to the serious overhead of local characteristic projection and WENO weights calculation for each species.

References

  • [1] N.A. Adams and Shariff K. A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems. J. Comput. Phys., 127:27–51, 1996.
  • [2] P. Colella and P.R. Woodward. The Piecewise Parabolic Method (PPM) for gas-dynamical simulations. J. Comput. Phys., 54(1):174–201, 1984.
  • [3] V. Daru and C. Tenaud. Evaluation of tvd high resolution schemes for unsteady viscous shocked flows. Computers & fluids, 30(1):89–113, 2000.
  • [4] G.S. Jiang and C.W. Shu. Efficient implementation of weighted ENO schemes. J. Comput. Phys, 126:202–228, 1996.
  • [5] E. Johnsen, J. Larsson, A.V. Bhagatwala, W. H. Cabot, P. Moin, B. J. Olson, P. S. Rawat, S. K. Shankar, B. Sjögreen, H.C. Yee, X. Zhong, and S. K. Lele. Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves. J. Comput. Phys., 229(4):1213–1237, 2010.
  • [6] S. Kawai and SK Lele. Localized artificial diffusivity scheme for discontinuity capturing on curvilinear meshes. J. Comput. Phys., 227(22):9498–9526, 2008.
  • [7] D. Kim and J.H. Kwon. A high-order accurate hybrid scheme using a central flux scheme and a weno scheme for compressible flowfield analysis. J. Comput. Phys., 210(2):554–583, 2005.
  • [8] K.H. Kim and C. Kim. Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows:: Part II: Multi-dimensional limiting process. J. Comput. Phys., 208(2):570–615, 2005.
  • [9] J. Larsson and B. Gustafsson. Stability criteria for hybrid difference methods. J. Comput. Phys., 227(5):2886–2898, 2008.
  • [10] G. Li and J. Qiu. Hybrid weighted essentially non-oscillatory schemes with different indicators. J. Comput. Phys., 229(21):8105–8129, 2010.
  • [11] S. Pirozzoli. Conservative hybrid compact-weno schemes for shock-turbulence interaction. J. Comput. Phys., 178(1):81–117, 2002.
  • [12] S. Pirozzoli. Numerical methods for high-speed flows. Ann. Rev. Fluid Mech., 43:163–194, 2011.
  • [13] Y.X. Ren, M. Liu, and H. Zhang. A characteristic-wise hybrid compact-weno scheme for solving hyperbolic conservation laws. J. Comput. Phys., 192(2):365–386, 2003.
  • [14] C.W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II. J. Comput. Phys., 83(1):32–78, 1989.
  • [15] B. Sjögreen and HC Yee. Grid convergence of high order methods for multiscale complex unsteady viscous compressible flows. J. Comput. Phys., 185(1):1–26, 2003.
  • [16] Zhen-Sheng Sun, Yu-Xin Ren, Cédric Larricq, Shi-ying Zhang, and Yue-cheng Yang. A class of finite difference schemes with low dispersion and controllable dissipation for dns of compressible turbulence. J. Comput. Phys., 230(12):4616–4635, 2011.
  • [17] Q. Zhou, F. He, and MY Shen. A family of efficient high-order hybrid finite difference schemes based on weno schemes. International Journal of Computational Fluid Dynamics, 26(4):205–229, 2012.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
84166
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description