## Abstract

We analyse the effect of second- and fourth-order accurate central finite-volume discretizations on the outcome of large eddy simulations of homogeneous, isotropic, decaying turbulence at an initial Taylor–Reynolds number *Re*_{λ}=100. We determine the implicit filter that is induced by the spatial discretization and show that a higher order discretization also induces a higher order filter, i.e. a low-pass filter that keeps a wider range of flow scales virtually unchanged. The effectiveness of the implicit filtering is correlated with the optimal refinement strategy as observed in an error-landscape analysis based on Smagorinsky's subfilter model. As a point of reference, a finite-volume method that is second-order accurate for both the convective and the viscous fluxes in the Navier–Stokes equations is used. We observe that changing to a fourth-order accurate convective discretization leads to a higher value of the Smagorinsky coefficient *C*_{S} required to achieve minimal total error at given resolution. Conversely, changing only the viscous flux discretization to fourth-order accuracy implies that optimal simulation results are obtained at lower values of *C*_{S}. Finally, a fully fourth-order discretization yields an optimal *C*_{S} that is slightly lower than the reference fully second-order method.

## 1. Introduction

Direct numerical simulation (DNS) and large eddy simulation are two important strategies for the numerical investigation of turbulent flows. Within the constraints of present-day computing infrastructure, the direct simulation approach is adopted for full resolution of flow problems of modest complexity, e.g. to underpin theoretical and modelling studies (Pope 2000). Instead, the focus in large eddy simulation is on a computationally more accessible coarsened flow description, which is obtained by low-pass spatial filtering (Sagaut 2001; Geurts 2003). This allows an external control over the required spatial resolution. However, low-pass filtering gives rise to the well-known closure problem for the turbulent stress, which represents the dynamic effects of the filtered-out small-scale turbulence on the retained flow structures. Viewed entirely from the PDE level corresponding to the spatially filtered Navier–Stokes equations, the remaining task is to close the system of equations by modelling these small-scale dynamic effects in terms of the resolved flow. Considerable effort has been put into construction, testing and tuning of such so-called subgrid models over recent years (Meneveau & Katz 2000). In such testing procedures, one frequently compares predictions from large eddy simulations with filtered data from DNS and/or experimental observations.

The above sketch of direct and large eddy simulations is incomplete in at least one important respect, as it does not contain the unavoidable spatial and temporal discretization steps. In fact, since the numerical representation often adopts only marginal subgrid resolution, a significant alteration of the resolved scales' dynamics may be introduced in the computational model (Rogallo & Moin 1984; Ghosal 1996; Salvetti & Beux 1998; Geurts & Fröhlich 2002; Gullbrand & Chow 2003; Meyers *et al*. 2003).

It is the purpose of this paper to determine the effect of second- and fourth-order accurate finite-volume discretization of the convective and viscous fluxes on the outcome of large eddy simulation of homogeneous, isotropic, decaying turbulence based on the Smagorinsky subgrid model (Smagorinsky 1963). Changing the spatial discretization will have an effect on the value of the Smagorinsky parameter *C*_{S} at which a suitably constructed multiscale measure for the total simulation error is minimal, at given spatial resolution (Meyers *et al*. 2007). Here, the word ‘multiscale’ refers to the fact that the simulation error is quantified in terms of simultaneous contributions at different length scales. This effect will be stronger on coarse meshes (Vreman *et al*. 1996; Geurts & Fröhlich 2002), as adopted frequently in large eddy simulation of turbulent flow of realistic complexity. The partial cancellation of simulation errors due to discretization and subfilter modelling (Vreman *et al*. 1996) will be optimally balanced at different values of *C*_{S}, depending on the precise spatial discretization that is adopted. This paper is dedicated to quantifying this effect and interpreting it in terms of the implicit spatial filter induced by the finite-volume discretization that is adopted (Rogallo & Moin 1984).

The relative importance of the turbulent subgrid stresses compared with contributions due to discretization errors depends strongly on the subfilter resolution *r*=*Δ*/*h*, with *Δ* the width of the spatial filter and *h* the grid spacing (Ghosal 1996; Kravchenko & Moin 1997). If *r* is sufficiently large, the grid-independent large eddy solution consistent with the assumed value of *Δ* may be accurately approximated. However, large eddy simulation of applications with a realistic complexity is typically associated with only marginal resolution corresponding to *r*=1−2. In that case, the numerically induced effects are comparable to or even larger than the turbulent stresses for popular finite-volume discretization methods (Geurts & van der Bos 2005). Thereby, the computational large eddy closure typically contains an important contribution, which is sensitive to the adopted spatial discretization. In this paper, we study the interacting error dynamics in terms of a so-called error landscape that provides a concise visualization of the induced errors (Meyers *et al*. 2003). Each point in this landscape is characterized by the simulation error, relative to filtered DNS results, which occurs at a given Smagorinsky coefficient *C*_{S} and spatial resolution *N*^{3}. At given *N*, a coefficient may be inferred, at which the chosen measure for the error is minimal. Determining the locus of such *C*_{S} values yields the ‘optimal refinement strategy’ for the particular spatial discretization method (Meyers *et al*. 2007). Determining the sensitivity of this locus to the chosen finite-volume method is a main element of this paper.

The organization of this paper is as follows. In §2, we introduce the governing equations and consider the filter that is induced by the spatial discretization. Section 3 is devoted to an error-landscape analysis of large eddy simulation of homogeneous, isotropic, decaying turbulence in a ‘Smagorinsky fluid’. We show that an increased formal order of accuracy for the discretization of the convective fluxes implies higher values of , while higher order accuracy for the viscous fluxes results in a lowering of . Concluding remarks are collected in §4.

## 2. Spatial discretization and induced filter

In this section, we discuss the subgrid closure and the finite-volume discretization, focusing on the implicit spatial filter that is induced.

Filtering the Navier–Stokes equations in the large eddy framework requires a low-pass spatial filter *L* (Germano 1992). Here, a top-hat filter is adopted which, in one spatial dimension, couples a filtered velocity to an unfiltered velocity *u*(2.1)where *Δ* denotes the width of the filter. For incompressible fluids, the spatially filtered mass and momentum conservation equations can be written as follows:(2.2)Here, ∂_{t} (resp. ∂_{j}) denotes partial differentiation with respect to time *t* (resp. spatial coordinate *x*_{j}). Summation over repeated indices is implied. The component of the filtered velocity in the *x*_{j} direction is and is the filtered pressure. Finally, *Re* denotes the Reynolds number of the flow. The closure problem on the PDE level is expressed by the divergence of the turbulent stress tensor . To close the filtered equations, we approximate *τ*_{ij} in terms of the classical Smagorinsky model (Smagorinsky 1963), which amounts to(2.3)where *C*_{S} is the Smagorinsky coefficient and is the magnitude of the filtered rate of strain tensor .

The discretization of the governing equations introduces a separate source of approximation in a large eddy simulation. Therefore, the numerical method should be fully included in the analysis (Geurts & van der Bos 2005). We consider the filtered convective flux in one dimension, which contains a numerical flux term that is directly available in the computational model and a remainder grouped in the closure problem. The discrete derivative operator can be expressed as , where denotes the implicit filter associated with the specific discretization method that was used (Rogallo & Moin 1984; Geurts & van der Bos 2005). We find(2.4)in which we introduced the numerical turbulent stress tensor(2.5)where denotes a high-pass filter associated with . The difference between *Ξ* for a given discretization method and *τ* depends on the subgrid resolution *r*=*Δ*/*h*, which characterizes the strength of the numerical filtering. If the subgrid resolution *r* is sufficiently large, the numerical filter operator approaches the identity operator for all length scales relevant to *τ*. Hence, as follows from (2.5), this implies that *Ξ*→*τ*. In practical situations, the grid spacing is chosen such that *r* assumes quite modest values, and the numerical filter component in the full closure problem needs to be explicitly accounted for.

For central finite-difference schemes, denoted by , the induced filter may readily be inferred in the following form:(2.6)where denotes the induced filter and {*d*_{j}} are the differencing weights. For convenience, we restrict ourselves to uniform grids. The filter may be expressed in terms of a weighted average involving the top-hat filter. As an example, for the second-order accurate central discretization, we have(2.7)which corresponds to *n*=1 and *d*_{1}=1. Hence, we find(2.8)We may write the induced filter in terms of the local top-hat filter(2.9)for which *Δ*=*b*−*a*. Hence, we observe that the second-order accurate central differencing method induces a spatial top-hat filter with filter width equal to twice the grid-spacing *h*. The fourth-order accurate central differencing method corresponds to *n*=2 and may be written as(2.10)from which we infer that *d*_{1}=4/3 and *d*_{2}=−1/3. Hence, the induced filter corresponding to may be written as follows:(2.11)

The discretization of second-order derivatives can be expressed analogously, leading in one dimension to the well-known second-order method(2.12)and the fourth-order accurate method(2.13)The modification of the viscous flux due to numerical approximations can be understood in terms of the induced spatial filter as well, which is equivalent to a modified wave number analysis (Geurts 2003).

The effect of the induced numerical filter may be concisely illustrated through its operation on *u*=sin(*kx*). Specifically, we may write . For the second- and fourth-order accurate central finite differencing of a first derivative, we find(2.14)where *Γ*(*z*)=sin(*z*)/*z* is the Fourier transform of the kernel of the top-hat filter. In figure 1, we plot these characteristic functions. We note that the filter corresponding to a higher order spatial discretization is itself a higher order filter; with increasing order, the Fourier transform of the kernel stays closer to a value of 1 for a wider range of wave numbers (Geurts 1997; Vasilyev *et al*. 1998). Each discretization method gives rise to a particular damping of the amplitude of individual modes, which induces a specific dynamic contribution in an actual large eddy simulation.

The full three-dimensional finite-volume methods as adopted in this paper will be introduced next, closely following Meyers *et al*. (2007). For the discretization of the convective fluxes, the adopted second-order discretization corresponds to(2.15)wherewithwhere *f*_{i,j,k} is the convective-flux vector at node (*i*, *j*, *k*). Here, it is understood that the coordinate direction *x* corresponds to index *i*, and similarly *y*, *z* correspond to *j* and *k*. One may recognize a second-order accurate central finite-difference scheme in (2.15). This scheme is applied to the intermediate field *r*_{i,j,k} that is obtained by interpolation of the flux-vector field *f*_{i,j,k} in directions perpendicular to the direction with respect to which the derivative is evaluated. This additional averaging over *j* and *k* increases the robustness of the scheme (Vreman *et al*. 1996, 1997). The fourth-order convective discretization corresponds to(2.16)wherewith

For the discretization of the viscous terms and the subgrid-scale model, we use discretization schemes that are based on two consecutive applications of first derivatives. One may view this as the application of an ‘inner’ and an ‘outer’ first derivative. For the second-order discretization, the inner derivative is(2.17)wherewithShifting this discrete operator over half a grid spacing yields the outer derivative(2.18)wherewithThe sequential application of these inner and outer schemes yields the total, second-order accurate viscous flux treatment. For the fourth-order accurate discretization of the viscous and subgrid fluxes, a combination of two fourth-order accurate first derivatives is used (Yzerman 2000; Meyers *et al*. 2003). We adopt as inner method(2.19)wherewithCompletion of the second derivatives follows from application of the corresponding scheme, but now operating on staggered field data and returning an approximate derivative at the collocated grid locations.

Based on these discretization schemes, we define four different finite volume methods. We use the acronyms 2–2, 2–4, 4–2 and 4–4 for the different combinations, in which the first index refers to the order of accuracy of the method used for the convective fluxes and, similarly, the second index labels the viscous flux treatment. Throughout, the Smagorinsky models are treated identical to the viscous fluxes. In §3, we investigate the influence of the spatial discretization method on the total simulation error arising in actual large eddy simulations.

## 3. Error-landscape analysis of spatial discretization

In this section, we first review the error-landscape approach, which is an effective tool to characterize errors in large eddy simulations (Meyers *et al*. 2003). A central element of an error landscape is the so-called optimal refinement strategy which, at given resolution *N*, defines the value of the Smagorinsky coefficient that yields minimal error in a simulated quantity. We investigate the influence of the spatial discretization methods on the locus and correlate this with properties of the induced implicit spatial filter.

A database of large eddy simulations of homogeneous, isotropic, decaying turbulence at Taylor–Reynolds number *Re*_{λ}=100 was generated at a variety of resolutions and Smagorinsky coefficients, using grids with *N*^{3} grid cells where *N*≤96. For the direct simulation, up to 384^{3} grid cells were used and the level of convergence for various flow quantities was established (Kuczaj & Geurts 2007). To measure the error that arises in a simulation, we compare large eddy results with filtered DNS results. Specifically, we consider errors sensitive to different ranges of scales by introducing(3.1)in terms of a time integral of the difference between the LES and DNS spectra *E*_{LES} and , respectively, weighed with the *q*th power of the wave number *k*, for scales up to the grid-scale *k*_{c}. Through the integration over time *t* in (3.1), the error arising in each separate LES is represented by a single number. Variation of *q* allows error contributions to be emphasized in different wave number ranges. In fact, as *q*=−1, a connection with an integral scale can be made, *q*=0 suggests a connection to the resolved total kinetic energy, while *q*=2 is reminiscent of enstrophy. Each of the error measures *ϵ*_{q} defines a corresponding error landscape with its own optimal refinement . To arrive at a measure that takes errors at all these scales into account simultaneously, we use the following weighing:(3.2)The error landscape based on provides a multiscale optimal refinement strategy , to which we turn next.

As an illustration, the relative error *ϵ*_{0} as arises with the 2–2 method at various (*C*_{S}, *N*) combinations is shown in figure 2. A striking non-monotonous dependence of the simulation error on *C*_{S} is observed due to a partial cancellation of modelling and discretization errors (Vreman *et al*. 1996). Optimal results at a given resolution are obtained only when a proper balance between these sources of error is struck. The error landscapes obtained at different ‘scale parameter’ *q*, and at different discretization methods, are qualitatively the same and are hence not included explicitly.

In figure 3, we collected the optimal refinement strategies as obtained with the different spatial discretization methods. We observe considerable differences in the amount of subgrid fluxes that should optimally be added to a given central discretization, in order to achieve the desired balancing of discretization and modelling errors. Compared with the reference 2–2 method, when a fourth-order accurate convective discretization is adopted, the nonlinear generation of small scales becomes stronger since the implicit filter is less effective on a wider range of scales. Hence, a higher *C*_{S} is needed to optimally balance the error components (4–2 method). Conversely, if a fourth-order accurate discretization is adopted for the viscous fluxes, the rate of strain tensor is better approximated and a correspondingly lower level of subgrid dissipation is required to reach the optimum (2–4 method). For the 4–4 method, the optimal refinement strategy is quite similar to that of the 2–2 method.

## 4. Concluding remarks

We determined the sensitivity of the multiscale optimal refinement strategy, arising from the use of different second- and fourth-order accurate finite-volume discretization methods. At coarse spatial resolutions, the implicit filter that is induced by the spatial discretization has a significant effect on the deviation of the simulation results from DNS results. A higher order spatial discretization induces a filter that is less effective in attenuating the smaller resolved scales in a numerical solution. This has a direct influence on the locus of the corresponding optimal refinement strategy.

A higher order discretization method for the convective fluxes will imply a stronger presence of smaller scales on a given grid. However, these smaller scales are also more prone to numerical contamination at fixed resolution. To nevertheless achieve accurate simulation results, this implies that some of the contaminated numerical flow structures need to be toned down, requiring a slightly larger value of the Smagorinsky coefficient to achieve minimal total error. Hence, the optimal refinement strategy shifts upward when the convective fluxes are treated with a fourth-order accurate method, instead of a second-order one. The converse is observed when the treatment of the viscous fluxes is based on a method of higher formal order of accuracy. The numerical viscous flux is larger with such a higher order method, and a correspondingly smaller amount of subgrid dissipation is required to achieve optimal error balancing at coarse resolution.

## Acknowledgments

The author gratefully acknowledges stimulating discussions with Johan Meyers (KU Leuven). Simulations were performed at SARA and made possible through grant SG-213 of the Dutch National Computing Foundation (NCF).

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Applied large eddy simulation’.

- © 2009 The Royal Society