Suresh Menon, Vaidyanathan Sankaran, and Christopher Stone

Georgia Institute of Technology, Atlanta, Georgia

Subgrid Combustion Modeling for the Next

Generation National Combustion Code

NASA/CR—2003-212202

April 2003

CCL–02–003

https://ntrs.nasa.gov/search.jsp?R=20030038812 2018-05-15T15:53:52+00:00Z
The NASA STI Program Office . . . 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 Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published 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’s counterpart of peerreviewed formal professional papers but

has 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. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, 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 Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
Suresh Menon, Vaidyanathan Sankaran, and Christopher Stone

Georgia Institute of Technology, Atlanta, Georgia

Subgrid Combustion Modeling for the Next

Generation National Combustion Code

NASA/CR—2003-212202

April 2003

National Aeronautics and

Space Administration

Glenn Research Center

Prepared under Grant NAG3–2660

CCL–02–003
Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

The Propulsion and Power Program at

NASA Glenn Research Center sponsored this work.

Available electronically at http://gltrs.grc.nasa.gov
NASA/CRó2003-212202 1 

Subgrid Combustion Modeling for the 

Next Generation National Combustion Code 

Suresh Menon, Vaidyanathan Sankaran, and Christopher Stone 

School of Aerospace Engineering 

Georgia Institute of Technology 

Atlanta, Georgia 30332 

1. Introduction

 Recent more stringent emission regulations have pushed for the development of 

more fuel-efficient and low-emission combustion systems for both ground-based and 

flight gas turbine applications. Experimental estimate of pollutant formation is limited 

due to the inability to probe non-intrusively, the highly three-dimensional flame zone. As 

a result, detailed understanding of where (and why) the pollutants (primarily, NOx, CO 

and unburned hydrocarbons, including soot) are forming is difficult, if not impossible to 

determine using experimental data. This is a major stumbling block in the design of the 

next generation, fuel-efficient gas turbine engines. Numerical prediction could help but 

most codes currently in use employ steady state methods and are unable to capture the 

unsteady dynamics that control the pollutant formation process. In addition to the 

inaccuracy in the prediction of NOx formation, CO emission is also not predicted very 

well under dynamic conditions. For example, recent studies (Bhargava, 2000) have 

shown that as the equivalence ratio is reduced, CO emission first decreases and then 

suddenly increases rapidly. This sudden increase occurs during lean combustion when 

perturbation to the flame can cause instability and also lead to incomplete combustion. 

The condition at which this sudden increase occurs is very important to predict and 

control since it impacts the performance and efficiency of the gas turbine engine. 

However, so far, no numerical simulation tool has been able to predict this behavior. 

 Accurate prediction of unsteady combustion requires a comprehensive model that 

can predict flame structure and propagation characteristics, pollutant formation and 

transport, and ignition/extinction phenomena over a wide range of flow conditions in 

high Reynolds (Re) number, three-dimensional (3D) flows as in a gas turbine combustor. 

Direct Numerical Simulations (DNS) are not practical since the resolution and 

computational resource requirements far exceed the present and even future capability. A 

simulation approach that has become popular in recent years is large-eddy simulation 

(LES). In LES, all scales larger than the grid resolution are numerically simulated using a 

space- and time-accurate scheme while the effect of scales below the grid resolution is 

modeled using a subgrid model. Although many LES studies have been reported, a 

validated LES approach for practical systems is yet to be demonstrated although the 

approaches discussed in this paper have shown a serious potential to achieve this 

capability. 
NASA/CRó2003-212202 2 

 Over the last few years, there has been a major effort at NASA-GRC to develop a 

unique computational capability called the Numerical Propulsion Simulation System 

(NPSS). In this approach, a simulation model that will allow simultaneous calculation of 

the flow in the entire engine (i.e., inlet, compressor, combustor, turbine and nozzle) is 

currently being pursued and developed. Within this structure, the simulation tool for the 

combustor, called the National Combustion Code (NCC) has reached a major state of 

maturity and is currently being evaluated for practical applications both at NASA and at 

industrial R&D sites (e.g., GEAEC). Various advanced models have been implemented 

to make NCC a state-of-the-art modeling tool and its application to various problems are 

being demonstrated and/or discussed in numerous recent papers (e.g., Liu, 2001, Raju, 

2001, Lannetti et al., 2001, Moder, 2001, Ajmani, 2001, Shih et al., 2001). 

 The NCC is based on Reynolds-Averaged Navier-Stokes (RANS) modeling 

approach and is the current choice since it is computationally efficient when complex 

systems have to be simulated. However, it has some inherent (and well-known) 

limitations. In RANS modeling, the effect of all turbulent length and time scales are 

modeled using global turbulence models such as the popular k-e model. Closure of the 

turbulent reaction rate and turbulent species flux are also based on global model that 

ignore some fundamental physics of turbulent mixing and combustion. For example, the 

dynamical interactions between flow, unsteady heat release, acoustic (pressure) 

oscillation and turbulence is lost since only the time-averaged flow field can be simulated 

with such an approach. These very important phenomena directly impact fuel-air mixing 

process and are controlled by processes that occur at the small-scales. Unfortunately, 

even unsteady RANS models cannot resolve these fine-scale interactions in a physically 

consistent manner since all these effects are modeled. 

 In this three-year research study, we plan to take advantage of two crucial 

developments: (a) the demonstrated capability of NCC to simulate flows in complex 

configurations using robust algorithms and (b) the recent demonstration of a unique 

subgrid combustion model for LES applications. We (the team from Georgia Tech and 

NASA/GRC) will combine these two capabilities into a new simulation model, denoted 

hereafter NCCLES that will contain all the capabilities of the original NCC but will now 

contain additional algorithms that will allow it to be used in a LES mode. It is envisioned 

that NCCLES will eventually become a specialized tool that will address issues that the 

original NCC cannot. 

2. Technical Objectives 

 The proposed research is carrying out both model development and model 

implementation simultaneously. The model development effort (which will primarily be 

at Georgia Tech) is addressing issues related to complex high-Re flows, including 

supersonic flows. Model implementation of the subgrid kinetic energy model and the 

subgrid mixing and combustion model within NCC is being carried out at NASA/GRC. 

Substantial progress in both the collaborative effort and in the new advances in the LES 

subgrid models has been achieved already. More details of the progress made to date are 

summarized below. 
NASA/CRó2003-212202 3 

The overall technical objectives of this three-year project are: 

1. Collaborate with the NCC team and provide algorithms that will allow them to 

convert the current NCC algorithm to a LES version.

It is worth noting that since many of the baseline algorithms are already 

developed, transition of this technology to the NCC team has already begun. More 

details of the current collaborative effort are summarized below. 

2. Further develop the subgrid models for momentum and scalar transport for 

implementation into NCCLES.

This is a major effort since NCCLES is expected to deal with many 

different scenarios and the subgrid models have to be able to simulate all such 

flows without compromising the physics. Therefore, the models currently 

operational need to be re-evaluated in this context. This effort has already begun 

and we have made substantial progress to-date. 

3. Collaborate with the NCC team to demonstrate and validate the next generation 

NCC for LES applications. 

All the models developed during this research will be implemented in our 

in-house code, LESLIE3D that has been extensively used for such studies. Once 

the validation of the new algorithms is completed, we will transition this 

technology in a timely manner to the NCCLES for further evaluation. Close 

collaboration with the NCC team will be maintained during this effort. 

3. The LESLIE3D Code 

 The LESLIE3D code is a fully compressible, 3D finite volume solver that has 

been used extensively in the past to carry out LES studies of complex turbulent reacting 

flows as in the gas turbine combustor. The code can be used for both DNS and LES and 

is second-order accurate in time and has an option to do either second or fourth-order 

accurate scheme in space. A constant coefficient model and a localized dynamic subgrid 

model (where the model coefficients are obtained as a part of the simulation and hence, in 

this case, the closure has no parameters to set or adjust) for the subgrid kinetic energy are 

options in the code. In addition, models for subgrid combustion are also implemented for 

both premixed and non-premixed combustion. For completeness, the LES model 

equations and closures are summarized here. 

3.1 The LES Equations 

The compressible LES equations are obtained by spatial filtering of the 

conservation equations of mass, motion, energy and species. The final set of LES 

equations are:
NASA/CRó2003-212202 4 

[ ]

[ ] ( )

m

sgs

i m

sgs

i m

i

m

m i m

i

m

sgs

ij

sgs

i i j ij i

i

sgs

i j ij ij ij

j

i

i

i

w

x

Y Y u D

t x

Y

E p u q u H

t x

E

u u p

t x

u

x

u

t

= ! 



 





+ Φ + Θ

∂

∂ − ∂

∂

+

∂

∂

+ + − + + = ∂

∂

+

∂

∂

+ − + = ∂

∂

+

∂

∂

= ∂

∂

+

∂

∂

, ,

~ ~ ~ ~

0 ~ ~ ~ ~

0 ~ ~ ~

0

~

ρ ρ ρ

ρ τ σ ρ

ρ δ τ τ ρ

ρ ρ

 (1) 

Here, the filtered total specific energy E

~ is given by: l l sgs k u u E = e + +

2

~ ~ ~ ~ , where 

sgs k is the subgrid turbulent kinetic energy given by ( ) 2 ~ ~ l l l lsgs k = u u − u u . The unclosed 

terms sgs

i m

sgs

i m

sgs

ij

sgs

i

sgs

ij H , , τ , ,σ ,Φ ,Θ represent the subgrid stress tensor, heat flux, viscous 

work, species mass flux and species diffusive mass flux, respectively, and are: 

( ) i j i j

sgs

ij u u u u ~ ~ τ = ρ −

( ) T j T j

sgs Hij H u H u

~ = ρ −

( ) i ji i ji

sgs

j σ = ρ u τ − u τ (2) 

( ) k j k j

sgs

k j Y u Y u ~ ~ Φ , = ρ −

( ) k k j k k j

sgs

k , j Y U , Y U ,

~ ~ Θ = ρ −

The pressure p is obtained from the equation of state for a thermally perfect gas, 

p = ρRT . The mixture gas constant is given by = ∑

m

R Ym Rm . A similar expression also 

applies for mixture specific heats Cp Cv and where Cp − Cv = R . 

3.2 Subgrid Closure for Convective Fluxes 

Closure of the subgrid flux terms is achieved using a transport equation for the 

subgrid kinetic energy per unit mass ( ) 



 



 = − 2

~

2 ~

2

1

k k

sgs k u u (Menon et al., 1996): 

( ) 















∂

∂

∂

∂ = − +

∂

∂

+

∂

∂

i

sgs

t

t

i

sgs sgs sgs

i

i

sgs

x

k

x

u k P D

t x

k

Pr

~ ν

ρ ρ ρ (3) 

Here, 

j

sgs i

ij

sgs

x

u P

∂

∂ = −

~

τ and = ( ) / ∆ 3 / 2 sgs sgs D C ρ k ε are the production and 

dissipation terms, respectively. Subgrid stress tensor is closed by an eddy viscosity 

hypothesis based on the subgrid kinetic energy. Thus, 

ij

sgs

t ij kk ij

sgs

ij τ ρν S S δ ρk δ

3

~ 2

3

~ 1 2  +



 



 = − − where the eddy viscosity is = ( ) ∆ 1/ 2 sgs

t C k ν ν . 
NASA/CRó2003-212202 5 

The two model coefficients, Cν and Cε , are 0.067 and 0.916, respectively for the 

constant coefficient model based on earlier studies (Chakravarthy and Menon, 2001a, 

2001b). A localized dynamic approach developed using a scale-similarity approach can 

also be used to obtain these coefficients model dynamically (Kim and Menon, 1999, Kim 

et al., 1999, Kim and Menon, 2000). Given the subgrid kinetic energy, the unclosed terms 

in LES equations can be modeled. For example, the sub-grid heat flux and the subgrid 

mass flux are: 

t i

sgs t

i

x

h H

∂

∂ = −

~

Pr

ν

ρ

i

m

t

sgs t

i m x

Y

Sc ∂

∂ Φ = −

~

,

ν

ρ

Where Sct and Prt are the turbulent Schmidt number and Prandtl number. At present, 

these are assumed to be unity. 

 LESLIE3D contains closure models for both premixed and non-premixed mixing 

and combustion process. These are briefly summarized below. 

3.3 Subgrid Flame Speed Model 

This approach is only applicable to premixed combustion in the flamelet regime. 

In this case, the flame is much thinner than the smallest scales of turbulence and 

therefore, turbulence can only wrinkle the flame but cannot affect the flame structure. 

Flame wrinkling and stretching increases the effective flame surface area and therefore, 

can increase the burning velocity. In the flamelet limit, a simplification can be carried out 

whereby the flame is modeled as an infinitely thin surface that separates the unburnt 

reactants from the burned products. In this approach, reaction kinetics need not be solved 

and the problem essentially reduces to the tracking of the flame front. This approach has 

been used quite successfully to simulate lean premixed combustion in full-scale gas 

turbine engine such as the General Electric LM6000. We will be using this version of the 

LESLIE3D initially to develop and evaluate the NCCLES. 

The model involves the transport solution of a progress variable G that lies in the 

range [0,1] with the flame located between the two limits. The G-field is convected by 

the flow field and also propagated normal to itself with the normal burning velocity. In 

the laminar case, this is the classical laminar burning velocity (which contains all the 

thermo-chemical effects of finite-rate kinetics) but in the turbulent field, a closure is 

needed to estimate the effective turbulent burning speed. The LES model for the flame 

front is obtained by the solution of the following filtered G-equation: 

) ~ ~ (

~ ~ ~

u G u G

x

S

x

u G

t

G

i i

i

sgs

i

i − ∂

∂ = − − ∂∂

+

∂

∂

ρ ρ ρ (4) 
NASA/CRó2003-212202 6 

Here, the terms of the right-hand-side require closure. Here, sgs S is the turbulent 

flame propagation term typically given as: S ST G sgs ~ = ∇ where ( ', ) ST = ST u SL is the 

turbulent flame speed that is modeled in terms of the laminar flame speed and the subgrid 

turbulence (obtained from subgrid kinetic energy: ( ) 



 



 = − 2

~

2 ~

2

1

k k

sgs k u u ). In the present 

study, a model proposed by Pocheau has been used to obtain the turbulent burning 

velocity. The second term ) ~ ~ G (ui

G uiG sgs

i = − is the subgrid convective flux that also 

requires closure. Various models have been proposed for this term: (a) localized dynamic 

eddy diffusion model (Kim and Menon, 2000), (b) flame curvature model (Peters, 2000), 

etc. We have evaluated all these models as reported below. 

Heat release is incorporated by including the heat of formation in the definition of 

the internal energy using the following relationship: ) ~( ~ H G G0 e C T e = v + f − where, 

H(Q) is the Heavy-side function and G0 is an arbitrarily chosen value in the range [0,1] 

where the flame is located. More details are given elsewhere (Kim et al., 1999, Kim and 

Menon, 2000) and therefore, avoided here, for brevity. Although this approach is 

computationally very efficient (since species equations and/or finite-rate kinetics do not 

have to be modeled), this approach cannot be used to capture species variations and/or be 

used to predict pollutant formation. To do this, one has to resort to the more complex 

finite-rate kinetics closure. 

3.4 Subgrid Scalar Closure using Eddy Diffusivity 

For both premixed and non-premixed combustion, if detailed kinetics is of interest 

multi-species equations have to be solved along with the LES equations. In this case, two 

terms, the sub-grid scalar flux, ] ~ ~ [ ,

~

, i,m m i m m

sgs

i m = V Y −V Y

$#"

θ ρ and the filtered reaction rate 

term, ω m! require closure. In a conventional approach, a subgrid gradient approximation: 

m

sgs t m

i m Sc

Y

~

,

∇ ≈ − ρν θ where Scm is a turbulent Schmidt number, is employed for the subgrid 

scalar flux. Note that, since large-scale motion is resolved in a LES, associated countergradient processes are resolved (even when a gradient closure is employed for sgs θ i,m ). The 

closure for ω m! s complicated due to its highly non-linear nature. Assumed PDF subgrid 

models, subgrid Eddy Breakup and subgrid Eddy Dissipation models, all of which are 

variants of the models used in RANS codes, have been proposed for this closure. 

 

However, at this time, the accurate closures for the reaction rate within the 

conventional LES context have yet to be established. In order to ensure accuracy, we 

have developed a new closure model based on a subgrid simulation approach. This is 

described in the next section. 
NASA/CRó2003-212202 7 

3.5 Subgrid Scalar Closure using LEM 

In the LEM closure, the species fields evolve within each LES cell due to 

localized, stochastic reaction-diffusion processes, and are then transported across the LES 

cells due to convective flux. The local sub-grid domain within each LES cell is resolved 

on a one-dimensional domain that resolves all scales of motion and thus, processes within 

this domain can be considered a localized 1-D Direct Numerical Simulation. Thus, all 

scales down to the Kolmogorov and Batchelor's scale are resolved in this simulation 

(Menon et al., 1993, Calhoon and Menon, 1996, Chakravarthy and Menon, 2001a, 

2001b). 

To describe this method, consider an ìexactî reaction-diffusion equation for a 

scalar φ: 

φ ω φ

φ φ φ + ! ∂

∂

∂

∂ = ∂

′ ∂

+ +

∂

∂ ) [ ] ~(

x

D

x x

u u

t i i

i i (5) 

Here, ui

~ and ui

′ are the resolved and unresolved velocities, respectively. Also, 

ui

′ is due to both the Lagrangian convection through the cell face ulag ′ , and turbulent 

convection at scales smaller than the resolved grid, ustir ′ 

The one dimensional linear-eddy domain, in every LES cell is initialized with a 

constant number of LEM cells. The number of linear-eddy cells that are needed to resolve 

all the sub-grid scales is computed based on a subgrid Reynolds number probability 

density, obtained over the whole of the computational domain. The length of the linear 

eddy domain is set equal to the LES filter width. Since the Lagrangian convection across 

the LES cell face is based on the mass, the volume of the linear-eddy cell is also carried 

as an additional scalar. At the start of the simulation, volume of each linear eddy cell is 

set to have an equal fraction of the LES cell volume. At each LES step, after updating the 

density, momentum and the energy field, LEM is used to evolve the scalar field. Note 

that, as opposed to the conventional scalar closure approach described in the previous 

section, no scalar equations are solved on the resolved grid. 

Using a fractional splitting technique, we can split the equation (5) as follows: 

Subgrid Processes: dt

x

D

x x

u

i i

t t

t

stir

LES

O

[ [ ] ] φ ω φ

φ φ + ! ∂

∂

∂

∂

+

∂

∂ ′ ∫

+ ∆

 (6) 

Supergrid Convection: 

i

i lag

LES

n

x

u u

t ∂

∂ = + ′ ∆

φ − φ φ ) ~(

*

 (7) 

Equation (6) represents the sub-grid reaction-diffusion processes that occur 

locally within each LES cell using the LEM model. Three processes occur within each 

LES cell: (a) turbulent convection at scales smaller than the resolved grid (which is 

modeled using stochastic stirring events called triplet maps, (b) molecular diffusion and 

finite-rate kinetics, and (c) volumetric heat release. Additional details are given in cited 
NASA/CRó2003-212202 8 

references. After the first two subgrid processes, volumetric gas expansion caused by 

heat release is implemented by expanding each linear-eddy cell in the domain by an 

amount equal to 1

1

, +

+ ∆ = n

i

n

n i VLEM i

ρ

ρ . 

Equation (7), which represents the advection (convective transport) of the scalar 

field by the resolved velocity field is modeled by a volume of fluid (VOF) approach 

(Chakravarthy and Menon, 2001a, 2001b) in which the linear eddy cells are convected 

across LES cell face at the advective time scale given by ( ) ( )min

i sgs

LES

conv u u

x

t CFL

+ ′

∆ ∆ = . 

Ideally, splicing in a compressible simulation can be performed every conv ∆t step. But, in 

order to provide a close coupling between the fluid dynamic and the flame front 

evolution, splicing is performed every LES time step. This considerably improves the 

numerical stability of the simulation. 

On a general three-dimensional grid the mass transfer from any control volume to 

any neighboring control volume is also predominantly 3D. As a result, different numbers 

of linear-eddy cells are transported in different spatial directions. In the numerical 

implementation of the scalar advection, the 3D advection operator is approximated by a 

sequence of three, one dimensional advection operators. So the order in which these 

operators act on the scalar field can have a significant effect on the scalar field evolution. 

In the new algorithm, the linear-eddy cells going out in the direction of largest out-flux 

are fluxed out first from the right end of the sub-grid domain. Similarly the largest influx 

is added first to the left end of linear eddy domain. This upwind type flux update has been 

used successfully in the past (Chakravarthy and Menon, 2001a, 2001b). 

Figure 1. Schematic illustrating the Splicing Algorithm used for scalar convection. 

3.5.1 New Splicing Algorithm Development under this Project

The earlier studies were limited to low-Mach number flows and therefore, some 

issues were substantially simplified and easy to implement. For example, in low-speed 

flows, the CFL time step is the convective time step, which is also the same time needed 

for the scalar transport via splicing. However, in the present study, we are eventually 

Arrows - Convection Direction 

Color Lines - Flame 

Blue Cell - LES Cells 

Black Cell - LEM Cells
NASA/CRó2003-212202 9 

interested in extending this approach to fully compressible flows, including supersonic 

flows. In this case, the CFL time-step is the acoustic time step whereas the scalar 

transport is still restricted by the convective time-step, which is at least an order of 

magnitude larger than the acoustic time-step. Thus, there is a possibility of de-coupling 

between the fluid dynamics (which evolves at the acoustic time) and the scalar transport 

(that occurs at the convective time). A substantial effort has been spent on understanding 

this issue and we have now developed a new convective transport model applicable for 

compressible flows. 

 The present splicing algorithm advects the sub-grid scalar field based on the 

mass, as opposed to the volume of the LEM cell that was used in the past. The advantage 

of using mass instead of volume is that, it ensures mass conservation at the resolved and 

unresolved scales at each instant of the calculation. Also in the earlier approach, partial 

LEM cells that have to be convected are accumulated, until they reach an integer value. 

In the new algorithm, every time a partial cell has to be convected out, a single cell of the 

same volume and mass, as that of the partial cell and with all the scalar properties is 

convected out. This potentially avoids the mass conservation errors and spurious 

numerical instabilities in the computation. 

Volumetric expansion due to heat release and splicing by itself, also lead to a nonuniform distribution of the volume of the linear-eddy cells. To avoid programming 

complexities in a parallel environment, it is desirable to have same number of linear eddy 

cells everywhere in the computational domain. So, to equalize the volume of the lineareddy cells and to retain the same number of linear eddy cells everywhere in the 

computational domain, the linear eddy domain is re-gridded to have cells of equal volume 

after each splicing. Errors introduced due to the spurious diffusion associated with regridding have not been evaluated yet. But it is expected to be minor, since at most only a 

couple of cells are affected. 

Furthermore, to reduce any biasing error introduced by the operator splitting, the 

order in which reaction, diffusion and stirring are performed is switched after every 

splicing. An additional advantage of this switching is that, any discontinuities introduced 

by splicing are smoothed quickly by turbulent stirring, thus preventing any spurious 

diffusion. 

3.5.2 Subgrid Simulation using LEM

The details of the sub-grid one-dimensional simulation using LEM are briefly 

summarized here. LEM is a stochastic model that treats reaction-diffusion and turbulent 

convection separately but concurrently. Reaction-diffusion processes evolve on a onedimensional (1D) domain in which all the characteristic length scales in the turbulent 

field (from the integral scale L to the Kolmogorov η) are fully resolved (6 cells are used 

to resolve η). The orientation of the 1D domain is in the direction of the scalar gradient 

(and thus, for premixed flame, is in the flame normal direction), and within this domain, 

the equations governing constant pressure and adiabatic flame propagation are: 

kstir

k k k k k F W

x

Y V

t

Y

+ +

∂

∂ = − ∂

∂

ρ

ρ ω

ρ

1 ! (8) 
NASA/CRó2003-212202 10

Tstir

N

k

k k k

p

N

k p

p k k k

p

h W F

x C

T

x C x

T C Y V

t C

T ∑ ∑ = =

− +

∂

∂

∂

∂

+

∂

∂ = − ∂∂

1 1

,

1 ( ) 1 1

ω

ρ

κ

ρ ρ ! (9) 

The equation of state for the scalar mixture is ∑=

=

N

k k

k u

W

Y R P T

1

ρ and the caloric relation is 

given by ∫ = ∆ +

T

T

hk hf k Cp k T dT

0 ,

0

, ( ') ' . Here, T is the temperature, P is the thermodynamic 

pressure, Ru is the universal gas constant and ρ is the mass density. Also, 

k k p k k k Vk Y ,W ,C , ,h , , ω! and 0 ∆hf ,k are respectively, mass fraction, molecular weight, specific 

heat at constant pressure, mass reaction rate, enthalpy, diffusion velocity and standard 

heat of formation (at standard temperature, To

) of the k-th species. The mixture-averaged 

specific heat at constant pressure and thermal conductivity are respectively, Cp and κ . 

Fickian diffusion velocity law: dx

dY

Y

D V k

k

k

k = − , where Dk is the mixture-averaged 

diffusivity of the k-th species, is used for molecular diffusion. 

The convective terms i k i u ∂Y ∂x and i i u ∂T ∂x in the species and the energy 

equations are represented as Fkstir and FTstir, respectively. These terms are implemented 

using stochastic re-arrangement events called triplet maps, each of which represents the 

action of a turbulent eddy on the scalar fields. It has been shown that this mapping can 

capture correctly the physical increase in scalar gradient (without affecting the mean 

scalar concentration) due to eddy motion. Three parameters are needed to implement 

these turbulent stirring events: the typical eddy size l, the eddy location within the 1D 

domain and the stirring frequency (event rate) λ. The eddy size in the range L to η is 

determined randomly from an eddy size distribution, f(l) which is obtained using inertial 

range scaling in three-dimensional turbulence: f(l) = (5/3) l-8/3 / (η-5/3 - L-5/3). Here, η is 

determined from inertial range scaling law η = Nη L Re-3/4 where Nη is an empirical 

constant. This constant reduces the effective range of scales between L and η but does not 

change the turbulent diffusivity, as described earlier. 

The event location is randomly chosen from a uniform distribution and the event 

rate (frequency per unit length) is determined by: 

[1 ( / ) ]

Re [( / ) 1]

5

54

4 / 3

5 / 3

3 L

L

C L η

ν η λ

λ −

− = (10) 

The time interval between events is then given as ∆tstir = 1/(λXLEM) where XLEM is the 

length of the 1D domain. 

The unique feature of the LEM model is that although the scalar evolution is 

simulated in 1D, the effect of turbulence on the scalar fields is modeled using 3D scaling 

laws. As a result, flame wrinkling occurs at spatial and temporal scales that mimics the 

effect of realistic 3D turbulent eddies on the flamelets. This formulation has two 

constants: Cλ and Nη both of which arise from the use of scaling laws. These parameters 

are obtained by comparing LEM predictions to experimental data in the flamelet regime 

(Smith and Menon, 1996). The present study uses these same values (Cλ = 15 and 

Nη = 4). 
NASA/CRó2003-212202 11

4. Progress to Date 

 

Here, we summarize the progress made to date. It is worth noting that the project 

officially started at Georgia Tech in late January 2002 and therefore, the progress 

noted below is for a period of 9 months. We report here primarily on the effort carried 

out at Georgia Tech over this period. There were two major areas of progress: (a) 

transition of LES technology to NCC team at NASA/GRC and (b) further development of 

LES model for more complex flows. 

 4.1 Collaboration with NCC Team at NASA/GRC 

 Discussions with the NASA/GRC team members were carried out to set up the 

work plan. It was determined that one area that needs to be explored is the manner by 

which LES is carried out since it differs substantially from RANS application. It was also 

decided that this “learning” process needed to be carried out prior to implementing the 

LES routines into NCC. In order to do this Georgia Tech provided the NCC team with 

our advanced LES code, LESLIE3D (Large-Eddy Simulations using Linear Eddy in 

3D). This code is the primary research tool employed at Georgia Tech to develop all the 

algorithms to be eventually implemented in NCC. Close collaboration has been 

established in order to “teach” the method for doing LES using our code. 

4.1.1 LES Simulation strategy using LESLIE3D 

 In order to establish an understanding of how to do LES, the first phase of 

collaboration was to use the LESLIE3D to simulate a test case that has also been 

simulated at Georgia Tech. The NCC researchers have been using the LESLIE3D code to 

carry out the baseline simulation of the GE LM6000. This test case was also simulated at 

Georgia Tech and the detailed analysis package was also provided to NASA/GRC. 

 The test case is a full-scale GE LM6000 that was experimentally studied at 

GEAEC. Figure 2 shows the test configuration. 

 

Figure 2. Schematic of the GE LM6000. 
NASA/CRó2003-212202 12

Premixed methane-air combustion is simulated under highly swirling conditions. 

The Reynolds number based on the inlet mean velocity and diameter is the 350,000, the 

inlet Swirl Number at the dump plane is 0.56, the inlet temperature is 688 Kelvin, the 

combustor mean pressure is approximately 6.5 ATM. We have employed relatively 

coarse grid of around 600,000 in the past to simulate this problem using the localized 

dynamic subgrid kinetic energy model and the dynamic flame speed model. A typical 

result is given in Fig. 3 which shows the centerline velocity decay. Various closure 

models for the subgrid diffusion term ) ~ ~ G (uiG uiG sgs = − were attempted but the results 

are relatively insensitive to them and all cases show reasonable agreement with data. 

Additional comparison of the velocity profiles also show good agreement. 

Figure 3. Centerline mean axial velocity variation. Symbols are the experimental 

data from GEAEC. The various lines correspond to the subgrid closure for the 

convective flux term ) ~ ~ G (uiG uiG sgs = − . All models show reasonable agreement 

with data. 

4.1.2 DNS of Decaying Isotropic Turbulence 

 Another area where collaboration has been established is in the area of inflow 

turbulence and initialization process. Typically, isotropic turbulence field is added to the 

mean flow. A version of the LESLIE3D was also developed and provided to NASA/GRC 

that can be used to generate isotropic turbulence. Since this solver has the ability to do 

both DNS and LES it was decided to simulate this sort of flow using both the second and 

fourth-order spatial schemes so that a baseline data base can be established to evaluate 

the ability of the NCC base scheme to simulate isotropic turbulence. 

One case chosen for this purpose is direct numerical simulation (DNS) of 

decaying isotropic turbulence. Since the physics of the isotropic flows (such as the 

energy decay rate, form of the spectrum, etc.) are well established from the numerous 

experimental and numerical studies from the past, this exercise can serve as a powerful 

test-bed for establishing the capabilities of the code. 
NASA/CRó2003-212202 13

Decaying isotropic turbulent flow experiments of Comte-Bellot and Corrsin 

(CBC) (1971) were chosen for validating the compressible solver. In the experiments, 

measurements were taken at three locations downstream of the turbulence-generating 

grid. The mesh size of the grid was 5.08 cm and the mean velocity ahead of the grid was 

10 m/s, giving a mesh Reynolds number of Re = 34000. Downstream of the grid, the 

wind tunnel had a slight contraction (1.27:1) that served to provide an isotropic condition. 

Conceptually, homogenous isotropic turbulence decays in time. However, in the wind 

tunnel environment, the turbulence decays in space as it evolves downstream. In order to 

convert a spatial evolution to the temporal evolution, Taylorís hypothesis is used and the 

turbulence in space is related to the decay in time of an imaginary box of homogenous 

turbulence through the relation: ∫ =

X

U x

dx

t

0 ( )

 , where U ( x) is the mean velocity of the 

flow at the streamwise location x. 

To design a computation for the simulation of the CBC experiment, appropriate 

initial conditions that has the same spectrum as in the experiment is needed. In general 

the initial velocity field for the isotropic turbulence must satisfy at least three conditions: 

(a) conservation of mass, (b) the generated velocity field must be a real function of space 

and time, and (c) the turbulent field must have a realistic energy spectrum. 

The first two conditions are required in order to obtain a stable numerical solution 

while the third conditions is needed to reduce the initial transient as the velocity field 

evolves to form realistic isotropic turbulence. In reality, not only does one need to 

specify the energy spectrum, but also the relative phases of the modes. In practice, 

however, this phase information is not available so that fields are usually constructed 

with random phase, but in such a way that they have a prescribed initial energy spectrum. 

Therefore, this initial condition will lead to a transient in which the phases adjust 

themselves to appropriate values. 

The governing conservation equations for compressible flows are solved using a 

finite volume scheme that is second order in time and second or fourth order in space. In 

general, for a given Reynolds number, the second-order scheme requires twice the 

resolution of the fourth-order scheme to maintain DNS and/or LES accuracy. This is an 

important issue for model and code validation. The results reported below are for DNS 

and LES using the fourth-order scheme. The computational grid used is 64 x 64 x 64 with 

uniform spacing (for DNS) and periodic conditions are imposed on all the boundaries. 

The parameters corresponding to the simulation and the experiments are given in Table 1. 

Table 1 

Station

L

Ut Re λ Simulation Grid 

1 42 71.6 LES 32 x 32 x 32 

2 98 65.3 DNS 64 x 64 x 64 

3 171 60.3 DNS 64 x 64 x 64 
NASA/CRó2003-212202 14

The DNS simulations are initialized with a velocity field corresponding to the 

Re λ=65.3 (station 2, L

Ut

=98) and is evolved for a time corresponding to the equivalent 

spatial distance between station 2 and station 3. Figure 4 shows the energy spectrum 

along with the experimental spectrum at two of the locations. It can be observed the 

present solver captures the energy cascade in the inertial range very well. The most 

significant aspect of this figure is that dissipation scales are resolved accurately. This is 

very important for any solver to demonstrate. If the numerical dissipation of the scheme 

is too high, dissipation range cannot be resolved very well. On the other hand, if the 

dissipation is not sufficient, the code becomes numerically unstable. The fact that the 

Kolmogorov cascade and the dissipation range of the spectrum are captured accurately is 

a proof that the code provides the correct amount of dissipation, while still retaining the 

stability of the solution. The deviation in the largest scales is because of the smaller size 

of the physical box used for the simulation. Since a bigger box size would require a 

higher resolution, slightly smaller box is used in the current simulation. But, as it is 

pointed out earlier, it is the smaller scales that are usually difficult to capture and hence, 

this should not be of any concern. 

Figure 5 shows the vortical structures at the initial time and at a time when the 

turbulence has fully evolved. The decay of the turbulence and the growth of the small 

scales is clearly apparent in this figure. 

4.1.3 LES of Decaying Isotropic Turbulence: 

Having established the accuracy of the numerical scheme using the DNS, it 

remains to establish the accuracy of the sub-grid models employed in the present 

simulation. For this purpose, a large eddy simulation with a higher Taylorís Reynolds 

number ( Re λ = 71.6) was carried out. The initial velocity field was initialized with a 

velocity field that would provide the correct resolved and subgrid energy (as that in the 

experiments) corresponding to the computational grid (32 x 32 x 32) employed and the 

Re λ . A one-equation sub-grid kinetic energy along with an eddy-viscosity type closure is 

used to model the sub-grid scale processes. The rationale behind this and the validity of 

this approach has been very well established and documented. 

Figure 6 shows the decay of the turbulent energy decay rate. As can be seen the 

current model accurately captures the energy decay rate from the large scales to the small 

scales. Sub-grid models are designed to provide dissipation from the resolved to the subgrid scales and it is apparent that the current subgrid model provides the right amount of 

dissipation since otherwise it would not be possible to capture the correct decay rate 

Figure 7 shows the energy spectrum from the present LES and the CBC 

experiment. Inertial range is captured very well despite the coarse grid used in the 

simulation. As expected, the energy spectrum deviates from the experiments at the small 

scales, since the dissipation range scales are not resolved explicitly in any LES. Finally, 

Figure 8 shows the evolution of the vorticity field from the current LES at the initial time 

and at a later time. 
NASA/CRó2003-212202 15

Figure 4. Energy Spectrum for Re λ = 65.3 corresponding to CBC experiment.

Figure 5. Vorticity Field corresponding to Re λ = 65.3 at t = 0 and t = 0.2 seconds. 
NASA/CRó2003-212202 16

 

Figure 6. Kinetic Energy Decay in time. 

Figure 7. Turbulent Kinetic Energy spectrum for Re = 71.6 λ . 

Figure 8. Vorticity Field in LES for Re = 71.6 λ at initial time and at a later time. 
NASA/CRó2003-212202 17

 4.2 Research Development at Georgia Tech 

 In addition to the collaboration with the NCC team reported in Section 4.1, 

research was undertaken at Georgia Tech to further develop the LESLIE3D code. Four 

areas have been under focus. These are (a) Flame-turbulence test problem for NCCLES 

validation, (b) LES of supersonic spatially evolving shear layers, (c) LES of combustion 

in the Trapped Vortex Combustor (TVC), and (d) DNS and LES of two-phase mixing 

layers. The reasons for the choice of these test cases are noted in the following sections. 

 4.2.1 Flame-Turbulence Test Problem for NCCLES Validation 

 The interaction of a premixed flame with turbulence is a well-established test 

problem used in many past studies to understand the impact of turbulence on flame 

dynamics. This problem is particularly suited for code validation since both DNS and 

LES have been carried out in the past (including detailed kinetics) and therefore, the 

physics of this problem is well established. This test case is also relatively simple to 

implement and test. This test problem can also be used for evaluating the algorithmís 

accuracy in reacting flows. Some results are reported below on the current effort. 

DNS and LES of the turbulent premixed stoichiometric methane-air flames with 

finite-rate effects have been performed. To reduce the cost of computations incurred due 

to the high resolution DNS and by the inclusion of finite rate chemistry, the simulations 

are restricted to two-dimensions. However, this is not a problem in this study since 

similar studies have been performed in the past (Echkeki and Chen, 1991) to understand 

flame-turbulence interactions. 

The governing conservation equations are solved using a finite-volume scheme 

that is second-order accurate in time and fourth-order in space. The DNS was initialized 

in a two-dimensional computational domain with a plane laminar premixed flame. The 

reactant composition is a stoichiometric methane-air mixture at atmospheric pressure and 

temperature. The initialization for the scalars (species & Temperature) is obtained from a 

one-dimensional steady unstrained computations using PREMIX. A single step global 

chemistry (Westbrook and Dryer, 1971) with 5 species is used. The boundary conditions 

are periodic in the transverse direction and non-reflecting at the inflow and outflow 

(Poinsot and Lele, 1991). A computational grid of 400 x 400 with uniform spacing is 

employed for DNS and 100 x 100 for LES. 

The initial turbulent field is generated from a specified energy spectrum, 

satisfying continuity. The spectrum is of the form:































 = ′ −

2

5

4

2 exp 2 2 ( ) 16

o o

o k

k

k

k E k u

π , 

where k is the wave number, o k (= 6) is wave number of the most energy containing 

scale and uo

′ is turbulent intensity which is equal to 10.57 m/s in the present simulation. 

The turbulence velocity field is super-imposed on the laminar flame. The ratio of the 

turbulence intensity uo

′ to the laminar flame speed, S L is 4.2 and the ratio of the integral 

scale to the laminar flame thickness, is 3.4. The turbulence Reynoldsí number based on 

the integral scale is 137 and the Taylorís micro-scale is 36. 
NASA/CRó2003-212202 18

Three sets of simulations were carried out to contrast and compare the 

performance of the models employed in the present computations. DNS, LES-LEM and 

another LES using Eddy-Break-up closure (EBU) for scalars were performed. LES-LEM 

and LES-EBU simulations typically employed a uniformly spaced 100 X 100 

computational grid. LES-LEM used 18 LEM cells per LES cell. 

Figure 9 shows the instantaneous temperature contours computed using the LESEBU, DNS and LES-LEM simulations, respectively. The flame topology indicates the 

presence of highly wrinkled flamelet due to flame-turbulence interactions. Note that LESEBU diffuses the flame zone considerably, whereas LES-LEM actually makes the flame 

sharper when compared to the DNS. This is not surprising since LES-LEM actually 

captures the flame within the subgrid with very high resolution (actually better than 

DNS). For example, where as the DNS is capturing the flame over approximately 15 grid 

points, LES-LEM is capturing the flame over approximately 3-4 LES grid points. 

However, since there are 18 LEM cells per LES cell, the actual resolution is over 

50 subgrid points. Thus, LES-LEM captures the flame very accurately. This is similar to 

what was observed earlier (Chakravarthy and Menon, 2001a, 2001b). This is 

demonstrated in Fig. 10 where the instantaneous vorticity and density contours for the 

three cases simulated are compared. (Contours not shown at the same time. Note, this is 

only a qualitative comparison). Vigorous interaction between the local vorticity and the 

flame are apparent from this figure. 

 (a) Eddy-Break-up (b) DNS (c) LES - LEM 

Figure 9. Instantaneous temperature contours during flame-turbulence interaction. DNS is 

using 400x400 where as the LES is using 100x100. LES-LEM captures the flame as a sharp 

front (even better than DNS) since the subgrid resolution is better than in the conventional 

DNS. LES-EBU on the other hand diffuses the flame and fails to resolve the wrinkles. 
NASA/CRó2003-212202 19

 

4.2.2 LES of Supersonic Mixing Layers 

One of the eventual goals of this project is to extend the NCCLES capability to 

supersonic combustion applications. However, in order to achieve this capability, the 

LES-LEM approach has to be extended to highly compressible flows. 

The goal of the present study in the first year is to extend LES formulation to the 

highly compressible flows and to investigate the effects of the compressibility on the 

flows of engineering interest. In the first stage of this study the dynamic subgrid scale 

models are formulated with explicit correction terms, which accounts for compressibility 

effects. Simulations are underway for spatially evolving supersonic mixing layer for 

which there is some experimental data for model validation. Also, this particular flow 

was investigated earlier (Nelson and Menon, 1998) using an earlier version of the 

LESLIE3D code. 

LES Equations for Compressible Turbulence: 

The governing equations are obtained from Navier-Stokes equations by using 

standard filtering instead of Favre filtering as done in the previous studies. While it 

results in somewhat more complicated form of the governing equations than using Favre 

filtering, it provides for explicit terms that can be identified as “compressible” terms. 

Therefore, the modeling of these terms can be considered independently of the modeling 

issues typically used for the closure of the subgrid shear stresses. 

(a) Eddy-Break-up (b) DNS (c) LES - LEM 

Figure 10. Instantaneous density and vorticity contours during flame-turbulence 

interactions. Comparison is not shown at the same time. LES-LEM approach captures the 

flame very sharply compared to LES-EBU, even though same number of grid points are 

used in both the approaches, 
NASA/CRó2003-212202 20

The filtered LES equations can be written as 

( ) i i

i i

i u u

x x

u

t ρ ρ ρ ρ − ∂

∂ − ∂

∂ = − ∂

∂ (11a) 

( ) ( ) ( ) ij ij i i

j

i j i j

ij j

ij

i i

i j i u u

t

t

x

u u u u

x x

t

x x

u u

t

u

ρ ρ τ ρ ρ ρ ρ ρ − ∂

∂ − − ∂∂ − +

∂

∂ − ∂

∂

+

∂

∂ − ∂

∂ = − ∂

∂ 

 (11b) 

( ) ( ) 2

1

( ) 1 ( ) ( ) ( )

( ) ( ) ( ) ( )

sgs

i i i i

i i i i i

j ij j ij

i

j ij

i

i i

i

i i

i

i

i

k

t

u u u u

t

T T

t

R

x

T k

x

T k

x x

T k

x

u u t

x

u t

x

u u

x

u E u E

x

E u

t x

E

ρ ρ ρ

ρ ρ

γ

τ

ρ ρ ρ ρ ρ ρ ρ

∂

∂ − +

∂

∂ −

− − ∂

∂

− − ∂

∂ − ∂∂

∂

∂

+

∂

∂

∂

∂ − +

∂

∂

+

+

∂

∂ − +

∂

∂ − − ∂

∂

+ − ∂

∂ = − ∂

∂

 

 (11c) 

The resolved viscous stress tensor can be written as ij

k

k

i

j

j

i

ij x

u

x

u

x

u

t μ δ

∂

∂ − ∂

∂

+

∂

∂ = 3

1 ( ) 21 2 (

The total energy is sgs

i i E = e + u u + k

2

1 , where subgrid kinetic energy is defined as 

( ) 2

1

i i i i

sgs k = u u − u u . The filtered equation of state ρ = ρRT + R(ρT − ρT) completes 

the set of the governing equations. 

Compressibility correction terms: 

The above LES equations contain additional terms that require closure. The 

density-velocity correlation terms in the continuity equation is purely a compressible 

term and this term is rewritten as the difference between Favre and conventionally 

filtered velocity: ) ~ c ui ui (u u sgs

i = ρ − ρ = ρ − . This term would vanish in the inviscid 

limit and it is generally negligible in the region of mild gradients (Chen et al., 1989). In 

the region of stronger gradients, it can be argued that the contribution of this term would 

be proportional to the mean density gradient, which lead to a gradient diffusion model for 

this term:

i

i i c

sgs

i

x

c u u

∂

∂ = − ≈ − ρ ρ ρ ν where νc

 is a “compressible viscosity’’ that must be 

determined. 

Since this term is expected to be significant only near strong density gradients, it 

is advantageous to formulate it in such a way that it would include a density gradient 

switch, avoiding therefore the excessive dissipation in the smooth mean flow regions. 

From dimensional consideration this term can be formulated as: ν c = α c S p | uknk | ∆ as 

“compressible viscosity correction term,” where 
NASA/CRó2003-212202 21





















≤

> −

−

=

∆

∆

∆

min

min

1 min

0

0 Re Re

) Re Re

(Re Re )

1 exp(

c K

α

α (12) 

Here, 

t

uknk

ν +ν

∆

∆ = | | Re , where ∆ the local grid spacing is the characteristic length scale 

and velocity normal to cell face is the velocity scale. The expression for the scaling 

coefficient α was obtained from the numerical simulations of 1-D non-linear Burger’s 

equation (Nelson and Menon, 1998). 

The turbulent eddy viscosity, used in the cell Reynolds number definition is 

defined as ≈ ∆ sgs

t c k ν ν (following Schumann, 1976) and used for modeling of 

‘incompressible’ terms of subgrid scale closure equations. 

The LES continuity equation can now be rewritten as 

( )

i

c

i i

i

x x x

u

t ∂

∂

∂

∂

+

∂

∂ = − ∂

∂ ρ ν ρ ρ (13) 

The subgrid term in the momentum equation can be rewritten as i j i j

sgs c τ ij = ρu u − ρu u ( ) . 

Again, as in continuity equation, compressibility effects are significant only near the 

regions of high gradients, and again the compressibility effects model should be triggered 

by the presence of high gradients (e.g., shocks). Since shocks were observed to have a 

dissipative effect on turbulence, the compressibility effect is modeled as an additional 

dissipation term ) 3

1 2 ( ( )

c ij kk ij

sgs c Mij ≈ − ρν S − S δ where resolved rate of strain tensor is 

( ) 2

1

i

j

j

i

ij x

u

x

u S

∂

∂

+

∂

∂ = . Thus, the modeled momentum equation can then be written as 

j

sgs c

ij

j

sgs i

ij

j i

ij

i

i j i

x

M

x x x

t

x

u u

t

u

∂

∂ − ∂

∂ − ∂

∂ − ∂

∂

+

∂

∂ = − ∂

∂ ( ) ( ) ρ ρ ρ τ

 (14) 

The modeled subgrid kinetic energy equation in the compressible case is written as: 

(15) 

This subgrid model contains two model coefficients which can be dynamically 

computed as before (Nelson and Menon, 1998). Finally, the compressibility effects in the 

energy subgrid model closure equation can be estimated as 

) ) Pr ( ) ( ) ( ) ((

i

v c

i

j ij

i i

sgs

e

i

i

i x

c T k

x

u t

x x

H

c k

x

E u

t x

E

∂

∂

+

∂

∂

+

∂

∂

+

∂

∂ ∆

∂

∂

+ +

∂

∂ = − ∂∂ γ ρν

ρ ρ ρ ρ

(16) 

∆ − ∂

∂ − ∂

∂

+

∂

∂ − ∂

∂ = − ∂

∂ 2

3

( ) ( ( ) )

sgs

i

sgs j

ij

i

sgs

t c

i i

i sgs

i

sgs k

c

x

u

x

k

x x

k u

t

k

ρ ε ρ ν ν τ ρ ρ
NASA/CRó2003-212202 22

This equation introduces yet another subgrid scale modeling coefficient, which 

can also be dynamically modeled as well. 

Supersonic Spatial Mixing Layers: For validation of the numerical scheme and 

compressibility correction terms the comparison with Samimy et al (1991) experimental 

data was conducted. The original experiments were conducted in a high Reynolds 

number aerodynamic tunnel. Two stream mixing layers at different Mach numbers was 

separated at the entrance to the test region by a splitter plate. The flow scheme is 

provided in Figure 11. The optical test section IS 500x80 mm and velocity measurement 

were conducted using two-component laser Doppler velocimetry. 

Figure 11. Spatially evolving supersonic mixing layer setup of Samimy et al. (1991). 

The incoming flow parameters are provided in the Table 1. 

T (K) P1 (kPa) M1 M2 Mc U1

(m/s) 

U2/U1 ρ2/ρ1 δ (mm) Re 

291.0 314.0 1.80 0.51 0.52 479.5 0.355 0.638 8.0 9.4x106 

The computational area starts immediately after the edge of the splitter plate and 

has a length of 400mm and height 152 mm. Both upper and lower walls are slip walls. 

The periodic boundary is used in the z-direction. The inflow is fully specified on the 

supersonic part of the flow while on the subsonic part a characteristic boundary condition 

is used. The outflow boundary conditions are calculated using characteristic form of the 

governing equations using Poinsot and Lele method (1991). The initial conditions for the 

incoming boundary layers are computed using a 2-D boundary layer solver based on the 

known experimental data. For better comparison with the experiment a pseudo-turbulent 

velocity fluctuation are added to the inflow boundary. The turbulence with the assumed 

spectrum is calculated in a box and convected into the inflow boundary, scaled according 

to boundary layer solution on inflow boundary. 

Figures 12 and 13 illustrates the instantaneous density, pressure and spanwise vorticity 

contours iso-surfaces for the developed mixing layer obtained in a recent simulation 

using the current LESLIE3D with the compressible correction. A grid 120x91x32 is 

employed for this simulation. This study is still underway and more results will be 

reported in the near future. Figure 14 shows the results obtained in an earlier study 

(Nelson and Menon, 1998), which compares the typical results for the shear growth rate 

and the normalized velocity profiles. 
NASA/CRó2003-212202 23

 

Figure 12. Iso-surfaces of density and pressure for the supersonic 

shear layer mixing. 

Figure 13. Spanwise vorticity in the supersonic shear layer. 

Figure 14. Comparison of the mixing layer growth and mean velocity profile with 

the experimental results. These results were obtained earlier (Nelson and Menon, 

1998) but are being revisited for LESLIE3D validation in the supersonic regime. 
NASA/CRó2003-212202 24

The next step would be to attempt to compare the numerical simulation results 

with the DNS data for mixing layers with significant compressibility effects. This will be 

an area of research over the next few months. 

 

4.2.4 DNS of Two-Phase Flows in Temporal Mixing Layers 

 Another area of research was in the development of DNS/LES methods for 

simulating single and two-phase mixing process in a temporal mixing layer. This test case 

is also chosen for NCCLES validation for some very important reasons: 

(1) The current NCC code has the capability to do periodic conditions in only 

one-direction. Thus, with this capability, the isotropic turbulence test case 

noted in Section 3 will be quite difficult to accomplish unless major revision 

of the code is carried out. On the other hand, the temporal mixing layer can be 

simulated using periodic in one-direction and slip boundary conditions in the 

other two directions. 

(2) The temporal mixing layer has been used extensively in the past for both DNS 

and LES studies and there is a large amount of data available for sort of flow. 

Also, temporal mixing layer is a subset mixing flow field of direct relevance 

to reacting flows. 

(3) The gas phase temporal mixing layer simulation capability was developed 

many years ago at Georgia Tech. Here, we have extended its capability to 

simulate two-phase flows (e.g., spray mixing). This is with the focus for the 

eventual application to spray combustion, which is the major area of interest 

in the application of NCCLES. 

In the following figures we show some samples of the DNS and LES studies using the 

temporal mixing layers. Results of the present DNS have been compared with DNS 

predictions in the past using a spectral code (Ling et al, 1998). Very good agreement with 

past studies is observed in all cases. The effect of particle size (as determined by the 

Stokes number) was the principal focus these studies. It can be seen that the simulations 

are able to capture all these effects accurately. The effect of using second-order and 

fourth-order spatially accurate schemes using the same grid resolution is also studied 

here. It can be seen that both schemes are reasonably accurate; however, the fourth-order 

scheme is better. In the present case, the second-order scheme is quite good since these 

simulations primarily focused on the rollup and pairing of large-scale coherent structures. 

These coherent structures are primarily 2D in their nature (and as initialized here) and 

therefore, the second-order scheme is reasonably accurate. However, if the mixing layer 

is initialized with 3D isotropic turbulence then it would become apparent that the fourthorder scheme would be superior for the same resolution. 
NASA/CRó2003-212202 25

Figure 15. Evolution of the temporal mixing layer (shown here is the vorticity) 

during vortex pairing in the presence of particles for Stokes number of 4. 

Figure 16. Particle motion and entrainment during vortex pairing (for a Stokes 

Number of 4). Top sequence is from DNS of Ling et al (1998) while the bottom 

sequence is from the present DNS study. 
NASA/CRó2003-212202 26

 

Figure 17. Comparison of past DNS (Ling et al, 1998 top series) prediction of 

particle motion with present DNS (bottom series) results. Figures from left to right 

are for Stokes Numbers of 1, 4 and 100, respectively at the same non-dimensional 

time. 

Figure 18. Modal energy and momentum thickness growth rates in the temporal 

mixing layer. The effect of second and fourth-order accuracy of the DNS is shown. 

Figure 19. RMS of the particle number density in a cell as a function of Stokes 

Number is shown. Comparison with past DNS studies is shown. 
NASA/CRó2003-212202 27

5. Conclusions and Future Plans 

The results reported here show that substantial progress has been made in the 

present research in this project. In addition to providing LESLIE3D codes and helping 

NASA/GRC researchers to understand LES methodology and its application, we have 

been developing additional features that will eventually find application in the NCCLES. 

In particular, three areas of development and validation are underway: (a) FlameTurbulence Interactions using premixed combustion, (b) Spatially evolving supersonic 

mixing layers, and (c) Temporal single and two-phase mixing layers. The configurations 

chosen are such that they can be implemented in NCCLES and used to evaluate the 

ability of the new code. 

 The studies planned in the near future for the first yearís effort is to complete the 

LES-LEM application in the flame-turbulence interactions and to further carry out LES 

of the supersonic spatially evolving shear layers. In addition, the extension of the LESLEM approach for spray combustion will be addressed. In the second year, we anticipate 

further extension of the LES approach to more complex supersonic flows including shock 

waves. 

 In addition, in the second year we will begin LES study of the Trapped-Vortex 

Combustor (TVC). This combustor was chosen as a final test case for NCCLES and 

General Electric Aircraft Engine Company team had agreed to collaborate with the 

NASA/Tech team in this research. Therefore, we will begin discussion with GEAEC to 

identify a generic TVC configuration of their interest and begin setting up the LESLIE3D 

for this configuration. It is anticipated that some new insight into the TVC flow field will 

be achieved during the second yearsí effort under this project. 

References 

Ajmani, K. and K.-H. Chen (2001), ìUnsteady-Flow Validation for the NCC,î AIAA 

2001-0972. 

Bhargava, A., Kendrick, D.W., Colket, M.B., Sowa, W.A., Casleton, K.H., and Maloney, 

D.J. (2000), ìPressure effect on NOx and CO emission from industrial gas turbines,î 

ASME, GT-97. 

Chakravarthy, V.K. and Menon, S. (1997), ìOn Large-Eddy Simulations of NonHomogeneous Flows,î AIAA 97-0652, 35th Aerospace Sciences Meeting, Reno, NV, 

January 6-9, 1997. 

Chakravarthy, V.K. and Menon, S. (2001a), ìModeling of Turbulent Premixed Flames in 

the Flamelet Regime,î Combustion, Science and Technology, Vol. 162, pp. 1-50, 

2001. 

Chakravarthy, V.K. and Menon, S. (2001b) ìSubgrid Modeling of Turbulent Premixed 

Flames in the Flamelet Regime,î Flow, Turbulence and Combustion, Vol. 65. 

Chen, C., Riley, J.J, and McMurthy, P.A. (1991) ìA study of Favre averaging in turbulent 

flows with chemical reaction,î Combustion and Flame, 87:257-277, 1991. 

Calhoon, W. H. and Menon, S. (1997) ìLinear Eddy Subgrid Modeling for Reacting 

Large-Eddy Simulations: Heat Release Effects,î AIAA-97-0368, 35th Aerospace 

Sciences Meeting, Reno, NV, January 6-9, 1997. 
NASA/CRó2003-212202 28

Comte-Bellot, G., and Corrsin, S., Simple Eulerian time correlation of full- and narrowband velocity signals in grid-generated, isotropic turbulence, Journal of Fluid 

Mechanics, 48:273-337, 1971. 

Echekki, T., and Chen, J.H. (1991) ìUnsteady Strain rate and curvature effects in 

turbulent premixed Methane-Air Flames,î Combustion and Flame, 106:184-202, 

1991. 

Goldin, G. and Menon, S. (1997) ìA Scalar PDF Construction Method for Steady State 

Turbulent Combustion,î Combustion Science and Technology, Vol. 125, pp. 47-72, 

1997. 

Kee, R.J., Rupley, F.M., and Miller, J.A., Sandia National Laboratories report SAND8215b, 1987. 

Kim, W.-W. and Menon, S. (1998) ìAn Unsteady Incompressible Navier-Stokes Solvers 

for Large-Eddy Simulations of Turbulent Flows,î International Journal of Numerical 

Methods in Fluids, Vol. 14, 1998. 

Kim, W.-W. and Menon, S. (2000) ìNumerical Modeling of Turbulent Premixed Flames 

in the Thin-Reaction-Zones Regime,î Vol. 160, pp. 17-48, Combustion Science and 

Technology, 2000. 

Kim, W.-W., Menon, S. and Mongia, H. (1999) ìNumerical Simulations of Reacting 

Flows in a Gas Turbine Combustor,î Combustion Science and Technology, Vol. 143, 

pp. 25-62, 1999. 

Lannetti, A., R. Tacina, S.-M. Jeng, and J. Cai (2001) ìTowards Accurate Prediction of 

Turbulent, Three-Dimensional Recirculating Flows with the NCC,î AIAA-2001-

0809. 

Liu, N.-S. (2001) “On the Comprehensive Modeling and Simulation of Combustion 

Systems,” AIAA-2001-0805, 38th Aerospace Sciences Meeting, Jan 2001.

Menon, S. (2000) ìOn Subgrid Modeling for Large-Eddy Simulations of Turbulent 

Combustion,î International Journal of Engine Research, Vol. 1, No. 2., pp. 12-33. 

Menon, S. and Calhoon, W. (1996) ìSubgrid Mixing and Molecular Transport Modeling 

for Large-Eddy Simulations of Turbulent Reacting Flows,î Proceedings of the 

Combustion Institute, 26, 59-66, 1996. 

Menon, S., McMurtry, P., and Kerstein, A. R., ìA Linear Eddy Mixing Model for LES of 

Turbulent Combustion,î in Large-Eddy Simulations of Complex Engineering and

Geophysical Flows, (B. Galperin and S. A. Orszag, Eds.), Cambridge University 

Press, pp. 287-314, 1993. 

Menon, S., Yeung, P.-K. and Kim, W.-W. (1996) ìThe Effect of Subgrid Models on the 

Computed Interscale Energy Transfer in Isotropic Turbulence,î Computer and Fluids, 

Vol. 25, pp. 165-186, 1996. 

Moder, J. (2001) ìValidation of the Radiation Heat Transfer Module Coupled with the 

Flow Solver in the NCC,î AIAA-2001-0810. 

Nelson, C. and Menon, S., ìUnsteady Simulations of Compressible Spatial Mixing 

Layers,î AIAA Paper No. 98-0786, 36th Aerospace Sciences Meeting, Reno, NV, 

January 12-15, 1998. 

Pannala, S. and Menon, S., (1998) ìLarge-Eddy Simulations of Two-Phase Turbulent 

Flows,î AIAA Paper No. 98-0163, 36th Aerospace Sciences Meeting, Reno, NV, 

January 12-15, 1998. 

Poinsot, T., and Lele, S. (1991) ìBoundary Conditions for Direct Simulations of 

Compressible Viscous Flow,î Journal of Computational Physics, 101:1, 1991. 
NASA/CRó2003-212202 29

Raju, M. (2001) ìA Validation Summary of Scalar Monte Carlo PDF Applied to Spray 

Flames,î AIAA-2001-0806. 

Samimy, M., Reeder, M.F., and Elliot, G. (1992) ìCompressibility effects on large 

structures in free shear flows,î Physics of Fluids A, 4(6):1251-1258, 1992. 

Sankaran, V. and Menon, S. (2000) ìThe Structure of Premixed Flames in the thinreaction zone,î Proceedings of the Combustion Institute, Vol. 28, 2000. 

Schumann, U. (1975) ìSubgrid scale model for finite difference simulations of turbulent 

flows in plane channels and annuli,î Journal of Computational Physics, 18:376-404, 

1975. 

Shih, T.-H., Norris, A., Lannetti, A. and Povinelli, L. (2001) ìA Study of a Hydrogen/Air 

Combustor using NCC,î AIAA-2001-0808. 

Smith, T. and Menon, S. (1996) ìModel Simulations of Freely Propagating Turbulent 

Premixed Flames,î Proceedings of the Combustion Institute, 26, 299-306, 1996. 

Sone, K. and Menon, S. (2001) ìLarge-Eddy Simulations of Fuel-air Mixing in Internal 

Combustion Engines,î AIAA 2001-0910, Reno 2001. 

Wei Ling, Chung, J.N., Troutt, T.R. and Crowe, C.T. (1998) ìDirect numerical 

simulation of a three-dimensional temporal mixing layer with particle dispersion,î 

Journal of Fluid Mechanics, 358:61-85, 1998. 

Westbrook, C.K. and Dryer, F.L. (1981) ìSimplified reaction mechanisms for the 

oxidation of hydrocarbon fuels in flames,î Combustion Science and Technology,

27:31-43, 1981. 
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

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 Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

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

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Final Contractor Report

Unclassified

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

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

National Aeronautics and Space Administration

Washington, DC 20546–0001

Available electronically at http://gltrs.grc.nasa.gov

April 2003

NASA CR—2003-212202

CCL–02–003

E–13805

WU–708–87–23–00

NAG3–2660

35

Subgrid Combustion Modeling for the Next Generation National

Combustion Code

Suresh Menon, Vaidyanathan Sankaran, and Christopher Stone

Subgrid mixing; Large eddy simulation; Spray combustion

Unclassified -Unlimited

Subject Categories: 01 and 07 Distribution: Nonstandard

Georgia Institute of Technology

790 Atlantic Drive

Atlanta, Georgia 30332

Project Manager, Nan Suey Liu, Turbomachinery and Propulsion Systems Division, NASA Glenn Research Center,

organization code 5830, 216–433–8722.

In the first year of this research, a subgrid turbulent mixing and combustion methodology developed earlier at Georgia Tech has been

provided to researchers at NASA/GRC for incorporation into the next generation National Combustion Code (called NCCLES

hereafter). A key feature of this approach is that scalar mixing and combustion processes are simulated within the LES grid using a

stochastic 1D model. The subgrid simulation approach recovers locally molecular diffusion and reaction kinetics exactly without

requiring closure and thus, provides an attractive feature to simulate complex, highly turbulent reacting flows of interest. Data

acquisition algorithms and statistical analysis strategies and routines to analyze NCCLES results have also been provided to

NASAJGRC. The overall goal of this research is to systematically develop and implement LES capability into the current NCC. For

this purpose, issues regarding initialization and running LES are also addressed in the collaborative effort. In parallel to this technology

transfer effort (that is continuously on going), research has also been underway at Georgia Tech to enhance the LES capability to tackle

more complex flows. In particular, subgrid scalar mixing and combustion method has been evaluated in three distinctly different flow

field in order to demonstrate its generality: (a) Flame-Turbulence Interactions using premixed combustion, (b) Spatially evolving

supersonic mixing layers, and (c) Temporal single and two-phase mixing layers. The configurations chosen are such that they can be

implemented in NCCLES and used to evaluate the ability of the new code. Future development and validation will be in spray

combustion in gas turbine engine and supersonic scalar mixing.