NASA/TM--2001-210762 AIAA-2001-0289

Development of a Hybrid RANS/LES Method

for Compressible Mixing Layer Simulations

Nicholas J. Georgiadis

Glenn Research Center, Cleveland, Ohio

J. Iwan D. Alexander and Eli Reshotko

Case Western Reserve University, Cleveland, Ohio

Prepared for the

39th Aerospace Sciences Meeting and Exhibit

sponsored by the American Institute of Aeronautics and Astronautics

Reno, Nevada, January 8-11, 2001

National Aeronautics and

Space Administration

Glenn Research Center

April 2001
NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

Available from

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

Available electronically at http://gltrs.grc.nasa.gov/GLTRS
Development of a Hybrid RANS/LES Method

for Compressible Mixing Layer Simulations

Nicholas J. Georgiadis *

NASA Glem_ Research Center, Clevcla_M, Ohio 44135

J. Iwan D. Alexander rand Eli Reshotko *

Case WesterT_ Re:scrt,e _riversity, Cle_,elaTM, Ohio _4106

A hybrid method has been developed for simulations of compressible turbulent mixing

layers. Such mixing layers dominate the flows in exhaust systems of modern day aircraft

and also those of hypersonic vehicles currently under development. The hybrid method

uses a Reynolds-averaged Navier-Stokes (RANS) procedure to calculate wall bounded

regions entering a mixing section, and a Large Eddy Simulation (LES) procedure to

calculate the mixing dominated regions. A numerical technique was developed to enable

the use of the hybrid RANS-LES method on stretched, non-Cartesian grids. The hybrid

RANS-LES method is applied to a benchmark compressible mixing layer experiment.

Preliminary two-dimensional calculations are used to investigate the effects of axial grid

density and boundary conditions. Actual LES calculations, performed ill three spatial

directions, indicated an initial vortex shedding followed by rapid transition to turbulence,

which is in agreement with experimental observations.

Introduction

The use of computational fluid dynamics (CFD) to

assist in the analysis and design of aerospace vehicles

and their components has substantially increased in recent years. For analyzing one particular class of flows,

that of aircraft engine exhaust nozzles, Reynoldsaveraged Navier-Stokes (RANS) codes have been used

extensively by government organizations (i.e. NASA)

and aerospace companies. Exhaust nozzles being developed for modern day commercial aircraft typically'

have nmltiple streams with a core flow and one or more

bypass streams which mix with the high energy core

flow before exiting the nozzle to lower jet noise while

maintaining high thrust levels. Similarly in-NASA's

recent High-Speed Research program, the engine exhaust systems for this proposed supersonic transport

were designed to be mixer-ejector nozzles, which entrain secondary air into the exhaust nozzle to mix with

the core engine stream, again with the goal of sinmltaneously lowering jet noise and maintaining sufficient

thrust. 1 Propulsion systems currently under development for use on hypersonic and reusable space launch

vehicles, such as the Turbine-Based Combined-Cycle

(TBCC) and l_ocket-Based Combined-Cycle (RBCC)

concepts also employ mixer ducts.

The flows in these nozzle systems all have compressible turbulent mixing as the dominant flow character-

*Aerospace Engineer, Nozzle Branch, Senior Member AIAA

lAssoc. Professor, Department of Mechanicaland Aerospace

Engineering, Member AIAA

_Kent H. Smith Professor Emeritus, Department of Mechanical and Aerospace Engineering , Fellow AIAA

istic. RANS codes used by research and development

engineers to analyze these nozzles have employed turbulence models to replace the unsteady turbulent motion with an effective eddy viscosity. Unfortunately, no

turbulence model has been developed to date which is

able to accurately represent the turbulent motion for

such nozzle flows. Validation studies -''3 have shown

that the "state of the art" turbulence models available in production-use RANS codes have major deficiencies in predicting turbulent, mixing in nozzle and

jet flows involving compressibility, high temperatures,

and three-dimensionality.

The known limitations of RANS techniques to calculate complex turbulent flows, coupled with continually increasing computing power, have led to interest in more sophisticated calculation techniques such

as direct numerical simulation (DNS) and large eddy'

simulation (LES). DNS is currently limited by computer hardware to very simple flows at low Reynolds

numbers, and LES, which directly solves for the large

turbulent, scales and limits empirical modeling to the

smallest scales, is becoming practical for more complex flows at, higher Reynolds numbers. Birch 4 and

Bradshaw 5 suggest, that LES techniques offer the best

prospects for improving the capability to calculate turbulent flows, particularly for flow regions not. including

wall boundary layers.

As a result, an LES-based technique is an attractive

option for calculating the mixing dominated regions

of nozzle flows. However, applying such an LES technique simultaneously to the wall bounded regions that

enter the mixing region (which are an important part.

NASA/TM--2001-210762 1
of nmlti-streamnozzlesthatshouldbecalculatedaccurately')will notbepracticalinthenearfuture.This

is becausecomputationalresourcesfar greaterthan

thoseavailabletodaywouldberequiredto capturethe

widerangeofturbulentimeandlengthscalesthatare

inq)ortantin sucha problem.Theseturbulentscales

rangefromverysmalleddies in the wall boundary layers to very large eddies in the developing mixing layer.

Nearly all LES simulations of jet and mixing layer

flows performed to date have placed the inflow of

the computational domain downstream of any wall

bounded regions and have either ignored the upstream

boundary layer effects or used some approximation to

initialize the turbulent mixiug layer. Several authors,

such as Ragab e,r and Hedges s have imposed hyperbolic tangent mean velocity profiles at the plane which

represents the end of the wall boundary layer regions

and the beginning of the mixing region. Ragab then

used the results of a linear stability analysis to generate

perturbations about the mean velocity profile located

at. the mixing plane. Hedges added small amplitude

perturbations to the vertical velocity component in

simulations of heated jet flows.

The difficulty in using an artificially generated inttow, such as that assuming a hyperbolic tangent velocity profile, is that the characteristics of the upstream

boundary layers, including velocity, temperature, and

turbulence profiles, art, not accurately represented.

This is a significant deficiency since the state of the incoming boundary layers have been shown to have substantial effects on the development of turbulent, mixing

layers in the experiments conducted by Bradshaw, _'

Browand and Lat, igo, 1° and Hussain and Ctark. 11

While RANS-based methods have major deficiencies

in predicting compressible mixing layers and inherently are not formulated for calculation of unsteady

turbulent flows, they have been shown to predict the

mean flow behavior of wall bounded regions quite well,

particularly in the absence of adverse pressure gradients. As a result, it would be desirable to combine

a RANS-based technique for the wall boundary layers

upstream of the mixing region with an LES-based technique for the downstream unsteady, turbulent mixing

region. The developmeut of such a hybrid R.ANS/LES

approach is the subject of this work.

The method developed here is proposed as an alternative computational technique to performing LES

calculations everywhere in the computational domain,

that. includes the mean flow characteristics of the incoming boundary layers and is also feasible when considering foreseeable computational resources.

Hybrid RANS-LES Methods

The realization that LES calculations of flows

in aerospace and industrial applications at realistic

Reynolds numbers will not be possible for some time

has led to interest into the development of hybrid

techniques. The objective of a hybrid method is t.o

retain the essential features of the LES method, but,

to employ a computatioually cheaper R ANS method in

regions where it is appropriate. As a result., nearly all

hybrid methods proposed to date apply a llANS approach to attached wall boundary layer regions and an

LES approach to regions of large scale separation. The

work detailed in this paper represents the first hybrid

method development for application to compressible

mixing layers.

The most widely publicized hybrid method to date

is the Detached Eddy Sinmlation (DES) method of

Spalart. 1"-'14 In the DES method, the wall bounded

regions are calculated using ttANS with the SpalartAllmaras 15 one equation turbulence model. (:onstantinescu and Squires 1'_ have applied Spalarl.'s DES

method to turbulent flow over a sphere, which is an

appropriate geolnetry for the method due to the large

scale separation in the wake of the sphere.

Speziale lr suggested an approach thai allows for

computations varying from RANS in the coarse grid

limit, through LES, and finally to DNS in the very fine

grid limit. A Reynolds Stress model is used to close

the turbulent stresses in the RANS limit, and provides

the basis for a subgrid model necessary m LES simulations. Batten et al. is also propose a hybrid model that.

employs a Reynolds-Stress model to (:lose the RANS

and LES equations. Lastly, Arunajatesan et al. 1_' have

applied a hybrid RANS-LES method to cavity flowfields. Their approach employed a two-equation k-kl

turbulence model to close the RANS equations and

a one-equation model solving for the filtered subgrid

kinetic energy to close the LES equations.

The hybrid RANS-LES method presented here is developed for application to configurations such as the

mixing layer shown in figure 1. This relatively simple

configuration is representative of more complex nozzle

geometries in that two wall bounded regions provide

isolated flows to a single region where compressible

mixing is the primary flow characteristic. Development of the hybrid method and preliminary assessment

of the method for a benchmark compressible mixing

layer configuration are the focus of this paper.

The hybrid method employs a RANS approach

to provide the mean flow characteristics of the wall

boundary layers entering the mixing region. The

downstream mixing layer is then calculated using LES.

The method developed here is intended for those nozzle and mixing layer problems in which a geometric

feature, such as the base region of a nozzle or splitter plate separating the upstream flows, provides the

dominant unsteady mechanism to drive the development of turbulence in the mixing layer. Although the

upstream RANS approach does not provide any unsteady turbulent information to the mixing layer, the

mean flow monmntum and thermal boundary layer effects are calculated and provided to the LES region.

NASA/TM--2001-210762 2
Wall Bounded Regions- RANS

/

Free Shear Layer Region - LES

Fig. 1 Schematic of mixing layer demonstrating the hybrid RANS/LES approach

Method Formulation

Both the RANS and LES equations are derived

starting from the general form of the Navier-Stokes

equations, written in tensor form. The expression

for conservation of mass, or the continuity equation

is written in tensor form as:

Op c_

_7 + _ (P_) = o (1)

Conservation of momentum is written:

0 0 OP Orij

0-7(P") + _ (pub,j) = -ox--7+ Ox-7 (2)

Conservation of energy is expressed as follows:

OEt 0 CO Oqj (3)

a--t- + _ (uj (Et + P)) = _ (uirij) - ax----_

Here, the variable Et represents the total energy (internal energy plus kinetic energy) per unit volume:

1

Et _-- pf "JC "_])tli lli (4)

The equation of state for an ideal gas is used to relate

the pressure, temperature, and density through:

P = pRT (5)

For the viscous stresses rij, it is assumed that the fluid

is a Newtonian fluid, and as a result, the viscous stress

is proportional to the rate of strain. This is written:

Ouj ,-ij= 2_,&j + A7:--5q (6)

(J*j

where the rate of strain tensor _%'ij is:

&J = -2\ O.j + Ox_] (7)

Using Stokes's assumption that the thermodynamic

and mechanical pressures are the same for a fluid undergoing and expansion or compression:

Equation (6) can then be rewritten as:

e Ouj_

7"ij = 2ptSij -- 5#_Oij

t1*j

(9)

The Sutherland model is used to calculate the viscosity. The heat flux qj is obtained from Fourier's law:

k OT

qj =- _ (10)

where k is the thermal conductivity. It is assumed

that the fluid is thermally perfect, such that the internal energy and enthalpy are only functions of the

temperature, and it is also assumed that the fluid is

calorically perfect such that the specific heats Cv and

Cp are constants. As a result, the internal energy and

enthalpy can be written as:

e = CvT

(]l)

h = CpT

Assuming that the air is of constant composition and

does not undergo any chemical reaction, the thermal

NASA/TM--2001-210762 3
con{tuctivity is only a function of t.emperature. Using

the specification of constant specific heats, the following expression is obtained ['or tile t.hernlal conductivity

as a flulction of tile const.anl pressure specific heat,

Prandll number, and the viscosity:

k -- IJ("P

Pr (12)

Mass-Weighted RANS Equations

tn the classical form of Reynolds averaging, the time

dependent form of the Navier-Stokes equations given

by (1) through (3) are averaged over a period of time

that is much larger than the period of turbulent fluctuations. Each of the dependent variables appearing

in these equations is replaced by the sum of mean and

fluctuating components. As an example, the velocity

would be given by':

ui = _i + u'i (la)

where the time averaged velocity _7i is given by:

= 1 /t+r uidt (14)

For the current work, where fluctuations in density are

important, a lnass (or density) weighting is employed

ill the averaging process, which makes the final form

of the RANS equations much more convenient to work

with. Tile dependent variables are again broken into

nlean and fltlct.uating conlpotlents:

ui = ui + ul' (15)

where the time averaged (using mass weighting) velocit,y ui is given by:

= 1 /t+r puidt (16)

ill _ ,It

This Inass-weighted Reynolds averaging process is frequently referred t.o as Favre averaging, and in general,

the Favre average of any variable f is defined by

pf

f = -- (lr)

P

Applying this averaging procedure to equations (1) -

(3) results in the following expressions for continuity,

momentum, and energy:

e¥ 0

0--7+ _ (_u,) = o

w

0 0 OP O_j

o-7(-_":")+ _ (>""') + 0.-7 - 0%-}

(18)

- o

c')x.i

(19)

. (.,)+ 0 0

0

+OTj (_j + qf) = 0

(20)

Spatially Filtered LES Equations

To derive the LES equations used in this work, tile

time det)endent form of the Navier-Stokes equations

given in equations (1) through (3) is again used as

tile starting point. Instead of time averaging lhese

equatioas, however, an approach similar to the work

of Ragab and Sheen _'r and Erlebacher et al "-'l_is used

that will filter out small scale fluctuations, and only

retain scales that are large enough to be resolved by'

a particular computational scheme and the coml)utational mesh. The filtering operation is defined on any

variable f by' the expression:

f(x,t) = _ G (x - _, A) f(5,l)da5 (21)

In equation (21), G is the filter function. D is the flow

domain, and A is the filter width. The filter width

A is usually taken to be the grid spacing, and is the

approach taken here. Note that the overbar used in

equation (21) indicates a filtered variable. This is in

contrast to the previous use of an overbar to indicate a

time averaged quantity in the previous section. As discussed by Nelson, 21 the exact form of the filter function

is not typically known. However, the filter function

must satisfy':

D6' (x -- {, _-_)d3{ = ]

(22)

In large eddy simulations of compressible flows, it is

common to use Favre-filtering which is defined as:

f = =Pf {23)

P

where a quantity f is decomposed into resolved and

unresolved (also referred to as sub-grid scale) components as:

f = ] + f' (24)

Equations (23) and (17) are very similar in appearance, but they refer to very different operations. Both

operations employ mass (density) weighting, but equation (17) defines a time averaged quantity, while equation (23) defines a spatially filtered quantity. As was

the case for the RANS equations, the density, pressure,

and heat. flux terms are not decomposed using massweighting. Again note that the overline represented a

time averaging process in the previous section, but it.

will refer to a spatial filtering operation in the current

section. In addition, Favre filtering differs from Favre

time averaging in that:

] ¢:/

and

]' # 0 (26)

NASAfrM--2001-210762 4
Afterapplyingthefilteringproceduretoequations(1)

- (3)theresultingLES expressions for continuity, momentum, and energy are:

0_ 0

0-7+ _ (_i) = 0 (27)

0 0 o-P o_ij ori'j""

0-7(-_a,) + _ (-_i%) + oxi m.j o_j - 0

(2s)

-- 7,-_gs

+_ (_j + qj_ ) ----O

(29)

Turbulence Modeling

Both the RANS and LES sets of equations require

a turbulence model to close the momentum and energy equations. In the RANS approach, all unsteady

turbulent motion is replaced by a turbulence model.

The resulting LES equations are very similar in appearance to the RANS equations, and also require a

model to close the momentum and energy equations.

The difference for the LES equations, however, is that

the terms replaced by a model are only the turbulent terms that are too small to be resolved using the

filtered LES equations. As a result, the large scale turbulent motion is directly calculated, and the effects of

the smallest scale turbulence are accounted for using

a subgrid turbulence model.

The turbulence model employed here to close the

RANS equations is the Cebeci-Smith algebraic turbulence model. 22'23 Since the RANS equations are only

used in this hybrid method to calculate wall boundary layer regions with no adverse pressure gradients,

the selection of a relatively simple algebraic model

such as the Cebeci-Smith formulation is appropriate.

The walt function technique of Ota and Goldberg 24

is used in conjunction with the Cebeci-Smith model

to enable use of a computational grid with the first

point off solid boundaries placed in the logarithmic

layer. This wall function approach is based upon the

compressible law of the wall formulation of White and

Christoph.25'26 The filtered LES equations are closed

using the Smagorinsky subgrid model. 27

hnplementation of the wall function technique is

critical to the development of this hybrid approach in

order to enable use of a single computational grid extending continuously from the RANS regions to the

LES regions. If a wall function approach were not

used, grids for the RANS regions would have to be

packed very tightly to the wall and use significant grid

stretching, while a separate grid which minimizes grid

stretching would need to be constructed for the LES

region. Use of such non-continuous grids for the RANS

and LES regions would require a scheme that would

likely introduce interpolation errors into the combined

hybrid method.

lTse of the Cebeci-Smith model to close the RANS

equations and the Smagorinsky subgrid model to close

the LES equations is desirable in terms of code implementation. While the function of the Cebeci-Smith

model to replace all of the turbulent stresses with a

model is quite different from that of the Smagorinsky

subgrid model, which only replaces the small subgrid

turbulent stresses, both are eddy viscosity models and

are derived at least in part from mixing-length theory.

The similar formulation of these two models enables

the RANS equations and LES equations to be solved

with a single solution scheme and computational grid,

as mentioned previously. For a compressible nozzle or

mixing layer flow, such as that depicted in figure 1.

the change from RANS regions to LES region occurs

at. the vertical plane passing through the trailing edge

of the splitter separating the wall bounded flows.

RANS Turbulence Model

The unclosed terms from the RANS momentum and

energy equations are the Reynolds stress and turbulent

heat flux. The Boussinesq approximation is used to

relate the turbulent Reynolds stress to the mean rate

of strain tensor through a turbulent (or eddy) viscosity.

This is directly analogous to equation (9) which relates

the viscous stress to the mean rate of strain through

the molecular (or laminar) viscosity. The turbulent

analogy to equation (9) is:

_a -P i j

= / 0/_k (30)

Similarly, the turbulent heat flux is related to the temperature gradient through a turbulent conductivity,

kr:

2u'h" q_f=P j

OT (31)

-- __k T --

Oxj

The turbulent Prandt] number, Pr r is used to relate

the turbulent viscosity to the turbulent conductivity:

Pr r _ prCp (32)

k _

The turbulent, Prandtl number is taken to be a constant here and equal to 0.9. Using equation (32) and

assuming a constant, turbulent Prandtl number enables

the turbulent heat flux to be expressed as a function

of the turbulent viscosity that is used to calculate the

Reynolds stress. The turbulent, heat flux becomes:

Cpp ToT

q_--- Pr r Oxj (33)

NASA/TM_2001-210762 5
The(?ebeci-Smithmodel,whichwaschosenhereto

{'losetheRANS equations in the wall boundary layer

regions, treats the wahl boundary layer as having inner and outer regions where the turbulent viscosity is

defined as:

= t,........ .u _<.u,,, (34)

tP_ !,I"L,, Y > Y-,

Ill equation (34), .q,. is defined as the smallest value

of y (the distance away from a wall) at, which p_o .... =

p[_,_ . Tile expressions for the inner and outer layer

turbulent viscosities are as follows:

hmer Layer:

with the mixing length t',.... is given by:

(m,_ = trY(1- e'-v+/A+)

(35)

(36)

Note that the specific form of equation (35) is for a

two-dimensional boundary layer.

Outer Layer:

(37)

In equation (37), the quantity (_,, is the velocity thickness, u_ is the boundary layer edge velocity, and F,_M,

is the Klebanoff intermittency function. Tile velocity

thickness is defined as:

5v = 1 - _ dy (38)

This velocity thickness is identical to the displacement

thickness for incompressible flows. Klebanoff ='s,''9 presented an expression for the intermittency of turbulence near the edge of a boundary, which has a functional form involving the complimentary error function. This original intermittency function is usually

approximated by the following fornmla, as indicated

by Cebeci: 2_

[ (_)6] -1 Fkleb = 1 + 5.5 (39)

The closure coefficients appearing in equations (36)

and (37) are

tc = 0.40 c_ = 0.0168 A + = 26

Wall Function Implementation: The Cebeci-Smith

model is usually integrated down to the wall, using a

computational grid with the first point off of the wall

placed well within the laminar sublayer, corresponding to y+ < 5. For the hybrid method developed in

this work, the objective is to place the first point off

of the wall in the logarithmic layer to enable the use

of computational grids that are not packed as tightly

to the wall. Removing the tight spacing requirement

will enable a continuous grid into the LES region. In

addition, because the allowable time step of the computations is directly proportional to the size of the

smallest grid cell. a less tightly packed grid enables

a larger tinle step for the solution scheme. The wall

function technique of Ola and Goldberg "-'4is one of the

more simple and effective methods currently in use,

and it. is the technique used in this work.

Wall functions have been iml)lemenled most fiequently in conjunction with two-equation k-( models.

The benefits of implementing a wall function for use

with a k-e model are the same as that for the CebeciSmith model used in this work including reducing grid

requirements, and increasing the permissible time step

of the computations.

The use of a wall-function approach is strictly only

valid in flow regions absent of adverse pressure gradients and separations, due to tile assumption that

the law of the wall holds. However, Avva et al. 3°

have shown results for separated flows in which wall

function methods perform no worse than methods integrating to the wall. The intention of the wall funct.ion implementation in this work is to only apply the

method to attached wall boundary layers where the

law of the wall is valid.

The Ota-Goldberg wall function employs the WhiteChristoph 2_'26 compressible law of the wall:

= (4o) t ,< Vo+)

where the compressibility parameter _t is given by:

ru_

_ - 2CpT_, (41)

In equation (40) u+ is the value of u + at the first, poim

off of the wall, V+ is the value of V+ at. the first, point

offof the wall, and 9+ = 0.1287. In equation (41), the

parameter r is the recovery factor, which is typically

1

taken to be Prg for turbulent boundary layers, and

T,_, is the wall temperature. An iteration procedure

is used with equations (40) and (41) to soh, e for u +,

from which the shear velocity u_ can be obtained:

t/q+ U9

. = -- (42)

U r

Finally, the shear velocity is used to compute the wall

shear stress through:

r,,, = pu 7 (43)

The wall shear stress calculated in equation (43) is

then used in the solution scheme for the momentum

and energy equations in the RANS regions.

NASAfrM--2001-210762 6
LES Subgrid Scale Model

The terms that must be closed for LES equations

are the subgrid-scale stress and the subgrid scale heat

flux. The earliest subgrid scale model for LES com-

• ')7 putations was developed by Smagornlsky.- Despite

significant, efforts to develop more sophisticated subgrid scale models, the Smagorinsky formulation is still

widely used, and is itself the foundation upon which

some of the more sophisticated models are derived.

The form of the model is very similar to the CebeciSmith model used for the HANS equations, in that. a

gradient-diffusion mixing-length approach is used.

The Smagorinsky expression for the subgrid scale

stress is:

T_7'= _ (_Tsj - _ _j )

(44)

The parameter 7r is defined:

_r : Sij S'ij (45)

The parameter A is the filter width and as a result, it.

is also used as the length scale that. is characteristic of

the subgrid turbulence. For use with a computational

method, A is usually taken to be the grid spacing. In

three dimensions, with a computational grid having

unequal spacing in the three directions, this subgrid

length scale is usually taken to be:

1

A = (AxAyAz)5 (46)

For computational grids with substantially different,

spacing in the three directions, an alternative form (see

Ragab _) is

i

A: [(Ax)2 -t- (AY)23 -t- (Az)2.] _ (47)

This second expression for the length scale is employed in this work. The constants Cs and CI have

been found to be dependent on the flow under investigation. Rogallo and Moin 31 suggest, a range for Cs

in the range 0.10 < Cs < 0.24. The upper limit on

Cs given by Rogallo and Moin is investigated for the

mixing layer in this work. The constant CI is usually

equal to 0.01, but several authors, including Ragab 7

and Choi et al. 32 mention that, the contribution of the

term involving C/ may not. be important and may be

neglected. This approach is taken in this work, and as

a result, the original expression for the Smagorinsky

subgrid scale stress in equation (44) may be rewritten

as follows:

2 ,_, Off/

tiT" = 2p"@ij -- 5P _xi°iJ (48)

where the subgrid scale turbulent viscosity is given by:

p'_' = -fi(Cs A )'-v/-_ (49)

Note the similar form of equation (49) to the expression for the Cebeci-Smith ilmer region turbulent

viscosity in equation (35). While the mixing length

defined for the Cebeci-Smith model is used to chara.cterize all of the turbulent motion, the length scale

defined here for the Smagorinsky model only characterizes the subgrid-scale motion. Finally, the subgrid

scale heat flux is modeled analogously to that done for

the turbulent heat flux of the BANS equations:

0_ (50)

8xj

where k _ is related to /t "_" through the turbulent

Prandtl number. As in the RANS regions, the turbulent Prandtl number is assumed to be constant in

the LES regions and equal to 0.9. The subgrid scale

heat flux becomes

Cpt_ 'g" ()r

qJ_'--- Pr T c)xj (51)

Solution Procedure

The similar form of the HANS and LES equations

derived here enable both sets to be solved with a single

computational method. In this work, the GottliebTurkel scheme 33 was used. In addition, the equations

are transformed to generalized coordinates and solved

using a set of metric terms that are consistent with the

order of the Gottlieb-Turkel scheme. This was done

to enable the use of the hybrid method on stretched,

non-Cartesian grids• The Gottlieb-Turkel scheme is

cited by several authors, such as Hudson and Long, 34

as providing second order accuracy in time and fourth

order accuracy in space. Bayliss et al. a_ indicated

that the scheme has fourth order accuracy only if At

is of the order (Ax) _, and Nelson 21 showed that the

spatial accuracy of the scheme is strictly only third

order for CFL numbers approaching 1. In the case of

simulations with grids employing significant stretching, however, the spatial accuracy is effectively fourth

order for most. of the computational domain.

Mixing Layer Simulations

Experimental Configuration

The hybrid RANS/LES method has been applied to

one of the benchmark compressible mixing layer experiments of Goebel and Dutton 36-3s in which two

isolated supersonic streams, separated by a splitter

plate, provide the flows to a constant area mixing

section. In particular, case 2 of their experiments is

examined here. A simplified schematic of their experimental configuration provided in figure 2 shows

NASA/TM--2001-210762 7
top flow bottom flow

MachNo. 1.91 1.36

t;(m/s) i"00 39.9

7;(K) 57s _95

T(K) 334 215

a(n_/s) 366 293

P(I,'Pa) 49 4!1

p(kg/m s) 0.51 0.79

_(mm) 2.9 2.5

a'* (ram) 0.90 {}.44

0(ram) 0.2!} 0.21

Table 1 Flow conditions for case 2 of the GoebelDutton experiments

that two isolated streams, in which boundary layers

develop over a splitter plate surface, are brought, together into a constant area mixing section. In all of

their experiments, the higher speed primary stream occurred over the top surface of the splitter plate. The

top stream enters the mixing section axially, while the

bottom stream enters the mixing section at an angle

of 2.5 degrees. The splitter plate thickness has a base

height of 0.5 mm at the trailing edge. Upstream of the

straight sections for the two isolated flows shown in figure 2. contoured nozzle blocks provided the supersonic

flows with nearly uniform exit flow conditions.

The mixing section height was 48 ram, and the overall length of the mixing section available for flowfield

measurentents was 500 ram. The width of the mixing

section was 96 ram, and as a result, the mean flow development could be considered two-dimensional. The

divergence angle of the lower and upper walls of the

mixing section were adjusted in each experiment with

two incoming supersonic flows, to account for boundary layer growth along these two surfaces and to effectively remove any streanlwise pressure gradient.

Single components LDV measurenmnts were used to

calculate the boundary layer, displacement, and momentunl thicknesses of the two streams as they entered

the mixing section. These quantities are provided

along with the other operating conditions of case 2 in

table 1. The documentation of the incoming boundary layer characteristics makes the Goebel-Dutton experiments one of the more thoroughly documented

benchmark data sets available for compressible mixing

layers. In the mixing region, a two-component LDV

system was used to measure the axial and transverse

velocities. In addition, a Schlieren system with a 20 ns

pulse duration was used to obtain nearly instantaneous

snapshots of the mixing layer.

RANS Validation

To begin the investigation of the hybrid method,

initial RANS calculations for the two streams feeding

the mixing layer were constructed. For the Mach 1.91

stream, wall-integration and wall-function grids were

generated. Both grids extended 300 mm axially and

150 nm_ vertically. In addition, both grids used 141

points in the axial direction and 141 points in the vertical direction. The wall-integration grid was packed

to the wall such that the first grid spacing was 0.(}06

mn_, corresponding to an average y+ of 2.5. The wall

function grid had the first point placed at 0.05 nm_,

corresponding to an average y+ of approximately 20.

This wall spacing was chosen for use with the wallfunction grid so that the initial grid spacing at the wall

was exactly 1/10th of the splitter plate thickness in the

experiment. This grid spacing coukt then be continued

into the mixing region, with 10 points spaced equally

in the vertical direction at the base of the splitter.

In the axial direction, both the wall-integration and

wall-function grids were packed to leading and trailing

edges with spacings at the two ends set to 0.10 ram.

corresponding to 1/Sth of the splitter base thickness.

For the Mach 1.36 case. grid construction was very

similar to that used for the Mach 1.91 stream, with the

exception that the axial dimension was shorter for the

Mach 1.36 calculations. This was done because the

Mach 1.3(5 boundary layer thickness in the DuttonGoebel experiment measured at the beginning of the

mixing section was smaller than the Mach 1.91 boundary layer. The grids extended 200 mm in the axial

direction and 150 mm vertically. Both grids used 141

points in the axial and vertical directions. The wallintegration grid was packed such that the first spacing

was 0.006 turn, corresponding to an average y+ of approximately 3.0. The wall function grid placed the first

point at 0.05 nmt, corresponding to an aw_rage y+ of

approximately 25.

For each of the cases discussed in this section, the

wall boundary condition for the calculations was set to

be an adiabatic no-slip surface and the extreme vertical boundary was set as a slip surface. The inflow

boundary was fixed to uniform supersonic flow and

the outflow boundary extrapolated all quantities.

Figure 3(a) compares velocity profiles obtained with

the two approaches to analytical predictions from the

compressible boundary layer method of Tucker a9 for

the Mach 1.36 boundary layer. Figure 3(b) provides

the same comparison for the Mach 1.91 boundary

layer. The solutions provided with and without the

wall-function method provide close agreement with the

Tucker theory. Figure 4 provides a comparison of the

nlomentum thickness growth along the two flat plates

and also indicates close agreement between the calculations and the method of Tucker. These results

verify that the RANS part of the current method

accurately predicts the mean flow behavior of the supersonic boundary layers entering the mixing section.

Two-Dimensional Mixing Layer Calculations

While true LES simulations require computations in

three spatial directions, initial two-dimensional hybrid

calculations of the entire mixing layer configuration

NASAfrM--2001-210762 8
U1 = 700 m/s

Tt._= 578 K

M1= 1.91 Mixingsection

U2 = 399 m/s

Tt.2= 295 K

M2= 1.36

_ SplitterPlate

,,]

48 mm

1

\\

SplitterPlateParameters:

T.E.thickness=0.50mm

Bottomsurfaceangle= 2.5o

not to scale

Fig. 2 Schematic of Goebel-Dutton mixing layer experiment and Case 2 operating conditions

were used to investigate effects of axial grid resolution

and the splitter plate treatment. No subgrid model

was used in these preliminary calculations.

The development of the computational model began by using the results of the Mach 1.91 and Mach

1.36 boundary layer sinmlations. Specifically, the wallfunction solutions were used, because the computational grids employed with the wall-function approach

enabled a continuous grid into the mixing region for

use with the hybrid RANS-LES solver. The objective

was to construct two RANS regions that would provide boundary layer quantities & 6=, and 0 that nearly

matched those measured in the experiment. Because

it would be virtually impossible t.o match all of three

quantities exactly, the momentum thickness 0 was chosen as the key parameter to match the computations

with the experiment. The momentum thickness represents the mean momentum deficit entering the mixing

section and is fundamental to the mixing layer behavior.

Examining the Mach 1.91 boundary layer first, table 1 indicates that the momentmn thickness for this

stream was measured as 0.29 rnrn at the trailing edge

of the splitter plate. For the wall function calculations discussed in the previous section, the momentum

thickness reached 0.29 mm at. a Reynolds number of

3,680,000 (see figure 4(b)), corresponding to an axial

position of 198 mm from the leading edge of the plate.

The next step was to construct a new computational

grid. The axial domain was shortened to 198 mrn while

retaining the same number (141) of axial grid points.

The grid stretching was modified to accommodate the

shorter domain while maintaining the initial and terminal grid spacings. The vertical domain was reduced

to 23.75 rmn to exactly match the height of the Math

1.91 stream in the experiment.. The vertical domain

of the original grid reached 23.75 mm at the 94tll grid

point, so the first 94 points from the original grid were

used in the modified grid. Calculations obtained with

this 141 x 94 point, 198 rnm by 23.75 rnrn grid provided boundary layer quantities identical to that of

the original 141 x 141 grid, further validating the grid

independent, characteristics of the RANS method.

A similar procedure was used for the Mach 1.36

boundary layer. Table 1 indicates that the momentum

thickness for this stream was 0.21 mm at the splitter trailing edge. Examination of the wall-function

solution obtained for the Math 1.36 boundary layer revealed that the momentum thickness became 0.21 ram

at an axial position of approximately 120 rnrn, corresponding to a plate Reynolds number of 2,720,000. A

modified grid was constructed using 141 axial and 94

vertical points for this Math 1.36 case, corresponding

to a physical domain of 120 mm by 23.75 ram. Calculations with this grid provided a solution identical to

that obtained with the original 141 x 141 grid.

The last step before constructing the entire RANSLES computational grid was to extract the grid points

from the modified wall-function grids just discussed,

extending from 81-141 in the axial domain. The solutions were then held constant at the 81st axial station

as the inflow of the final hybrid RANS-LES grid. This

was done to reduce the grid requirements for the subsequent RANS-LES computations. Boundary layer calculations subsequently obtained with these two shortened axial grids using 61 axial points and 94 vertical

NASA/TM--2001-210762 9
y/h

y/a

Fig. 3

lations

1,2 . . .

1

0.8

0.6

0.4

0.2

0

0

! - . . ! - . . ! ......

-- Tucker Correlation (1/6.37 profile) ]

..... C-S model, no wall function

......... C-S model, O-G wall function 1

, , , i

0.2 0.4 0.6 0.8

U/U

a) Mach 1.36

1.2

• " " I " " • t " " • 1 • • • i • • •

Tucker Correlation (1/6.275 profile) I

..... C-S model, no watlfunction I

1 ......... C-S model, O-G wall function

o.6 _l i

0.4 i i .. .

0 • • • ,

0 0.2 0.4 0.6 0.8

U/U=

b) Mach 1.91

Velocity profiles for boundary layer caleupoints again returned the same boundary layer quantities and velocity profiles as the other solutions. These

were the RANS grids that were used to join with the

LES mixing region computational domain.

Three axial grid spacings (with 200, 400, and 800

points in the mixing section respectively) were examined in the initial two-dimensional hybrid calculations.

In all cases, however, the vertical spacing from the

two RANS regions was continued throughout the entire LES region. With the tightest vertical spacing of

tile wall-function boundary layer solutions set. to 0.05

ram. 10 grid spacings were used vertically across the

splitter base. As a result, all of the hybrid grids used

197 vertical points in the mixing region.

Fixed inflow boundary conditions were used at the

RANS inflows. The fixed inflow for the Mach 1.91

upper stream was placed at an axial position 67 mm

upstream of the splitter plate trailing edge and the

Mach 1.36 lower stream inflow was placed 42 rnm up-

(mm)

0.35

0.3

0.25

0.2

0.15

0.1

0.05

0

.... i .... i ........

-- Tucker Correlation ] I

..... C-S model, no wall function I

......... C-S r_xle,, O-G wall funct on t ................. t ....

..,s>" i -

1 10 _ 2 10 _ 3 10 s 4 106

Re x

a) Mach 1.36

I:)(ram)

0.5

0.4

0.3

-- Tucker Correlation I

..... C-S modal, no walt function I

. ........ C-S model, O-G wall funcl_on [

B

o.2j

0.1

0 .... I .... i .

0 1 106 2 10 e 3 10 e

F ....

J

I "--£--

i

1,". _ 1

!

4 10 _

t . .

10 s

Re x

b) Mach 1.91

Fig. 4 Momentum thickness growth along supersonic fiat plates

stream of splitter tip. At the outflow of the mixing

section, corresponding to an axial distance of 300 mm

from the trailing edge of the splitter plate, an extrapolation boundary condition was used, which is appropriate for the mixing supersonic flow exiting the axial

domain. The top and bottom walls of the mixing section were approximated as slip walls, and no attempts

were made to simulate boundary layers developing on

these two surfaces. The Goebe]-Dutton experiment

configuration was specifically designed with adjustable

divergence angles for these walls to account for boundary layer growth, and to thereby provide data that. may

be directly compared to calculations which do not include the mixing section boundary layers.

The first two-dimensional simulation investigating

grid density effects was for the coarsest grid using 200

axial points. A vortex shedding pattern was observed

to originate from the trailing edge of the splitter plate

due to the separation of the two flows leaving the wall

bounded RANS regions and entering the LES mixing

NASA/TM--2001-210762 10
section. Tile vortex shedding quickly dissipates and

the flow appeared to be mostly laminar throughout

the rest, of the mixing section.

Tile second two-dimensional sinmlation was for tile

computational grid using 400 axial points. Tile characterist.ics of this flow were substantially different from

those of the 200 grid point, case. Instantaneous density

contours are provided in figure 5. A stronger vortex

shedding is evident for this case. although the vortex

strength gradually' dissipates back to a laminar state

approximately 40 mm downstream of the splitter trailing edge. An instability fornls at an axial position of

80 mm which in turn slightly dissipates before resuming growth at a position 180 mm downstream of the

splitter plate trailing edge. A turbulent, pattern then

grows fronl this location to the exit at. x =300 ram.

The 400 point axial grid consisted of smaller axial

grid steps than the 200 axial point, grid and had a

significantly reduced stretching factor relative to the

coarse grid. As a result, the truncation error term in

the Gottlieb-Turkel scheme that effectively smooths

discontinuities is substantially reduced for the 400 axial point, grid, and the capability to resolve unsteady

flow behavior was improved.

The third computational grid investigated had 800

axial points. A substantially different, flow behavior

was also observed for this case compared to the solutions obtained with 200 and 400 axial points. The

density, contours in figure 6 again indicate a vortex

shedding pattern that. originates from the trailing edge

of tile splitter plate, but. unlike the 400 point case, the

solution does not. return to a laminar state before transitioning over to a turbulent-like pattern at m = 80

ram. The very tight axial spacing for this case is sufficient, to minimize the truncation error damping effects

on the unsteady flow developnlent.

An additional two-dimensional grid with 800 axial

points was constructed for a modified splitter geometry

in which the splitter trailing edge is reduced to a sharp

tip, and as a result, the flow separation and vortex

shedding is removed. Figure 7 provide instantaneous

density contours for the case with a sharp trailing edge.

As expected, the vortex shedding evident, in the baseline case was removed in the current case with the

sharp tip. The lack of a separation region results in an

initially laminar mixing layer through the beginning

of the mixing section. Interestingly, the laminar flow

begins to transition to a more turbulent structure at,

nearly the same position observed for the baseline case,

at. approximately x = 80 ram. The structures remain

relatively small until x = 150 mm where large scale

turbulence forms. These structures are more similar

to the Brown-Roshko organized structures than were

those of the baseline case with 800 axial points.

Three-Dimensional Mixing Layer Calculations

In the previous section, two-dimensional calculations were used to construct, the initial computational

model of the mixing layer and to examine preliminary

effects of grid resolution and splitter plate treatment.

To correctly investigate the capability of the hybrid

method, however, LES calculations obtained in three

spatial directions with the use of a subgrid scale model

were performed.

The grid topologies and boundary conditions used

for the three-dimensional simulations were very similar to those used for the t.wo-dimensional simulations.

The two-dinlensional computational grids with 200,

400, and 800 axial point grids were used to construct tile three-dinlensional grids used here. To add

the third computational direction in each case, the

two-dimensional planar grid was copied to provide 11

points in the third (or z) direction. The grid spacing in the z direction was uniform and set. equal to

the axial spacing at the splitter trailing edge, or ,Xz =

L__,X1 _--- 0.10 _'nrtl. Because of the very small number of

grid points used in the z direction and the small physical space that is represented, only very' snlall wave

components in this direction could be simulated, and

a periodic boundary condition was used in this direction.

All of the other boundary conditions and the solution procedure are identical to that used for the

two-dimensional nlixing layer calculations, with the

one exception being that the Smagorinsky subgrid

scale model was used in these three-dimensional simulations. The switch from the RANS regions to the

LES region at. the mixing plane (corresponding to a

vertical plane drawn through the trailing edge of the

splitter plate) was accomplished by' changing the eddy

viscosity used in the flow solver from the Cebeci-Snlith

turbulence model to the Smagorinsky subgrid scale

model. As a result, the effect of the eddy viscosity

changes from that of replacing all of the turbulent

stresses in the RANS regions to that of only replacing the subgrid stresses in the LES regions.

The first, three-dimensional simulation was obtained

using the computational grid with 200 axial points.

The behavior of the flow just downstream of the splitter trailing edge is very similar to the corresponding

two-dimensional calculation. The lack of adequate axial grid resolution results in a rapid dissipation of the

initial vortex pattern, and is even more rapid in the

three dimensional case due to the dissipative nature of

the Smagorinsky subgrid scale model.

The three-dimensional grid using 400 axial points

was utilized next.. The density contours shown in figure 8 indicate a fundamentally different flow structure

than observed with any of the two-dimensional calculations. In particular, the vortex shedding is observed

to transition into a turbulent, pattern much closer to

the trailing edge. While the Schlieren photographs

NASA/TM--2001-210762 11
a) Beginning of mixing section

F ........... .,__(_;%,,_, _LYZ-_ ....... : _ .....

/ ; ,?- " ::_<L%_.*_-, ._ _-0.;_,. ' I " .- .. ', ", f? ; i , ' , _-" ,'1

I I I ...... _ ffl-I "1 Ih I/ I I I " I " I i I [ / I ' I (-i i I

b) Entire mixing section

Fig. 5 Instantaneous density contours for the 400 axial grid point case (2D)

! ""._, ">A_\\_ \

i " " \ -" '?_,X

// .''i s' _',_q

• i( i" ,....

a) Beginning of mixing section

Fig. 6 Instantaneous density contours for the baseline 800 axial grid point case (2D)

NASA/TM--2001-210762 12
a) Beginning of mixing section

.::_ ........ .:::.:L ..... ", _........ ;_ .,, ......... :> .-- +'- -< " ' -3"- • _'<,_ " " _"

b) Entire nlixing section

Fig. 7 Instantaneous density contours for the 800 axial grid point case with a sharp trailing edge for the

splitter (2D)

- . , , ,.

. - • , • , ; ,A  ̧ ": _ • • . •

-. _ '/. . .. . -.. ,. _, . . ... _ _. •

..._ ,. ,,. ...... , ......... h:_<_"..,--; ;-',_:":"_'"-_-..... <

a) Beginning of mixing section

b) Entire mixing section

Fig. 8 Instantaneous density contours for the 3D case using 400 axial points

NASA/TM--2001-210762 13
ylH

0.4

o Goet_H-Dutton Data, 150 mm i

0.3 _D Goebel-Du_on Data, 200 mm i

0.2 :ill' _.......................................................................................................

: i I .. _ O i •

i b , _ i

_o i ! ! I !

• _ i _

-0.3 _

-OA

..0.2 0 0.2 0.4 0.6 0.8 1 1,2

(U -U2) / AU

Fig. 9 Velocity profiles in the mixing section

from the Goebel-Dutton experiments did not probe the

flow details near the splitter plate trailing edge, such

details were examined by Clemens and Mungal 4° for

a similar experiment configuration and flow operating

conditions. A Schlieren photograph from the ClemensMungal experiments indicates an initial flow structure

similar to that observed in the three-dimensional calculation with initial vortex shedding from the trailing

edge of a splitter plate followed by a rapid transition

into turbulence. The turbulent structure in the experiment was observed to be of primarily small scales,

while the LES calculations, by definition, only caplure the large scale structures. While the length of

the organized vortex structure in the calculations does

not exactly match that of the Clemens-Mungal experiment, the three-dimensional calculations predict

the transition location nmch more accurately than the

previous two-dimensional results. In addition, these

results verify" that LES calculations must be run in

three dimensions to allow instabilities to form in all

three directions.

Mean axial profiles obtained from this threedimensional calculation are compared to experimental

data of Goebel and Dutton in figure 9. The calculation

indicates greater shear layer spreading than shown by

the experimental data. The profiles of u_,n_ in figure

10(a) and v_,_, in figure 10(b) generally indicate overpredictions from the calculation, which corresponds

to the wider axial velocity profiles in figure 9. Liou

et al. 41 and Inoue 42 also reported large overpredictions in the turbulence intensities for planar mixing

layers. With the confinement of the current threedimensional calculations to a very small domain in

the z direction, the inability to calculate large scale

fluctuations in this direction may be responsible for

l.he overpredictions of u,,_, and ¢_,,,. In addition, the

Goebel-Dutton Schlieren photographs indicated a very

fine turbulent structure contained within the larger

scale development., while the LES calculations inhery/b

y/b

Fig. 10

tion

2

1.5

1

0.5

0

-0.5

-1

-1.5

-2

2

1.5

1

0.5

0

-0.5

-1

-1.5

-2

= i1! I "

O_ _ i / ° Goebel-Dutton Data, 200 mm

---¢_-_z---F 4 -- Hybrid CalculatiOn. 150 mm

q' i i i

ooi = ! I --

i° o l -i

o° ,, i :

.._ i ._, i ....

0,05 0.1 0.15 0.2

Urm s / AU

a) _lrms

i oao= T"__ 1

i °o ../

! o o • ..7 __i

=i" 7

_ °,, i ____ i

z _ i- - i

. • . . i .... i . . . . i ....

0.05 0.1 0.15 0.2

VrmJ A U

b) _',,_ _

Turbulence intensities in the mixing seeently can resolve only the larger turbulent scales.

The last, three-dimensional case was run using 800

axial points in the mixing section. The finer grid used

reduced the permissible time step by nearly a factor of

two relative to the 400 axial point case. This reduced

time step and doubling the number of grid points in the

mixing section would require a factor of four increase

in the computer CPU time requirements to run this

case to completion, relative to the case with 400 axial

points. Considering that, the 400 axial point case required 500 CPU hours on a Cray C90 computer, 2000

Cray C90 hours would be required to complete the

case with 800 axial grid points. As a result, this last,

three-dimensional case was run long enough to allow

the flow to fully develop in the mixing section, but not

long enough to enable time averaging of the turbulent

statistics.

In contrast to the grid refinement studies performed

for the two-dimensional cases, the large scale turbulent

NASA/TM--2001-210762 14
development did not, change significantly when increasing the number of axial points in the three-dimensional

computations from 400 to 800. The instantaneous

density contours in figures 11 indicate a large scale

turbulent structure which closely resembles those of

tile 400 axial grid point ease. In particular, the breakdown of the organized vortex structure to turbulence

is very similar to those previously shown for the 400

axial grid point case. Within the large scale structures, more fine scale turbulence is evident in the 800

axial point ease. This behavior corresponds directly

with the philosophy of LES in that as the computational grid is refined, smaller structures are able to

be resolved and tile role of the subgrid scale model is

reduced. In the idealized limit of a grid which is sufficiently fine to resolve all turbulent length scales, a

direct numerical sinmlation is obtained.

Conclusions

The work described here represents the initial efforts

to develop and evaluate a hybrid RANS-LES method

for compressible mixing layer simulations. Although

the RANS approach does not provide any unsteady

turbulent information to the LES region, the mean

flow boundary layer characteristics are provided. The

hybrid method was developed for the analysis of nozzle

and mixing layer configurations in which a dominant

structural feature, such as the base region of a nozzle

or splitter plate separating the upstream flows, will

provide the dominant unsteady mechanism to drive

the development of turbulence in the mixing layer.

The Cebeci-Smith turbulence model, despite its relatively simple form, was demonstrated to provide accurate calculations of boundary layer flows in the

RANS regions that are free of adverse pressure gradients or separations. Further, the use of the CebeciSmith model in conjunction with the Ota-Goldberg

wall function enabled calculations of supersonic wall

boundary layers to nearly the same accuracy as that of

the standard approach of integrating the Cebeci-Smith

model through the viscous sublayer, while enabling a

significantly larger vertical grid spacing near tile wall.

The wall function approach enabled a continuous computational grid to be used from the RANS to the LES

regions, and the method thereby avoided the use of

discontinuous grid zones that would have required an

interpolation scheme at the RANS-LES interface. In

addition, the origins of the Cebeci-Smith RANS turbulence model and the Smagorinsky LES subgrid scale

model are both in mixing length theory, and this similar form of the models assisted in code implementation.

As a result, the use of a more sophisticated turbulence

model to close the RANS equations was found to be

unnecessary, provided the RANS regions are restricted

to attached, zero pressure gradient wall boundary layer

regions.

While true LES calculations require computations in

three spatial directions, two-dimensional sinmlations

of a benchmark nfixing layer experiment were considered first to address effects of axial grid resolution and

splitter plate treatment. The parametric study of axial grid resolution indicated more realistic turbulent

development with increasing axial grid density. For all

of the cases examined, a vortex shedding was found

t.o originate from the base region of a splitter plate

separating the upst.ream wall bounded regions. For

the finest grid examined, the unsteady vortex pattern

eventually trausitioned to a turbulent structure. The

location of this transition, however, was much further downst.ream than observed in the experiments.

Calculations obtained for the case in which the finite

thickness splitter base was changed to a sharp tip indicated that the vortex shedding was removed.

Three-dimensional calculations were obtained next

for grids constructed by copying the two-dimensional

planar grids to locations in the third computational

direction. Only a small domain was modeled in

this third direction, and periodic boundary conditions

were employed along the extreme boundaries. For

the coarse three-dimensional grid, again no turbulent

flow development was observed. For the intermediate grid, the vortex shedding found previously in the

two-dimensional simulations was also observed in the

three-dimensional calculations. However, the organized vortical structure rapidly disintegrated into a

significantly more realistic turbulent flow structure.

This rapid transition to turbulent, flow was nearly identical to that found in experimental investigations of a

similar mixing layer configuration. Although the extent of the third dimension in these calculations was

very small, an unsteady mechanism by which disturbances could develop in all three directions and result

in a rapid transition to turbulent flow was enabled by

the three dimensional calculations. In contrast,, a twodimensional approach, by definition, does not. allow for

such three-dimensional disturbances to develop. The

results of these calculations verified that. LES simulations must be performed in three dimensions. A final

three-dimensional calculation was investigated using a

computational grid constructed from the most densely

packed two-dimensional grid. Because a prohibitively

long run time would be required to complete this

solution for turbulent statistics purposes, the calculations were run only long enough to allow the flow

to fully develop. The initial transition to turbulence

and large scale turbulent structures evident, for this

case were very similar to those for the intermediate

three-dimensional case. More resolution of the finer

turbulent scales contained with the larger structures

was observed for the fine grid case, which is in line

with the philosophy of LES to resolve finer scales as

the grid density is increased.

Comparisons of time-averaged axial velocities and

turbulence intensities from the calculations to experiNASA/TM--2001-210762 15
a) Beginning of mixing section

:,,,,!3.!='_. ;:_: _._>?._,9_'_*•!-:'-:_r_;i" "_.... " !"-

I............................................. .... . _..../ " ........ ......

b) Entire mixing section

Fig. 11 Instantaneous density contours for the 3D case using 800 axial points

mental data indicated reasonable agreement, with the

solutions indicating somewhat higher levels of turbulent mixing. A major source of discrepancy between

the calculations and experiment is believed to be the

lack of adequate grid resolution to resolve the small

turbulent scales contained within the larger turbulent

structures. Another source of discrepancy was the very

small domain used in the third computational direction. Despite these limitations, the three-dimensional

calculations demonstrated the success of the hybrid

method to capture the dominant characteristics of the

mixing layer, and in particular, the rapid transition of

the organized vortex structure to a turbulent mixing

layer structure. It is expected that improvements in

the fidelity of the solution scheme, and more importantly, improvements in computing power, will enable

better predictions of the turbulent statistics.

References

t J.M. Seiner and E.A. Krejsa. Supersonic Jet Noise and the

High Speed Civil Transport. AIAA Paper 89-2358, 1989.

2T.J. Barber, L.M. Chiappetta, J.R. DeBonis, N.J. Georgiadis, and D.A. Yoder. An Assessment of Parameters Influencing the Prediction of Shear Layer Mixing. A1AA Journal of

Propulsion and Power, 15(1):45-53, 1999.

3N.J. Georgiadis, D.A. Yoder, and J.R. DeBonis. A Comparison of Three Navier-Stokes Solvers for Exhaust Nozzle Flowfields. AIAA Paper 99-0748, 1999.

4S.F. Birch. The Role of Structure in Turbulent Mixing.

AIAA Paper 97-2636. 1997.

sP. Bradshaw. Turbulence: The Chief Outstanding Difficulty of Our Subject• In The Fifth Symposium on Numerical

and Physical Aspects o,f Aerodynamic Flows. California State

University, 1993. Stewartson Memorial Lecture.

SS.A. Ragab and S.-C. Sheen. An Investigation of FiniteDifference Methods for Large-Eddy Simulation of a Mixing

Layer. AIAA Paper 92-0554, 1992.

7S.A. Ragab and S.-C. Sheen. Large Eddy Simulation

of Mixing Layers. In B. Galerpin and S.A. Orszag, editors,

Large Eddy Simulation o.f Complex Engineering and Geophysical Flows, pages 255 285. Cambridge University Press, 1993.

sL.S. Hedges and S. Eberhardt. Numerical Simulation of

Total Temperature Separation in Jets. AIAA Paper 92-0535,

1992.

9p. Bradshaw. The Effect of Initial Conditions on the Development of a Free Shear Layer. Journal of Fluid l"vlechanics,

26:225-236, 1966.

I°F.K. Browand and B.O. Latigo. Growth of the TwoDimensional Mixing Layer from a Turbulent and Nonturbulent

Boundary Layer. Physics of Fluids, 22{6):1011-1019, 1979.

I_A.K.M.F. Hussain and A. R. Clark. lTpstream Influence

on the Near Field of a Plane Turbulent Jet. Physics of Fluids,

20(9):1416-1425, 1977.

_2P.R. Spalart, W.H. Jou, M. Strelets, and S.R. Allmaras.

Comments on the Feasibility of LES for Wings, and on a Hybrid

RANS/LES Approach. In First AFOSR International Conference or_ DNS/LES, 1997.

_P.R. Spalart. Trends in Turbulence Treatments. AIAA

Paper 2000-2306, 2000.

14p.R. Spalart. Strategies for Turbulence Modeling and

Simulations. International Journal oJ Heat and Fluid Flow,

21(31:252 263,2000.

15p.R. Spalart and S.R. Allmaras. A One-Equation Turbulence Model for Aerodynamic Flows. La Recherche Aerospatiale,

1(1):5-21, 1994.

16G.S. Constantinescu and K.D. Squires. LES and DES Investigations of Turbulent Flow over a Sphere. AIAA Paper

2000-0540, 2000.

ITC.G. Speziale. Turbulence Modeling for Time-Dependent

RANS and VLES: A Review. AIAA Journal, 36(2):173-184,

1998.

_sp. Batten, U. Goldberg, and S. Chaknavarthy. Sub-Grid

NASA/TM----2001-210762 16
Turbulence Modeling for Unsteady Flow with Acoustic ttesonance. AIAA Paper 2000-0473. 2000.

i9S. Arunajatesan, N. Sinha, and S. Menon. Toward Hybrid

LES-RANS Computations of Cavity Flowfields. AIAA Paper

2000-0401, 2000.

2°G. Erlebacher, M.Y. Hussaini, C.G. Speziale, and T. A.

Zang. Toward the large-Eddy Sinmlation of Compressible Turbulent Flows. NASA CR 187460, 1990. see also ICASE Report

90- 76.

21C.C. Nelson. Simulations of Spatially Evolving Compressible Turbulence Using a Local Dynamic 5"ubgrid Modal. PhD

thesis. Georgia Institute of Technology, 1997.

22T. Cebeci. Calculation of Compressible Turbulent Boundary Layers with Heat and Mass Transfer. AIAA Paper 70-741,

1970.

23T. Cebeci and A.M.O. Smith. Analysis of Turbulerlt

Boundary Layers. Academic Publishers, 1974.

24D.K. Ota and U.C. Goldberg. Law of the Wall with Consistent No-Slip Limit. AIAA Paper 97-2246, 1997.

25F.M. White and G.H. Christoph. A Simple Analysis of

Compressible Turbulent Two-Dimensional Skin Friction Under

Arbitrary Conditions. AFFDL Technical Report 70-133, 1970.

26F.M. White. Viscous Fluid Flow. McGraw-Hill, 1974.

273. Smagorinsky. General Circulation Experiments with

the Primitive Equations. ,_onthly Weather Review, 91:99 164.

1963.

2sp.s. Klebanoff. Characteristics of Turbulence in a Boundary Layer with Zero Pressure Gradient. NACA Technical Note

3178, 1954.

29p.s. Klebanoff. Characteristics of Turbulence in a Boundary Layer with Zero Pressure Gradient. NACA Report 1247,

1955.

3°lR. Avva, C. Smith, and A. Singhal. Comparative Study of

High and Low Reynolds Number Version of k-e Models. AIAA

Paper 90-0246, 1990.

3Ltt.S. 1Rogallo and P. Moin. Numerical Simulation of Turbulent Flows. Annual Review of Fluid Mechanics, 17:99 137,

1984.

32D. Choi, T.J. Barber, and L.M. Chiappetta. Large Eddy

Simulation of High-Reynolds Number Jet Flows. AIAA Paper

99-0230, 1999.

33D. Gottlieb and E. Turkel. Dissipative Two-Four Methods

for Time-Dependent Problems. Mathematics of Computation,

30( 136):70.3-723. 1976.

a4D.A. Hudson and L.N. Long. Time Accurate Application

of the MacCormack 2-4 Scheme on Massively Parallel Computers. In ICASE/LaRC Workshop on Be rachrnark Problems in

Computational Aeroacoustics, pages 209-215, 1995.

35A. Bayliss, P. Parikh, L. Maestrello, and E. Turkel. A

Fourth-Order Scheme for the Unsteady Compressible NavierStokes Equations. NASA CB 177994, 1985.

36S. G. Goebel. An Experimental Investigation of Compressible Turbulent MiJ'in9 Layers. PhD thesis, University of Illinois

at Urbana-Champaign. 1990.

37S. G. Goebeland J.C. Duttou. Velocity Measurements of

Compressible Turbulent Mixing Layers. AIAA Paper 90-0709,

1990.

38S. G. Goebel and J.C. Dutton. All Experimental Stud.,,'

of Turbulent Compressible Mixing Layers. AIAA Journal,

29(4):538-546, 1991.

39M. Tucker. Approximate Calculation of Turbulent

Boundary-Layer Development in Compressible Flow. NASA

Technical Note 2337. 1951.

4°N.T. Clemens and M.G. Mungal. Two- and ThreeDimensional Effects in the Supersonic Mixing Layer. AIAA

Journal, 30(4):973 981, 1992.

4lT.M. Liou. W.Y. Lien, and P.W. Hwang. Compressibility

Effects and Mixing Enhancement in Turbulent Free Shear Flows.

AIAA Journal. 33(12):2332-2338, 1995.

420. Inoue. Vortex Simulation of a Turbulent Mixing Layer

AIAA Journal, 23(3):367-373, 1985.

NASA/TM--2001-210762 17
REPORT DOCUMENTATION PAGE FormApproved

OMB No. 0704-0188

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.

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

April 2001 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

Development of a Hybrid RANS/LES Method for Compressible

Mixing Layer Simulations

6. AUTHOR(S) WU-708-90--43-00

Nicholas J. Georgiadis, J. Iwan D. Alexander, and Eli Reshotko

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

REPORT NUMBER

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field E-12725

Cleveland, Ohio 44135 -3191

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

AGENCY REPORT NUMBER

National Aeronautics and Space Administration

Washington, DC 20546-0001 NASA TM--2001-210762

AIAA-2001-0289

11. SUPPLEMENTARY NOTES

Prepared for the 39th Aerospace Sciences Meeting and Exhibit sponsored by the American Institute of Aeronautics and

Astronautics, Reno, Nevada, January 8-11, 2001. Nicholas J. Georgiadis, NASA Glenn Research Center; J. Iwan D.

Alexander and Eli Reshotko, Case Western Reserve University, Cleveland, Ohio 44106. Responsible person, Nicholas J.

Georgiadis, organization code 5860, 216--433-3958.

12a. DI._T_tBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

Unclassified - Unlimited

Subject Categories: 02 and 34 Distribution: Nonstandard

Available electronically at htto://_ltrs._c.nasa.gov/GLTR$

This publication is available from the NASA Center for AeroSpace Information, 301-621-0390.

13. ABSTRACT (Maximum 200 words)

A hybrid method has been developed for simulations of compressible turbulent mixing layers. Such mixing layers

dominate the flows in exhaust systems of modern day aircraft and also those of hypersonic vehicles currently under

development. The hybrid method uses a Reynolds-averaged Navier-Stokes (RANS) procedure to calculate wall bounded

regions entering a mixing section, and a Large Eddy Simulation (LES) procedure to calculate the mixing dominated

regions. A numerical technique was developed to enable the use of the hybrid RANS-LES method on stretched, nonCartesian grids. The hybrid RANS-LES method is applied to a benchmark compressible mixing layer experiment.

Preliminary two-dimensional calculations are used to investigate the effects of axial grid density and boundary conditions.

Actual LES calculations, performed in three spatial directions, indicated an initial vortex shedding followed by rapid

transition to turbulence, which is in agreement with experimental observations.

14. SUBJECT TERMS

Turbulent mixing; Nozzle; Compressibility; Large eddy simulation; Navier-Stokes

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

19. SECURITY CLASSIRCATION

OF ABSTRACT

Unclassified

NSN 7540-01-280-5500

15. NUMBER OF PAGES

23

16. PRICE CODE

20. LIMITATION OF ABSTRACT

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102