James R. DeBonis

Glenn Research Center, Cleveland, Ohio

A High-Resolution Capability for Large-Eddy

Simulation of Jet Flows

NASA/TM—2011-216833

March 2011

AIAA–2010–5023
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c 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 Offi cer. 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 Reports 

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 signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings 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. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

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 via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
James R. DeBonis

Glenn Research Center, Cleveland, Ohio

A High-Resolution Capability for Large-Eddy

Simulation of Jet Flows

NASA/TM—2011-216833

March 2011

AIAA–2010–5023

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

40th Fluid Dynamics Conference and Exhibit

sponsored by the American Institute of Aeronautics and Astronautics

Chicago, Illinois, June 28 to July 1, 2010
Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

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

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 

This report is a formal draft or working 

paper, intended to solicit comments and 

ideas from a technical peer group.

This report contains preliminary fi ndings, 

subject to revision as analysis proceeds.
A large-eddy simulation (LES) code that utilizes high-resolution numerical schemes is

described and applied to a compressible jet flow. The code is written in a general manner

such that the accuracy/resolution of the simulation can be selected by the user. Time discretization is performed using a family of low-dispersion Runge-Kutta schemes, selectable

from first- to fourth-order. Spatial discretization is performed using central differencing

schemes. Both standard schemes, second- to twelfth-order (3 to 13 point stencils) and

Dispersion Relation Preserving schemes from 7 to 13 point stencils are available. The code

is written in Fortran 90 and uses hybrid MPI/OpenMP parallelization. The code is applied

to the simulation of a Mach 0.9 jet flow. Four-stage third-order Runge-Kutta time stepping

and the 13 point DRP spatial discretization scheme of Bogey and Bailly are used. The high

resolution numerics used allows for the use of relatively sparse grids. Three levels of grid

resolution are examined, 3.5, 6.5 and 9.2 million points. Mean flow, first-order turbulent

statistics and turbulent spectra are reported. Good agreement with experimental data for

mean flow and first-order turbulent statistics is shown.

Nomenclature

an finite difference stencil coefficient

bn filter stencil coefficient

b shear layer thickness

c phase speed

c

∗ numerical phase speed

et total energy

f arbitrary function

i, j grid indices

k turbulent kinetic energy

p pressure

qi heat flux vector

r radial distance from centerline

t time

ui velocity tensor

u, v, w axial, radial and azimuthal velocity components

u

0

, v

0

, w

0 axial, radial and azimuthal turbulence intensities

xi coordinate tensor

x, y, z cartesian coordinates

Dj jet diameter

D spatial discretization operator

D(κ∆x) filter damping function

E(κ) energy spectral density

M Mach number, number of Runge-Kutta stages

N stencil half width

Re jet Reynolds number

A High-Resolution Capability for Large-Eddy 

Simulation of Jet Flows 

James R. DeBonis 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 

NASA/TM—2011-216833 1
T temperature

Uj ideal jet velocity

∆t time step

∆x grid spacing

α angle of attack

αm, βm Runge Kutta scheme coefficients

κ wave number

μ viscosity

ρ density

σ filter strength coefficient

τij stress tensor

subscripts

jet jet exit condition

0 stagnation/total condition

∞ freestream condition

superscripts

spatial filtered

e favre filtered

b filtered

sgs sub-grid scale

I. Introduction

The prediction of the flowfield of turbulent jets using computational fluid dynamics (CFD) is an important

research area. Turbulent jet flows are a critical component in propulsion, aeroacoustics, combustion and

low-observable (reducing the infrared signature of the exhaust) technologies. A better understanding of the

flow physics and improved predictive tools could lead to quieter, more efficient and less vulnerable aircraft.

The standard practice in the aerospace community for predicting these flows is using Reynolds Averaged

Navier-Stokes (RANS) based CFD codes. RANS methods are based on the time-averaged Navier-Stokes

equations and rely on a model to represent the effects of the turbulence. Since jet flows are dominated by

turbulent mixing, the turbulence model in the RANS solver is the critical component in the simulation.1

Despite decades of development, to date RANS has failed to produce reliable predictions of the jet flowfield.2

Standard models such as Menter’s SST model3 and k −

4

fail to accurately predict the potential core length

(the inviscid portion of the jet not effected by the shear layer) and shear layer growth rate. Turbulence

models specially formulated for round jets produce better results, but these models can not be generalized

for other types of jets.5In addition, because RANS produces a steady/time-averaged solution, no detailed

information on the turbulent structures is generated.

Large-eddy simulation (LES) is a CFD technique that offers the promise of improved predictions and

increased information for turbulent flows. LES is an unsteady solution technique that directly computes

the large-scale turbulent structures in the flow. Because the large structures carry the vast majority of the

turbulent kinetic energy, directly computing them should improve predictions. The effect of structures that

are too small be be resolved by the computation is modeled or neglected. While LES has the potential

for improved accuracy, up to this point RANS methods perform as well or better at a much lower cost for

compressible jets.6 A recent survey comparing several different LES methodologies for the prediction of jet

flows was done by Bodony and Lele.7In general the methods do a poor job of predicting the end of the

potential core, centerline velocity profiles and centerline turbulence intensities.

The key to an accurate LES is the proper resolution of the large-scale turbulent structures which dominate

the flow. The resolution of these structures is determined by a combination of the numerical scheme and

computational grid used. Improved resolution can be achieved by adding additional grid points or increasing

the resolution of the scheme. To quantify the ability of numerical schemes to resolve and convect wavelike

structures many researchers have used a Fourier analysis based technique.8 This technique has been a

valuable tool in the development and characterization of numerical schemes. Recent work has attempted to

extend this work to jet LES.9

NASA/TM—2011-216833 2
High resolution schemes are defined here as having the ability to resolve structures containing large wave

numbers with fewer grid points than the typically used second- and third-order schemes. High resolution

schemes may, but do not necessarily, have a high order of accuracy. Commonly used high resolution schemes

are compact schemes10 and Dispersion Relation Preserving (DRP) schemes.11, 12 These schemes are more

efficient for turbulent simulations; they can achieve a given level resolution with less computer resources than

low-order schemes. A drawback of high resolution schemes is the increased complexity required in the flow

solver. Large difference stencil sizes make it difficult to implement the numerical scheme near computational

boundaries and more information must be passed between grid interfaces. This makes using high resolution

schemes difficult for realistic problems with complex shapes.

This paper describes a new LES code and its application to jet flows. The code incorporates high

resolution numerical schemes for efficient simulations of turbulent flows. Details on the numerical schemes,

boundary conditions and parallel processing options are provided. Several features that facilitate the solution

of flows around complex geometries, generalized curvilinear coordinates, internal grid boundaries, hole cutting

and grid blocking, are included and their implementation with respect to the spatial differencing schemes

is discussed. The code is applied to a Mach 0.9 turbulent compressible jet, a standard test case, and its

solutions are compared against experimental data. Mean velocities and turbulence intensities are compared

to the experiment with good results.

II. Code Description

The code used in this study, WRLES (Wave Resolving Large-Eddy Simulation), is a special purpose largeeddy simulation code that uses high-resolution temporal and spatial discretization schemes to accurately

simulate the convection of turbulent structures. The code solves the compressible Favre-filtered NavierStokes equations. The code is written entirely in Fortran 90 and utilizes both Message Passing Interface

(MPI) libraries13 and OpenMP compiler directives.14 It was developed at the NASA Glenn Research Center

for the study of turbulent jets, but can be applied to other flows.

A. Governing Equations

The basis of large-eddy simulation is the separation of the large and small scale turbulent fluctuations. The

large-scale turbulence is computed directly using the Navier-Stokes equations. The small scale turbulence is

modeled. Since the large-scales carry the vast majority of the turbulent kinetic energy the technique offers

promise for accurate turbulent simulations. Since the small scale turbulence serves to dissipate the large

scales and is isotropic, a simple model can be constructed to impose the effect of the small scales on the rest

of the flow.

A filtering process is applied to the continuity, momentum, and energy equations. The resulting equations are comprised of resolved and unresolved terms. The resolved terms in the filtered equations directly

correspond, in form, to the unfiltered equations. The resulting LES expressions for conservation of mass,

momentum, and energy are:

∂ρ

∂t +

∂

∂xi

(ρu ̃i) = 0 (1)

∂

∂t(ρu ̃i) + ∂∂xj

(ρu ̃iu ̃j + p) = ∂

∂xj



τ ij + τ

sgs

ij 

(2)

∂

∂t (ρe ̃t) + ∂∂xi

(ρu ̃ie ̃t +  ̃uip) = ∂

∂xi

u ̃j τ ij +  ̃uj τ

sgs

ij 

− (qi + q

sgs

i

)



(3)

The overbar,  ̄( ) represents a spatially filtered quantity and the tilde,  ̃( ) represents a Favre-averaged

quantity;  ̃f = (ρf)/ρ ̄

The unclosed terms from the LES momentum and energy equations are the sub-grid scale stress, τ

sgs

ij , and

the sub-grid scale heat flux, q

sgs. In a traditional LES approach a sub-grid model would be used to compute

τ

sgs

ij and q

sgs. The Smagorinsky sub-grid model is included in the code.15 However, many practitioners

forego sub-grid modeling, relying on the dissipation implicit in the numerical scheme to dissipate the energy

at the small-scales. This approach is sometimes called Implicit LES (ILES).

NASA/TM—2011-216833 3
B. Numerical Method

The Favre filtered Navier-Stokes equations are discretized and solved using the explicit numerical methods

described below.

1. Time Discretization

The code uses a family of explicit Runge-Kutta time stepping schemes written in a general M-stage 2-N

storage formulation.16

dfm = αmdfm−1 + ∆tD (fm−1) (4a)

fm = fm−1 + βmdfm (4b)

for m = 1 . . . M, and where the initial value at the start of the iteration is f0 = f

n and the final value is

fM = f

n+1. The coefficient α1 is typically set to zero for the algorithm to be self-starting.

The operator D is the spatial finite difference operator. One-stage is required for each order of accuracy

desired. Additional stages can be used to increase accuracy or alternatively, to reduce dispersion error. The

temporal discretization can be varied by simply changing the number of stages and the coefficients, αm and

βm. Table 1 lists the schemes that are currently implemented in the code.

Scheme Number of Stages Order of Accuracy

Euler forward 1 1

Second order 2 2

Williamson16 3 3

Gottlieb & Chu, TVD17 3 3

Carpenter & Kennedy18 4 3

Carpenter & Kennedy18 5 4

Stanescu & Habashi19 5 4

Stanescu & Habashi19 6 4

Table 1. Low-dispersion Runge-Kutta schemes implemented in the code

2. Spatial Discretization

Central differencing is used for the spatial discretization because of its non-dissipative properties. This helps

ensure the accurate convection of turbulent structures. The central difference stencil can be written for an

arbitrary stencil size

∂f

∂x









i

=

X

N

n=1

1

∆x

[an(fi+n − fi−n)] (5)

where the half width of the stencil, the number of points on either side of the central point, is N. The total

number of points in the stencil is 2N + 1. The spatial discretization can be varied by simply changing the

width of the stencil and the coefficients, an. Standard stencils from 2nd- to 12th-order are included. In

addition Dispersion Relation Preserving (DRP) stencils from Tam,11 7-point, and Bogey & Bailly,12 7-, 9-

and 13-point, are included in the code.

Standard central difference stencils are designed to minimize the truncation error for a given stencil

size. In DRP schemes, the coefficients in the stencil are chosen to minimize the dispersion error of the

scheme, rather than the truncation error. DRP schemes are designed to accurately capture and convect

large structures, which are well resolved by the grid, and are generally a better choice for LES.

3. Filtering

The lack of dissipation in central differencing makes the schemes unstable. To ensure a stable solution

without adversely affecting the resolving properties of the scheme, solution filtering is used. This is a

NASA/TM—2011-216833 4
low-pass filter that leaves the low-wavenumber well-resolved structures untouched and removes the highwavenumber unresolved structures that can cause instability. The filter must be properly matched to the

differencing scheme, so that the filter removes only those waves that are not properly computed.

ˆfi = fi + σ

X

N

n=−N

bnfi+n (6)

For the standard central difference schemes the filters of Kennedy and Carpenter20 are implemented in the

code. Bogey and Bailly developed filters to match their DRP stencils and are included.

4. Properties of the Numerical Scheme

Fourier analysis has been used to evaluate the ability of numerical schemes to resolve wave motion.8 The

technique analyzes the numerical scheme on the one-dimensional convection equation and can provide the

scheme’s behavior in terms of dissipative and dispersive errors. The effect of the temporal discretization is

not considered here.

Fourier analysis shows that central difference schemes have no inherent dissipation, and thus are ideally

suited for LES calculations. They do however produce dispersive errors. Figure 1(a) shows the phase error,

written as the ratio of numerical phase speed to actual phase speed, c

∗/c, versus wavenumber per grid

spacing, (κ∆x), for both standard schemes and the DRP schemes. If one selects a maximum acceptable

phase error for a simulation then the maximum resolvable wave number per grid spacing, (κ∆x)max, is

determined for each stencil. The maximum acceptable error used was 5 · 10−4.

The damping function of the filter can also be expressed in terms of wave number per grid spacing. The

damping functions, D(κ∆x), are shown in figure 1(b) for the standard filters of Kennedy and Carpenter20

and the DRP filters derived by Bogey and Bailly to match their difference stencils. The damping is zero

until it reaches a “cutoff” wave number, where it rapidly increases, effectively removing all the structures of

higher waves numbers. Since filters typically also utilize a coefficient, σ, to regulate the amount of damping,

the total amount of dissipation is σD(κ∆x). The cutoff wave number, (κ∆x)cut is the wave number where

the total amount of dissipation exceeds 5 · 10−4, the same error level used for the differencing schemes. A

typical damping coefficient of 0.2 is used here.

!"x

c

*

/c

0 0.5 1 1.5 2 2.5 3 0

0.2

0.4

0.6

0.8

1

2nd order

4th order

6th order

8th order

10th order

12th order

Tam & Shen

B & B, 9 pt.

B & B, 11 pt.

B & B, 13 pt.

(a) Error in phase speed for the finite difference stencils

!"x

D(!"x)

0

0.2

0.4

0.6

0.8

1 2nd order

4th order

6th order

8th order

10th order

12th order

B&B 9pt.

B&B 11pt.

B&B 13pt.

#/2 #

(b) Damping functions of the explict filters

Figure 1. Properties of the numerical scheme

The maximum resolvable wave number for the 8th, 10th, and 12th order standard schemes and the 9-,

11- and 13-point DRP schemes is plotted against the stencil size as solid lines in figure 2. The cutoff wave

number of the corresponding filter is plotted as a dashed line. For a stable solution, the cutoff wave number

of the filter should be less than the maximum resolvable wave number of the numerical scheme. This ensures

NASA/TM—2011-216833 5
that all of the structures within the domain are properly resolved. For the standard schemes, the filter’s

cutoff is too high, leaving improperly resolved structures in the flowfield. For these schemes, choosing the

next lower order filter effectively shifts the filter curve to the right and should insure a stable solution. Figure

2 clearly illustrates superior resolution offered by the DRP schemes.

Stencil Size

!"

x

3 5 7 9 11 13 0

0.4

0.8

1.2

1.6

2

Std. Diff.

Std. Filter

Shifted Filter

DRP Diff.

DRP Filter

Figure 2. Cutoff wave numbers for schemes and filters

C. Geometry Handling

Many high resolution/high-order of accuracy CFD solvers are limited to canonical problems with very simple

geometries. The WRLES code contains several features that allow the solution of flows about moderately

complex geometries.

1. Generalized Curvilinear Coordinates

The governing equations are cast in generalized curvilinear coordinates allowing the computational grids to

stretch and deform, fitting the geometric surface. The code has options to solve these equations in either the

strong conservation form or the chain rule form.21 The strong conservation form is used in most codes.22

However, Hixon et al.23 recently showed that in practice, the chain rule form of the equations is more

accurate than the strong conservation form when the metric terms are computed numerically. There are

two practical considerations when choosing the form of the equations to solve; 1) the strong form is more

computationally efficient, but 2) the chain rule form is more stable when using the centerline differencing

condition for o-grid singularities. Grid metrics for the coordinate transformation are computed using the

same differencing scheme that is selected for the flow solver.

2. Internal Boundaries and Holes

To allow greater flexibility in generating grids, boundary conditions can be specified on any portion of

any grid surface in the domain. In addition, holes within the computational domain can be created. The

combination of these features allows for complex geometric features to be embedded within a grid block.

3. Near Boundary Stencils

As you approach the boundaries, the full width of the difference stencil can not be maintained. To preserve the

non-dissipative and non-dispersive characteristics of the scheme, the stencil width is reduced in the direction

of the boundary and central differencing stencils are maintained. This approach reduces the resolution of the

scheme near the boundaries and care must be taken to not adversely affect the solution. Wherever possible

the boundaries should be placed as far from the region of interest as possible. For the jet problem, only the

NASA/TM—2011-216833 6
very start of the mixing layer will be affected by the smaller stencils. Figure 3 illustrates the reduction of

the stencil half widths, N, near an internal boundary representing and embedded object. The boundary is

represented by the red line and the hole points are shown by the x’s.

6 

6 

6 

6 

6 

6

X  X  X  X X 

X  X  X  X X 

6 

6 

6 

6 

6 

6

6 

6 

6 

6 

6 

6

6 

6 

6 

6 

6 

6

6 

6 

6 

6 

6 

6

6 

6 

6 

6 

6 

6

X  X  X  X X 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

1  2  3  4 5  6 

1  2  3  4 5  6 

1  2  3  4 5  6 

1  2  3  4 5  6 

i

j

(a) i-direction

3 

2 

1 

4 

5 

6

X  X  X  X X 

X  X  X  X X 

3 

2 

1 

4 

5 

6

3 

2 

1 

4 

5 

6

3 

2 

1 

4 

5 

6

3 

2 

1 

4 

5 

6

3 

2 

1 

4 

5 

6

X  X  X  X X 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

6  6  6  6 6  6 

i

j

(b) j-direction

Figure 3. Stencil half widths, N, near an embedded object

4. Grid Blocking

The code is also capable of computing grids made up of multiple grid blocks. In order to maintain consistent

numerical resolution across the block boundaries, the blocks must be point to point matched and share

several common planes. The number of common planes required is 2N, one less than the stencil width. This

restriction ensures that every grid point in the domain is computed using a full central difference stencil.

Block interfaces for stencil widths of 5 and 13 are currently included in the code. A schematic showing the

overlap requirement and the resulting stencil half widths at each point for a 13-point stencil is given in figure

4.

1  2  3 4  5  6  6  6 6  6  6 6 

6  6  6  6  6  6  5  4  3  2  1 

0 

6  0 

6 

6 

Figure 4. Stencil half widths, N, at overlapping block boundary

D. Boundary Conditions

Boundary conditions are available to simulate a wide variety of flows. Because the width of the central

difference stencil is reduced as the boundary is approached, reducing the resolution, most of the boundary

conditions do not attempt to maintain high resolution. This is true for wall, inflow and outflow boundaries.

High resolution is maintained for periodic, cylindrical centerline and block interface boundaries. The formulation of the boundary conditions is standard within the aerospace community. Table 2 lists the available

boundary conditions.

NASA/TM—2011-216833 7
Boundary Condition Condition(s) enforced

subsonic inflow p0, T0

supersonic inflow ρ, ρu, ρv, ρw, ρet

characteristic boundary M∞, α

subsonic outflow p∞

supersonic outflow none

inviscid (slip) wall ∂ui/∂n = 0

viscous (no-slip) wall ui = 0, ∂q/∂n = 0

singular axis

periodic

block interface

centerline differencing

exit zone

Table 2. Standard boundary conditions

An o-grid that forms a cylindrical coordinate system is the natural choice for a round jet problem.

However this grid type creates a singularity on the centerline causing two problems; 1) a boundary condition

must be applied in the center of the domain where the flow should be computed and 2) the small grid cells at

the singularity create numerical instabilities. The centerline differencing condition removes both problems

by removing the singularity, increasing the cell size, and creating a difference stencil across the centerline so

that the equations can be solved there. The formulation is taken from Hixon24 and the concept is illustrated

in figure 5. The grid is constructed with a cylindrical void on the centerline to avoid a collapsed surface.

The cell spacing is set such that the diameter of the void is equal to the adjacent radial grid spacing. Then

a stencil is constructed in such a manner that the points on the opposite side of the singularity are used.

j,k

j+1,k 

j+2,k 

j+3,k 

j-1 

j-1, k+kmax/2 

j, k+kmax/2 

Figure 5. Sample differencing stencil across the centerline

The exit zone condition applies additional damping near an outflow boundary to remove outgoing waves

that would otherwise reflect and contaminate the solution. The damping is introduced through grid stretching

and additional low-order filtering whose strength is gradually increased as the boundary is approached.

NASA/TM—2011-216833 8
E. Parallel Processing

WRLES is written using hybrid parallel constructs to take advantage of multi-processor and multi-core Linux

cluster systems. The hybrid parallel implementation uses both Message Passing Interface (MPI) libraries and

OpenMP compiler directives. The grid can be divided into multiple blocks and computation of each block

can be distributed to the nodes of the Linux cluster using the MPI libraries. OpenMP compiler directives

parallelize the computation at the loop level across the multiple processors and multiple cores within each

computer node. The requirement of grid overlap adds additional grid points to the computation when a

domain is subdivided for parallelization. This limits the parallel scalability for large numbers of processors.

III. Simulation of a Mach 0.9 Jet

The WRLES code was used to obtain flowfields for a Mach 0.9 turbulent cold jet, a standard test case

for jet LES. The configuration used for this study was the 2 inch diameter Acoustic Reference Nozzle (ARN)

tested at the NASA Glenn Research Center by Bridges and Wernet.25 The experimental data was obtained

using Particle Image Velocimetry (PIV). The data was obtained on a streamwise plane through the jet

centerline downstream of the nozzle exit from 1 ≤ x/Dj ≥ 25 and −1.4 ≤ r/Dj ≥ 1.4. Three components of

velocity and turbulence intensities are used for comparison to the calculations. The jet Reynolds number,

Rej = ρjUjDj/μj , from the experiment was approximately 2 million. In order to reduce the range of

turbulent scales in the simulation, the Reynolds number of the simulation was reduced to 50,000.

A. Computational Approach

1. Computational Grid

A blocked structured grid was created using the Gridgen software package.26 The computational grid (figure

6) models the internal and external nozzle geometry from plenum to exit. The majority of grid points are

placed in the area of interest, downstream of the nozzle exit in the plume region. The grid extends 40 jet

diameters, Dj , downstream of the nozzle exit and 30 diameters radially from the centerline. The nozzle lip

thickness was 0.015Dj (0.03 inches) and 30 grid points were used to model it. For the viscous walls, the grid

spacing at the wall was set to 3 · 10−4Dj (6 · 10−4inches).

x/Dj

y/Dj

0 5 10 15 20 25

-4

-2

0

2

4

Figure 6. Computational grid, xy-plane (every other grid line removed)

The grid is constructed in 3 grid blocks; interior nozzle flowpath, exterior nozzle flowpath and nozzle

plume. To ensure that each grid point is computed using the maximum stencil width possible the nozzle

plume grid block is overlapped into the interior and exterior nozzle blocks. The grid points that fall in

between the two upstream blocks are designated as hole points and not computed. Figure 7 shows a detail of

this region. The overlapped region, hole points (in purple) and coincident points can be seen in the exploded

view.

NASA/TM—2011-216833 9
x/Dj

y/Dj

-0.5 0 0.5 1

0.5

1

1.5

(a) 3-block grid around nozzle lip

x/Dj

y/Dj

-0.5 0 0.5 1

0.5

1

1.5

(b) exploded view, showing overlap and i-blanked region

Figure 7. Modeling of nozzle lip

Three different grid densities were examined (table 3). The initial coarse grid was constructed based on

best practices for a Reynolds Averaged Navier-Stokes (RANS) simulation. For this grid, the smallest grid

spacings were in the radial direction, due to the clustering near solid walls, and the largest were in the axial

direction. The ability to resolve a three-dimensional turbulent eddy is limited by the largest grid spacing in

each of the three directions. Therefore the first refined grid maintained the radial grid spacing and added

points in the axial and azimuthal directions. The second refinement added points only in the axial direction.

Grid Interior of Nozzle Exterior of Nozzle Nozzle Plume Total

coarse 45 x 55 x 102 45 x 65 x 102 196 x 148 x 102 3,509,616

medium 45 x 55 x 132 45 x 65 x 132 294 x 148 x 132 6,456,384

fine 45 x 55 x 132 45 x 65 x 132 435 x 148 x 132 9,210,960

Table 3. Dimensions of the computational grid

2. Boundary Conditions

Total pressure and total temperature, consistent with a Mach 0.9 cold jet, were specified at the inflow of the

nozzle plenum. Freestream conditions were specified using a farfield characteristic boundary condition. The

experiment had no freestream velocity, but a very low Mach number, 0.05, was specified at the computation’s

farfield boundary to maintain stability. Static pressure was specified at the outflow and an exit zone was

initiated at x/Dj > 30; the grid spacing is gradually increased and a sixth-order filter is applied to remove

spurious reflections. The solution was computed at the jet centerline using the singularity differencing

boundary condition.

3. Solution Parameters

Temporal discretization was done using the four-stage third-order low-dispersion Runge-Kutta scheme developed by Carpenter and Kennedy.18 Spatial discretization was done using the 13-point DRP scheme

developed by Bogey and Bailly.12 The filter corresponding to this scheme was used for damping with a

coefficient of 0.5.

NASA/TM—2011-216833 10
For these calculations the implicit large-eddy simulation (ILES) approach was used. In this approach,

no model is used for the sub-grid scale turbulent stresses. The solution relies on the dissipation of the filter

to represent the sub-grid scales.

The solution does not utilize any artificial means to generate or initialize turbulence in the flowfield.

The grid for the nozzle boundary layers is not resolved enough to sustain turbulent structures, resulting in

laminar flow at the nozzle exit. At the nozzle exit, vortex shedding from the nozzle lip introduces large-scale

unsteadiness into the jet mixing layer, which rapidly transitions into turbulence.

The initial 3 block grid was subdivided into 16 blocks for parallel processing. The code was run on a

linux cluster using one master node and 16 dual-processor dual-core worker nodes. The block interface data

was shared between nodes via MPI libraries over an Infiniband network. The computational work on each

node was parallelized across the four processor cores using OpenMP.

4. Time Advancement

The solution was run by specifying a CFL number. The time step at each grid point is computed using

that CFL, and the minimum value for the entire domain is used to advance the solution everywhere. The

maximum stable CFL for this calculation was 0.7, which corresponds to a physical time step of approximately

2.45 · 10−8seconds. The solution was saved at a Nyquist rate of 2.5 · 10−5seconds, yielding a maximum

frequency of 20kHz. A total of 2,048 solutions, representing 0.0512 seconds of physical time (7.7 flow through

times) were saved for post processing.

B. Results

Instantaneous Mach number contours for the medium grid are shown in figure 8. Upstream of the nozzle

exit the flow is steady. Downstream, the flow is dominated by turbulent mixing. Immediately downstream

of the nozzle lip relatively large structures are seen; the structures are on the order of the the nozzle lip

thickness. In general the size of the structures increases with distance from the nozzle exit and grow to be

on the order of the nozzle exit diameter. The thin void along the centerline of the domain is an artifact of

the centerline boundary condition.

Figure 8. Instantaneous Mach contours, M

1. Prediction of the Mean Flowfield

The prediction of the mean velocity along the jet centerline for all three grids is compared to experiment in

figure 9. Overall, agreement with experiment is very good for all three grid levels. The data shows that the

predicted end of the potential core (the region of jet flow unaffected by viscous effects) shifts downstream

with the initial grid refinement, more closely matching the data. Further refinement shifts the end of the

potential core slightly back downstream. Downstream of the potential core, the decay rate of the centerline

velocity increases with grid refinement. This is most likely due to the additional turbulent structures that

are resolved by the finer grids.

NASA/TM—2011-216833 11
x/Dj

u/Uj

0 5 10 15 20 0.2

0.4

0.6

0.8

1

1.2

Coarse Grid

Medium Grid

Fine Grid

Experiment

Figure 9. Centerline velocity profile

Radial profiles of mean velocity at four axial locations, x/Dj = 1, 4, 7 and 10, are shown in figure 10.

Again, agreement with experiment is very good. The profiles indicate that the computations predict a thicker

initial shear layer thickness than the experiment. Grid refinement appears to reduce the thickness, moving

the predictions closer to the data.

u/Uj

r/Dj

0 1.2 2.4 3.6 4.8 0

0.2

0.4

0.6

0.8

1

1.2

1.4

Coarse Grid Medium Grid Fine Grid Experiment

x/Dj = 1 x/Dj = 4 x/Dj = 7 x/Dj = 10

Figure 10. Radial profiles of axial velocity

To help quantify jet mixing a shear layer thickness parameter was defined. The inner edge of the shear

layer was defined as the point where the velocity was 90 percent of the jet velocity and the outer edge was

defined as the point where the velocity was 10 percent of the jet velocity.

b = r







u=0.1Uj

− r







u=0.9Uj

(7)

The thickness, b, based on the edge points was computed from the jet exit to the end of the potential core

and is shown in figure 11. The experiment indicates a near linear growth rate from the jet exit. All three

computations show that the shear layer grows rapidly for the first diameter and then relaxes to a growth

NASA/TM—2011-216833 12
rate similar to the experiment. The initial thickness does not improve with grid resolution. This could

be caused by; 1) the the lack of turbulence entering the shear layer from the nozzle boundary layers, 2)

the fact that even the finest grid can not resolve the small scale turbulence here, or likely a combination

of the two. Downstream the growth rates improve with grid resolution and the fine grid solution growth

rate approaches the experimental data. For the coarse grid, the error in the initial shear layer thickness

and the slow downstream growth rate compensate for each other. This is the reason the mean centerline

velocity profile appears correct for this case. For the fine grid the initial shear layer is still too thick, but the

downstream growth rate matches the experimental growth rate closely. The poor prediction of the initial

thickness is the reason the fine grid centerline velocity decays too early.

x/Dj

b/Dj

0 2 4 6 8 0

0.2

0.4

0.6

0.8

1

1.2

Coarse Grid

Medium Grid

Fine Grid

Experiment

Figure 11. Shear layer thickness

2. Prediction of Turbulent Statistics

Contours of the three components of turbulence intensity are shown in figures 12 - 14. The three levels of

grid resolution are compared to experimental data. For all components of intensity and all levels of grid

resolution the computations show a high level of turbulence just downstream of the nozzle lip and persisting

for 2 nozzle diameters. The experimentalists were unable to take accurate measurements in the first jet

diameter downstream of the nozzle exit. Nevertheless based on the existing data one can infer that the high

predicted levels of turbulence near the nozzle lip are nonphysical. Because the nozzle boundary layers are

not properly resolved, large laminar-like vortices are shed. The situation is exacerbated by the lack of grid

resolution in this region necessary to resolve small turbulent structures. These large vortices serve to disturb

the flow and trigger the flow to transition to a realistic turbulent state.

The contours show that the predictions of turbulence intensity are greatly improved with increasing grid

resolution. For all three components of turbulence intensity the fine grid predictions beyond x/Dj = 5 are

very good.

The discrepancy near the nozzle lip can also be seen on line plots of turbulence intensity along the jet lip

line, r/Dj = 0.5 (figure 15). A sharp peak of about twice the expected intensity in the downstream region is

present near x/Dj = 0.5. For the coarse grid, levels of axial intensity are overpredicted and the radial and

azimuthal values are underpredicted. Grid refinement reduces axial levels and increases radial and azimuthal

levels improving agreement with experiment for all components.

Axial, u

0/Uj , and radial, v0/Uj , turbulence intensities on the jet centerline are shown in figure 16. On the

baseline grid, the axial intensity levels are over predicted and the radial intensity levels are under predicted.

Grid refinement improves both predictions, decreasing the axial intensity and increasing the radial intensity.

Even though the flow at the centerline upstream of the end of the potential core is inviscid, the effect of

the large structures at the nozzle lip are felt on the centerline as evidenced by a peak near x/Dj = 0. The

disturbances from the jet lip are propagated to the centerline through pressure fluctuations.

NASA/TM—2011-216833 13
(a) coarse grid

(b) medium grid

(c) fine grid

(d) experiment

Figure 12. Contours of axial turbulence intensity, u

0/Uj

NASA/TM—2011-216833 14
(a) coarse grid

(b) medium grid

(c) fine grid

(d) experiment

Figure 13. Contours of radial turbulence intensity, v

0/Uj

NASA/TM—2011-216833 15
(a) coarse grid

(b) medium grid

(c) fine grid

(d) experiment

Figure 14. Contours of azimuthal turbulence intensity, w0/Uj

NASA/TM—2011-216833 16
x/Dj

u’/Uj

0 5 10 15 20 0

0.05

0.1

0.15

0.2

0.25

0.3

Coarse Grid

Medium Grid

Fine Grid

Experiment

(a) axial intensity

x/Dj

v’/Uj

0 5 10 15 20 0

0.05

0.1

0.15

0.2

0.25

0.3

Coarse Grid

Medium Grid

Fine Grid

Experiment

(b) radial intensity

x/Dj

w’/Uj

0 5 10 15 20 0

0.05

0.1

0.15

0.2

0.25

0.3

Coarse Grid

Medium Grid

Fine Grid

Experiment

(c) azimuthal intensity

Figure 15. Lip-line turbulence intensities

NASA/TM—2011-216833 17
x/Dj

u’/Uj

0 5 10 15 20 0

0.05

0.1

0.15

0.2 Coarse Grid

Medium Grid

Fine Grid

Experiment

(a) axial intensity

x/Dj

v’/Uj

0 5 10 15 20 0

0.05

0.1

0.15

0.2 Coarse Grid

Medium Grid

Fine Grid

Experiment

(b) radial intensity

Figure 16. Centerline turbulence intensities

Radial profiles of turbulence intensity at x/Dj = 1, 4, 7 and 10 are shown in figure 17. The trends seen

in figures 12 - 16 are further quantified here: 1) the LES overpredicts the turbulence levels near the nozzle

lip, 2) further downstream the turbulence levels are slightly underpredicted and 3) the predictions improve

universally with grid refinement.

The turbulence intensity data clearly show that grid refinement greatly improves the LES predictions.

This fact was not evident in the mean velocity data, and the opposite conclusion could have been drawn.

The effect of two competing sources of error, the poor resolution of the initial shear layer and the predicted

growth rate downstream compensated for each other in the mean flow data.

It is also important to note that the prediction of all components of turbulence intensity were improved by

refining the axial and azimuthal grid spacing for the first refinement and only the axial spacing for the second

refinement. In other words, improving the grid resolution in the coarsest direction has a significant effect.

This illustrates the importance of resolving the turbulent structures in all three computational directions

and suggests that isotropic grid cells are ideal for LES. This is a significant shift from the practices of RANS

analyses where grid refinement is done to resolve mean flow gradients.

3. Prediction of Turbulent Spectra

Turbulent spectra were computed at three locations in the jet mixing layer. The points were located on the jet

lip line, r/Dj = 0.5 at three axial locations, x/Dj = 4, 7, and 10. At each of these points, the three velocity

components for each saved solution were extracted, yielding a time history of the velocity. The data was

converted to instantaneous turbulent kinetic energy and was processed using a Fast Fourier Transform. This

yielded a turbulent spectrum in the frequency domain. The flow was assumed to be statistically stationary

and Taylor’s hypothesis27 was used to convert the frequency domain to the spatial domain.

Plots of the turbulent spectra for both grids are shown in figure 18. The majority of the turbulent energy

is contained in the small wave number/large-scale structures. A dashed line with a −5/3 slope indicating

the inertial sub-range is fitted to each spectra. The inertial sub-range is not well resolved (some may argue

whether or not it is present in the simulation at all). Beyond this point the energy decays very rapidly, most

likely the effect of the filter. The spectra show that the energy in the larger scales increases with downstream

location. This is consistent with the fact that the turbulent structures grow in size as they move downstream.

Grid refinement increases the energy in the small scales and shifts the start of the inertial sub-range, the

−5/3 slope, to higher wave numbers. The apparent extent of the inertial sub-range increases slightly with

grid refinement. It is interesting to note that the turbulence field is well predicted on the fine grid despite a

lack of resolution of the turbulent inertial sub-range.

NASA/TM—2011-216833 18
u’/Uj

r/Dj

0 0.2 0.4 0.6 0.8 0

0.5

1

1.5

x/Dj = 1 x/Dj = 4 x/Dj = 7 x/Dj = 10

Coarse Grid Experiment Medium Grid Fine Grid

(a) axial intensity

v’/Uj

r/Dj

0 0.2 0.4 0.6 0.8 0

0.5

1

1.5

x/Dj = 1 x/Dj = 4 x/Dj = 7 x/Dj = 10

Coarse Grid Experiment Medium Grid Fine Grid

(b) radial intensity

w’/Uj

r/Dj

0 0.2 0.4 0.6 0.8 0

0.2

0.4

0.6

0.8

1

1.2

1.4 x/Dj = 1 x/Dj = 4 x/Dj = 7 x/Dj = 10

Coarse Grid Experiment Medium Grid Fine Grid

(c) azimuthal intensity

Figure 17. Radial profiles of turbulence intensity

NASA/TM—2011-216833 19
!Dj

E(k)/(D

Uj

2j )

10-2 10-1 100 101 10-3

10-2

10-1

100

Coarse Grid

Medium Grid

Fine Grid

-5/3 Slope

(a) x/Dj = 4

!Dj

E(k)/(D

Uj

2j )

10-2 10-1 100 101 10-3

10-2

10-1

100

Coarse Grid

Medium Grid

Fine Grid

-5/3 Slope

(b) x/Dj = 7

!Dj

E(k)/(D

Uj

2j )

10-2 10-1 100 101 10-3

10-2

10-1

100

Coarse Grid

Medium Grid

Fine Grid

-5/3 Slope

(c) x/Dj = 10

Figure 18. Turbulent kinetic energy spectra

NASA/TM—2011-216833 20
IV. Summary and Conclusions

A code for large-eddy simulation of turbulent jet flows was developed. The code employs low-dispersion

Runge-Kutta time discretization and central differencing spatial differencing. It is coded in a general manner

such that any number of temporal and spatial schemes that can be written in the given form can be used.

Fourier analysis of the spatial discretization scheme was used to characterize the ability of the scheme to

resolve waves on a computational grid and it was shown that the Dispersion Relation Preserving schemes

are superior to standard central differencing schemes for the resolution of turbulent structures. Details on

the implementation of the numerical scheme for generalized curvilinear coordinates, blocked structured grids

and internal boundaries are provided.

The code was applied to the flowfield of a Mach 0.9 cold jet. A four-stage third-order low-dispersion

Runge-Kutta scheme was used for the temporal discretization and the 13-point Dispersion Relation Preserving scheme of Bogey and Bailley12 was used for the spatial discretization. A filter designed to match the

spatial discretization scheme was used to maintain stability and to represent the sub-grid scale dissipation.

Three different grids of 3.5, 6.5 and 9.2 million points were examined. The results of the computation were

compared to the experimental data of Bridges and Wernet.25 Excellent agreement was obtained for the mean

velocities at all grid resolutions. Good agreement was achieved for the turbulence intensities for the baseline

grid and the agreement improved significantly with grid refinement. The results represent an improvement

over the previously cited RANS2 and LES7results. The improvement with selected grid refinement suggests

that grid cells that tend toward isotropic are ideal for predicting turbulent quantities. Results indicated that

the prediction of the initial portion of the shear layer 0 ≤ x/Dj ≤ 2 is a limiting factor in the accuracy of

the simulation and further research is needed to improve predictions in this area.

Turbulent energy spectra were obtained at three points in the flowfield for both grids. The spectra

appear reasonable, but lack a well defined inertial sub-range. A −5/3 slope turbulent decay can be seen for a

small range of wave numbers. Despite a lack of resolution of the turbulent inertial sub-range, the turbulence

intensities were well represented on the fine grid.

References

1Barber, T. J., Chiappetta, L. M., DeBonis, J. R., Georgiadis, N. J., and Yoder, D. A., “An Assessment of Parameters

Influencing the Prediction of Shear Layer Mixing,” Journal of Propulsion and Power, Vol. 15, No. 1, 1999, pp. 45–53.

2Georgiadis, N. J. and DeBonis, J. R., “Navier-Stokes Analysis Methods for Turbulent Jet Flows with Application to

Aircraft Exhaust Nozzles,” Progress in the Aerospace Sciences, Vol. 42, 2006, pp. 377–418.

3Menter, F. R., “Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications,” AIAA Journal, Vol. 32,

No. 8, 1994, pp. 1598–1605.

4Jones, W. P. and Launder, B. E., “The Prediction of Laminarization with a Two-Equation Model of Turbulence,”

International Journal of Heat and Mass Transfer, Vol. 15, 1972, pp. 301–314.

5Georgiadis, N. J., Yoder, D. A., and Engblom, W. A., “Evaluation of Modified Two-Equation Turbulence Models for Jet

Flow Prediction,” AIAA Journal, Vol. 44, No. 12, 2006, pp. 3107–3114.

6DeBonis, J. R., “Progress Toward Large-Eddy Simulations for Prediction of Realistic Nozzle Systems,” Journal of Propulsion and Power, Vol. 23, No. 5, 2007, pp. 971–980.

7Bodony, D. J. and Lele, S. K., “Review of the Current Status of Jet Noise Predictions Using Large-Eddy Simulations,”

AIAA Paper 2006-0468, 2006.

8Vichnevetsky, R. and Bowles, J. B., Fourier Analysis of Numerical Approximations of Hyperbolic Equations, Society for

Applied and Industrial Mathematics, 1982.

9DeBonis, J. R., “An Examination of the Spatial Resolution Requirements for LES of a Compressible Jet,” Quality and

Reliability of Large-Eddy Simulation II , Springer, 2010, in process.

10Lele, S. K., “Compact Finite Difference Schemes with Spectral-Like Resolution,” Journal of Computational Physics,

Vol. 103, 1992, pp. 16–42.

11Tam, C. K. W. and Webb, J. C., “Dispersion Relation-Preserving Finite Difference Schemes for Computational Aeroacoustics,” Journal of Computational Physics, Vol. 107, 1993, pp. 262–281.

12Bogey, C. and Bailly, C., “A Family of Low Dispersive and Low Dissipative Explicit Schemes for Flow and Noise

Computations,” Journal of Computational Physics, Vol. 194, 2004, pp. 194–214.

13Forum, M. P. I., MPI: A Message-Passing Interface Standard, Version 2.2 , High Performance Computing Center

Stuttgart (HLRS), 2009.

14Chapman, B., Jost, G., and van der Pas, R., Using OpenMP: Portable Shared Memory Parallel Programming, The MIT

Press, 2007.

15Smagorinsky, J., “General Circulation Experiments with the Primitive Equations, Part I: The Basic Experiment,” Monthly

Weather Review, Vol. 91, 1963, pp. 99–152.

16Williamson, J. H., “Low-Storage Runge-Kutta Schemes,” Journal of Computational Physics, Vol. 35, 1980, pp. 48.

17Gottlieb, S. and Chu, C. W., “Total Variation Diminishing Runge-Kutta Schemes,” NASA CR 201591, 1996.

NASA/TM—2011-216833 21
18Carpenter, M. H. and Kennedy, C. A., “Fourth-Order 2N-Storage Runge-Kutta Schemes,” NASA TM 109112, 1994.

19Stanescu, D. and Habashi, W. G., “2N-Storage Low Dissipation and Dispersion Runge-Kutta Schemes for Computational

Acoustics,” Journal of Computational Physics, Vol. 143, No. 12, 1998, pp. 674–681.

20Kennedy, C. A. and Carpenter, M. H., “Comparison of Several Numerical Methods for Simulation of Compressible Shear

Layers,” NASA TP 3484, 1997.

21Tannehill, J. C., Anderson, D. A., and Pletcher, R. H., Computational Fluid Mechanics and Heat Transfer, Taylor &

Francis, 2nd ed., 1997, pp. 338–341.

22Vinokur, M., “Conservation Equations of Gas-Dynamics in Curvilinear Coordinate Systems,” Journal of Computational

Physics, Vol. 14, 1974, pp. 105–125.

23Hixon, R., Shih, S. H., Dong, T., and Mankbadi, R. R., “Evaluation of Generalized Curvilinear Coordinate Transformations Applied to High-Accuracy Finite-Difference Schemes,” AIAA Paper 98-0370, 1998.

24Hixon, R., Shih, S.-H., and Mankbadi, R. R., “Numerical treatment of cylindrical coordinate centerline singularities,”

AIAA Paper 1999-2391, 1999.

25Bridges, J. E. and Wernet, M. P., “Measurements of the Aeroacoustic Sound Source in Hot Jets,” AIAA Paper 2003-3130,

2003.

26Pointwise, Inc., http://www.pointwise.com/gridgen/.

27Pope, S. B., Turbulent Flows, Cambridge University Press, 2000, pp. 223–224.

NASA/TM—2011-216833 22
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

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-03-2011 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

A High-Resolution Capability for Large-Eddy Simulation of Jet Flows 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

DeBonis, James, R. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 984754.02.07.03.17.04 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17463 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2011-216833 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 34, 07, and 64 

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

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A large-eddy simulation (LES) code that utilizes high-resolution numerical schemes is described and applied to a compressible jet flow. The 

code is written in a general manner such that the accuracy/resolution of the simulation can be selected by the user. Time discretization is 

performed using a family of low-dispersion Runge-Kutta schemes, selectable from first- to fourth-order. Spatial discretization is performed 

using central differencing schemes. Both standard schemes, second- to twelfth-order (3 to 13 point stencils) and Dispersion Relation 

Preserving schemes from 7 to 13 point stencils are available. The code is written in Fortran 90 and uses hybrid MPI/OpenMP parallelization. 

The code is applied to the simulation of a Mach 0.9 jet flow. Four-stage third-order Runge-Kutta time stepping and the 13 point DRP spatial 

discretization scheme of Bogey and Bailly are used. The high resolution numerics used allows for the use of relatively sparse grids. Three 

levels of grid resolution are examined, 3.5, 6.5, and 9.2 million points. Mean flow, first-order turbulent statistics and turbulent spectra are 

reported. Good agreement with experimental data for mean flow and first-order turbulent statistics is shown. 

15. SUBJECT TERMS

Large-eddy simulation (LES); Turbulent flow; Jet flow 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

28 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18