## Abstract

In this article, we give a brief overview on high-order accurate shock capturing schemes with the aim of applications in compressible turbulence simulations. The emphasis is on the basic methodology and recent algorithm developments for two classes of high-order methods: the weighted essentially non-oscillatory and discontinuous Galerkin methods.

## 1. Introduction

In this overview article, we are concerned with the numerical simulation of compressible turbulence. This includes both direct numerical simulation (DNS), for which the aim is to resolve the compressible Navier–Stokes equations to all the relevant scales, and various turbulence models, for which the smaller scales are not resolved but rather are modelled.

The numerical methods we are discussing in this article are the high-order accurate methods. High-order accuracy refers to the fact that the local truncation error, which is the error incurred when using the scheme to evolve from the exact solution of the partial differential equation (PDE) for one time step, is proportional to *N*^{−p} for large *p*, where *N* is the total number of degrees of freedom (e.g. the total number of mesh points for a finite difference scheme). Usually, this can be achieved only if the exact solution of the PDE is smooth. That is, the constant coefficient to the *N*^{−p} term in the local truncation error is usually related to some high-order derivatives of the exact solution. If we are concerned with the inviscid Euler equations of compressible gas dynamics, the exact solutions will be no longer smooth, since shocks, contact discontinuities, etc., will appear. However, the solutions are still in general *piecewise smooth*, hence the local truncation error will still be of small size proportional to *N*^{−p} in the smooth region. Therefore, such numerical methods are still called high-order schemes.

Spectral methods [1–3] are the ultimate high-order schemes, with truncation errors smaller than *N*^{−p} for any *p* as long as the exact solution has enough smooth derivatives. In fact, the truncation errors of a spectral method is exponentially small, e^{−αN} for some *α*>0, if the exact solution of the PDE is analytic. Therefore, spectral methods are widely used in turbulence simulations, especially in incompressible turbulence simulations for which it is still possible to resolve all the scales for reasonable Reynolds numbers. There are however some limitations of spectral methods:

Spectral methods for non-periodic problems or for problems in complex geometry are more difficult to design.

Spectral methods are global methods. Therefore, if the solution is not smooth in some part of the computational domain, the accuracy is degenerated throughout the complete domain. The method may also become unstable in such a situation.

The first limitation can be mitigated by the design of polynomial (rather than trigonometric polynomial) spectral methods, e.g. Chebyshev method, which involve non-uniform collocation points, and the spectral element method or multi-domain spectral method [1,2,4]. The second limitation is more difficult to overcome. A typical approach is to add filtering or spectral vanishing viscosities, which would enhance stability while improving accuracy in smooth regions [3].

When the spectral method is not convenient to use (non-uniform mesh, small time step for stability, etc.), a good alternative is the high-order compact schemes [5], which can give ‘spectral-like’ resolution for many wave problems. Compact schemes are very popular choices for turbulence simulations.

Both the spectral method and the compact schemes are linear schemes, that is, they are linear when applied to a linear wave PDE
1.1
where *a* is a constant. Godunov [6] proved that *all* linear schemes are either first-order accurate or oscillatory (that is, it does not satisfy a maximum principle when applied to the wave equation (1.1)). Therefore, the spectral method and compact schemes are all oscillatory when applied to problems with discontinuities (or with sharp gradients). This is not only unpleasant visually, but also leads to global deterioration of errors and even to nonlinear instability (blowup of the numerical solution). For this reason, various nonlinear schemes (which are nonlinear even for the linear PDE (1.1)) have been designed in the literature. In this article, we mainly discuss two classes of such methods: the weighted essentially non-oscillatory (WENO) finite difference and finite volume methods, and the Runge–Kutta discontinuous Galerkin (DG) finite-element methods with nonlinear limiters, in the perspective of application to compressible turbulence simulations.

Compressible Navier–Stokes equations with high Reynolds numbers belong to the class of convection-dominated problems. The main challenge to the design of numerical schemes for such problems, especially for hyperbolic conservation laws, including the Euler equations, is the presence of discontinuities (such as shocks and contact discontinuities) or sharp gradient regions. Among the successful numerical methods in such a situation, we mention the first-order schemes such as the Godunov scheme [6] or the Roe scheme [7]. These schemes can resolve the discontinuities monotonically without spurious numerical oscillations; however, they contain relatively large numerical dissipation in the smooth part of the solution, hence are not suitable for turbulence simulations. These first-order schemes are, however, used as building blocks for the design of higher-order schemes. The first nonlinear, non-oscillatory scheme higher than first-order accuracy seems to be that designed by Kolgan in 1972; see [8]. The second-order ‘high-resolution’ schemes designed in the 1970s and 1980s, represented by the MUSCL schemes [9], total variation diminishing (TVD) schemes [10] and piecewise parabolic method (PPM) schemes [11], are second-order accurate in smooth regions (and can be higher order in smooth, monotone regions away from smooth extrema), and can resolve discontinuities monotonically with a sharper transition than first-order schemes. These high-resolution schemes are very popular in applications, including in compressible turbulence simulations, especially for various turbulence models. However, in some of the situations, for example DNS of turbulence, it is desirable to have non-oscillatory schemes that are higher-order accurate. The WENO finite difference and finite volume schemes and the Runge–Kutta DG finite-element methods with nonlinear limiters belong to this class.

## 2. Time discretization

The discussion in this paper is mainly concentrated on spatial discretization. However, it is equally important to obtain stable, high-order time discretization. In this section, we briefly discuss the issue of time discretization.

There are mainly two approaches to discretize the time variable. The more popular method is to first discretize in space, obtain a large ordinary differential equation (ODE) in the time variable, then discretize this ODE using one of the standard ODE solvers. This approach is usually referred to as the method of lines. While in principle any ODE solver can be used, there are studies in the literature to obtain ODE solvers with special properties when coupled with suitable spatial discretizations. For example, the so-called TVD time discretization, first developed in Shu [12] for the multi-step methods and in Shu & Osher [13] for the Runge–Kutta methods, is designed to maintain stability in the total variation semi-norm, or any other norm or semi-norm, of the first-order Euler forward method with the same spatial discretization. The most popular TVD time discretization is the third-order Runge–Kutta method in Shu & Osher [13]:
2.1
where *L*(*u*) is the spatial discretization. Higher-order versions of TVD Runge–Kutta or multi-step types are also available. These types of time discretizations are also known as the strong stability preserving (SSP) methods in the later literature. We refer to Gottlieb *et al.* [14–16] for more details.

While Runge–Kutta or multi-step methods are easy to implement, they are more difficult to obtain for very high-order accuracy. This is especially the case for SSP time discretizations, which have various order barriers [15,16]. For turbulence simulations, one would in principle be interested in time discretizations that are of the same high-order accuracy as the spatial discretization. This leads to the second type of time discretization, which is referred to as the Lax–Wendroff type time discretization. In this approach, one first performs a Taylor expansion in time to the desired accuracy, then one uses the PDE repeatedly to replace all the time derivatives by spatial derivatives, and then one discretizes these spatial derivatives to suitable order of accuracy. The traditional second-order Lax–Wendroff scheme [17] uses this type of time discretization to obtain second-order accuracy both in space and in time. The first essentially non-oscillatory (ENO) scheme of Harten *et al.* [18] also uses this type of time discretization. Later studies on this type of time discretization [19–22] discuss both the WENO and DG spatial discretizations. This type of time discretization has a good potential in future applications of turbulence simulation.

## 3. Weighted essentially non-oscillatory finite difference and finite volume schemes

In this section, we give an overview of the WENO finite difference and finite volume schemes, aiming at numerical simulation of compressible turbulence. We start with a brief introduction to ENO schemes, from which WENO schemes were developed.

### (a) A brief introduction to essentially non-oscillatory schemes

In most numerical simulations, including turbulence simulations, the generic solutions are piecewise smooth. That is, there are isolated curves (in two dimensions) or surfaces (in three dimensions) on which the solutions are discontinuous; however, between them the solutions are smooth (having bounded high-order derivatives).

High-order polynomial interpolations or reconstructions may have difficulties approximating such functions if care is not taken. For example, consider the simple problem that we know the point values *u*_{i} of a piecewise smooth function on the uniform grids *x*=*x*_{i}, and we would like to obtain approximations to the function at middle points , or the derivative of the function at *x*=*x*_{i}. The typical procedure is to take a ‘stencil’ containing several grid points, for example {*x*_{i−1},*x*_{i},*x*_{i+1},*x*_{i+2}}, and obtain an interpolation polynomial *p*(*x*) of suitable degree (in this case of degree at most three) such that *p*(*x*_{j})=*u*_{j} at *j*=*i*−1,*i*,*i*+1,*i*+2. Then, we can take the approximation to the function at the middle point as *p*(*x*_{i+1/2}), or the approximation to the derivative of the function at *x*=*x*_{i} as *p*′(*x*_{i}). If the function is completely smooth in the stencil (the region covered by all the grid points in the stencil), such approximations are very accurate. However, if there is a discontinuity of *u*(*x*) within the stencil, say between *x*_{i+1} and *x*_{i+2}, then the interpolation polynomial *p*(*x*) is a very poor approximation to the function *u*(*x*). The idea of ENO approximation [18] is to choose the stencil adaptively. That is, in the situation that there is a discontinuity of *u*(*x*) between *x*_{i+1} and *x*_{i+2}, we will choose a stencil that avoids this discontinuity, for example, we could choose {*x*_{i−2},*x*_{i−1},*x*_{i},*x*_{i+1}}. Now the interpolation polynomial will be accurate since the stencil does not contain any discontinuity. Of course, in practice, there is no sharp boundary between a discontinuity and a sharp gradient region, hence the ENO procedure in Harten *et al.* [18] avoids explicit determination of discontinuities. Rather, it uses the information of the given data {*u*_{i}} and their divided differences to automatically choose a locally ‘smoothest’ stencil at every location where interpolation or reconstruction is needed. The ENO procedure is uniformly high-order accurate and essentially non-oscillatory for piecewise smooth functions. Even though theory is scarce for ENO schemes, in the recent work [23] there are very interesting results of ENO schemes being sign-preserving and stable in the sense of controlling the ratio of reconstruction jumps and jumps in the given data.

### (b) History and overview of weighted essentially non-oscillatory schemes

The first WENO scheme was introduced in 1994 by Liu, Osher & Chan [24] in their pioneering paper, in which a third-order accurate finite volume WENO scheme was designed. In 1996, Jiang & Shu [25] provided a general framework to construct arbitrary-order accurate finite difference WENO schemes, which are more efficient for multi-dimensional calculations. Most applications of WENO schemes in the literature use the fifth-order accurate WENO scheme designed in Jiang & Shu [25].

We will briefly describe the fifth-order finite difference WENO scheme in Jiang & Shu [25] for the scalar one-dimensional conservation law
3.1
We further assume the wind direction is from the left to right, namely *f*′(*u*)≥0, to simplify the presentation. The spatial derivative term in (3.1) is approximated by a conservative flux difference
3.2
where the numerical flux is obtained by the procedure given below. We use the shorthand notation *f*_{i}=*f*(*u*_{i}) where *u*_{i} is the approximation to the grid value of the solution *u*(*x*_{i},*t*). First we denote the following three third-order numerical fluxes:
The numerical flux in (3.2) is obtained by a convex combination of these third-order fluxes
3.3
where the nonlinear weights *ω*_{m} are given by
3.4
Here *ϵ* is a parameter to avoid the denominator becoming zero and is usually taken as *ϵ*=10^{−6}. The linear weights *γ*_{l} are given by
3.5
and the smoothness indicators *β*_{l} are given by
These smoothness indicators measure the smoothness of the function in the relevant stencils. In smooth regions *β*_{l} are of the size *O*(*Δx*^{2}); however, if the stencil contains a discontinuity, *β*_{l} can be *O*(1). This in turn implies that the nonlinear weights *ω*_{m} defined in (3.4) will be much smaller for the stencil containing a discontinuity when compared with the one for the stencil in which the solution is smooth. The essentially non-oscillatory property is a consequence of this fact. We can also show that uniform high-order accuracy is maintained in smooth regions.

After the spatial discretization, we obtain a method of lines ordinary differential equation (ODE) 3.6 This ODE system can be discretized in time by the TVD Runge–Kutta time discretization methods [13] described in §2, for example, the third-order version (2.1). Of course, other types of time discretization, such as multi-step and Lax–Wendroff type time discretizations discussed in §2 can also be applied.

The finite difference WENO schemes can be easily generalized to multiple dimensions and to systems, and they perform well for problems containing both strong discontinuities and complex smooth regions features. They are therefore good candidates for compressible turbulence simulations. Higher-order versions of the finite difference WENO schemes (order 7 to 11) are derived in Balsara & Shu [26]. For more details of WENO schemes, we refer to the Lecture Notes [27] and the review paper [28].

### (c) Development of weighted essentially non-oscillatory schemes tailored to compressible turbulence simulation

When applied to turbulence simulations, one would like to use very high-order WENO schemes, for example those in Balsara & Shu [26]. Also, we mention research performed in the following two aspects.

#### (i) Dissipation and dispersion optimized weighted essentially non-oscillatory schemes

Turbulence simulation involves wave propagation over a long time with a relatively coarse mesh (that is, the number of grid points per wave is not large). Therefore, dissipation and dispersion errors in such situations could be very large, even for high-order schemes. There is a systematic approach [29] to modify the coefficients of the finite difference approximation, with the objective of increasing the amplitude or phase resolution at the price of lowering the formal order of accuracy for the same stencil. Even though such dissipation or dispersion optimized finite difference schemes are less accurate for asymptotically small mesh sizes (that is, when the number of points per wave is large), they are better than the usual finite difference schemes on the same stencil for marginally resolved waves (that is, when the number of points per wave is small).

Even though WENO schemes are nonlinear, they can also be adapted by such dissipation and dispersion optimized techniques to enhance their resolution power for marginally resolved waves. Such WENO schemes are therefore more suitable for compressible turbulence simulations. For efforts in this direction, we mention, for example, the work in [30–33]. Other efforts to improve the performance of WENO schemes and their simulation for turbulence simulations can be found in the studies of Taylor & Martin [34,35] and Chaudhuri *et al.* [36].

#### (ii) Weighted compact schemes

As mentioned in §1, compact schemes are very popular in turbulence simulations because of their relatively small dissipation and dispersion errors. It would be desirable to combine the advantages of compact schemes in their small dissipation and dispersion errors and of the WENO schemes in their shock capturing capability. This leads to various efforts in designing compact WENO schemes [37–40].

### (d) Weighted essentially non-oscillatory schemes for simulating shock–vortex interactions

Shock waves and vortices are two basic elements of compressible flow. The interaction between them is a common phenomenon and is very important in many applications, such as supersonic mixing layers, supersonic jets and combustion instability. In particular, because a number of shock waves and vortices coexist in the complicated supersonic turbulence flow, the interaction of a shock wave and vortices can be seen as a simplified model of shock–turbulence interaction, which is one of the major sources of noise and has received increasing attention.

WENO schemes are particularly suited for the simulation of shock–vortex interactions because of their simultaneous capability of high-order solution and non-oscillatory shock transition. Examples of such applications can be found in earlier studies [41–47].

## 4. Runge–Kutta discontinuous Galerkin methods with nonlinear limiters

In this section, we give an overview of the DG finite-element method, aiming at numerical simulation of compressible turbulence.

DG methods are a class of finite-element methods using completely discontinuous basis functions, which are usually chosen as piecewise polynomials. Since the basis functions can be completely discontinuous, these methods have a flexibility that is not shared by typical finite-element methods, such as the allowance of arbitrary triangulation with hanging nodes, complete freedom in changing the polynomial degrees in each element independent of that in the neighbours (*p* adaptivity), and extremely local data structure (elements only communicate with immediate neighbours regardless of the order of accuracy of the scheme) and the resulting embarrassingly high parallel efficiency (usually more than 99% for a fixed mesh, and more than 80% for a dynamic load balancing with adaptive meshes that change often during time evolution) [48,49].

The first DG method was introduced in 1973 by Reed & Hill [50], in the framework of neutron transport, i.e. a time-independent linear hyperbolic equation. A major development of the DG method was carried out by Cockburn *et al.* in a series of papers [51–54], in which they have established a framework to easily solve *nonlinear* time-dependent hyperbolic conservation laws, such as the Euler equations of gas dynamics, using explicit, nonlinearly stable high-order Runge–Kutta time discretizations [13] and DG discretization in space with exact or approximate Riemann solvers as interface fluxes and total variation bounded (TVB) nonlinear limiters [55] to achieve non-oscillatory properties for strong shocks. The viscous terms (e.g. those in the Navier–Stokes equations) can be approximated within the DG framework, e.g. by the local DG methods [56].

We again use the following one-dimensional conservation law
4.1
to describe the general idea of the DG method. We assume the following mesh to cover the computational domain, consisting of cells *I*_{i}=[*x*_{i−1/2},*x*_{i+1/2}], for 1≤*i*≤*N*. We denote
The mesh does not need to be uniform but is usually assumed to be regular, namely there is a constant *c*>0 independent of *h* such that
We define a finite-element space consisting of piecewise polynomials
4.2
where *P*^{k}(*I*_{i}) denotes the set of polynomials of degree up to *k* defined on the cell *I*_{i}. The semi-discrete DG method for solving (4.1) is defined as follows: find the unique function such that, for all test functions and all 1≤*i*≤*N*, we have
4.3
Here, is the numerical flux, which is a single-valued function defined at the cell interfaces and in general depends on the values of the numerical solution *u*_{h} from both sides of the interface
We use the so-called monotone fluxes from finite difference and finite volume schemes for solving conservation laws. For example, if the wind is always from left to right, or *f*′(*u*)≥0, we can take .

After the spatial discretization, we obtain a method of lines ODE in the form of (3.6), with the only difference that *u*_{j} is now a vector containing all the degrees of freedom in the cell *I*_{j} (e.g. the coefficients when the solution polynomial is written out using a local basis). This ODE can again be solved by the TVD or SSP time discretization techniques, such as the third-order Runge–Kutta method (2.1). The time-step restriction (CFL numbers) can be found in the review paper [57]. The CFL numbers are usually much smaller than those for finite difference or finite volume schemes. However, since there are many degrees of freedom in each cell, the DG method can typically use much coarser meshes to obtain the same error magnitudes when compared with finite volume schemes of the same order of accuracy. See Zhou *et al.* [58] for some comparison results.

The DG method can be easily generalized to multiple dimensions and to systems, over arbitrary unstructured meshes. Unlike the finite difference or finite volume methods described in §3, the DG method has a provable cell entropy inequality for the square entropy, for all nonlinear scalar equations or nonlinear symmetric hyperbolic systems, in any space dimensions on any triangulations, for any order of accuracy, without any assumption on the smoothness of the solution [59,60]. This cell entropy inequality implies *L*^{2} stability. The DG method therefore works well for smooth problems, or for problems containing weak shocks. However, if the problem contains strong shocks, the *L*^{2} stability is not enough to control the spurious numerical oscillations. Some nonlinear limiters would be needed, for example, the TVB limiter [55]. So far, the study of limiters that are easy to implement, robust and accurate is still ongoing. This can be considered as a bottleneck for the DG method when applied to compressible flow simulations with strong shocks. Recently, there has been effort to use the WENO methodology as limiters for the DG methods [61].

For more details about the DG methods, we refer to the survey paper [62], other papers in those Lecture Notes, the review paper [57], and the three recent special journal issues on DG methods [63–65]. Recently, there are at least four books [66–69] covering different aspects of the DG method. See also the lecture notes [70].

Efforts in the applications of DG methods in compressible turbulence simulations can be seen in earlier studies [71–73]. There is an excellent potential for more applications of the DG methods to compressible turbulence simulations.

## 5. Conclusions and future work

In this short overview paper, we have summarized the basic algorithm methodology of the high-order accurate WENO schemes and the DG finite-element methods, with the aim of applications in compressible turbulence simulations. Current activities in such applications are briefly reviewed.

Both the WENO schemes and the DG methods are good choices for compressible turbulence simulations. For problems on simple geometry and with possibly strong shocks, the high-order finite difference WENO schemes are suitable choices. Improvement of such schemes in terms of a further reduction in dissipation and dispersion errors would be interesting future work. For problems on complex geometry, or problems with local structures requiring adaptive meshes, the high-order DG methods are suitable choices. Further investigation on finding more robust, efficient and accurate nonlinear limiters for DG methods in the presence of strong shocks, especially for highly stretched meshes, would be interesting future work.

## Acknowledgements

Research is supported in part by ARO grant nos W911NF-08-1-0520, W911NF-11-1-0091 and NSF grant no. DMS-1112700.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales I’.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.