August 2012 

NASA/TM–2012-217593 

Parametric Study of Decay of Homogeneous 

Isotropic Turbulence Using Large Eddy 

Simulation 

R. C. Swanson, Christopher L. Rumsey, Robert Rubinstein, Ponnampalam Balakumar, and 

Thomas A. Zang 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

 TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

 CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

 TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

 Access the NASA STI program home page 

at http://www.sti.nasa.gov

 E-mail your question to help@sti.nasa.gov

 Fax your question to the NASA STI 

Information Desk at 443-757-5803 

 Phone the NASA STI Information Desk at 

443-757-5802 

 Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

August 2012 

NASA/TM–2012-217593 

Parametric Study of Decay of Homogeneous 

Isotropic Turbulence Using Large Eddy 

Simulation 

R. C. Swanson, Christopher L. Rumsey, Robert Rubinstein, Ponnampalam Balakumar, and 

Thomas A. Zang 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 

and Space Administration.
Abstract

Numerical simulations of decaying homogeneous isotropic turbulence are performed with both loworder and high-order spatial discretization schemes. The turbulent Mach and Reynolds numbers

for the simulations are 0.2 and 250, respectively. For the low-order schemes we use either secondorder central or third-order upwind biased differencing. For higher order approximations we apply weighted essentially non-oscillatory (WENO) schemes, both with linear and nonlinear weights.

There are two objectives in this preliminary effort to investigate possible schemes for large eddy

simulation (LES). One is to explore the capability of a widely used low-order computational fluid

dynamics (CFD) code to perform LES computations. The other is to determine the effect of higher

order accuracy (fifth, seventh, and ninth order) achieved with high-order upwind biased WENObased schemes. Turbulence statistics, such as kinetic energy, dissipation, and skewness, along with

the energy spectra from simulations of the decaying turbulence problem are used to assess and

compare the various numerical schemes. In addition, results from the best performing schemes are

compared with those from a spectral scheme. The effects of grid density, ranging from 32 cubed to

192 cubed, on the computations are also examined. The fifth-order WENO-based scheme is found

to be too dissipative, especially on the coarser grids. However, with the seventh-order and ninthorder WENO-based schemes we observe a significant improvement in accuracy relative to the lower

order LES schemes, as revealed by the computed peak in the energy dissipation and by the energy

spectrum.

Nomenclature

A = Jacobian matrix

ADM = Approximate Deconvolution Model

CFD = Computational Fluid Dynamics

DIA = Direct Interaction Approximation

DNS = Direct Numerical Simulation

DSM = Differential Stress Model

E = specific total energy

E(κ) = energy as a function of wave number in Fourier space

ENO = Essentially Non-Oscillatory (numerical scheme)

FCT = Flux-Corrected Transport (numerical scheme)

G = LES convolution kernel

Hi = sum of heat flux plus work done by stresses

L = length scale

LED = Local Extremum Diminishing

LES = Large Eddy Simulation

ILES = Implicit LES

M = Mach number

MILES = Monotone Implicit LES

MUSCL = Monotone Upstream-centered Schemes for Conservation Laws (numerical scheme)

PPM = Piecewise Parabolic Method (numerical scheme)

Pr = Prandtl number

Qi = subgrid scale heat flux

R = gas constant

Re = Reynolds number

Reλ = Taylor microscale Reynolds number

RK = Runge-Kutta

S = (2SijSij )

1/2

Sij = strain rate tensor

1
Sk = velocity derivative skewness factor

SGS = Sub-Grid Scale

SSM = Scale-Similarity Model

T = temperature

Tij = total stress tensor

TVD = Total Variation Diminishing (numerical scheme)

U = velocity

Ui = state vector

WENO = Weighted ENO (numerical scheme)

WENO-L = Linear version of WENO (numerical scheme)

cp = specific heat constant

cs = constant in Smagorinsky SGS model

cv = specific heat at constant volume

e = internal energy per unit volume

f = represents general variable; also represents numerical flux

k = filtered stress kinetic energy; also used as an index counter

p = pressure

t = time

t = time

t

∗ = time, nondimensionalized by Uref and Lref

ui = velocity component u, v, w

xi = Cartesian component x, y, z

Greek:

∆ = LES filter width; also represents a spatial difference

∆ = LES grid filter

βk = smoothness measure of WENO scheme

δij = Kronecker delta

 = turbulent dissipation rate

εw = parameter to prevent division by zero in WENO weights

γ = ratio of specific heats

η = Kolmogorov length scale

κ = wave number in Fourier space

κ

0 = coefficient of thermal conductivity

κ = variable in the MUSCL scheme

λ = Taylor microscale

μ = dynamic viscosity

ν = kinematic viscosity

ρ = density

σij = laminar part of the stress tensor

τ = Kolmogorov time scale

τij = subgrid scale stress

χ = ratio of compressible kinetic energy to total kinetic energy

ωk = weights of WENO scheme

φ = represents general space-time variable

Subscripts, Superscripts, and Special Symbols:

overbar = resolved part; also spatially-filtered variable

tilde = Favre-averaged variable

angled brackets = average over all cells in the domain

C = compressible

I = incompressible

2
L = left

R = right

ref = reference

rms = root-mean square

t = turbulent

x, y, z = in x, y, z directions

0 = initial

∇ = mathematical “del” differential operator

1 Introduction

1.1 General Description of DNS and LES

Advancements in computer power, both processing speed and memory, in recent years and the

need for improved representation of the effects of turbulence in practical engineering computations

has encouraged increased interest in direct numerical simulation (DNS) and large eddy simulation

(LES). In DNS the full Navier-Stokes equations are solved, and all the relevant scales, down to

the Kolmogorov scale, of a physical problem are resolved. Unfortunately, even with the current

computer capability, the computational requirements for DNS are too demanding for high Reynolds

number (Re) problems since the number of grid points required is proportional to Re9/4. Furthermore, as indicated by Piomelli [1], the computational cost scales like Re3

. For these reasons DNS

has generally been applied to flow problems with relatively low Re and simple geometries. Instead

of resolving all the relevant scales of the turbulent motion, LES relies upon the assumption that the

smaller scales of the turbulent motion have a more universal behavior that allows their effect (due

to energy dissipation) on the larger scales to be represented by appropriate modeling. By modeling

the physics of the smallest scales there is a significant reduction in the computational requirements

relative to DNS, allowing the possibility to compute high Re turbulent flows over complex geometric configurations. There are two general approaches to LES. One involves the filtering of the

Navier-Stokes equations in order to eliminate scales smaller than a prescribed value. The filtering

introduces additional stresses and heat fluxes that are determined by the unresolved small scale (subgrid scale) eddies. This approach relies on direct (explicit) modeling of the additional terms, and it

is often called conventional or explicit LES. The second LES formulation depends on the numerical

dissipation associated with the spatial discretization of the flow equations to emulate the physics of

the smallest scales. This approach is called integrated or implicit LES (ILES).

1.2 Conventional LES Methods

In conventional LES the filtering to separate the large energy carrying eddies from the unresolved

small scales is designed to remove the scales equal to or smaller than some filter width. This filter

width is frequently represented by ∆, and it defines the characteristic resolution length. For the LES

∆  η, with η = O(Re−3/4) denoting the Kolmogorov length scale, whereas for DNS there is no

“filter” so it is equivalent to LES with ∆ ≈ η. The ultimate success in the LES computation depends

upon how well the effects of the filtered scales on the larger resolved scales are represented. The explicit subgrid-scale (SGS) models of conventional LES have two general classifications: functional

models and structural models. There are also mixtures of these two types (e.g., functional models

with insufficient dissipation of the small scales also require some type of additional dissipation).

Functional models represent the subgrid terms (i.e., stress tensor and heat flux vector) indirectly

by modeling the effect of the dissipative action of the small scales on the larger scales, whereas

structural models represent the subgrid terms directly by approximating the unfiltered flow variables.

One example of a functional model of the subgrid scale stress tensor is the Smagorinsky model [2].

With this algebraic eddy viscosity model the SGS stress tensor is a function of the mean strain rate

3
tensor and the square of the filter width, which is usually taken as the mesh cell size. There is also a

constant that must be specified. An improvement of this methodology is the dynamic Smagorinsky

model [3], which replaces the square of the constant by a coefficient that is dynamically adjusted

according to information available at the smallest resolved scales. In general, eddy-viscosity models

can reproduce SGS dissipation but not the subgrid forces in the filtered momentum equations.

In the case of the structural models the objective is to reconstruct some of the SGS information

directly by generating a new solution field from the filtered one, resulting in an improved approximation of the unfiltered field. Structural models are generally superior to functional models because

they do not have the isotropic assumption that most functional models have. Some examples of

structural models are the scale-similarity model (SSM) and mixed model of Bardina et al. [4], the

differential stress equation model (DSM) of Deardorff [5], and the approximate deconvolution model

(ADM) of Stolz and Adams [6, 7]. In essence, the scale similarity hypothesis [4] assumes that the

effect of the smallest resolved scales on the unresolved scales is similar to the effect of the largest

resolved scales on the smallest resolved scales. Based on this similarity assumption, the SSM approximates the SGS stress tensor using the resolved velocity field and a second application of the

spatial filter. Due to the filtering process this model allows for backscatter (i.e., energy transfer

from the subgrid scales to the resolved scales). The DSM [5] uses a transport equation to model

the subgrid stress tensor. For the ADM [6, 7] the unfiltered velocity is estimated by an approximate

deconvolution operation, which uses an approximation of the inverse filter kernel. Like the SSM, the

ADM does not have adequate dissipation of energy. In order to represent the energy transfer from

the resolved to unresolved scales, an additional relaxation is included in the discrete equations. This

is equivalent to an explicit filtering of the resolved velocity field at each time step in the evolution

of the solution. As pointed out by Grinstein and Fureby [8], a principal difficulty with explicit SGS

models is the need to separate the effects of explicit filtering and SGS reconstruction models from

those due to discretization.

1.3 ILES Methods

In general, one could classify any LES method that uses numerical dissipation rather than an explicit

model to represent the effect of the subgrid scales as an ILES method. According to Sagaut [9] any

ILES approach is implicitly based on the following hypothesis: The action of subgrid scales on the

resolved scales is equivalent to a strictly dissipative action. The numerical dissipative terms are

produced when approximating the convective terms in the fluid dynamic equations. These terms

are usually either introduced implicitly through an upwind scheme framework or explicitly when

using a central difference scheme. The numerical dissipation plays several roles. One role is to

ensure numerical stability. In the linear sense, assuming well-posedness and consistency, stability

is a necessary and sufficient condition for convergence (Lax equivalence theorem, see Richtmyer

and Morton [10]). It is not a sufficient condition to guarantee physically realizable solutions. To

ensure a unique and physically meaningful solution the dissipation must also allow satisfaction of

an entropy condition. Another role for dissipation is to prevent oscillations in the solution when

extrema occur due to physical phenomena such as shock waves and contact discontinuities. There is,

in general, a reduction in the approximation order (e.g., second-order schemes become first order) for

the convective terms at discontinuities in order to enable discontinuity capturing without oscillations.

Frequently, the distinction between upwind and central difference methods depends primarily on

the nature of the limiting (or switching) that is imposed to change the order of approximation at

discontinuities.

A subclass of ILES algorithms is the monotone integrated LES (MILES) schemes. To understand

what methods can be placed in this category it is appropriate to provide some historical background

and clarification. The MILES schemes are defined to be at least second-order accurate in smooth

regions of the flow field. The classification of MILES was created in 1990 by Boris [11] in reference

to the flux corrected transport (FCT) scheme introduced by Boris [12] in 1971 and generalized by

Boris and Book [13] – [15]. It stems from the fact that the FCT scheme has a predictor-corrector

4
structure with the first approximate solution determined with a first-order monotone scheme. This

initial approximation is then corrected by removing the large dissipative error with a limited antidiffusion. Through the limiting process monotonicity and positivity conditions are imposed. The

monotonicity constraint, which is the underlying FCT principle given in [12] and [13], ensures

that no new extrema are created and that there is no growth of existing extrema. The positivity

condition prevents positive definite physical quantities (e.g., density) from becoming negative. The

resulting algorithm is at least second-order accurate in smooth regions, and it satisfies monotonicitypreserving and positivity-preserving properties. In addition, the algorithm is designed to satisfy the

physical properties of global conservation, causality, and locality. Causality simply means that a

fluid passing between two locations in a discrete domain must pass through all the cells in the flow

path between those locations. Since the FCT scheme is conservative, causality must be satisfied.

Locality requires the approximations to derivatives to be compact.

Boris classifies the FCT scheme as a monotone fluid dynamics algorithm. In [16] Boris lists

other examples of monotone algorithms: monotone upstream-centered scheme for conservation

laws (MUSCL) of van Leer [17] – [21]; piecewise parabolic method (PPM) of Woodward and

Colella [22, 23]; total variation diminishing (TVD) schemes of Harten [24, 25], (TVD) schemes

of Harten [24, 25], Sweby [26], Yee et al. [27]; and second-order Godunov methods of Colella [28]

and Colella and Glaz [29]. Monotone schemes are monotonicity preserving and guarantee a monotone solution if the scheme converges. In addition, convergent monotone methods generate solutions

that automatically satisfy an entropy condition. This means that if the scheme is consistent and conservative, then it not only generates a weak solution (Lax-Wendroff theorem [30]) but also a unique

solution. Consequently, both smooth and discontinuous flow solutions are necessarily the physically

correct solutions. We should also point out that these monotonicity-preserving schemes impose the

total variation bounded (TVB) condition, which is a weaker condition than TVD. This is because

TVD schemes are at most first-order accurate (see theorems of Godunov [31], Harten [32]).

There are several approaches to building schemes that preserve monotonicity. One is to start with

a low-order monotonic scheme and introduce anti-diffusion that is limited to ensure local monotonicity (e.g., FCT scheme). Another approach is to use a monotone interpolation method (e.g., MUSCL)

to determine the primitive variables for reconstruction of the fluxes at the cell interfaces. Limiting is

applied to the slopes used in the interpolation to impose monotonicity. A third way is to blend loworder and high-order diffusive terms and choose appropriate coefficients to ensure that the scheme

is positive; and thus, the scheme is local extremum diminishing (LED). The concept of LED was

introduced by Jameson [33, 34] to provide a theoretical foundation for extending the non-oscillatory

design principle (i.e., principle of nonincreasing maxima and nondecreasing minima) to multiple dimensions, since the TVD property is a one-dimensional concept. In all of these methods the limiter

must satisfy certain properties. In general, these are natural averaging properties and the property

to ensure satisfaction of local monotonicity (e.g., LED). These schemes are also designed to preserve positivity. In addition, by relaxing the LED property (i.e., essentially LED), these schemes

can be made higher order even at smooth extrema. Although most of the theory for non-oscillatory

algorithms has been developed using scalar equations and constant coefficient one-dimensional systems, these nonlinear high-resolution schemes have proved to be effective for the multidimensional

nonlinear systems of fluid dynamics.

So far we have discussed primarily second-order non-oscillatory methods, with some limited

reference to the third-order scheme called PPM, which is based on a piecewise quadratic representation of the solution. In general, higher order methods (order greater than two) rely upon techniques

such as filtering and/or nonlinear coefficients (e.g., essentially non-oscillatory schemes) to prevent

oscillations at extrema and to ensure stability (i.e., to prevent unbounded energy growth due to highfrequency error components). For example, Yee and Sj ̈ogreen [35] developed a class of high-order

central difference schemes that apply a multistep filter method that includes high-order linear and

nonlinear filters. The linear filter consists of a product of a sensor and either a spectral-like filter or

a high-order centered dissipative operator. The nonlinear filter consists of a product of a different

sensor and either an artificial compression method indicator or a wavelet sensor and the dissipative

5
part of a high-resolution shock-capturing algorithm. Another high-order approach for preventing

spurious oscillations at extrema includes the essentially non-oscillatory (ENO) and weighted essentially non-oscillatory (WENO) schemes (see Shu [36]). The ENO schemes are Godunov-type

schemes that achieve high-order spatial accuracy in smooth regions by a piecewise polynomial approximation of the solution based on cell-averaged values. A nonlinear adaptive procedure is used

to automatically choose the smoothest local approximation stencil, and thus, avoid the Gibbs phenomenon at discontinuities (i.e., generation of O(1) spurious oscillations). Both the reconstruction

and the resulting scheme are essentially non-oscillatory (i.e., almost TVD). The ENO schemes are

essentially non-oscillatory because they allow small spurious oscillations that have size on the order

of h

r

, where h is the mesh spacing and r is the approximation order of the scheme. Instead of using

the optimal stencil, the WENO schemes use a convex combination of the candidate stencils with

nonlinear adaptive coefficients (weights). In practice, both the ENO and WENO schemes generally

produce solutions without oscillations.

There are also hybrid high-order algorithms that use a WENO algorithm in the neighborhood of

discontinuities and non-dissipative algorithms such as a tuned center difference method or spectral

method (exponential convergence method) for smooth regions (e.g., Hill and Pullin [37], Costa and

Don et al. [38]). Another type of high-order scheme that is frequently used is the compact scheme,

which is based upon a Pade differencing operator. Examples of such schemes are the centered

sixth-order and fourth-order schemes using five and three points, respectively, employed by Visbal

and Rizzeta [39] and Visbal and Gaitonde [40]. They use a high-order low-pass spatial filtering

technique to provide dissipation at the high modified wavenumbers (preventing energy build up at

the high frequencies) when there are significant dispersion errors due to the spatial discretization.

At discontinuities they use a Roe scheme with the MUSCL approach for reconstruction, which

produces a third-order upwind-biased scheme in the case of uniform grids. Various types of detection

mechanisms have been considered with these compact schemes, and they include: (1) Jameson-type

pressure switch, (2) the artificial compression method, (3) multiresolution wavelets, and (4) WENOtype smoothness criteria.

At this point we have simply identified different types of ILES schemes. Now, we need to be

more specific about which ILES schemes are appropriate for indirectly introducing a subgrid scale

model. Grinstein and Fureby [8] have demonstrated that the leading order truncation error terms for

a certain class of ILES methods (i.e., methods having a flux function determined by blending of loworder and high-order fluxes with an appropriate limiter function) provide implicit SGS models that

are similar in form to conventional (explicit) mixed SGS models. While a general theory for the ILES

approach is yet to be established, the results obtained for fundamental turbulence problems such as

isotropic turbulence decay and Taylor-Green vortex problems provide some support to the validity

of the approach. Nevertheless, general validity can only be established by directly connecting the

ILES methodology to the physical behavior of small scale turbulence.

For both conventional LES (with explicit SGS models) and ILES (for which the SGS model is

“implicitly” included through the numerical scheme), the magnitude of the numerical errors relative

to the magnitude of the subgrid terms is critical in determining the accuracy of LES computations.

There are two fundamental types of numerical errors: discretization (truncation) errors and aliasing errors. The truncation errors are determined by the order of approximation of the terms in the

governing flow equations. Aliasing errors are caused by the interaction of nonlinear terms, resulting in high-frequency errors reappearing as low-frequency errors. For central differencing and pure

upwind differencing the leading truncation error terms are dispersive and dissipative, respectively.

Upwind biased schemes have both types of errors. In order to accurately represent the turbulent kinetic energy spectrum the difference approximation of convective terms must have low dissipation.

At the same time, the dissipation should be sufficiently large to suppress oscillations and ensure

numerical stability. In addition, there is a need to have appropriate dissipation for shock capturing.

There is also a need to have low dispersion errors. For any scheme being used for LES, these numerical errors must be monitored and controlled to ensure that they are sufficiently small so that the

dominant errors in the simulation are due to SGS modeling. Only then will it be possible to effec6
tively improve the models and obtain a simulation capability for a broad range of flow conditions.

Since many existing schemes use upwind biased spatial differencing for shock capturing, there is

clearly a need to investigate and validate these schemes for LES computations.

1.4 Related Work and Current Contribution

Garnier et al. [41] investigated six shock capturing schemes for LES computations. These methods

included the following: Jameson scheme with scalar dissipation, Roe scheme with the MUSCL extrapolation procedure and two minmod limiter variations, and three different types of ENO schemes.

Their spatial accuracy range was between second and fifth order, and either second-order or thirdorder TVD Runge-Kutta (RK) schemes were used for time integration. Each scheme was applied to

decaying isotropic turbulence problems ignoring viscous effects, and Mach numbers were between

0.2 and 1.0. Simulations were performed for mesh resolutions of 643and 1283. The effects of both

implicit and explicit SGS models were discussed. Garnier et al. [41] considered a scheme suitable

for LES if it satisfied one of the following two conditions: (1) numerical dissipation is much smaller

than the physical SGS dissipation, or (2) numerical dissipation mimics that of a SGS model. They

concluded that none of the schemes tested satisfied either of these two conditions.

Another study concerning the behavior of high-resolution and high-order methods for implicit

LES was performed by Thornber et al. [42]. Their investigation, at low Mach number, was also done

for homogeneous decaying turbulence. Computing on mesh densities between 323and 2563, they

considered seven different schemes with spatial accuracy ranging from second to ninth order. In

all simulations the Euler equations were solved with finite-volume methods. For higher order spatial accuracy, cell interface variables were determined by either the MUSCL or WENO approach.

Third-order RK schemes were used to advance the solution in time. With grids exceeding 323 or

with numerical methods having spatial accuracy higher than third order, their results indicated that

the growth of large scales and the dissipation of kinetic energy could be captured. Based on energy spectra, Thornber et al. also concluded that all the upwind-biased schemes evaluated were too

dissipative at the high wavenumbers, resulting in an early onset of decay in the inertial range. We

should point out that with this approach the question does arise concerning an appropriate representation of an inertial subrange in the absence of physical viscous effects (i.e., the physical effects

when Re = ∞ are clearly not the same as for high Re). Thornber et al. also observed that due to

differences in the effective cut-off wavenumber, methods with fifth-order or higher spatial accuracy

performed better on a 1283 grid than those with second-order accuracy on a 2563 grid.

The primary objective of this report is to compare the behavior of representative second-order

finite-volume schemes and high-order finite-difference WENO-based schemes when applied to the

problem of isotropic decay of turbulence. The second-order schemes have either central or thirdorder upwind biased spatial differencing, and for each scheme second-order temporal differencing is used. Both explicit and implicit LES computations are performed with these methods. For

the WENO schemes we investigate fifth-order, seventh-order, and ninth-order spatial accuracy for

ILES. A RK scheme, which is third-order TVD, is used for time advancement. There is a recognition (Martin et al. [43], Taylor and Martin [44]) that the standard fifth-order WENO scheme with

Lax-Friedrichs flux splitting is too dissipative at the higher wavenumbers, as clearly revealed in the

energy spectrum. Here, the effect of increasing the order of approximation is discussed. In the simulations, the filtered Navier-Stokes equations are solved. Since the Kolmogorov scale is significantly

smaller than the mesh size, Garnier et al. [41] and Thornber et al. [42] indicated either indirectly

or directly that the viscous effects are negligible. However, in any practical simulation the viscous

terms and their effects are present, regardless of the spatial resolution. This certainly suggests that

including these terms is necessary to obtain an effective measure of the wave cut-off number in the

energy spectrum. Moreover, as discussed by Ishihara [45], the turbulence statistics have significant

Re dependencies, even at high values of Re.

In the initial part of the paper the governing equations are presented. Then the numerical methods

are described. This is followed by the results section, which begins by describing the formulation of

7
the initial solution and providing definitions of the turbulence statistics. In the next part the results

from LES simulations are presented. The turbulence results include kinetic energy, dissipation of

energy, and skewness. Vorticity and Q-invariant contours are also shown. The initial peak of the

energy dissipation during the energy redistribution process clearly reveals the improved accuracy

achieved with the higher order WENO schemes. The computed energy spectrum is used to determine

the dissipative behavior of the numerical methods over the frequency range. Results with the various

numerical schemes are compared with each other. In addition, there are some comparisons with

results from a spectral scheme.

2 Governing Equations

2.1 Direct Numerical Simulation (DNS)

Written in strong conservation form, the compressible form of the Navier-Stokes equations is:

∂ρ

∂t +

∂ρuk

∂xk

= 0 (1)

∂ρui

∂t +

∂ρuiuk

∂xk

= −

∂p

∂xi

+

∂

∂xk



μ



∂ui

∂xk

+

∂uk

∂xi

−

2

3

δki

∂uj

∂xj

 (2)

∂ρE

∂t +

∂[(ρE + p)uk]

∂xk

=

∂

∂xk



κ

0 ∂T

∂xk

+



μ



∂ui

∂xk

+

∂uk

∂xi

−

2

3

δki

∂uj

∂xj

 ui



(3)

where the coefficient of thermal conductivity κ

0

is taken to be κ

0 = cpμ/Pr, and the variable E

represents the specific total energy E = e + (u

2 + v2 + w2

)/2. For a perfect gas, e = cvT,

cv = R/(γ − 1), γ = cp/cv, and the equation of state is p = ρRT, or:

p = (γ − 1) ρE −

1

2

ρ



u

2 + v2 + w2





(4)

By definition, DNS is a numerical computation of the Navier-Stokes equations that resolves all

relevant spatial and temporal scales in a flow field. This means that the grid needs to be fine enough

to resolve features on the order of the Kolmogorov dissipation length scale:

η =



ν

3

/ε 1

4

(5)

where ν is the kinematic viscosity and ε is the dissipation rate. Similarly, the time step needs to be

fine enough to resolve temporal dynamics on the order of the Kolmogorov time scale:

τ = (ν/ε)

1

2

(6)

It is impossible today to achieve these resolutions at reasonably high Reynolds numbers. Successful

low Reynolds number simulations require very low-dissipation schemes to resolve the smallestscale features with good fidelity. The spatial order of accuracy of the scheme also influences grid

requirements. Additional dissipation is required to stabilize some numerical schemes; for example,

compact schemes are often stabilized through the use of additional filtering [46]. Note that an underresolved DNS simulation can exhibit time-dependent features that appear to the eye like an adequate

DNS simulation, but if the smallest space and time scales are not resolved, some flow physics are

being omitted and its validity remains unclear.

8
2.2 Large Eddy Simulation (LES)

In LES, the Navier-Stokes equations are filtered via a low-pass filter defined as a convolution product:

φ(~x, t) = Z +∞

−∞

Z +∞

−∞

φ(

~ξ, t0

)G(~x − ~ξ, t − t

0

)dt0d

3~ξ (7)

where φ(~x, t) is the resolved part of a space-time variable φ(~x, t). For the compressible NavierStokes equations, the equations contain both ordinary and Favre-filtered variables (see, e.g., Knight

et al. [47]). The Favre filter is defined by:

 ̃f =

ρf

ρ

(8)

where the overbar denotes the ordinary spatially filtered variable, and the tilde denotes the Favrefiltered variable. After filtering, the compressible governing equations (Eqs. (1) - (3)) become:

∂ρ

∂t +

∂ρu ̃k

∂xk

= 0 (9)

∂ρu ̃i

∂t +

∂ρu ̃iu ̃k

∂xk

= −

∂p

∂xi

−

∂Tik

∂xk

(10)

∂ρE ̃

∂t +

∂[(ρE ̃ + p) ̃uk]

∂xk

=

∂Hk

∂xk

(11)

where the total stress tensor is Tik ≡ ρτik − σik, and the sum of the heat flux plus work done

by the stresses is Hk ≡ Qk + κ

0∂T /∂x  ̃

k + Tiku ̃i with κ

0

typically expressed as cpμ/Pr, and σik

approximated as:

σik = μ



∂u ̃i

∂xk

+

∂u ̃k

∂xi

−

2

3

δik

∂u ̃j

∂xj



(12)

As a result of the filtering, the effect of the filtered stress kinetic energy (k ≡ τii/2) should be

included in the definition of the total energy (Knight et al. [47]), that is, ρE ̃ = ρe ̃ + ρ( ̃u

2 +  ̃v2 +

w ̃

2

)/2 + ρk. Then, the equation of state becomes:

p = (γ − 1) ρE ̃ −

1

2

ρ



u ̃

2 +  ̃v2 +  ̃w2



− ρk(13)

However, for lower speed flows, the effect of k is usually negligible, and many researchers ignore it

in the definition of energy and equation of state.

As a result of the filtering, two unknown quantities emerge: the subgrid scale stress τik and the

heat flux Qk:

τik = ugiuk − u ̃iu ̃k (14)

Qk = −cpρ



ugkT − u ̃kT ̃



(15)

The quantities in Eqs. (14) and (15) must be modeled. As discussed in Sagaut [9], there are two

basic methods: explicit modeling and implicit modeling. In explicit modeling, a subgrid model

is explicitly introduced. In implicit modeling (ILES), no additional eddy viscosity is introduced

into the governing equations, but the numerical method is selected such that the numerical error

fulfills desired properties and effectively acts like a subgrid model. With no explicit model present,

ILES appears to be functionally similar (from an implementation point of view) to under-resolved

9
DNS; neither method resolves the finest scales, and both require enough inherent dissipation in the

numerical scheme to prevent the non-physical build-up of energy at the smallest resolved scales.

The commonly used Smagorinsky subgrid scale model [2] is an explicit SGS model, which

utilizes the the Boussinesq assumption for the SGS stress tensor:

τik = −2νt



Sik −

1

3

δikSjj+

2

3

δikk (16)

where the νt in the Smagorinsky model is given by:

νt =



cs∆

2

|S| (17)

with |S| =

q

2SijSij , Sij = (∂u ̃i/∂xj +∂u ̃j/∂xi)/2, and ∆ is typically defined by some measure

of the local grid spacing, such as ∆ = (∆x∆y∆z)

(1/3). The cs term is a constant for the original

Smagorinsky model, commonly set to approximately 0.1.

Note that in most subgrid scale models, neither the formal filter function nor the filter width are

explicitly defined. Instead, the models have a built-in filter, related to the local grid spacing: for the

Smagorinsky model it is cs∆. This built-in filter controls the size of the smallest locally resolved

flow structures. It can be difficult to choose an optimum built-in filter width, which has been noted to

be flow dependent. This difficulty has been one reason for the development and subsequent success

of the widely used dynamic model of Germano et al. [3]. This model dynamically computes a

variable cs term, rather than setting it to a constant value.

For compressible flows, an SGS viscosity model for the subgrid-scale heat flux is (see Urbin and

Knight [48]):

Qk =

ρcpνt

Prt

∂T ̃

∂xk

(18)

3 Numerical Methods

For the current problem of decay of isotropic turbulence, two compressible CFD codes were employed, the NASA production computer code CFL3D [49] and an in-house NASA research code that

employs the weighted essentially non-oscillatory (WENO) scheme [36]. Some of the characteristics

of these numerical methods will be described next.

3.1 Upwinding and Central Differencing in CFL3D

The CFL3D code is a structured-grid upwind multi-zone CFD code that solves the generalized thinlayer or full Navier-Stokes equations [49]. The code can employ local time-step scaling, grid sequencing and multigrid to accelerate convergence to steady state. A second-order time-accurate

mode (using standard backward differentiation formula (BDF)) is used in this study with pseudotime stepping plus multigrid, and an implicit approximate factorization method to advance the solution to the next physical time step. For the current study, a 3-level W-cycle multigrid strategy was

employed, along with 30 sub-iterations (pseudo-time steps) per time step. This number of subiterations was enough to drive the L2-norm of the subiteration density residual down approximately 14

orders of magnitude, to near machine zero for the current problem.

The CFL3D code employs a cell-centered finite-volume method. It uses the so-called “kappa

scheme” in a MUSCL approach [50] for its spatial differencing for the convective and pressure

terms, and second-order central differencing for the viscous terms. Roe’s flux difference-splitting

(FDS) method [51] is used to obtain fluxes at the cell faces. In the kappa scheme, the reconstruction

10
at cell interfaces is accomplished via:

UL = Ui +

1

4

[(1 − κ)(Ui − Ui−1) + (1 + κ)(Ui+1 − Ui)] (19)

UR = Ui+1 −

1

4

[(1 + κ)(Ui+1 − Ui) + (1 − κ)(Ui+2 − Ui+1)] (20)

where Ui

is the state vector at cell center i, and the subscripts L and R refer to the left and right

states for a local approximate Riemann solver. The “standard” scheme for CFL3D in RANS applications uses κ = 1/3, which is a third-order upwind-biased stencil on a uniformly spaced Cartesian

grid [52]. For the LES applications in the current paper, upwinding was believed to be too dissipative, so for most of the applications the central-differencing option κ = 1 was used. When κ = 1,

UL = UR =

1

2

(Ui + Ui+1) (21)

The Roe-type numerical flux is written as

fi+1/2 =

1

2

(fL + fR) −

1

2

|A|

i+1/2

(UR − UL) (22)

Since UL = UR, the dissipative term in Eq. (22) vanishes, and the flux balance for a given cell

yields a standard central difference approximation with second-order spatial accuracy. In general,

central differencing with no additional dissipation can yield numerical problems such as odd-even

point decoupling. However, for the particular case of isotropic turbulence decay considered here,

the solutions on the grid sizes considered were stabilized sufficiently by the viscosity of the fluid

(although not shown, computations attempted on a very fine grid with 2563cells did yield odd-even

point decoupling behavior).

The standard Smagorinsky model has been coded into CFL3D, with the constant coefficient cs

specified by the user. When cs is set to 0, the method becomes an ILES type, with no additional

subgrid eddy viscosity model added.

3.2 Weighted Essentially Non-Oscillatory (WENO) Scheme

In this section, for completeness, we present the elements of the weighted essentially non-oscillatory

(WENO) scheme and the implementation of this scheme in the present work. To simplify the description of the scheme, we consider a one-dimensional scalar conservation equation

∂u

∂t +

∂f(u)

∂x =

∂f v(ux)

∂x (23)

where u = u(x, t) is a conserved quantity, f(u) and f

v

(ux) are inviscid and viscous flux functions

with ux denoting the derivative ∂u/∂x, and the independent variables x and t represent space and

time, respectively. By taking this approach there is no compromise of the generality of the discussion since the same spatial approximation techniques can be applied in each coordinate direction of

a multidimensional fluid dynamics problem. We obtain a semi-discrete form of Eq. (23) by considering a one-dimensional domain with uniform spacing ∆x. The unknowns are located at the points

xi = i∆x, i = 0, . . . , N, and cell boundaries are at xi+1/2 = xi + ∆x/2. As seen in Fig. 1, the

grid points are also cell centers. The semi-discrete form of Eq. (23) is given by

dui(t)

dt = −

∂f

∂x









x=xi

+

∂f v

∂x









x=xi

(24)

where ui(t) is a numerical approximation of the point value u(xi, t). To construct a conservative

difference scheme for this system of ordinary differential equations, the fluxes must be evaluated

11
at the interfaces of the cells to the desired order of accuracy. Since the fifth-order WENO scheme

is frequently used in applications, we present the details of the discretization for the derivative of

f(u) using a fifth-order finite-difference form of this scheme, following the work of Shu [36], Jiang

and Shu [53], Henrick et al. [54], and Borges et al. [55]. After discussing the approximation of

convective term, we then consider different methods for discretizing the viscous (i.e., diffusion)

term, giving the advantages and disadvantages of each method. Next, we provide a brief discussion

of the extension to even higher order based on the paper by Balsara and Shu [56]. In the last part of

this subsection we describe the time integration method.

3.2.1 Approximating the Inviscid Term

This section contain three parts. In the first part we consider the linear form of the discretization for

the derivative of the convective (inviscid) flux function. This form is dominant in smooth regions

of the domain, and it provides a natural foundation for the development of the nonlinear WENO

scheme. The final part of this section discusses the flux splitting used to account for different signs

of the flux Jacobian.

3.2.1a Linear Form of 5th Order WENO Scheme

Consider the inviscid flux function. One can satisfy the conservative property by defining a numerical flux function h(x) implicitly as

f(x) = 1

∆x

Z x+∆x/2

x−∆x/2

h(ξ)dξ (25)

so that the spatial derivative of f(x) is exact; and thus,

∂f

∂x









x=xi

=

1

∆x



hi+1/2 − hi−1/2



(26)

where hi±1/2 = h(xi±1/2). The numerical flux function hi±1/2 can be approximated by a flux

function ˆfi±1/2 that is defined by a polynomial of sufficient degree to obtain the desired order of

accuracy. Then

∂f

∂x









x=xi

≈

1

∆x

(

ˆfi+1/2 − ˆfi−1/2) (27)

and

ˆfi+1/2 = h(xi+1/2) + O(∆x

r

) (28)

where r is the desired order of accuracy. The numerical flux function ˆfi+1/2 depends on the discrete

values of the physical flux f at the points which determine the stencil of the approximation.

For the WENO scheme a fifth-order accurate flux ˆfi+1/2 is computed as a convex combination

(i.e., weighted average) of possible third-order approximations. In Fig. 2 the three candidate stencils

are shown. The weighted average of the third-order stencils can be expressed as

hi+1/2 ≈ ˆfi+1/2 =

Xr−1

k=0

ωk

ˆf

k

i+1/2

(29)

where r = 3 and ωk is the weight corresponding to stencil Sk. To obtain the third-order approximations for each ˆf

k

(x), we represent each approximate numerical flux function with a second-degree

polynomial, and thus

h(x) ≈ ˆf(x) = a0 + a1x + a2x

2

(30)

12
with the undetermined coefficients ak, k = 0, 1, 2. Substituting ˆf(x) for h(x) in Eq. (25) and

integrating we obtain

ˆf(x) = a0 + a1x + a2



x

2 +

∆x

2

12 

(31)

The coefficients ak of the polynomial are determined from values of f(x) at nodes of the grid. Consider the grid cell with boundaries at xi−1/2 and xi+1/2. To obtain the ak so that the flux of Eq. (29)

produces a fifth-order upstream-centered scheme at the interface xi+1/2, the values of f(x) are

required at the following three sets of nodes: {xi−2, xi−1, xi}, {xi−1, xi

, xi+1}, {xi, xi+1, xi+2}.

Evaluating the polynomial of Eq. (31) for each set of nodes gives three algebraic equations that must

be solved to obtain the coefficients of the interpolating polynomial. For the evaluations the value of

xiis zero. The numerical flux functions for the three sets of points are

ˆf

0

(x) = −fi−2 + 2fi−1 + 23fi

24

+



fi−2 − 4fi−1 + 3fi

2∆x



x +



fi−2 − 2fi−1 + fi

2∆x

2



x

2

ˆf

1

(x) = −fi−1 + 26fi + fi+1

24

+



fi+1 − fi−1

2∆x



x +



fi−1 − 2fi + fi+1

2∆x

2



x

2

(32)

ˆf

2

(x) = (23fi + 2fi+1 − fi+2)

24

+



−3fi + 4fi+1 − fi+2

2∆x



x +



fi − 2fi+1 + fi+2

2∆x

2



x

2

At i + 1/2 (xi+1/2 = ∆x/2) the fluxes of Eq. (32) become

ˆf

0

i+1/2 =

1

6

(2fi−2 − 7fi−1 + 11fi)

ˆf

1

i+1/2 =

1

6

(−fi−1 + 5fi + 2fi+1) (33)

ˆf

2

i+1/2 =

1

6

(2fi + 5fi+1 − fi+2)

One can easily show (see Henrick et al. [54]) that the fifth-order flux at xi+1/2 for the upstreamcentered scheme is given by

ˆfi+1/2 =

1

60

(2fi−2 − 13fi−1 + 47fi + 27fi+1 − 3fi+2) (34)

Remark: The third-order fluxes of Eq. (33) and the fifth-order flux of Eq. (34) can also be determined

by using Table 2.1 in the 1997 paper by Shu [36].

By substituting the third-order approximations of Eq. (33) into Eq. (29) and equating corresponding terms, we obtain the weights

ω0 = 1/10, ω1 = 6/10, ω2 = 3/10 (35)

which give the desired fifth-order flux function in regions where the solution is smooth. The weights

sum to unity, and they are called the ideal weights. To obtain the flux ˆfi−1/2 for approximating the

first derivative of the inviscid flux function f, the stencil for ˆfi+1/2 is shifted to the left by one mesh

point. Substituting for the fluxes in Eq. (27) gives

∂f

∂x









x=xi

≈

1

60∆x

(−2fi−3 + 15fi−2 − 60fi−1 + 20fi + 30fi+1 − 3fi+2) (36)

= f

0

(x) + O(∆x

5

)

Even though the fluxes at xi−1/2 and xi+1/2 are only fifth-order accurate and the flux balance of

Eq. (27) is divided by ∆x, the resulting flux derivative is fifth-order accurate because the two stencils

approximating the fluxes only differ by a grid point shift, resulting in a canceling of the lowest order

truncation error terms.

13
3.2.1b 5th Order Non-Oscillatory Weights

So far we have only considered a linear scheme. The approximation of the inviscid flux derivative

given by Eq. (36) requires a smooth solution for high-order accuracy. If a discontinuity such as a

shock wave occurs in the flow field, oscillations will occur in the discrete solution. The full nonlinear

WENO scheme can be designed to give comparable high-order accuracy in smooth regions to that

of the linear scheme and prevent spurious oscillations at discontinuities. This is accomplished by

introducing appropriate nonlinear weights in Eq. (29) such that the weight of each stencil depends

upon its smoothness. In the nonlinear scheme, stencils that cross a shock wave are automatically

assigned an essentially zero weight, producing a one-sided upwind difference approximation. Such

an approach is similar to what is done with second-order schemes, which use a limiter for shock

capturing that produces a first-order upwind difference at shock waves. One can view the nonlinear

weights of the WENO scheme as playing the role of a limiter, and by their action, the scheme is

essentially non-oscillatory. In the neighborhood of a discontinuity the order of accuracy is reduced.

That is, if r is the order of accuracy of the weighted fluxes, then in regions where the solution is

smooth the order of accuracy of the scheme is 2r − 1 and in the vicinity of a shock wave the order

of accuracy is reduced to r.

Jiang and Shu [53] define the non-oscillatory weights as

ωk =

αk

Pr−1

k=0 αk

, αk =

dk

(βk + εw)

p

(37)

where the αk are normalized weights, dk are the ideal weights, and the parameter εw is a small

positive number used to bound the weights (i.e., to prevent division by zero). For r = 3 the exponent

p = 2 is chosen, since it proves to be adequate to obtain essentially non-oscillatory approximations.

The βk is a smoothness measure of the flux function on the stencil Sk, and it is defined by

βk =

Xr−1

l=1

∆x

2l−1

Z xi+1/2

xi−1/2

 

d

l ˆfk

dxl

!2

dx (38)

Since ˆf

k = a0 + a1x + a2x2

, the equation for βk becomes

βk = a

2

1∆x

2 +

13

3

a

2

2∆x

4

(39)

when the stencils are centered around xi = 0. Substituting for the coefficients a1 and a2 corresponding to each value of k we obtain the smoothness indicators

β0 =

13

12

(fi−2 − 2fi−1 + fi)

2 +

1

4

(fi−2 − 4fi−1 + 3fi)

2

β1 =

13

12

(fi−1 − 2fi + fi+1)

2 +

1

4

(fi+1 − fi−1)

2

(40)

β2 =

13

12

(fi − 2fi+1 + fi+2)

2 +

1

4

(3fi − 4fi+1 + fi+2)

2

which have Taylor series expansions at xi given by

β0 = f

0

i∆x

2 +



13

12

f

00

i

2 −

2

3

f

0

i

f

000

i



∆x

4 −



13

6

f

00

i

f

000

i −

2

2

f

0

i

f

0000

i



∆x

5 + O(∆x6

)

β1 = f

0

i∆x

2 +



13

12

f

00

i

2 +

2

3

f

0

i

f

000

i



∆x

4 + O(∆x6

) (41)

β2 = f

0

i∆x

2 +



13

12

f

00

i

2 −

2

3

f

0

i

f

000

i



∆x

4 +



13

6

f

00

i

f

000

i −

2

2

f

0

i

f

0000

i



∆x

5 + O(∆x6

)

14
with the prime denoting differentiation (e.g., f

000 = d3f /dx3

). Appropriate smoothness indicators

have the following behavior. In regions where the solution is smooth, the βk are small and about the

same size, resulting in weights ωk that are close to the ideal weights dk. If the approximation stencil

contains a discontinuity, then the associated smoothness measure βk is O(1), and the corresponding

weight ωk is small relative to the other weights.

When the solution is smooth, a requirement in choosing the appropriate βk’s is that the corresponding ωk’s approach the ideal weights at a sufficiently fast convergence rate. Henrick et al. [54]

and Borges et al. [55] derive the necessary and sufficient conditions to ensure the design order of

accuracy in smooth regions. These conditions are as follows:

Xr−1

l=1

(ω

±

k − dk) = O(∆x

2r

) (42)

Xr−1

l=1

Ak(ω

+

k − ω

−

k

) = O(∆x

r

) (43)

ω

±

k − dk = O(∆x

r−1

) (44)

where Ak is the coefficient (independent of ∆x) of the leading term in a Taylor series expansion

for ˆf

k

i±1/2

, and the superscript ± refers to this flux function. As indicated by Borges et al. the first

condition is always satisfied, since Pr−1

k=0 ω

±

k =

Pr−1

k=0 dk. A sufficient condition for fifth-order

convergence is given by

ω

±

k − dk = O(∆x

3

) (45)

Following the analysis of Henrick et al. for the smoothness indicator defined by Eq. (40), one can

show that if βk can be expressed as

βk = D[1 + O(∆x

r−1

)] (46)

where D is a non-zero constant independent of k, then ωk = dk + O(∆x

r−1

). When r = 3 we

readily see from the Taylor series expansions for the βk’s given in Eq. (41) that Eq. (46) is satisfied,

and thus, ωk = dk + O(∆x

2

). Borges et al. emphasize that the condition of Eq. (43) must be

satisfied to achieve the design order of convergence.

At critical points, which are points where the first or higher derivatives of the function f vanish,

these conditions are not necessarily satisfied using the βk given in Eq. (40). For example, if f

0 = 0,

we find from the Taylor series expansions of Eq. (41) that βk = D[1 + O(∆x)] for k = 0, 2, and

thus, ωk = dk+O(∆x), resulting in third-order convergence. In order to overcome this deterioration

in the convergence of the scheme due to critical points, Borges et al. introduced a different definition

of the measures of smoothness βk given by

β

z

k =



βk + εw

βk + τ5 + εw



, k = 0, 1, 2 (47)

where τ5 = |β0 −β2|, and the superscript z indicates that the quantity, a smoothness measure in this

case, is associated with the WENO-Z scheme of Borges et al. The corresponding weights ω

z

k

are

defined as

ω

z

k =

α

z

P

k

r−1

k=0 α

z

k

, αz

k =

dk

β

z

k

= dk



1 + 

τ5

βk + εw

q

, k = 0, 1, 2 (48)

The exponent q ≥ 1 increases the rate of convergence to the ideal weights for smooth regions of the

solution, and it accelerates the convergence rate to zero in the vicinity of discontinuities.

In practice, the question arises as to whether or not the ωk are sufficiently good approximations

of the ideal weights for smooth parts of the solution. The appropriate scaling of the smoothness

15
indicators βk is essential to ensure that they do not dominate the parameter εw in smooth regions and

produce weights that will result in too much numerical dissipation. This suggests that a given scaling

of βk can have an effect on the proper choice of εw. At the same time, caution must be exercised to

make certain that the choice of εw does not compromise the non-oscillatory behavior of the WENO

scheme. In numerical computations we have observed that if the parameter εw is too small relative

to the smoothness measures defined by either Eq. (40) or Eq. (47), the WENO scheme will have too

much numerical dissipation. Other difficulties can also occur when εw is not sufficiently large. Shen

and Zha [57] found that when εw is too small (e.g., 10−6) relative to βk, convergence and stability

problems can occur due to small disturbances in the βk. These disturbances result in oscillations of

the nonlinear weights which prevent convergence to the ideal (optimal) weights in smooth regions.

Thus, there are increases in the numerical dissipation. Shen and Zha considered a seventh-order

accurate scheme with the smoothness indicators of Borges et al. [55] and solved either the twodimensional Euler or Navier-Stokes equations for subsonic and transonic steady flows. By setting

εw = 10−2, Shen and Zha prevented convergence stall and showed improved accuracy.

Oscillations in the weights are not only caused by flow discontinuities. As pointed out by Shen

and Zha [58], strong variations in the smoothness measures βk can occur due to nonuniformities in

the computational mesh. If εw is very small (e.g., 10−6

), the βk will not always be less than εw,

even in smooth regions of the flow field. Here again, by increasing the magnitude of εw one can

avoid oscillations in the nonlinear weights.

In the present work we apply the weights of Eq. (47). For comparison purposes we also consider

the weights of Eq. (37). The influence of the parameter εw on accuracy is shown in the results

section.

3.2.1c Flux Splitting

Throughout the previous discussion we have assumed that the flux Jacobian (∂f /∂u) is positive.

In general, there is a need to allow for a sign change in the flux Jacobian. This is usually done by

splitting the flux into positive and negative parts. In the present formulation for the WENO scheme

a global Lax-Friedrichs splitting is used. Thus,

f

±(u) = 1

2

(f(u) ± αu) (49)

where α = max|∂f /∂u|, and the maximum is taken over the complete solution field. A local

splitting has also been considered. For the decay of turbulence problem considered herein, the local

splitting did not have any significant effect on the results.

3.2.2 Approximating the Viscous Term

Here we discuss the approximation of the physical diffusion terms in the Navier-Stokes equations

by considering the derivative of the viscous flux function f

v

(ux) in Eq. (23). The derivative of the

viscous flux function is defined by

df v(x)

dx =

d

dx 

a(x)

du(x)

dx 

(50)

where in general, the coefficient a(x) = a[u(x), x]. Since we want to construct a conservative

difference scheme, this derivative is approximated by

df v(x)

dx ≈

f

v

i+1/2 − f

v

i−1/2

∆x

(51)

Thus, one needs to determine the interface approximations of a(x) and u(x) at xi+1/2 and xi−1/2.

With the fifth-order WENO scheme, the diffusion term is approximated with a conservative fourth16
order finite difference. When evaluating possible ways to determine these discrete interface functions, consideration must be given to not only accuracy and stability but also damping of all error

modes. We now consider two ways to compute ai±1/2 and ui±1/2.

One possible method of discretizing the diffusion term is to first compute the discrete f

v

(x) at

the solution points. A fourth-order approximation of f

v

(x) at xi, an interior point of the domain, is

given by

(f

v

)i =

ai

12∆x

[ui−2 − 8ui−1 + 8ui+1 − ui+2] (52)

Then the approximation for the derivative of f

v

(x)is evaluated when constructing the approximation

for the derivative of f(x). For the fifth-order linear form of the WENO scheme (i.e., ideal weights)

we obtain

df v

dx









x=xi

≈

1

60∆x

[−2(f

v

)i−3 + 15(f

v

)i−2 − 60(f

v

)i−1 + 20(f

v

)i (53)

+ 30(f

v

)i+1 − 3(f

v

)i+2]

=

1

720∆x

2

[−2ai−3ui−5 + (16ai−3 + 15ai−2)ui−4

− (120ai−2 + 60ai−1)ui−3 + (−16ai−3 + 480ai−1 + 20ai)ui−2

+ (2ai−3 + 120ai−2 − 160ai + 30ai+1)ui−1

− (15ai−2 + 480ai−1 + 240ai+1 + 3ai+2)ui

+ (60ai−1 + 160ai + 24ai+2)ui+1 + (−20ai + 240ai+1)ui+2

− (30ai+1 + 24ai+2)ui+3 + 3ai+2ui+4]

If a(x) and u(x) are smooth functions, then the approximation of Eq. (53) is fourth-order accurate.

While this method for differencing the derivative is convenient for implementation and operation

count, it does have the undesirable consequence that the π-mode is not damped. One can easily

show this by treating the discrete product of a(x) and u(x) as a single variable, Fourier transforming

the derivative, and evaluating the transformed discrete operator at a Fourier angle equal to π.

An alternative method that provides a compact difference stencil and ensures damping of the

π-mode is given by Kamakoti and Pantano [59]. This approximation enforces discrete conservation.

Let B be a matrix of stencil coefficients for approximating the viscous term. Then, the flux form is

given by

df v(x)

dx









x=xi

≈

1

∆x

2

Xp

m=−p

Xp

n=−p

Bm+p+1,n+p+1 ai+mui+n =

f

v

i+1/2 − f

v

i−1/2

∆x

(54)

where

f

v

i+1/2 =

1

∆x

Xp

m=−p+1

Xp

n=−p+1

Mm+p,n+p ai+mui+n (55)

with M denoting the matrix of the unknown coefficients for the viscous flux to construct the difference stencil. The widths of the flux stencil and the stencil for the variable coefficient second-order

derivative are 2p and 2p + 1, respectively. One uses Taylor series expansions to determine the

elements of the matrix M.

If we expand Eq. (55) for a fourth-order approximation, we obtain after rearranging terms

f

v

i+1/2 =

1

∆x

[(M11ai−1 + M21ai + M31ai+1 + M41ai+2)ui−1 + (56)

(M12ai−1 + M22ai + M32ai+1 + M42ai+2)ui +

(M13ai−1 + M23ai + M33ai+1 + M43ai+2)ui+1 +

(M14ai−1 + M24ai + M34ai+1 + M44ai+2)ui+2]

17
The stencil matrix M is given by

M =





















1

8 −

1

6

1

24 0

−

1

6 −

3

8

2

3 −

1

8

1

8 −

2

3

3

8

1

6

0 −

1

24

1

6 −

1

8





















(57)

The resulting approximation of Eq. (54) for the viscous term is compact, requiring just five points

for fourth-order accuracy. Thus, it has a reduced error relative to the approximation of Eq. (53), and

it has the advantage that the π-mode is damped.

The discretization method used to obtain Eq. (53) was employed in the original ENO (Essentially

Non-oscillatory) code of Atkins [60] that was modified by Balakumar et al. [61] to create a WENO

version. While the π-mode will be damped when there is deviation between the WENO weights and

the ideal weights, this mode may not be adequately damped in regions where the diffusion terms

dominate or the ideal weights are recovered. For this reason, along with reducing the truncation

error and ensuring that this error is purely dissipative, the discretization method for Eq. (53) should

be replaced by that for Eq. (54).

In this initial work we have retained the discretization method for the diffusion terms used by

Atkins [60] and Balakumar et al. [61]. For the WENO scheme the polynomial that interpolates

f

v

(x) at the solution points to the cell interfaces requires the smoothness indicators of Eq. (40) to

determine the appropriate non-oscillatory stencil weights. We should mention that periodic boundary conditions are imposed for the isotropic turbulence problem considered in this paper. This is

accomplished in the WENO code by introducing additional points outside the physical domain.

3.2.3 Higher Order Schemes

Higher order WENO schemes can be derived by following the steps presented in the first two subsections of Section 3.2.1 for the fifth-order scheme, starting with Eq. (29) and the appropriate r. For the

seventh-order and ninth-order schemes the fluxes in the summation are fourth-order and fifth-order,

respectively. Thus, each flux ˆf

k

i+1/2

for the seventh-order and ninth-order schemes is determined

by third-degree and fourth-degree polynomials, respectively. The stencils for theses schemes that

define the corresponding WENO flux ˆfi+1/2, the ideal weights dk, and the smoothness indicators

βk are presented in the paper by Balsara and Shu [56]. For completeness, these quantities are included in Appendix A. It should be noted that the smoothness indicators in Appendix A include the

scaling that arises when using Eq. (38), whereas those of Balsara and Shu do not. For extremely

small values of the εw in Eq. (37) and Eq. (48) this does not matter. However, it does have an effect

for larger values of εw such as 10−2.

For the seventh-order and ninth-order schemes the f

v

(x) of the diffusion term is discretized at

interior points with sixth-order and eighth-order central differences, respectively. Stencils for these

approximations can be obtained using Table 2.1 of Shu [36], which is supplemented by Table A1 in

Appendix A to include eighth-order. The sixth-order approximation is given by

(f

v

)i =

ai

60∆x

(−ui−3 + 9ui−2 − 45ui−1 + 45ui+1 − 9ui+2 + ui+3) (58)

and the eighth-order one is given by

(f

v

)i =

ai

840∆x

(3ui−4 − 32ui−3 + 168ui−2 − 672ui−1 + 672ui+1 (59)

− 168ui+2 + 32ui+3 − 3ui+4)

The interpolating polynomials for obtaining the viscous fluxes at the cell interfaces are constructed

with the fluxes and WENO weights given in Appendix A.

18
3.2.4 Time-Stepping Scheme

To advance the solution in time, a third-order Runge-Kutta (RK) scheme that satisfies the total

variation diminishing (TVD) property is used. A scheme is said to be TVD if

T V (u

n+1) ≤ T V (un

), T V =

X

i

|ui+1 − ui| (60)

The RK scheme has three stages given by

u

(1) = un + ∆tLun

,

u

(2) =

1

4

(3u

n + u(1) + ∆tLu(1)), (61)

u

(3) =

1

3

(u

n + 2u(2) + 2∆tLu(2)),

u

n+1 = u(3) (62)

where L is the discrete spatial operator, and the superscript n denotes time level (time t = n∆t).

4 Numerical Results

4.1 Problem Definition

4.1.1 Decay of Isotropic Turbulence (DOIT)

In the classical problem of decay of isotropic turbulence, an initial turbulent field is allowed to

decay in a “box-like” domain with periodicity imposed on all six of its sides. The turbulent statistics

are collected and compared with experiments or with other numerical studies. Particular focus is

often placed on the decay of the total turbulent kinetic energy with time, the initial growth followed

by subsequent decay of the total turbulent dissipation with time, and the energy spectrum at given

instants in time.

This classical problem is a good discriminator for SGS models and numerical schemes: poor

models or overly-dissipative schemes tend to diffuse the energy-containing structures, and consequently miss many of the characteristic behaviors that can be accurately captured by well-resolved

DNS. See, for example, the experimental results of Compte-Bellot and Corrsin [62], and the numerical studies of Clark et al. [63], Vreman et al. [64], Mansour and Wray [65], Spyropoulos and

Blaisdell [66], Seror et al. [67], Hughes et al. [68], Hansen et al. [69], and Thornber et al. [42].

In the current study, the domain uses side lengths of 2π. Various Cartesian grid sizes are employed, ranging from the coarsest with 333 grid points (323 grid cells) to the finest with 1933 grid

points (1923 grid cells).

To evaluate the numerical solutions to this problem, we use DNS results computed with a spectral scheme. Spectral solutions provide a good measure of accuracy due to their low discretization

errors. They have the additional advantage of much faster convergence than finite-difference or

finite-volume methods (i.e., exhibiting exponential rather than algebraic convergence with mesh refinement). These computations were performed with the same spectral computer code used to obtain

the results presented in the paper by Zang et al. [70]. The conditions of the problem considered in

the present work are for one of the cases included in the paper of Zang et al. In the DNS computer

code, the spatial discretization is Fourier collocation with evenly spaced points in all coordinate

directions. The code employs a time-step splitting technique that allows sound waves to be integrated exactly in time [71]. The other terms are advanced with a third-order Runge-Kutta explicit

time-stepping scheme. An isotropic truncation is applied at the end of each Runge-Kutta stage. This

ensures that only the resolvable modes (i.e., |κ| ≤ N/2 where κ is the Fourier wavenumber vector)

are nonzero. The semi-implicit scheme provides important computational advantages. For example,

19
at small fluctuating Mach numbers it can allow the use of time steps several times larger than those

possible with a fully explicit scheme.

4.1.2 Initial Solution

The method for creating the initial turbulent flowfield condition is taken directly from Erlebacher

et al. [71, 72], and is similar to a method introduced by Passot and Pouquet [73]. Random fields in

physical space are generated for density, velocity, and temperature (ρ, u ̃, and T ̃). Subsequently the

velocity field is decomposed in Fourier space into incompressible and compressible components:

u ̃0(x) =  ̃u

I

0

(x) +  ̃u

C

0

(x) (63)

The u ̃

C

0

(x) has mean value of zero by construction, and u ̃

I

0

(x) is divergence free. These components

are transformed back to physical space and re-scaled in order to impose a prescribed autocorrelation

spectrum. The spectrum

E(κ) = (κ + 1)4exp −2

(κ + 1)2

κ

2

0



, κ ≥ 1 (64)

is imposed on each of the variables ρ, u ̃

I

0

, u ̃

C

0

, and T ̃ (for κ = 0 E(κ) = 0). Also, the initial value

for the ratio of the compressible kinetic energy to the total kinetic energy, which is often called χ, is

set to 0.2 (see Zang et al. [70]). Note that Eq. (64) is slightly different than that reported in Zang et

al. [70]. The reference pressure is chosen (based on a reference density, ρref , and reference velocity,

U ̃

ref ) to be pref = ρrefU ̃ 2

ref , so the initial nondimensional pressure is:

p =

ρT ̃

γM2

ref

(65)

The reference Mach number

Mref =

U ̃

q

ref

γRT ̃

ref

(66)

is taken to be Mref = 0.2 for all the current simulations. The Reynolds number, which is defined as

Re = (ρrefU ̃

refLref )/μref , is taken to be 250, and the Prandtl number and ratio of specific heats

have the standard values of P r = 0.72 and γ = 1.4, respectively. The reference length Lref in this

case is unity.

The Taylor microscale Reynolds number is defined by:

Reλ =

hρi

hμi

λU ̃

rms (67)

where the angle brackets indicate an average over all cells in the domain, and

λ =

1

3

(λx + λy + λx) (68)

λx =

s

hu ̃

2i

h(

∂u ̃

∂x )

2i

λy =

s

hv ̃

2i

h(

∂v ̃

∂y )

2i

λz =

s

hw ̃

2i

h(

∂w ̃

∂z )

2i

(69)

U ̃

rms =

p

hu ̃

2 +  ̃v2 +  ̃w2i (70)

For all the cases herein, the initial Reλ is approximately 125. The x-, y-, or z-direction Taylor

microscale Reynolds numbers can also be found. For example, the x-direction Taylor microscale

Reynolds number is defined by:

Reλ,x =

hρi

hμi

λx

p

hu ̃

2i (71)

The initial Reλ,x = Reλ,y = Reλ,z ≈ 75.

20
4.1.3 Turbulence Statistics

This section defines some of the turbulence statistics. The total turbulent kinetic energy is:

k =

1

2

X

cells



u ̃

2 +  ̃v2 +  ̃w2



(72)

The total turbulent dissipation is:

ε = 2

hμi

hρi

X

cells

SijSij (73)

The skewness factor of the velocity derivative, which supplies a measure of the nonlinear vortex

stretching, is given by:

Sk =

1

3

"

h(

∂u ̃

∂x )

3

i

[h(

∂u ̃

∂x )

2i]3/2

+

h(

∂v ̃

∂y )

3

i

[h(

∂v ̃

∂y )

2i]3/2

+

h(

∂w ̃

∂z )

3

i

[h(

∂w ̃

∂z )

2i]3/2

#

(74)

Some discussion about the skewness factor and its relationship to the numerical scheme is given in

Appendix B.

The enstrophy is defined by:

enstrophy = 1

2

h|∇ × u|

2

i (75)

The turbulence spectra is computed as follows. Given (at an instant in time) the velocities everywhere in the field, a 3-D fast Fourier transform (FFT) is performed on each velocity component. The

energy E0(k) within each energy “bin” is computed by summing over the squares of all the Fourier

indices that contribute to that energy level. Then, the result is normalized via:

E(κ) = E

0

(κ)/

Z ∞

0

E

0

(κ)dκ(76)

4.2 Results with Upwinding and Central Differencing

Table 1 summarizes the computer runs performed using CFL3D. Both upwinding and central differencing methods were employed, and the effects of grid size, time step, and Smagorinsky filter constant cs were explored. The time step, ∆t

∗

, was nondimensionalized using ∆t

∗ = ∆tUref /Lref ,

where Uref and Lref are the reference velocity scale and reference length scale, respectively. This

time scale is a measure of “eddy turnover time” in terms of mean flow quantities. Most of the runs

were performed at a nondimensional time step of ∆t

∗ = 0.002. As can be seen in the table, one run

was performed at a time step that was one fourth the nominal time step. Although not shown, results

using the two different time steps were essentially identical, so the larger time step, ∆t

∗ = 0.002,

was deemed to be fine enough to accurately capture the temporal development with the current

scheme in CFL3D.

4.2.1 Implicit LES (No Subgrid Model)

Results using CFL3D with the upwind-biased κ = 1/3 scheme are shown in Figs. 3(a)–(c) for four

different grid cell sizes. Note that the k and ε were normalized by their initial values. All cases

except the 1923 were carried out to total time of t

∗ = 10. The main point to be made is that on the

coarsest grid (323) the upwind result is very poor, particularly in the fact that it shows no increase

at all in ε for 0 < t∗ <≈ 1.5. As will be shown later, ε should increase to a level of nearly 2 (near

t

∗ = 1.3) before dropping off. This rise in ε occurs during the energy redistribution phase of the

21
Table 1. Summary of CFL3D Cases

κ Grid size (# cells) ∆t

∗

cs

1/3 (upwind) 323 0.002 0

1/3 (upwind) 643 0.002 0

1/3 (upwind) 1283 0.002 0

1/3 (upwind) 1923 0.002 0

1 (central) 323 0.002 0

1 (central) 643 0.002 0

1 (central) 1283 0.002 0

1 (central) 1923 0.002 0

1 (central) 323 0.002 0.2

1 (central) 643 0.002 0.2

1 (central) 1283 0.002 0.2

1 (central) 323 0.002 0.1

1 (central) 323 0.002 0.4

1 (central) 323 0.0005 0.2

simulation. The skewness using upwind with no model ends up between about −0.2 and −0.5. For

incompressible flows the experimental value for isotropic turbulence is around −0.4 (Lesieur [74]),

and the DNS value at Reλ = 40 is −0.5 (Erlebacher [71]).

Isocontours of vorticity are shown at t

∗ = 5 on three of the grids in Fig. 4. (Contours on the

1923 grid are not shown; visually there are more small scale structures than for the 1283 grid, as

expected. The entire cube, 2π on a side, is shown here along with the subdivisions used for parallel

processing.) However, in the literature, the so-called Q-invariant or Q-criterion is often preferred for

visualizing LES results. See, for example, Piomelli et al. [75] and Slomski and Chang [76]. The Q

variable is the second invariant of the velocity gradient tensor:

Q =

1

2

(ΩijΩij − SijSij ) = −

1

2



∂ui

∂xj

∂uj

∂xi



(77)

The Q-invariant is useful because it tends to highlight the smaller-scale structures. Figure 5 shows

the Q-invariant plotted for the upwind results. It is evident that using upwind on the 323 grid there

is little eddy content, whereas on finer grids there are additional smaller scale structures present.

Results using CFL3D with second-order central differencing (κ = 1) and no SGS model (cs =

0) are shown in Figs. 6(a) – (c) for four different grid sizes. These results can be thought of as

under-resolved DNS computations because there is no numerical or SGS model dissipation, and the

grids are too coarse for a second-order finite-volume scheme to resolve down to the Kolmogorov

dissipation length scale. Comparing Fig. 6(b) with Fig. 3(b), the central differencing yielded better

results for temporal development of ε on a given grid size (the peak level was closer to the expected

value of approximately 2). This improvement is to be expected because upwinding is inherently

more dissipative. For example, using central differencing, the peak normalized ε reached nearly 1.4

on the 643 grid, which was almost as good as the upwind result on the 1283 grid. However, note

that the skewness tended to be smaller in magnitude using central differencing compared to upwind.

The reason for this underprediction is not known.

Isocontours of the Q-invariant using central differencing are shown at t

∗ = 5 on three of the

grids in Fig. 7. Comparing to Fig. 5, it is easy to see that central differencing yielded greater eddy

content than upwind differencing on a given grid size.

Figure 8 shows energy spectra for eight different results at t

∗ = 5. The theoretical slope of

κ

−5/3

from Kolmogorov theory (known as the Kolmogorov −5/3 law) is also indicated in the plot.

22
In the inertial subrange where the energy cascade transfers energy from large to small scales, the

spectrum should follow this slope. However, very high Reynolds numbers are required in order to

see evidence of the inertial subrange. The current Reynolds number of 250 is far too low for a region

with constant slope to emerge. See Appendix C for additional discussion. Nonetheless, in this and

subsequent plots of the energy spectrum, the κ

−5/3

line will continue to be shown to help provide a

reference for the spectrum. In all cases shown in Fig. 8, results with central difference were better

than for upwind (less rapid drop-off in energy with increasing κ). The greater decay in the upwind

result is a consequence of the increased numerical dissipation with the upwind differencing, which

results in greater smearing of the small scale structures associated with the higher wavenumbers.

However, the finer the grid, the closer the central and upwind results were to each other. Figure 9

shows the progression of the energy spectrum with eddy turnover time for central differencing (no

SGS model) on the 1283 grid.

4.2.2 Explicit LES with Smagorinsky Model

Central differencing was also used in conjunction with the Smagorinsky model, using CFL3D.

Figs. 10(a) – (c) show the statistical information over time. Generally, the higher cs, the worse

the result (on the 323 grid). The best result was achieved using no model at all (cs = 0). This

conclusion is also reflected in the spectrum results at t

∗ = 5, which is shown in Fig. 11. The higher

the cs coefficient, the more the energy dropped off at the higher wavenumbers.

On finer grids, the effect of cs was not as great, as seen in Figs. 12(a) – (c) and Fig. 13. For

example, on the 1283 grid, the energy spectrum result using cs = 0.2 is only slightly worse than the

result with cs = 0. Isocontours of the Q-invariant using central differencing with the Smagorinsky

model and cs = 0.2 are shown at t

∗ = 5 in Fig. 14.

Upwind was not run in conjunction with the Smagorinsky model (cs > 0), because the combination would yield even more dissipative results than upwind ILES (cs = 0).

4.3 Implicit LES Results with WENO and WENO-L

The focus in this section is on the effect of higher order accuracy on the turbulence statistics. Fifthorder, seventh-order, and ninth-order schemes are all considered. To ensure that the best possible

results are obtained with the higher order approximations, we have switched off the nonlinear part

of the WENO scheme for the initial set of computations. The linear scheme is designated WENO-L.

In the second part of this section we investigate the nonlinear scheme, which is simply designated

as WENO. The two ways of computing the nonlinear weights, ωk, and smoothness measure, βk, described in Section 3.2.1, are considered. In addition, the effect of the parameter εw is also evaluated

and discussed. It is demonstrated by numerical testing that nearly the same solutions can be obtained

with the nonlinear WENO scheme as for the WENO-L scheme if a sufficiently large value for εw

in the non-oscillatory weights is chosen. For all high-order schemes the effect of mesh resolution

on the computed solutions is examined. The grid density varied from 323to 1923; these are the

same grids used for all results in the previous section. Results obtained with the WENO schemes

are also compared with the second-order central (cs = 0 and cs = 2) and third-order upwind biased

schemes. In addition, they are compared with those from the spectral scheme of Zang et al. [70]. For

all simulations a ∆t

∗ = 0.002 was used. Reducing the time step by a factor of four, which reduces

temporal error by a factor of 64, produced essentially the same numerical results. All of the WENO

and WENO-L calculations were performed without an explicit SGS model.

4.3.1 Linear WENO Scheme

Figures 15–17 show the turbulence quantities computed using the WENO-L scheme with different

orders of spatial accuracy. Comparing the decay in the turbulent kinetic energy of the dissipative

higher order schemes with that of the nondissipative central scheme, which was shown in Fig. 6(a),

23
Table 2. Normalized Dissipation Peak

Scheme Order Grid size (# cells) Dissipation Peak

Central 2 323 1.10

Central 2 643 1.40

Central 2 1283 1.64

Central 2 1923 1.75

Upwind 3 323 1.00

Upwind 3 643 1.10

Upwind 3 1283 1.46

Upwind 3 1923 1.65

WENO-L 5 323 1.00

WENO-L 5 643 1.37

WENO-L 5 1283 1.75

WENO-L 5 1923 1.87

WENO-L 7 323 1.08

WENO-L 7 643 1.56

WENO-L 7 1283 1.84

WENO-L 7 1923 1.91

WENO-L 9 323 1.17

WENO-L 9 643 1.69

WENO-L 9 1283 1.96

WENO-L 9 1923 2.02

we observe that the dissipative schemes approach the solution on the 1923 grid from below while

the nondissipative scheme approaches from above, as expected. Moreover, it is evident that the

fifth-order scheme is highly dissipative on the coarser meshes, especially on the 323 mesh. In (b)

of Figs. 15 – 17, we see that the peak value of the normalized dissipation, which occurs during the

energy redistribution phase of the simulation, increases both with grid refinement and with order of

accuracy. The peak values for the different approximation orders and the range of mesh densities

are given in Table 2, along with results from the earlier central (cs = 0) and upwind computations.

These results are also shown in Fig. 18. The largest energy dissipation (ε) peaks, obtained with the

seventh-order and ninth-order schemes on mesh densities 1283and 1923, are roughly in the range

of 1.8 to 2.0 and occur at t

∗ ≈ 1.30. On the 1283 grid, the peak value in ε is about 20% larger

than the one computed with second-order central differencing. Like the third-order upwind biased

scheme, the fifth-order WENO-L scheme does not exhibit any dissipation peak on the 323 grid, even

though its dissipation is somewhat smaller (as seen by comparing the turbulent kinetic energy levels

in Figs. 3(a) and 15(a)). Even the seventh-order scheme only has a very small peak on this grid, and

this peak is about the same size as that for the central difference scheme. The peak in dissipation for

the fifth-order scheme exceeds 1.8 on the 1923 grid, whereas for the seventh-order and ninth-order

schemes it is greater than 1.8 on the 1283 grid. These results indicate that the resolution on the 323

grid is not adequate for any of the schemes. The numerical viscosity is not sufficiently small on such

a coarse grid.

We now consider the skewness factor (Sk) plots in Figs. 15 – 17. On the finest mesh the quasiequilibrium value is between -0.5 and -0.6 for the different approximation orders of the WENO-L

scheme. It is nearly -0.6 using the ninth-order scheme on the 1923 grid. This is about the same Sk

24
obtained with the third-order upwind biased scheme on the finest mesh.

In Fig. 19, the energy spectra are shown for computations with different orders of accuracy and

grid sizes. The fifth-order results are comparable to those obtained with the third-order upwind

biased scheme, although the fifth-order results are more dissipative at the highest frequencies. In

comparing the spectra for the higher order schemes with the fifth-order scheme, we observe roughly

an order of magnitude less drop in the energy decay at the the cut-off wavenumber (N/2, where N is

the number of mesh points), reflecting the reduced dissipation of the small scales. Figure 20 clearly

shows how the high-frequency content improves with increasing mesh density for the ninth-order

scheme. In Fig. 21, we see that the spectra obtained on the 1283 grid with the ninth-order WENO

and the second-order central (cs = 0) schemes are nearly the same. The fifth-order WENO and

third-order upwind spectra are also similar to each other

Figure 22 displays the isocontours of the nondimensional Q-invariant for the ninth-order WENOL scheme computations for various grid sizes. Comparing with the corresponding contours obtained

with second-order central differencing (cs = 0) in Fig. 7, we note that in these plots the two schemes

appear to have similar eddy content on the 1283 grid. However, comparing with the explicit LES

scheme with central differencing and cs = 0.2 (see Fig. 14), we observe that there is a significant

increase in the eddy content for ninth-order over the range of meshes.

4.3.2 Nonlinear WENO Scheme

Appropriate determination of the nonlinear weights is essential since these weights allow nonoscillatory shock capturing and determine the rate at which the nonlinear WENO scheme returns

to a linear form. We first evaluate the two methods for computing the nonlinear weights, ωk, and

smoothness measures, βk. Figure 23 compares the two methods discussed in Section 3.2.1. In the

figure, α1 refers to the method of Jiang and Shu (Eqs. (37) and (38)), and α2 refers to the method

of Borges et al., Eq. (48)). For each method, two values of εw, the parameter for ensuring that

βk is bounded, are considered. The computational results are for the fifth-order and seventh-order

schemes on the 1283 grid. We observe in the figure a slight improvement in the prediction of the

peak value of ε with the method of Borges et al. (α2). For the larger εw = 10−2, the energy spectrum

shows that the WENO scheme is dissipative at the highest frequencies, whereas there is some energy

buildup at the highest frequencies for the smaller εw = 10−6. This behavior is most prominent for

the fifth-order scheme.

In Fig. 24, the method of Borges et al. is considered, and the effect of varying εw on the turbulence statistics using a 643 grid is shown. The WENO and WENO-L results are nearly the same

when εw = 10−1. However, such a large value could easily adversely affect the nonlinear behavior

(i.e., shock capturing properties) of the WENO scheme. Shen and Zha [57, 58] have demonstrated

for transonic fluid dynamics problems with shock waves that a value of εw = 10−2

is satisfactory.

It still remains, however, to show that this value of εw is suitable for a wide range of flow problems

with various levels of complexities and Mach numbers. To reduce the possibility of energy buildup

at the highest frequencies that can result in instability, the method of Borges et al. with εw = 10−2

will be used for all subsequent results.

Figure 25 compares the dissipation computed with the WENO-L and WENO schemes for different grids and orders of discrete approximation. In general, the linear and nonlinear results are

quite close, with the greatest difference occurring on the coarsest (323) grid. A similar comparison

is shown in Fig. 26 for the energy spectrum. There is no discernible difference between the linear

and nonlinear results.

We now consider how the ninth-order WENO scheme compares with other previously presented

LES schemes. Figure 27 displays the turbulence statistics for the ninth-order WENO ILES scheme,

second-order central scheme with the Smagorinsky SGS model (cs = 0.2), and the third-order

upwind biased ILES scheme. Results are given for the 643and 1283 grids. The largest differences

in the turbulent kinetic energy (k) decay for the two grids occur with the upwind scheme. With the

WENO scheme the k decay for the 643 grid closely tracks that for the 1283 grid. When comparing

25
the dissipation (ε) peaks, we see that the peak on the 1283 grid with WENO scheme is significantly

larger than the peaks for the lower-order LES schemes. Furthermore, the ε peak on the 643 grid

with the WENO scheme is approximately 8% and 15% higher than the peaks on the 1283 grid with

the central and upwind schemes, respectively. These results indicate that more than an order of

magnitude increase in the grid density would be required to produce similar ε peaks with the lower

order LES schemes. In Fig. 27(c) a comparison of the LES schemes is given for the energy spectrum

at t

∗ = 5. The WENO result for the 643 grid follows the fine grid energy spectrum until between

κ = 15 and κ = 20, where it begins to show a more rapid decrease in the energy with increasing κ.

As expected, the lower order LES schemes exhibit a much earlier departure (roughly κ = 5) from

the fine grid result with WENO.

A comparison of the turbulence statistics at t

∗ = 2 for the ninth-order WENO ILES scheme,

second-order central scheme (cs = 0), and a spectral scheme is presented in Fig. 28. In the absence

of numerical dissipation (either through discretization or filtering) and an explicit SGS model, the

central scheme relies on the dissipative effect of the physical viscous terms for stability. Thus, the

computation with the central scheme is essentially an under-resolved DNS. One could view this

solution as representing a bound on the best possible result with a second-order scheme for an

explicit SGS model (i.e., in the limit as cs → 0). The results for the spectral scheme are from DNS

computations with the computer code of Zang et al. [70]. In Fig. 28a, we observe a similar decay

rate of the turbulent kinetic energy for all schemes. Figure 28(b) shows that the peak dissipation

for the ninth-order WENO scheme on the 1923 grid is approximately 4% high compared with the

DNS, whereas the second-order central scheme is approximately 10% low. Although not shown, a

simulation on a 3843 grid was performed to determine if there is any further increase in the peak

dissipation with the WENO scheme. Essentially the same peak value was obtained on the finer grid.

In Fig. 28c, we observe fairly good agreement between the energy spectra computed with the WENO

and spectral schemes, except at the highest frequencies where the WENO scheme exhibits a faster

decay in κ due to dissipative effects. There is also some difference at the lowest frequencies. The

source of this difference appeared to be a consequence of variations in the methods used for the two

schemes to determine the spectrum from the simulation data. The second-order central scheme on

the 1923 grid exhibits less energy content at the highest wave numbers than the ninth-order WENO

scheme, and thus, has poorer agreement with the DNS spectrum.

The computational cost of the WENO scheme with increasing approximation order is especially important due to the high operation count of the scheme. A comparison of the results for the

ninth-order and fifth-order schemes indicates that the ninth-order scheme requires half as many grid

points in each coordinate direction as the fifth-order scheme to obtain comparable resolution. We

compared the computational (CPU) times of these schemes for the isotropic turbulence decay simulation using the 643 grid for the ninth-order scheme and the 1283 grid for the fifth-order scheme.

All of the simulations were performed with parallel processing. To have approximately the same

parallel efficiency for both schemes, we used more processors for the fifth-order scheme on the

1283 grid than for the ninth-order scheme on the 643 grid (i.e., 64 and 16 processors, respectively).

This comparison showed that the ninth-order scheme required about 17% less computing time than

the fifth-order scheme for comparable resolution. We should point out that Shi et al. [77] showed

that for problems with discontinuities and complex solution structures (e.g., double Mach reflection

and Rayleigh-Taylor instability problems) that a ninth-order WENO scheme required half as many

mesh points in each coordinate direction to obtain comparable resolution (i.e., similar numerical

viscosity) as the fifth-order scheme. Furthermore, they showed that with parallel processing that the

ninth-order scheme could reduce the CPU time between 66% and 80% with their implementation

of the WENO scheme. Their results suggest that it may be possible to significantly improve the

present implementation of the WENO scheme, especially since no attempt was made to optimize it.

In addition, their results emphasize the capability of the high-order WENO schemes to capture small

scale structures of complicated flows that would require much finer grids with lower order schemes.

This is especially important in the simulation of turbulence.

26
4.3.3 Future Work: Improved WENO Scheme

As we have seen in the previous section, the conventional fifth-order WENO scheme of Shu [36]

with either the weights of Jiang and Shu [53] or Borges et al. [55] is too dissipative. In addition,

there is no energy stability proof (guaranteeing all eigenvalues lie in the left half of the complex

plane) for the conventional WENO scheme. To improve this situation Yamaleev and Carpenter [78]

– [80] have introduced energy stable modifications for existing WENO schemes. Such modified

schemes are called energy stable WENO (ESWENO) schemes. By modifying the nonlinear weights

of the conventional WENO schemes, Yamaleev and Carpenter achieved a reduction of the nonlinear

dissipation in the neighborhood of discontinuities, while retaining an energy stability of the spatial discretization. With their systematic methodology (see [78] and [80]) they have constructed

ESWENO schemes for both upwind biased and centered finite differences. Thus, stable higher

order central difference schemes can be used in smooth regions of the flow field, eliminating the

dissipation in these regions.

An essential constraint in minimizing numerical errors produced by dissipation is to enhance

the stability of the spatial discretization. For a given order of approximation, one approach that

can lead to not only reduction of numerical dissipation but also minimization of aliasing error is to

use a skew-symmetric form of the convective operator (e.g., see Kravchenko and Moin [81]). This

form of the convective operator is defined by a linear combination of the conservative (divergence)

and non-conservative operators using a central difference discretization. Such a discrete operator

is conservative (i.e., can be expressed in terms of telescopic fluxes), and in the absence of external

forces and dissipation, it produces a kinetic-energy preserving scheme. With kinetic energy conservation, one can show that the scheme is energy stable, which means that without dissipation the

energy of the discrete system cannot increase (with dissipation the energy decreases). Moreover,

in conserving kinetic energy there is no spurious production or dissipation of kinetic energy due to

the discretization of the nonlinear convection. Some examples of papers with schemes based on a

skew-symmetric operator are those by Lilly [82], Ducros et al. [83], Verstappen and Veldman [84],

Feiereisen et al. [85], and Kok [86].

While the ESWENO formulation and alternative convective operators can significantly reduce

numerical dissipation, and thus, improve the representation of the entire energy spectrum, the dispersion errors still need to be controlled. As Ghosal [87] has demonstrated, even a centered scheme

with eighth-order accuracy requires some technique such as prefiltering to effectively remove such

errors so that the modeling errors dominate. One possible approach for minimizing the numerical

dispersion errors at high wavenumbers is to construct a central difference approximation with the

dispersion-relation preserving (DRP) scheme of Tam and Webb [88]. The DRP scheme uses a larger

stencil than necessary for a given order of approximation to create an additional degree of freedom

to allow minimization of dispersion.

In future work, we plan to consider both an ESWENO scheme and the skew-symmetric form

of the convective operator. For the viscous terms we will replace the current discretization with the

compact approximation of Kamakoti and Pantano [59], which ensures damping of the π-mode. We

will also consider minimization of dispersion error.

5 Concluding Remarks

Numerical simulations of decaying homogeneous isotropic turbulence have been performed with

both low-order and high-order spatial discretization schemes. Second-order central and third-order

upwind biased differencing have been used for the low-order schemes. The WENO framework

has been applied to obtain fifth-order, seventh-order, and ninth-order discretizations. Large eddy

simulation results have been presented for the various schemes on four grids, ranging from 323

to 1923. Turbulence statistics, such as kinetic energy, dissipation, and skewness, along with the

energy spectra from simulations of the decaying turbulence problem have been used to examine and

27
compare the schemes.

Through the simulation results for the fifth-order WENO scheme and the second-order central

scheme, with or without the Smagorinsky SGS model (cs ≤ 0.2), we have shown that the fifthorder WENO is too dissipative on the coarsest grid considered (323

). While the fifth-order WENO

was still too dissipative on the on the 643 grid when compared with the central scheme without a

SGS model, this was not the case when compared with the central scheme with a SGS model and

cs ≥ 0.2. On the coarser grids, the fifth-order WENO scheme was somewhat less dissipative than

the third-order upwind biased scheme. Comparisons of the fifth-order WENO scheme with the lower

order LES schemes have indicated that the fifth-order WENO has reasonable levels of dissipation

on the finer grids (1283and 1923).

With the seventh-order and ninth-order WENO schemes, we have observed a significant improvement in accuracy relative to the lower order schemes, as revealed by the computed peak in the

energy dissipation variation with eddy turnover time. The ninth-order WENO scheme has shown the

best overall results when comparing with the other LES schemes (i.e., second-order central scheme

with Smagorinsky SGS model and third-order upwind biased scheme). Moreover, it was indicated

(comparing results from 643and 1283 grids) that the other LES schemes would require more than

an order of magnitude increase in mesh density to produce a similar dissipation peak to that obtained with the ninth-order WENO Scheme. The ninth-order WENO scheme, which is an ILES

scheme, has also been compared with a spectral scheme for DNS. Similar results were obtained for

the turbulent kinetic energy and turbulent dissipation. Generally, the energy spectrum of the WENO

scheme closely tracked that of the DNS spectral scheme, except at the highest frequencies, where

the WENO scheme exhibited some numerical dissipation effects.

The computational cost of the WENO scheme is high. However, current simulation results have

provided indications that for the ninth-order scheme this additional cost relative to conventional

lower order LES schemes (such as those based on central differencing with the Smagorinsky model

and third-order upwind biased differencing) as well as fifth-order WENO is compensated by improved resolution. In addition, as revealed by Shi et al. [77], the ninth-order WENO scheme has

the additional advantage of capturing small scale structures of complex flows that can include shock

waves with significantly coarser meshes than required by much lower order schemes.

Although we have observed improved accuracy with the upwind biased higher order schemes,

there is still a need to reduce the numerical dissipation, as revealed by the computed energy spectra.

Recent developments for WENO schemes such as ESWENO suggest that this can be done while

retaining appropriate shock capturing capability. In addition, these schemes could possibly be designed so that they have low dispersion errors and preserve properties such as kinetic energy. Such

schemes would then have sufficiently low numerical errors that the dominant errors in a simulation

would be due to modeling, allowing a better understanding of SGS modeling effects that can lead to

improved LES results.

28
References

1. Piomelli, U., “Large Eddy and Direct Simulation of Turbulent Flows,” von K`arm`an Institute for

Fluid Dynamics Lecture Series, Introduction to Turbulence Modeling, March 22–26, 2004.

2. Smagorinsky, J., “General Circulation Experiments with the Primitive Equations,” Monthly

Weather Review, Vol. 91, No. 3, 1963, pp. 99–164.

3. Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., “A Dynamic Subgrid-Scale EddyViscosity Model,” Physics of Fluids A, Vol. 3, No. 7, 1991, pp. 1760–1765.

4. Bardina, J., Ferziger, J. H., and Reynolds, W. C., “Improved Subgrid Scale Models for Large

Eddy Simulations,” AIAA Paper 80-1357, Jan. 1980.

5. Deardorff, J. W., “The Use of Subgrid Transport Equations in a Three-Dimensional Model of

Atmospherical Turbulence,” Journal of Fluids Engineering Transactions, Vol. 95, 1973, pp. 429–

438.

6. Adams, N. A. and Stolz, S., “Deconvolution Methods for Subgrid Scale Approximation in LES,”

Modern Simulation Strategies for Turbulent Flow, editor: B. J. Geurts, R. T. Edwards, Inc.,

Philadelphia, 2001, pp. 21–44.

7. Stolz, S. and Adams, N. A., “On the Approximate Deconvolution Procedure for LES,” Physics

of Fluids, Vol. 11, 1999, pp. 1699–1701.

8. Grinstein, F. F. and Fureby, C., “Recent Progress on Flux-Limiting Based Implicit Large Eddy

Simulation,” European Conference on Computational Fluid Dynamics (ECCOMAS CFD) 2006,

editors: P. Wesseling, E. O ̃nnate, J. P`eriaux, TU Delft, Delft The Netherlands, 2006.

9. Sagaut, P., Large Eddy Simulation for Incompressible Flows, Third Edition, Springer, Berlin,

2006.

10. Richtmyer, R. D. and Morton, K. W., Difference Methods for Initial Value Problems, John Wiley

& Sons, 1967.

11. Boris, J. P., “On Large Eddy Simulation Using Subgrid Turbulence Models,” Whither Turbulence? Turbulence at the Crossroads, Ed.: J. L. Lumley, New York: Springer-Verlag, 1990, pp.

344–353.

12. Boris, J. P., “A Fluid Transport Algorithm that Works,” Computing as a Language of Physics,

Vienna: International Atomic Energy Agency, 1971, pp. 171–189.

13. Boris, J. P. and Book, D. L. “Flux-Corrected Transport I: SHASTA – A Fluid Transport Algorithm that Works,” Journal of Computational Physics, Vol. 11, 1973, pp. 38–69.

14. Boris, J. P., Book, D. L., and Hain, K. “Flux-Corrected Transport II: Generalizations of the

Method,” Journal of Computational Physics, Vol. 18, 1975, pp. 248–283.

15. Boris, J. P. and Book, D. L. “Flux-Corrected Transport III: Minimal-Error FCT Algorithms,”

Journal of Computational Physics, Vol. 20, 1976, pp. 397–431.

16. Grinstein, F. F., Margolin, L. G., and Rider, W. J. (editors), Implicit Large Eddy Simulation,

Cambridge University Press, 2007.

17. van Leer, B. “Towards the Ultimate Conservative Difference Scheme. I. The Quest of Monotonicity,” Lecture Notes in Physics, Vol. 18, editors: E. Cabannes and R. Temam, Springer-Verlag,

1973, pp. 163–168.

29
18. van Leer, B. “Towards the Ultimate Conservative Difference Scheme II. Monotonicity and Conservation Combined in a Second-Order Scheme,” Journal of Computational Physics, Vol. 14,

1974, pp. 361–370.

19. van Leer, B. “Towards the Ultimate Conservative Difference Scheme III. Upstream-Centered

Finite-Difference Schemes for Ideal Compressible Flow,” Journal of Computational Physics, Vol.

23, 1977, pp. 263–275.

20. van Leer, B. “Towards the Ultimate Conservative Difference Scheme IV. A New Approach to

Numerical Convection,” Journal of Computational Physics, Vol. 23, 1977, pp. 276–299.

21. van Leer, B. “Towards the Ultimate Conservative Difference Scheme V. A Second-Order Sequel

to Godunov’s Method,” Journal of Computational Physics, Vol. 32, 1979, pp. 101–136.

22. Colella, P. and Woodward, P. R., “The Piecewise Parabolic Method (PPM) for Gas-Dynamical

Simulations,” Journal of Computational Physics, Vol. 54, 1984, pp. 174–201.

23. Woodward, P. R. and Colella, P., “Numerical Simulation of Two-Dimensional Fluid Flow with

Strong Shocks,” Journal of Computational Physics, Vol. 54, 1984, pp. 115–173.

24. Harten, A., “High Resolution Schemes for Hyperbolic Conservation Laws,” Journal of Computational Physics, Vol. 49, 1983, pp. 357-393.

25. Harten, A., Engquist, B., Osher, S., and Chakravarthy, S., “Uniformly High-Order Accurate

Essentially Non-oscillatory Schemes, III,” Journal of Computational Physics, Vol. 71, 1987, pp.

231–303.

26. Sweby, P. K., “High-Resolution Schemes Using Flux Limiters for Hyperbolic Conservation

Laws,” SIAM Journal of Numerical Analysis, Vol. 21, 1984, pp. 995–1011.

27. Yee, H., Warming, R. F., and Harten, A., “Implicit Total Variation Diminishing (TVD) Schemes

for Steady-State Calculations,” Journal of Computational Physics, Vol. 57, 1985, pp. 327–360.

28. Colella, P., “A Direct Eulerian MUSCL Scheme for Gas Dynamics,” SIAM Journal of Scientific

Statistical Computing, Vol. 6, 1985, pp. 104–117.

29. Colella, P. and Glaz, H. M., “Efficient Solution Algorithms for the Riemann Problems for Real

Gases,” Journal of Computational Physics, Vol. 59, 1985, pp. 264–289.

30. Lax, P. and Wendroff, B., “Systems of Conservation Laws,” Communications on Pure and Applied Mathematics, Vol. 13, 1960, pp. 217–237.

31. Godunov, S. K., “A Finite-Difference Method for the Numerical Computation and Discontinuous Solutions of the Equations of Fluid Dynamics,” Matematicheskii Sbornik, Vol. 47, 1959, pp.

271–306.

32. Harten, A., Hyman, J. M., and Lax, P. D. (with Appendix by B. Keyfitz), “On Finite-Difference

Approximations and Entropy Conditions for Shocks,” Communications on Pure and Applied

Mathematics, Vol. 29, 1976, pp. 297–322.

33. Jameson, A., “Analysis and Design of Numerical Schemes for Gas Dynamics 1: Artificial Diffusion, Upwind Biasing, Limiters and Their Effect on Accuracy and Multigrid Convergence,”

International Journal of Computational Fluid Dynamics, Vol. 4, 1995, pp. 171–218.

34. Jameson, A., “Analysis and Design of Numerical Schemes for Gas Dynamics 2: Artificial Diffusion and Discrete Shock Structure,” International Journal of Computational Fluid Dynamics,

Vol. 5, 1995, pp. 1–38.

30
35. Yee, H. C. and Sj ̈ogreen, B., “Development of Low Dissipative High Order Filter Schemes for

Multiscale Navier-Stokes/MHD Systems,” Journal of Computational Physics, Vol. 225, 2007,

pp. 910–934.

36. Shu, C.-W., “Essentially Non-oscillatory and Weighted Essentially Non-oscillatory Schemes for

Hyperbolic Conservation Laws,” Institute for Computer Applications in Science and Engineering

(ICASE) Report No. 97-65 (NASA/CR-97-206253), Nov. 1997.

37. Hill, D. J. and Pullin, D. I., “Hybrid Tuned Center-Difference-WENO Method for Large Eddy

Simulations in the Presence of Strong Shocks,” Journal of Computational Physics, Vol. 194,

2004, pp. 435–450.

38. Costa, B. and Don, W. S., “Multi-domain Hybrid Spectral-WENO Methods for Hyperbolic

Conservation Laws,” Journal of Computational Physics, Vol. 224, 2007, pp. 970–991.

39. Visbal, M. R. and Rizzeta, D. P., “Large-Eddy Simulation on Curvilinear Grids Using Compact

Differencing and Filtering Schemes,” Journal of Fluids Engineering, Vol. 124, 2002, pp. 836–

847.

40. Visbal, M. R. and Gaitonde, D. V., “Shock Capturing Using Compact-Differencing-Based Methods,” AIAA Paper 2005-1265, Jan. 2005.

41. Garnier, E., Mossi, M., Saugaut, P., Comte, P., and Deville, M., “On the Use of Shock-Capturing

Schemes for Large-Eddy Simulation,” Journal of Computational Physics, Vol. 153, 1999, pp.

273–311.

42. Thornber, B., Mosedale, A., and Drikakis, D., “On the Implicit Large Eddy Simulations of

Homogeneous Decaying Turbulence,” Journal of Computational Physics, Vol. 226, 2007, pp.

1902–1929.

43. Martin, M. P., Taylor, E. M., Wu, M., and Weirs, V. G., “A Bandwidth-Optimized WENO

Scheme for the Effective Direct Numerical Simulation of Compressible Turbulence,” Journal of

Computational Physics, Vol. 220, 2006, pp. 270–289.

44. Taylor, E. M. and Martin, M. P., “Stencil Adaptation Properties of a WENO Scheme in Direct

Numerical Simulations of Compressible Turbulence,” Journal of Scientific Computing, Vol. 30,

No. 3, 2007, pp. 533–554.

45. Ishihara, T., Gotoh, T., and Kaneda, Y., “Study of High-Reynolds Number Isotropic Turbulence

by Direct Numerical Simulation,” Annual Review of Fluid Mechanics, Vol. 41, 2009, pp. 165–

180.

46. Gaitonde, D., Shang, J., and Young, J., “Practical Aspects of Higher-Order Numerical Schemes

for Wave Propagation Phenomena,” International Journal for Numerical Methods in Engineering, Vol. 45, 1999, pp. 1849–1869.

47. Knight, D., Yan, H., and Zheltovodov, A., “Large Eddy Simulation of Supersonic Turbulent

Flow in Expansion-Compression Corner,” DNS/LES - Progress and Challenges, Proceedings of

the Third AFOSR International Conference on DNS/LES, edited by C. Liu, L. Sakell, and T.

Beutner, Greyden Press, Columbus, 2001, pp. 183–194.

48. Urbin, G. and Knight, D., “Large-Eddy Simulation of a Supersonic Boundary Layer Using an

Unstructured Grid,” AIAA Journal, Vol. 39, No. 7, 2001, pp. 1288–1295.

49. Krist, S. L., Biedron, R. T., and Rumsey, C. L., “CFL3D User’s Manual (Version 5.0),” NASA

TM-1998-208444, June 1998.

31
50. Anderson, W. K., Thomas, J. L., and van Leer, B., “Comparison of Finite Volume Flux Vector

Splittings for the Euler Equations,” AIAA Journal, Vol. 24, No. 9, 1986, pp. 1453–1460.

51. Roe, P. L., “Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes,” Journal of Computational Physics, Vol. 43, 1981, pp. 357–372.

52. Koren, B., “Upwind Schemes for the Navier-Stokes Equations,” Proceedings of the Second

International Conference on Hyperbolic Problems, Aachen. Vieweg: Braunscheweig, 1988.

53. Jiang, G. S. and Shu, C.-W., “Efficient Implementation of Weighted ENO Schemes,” Journal of

Computational Physics, Vol. 126, 1996, pp. 202–228.

54. Henrick, A. K., Aslam, T. D., and Powers, J. M., “Mapped Weighted Essentially Non-oscillatory

Schemes: Achieving Optimal Order near Critical Points,” Journal of Computational Physics, Vol.

207, 2005, pp. 542–567.

55. Borges, R., Carmona, M., Costa, B., and Don, W. S., “An Improved Weighted Essentially Nonoscillatory Scheme for Hyperbolic Conservation Laws,” Journal of Computational Physics, Vol.

227, 2008, pp. 3191–3211.

56. Balsara, D. S. and Shu, C.-W., “Monotonicity Preserving Weighted Essentially Non-oscillatory

Schemes with Increasingly High Order of Accuracy,” Journal of Computational Physics, Vol.

160, 2000, pp. 405–452.

57. Shen, Y. and Zha, G., “Improved Seventh-order WENO Scheme,” AIAA Paper 2010-1451, Jan.

2010.

58. Shen, Y. and Zha, G., “Improvement of Stability and Accuracy for Weighted Essentially

Nonoscillatory Scheme,” AIAA Journal, Vol. 47, No. 2, 2009, pp. 331–344.

59. Kamakoti, R. and Pantano, C., “High-Order Narrow Stencil Finite-Difference Approximations

of Second-Order Derivatives Involving Variable Coefficients,” SIAM Journal of Scientific Computing, Vol. 31, No. 6, 2009, pp. 4222–4243.

60. Atkins, H. L., “High-Order ENO Methods for the Unsteady Compressible Navier-Stokes Equations,” AIAA Paper 1991-1557, June 1991.

61. Balakumar, P., Zhao, H., and Atkins, H., “Stability of Hypersonic Boundary Layers over a

Compression Corner,” AIAA Paper 2002-2848, June 2002.

62. Compte-Bellot, G. and Corrsin, S., “Simple Eulerian Time Correlation of Full- and NarrowBand Velocity Signals in Grid-Generated, Isotropic Turbulence,” Journal of Fluid Mechanics,

Vol. 48, Part 2, 1971, pp. 273–337.

63. Clark, R. A., Ferziger, J. H., and Reynolds, W. C., “Evaluation of Subgrid-Scale Models Using

an Accurately Simulated Turbulent Flow,” Journal of Fluid Mechanics, Vol. 91, Part 1, 1979, pp.

1–16.

64. Vreman, A. W., Geurts, B. J., Kuerten, J. G. M., and Zandbergen, P. J., “A Finite Volume

Approach to Large Eddy Simulation of Compressible, Homogeneous, Isotropic, Decaying Turbulence,” International Journal for Numerical Methods in Fluids, Vol. 15, 1992, pp. 799-816.

65. Mansour, N. N. and Wray, A. A., “Decay of Isotropic Turbulence at Low Reynolds Number,”

Physics of Fluids, Vol. 6, No. 2, 1994, pp. 808–814.

66. Spyropoulos, E. T. and Blaisdell, G. A., “Evaluation of the Dynamic Model for Simulations of

Compressible Decaying Isotropic Turbulence,” AIAA Journal, Vol. 34, No. 5, 1996, pp. 990–998.

32
67. Seror, C., Sagaut, P., Bailly, C., and Juve, D., “Subgrid-Scale Contribution to Noise Production

in Decaying Isotropic Turbulence,” AIAA Journal, Vol. 38, No. 10, 2000, pp. 1795–1803.

68. Hughes, T. J. R., Mazzei, L., and Oberai, A. A., “The Multiscale Formulation of Large Eddy

Simulation: Decay of Homogeneous Isotropic Turbulence,” Physics of Fluids, Vol. 13, No. 2,

2001, pp. 505–512.

69. Hansen, A., Sorensen, N. N., Johansen, J., and Michelsen, J. A., “Detached-Eddy Simulation of

Decaying Homogeneous Isotropic Turbulence,” AIAA Paper 2005-885, Jan. 2005.

70. Zang, T. A., Dahlburg, R. B., and Dahlburg, J. P., “Direct and Large-Eddy Simulations of ThreeDimensional Compressible Navier-Stokes Turbulence,” Physics of Fluids A, Vol. 4, No. 1, 1992,

pp. 127–140.

71. Erlebacher, G., Hussaini, M. Y., Speziale, C. G., and Zang, T. A., “Toward the Large-eddy

Simulation of Compressible Turbulent Flows,” Journal of Fluid Mechanics, Vol. 238, 1992, pp.

155–185.

72. Erlebacher, G., Hussaini, M. Y., Kreiss, H. O., and Sarkar, S., “The Analysis and Simulation of

Compressible Turbulence,” Theoretical Computational Fluid Dynamics, Vol. 2, 1990, pp. 73–95.

73. Passot, T. and Pouquet, A., “Numerical Simulation of Compressible Homogeneous Flows in the

Turbulent Regime,” Journal of Fluid Mechanics, Vol. 181, 1987, pp. 441–466.

74. Lesieur, M., Turbulence in Fluids, 2nd ed., Kluwer Academic, Norwell, MA, 1990.

75. Piomelli, U., Radhakrishnan, S., and de Prisco, G., “Turbulent Eddies in the RANS/LES Transition Region,” in Advances in Hybrid RANS-LES Modeling, editors: S.-H. Peng and W. Haase,

Springer, Berin, 2008, pp. 21–36.

76. Slomski, J. F. and Chang, P. A., “Synthesis and Advection of Isotropic Turbulence to a Target

Intensity and Spectrum,” AIAA Paper 2008-547, Jan. 2008.

77. Shi, J., Zhang, Y.-T., and Shu, C.-W., “Resolution of High Order WENO schemes for complicated flow structures,” Journal of Computational Physics, Vol. 186, 2003, pp. 690–696.

78. Yamaleev, N. K. and Carpenter, M. H., “High-Order Energy Stable WENO Schemes,” AIAA

Paper 2009-1135, Jan. 2009.

79. Yamaleev, N. K. and Carpenter, M. H., “Third-Order Energy Stable WENO Scheme,” Journal

of Computational Physics, Vol. 228, No. 8, 2009, pp. 3025–3047.

80. Yamaleev, N. K. and Carpenter, M. H., “A Systematic Methodology for Constructing HighOrder Energy Stable WENO Schemes,” Journal of Computational Physics, Vol. 228, 2009, pp.

4248–4272.

81. Kravchenko, A. G. and Moin, P., “On the Effect of Numerical Errors in Large Eddy Simulations

of Turbulent Flows,” Journal of Computational Physics, Vol. 131, 1997, pp. 310–322.

82. Lilly, D. K., “On the Computational Stability of Numerical Solutions of Time-Dependent NonLinear Geophysical Fluid Dynamics Problems,” Monthly Weather Review, Vol. 93, No. 1, 1965,

pp. 11–26.

83. Ducros, F., Laporte, F., Soule`res, T., Guinot, V., Moinat, P., and Caruelle, B., “High-Order

Fluxes for Conservative Skew-Symmetric-Like Schemes in Structured Meshes: Application to

Compressible Flows,” Journal of Computational Physics, Vol. 161, 2000, pp. 114–139.

33
84. Verstappen, R. W. C. P. and Veldman, A. E. P., “Symmetry-Preserving Discretization of Turbulent Flow,” Journal of Computational Physics, Vol. 187, 2003, pp. 343–368.

85. Feiereisen, W. J., Reynolds, W. C., and Ferziger, J. H., “Numerical Simulation of a Compressible Homogeneous, Turbulent Shear Flow,” Rep. TF-13, Thermosci. Div., Department of Mech.

Eng., Stanford University, 1981.

86. Kok, J. C., “A High-Order Low-Dispersion Symmetry-Preserving Finite-Volume Method for

Compressible Flow on Curvilinear Grids,” Journal of Computational Physics, Vol. 228, 2009,

pp. 6811–6832.

87. Ghosal, S., “An Analysis of Numerical Errors in Large-Eddy Simulations of Turbulence,” Journal of Computational Physics, Vol. 125, 1996, pp. 187–206.

88. Tam, C. K. W. and Webb, J. C., “Dispersion-Relation-Preserving Finite Difference Schemes for

Computational Acoustics,” Journal of Computational Physics, Vol. 107, 1993, pp. 262–281.

89. Batchelor, G. K., The Theory of Homogeneous Turbulence, Cambridge University Press, Cambridge, 1953.

90. Tsinober, A., An Informal Introduction to Turbulence, Kluwer Academic Publishers, Dordrecht,

2001.

91. Reid, W. H. and Harris, D. L., “Similarity Spectra in Isotropic Turbulence,” Physics of Fluids,

Vol. 2, 1959, pp. 139–146.

92. Kraichnan, R. H., “The Structure of Isotropic Turbulence at Very High Reynolds Numbers,”

Journal of Fluid Mechanics, Vol. 5, 1959, pp. 497–543.

93. Smith, L. M. and Reynolds, W. C., “The Dissipation Range Spectrum and the Velocity Derivative Skewness in Turbulent Flows,” Physics of Fluids A, Vol. 3, 1991, pp. 992-994.

34
xi+1/2 ∆x xi-1/2

xi xi+1 x xi-2 xi-1 i+2

Figure 1. Uniform computational grid for WENO scheme. f

i+1/2

S1

xi xi+1 x xi-2 xi-1 i+2

xi+1/2 S0

S2

Figure 2. The five point WENO stencil constructed from three-point stencils S0, S1, and S2.

35
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1923, upwind, ILES

1283, upwind, ILES

643, upwind, ILES

323, upwind, ILES

(a) (b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

1923, upwind, ILES

1283, upwind, ILES

643, upwind, ILES

323, upwind, ILES

(c)

Figure 3. Time history of various statistics for upwind ILES for various grid sizes; (a) normalized

kinetic energy, (b) normalized dissipation, (c) skewness.

36
Figure 4. Isocontours of nondimensional vorticity for upwind ILES for various grid sizes, at t

∗ = 5.

37
Figure 5. Isocontours of nondimensional Q-invariant for upwind ILES for various grid sizes, at

t

∗ = 5.

38
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1923, central, cs=0

1283, central, cs=0

643, central, cs=0

323, central, cs=0

(a) (b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

1923, central, cs=0

1283, central, cs=0

643, central, cs=0

323, central, cs=0

(c)

Figure 6. Time history of various statistics for central (cs = 0) for various grid sizes; (a) normalized

kinetic energy, (b) normalized dissipation, (c) skewness.

39
Figure 7. Isocontours of nondimensional Q-invariant for central (cs = 0) for various grid sizes, at

t

∗ = 5.

40
Z Z Z

Z

Z

Z

Z

Z Z

Z

Z

ZZ

ZZ

ZZ

ZZZZZZ

ZZ

ZZ

ZZZZZZZZZZZZZZZZZZZZZZZZZZ

ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ

κ

E(

κ

)

100 101 102 103 10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

κ

-5/3

1923, U

1283, U

643, C

643, U

323, U

323, C

1283, C

1923, C

Figure 8. Energy spectrum for upwind ILES and central (cs = 0), at t

∗ = 5; U indicates upwind, C

indicates central.

41
κ

E(

κ

)

100 101 102 103 10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

t*=1.5

t*=5

t*=10

κ

-5/3

Figure 9. Energy spectrum at three different times for central (cs = 0), 1283 grid.

42
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

323, central, cs = 0

323, central, cs = 0.1

323, central, cs = 0.2

323, central, cs = 0.4

(a) (b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

323, central, cs = 0

323, central, cs = 0.1

323, central, cs = 0.2

323, central, cs = 0.4

(c)

Figure 10. Time history of various statistics for central differencing with the Smagorinsky model on

323for various cs values; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness.

43
Figure 11. Energy spectrum for the Smagorinsky model with different values of cs, 323

, central

differencing, at t

∗ = 5.

44
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1283, central, cs = 0

1283, central, cs = 0.2

643, central, cs = 0

643, central, cs = 0.2

(a) (b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

1283, central, cs = 0

1283, central, cs = 0.2

643, central, cs = 0

643, central, cs = 0.2

(c)

Figure 12. Time history of various statistics for central differencing with the Smagorinsky model

with different values of cs on 643and 1283; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness.

45
Figure 13. Energy spectrum for the Smagorinsky model with different values of cs on two finer grid

sizes, central differencing, at t

∗ = 5.

46
Figure 14. Isocontours of nondimensional Q-invariant for central differencing (Smagorinsky model

with cs = 0.2) for various grid sizes, at t

∗ = 5.

47
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1923, 5th order

1283, 5th order

643, 5th order

323, 5th order

(a)

t

*

ε

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

1923, 5th order

1283, 5th order

643, 5th order

323, 5th order

(b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

1923, 5th order

1283, 5th order

643, 5th order

323, 5th order

(c)

Figure 15. Time history of various statistics on different grid sizes using 5th order WENO-L ILES

scheme; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness.

48
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1923, 7th order

1283, 7th order

643, 7th order

323, 7th order

(a)

t

*

ε

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

1923, 7th order

1283, 7th order

643, 7th order

323, 7th order

(b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

1923, 7th order

1283, 7th order

643, 7th order

323, 7th order

(c)

Figure 16. Time history of various statistics on different grid sizes using 7th order WENO-L ILES

scheme; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness.

49
t

*

k

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1923, 9th order

1283, 9th order

643, 9th order

323, 9th order

(a)

t

*

ε

0 2 4 6 8 10 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

1923, 9th order

1283, 9th order

643, 9th order

323, 9th order

(b)

t

*

S

k

0 2 4 6 8 10 -1.5

-1

-0.5

0

0.5

1

1923, 9th order

1283, 9th order

643, 9th order

323, 9th order

(c)

Figure 17. Time history of various statistics on different grid sizes using 9th order WENO-L ILES

scheme; (a) normalized kinetic energy, (b) normalized dissipation, (c) skewness.

50
N

-2/3

normalized peak

ε

0 0.0002 0.0004 0.0006 0.0008 0.001 1

1.2

1.4

1.6

1.8

2

2.2

central-2nd

upwind-3rd

weno-5th

weno-7th

weno-9th

643 grid

323 grid

1283 grid

1923 grid

Figure 18. Normalized peak ε as a function of grid spacing squared.

51
(a) (b)

(c) (d)

Figure 19. Energy spectrum for 5th-order, 7th-order, and 9th-order WENO-L ILES on four different

grids at t

∗ = 5; (a) 323

, (b) 643, (c) 1283, (d) 1923.

52
κ

E(

κ

)

100 101 102 103 10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

1923, 9th order

1283, 9th order

643, 9th order

323, 9th order

κ

-5/3

Figure 20. Energy spectrum for 9th-order WENO-L ILES at t

∗ = 5

.

53
κ

E(

κ

)

100 101 102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

κ

-5/3

1283, C2

1283, W9

1283, U3

1283, W5

Figure 21. Energy spectrum for C2 = central 2nd order, cs = 0 (CFL3D), U3 = upwind 3rd order

ILES (CFL3D), W9 = WENO-L 9th order, and W5 = WENO-L 5th order ILES on same 1283 grid,

at t

∗ = 5.

54
Figure 22. Isocontours of nondimensional Q-invariant for 9th-order WENO-L for various grid sizes,

at t

∗ = 5.

55
t

*

ε

0 0.5 1 1.5 2 2.5 3

1

1.2

1.4

1.6

1.8

2

WENO-L

WENO, α1, εw = 10-6

WENO, α1, εw = 10-2

WENO, α2, εw = 10-6

WENO, α2, εw = 10-2

1283, 5th order (a)

κ

E(

κ

)

20 40 60 80 100

10-10

10-9

10-8

10-7

10-6

10-5

10-4

10-3

WENO-L

WENO, α1, εw = 1 x 10-6

WENO, α1, εw = 1 x 10-2

WENO, α2, εw = 1 x 10-6

WENO, α2, εw = 1 x 10-2

1283, 5th order

(b)

t

*

ε

0 0.5 1 1.5 2 2.5 3

1

1.2

1.4

1.6

1.8

2

WENO-L

WENO, α1, εw = 10-6

WENO, α1, εw = 10-2

WENO, α2, εw = 10-6

WENO, α2, εw = 10-2

1283, 7th order

(c)

κ

E(

κ

)

20 40 60 80 100

10-10

10-9

10-8

10-7

10-6

10-5

10-4

10-3

WENO-L

WENO, α1, εw = 1 x 10-6

WENO, α1, εw = 1 x 10-2

WENO, α2, εw = 1 x 10-6

WENO, α2, εw = 1 x 10-2

1283, 7th order

(d)

Figure 23. Effect of two different ways of computing nonlinear WENO weights on the normalized

dissipation and energy spectrum (1283); (a) dissipation, 5th order, (b) spectrum at t

∗ = 2, 5th order

(c) dissipation, 7th order (d) spectrum at t

∗ = 2, 7th order.

56
t

*

k

0 1 2 3

0.6

0.8

1

1.2

643

643, εw = 0.1

643, εw = 0.01

643, εw = 0.001

(a)

t

*

ε

0 1 2 3

0.8

1

1.2

1.4

1.6

643

643, εw = 0.1

643, εw = 0.01

643, εw = 0.001

(b)

κ

E(

κ

)

100 101 102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

643

643, εw = 0.1

643, εw = 0.01

643, εw = 0.001

κ

-5/3

(c)

Figure 24. Time history of statistics and energy spectrum comparison for 5th-order WENO (643)

with εw = 0.1, εw = 0.01, εw = 0.001; (a) normalized kinetic energy, (b) normalized dissipation,

(c) energy spectrum at t

∗ = 1.5.

57
t

*

ε

0 1 2 3

0.8

1

1.2

1.4

1.6

1.8

2

2.2

1923, WENO-L

1923, WENO

1283, WENO-L

1283, WENO

643, WENO-L

643, WENO

323, WENO-L

323, WENO

5th order

(a)

t

*

ε

0 1 2 3

0.8

1

1.2

1.4

1.6

1.8

2

2.2

1923, WENO-L

1923, WENO

1283, WENO-L

1283, WENO

643, WENO-L

643, WENO

323, WENO-L

323, WENO

7th order

(b)

t

*

ε

0 1 2 3

0.8

1

1.2

1.4

1.6

1.8

2

2.2

1923, WENO-L

1923, WENO

1283, WENO-L

1283, WENO

643, WENO-L

643, WENO

323, WENO-L

323, WENO

9th order

(c)

Figure 25. Comparison of normalized dissipation computed with WENO-L and WENO; (a) 5th

order, b) 7th order, (c) 9th order.

58
κ

E(

κ

)

100 101 102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

1923, WENO-L

1923, WENO

1283, WENO-L

1283, WENO

643, WENO-L

643, WENO

323, WENO-L

323, WENO

κ

-5/3

5th order

(a)

κ

E(

κ

)

100 101 102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

1923, WENO-L

1923, WENO

1283, WENO-L

1283, WENO

643, WENO-L

643, WENO

323, WENO-L

323, WENO

κ

-5/3

7th order

(b)

κ

E(

κ

)

100 101 102 103

10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

102

104

1923, WENO-L

1923, WENO

1283, WENO-L

1283, WENO

643, WENO-L

643, WENO

323, WENO-L

323, WENO

κ

-5/3

9th order

(c)

Figure 26. Comparison of energy spectrum computed with WENO-L and WENO, t

∗ = 2; (a) 5th

order, b) 7th order, (c) 9th order.

59
t

*

k

0 1 2 3 4 5

0.4

0.6

0.8

1

1.2 643

, WENO, 9th order

1283, WENO, 9th order

643, central, 2nd order, cs=0.2

1283, central, 2nd order, cs=0.2

643, upwind, 3rd order

1283, upwind, 3rd order

(a) (b)

(c)

Figure 27. Time history of statistics and energy spectrum comparison for 9th-order WENO, 2ndorder central with Smagorinsky SGS model (cs = 0.2), and 3rd-order upwind biased; (a) normalized

kinetic energy, (b) normalized dissipation, (c) energy spectrum at t

∗ = 5.

60
t

*

k

0 1 2 3

0.4

0.6

0.8

1

1.2

1.4

1923, DNS, Spectral

1923, WENO, 9th order

1923, central, 2nd order

(a)

t

*

ε

0 1 2 3

0.8

1

1.2

1.4

1.6

1.8

2

2.2

2.4 1923

, DNS, Spectral

1923, WENO, 9th order

1923, central, 2nd order

(b)

κ

E(

κ

)

100 101 102

10-10

10-8

10-6

10-4

10-2

100

102

104 1923

, DNS, Spectral

1923, WENO, 9th order

1923, central, 2nd order

κ

-5/3

(c)

Figure 28. Time history of statistics and energy spectrum comparison for 9th-order WENO, 2ndorder central (cs = 0), and spectral DNS; (a) normalized kinetic energy, (b) normalized dissipation,

(c) energy spectrum at t

∗ = 2.

61
Appendix A

For the seventh-order WENO scheme, the fourth-order fluxes at xi+1/2 are as follows:

ˆf

0

i+1/2 =

1

12

(−3fi−3 + 13fi−2 − 23fi−1 + 25fi),

ˆf

1

i+1/2 =

1

12

(fi−2 − 5fi−1 + 13fi + 3fi+1),

ˆf

2

i+1/2 =

1

12

(−fi−1 + 7fi + 7fi+1 − fi+2),

ˆf

3

i+1/2 =

1

12

(3fi + 13fi+1 − 5fi+2 + fi+3),

The corresponding ideal weights are as follows:

ω0 = 1/35, ω1 = 12/35, ω2 = 18/35, ω3 = 4/35. (A1)

Then, substituting the ˆf

k

i+1/2

and the ideal weights into Eq. (29), we obtain the seventh-order flux at

xi+1/2, which is given by

ˆfi+1/2 =

1

420

(−3fi−3 + 25fi−2 − 101fi−1 + 319fi + 214fi+1 − 38fi+2 + 4fi+3). (A2)

Using Eq. (38) and the coefficients of the third-degree interpolating polynomials, the smoothness

indicators for the nonlinear weights of the seventh-order scheme are found to be

β0 = fi−3(547fi−3 − 3882fi−2 + 4642fi−1 − 1854fi)

− fi−2(7042fi−2 − 17246fi−1 + 7042fi)

+ fi−1(11003fi−1 − 9042fi) + 2107f

2

i

β1 = fi−2(267fi−2 − 1642fi−1 + 1602fi − 494fi+1)

− fi−1(2843fi−1 − 5966fi + 1922fi+1)

+ fi(3443fi − 2522fi+1) + 547f

2

i+1

(A3)

β2 = fi−1(547fi−1 − 2522fi + 1922fi+1 − 494fi+2)

− fi(3443fi − 5966fi+1 + 1602fi+2)

+ fi+1(2843fi+1 − 1642fi+2) + 267f

2

i+2

β3 = fi(2107fi − 9402fi+1 + 7042fi+2 − 1854fi+3)

− fi+1(11003fi+1 − 17246fi+2 + 4642fi+3)

+ fi+2(7043fi+2 − 3882fi+3) + 547f

2

i+3

and βk = βk/240 for k = 0, . . . , 3.

For the ninth-order WENO scheme, the fifth-order fluxes at xi+1/2 are as follows:

62
ˆf

0

i+1/2 =

1

60

(12fi−4 − 63fi−3 + 137fi−2 − 163fi−1 + 137fi),

ˆf

1

i+1/2 =

1

60

(−3fi−3 + 17fi−2 − 43fi−1 + 77fi + 12fi+1),

ˆf

2

i+1/2 =

1

60

(2fi−2 − 13fi−1 + 47fi + 27fi+1 − 3fi+2), (A4)

ˆf

3

i+1/2 =

1

60

(−3fi−1 + 27fi + 47fi+1 − 13fi+2 + 2fi+3),

ˆf

4

i+1/2 =

1

60

(12fi + 77fi+1 − 43fi+2 + 17fi+3 − 3fi+4),

The corresponding ideal weights are as follows:

ω0 = 1/126, ω1 = 10/63, ω2 = 10/21, ω3 = 20/63, ω4 = 5/126. (A5)

Then, substituting the ˆf

k

i+1/2

and the ideal weights into Eq. (29), we obtain the ninth-order flux at

xi+1/2, which is given by

ˆfi+1/2 =

1

7560

[12fi−4 − 123fi−3 + 597fi−2 − 1923fi−1 + 5637fi (A6)

+ 4125fi+1 − 915fi+2 + 165fi+3 − fi+4]

Using Eq. (38) and the coefficients of the fourth-degree interpolating polynomials, the smoothness

indicators for the nonlinear weights of the ninth-order scheme are found to be

β0 = fi−4(22658fi−4 − 20850fi−3 + 364863fi−2 − 288007fi−1 + 86329fi)

− fi−3(482963fi−3 − 1704396fi−2 + 1358458fi−1 − 411487fi)

+ fi−2(1521393fi−2 − 2462076fi−1 + 758823fi)

+ fi−1(1020563fi−1 − 649501fi) + 107918f

2

i

β1 = fi−3(6908fi−3 − 60871fi−2 + 99213fi−1 − 70237fi + 18079fi+1)

− fi−2(138563fi−2 − 464976fi−1 + 337018fi − 88297fi+1)

+ fi−1(406293fi−1 − 611976fi + 165153fi+1)

+ fi(242723fi − 140251fi+1) + 22658f

2

i+1

β2 = fi−2(6908fi−2 − 51001fi−1 + 67923fi − 38947fi+1 + 8209fi+2) (A7)

− fi−1(104963fi−1 − 299076fi + 179098fi+1 − 38947fi+2)

+ fi(231153fi − 299076fi+1 + 67923fi+2)

+ fi+1(104963fi+1 − 51001fi+2) + 6908f

2

i+2

β3 = fi−1(22658fi−1 − 140251fi + 165153fi+1 − 88297fi+2 + 18079fi+3)

− fi(242723fi − 611976fi+1 + 337018fi+2 − 70237fi+3)

+ fi+1(406293fi+1 − 464976fi+2 + 99213fi+3)

+ fi+2(138563fi+2 − 60871fi+3) + 6908f

2

i+3

β4 = fi(107918fi − 649501fi+1 + 758823fi+2 − 411487fi+3 + 86329fi+4)

− fi+1(1020563fi+1 − 2462076fi+2 + 1358458fi+3 − 233007fi+4)

+ fi+2(1521393fi+2 − 1704396fi+3 + 364863fi+4)

+ fi+3(482963fi+3 − 208501fi+4) + 22658f

2

i+4

and βk = βk/5040 for k = 0, . . . , 4.

63
Table A1. The constants crj for 8th order.

r j = 0 j = 1 j = 2 j = 3 j = 4 j = 5 j = 6 j = 7

-1 761/280 −1479/280 2441/280 −8357/840 6343/840 −613/168 171/168 −1/8

0 1/8 481/280 −499/280 481/280 −1007/840 463/840 −25/168 1/56

1 −1/56 15/56 341/280 −219/280 131/280 −167/840 43/840 −1/168

2 1/168 −11/168 73/168 743/840 −307/840 113/840 −9/280 1/280

3 −1/280 29/840 −139/840 533/840 533/840 −139/840 29/840 −1/280

4 1/280 −9/280 113/840 −307/840 743/840 73/168 −11/168 1/168

5 −1/168 43/840 −167/840 131/280 −219/280 341/280 15/56 −1/56

6 1/56 −25/168 463/840 −1007/840 481/280 −499/280 481/280 1/8

7 −1/8 171/168 −613/168 6343/840 −8357/840 2441/280 −1479/280 761/280

This eighth-order table is an amendment to Table 2.1 of Shu [36]. As in Table 2.1, the crj are

weights of interpolating polynomials; k and r denote the order of approximation and the number of

grid points to the left of point xi

, respectively.

64
Appendix B

Measurements generally give values for the velocity derivative skewness of Sk ≈ −0.5 (note that

sometimes the skewness is defined with opposite sign, giving Sk ≈ +0.5). Whether or not the

skewness depends on Reynolds number is still debated. Velocity derivative skewness is related to

vortex stretching, the basic mechanism of production of small scales in turbulence is as follows: The

spectral evolution equation for homogeneous isotropic turbulence [89] is

∂E(κ, t)/∂t = P(κ, t) − T(κ, t) − D(κ, t) (B1)

where E(κ, t) is the energy spectrum, P(κ, t) is a production spectrum, T(κ, t) is the nonlinear

transfer, and D(κ, t) = 2νκ2E(κ, t) is the dissipation spectrum. The second moment of this equation is the enstrophy balance [89]

Z ∞

0

dκ κ2∂E(κ, t)/∂t =

Z ∞

0

dκ κ2[P(κ, t) − T(κ, t) − D(κ, t)]

If, as is usual, production is concentrated at large scales, then

Z ∞

0

dκ κ2P(κ, t) ≈ 0

Writing the energy transfer as T(κ, t) = ∂F(κ, t)/∂κ, where F is the energy flux,

−

Z ∞

0

dκ κ2T(κ) = Z ∞

0

dκ κF(κ)

so that −

Z ∞

0

dκ κ2T(κ) ≥ 0 if the energy flux is from large to small scales. Thus, in the approximate enstrophy balance

Z ∞

0

dκ κ2∂E(κ, t)/∂t = −

Z ∞

0

dκ [κ

2T(κ, t) + κ2D(κ, t)] (B2)

The first term on the right side is a source term that generates enstrophy. Elementary considerations

relate this term to vortex stretching [89]. B1 The Kolmogorov theory asserts that the small scales

in any turbulent flow are approximately isotropic and homogeneous, thus, this special analysis for

homogeneous isotropic turbulence is believed to apply in much greater generality.

One of the first results of Taylor’s statistical theory of turbulence [91] is

Sk =

3

14

√

30

Z ∞

0

dκ κ2T(κ)

Z ∞

0

dκ κ2E(κ)

3/2

(B3)

Then, in view of Eq. (B2), velocity derivative skewness Sk is a dimensionless measure of the production of enstrophy, hence of the generation of small scales. The skewness also measures the strength

of the nonlinearity: For example, it vanishes in a Gaussian random field. The Batchelor skewness

relation

−

Z ∞

0

dκ κ2T(κ, t) = Z ∞

0

dκ 2νκ4E(κ, t)

B1That vortex stretching is in fact positive was demonstrated experimentally by Taylor, who doubted von Karm ́ an’s sug-  ́

gestion that vortex stretching vanishes [90].

65
permits an evaluation of the skewness in terms of integrals over the energy spectrum alone: since

Z ∞

0

dκ κ2E(κ) = 

ν

,

Sk = −

3

14

√

30

Z ∞

0

dκ 2νκ4E(κ)

Z ∞

0

dκ k2E(κ)

3/2

= −

3

14

√

30 ν



3/2

Z ∞

0

dκ 2νκ4E(κ) (B4)

Explicit evaluation requires knowing the exact energy spectrum, but the trends can be assessed

by assuming a model spectrum of the form

E(κ) = CK





2/3κ−5/3 κ ≤ κd



2/3κ

−5/3

d

f(κ/κd) κ ≥ κd

(B5)

where κd = α(/ν3)

1/4

is proportional to the Kolmogorov microscale and CK ≈ 1.5 is the Kolmogorov constant. Experimental and numerical data support the dissipation range form

f(κ/κd) = 

κ

κd

b

exp[a(1 − κ/κd)]

Direct interaction approximation (DIA) predicts this form with b = 2 and a ≈ 4. [92]. Assuming

the general expression Eq. (B5),

Sk = −

3

14

√

30α

10/3CK



3

5

+ 2e

a

Z ∞

1

dx x4+be

−ax

(B6)

Ignoring the dissipation range contribution would result in Sk ≈ −α

10/3

.

The most important conclusion is that the result is dominated by the small scales because of the

dependence on the parameter α, which determines the scale at which dissipation becomes significant, and the spectrum in the dissipation range.B2 However, the value of κd itself does not enter

Eq. (B6). In this analysis, then, the skewness is independent of Reynolds number; while it remains

uncertain whether this conclusion is supported by experimental data, we note that Reynolds number

dependence would require a departure from the Kolmogorov spectrum.

The value of the velocity derivative skewness is often taken as a figure of merit to validate numerical schemes; the connection of the skewness to vortex stretching and generation of small scales

certainly makes this reasonable. The previous analysis is easily applied to the case of the “numerical

turbulence” generated by a finite-difference scheme on a lattice with uniform isotropic spacing h.

Following Vreman et al. [64], such numerical turbulence is characterized by the multipliers lj that

express the action of the finite-difference operator Dj on “plane wave” exponential functions,

Dj e

iκ·x = ilj eiκ·x

Vreman et al. gives several examples. We cite the simplest one: the “vertex centered” operator for

which

lj =

sin(κjh)

h

The previous analysis can be repeated essentially verbatim, replacing wavevector κ and its corresponding wavenumber magnitude κ = |κ| by l and l = |l|. The result is a type of regularization

of the Navier-Stokes equation in which the lattice spacing h enters as a new parameter. The spectral evolution equation remains valid, but the constant flux Kolmogorov spectrum will be modified

B2Compare Smith and Reynolds [93], although their conclusion that the dissipation range should depend on e−κ

2

instead

of e−κ is not supported by data or theory.

66
by the regularization for scales k near the lattice spacing h: how near depends on the difference

scheme. The value of the skewness will depend on the behavior of the numerical spectrum, and

a value close to Sk ≈ −0.5 must not be expected unless the Kolmogorov spectrum can extend to

scales significantly larger than h. On the other hand, the sensitive dependence of the skewness to

small scale details justifies caution in making positive or negative claims for difference schemes on

this basis alone.

67
Appendix C

Inertial Subrange

A turbulent flow field contains (vortical) structures, which are generally called eddies. There are

multiple length scales associated with the eddies of various sizes, starting with the largest eddies

having a length scale proportional to the characteristic length of the field and terminating with the

small eddies with a length scale proportional to the Kolmogorov length scale. Nonlinear interactions

among eddies with a comparably large length scale (i.e., small wavenumber κ) are approximately

inviscid in nature and can result in eddies with smaller length scales (i.e., higher wavenumbers).

One may view this process as large scale structures breaking up into smaller scale structures.

There is a fluctuating velocity associated with each eddy; and thus, there is a certain kinetic

energy corresponding to each eddy. Each fluctuating velocity has a frequency. The largest eddies

correspond to low frequencies, and the smallest correspond to high frequencies. Consider a turbulent

flow field that is statistically homogeneous with respect to time. Then, the turbulent kinetic energy

of an eddy can be represented by a full range of wavenumbers. An ideal energy spectrum for

homogeneous isotropic turbulence is displayed in Fig. C1. The spectrum shows that during the

initial increase of the wavenumber κ the energy increases, due to the largest eddies in the domain,

to a peak level, where the eddies have an integral length scale that is an O(1) fraction of the large

scales. Interactions between the eddies with an integral scale start an energy cascade process. The

energy decreases as the frequencies of the spectrum increases. It decreases with a slope of −5/3, and

this part of the energy spectrum which is characterized by the Taylor microscale is often called the

Kolmogorov inertial subrange. To have a well defined inertial subrange there must be a sufficiently

large Taylor microscale Reynolds number. This regime of the energy spectrum is followed by a

region where the smallest eddies are dissipated.

It is difficult to establish an inertial subrange with κ

−5/3 behavior in decaying isotropic turbulence. As indicated by Ishihara [45], a Taylor microscale Reynolds number Reλ ≥ 200 is required

to have an inertial subrange. In order to produce a significant inertial subrange and to investigate the

physics of the small scales of isotropic turbulence, an alternative simulation approach is frequently

considered. With this approach, a forcing function is introduced to emulate the energy transfer to

the inertial subrange from the larger energy containing eddies. Thus, in a simulation the large scales

are modeled and the small scales are resolved, and higher Reλ can be realized, which can result in

a sufficiently wide inertial subrange. Due to the steady forcing the turbulence becomes statistically

steady.

Since the Taylor microscale Reynolds number for the case being considered is relatively low,

and there is no forcing function, we should observe at most a very small inertial subrange. Nevertheless, as a reference, we include the inertial subrange curve which has a -5/3 slope. This facilitates

comparing with results from other investigators who include this curve.

68
κT κ log(κ) I κk

log(E(

κ))

integral scale

large scale

inertial

subrange

dissipation

scale

Figure C1. Sketch of Kolmogorov energy spectrum indicating length scales corresponding to different parts of the spectrum.

69
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

Parametric Study of Decay of Homogeneous Isotropic Turbulence Using 

Large Eddy Simulation

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Swanson, R. C.; Rumsey, Christopher L.; Rubinstein, Robert; Balakumar, 

Ponnampalam; Zang, Thomas A.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20152

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 02

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

Numerical simulations of decaying homogeneous isotropic turbulence are performed with both low-order and high-order spatial discretization schemes.

The turbulent Mach and Reynolds numbers for the simulations are 0.2 and 250, respectively. For the low-order schemes we use either second-order central or 

third-order upwind biased differencing. For higher order approximations we apply weighted essentially non-oscillatory (WENO) schemes, both with linear and 

nonlinear weights. There are two objectives in this preliminary effort to investigate possible schemes for large eddy simulation (LES). One

is to explore the capability of a widely used low order computational fluid dynamics (CFD) code to perform LES computations. The other is to determine

the effect of higher order accuracy (fifth, seventh, and ninth order) achieved with high-order upwind biased WENO-based schemes. Turbulence statistics, such 

as kinetic energy, dissipation, and skewness, along with the energy spectra from simulations of the decaying turbulence problem are used to assess and

compare the various numerical schemes. In addition, results from the best performing schemes are compared with those from a spectral scheme. The effects of 

grid density, ranging from 32 cubed to 192 cubed, on the computations are also examined. The fifth-order WENO-based scheme is found to be too dissipative, 

especially on the coarser grids. However, with the seventh-order and ninth-order WENO-based schemes we observe a significant improvement in accuracy 

relative to the lower order LES schemes, as revealed by the computed peak in the energy dissipation and by the energy spectrum.

15. SUBJECT TERMS

Isotropic decay; Large eddy simulation; Numerical methods; Turbulence

18. NUMBER

 OF 

 PAGES

74

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

561581.02.08.07.46.04 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2012-217593

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 08 - 2012