NASA/TM--2001-210716

The Numerical Analysis of a

Turbulent Compressible Jet

James R. DeBonis

Glenn Research Center, Cleveland, Ohio

National Aeronautics and

Space Administration

Glenn Research Center

May 2001
Acknowledgments

This dissertation would not have been possible without the efforts and support of my teachers, colleagues, and friends.

I would like to thank my advisor Dr. Jim Scott for both his guidance and encouragement during my studies and also

for making this process an enjoyable one. I would also like to thank Dr. Richard Bodonyi, Dr. Michael Foster, and

Dr. Mohammad Samimy for reviewing this document. They have made my studies at Ohio State very rewarding.

This work was done under the doctoral study program at the NASA Glenn Research Center. I am very grateful

for the support of my supervisors, Mr. Bernard Blaha and Dr. Rickey Shyne, and project managers,

Ms. Mary Jo Long-Davis and Mr. Waldo Acosta. During my research, I have sought the help and

advice of many colleagues. In particular I would like to thank Mr. Ibraheem AI-Qadi at Ohio State,

Mr. Dennis Yoder at NASA Glenn, Dr. Ray Hixon at the Ohio Aerospace Institute, and

Dr. Tom Barber at UTRC. I am grateful to Dr. Jay Panda for his high quality experimental

data and for the many informative discussions on jet flows. A special note of thanks also

goes to my NASA office-mate and fellow Ph.D. candidate Nick Georgiadis. Nick and I

have spent many, many hours "discussing" the finer points of Large-Eddy Simulation.

My research is the better for it.

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.eTc.nasa.gov/GLTRS
TABLE OF CONTENTS

Page

List of Tables ....................................

List of Figures ...................................

List of Syml)ols ...................................

('hapters:

1. Introduction ..................................

1.1 Motivation ...............................

1.2 Physics of Jets .............................

1.:3 Computational Methods ........................

1.3.1 Reynolds Averaged Navier-Stokes Simulations ........

1.3.:? Direct Numerical Simulations .................

1.a.a Large-Eddy Simulations ....................

1.4 Discussion of the Present Work ....................

'2. Governing Equations .............................

"2.1 Tile Navier-Slokes Equations .....................

'2.'2 Axisymmetric Formulation .......................

2.a Flux \/ector Forms ...........................

').3.1 Three-dimensional .......................

2.3.2 Axisyn_metric ..........................

vii

ix

XV

1

,)

a

4

5

6

8

11

11

la

14

14

15

NASA/TM--2001-210716 iii
3.

4.

.

.Numerical Methods .............................. 1,

3.1 I)redictor-('orrector Schenms ..................... 17

3.1. l Mac('ormack 2-2 Scheme ................... 18

3.1.2 Gottlieb-Turkel 2-4 Scheme .................. 18

:3.2 Runge-Kutta S('hen3es ......................... 19

3.2.1 Standard Scheme ........................ 20

:3.2.2 Low Dispersion Scheme .................... 20

3.3 Fourier Analysis of Numerical Schemes ................ 22

:3.4 Filtering ................................. 24

1 D Error Analysis .............................. 31

4.1 Problem Descril)tions ......................... 32

4.1.1 Problen_ 1 ............................ :32

.1.1,2 Problem ") 32

4.1.3 Prol)lem 3 ............................ :13

4.2 Results ................................. 33

4.2.1 Problem 1 ............................ 3:3

4.2.2 Problem ') 35

4.2.:3 Problem 3 ............................ :35

Large-Eddy Simulation ............................ 49

5.1 Filtered Equations ........................... 50

.5.1.1 Continuity Equation ....................... 51

,5,1.2 Momentum Equation ...................... ,51

.5.1.3 Energy' Equation ........................ 52

5.2 Sub-Grid Scale Modeling ........................ .5:3

5.2.1 Momentum Equation ...................... 53

5.2.2 Energy Equation ......................... 5:1

5.3 hnplementation ............................. 55

Flow Solver . ................................. 57

6.1 Generalized Curvilinear (',oordinates ................. 57

6.1.1 ( :oordinate Transformation .................. 58

6.1.2 Chain 1Rule Form of the Governing Equations ........ 59

6.2 Time Stepping ............................. 60

6.:3 Treatment of the Viscous Terms .................... 60

6.4 Boundary Conditions .......................... 61

NASA/TM--2001-210716 iv
.

,

(;.5

6.6

6.7

6.4.1 Inflow and Outflow Boundaries ................ 61

6.4.2 Solid Surfaces .......................... 65

6.4.3 Other .............................. 67

Comlmter Resource Requiremenls .................. 71

Parallel hnplementation ........................ 71

Validation ................................ 75

6.7.1 Laminar Flat Plate. ...................... 76

6.7.2 Supersonic Wedge . ...................... 77

•6.7.:_ Supersonic Cone . ....................... 78

Computation of a Nozzle Flowfield ..................... 95

7.1 Description of the Nozzle ....................... 95

7.2 FlowIield Statistics ........................... 96

7.2.1 Time Averaging ........................ .%

7.2.2 Turbulent Statistics ...................... 97

7.2.3 Two Point Correlalions .................... 98

7.3 Axisvmmetric Solutions ........................ 100

7.3.1 Grid Generation ........................ 101

7.4

7.3.2

7.3.3

7.3.4

7.3.5

7.3.6

7.3.7

Specification of Boundary (!onditions ............. 102

Effect of Grid Resolution ................... 103

Comparison of Numerical Schemes .............. 104

Evahlation of Exit Zone Boundary ('ondition ........ 105

Effect of the Sub-Grid Scale Model .............. 106

Evaluation of tile Axisvmmetric LES Sohltion ........ 107

Three-Dimensional Solutions .....................

7.4.1

7.4.2

7.4.3

7.4.4

7.4.5

109

Grid Generation ........................ 110

Specification of Boundary Conditions ............. 110

Presentation of the Time Averaged Flowfield ........ 111

Presentation of the Instantaneous Flowfield and Turbulent

Statistics ............................ 113

Presentation of the Two Point Correlation Data ....... 115

Conclusions and Recommendations ..................... 161

Appendices:

A. Derivation of the Filtered Equations .................... 165

A.1 The Filter ................................ 165

A.2 Continuity Equation .......................... 166

NASA/TM--2001-210716 v
:\.3 ,_lonlentmn Equation .......................... 167

A.I Energy Equation ............................ 168

A.5 Determination of Pressure ....................... 171

Flibliography . ................................... 173

NASA/TM--2001-210716 vi
LIST OF TABLES

Table Page

3.1 Coefficients for fourth-order low-dispersion Runge-I(utta scheme . . 21

6.1 ('oincident grid points in overlap region ................. 70

6.2 Memory required ............................. 71

6.3 ('pIT time required ............................ 72

7.1 Nozzle operating conditions ....................... 96

7.2 Two point correlation location information ............... 100

".3 Effect of grid resolution on axisymmetric solution ........... 103

".4 Comparison of grid spacing to the Kolmogorov scale .......... 104

_.5 Effect of numerical scheme on axisvmmetric solution ......... 105

_.6 Effect of exit zone boundary condition on axisymmetric solution 106

-.7 Effect of sub-grid scale model on axisymmetric solution ........ 107

-.8 \,'elocity statistics at the two point correlation location ........ 115

".9 Curve fit coefficients for two point space correlations ......... 116

_.10 Turbulent length scales .......................... 116

".11 Curve fit coefficients for two point space-time correlations ....... 11S

NASA/TM--2001-210716 vii

LIST OF FIGURES

Figure Page

10 1.1 Jet schematic ...............................

3.1 Error in wave speed ............................ 27

3.2 Error in wave number ........................... 27

28 3.3 Filter response ..............................

3.4 Effect of grid resolution on filter response ................ 29

3.5 Effect of filter order on a five point wave ................ 30

4.1 Solution to linear convection problem .................. 37

4.2 Error of Gottlieb-Turkel scheme for the linear problem ........ 38

4.3 Improved spatial error of Gottlieb-Turkel scheme for the linear problen_ 39

4.4 Error of Runge-Kutta scheme for the linear problem .......... 40

4.5 Efficiency of tile numerical schemes for the linear convection problem 41

4.6 Solution to smooth nonlinear convection problem ........... 41

4.7 Error of Gottlieb-Turkel scheme for the smooth nonlinear problem . 42

4.8 Improved spatial error of Gottlieb-Turkel scheme for the smooth nonlinear problem .............................. 43

4.9 Error of Runge-Kutta scheme for the smooth nonlinear problem 44

NASA/TM--2001-210716 ix
-1.1(}Efiiciencvof tile numerical schemesfor tile linear convectiont)robleln .15

4.l 1 Solution to nonlinear convectionproblem with shock ......... 46

4.12 Error of Gottlieb-murkel schemefor the nonlinear problem with shock 47

4.1:1Error of Runge-I(utta scheme for the nonlinear prol)lem with shock 48

6.1 Modeling an internal objecl using hole points ............. 80

6.2 Pole boundary condition ......................... 81

6.3 Overlap boundary condition ....................... 82

6.4 Parallel processing speedup ....................... 83

6.5 Parallel processing efficiency . ...................... 83

6.6 2D lanfinar fiat plate ........................... 84

6.7 2D grid for laminar fiat plate calculations ................ 84

6.8 Skin fl'iction coefficient for laminar fiat plate .............. 85

6.9 Velocity profiles for laminar fiat plate .................. 86

6.10 Comparison of 2D and 3D Gottlieb-Turkel schenles for laminar fiat plate 87

6.11 (:omparison of 2D and 3I) Runge-Kutta schemes for laminar fiat plate 88

6.12 Mach 2 flow over a 15 degree wedge ................... 89

6.13 2D grid for wedge and cone calculations ................ 90

6.14 Pressure distributkm on wedge surface ................. 91

6.15 Comparison of 2D and 3D codes for wedge calculation ......... 92

6.16 Math 2 flow over a 15 degree cone .................... 93

NASA/TM--2001-210716 x
6.17 Pressuredistribution Oll conesurface ..................

6.18 Comparisonof 2D and 3D codesfor COilecalculation .........

7.1

7.2

7.3

".5

".6

93

94

Two point correlalion schematic . .................... 119

Axisxqnmetriccoml)utational grid .................... 11!)

Boundary conditions fox" axisvmmetric calculations .......... 120

Entropy contours for axisymmetric grid study . ............ 121

Turbulent statistics fox axisymmetric grid study . ........... 122

Entropy contours for axisymmetric scheme comparison ........ 123

-.7 Axisymmetric coml)utatiotml grid with exit zone ............ 123

".8 Entropy contours for exit zone boundary condition comparison .... 124

".9 Entropy contours showing effect of the sub-grid model ......... 125

".10 Density contours for the axisvmmetric solution ............. 126

_.11 Axial velocity contours for the axisymmetric solution ......... 126

".12 Radial velocity contours for the axisymmetric solution ......... 127

7.13 Root mean square density contours for the axisvmmetric solution 127

7.14 Root mean square axial velocity contours fox" the axisymmetric solution 128

7.15 Root mean square radial velocity contours for the axisymmetric solution 128

7.16 Turbulent kinetic energy contours for the axisymmetric solution 129

7.17 Nozzle exit detail for axisymmetric solution .............. 130

7.18 Time averaged centerline velocity profile for the axisymmetric LES

solution .................................. 131

NASA/TM--2001-210716 xi
7.19 Time averagedradial profiles of axial velocity for the axisvmmetric

I,ES solution ................................ 132

7.20 Three-dinlensionalgrid .......................... 135

7.21 Axial velocity conlours for 3D I,ES solution .............. 136

7.22 Radial velocity contours for 3D LES solution .............. 1:37

-.23 Azimuthal velocity contours for 3D LES solulion ............ 138

".2-4 Total velocity contours at ,r/l),i = 3 for :31) LES solution ....... 139

-.25 Total velocity contours at :r/Dj = 6 for 3D LES solution ....... 1:19

-.26 Total velocity contours at :c//Dj _- 9 for 3D LES solution ....... 140

".27 Nozzle exit detail tbr 3D LES solution ................. 141

".28 Dilatation contours for 3D LES solution ................ 142

7.29 Time averaged centerline velocity profile for the aD LES solution . 142

-.30 Time averaged radial profiles of axial velocilv for the 3D LES solution 143

-.al Comparison of sixth- and eighth- order filters on aD LES solution . . 146

-.a2 Instantaneous axial velocity contours for the aI) LES solution .... l-_T

".33 Instantaneous radial velocity contours for the 3D LES solution .... 148

-.a4 Instantaneous azimuthal velocity contours for the 3D LES solution . . 149

-.35 Instantaneous total velocity contours at :r/D; = 3 for the aD LES

solution .................................. 150

7.36 Instantaneous total velocity contours at z/Dj = 6 for the 3D LES

solution .................................. 1;51

Instantaneous total velocity contours at .r/Dj = 9 for the 3D LES

solution .................................. 152

7.:17

NASA/TM--2001-210716 xii
".38 Turl)ulent intensity contom'sfor lhc 3D LES solution ......... 153

".39 Turbulent kinelic energycontours for tile 3l) LES solulion ...... 15-4

".40 Ratio of axial to radial _t_,'}_Je_ll h_*eJ}sitv for l]_e 3I) i,ES so)_li(m 155

-.41 \,%locity history for two point covvelatiotls ............... 1.56

-.42 Two point space correlalion coetlicienl .................. 157

-.4a Two point space-time correlation coefficient ............... 158

-.44 Convection velocity . ........................... 159

NASA/TM--2001-210716 xiii

LIST OF SYMBOLS

a

c

(

_t

k

k

[

71i

P

qi

t

U i

Xi

y+

El_U, U_

X,y.2

,l" ,F.O

speed of sound

wave speed

specific hea,t a,t constan( volume

int ernal energy

total energy

thermal conductivity

turI)ulent kinetic energy

turbulent length scale

surface normal vector

pressure

heat flux vector

time

velocity vector

coordinate vector

inner variable distance

velocity components

cartesian coordinates

cylindrical coordinates

NASA/TM--2001-210716 xv
A+

('

('I

('!

('2

('y

( 'p

Di

l),i

D

[)

E.F.G

tt

J

Pr

Pl" t

Q

R

"R

Rc 5

tfet

Van l)riest daml)ing func!ion constant

Smagorinsky model coefficient

Smagorinsky model coefficient

Sut herland's viscosity law coefficient

Sutherland's viscosity law coefficient

sub-grid scale turbulent dissii)ation rate model coefficient

skin friction coefficient,

l)ressure coefficient

sub-grid scale turbulent diffusion vector

.jet diameter

spatial finite difference operator

eigenvalues of spat ial finite difference operator

flux veclors

axisvmmetric source term vector

Jacobian

Prandtl number

turbulent Prandtl number

conservation variable vector

sub-grid scale heat flux vector

specific gas constant

correlation coefficient

Revnolds number base on jet diameter

Reynolds number based on turbulent length scale

Reynolds number based on axial distance

NASA/TM--2001-210716 xvi
,q'ii

7'

[;

l;

OD

1]

P

oij

T

rij

Oi

5,1,(

strain rale tensor

temperature

velocilv

convective velocity

velocity normal to lhe surface

velocilv at angle 0 to the jet, axis

IRunge-K ut ta coefficienl

solution filter coefficient

gunge-I{ul ta coefficienl

ratio of specific heats

Kronecker delta

turbulent dissil)ation rate

similarity variable

Kolmogorov length scale

dynamic viscosity

kinematic viscosity

density

stress tensor

time delay

sub-grid scale stress tensor

two t)oint correlation location vector

two point, correlation separation vector

wave FlUnl t)el"

computational coordinates

NASA/TM--2001-210716 xvii
__.k filter width

_! tillwstep

__k.r spal ia,1 step

subscripts

a_'g mass averaged

,v exact

j jet

n_ (I,F tllaxilnuln

thin Ininimum

rm._ root mean square

sgs sub-grid scale

U V]SCOllS

0 stagnation/total condition

:x: freestream

superscripts

spatially filtered quantity

Favre averaged quantity

^ resolved quantity

time averaged quantity

' unresolved quantity

" perturbation from time averaged value

: turbulent intensity'

NASA/TM_2001-210716 xviii
CHAPTER 1

INTRODUCTION

The analysis of turbulent compressit)lo jets is a critical lechnology for improving

the performance of advanced aerospace propulsion systems. In bolh the militars' and

commercial sectors a need exists for a be_ter understanding and prediction of jel

engine exhaust plumes and fuel injector flows. This dissertation addresses this need

through the development and application of a. computational method to analyze these

flows, providing a better understanding of the underlying physics.

1.1 Motivation

The survivability of military aircraft depends on their ability to evade the enemy's

defenses. Anti-aircraft svstems employ infrared detectors and heat seeking missiles.

which rely on the high temperature exhaust from the aircraft's engines. The military

has a great interest in low-observable technology to reduce the size and intensity of

this exhaust plume in order to increase the aircraft's survival rate.

In the commercial aircraft industry, reducing the noise generated by the aircraft

has become an important, focus, as increasingly stringent restrictions on noise levels

are imposed near airports. The primary contribution to community noise is jet noise

from the engine at takeoff. The mechanisms by which the ,jet generates sound are

NASA/TM--2001-210716 1
nol well understood. Expanding lhe knowledge of jet aeroacoustics has t)econle lhe

tbcus of sew_ral recent nat, ional progr_mls including NASA's High Speed l_esearch aud

Advanced Subsonic Technology programs.

Both die military and conmwrcial sectors are interested in high speed ramjet and

scralnjeI propulsion systems. These systems are being considered for fasl response

weal)on syslems, high speed aircraft and efficient airt)reathing launch vehicles. The

COml_lete and rapid mixing of the fuel jet with the high speed airstremn is key to

their operation. A better prediclor of the performance of fllel iujectors and their

interaction with the ffeestream flowfield is necessary.

1.2 Physics of Jets

The physics of jet flows is dominated by turbulent motion. Turbulent flows by

definition are unsteady and randomly varying. Turbulence is lhree-dimensional and

rotational with many vortical structures. The scales of the structures vary from the

Kohnogorov microscales [1] to scales nearly' on the order of the jet diameter. The

large scales contain most of the turbulent energy and transport the majority of the

momentum and energy. The energy is cascaded from the largest scales to the smallest.

The smallest scales then dissipate the energy and are isotropic. The Kolmogorov

scales, _. are extremely small. Their ratio to the scales of the largest eddies is

= _673/4 ( 1.1 )

t

where g is the size of the largest scales and Rc_ is the turbulent Reynolds number.

However, Wilcox shows for a typical turbulent flow the Kolmogorov scales are approximately seventy two times the mean free path of the molecules [2]. Thus, the

continuum assumption can be used when modeling the flow at this level.

NASA/TM--2001-210716 2
The t,urlmlent motion in tile jet flowfield is confinedto the lnixing layer, which

is the interface betweenthe flow from the nozzleand the ambient air (figure 1.1 ).

The mixing laver begins at the jet lip and spreads radially with increasing dislanco

downstream. The nozzle flow that has not been affected by the mixing laver maintains

an inviscid character and is termed the potenlial core.

1.3 Computational Methods

An accurate and inexpensive method t,o predict the complex physics of jet flows

would be a great benefit to low-observable, aeroacoustic, and fuel injector technology.

Experimental studies of,jets are an expensive and difficult undertaking and provide

a limited amount of data. Analytical methods in this area are very limited in 1heir

applicability due to the steady-state and low l_eynolds number assumptions that,

usually accompany them.

Computational fluid dynamics (CFD) offers an excellent alternative for jet analyses. In CFD. solutions to the governing equations of fluid motion are obtained using

numerical methods. In principle, it does not require major simplifications to the

equations that limit the applicability and accuracy of the simpler analytical methods.

CFD provides a complete description of the flowfield at specified discretized points.

In this wav it is superior to experimental methods that provide a, limited amount of

data. In general, a CFD analysis is also less expensive than an experimental program.

Two sources of error limit the accuracy of CFD. The first error is the error introduced into the solution by the discretization of the equations. This discretization.

or truncatiom error is a function of both the numerical scheme used to solve the

equations and the computational grid that specifies the discrete locations at which

NASA/TM--2001-210716 3
the equations are solved. The second source of error is modeling error. This error

is tile result of assumplions lnade t,o model the particular geometry or physics of

the l)roblem. ,Modeling errors are introduced in grid generation, boundary condition

specification and within problem formulation and solution process. In fact some of

the largesl and most widespread sources of modeling error occur in formulating and

solving a turbulent flowfield. Three common methods are used to simulate turbulent

flows. They are outlined below.

1.3.1 Reynolds Averaged Navier-Stokes Simulations

The most common approach to simulate a turl)ulent flow is to solve the Reynolds

Averaged Navier-Stokes (HANS) equations. The HANS equations are obtained by

time averaging the Navier-Stokes equations. The contribution of the unsteady terms

in the equations is averaged out and the effect on the flow is replaced bv the Reynolds

stress tensor. In practice, this tensor is modeled using the Boussinesq approximation:

it is replaced by the product of an eddv viscosity and the strain rate tensor. To

further simpli_', isotropic turbulence is typically." assumed. The process of calculating

the eddy viscosity is comlnonly referred to as turbulence modeling.

Turbulence modeling has been an active area of research for many years and a great

number approaches exists. But due to the approximations inherent in the method.

these approaches have failed to produce an adequate simulation of a turbulent jet.

The most widely used turbulence model for HANS simulations of jet flows is the k

model [3]. In this approach, two additional partial differential equations, transport

equations for the turbulent kinetic energy, k, and the turbulent dissipation rate. e.

are solved arid the eddy viscosity is computed from these quantities.

NASA/TM--2001-210716 4
(:orrections to the k-e model to improvojel flowpredictions havebeendeveloped.

Sarkar[.1]and Zeman[5] both introduced corrections to account for compressibility

effects which retard the shear laver mixing in high speed jets. l)ol)e [6] developed a

vortex stretching correction that improved lhe prediclion of round .jets.

In a cooperative efforl among the aerospace commtmity. Barber et. al. [7] conlputed several jet flows using RANS techniques. They found lhat the solutions were

highly dependent on the formulation of the turbulence model and the corrections

used. They concluded that no one model provided adequate predictions over a range

of jet conditions.

RANS simulations are the leasl CPI; intensive method for computing a ttHbu]ent

flow. Typically, the discret.ized equations are marched in t.ime to convergence al

a steady state. Time accuracy in intermediate steps is not necessary and is often

sacrificed for computational speed. In general, the results obtained represent only

the time average of the flowfield. Some unsteady information may be available from

the turbulence model itself (i.e. turbulent kinetic energy is computed in the k e

model).

1.3.2 Direct Numerical Simulations

In theory the simplest and most straightforward way t.o compute a turbulent

flowfield is by performing a Direct Numerical Simulation (DNS). DNS methods solve

the Navier-Stokes equations in a time accurate ma.nner without approximation. In

order to obtain an accurate representation of the turbulence, the turbulent motion

down to the Kolmogorov scales must be accurately resolved in both time and space.

In other words, the grid spacing must be no larger than the Kolmogorov scale, a. The

NASA/TM---2001-210716 5
I_olmogorovscalesvary with tile lurbulent 1Reynolds numl)er as shown ill equation

(1.1). In order to resolve just one large scale eddy we would need ± grid points in

t¢

each direclion, a total of R( 9/4 points. :\nd since the time slep is linearly related lo

the grid size through the ('ourant Friedrichs Lewev (CFL) numl)er [8]. the cosl of the

simulation is on the order of R:_. The amount of computer memory available limils

the size of the comt)utationa] grid that can be used and the computer's processing

speed gives a practical limit to the number of time steps possible in a simulation.

Therefore. based on available computational resources. DNS is limited to very low

t¢evnolds nunll)er flows.

One way to hell) alleviate this iRevnolds number limitation on I)NS is the use

of high order numerical methods [9 15]. These methods have the ability to accuratelv capture small scale structures with fewer grid points than possible with traditional second-order accurate codes. This increased accuracy does come at the price

of increased computational expense and the trade-off between the two has not been

sufficiently investigated.

Freund [16-18] has computed Math 0.8.0.9, and 1.92 jets using a DNS technique.

These jets had Reynolds numbers of S00, 3600, and 2000 respectively. This work has

produced excellent agreement with experimental data and provided great insight into

the flow behavior. But the Reynolds numbers of the jets are orders of magnitude

smaller than any nozzle of practical interest.

1.3.3 Large-Eddy Simulations

A compromise between the approximations necessary in RANS and the computa,

tional limitations of DNS is Large-Eddy Simulation (LES) [19-23]. In LES the large

NASA/TM--2001-210716 6
scalesof turbulent mot.ionare simulateddirectls' in the Navier-Stokesequations,but

the small scalesare modeled ill a tnanner similar to t:{ANS turbulence modeling.

Because the larger scales carry mosl of the momentum and energy, computing lhem

directly should increase the simulation's accuracy. And. since the small scales are dissipative and isotropic, modeling lhem using a simple eddy viscosity approach appears

viable. The scales are separated by spatially filtering the Navier-Stokes equal ions.

This filtering process replaces the equalions wilh a set of resolved (large scale) equations of motion that contain additional tmresolved (small scale) terms that musl be

modeled. The size of the scales that are resolved and that are modeled is determined

by the width of the filter, ._k. which is on the order of the grid cell size.

The key to an accurate LES computation is the model used to ai)proximate the

unresolved or sub-grid scale t.erms. There are many forms of these sub-grid scale

models. The simplest and most popular model, the Smagorinsky model [24] is similar

to Prandtl's mixing length theory [25]. Germane et. al. [26] developed a very successN1 model which dynamically adjusts this constant based on the local flow conditions.

Others [27 30] have modified and improved on these basic ideas.

As with DNS calculations, high order schemes are important for proper and efficient resolution of the large scale structures. In addition, most sub-grid models scale

the eddy viscosity by the square of the filter width, which is directly related to the

grid size. In second-order accurate numerical schemes the truncation error is also a

function of the square of the cell size. LES performed with a second-order accurate

code would have a truncation error and a sub-grid model of the same magnitude.

And because they both scale with A 2 it would be impossible to separate the two by

NASA/TM--2001-210716 7
performhlg grid refinement. Therefore it is crilical lo usea higher order ac('urale

mNnericalmethod when performing LES computations.

Previous stltdies [:/1 3:3] have applied LES codes to high llevnolds lmmber jel

flows. In each of the cited cases, agreement with experimental data has been poor.

:\ common problem wilh these studies is the overprediclion of the length of lhe jet's

potential core. Each solution underpredicts the large scale structures and l'esultanl

turbulent mixing. In a,1] three sludies the computational domain began at lhe nozzle exit and a "'top-hat" inflow velocity profile was specified. The nozzle geometry

and influence of the ,jet lip was not modeled. Mankbadi et. al. [31] also artificially

perturbed the inflow conditions to in an attempt to excite the,jet and increase mixing.

1.4 Discussion of the Present Work

The objective of the present work is to accurately predict a high Reynolds numl)er

turbulent nozzle flowfield using computational methods and gain an improved understanding of the ,jet's behavior. A completely new analysis code has been developed

for this purpose. To meet the goal of accuracy, the sources of error present in any

CFD analysis must be minimized. These errors, discretization error and modeling

error are thoroughly investigated.

Discretization error is addressed through the examination of an existing high-order

numerical scheme and the development of a new scheme. (:are is taken to consider

both temporal and spatial accuracy conjointly. Both schemes are then examined

analytically through their truncation error and experimentally through numerical

experiments to ascertain their performance in solving representative problems. Grid

refinement is also examined to show how increased resolution affects the solution.

NASAfrM--2001-210716 8
The largest contribution to tho modeling error, modelingof tile turbulent slructures, is examined fully. The La.rge-Eddy simulatio, lecb1_ique was selecled as the

most I)rOlnising method to use. The LES equatiotis are derived in detail and the

resulting sub-grid terms are examined. Models fovthe dominate sul)-grid terms are

implemented in the flow solver. Boundary co,tilt ions and other modeling issues along

wilh correct code ilnplen_eiltalion are tested using simple validation cases and a representat.ive jet problem.

[:sing the knowledge gained from testing the codes on simpler problelns, a turbulent compressible round .]el is simulated. The jet has a 1Reynolds number of 1.2

million and an exil. Mach number of 1.4. A complete description of the jel flowfield

including turbulence information is presented. The time averaged flowfield is compared to experimental data to ascertain the accuracy of the simulation. Instantaneous

flowfield data and turbulent statistics lend insight in{o the complex hehavior of the

jet. Correlation of velocity. _ signals in the jet's mixing laver help quantify the large

scale structures present.

NASA/TM--2001-210716 9
Potential Core ______------_.....

.............................................. Z'2 --: :: ...............

 ozzle

Mixing Layer

Figure 1.1: Jet schemat, ic

NASA/TM--2001-210716 10
CHAPTER 2

GOVERNING EQUATIONS

The equations governing the flow of the unsteady compressil)le jet are presenled

below. The fluid is assumed t,o be a conl innunl, which implies that the smallest scales

of interest are much larger than the scales of inolecular motion.

2.1 The Navier-Stokes Equations

The most general set of equations considered are the three-dimellsional NavierStokes equations. As presented below, they express the conservation of mass. moment.urn, and energy for an unsteady compressible fluid in tensor form using cartesian

coordinates (x,y,z).

The continuity equation expresses the conservation of mass.

Op Opu,

0-t + 0xi -- 0 (2.1)

The momentum equation expresses Newton's Second Law. It. relates the time rate

of change of momentum to the forces applied.

Opui Op'ttiuj Op OcYij

o--7-+ ox----/-+ o - j 12.2

NASA/TM--2001-210716 11
]'he stress l ellsor iS defined as

and the strain rate l ensor is

cs,.j= -_t,_L;5,'_._, + 21t,q,j (2.3)

5',, = _- _ + O,r.;/ (2.4)

Sulherland's I,aw is used to model the viscosity

(71 T _

P = ('2 +

The constants for p expressed in English units ( ft.s ' are (7'1

C_ 2 = 19t';.6. -5 •

(2..'3)

= 2.27 x [0 -s and

('onservation of energy expresses the first law of t hermodvnanfics. It relates the

time ra,t,e of change of energy to the a.mount of heat added and the work clone.

Ope:t Ltpuie¢ Ouip Ottjdri) Oqi

- + -- + - (2.6)

Ot O:C i O3"i (-):Tri _),C i

The total energy of the fluid, _t is defined from the internal energy, _ = c,T

e, = _ + 1⁄2","i (2.7)

The heat flux is represented by

q_ = (2.8)

():2" i

The system of equations is closed using the equation of state for a perfed gas.

p = pRT (2.9)

NASA/TM--2001-210716 12
2.2 AMsymmetric Formulation

A solut ion to the t hree-dimensional Navier-Stokes equal ions requires a large amounl

of compuler xnemorv, storage and ('PU lime. Taking advant.age of the symmetry of

a 1)roblem is a conlmon technique used in ('FD 1o reduce the COml)uler resollrCes

required for an analysis. Since the geometry of the round nozzle is symmetric aboul

its ceid.erline, it may be possible to model the 1)roblem using the axisvmmetric form

of the Navier-Stokes equations. In order for this assumption Io be valid the flow field

must also be symmetric about the cenlerline.

The axisvmmetric form of the Navier-Stokes equations assumes that there are

no gradients or velocity components in the circumferential direction of a cylindrical

coordinate system. The symbols x. r. and 0 represent the axial, radial, and circumferential directions. The corresponding velocities are represented by u, t,, and w. The

resulting equations are presented below.

Continuity

Axial momentum

i)p _)pt_ I t)prt,

0--7÷ _-_x + r 0r -0 (2.10)

Opu Opu 2 10pruv Op i')_.x 0_,, 2 r 0 { _,

O---i-+--g-f-x+ + - + + crx_ - ,p-, (911) ," 0," O,r 0. _ _ , r, ""

Radial momentum

Opv Oput, 1 ?)prt, 2 i)p Oar_

--+--+ + -

Ot Ox r Or &" Oa,

+ Oer,._ 2 v "2 ?) l_g v7 \) -aT + - - -

('2.12)

NASA/TM--2001-210716 13
']'h¢' com[:)onent,_ of the stress teilsor are

o'.,..,. (40, 2 0'_, ) -- I* :30a" :] 0/'

( 4 7)v "20 u ) o,.,, = t/ :3 Oz' 3 Ox

o-,.,.= 11 ,_r' + O.r/

[ (o,,

(2. l:}a)

(2.1:lb)

(2.13c)

(2.1:/(t)

Energy

t)p_ t Opt t" I Op_:trc

--+--+

Ot Oa" r

Opu 1 Opr'v O

Ol' + _ + _" 0_" - O.r (uo_._.+ c,_.,. - q_.)

O

+ Or (uo>,. + vo,.,. - q,.) + '.o'x,. + vo',.,.

- q"- ._'7 - "_ U'7 J - "_ ,,:3 ,. /

(2.14)

The corresponding hea.t fluxes a.re represented by

I, OT

q:r _ -- _ -.

Oil_"

(2.15a)

2.3 Flux Vector Forms

OT

= -- (2.15b) % -k Or

For CFD, it is helpful to recast the equations in flux vector form. This puts all

the equations the salne form for easy discretization and solution.

2.3.1 Three-dimensional

The flux vector form of the three-dimensional Xavier-Stokes equations is

OQ OE OF 0(7 OE_. OF_. OG,, (2.16) o--f+_+_ + o:- ox +_+ o---2-_

NASA/TM--2001-210716 14
whPrc

P

p_

p_

o(

E

F

pu 2 + P

DUP

D u tu

pz'

pc 2 + P

DUll'

(p_ t + P)t,

f)U'

_U'U

pw_'

pu ,2 + P

(p_t + P) "'

0

('_ j. g.

(Ta,y

(Tr=

ttcTr,r 4- I!O'x_j -[-WO'3"= -- q,r

0

O'y_.

(Tyy

O'g.:

_[lO'yx _l- I'O'yy "]- II'(T,q.:. -- ([y

0

0".-:.3.

<7:y

0"::

2.3.2 Axisymmetric

The axisymmetric equations are similar to the three-dimensional form with the

addition of a source term. When the source term H is omitted, the equations simplify

to the two-dimensional/planar form.

0(2 OE OF OE,, OF,

o-7+ +N +H- ox

('_).17)

where

Q

NASA/TM---2001-210716 15
F

pu," + p ]

I

,,,, ]

pt'u [

(t'_ t + l') vJ

Et I " ]

('_,r,i-

('_3"7-

IIda,J, q:- t'O',rr -- qa'

[ ° 1

O'r. r

I_FF

UO"r. r .-J-- 1,0".. r -- q,..

tt = I

pl, '2

H_

NASA/TM--2001-210716 16
CHAPTER 3

NUMERICAL METHODS

The governing equations are solved using two explicit finile difference methods.

a MacCormack tyl)e predictor-corrector method and a Runge-l(utta method. Both

types of schemes are in wide use today for fluid c[vnamic and acoustic analyses of

nozzle flows. Each scheme ha.s an associated set, of strengths and wea,knesses which

will be discussed later in this document.

The numerical schemes will be 1)resented in terms of a one-dimensional model

equation

Oq Of

0_ + _ = 0 (:3._)

where f = .f (q). Extending the schemes to the Navier-Stokes equations (2.16 & 2.17)

is a straight-forward matter.

3.1 Predictor-Corrector Schemes

MacCormack developed a two-step explicit finite difference method [34] which is a

variant of the Lax-Wendroff scheme [35]. This method is very easy to understand and

implement and hence has become very popular in the CFD comnmnity. The scheme

is very robust and requires only two storage locations for each dependent variable

NASA/TM--2001-210716 17
(2-.Nstorage).As a consequenceseveralvariants of the MacCormacktechniquehave

])eelldeveloped.The original techniqueand a higller order varianl will be discussed.

3.1.1 MacCormack 2-2 Scheme

Mac('ormack's original schenle is second order accurate in both tinge and space.

That ix to say the truncation error of the scheme is proportional to both the time

step a,ld spatial step to the second power.

The first stage, or predictor step. COmlmt.es the solution at an intermediate time

based on the solut, ion a.t the previous time step. It. uses a one-sided forward difference

for the spatial derivative.

,, At

q7 = q,, _X:c 't''__+_- f/.) (:l._a)

The second stage, or corrector step. COlnputes the solution at the end of the time

step based on the solution at the intermediate time. It uses a one-sided backward

difference for the spatial derivat;ive.

q_+, = 7)1 [.(ll _- q7 -- ._kt '.__i-"'7 (f7 --fi--1 . )] (:],21,)

The leading truncation error terms for MacCormack's scheme are

(At) "203q _ (A,r) 2 Oaf _ At (Aa,)2A04 f

(:_.:_)

6 OU 6 Ox 3 24 Oa"_

where f has been linearized by f = Aq. The error associated with the temporal and

spatial terms are dispersive in nature, while the error associated with the cross term

is dissipative. This dissipative term is scaled by the time step.

3.1.2 Gottlieb-Turkel 2-4 Scheme

Gottlieb and Turkel [9] modilied .\lac(?ormack's scheme to be fourth-order accurate in space while retaining the second-order time accuracy. They simply modified

NASAITM--2001-210716 18
tile differencestencil of the spatial derivative to achievehigher accuracy.

. = ..:xt .,, st,, _ 71:') _1, qV (i A.r (-.1i+2 + ..,,+_

(3Aa)

l_l (TF- sf_L_ + ,tL_)] 6 A,r '

(3.4 b )

The leading terms in tile truncation error for the (k)ttliet)-Turkel scheme are

(At)2 ?)3q (-kx)40_q -kt(A'r)2 O4f (3.5)

6 i)t 3 + 30 0x '_ 18 "40x4

In addition to retaining the second-order time accurac v, il is important to note that

the cross term is dissipative and is scaled by .5t (_kx) 2. Bayliss et. al. [:}6] extended

this scheme to the solution of the Navier-Stokes equations.

3.2 Runge-Kutta Schemes

As seen with the Gottlieb-Turkel scheme, all higher order variants of MacCormack's method maintain second-order time accuracy. It will be shown later that one

cannot separate the temporal and spatial accuracy. They are clearly equated through

equation (3.1).

A true fourth-order accurate scheme in both time and space has been developed.

A fourth-order central difference for the spatial derivative is combined with a fourthorder Runge-Kutta time stepping scheme.

Runge-Kutta schemes are a popular family of numerical schemes with higher order

temporal accuracy. These multi-stage schemes can be formulated for any order of

accuracy. The number of stages in the scheme is equal to or greater than the desired

order of accuracy. Two fourth-order :Runge-Kutta schemes are investigated.

NASA/TM--2001-210716 19
3.2.1 Standard Scheme

Tile standaM four-stage fOUllh-order scheme as giv(,n I_v .Jameson [37] is

qO _ (/N

1

qz =q'_ -_AtD(qo)

q_ = q,_ _l_-ktD(q_)

3

1

q:3 = (I'_ --AtD (q_) 2

q,_+l = q'_ -_kiD (q:3) (3.6)

Tile operator D is the spatial finite difference operator. For equat ion 3.1, D (q) would

l)e a fourth-order finite difference stencil for -"_;,. A central difference stencil for .2-/-/,.,.

is used here.

Of -f,+2 + 8.1,+1 - 8f,-1 + fi-2

_ (a.r)

O.r 12(_ka')

This scheme requires two storage locations for' each dependent, variable (2-N storage).

The leading terins in the t,ruIlca,tion error for' this scheme are

120 Ot 5 30 0x 5 30 . Ox 6 (3.8)

3.2.2 Low Dispersion Scheme

Several researchers [38 41] have developed alternative Runge-Kutta schemes that.

have a lower dist)ersion error than the standard scheme leading to greater stability

and accuracy. To accomplish this, additional stages are required. The additional

stages provide a. means to impose the additional constraints necessary to minimize

the error. All of the schemes are based on a general M-stage 2-N storage formulation

NASA/TM--2001-210716 20
given by

dq._ = c_.,dq .... l + AID(q,._1) (3.!_a)

q,,, = q,,-1 + :J,, dq,, (3.9b)

for m = 1... M. ail([ where q0 = qn all([ q.U = qn+l. The coefficient o j is tyl)ically

set to zero for the algorithm lo I)e self-slarling. Again, the operator D is tile spatial

finite difference operator.

Carpenter and t,_ennedy's five-stage fourth-order scheme [38] was chosen for ils

fourth-order a,ccura.cy, low number of stages, and ease of programnfing. The coefficients tot the schelne are given in lable 3.1

stage (Tn) a,, ,(3,,_

] 0.00000000000 0.1496590219993

2 -0.41789047450 0.3792103129999

3 -1.192151694643 0.8229550293869

4 q.697784692471 0.6994504559488

5 -i.514183444257 0.1530572479681

Table 3.1: Coefficients for fourth-order low-dispersion Runge-Kutta scheme

The leading terms in the truncation error for this scheme are

300 0{" + 30 0a "5 :30 ?):r6

The dispersive error term due to the thne step (the first term) is two and one-half

times smaller than the standa,rd fourth-order scheme's error (3.8).

NASA/TM_2001-210716 21
3.3 Fourier Analysis of Numerical Schenms

In lhe field of high order numerical methods for [)NS and I,ES. the proper resolution of tlw waves ])reselll, in tile flowtield is critical. :\ popular method to assess

lhe abitily of a numerical scheme to resolve waves is the Fourier analysis of lhe semidiscrelized equation (spatial discret.ization only) [42].

The semi-discretization of equataion 3.1 is

Oq,

-D(%) (3.11)

Ot

where q,, is tile discrete solution.

number _,,. which has the form

\Ve choose a, sinusoidal trial solution with wave

q.(_:.t) = _,(_.,,t)¢ _'_'' (3.12)

\Ve substilute 3.1'2 into 3.11 and solve for v(_,'.l). The resulting solution of the

equalion is

q,_(w,/) = v(_.', 0)eR°[D_'lqe i_I'''-C'') (3.13a)

with

c'= -!Im [/)(_:)l]

where D(w) are the eigenvalues of D. Tile exact solution t,o equation 3.1 is

q(w', t) = *,(_', 0)e i_'('-_'}

(3.13b)

(3.14)

Bv comparing the numerical and exact solutions we can see how the numerical schelne

affects the propagation of the wave. The exact solution (equation (3.14)) shows that

NASA/TM--2001-210716 22
amplitude of tlle waveshould not change. However,in tlle ntunerical solution t ll(,

amplitude is a function of time (equation (3.13a)). Tlierefore.

if Re [/)(_,)] = O, then the scheme is conser,'aiive

if Re [b(_')] < O. then the scheme is dissipaiive

if Re [D(.:)] > O. then the scheme is unstal)le

Similarly. comparing the numerical wave speed c" to tlw exact wave speed c we find

if -- = 1, then the scheme ha.s no phase error

C

if' c[ # 1, then lhe scheme inlroduces a phase error

C

To achie,'e a conser,'ative scheme (Re [b(_')] = O) the matrix D llltlst ])e antisymmetric. Central difference schemes, such as those used here, satisfy this criteria.

I.Tpwind schemes, which are dissipative, do not have an anti-symmetric matrix and

were not considered suitable for this study. The waw" speed relative to the exact

speed for each scheme considered here is

c" sin( A.r)

C*

c .JAx

, for second-order spatial accuracy

(3.t.sa)

, for fourth-order spatial accuracy

(3.]5b)

Figure 3.1 compares the error in the wave speeds for second- and tburth-order

accurate central difference operators, D. The error is a function of the wave numl)er,

,_,, and grid resolution, Am. From this Fourier analysis perspective the numerical

scheme can be seen to act as a spectral filter of the exact solution. ]'his error is

sometimes plotted as numerical wave number versus exact, wave lmmber (_,,'.Xa" versus

NASA/TM--2001-210716 23
_'_Xx.tigure3.2). It is clearthat higherorder schemescanresolvehigher wavenunl])er

waveswith lessgrid resolution. But, heir,her scheme('all adequatelyresolvevery high

Wave lluin])ers, evell with very fine grid resolution.

If one determines an acceptable error in wave sl)eed for their calculation, the

maximum value of ,'.Sx that provides this error can be related to the number of grid

points per wavelength necessary to accurately resolve a wave [43].

27r

N - (3.16)

The wave will be properly resoh'ed if it is captured by at least N points. This "'points

per wavelength" has become a popular measure of a schenle's "'goodness". But, this

ineasure does not consider the error due to time discretization and the computalional

cost of the scheme.

3.4 Filtering

Both numerical methods introduce a. dispersive error into the solution. This error manifests itself in two ways. First, as shown above, waves that are that are not

resolved with enough grid points have their speeds altered. Secondly, new high wave

number waves are introduced into the solution. If unchecked, these errors can grow

arid cause the analysis to become unstable. To damp out these waves, an artificial

dissipation term is typically added to the equations [44 47]. All these "classical"

techniques involve the use of a user specified constant that scales the effect of the dissipation. The value of this constant is somewhat arbitrary and varies widely between

cases.

In this study a solution filtering technique is used. This technique can be regarded

in two ways. In the classical sense, solution filtering simply adds a dissipative term to

NASA/TM--2001-210716 24
the equationsto damI) the unwantedwaves.Froma Fourier analysisperspective,this

techniquefilters the high wavenumber componentsout of the solution. A prol)erly

selectedfilter removesboth tile numericallyintroduced wavesand the poorly resolved

waves,leaving only the portions of the flowfield that are accurately modeled. Several researchers have developed solution filters for CFI) including [10,42.48.49]. The

explicit filters of Kennedy and Carpenter [48] were chosen for their ease of implementation and clear relation to the unfiltered solution, The filtered solution is simply the

unfiltered solution plus a dissipative term.

cti = qi - oD(-Xa') 2_ O:r_,,

(-ll" and J_ :

_,vhel'e OD = 22,---'7--

(:/.17)

1.') I{ennedv and (a.zpent I developed a family

of filters with corresponding boundary stencils for tt = 1... 7. These filters are

implemented in the analysis code. By choosing the order of the filter. 2'tt, to be larger

than the order of accuracy of the numerical scheme, one can insure that the filter

does Ilot overly influence the numerical solution. The response of the filters is shown

in figure 3.3. As the order of the filter is increases the cut-off wave number of the

filter is increased and a greater portion of the domain remains unmodified.

Examples of the effect, of solution filtering are shown in figures 3.4 and :t.5. A

sine wave on the domain 0 _< ,r < rr was filtered 100,000 times (a typical number of

iterations for a jet calculation). To illustrate how increased grid resolution reduces the

effect of the filter, sine waves made up of 5, 15, 95, 35, and 45 points (_'Aa" = 1.257,

0.4189, 0.2513, 0.1795, and 0.1;196) were filtered using fourth-, sixth-, and eighthorder fil_ers. Figure 3.4 shows that increasing the resolution of the wave reduces the

NASA/TM--2001-210716 25
etfecl of the filter. Poorly resolvedwavescan be completelyeliminated (lhe fourlhorder filter completely removesthe ,5 point resolved wave) Figure 3.5 shows how ihe

resolved wave number increases wilh the order of the filter. :\ wave resolved by five

points was filtered with fourth-, sixth-, eighth-, lenth-, twelfth-, and fourteenth-order

filters. The wave (_:__X,r= 1.2:57) is completely removed by the fourth order filter, but

damped less than 10 percent by the tenth-order filter, and maintained exactly by lhe

fourteenth-order filter.

NASA/TM--2001-210716 26
1.2

0.8

0.6

0.4

0.2

0

0

.... I .... P ........ i ........ i ....

,'" i I--exac' .I

....................._,....................._._......... "_-, *, :. ..... 2nd order ! ........

.................?.............. - .... 4th order |

N_ ", i i

"< i

L ", i_. i

..................................... _.................... ,................ ,_............ _._,............. +....................

! I i " "x

.... i .... i .... _ .... i .... i .... :",,...

0.5 1 1.5 2 2.5 3 3.5

Figure 3.1" Error in wave spee(l

3

3.5

2.5

2

1.5

0.5

"'''l ................ ! .... ! ....

-- e_act I i ;

..... 2nd order _ .... "? .... _ " " --7_'_

.,,oral,,/ i

t I i i Ji '

.....i, ' "'. Ct-----

0.5 1 1.5 2 2.5 3 3.5

(o2_x

Figure 3.2: Error in wave number

NAS A/TM--2001-210716 27
1.2

1

0.8

0,6

0.4

0.2

0

................ ! .... ! .... I ....

i \ _,"_\\\ '_ ......'

0.5 1 1.5 2 2.5 3 3.5

_ZkX

Figure 3.3: Filter response

NASA/TM--2001-210716 28
o"

0.8

0,6

0.4

°I

.o2

0 0,5 1 1.5 2 2.5 3

(a,) fourth-order filter

o-

/

-- Exact •

--o-- 5 Points .

1 l- ....................... ........... _ ............. _15Points _

t _ - _ _25 Points .

[ .,_i'_-_-_"_'-._"_"l_ _ 35 Points

_f_'_ _"_-'_ --+-- 45 Points "

0.8 .............................................. ............... -

0.6 .....................-., ................. _..._

0.4 ...................

°I

-0.2

0 0.5 1 1.5 2 2,5 3

(b) sixth-order filter

Figure 3.4: Effect of grid resolution on filter response

continued

NASA/TM_2001-210716 29
Figure 3.4: continued

1.2

1

0.8

0.6

0.4

0.2

O,

-0.2

0

I +5 Points |J

............................................. i...................... _ ....... ] _ 15 Points U

.................. / + 25 Points

_ | _ 35 Points |7

i

0,5 1 1.5 2 2.5 3

(c) eight-order filter

1.2 L' '.........'...... '" ,-,--r-,-, ................i'_',_' r" ' ' ..........r'"" """T

' -- Exact

_- ---o--- 4th order

l- _ 6th order

1 _-............................................ _.............._'"-'"t ........... ---e--- 8th order

I" f j"_'._. _ + lOth order

I- J .J'j ._ "% _ 12th order

F // _Jl j _'_ N _ 14th order

o._ o.6 F........................ ................... i..-..___ ................................ i.....

o.24 I_io. ................................ ;.-

oT o " " i

-0.2

0 0.5 1 1.5 2 2.5 3

Figure 3.,5: Effect of filter order on a five point wave

NASA/TM--2001-210716 30
CHAPTER 4

1-D ERROR ANALYSIS

The t)ehavior of the numerical schemes presenled in ('hapter 3 is examined for

three one-dimensional model problen>. Truncation error, efficiency, and the consequence of disparate temi)oral and spatial accuracy are discussed. The schemes are

used t,o solve the one-dimensional inviscid convection of a gaussian pulse.

The accuracy of these schemes is commonly expressed separately in terms of spatim and temporal accuracy. A great deal of effort, has focused on increasing the spatial

accuracy of a numerical scheme without regard to the temporal accuracy,. This chapter examines the behavior of two commonly used schemes in terms of truncation error

and computational cost on a model equation. Particular attention is paid t,o the truncation error of the schemes and how the spatial and temporal errors affect the overall

order of accuracy of t}le scheme.

Three one-dimensional problems were used to test tire accuracy and efficiency of

the schemes. All three are based on the model equation (3.1).

NASA/TM---2001-210716 31
4.1 Problem Descriptions

4.1.1 Problem 1

The first problen] is the linear convection of a gaussian pulso, whero q = f = ,.

(i'll (-)?1

0-7+7=0 (4.1)

Tile initial condition is given by

.)2 .(.,,.0) = .0(_')= ,4,-_12t(5 (4.2)

Tile domain is -20 < x < 450 and the solution is run for 0 < t < ,i00.

solution to this problem exists and is given by

AJI (_xacl

u,o-(:r,t) = .o(:r - t)

4.1.2 Problem 2

(4.a)

The second problem is nonlinear where q = , and f= _t_l2.

&' 0-7+_ ()( 7v2 ) =0 (4.4)

The initial conditions and sinmtation time are set so that a smooth final solution is

obtained. The initial condition is

(4..5)

The domain is -50 < x < 50 and the solution is run tbr 0 < t < 100. A numerical

approximation to the exact solution is obtained using the method of characteristics.

NASA/TM--2001-210716 32
4.1.3 Problem 3

The third problenl is also nonlinear (equation (4.4)). but the initial condilions and

time of 1he simulation are set so lhat the pulse is allowed to coalesce into a shock.

,,,(_,,,o) = ,o(:,, ) = '-_-'_)( _/ (4._;)

2

The domain is -50 < a,' < 50 and the solution is run for 0 _< t _< 200. lake pl'oblena

'2, a. numerica.l api)roximat:ion to t,he exact solution is obtgined using the met.hod of

characteristics. The location of lhe shock is then fitted in the exact solution using

Whitham's area rule [501.

4.2 Results

All the one-dimensional calculations were run on an Apple Macintosh Powerbook

G3 computer with a 266MHz PowerP(' G3 processor.

The standard scheme was used for all the Runge-Kutta results in this section.

The error of the numerical scheme was measured by the 12 norm, which is computed

as follows

1

12 = x,,,_-- z,_/,_ }-_(u,:- u_.) 2 : (4.7)

i= 1

For each problem, both numerical schemes were run to determine the maximum

stable time step. Then, each scheme was run at a number of different x'a.lues of at as

the spatial step was halved.

4.2.1 Problem 1

A sample solution to the linear problem for both numerica.l schemes is presented in

figure 4.1. The Runge-Kutta scheme shows better agreement with the exact solution

NASA/TM--2001-210716 33
lhau the (-;ot,t lieb-Turkel schenle. Error from lhe (;ot tlieb-Iurkel scheme is plotted

wwsus spalaial step and time step (figure 4.2). The l,ernfinal slope of the line. p

indicates the order of accuracy of the scheme and is listed in the plol legend. The

slopes in figure 4.2(a) show that the error is of second-order accuracy for the varyhlg

spatial slep. Only where the lime step is much smaller than tile spatial slep does the

schelne behave as a tburth-order schelne. Also, the error is dependent on the time

step chosen. This is due to the third term ill tile truncation error (equal, ion (:{.5)).

('learly, this dissipalive error is signiticanl. Figure 4.2(b), sllows tile error 1)ehavior

with the time step. Again, the scheme is only second-order accurate except when the

t,ime step is ve|'y small coral)areal to the spa, tia[ step. Iu order to obtain fourth-order

behavior fi'om the Gotllieb-Turkel scheme, the error due to the time step must be

reduced at the same rate as the error due to tile spatial step. In order to accomplish

1

this, lhe time step must be reduced by a factor of 7? as the spatial step is reduced by

}. Figure 4.:3 verifies that when adjusting the time step in this manner fourth-order

accuracy is obtained.

[{esults for 1Runge-Kutta scheme, on the linear prob]enl (figure 4.4); indicate that

the scheme is truly a fourth-order accurate scheme. In addition, except for the highest

value of _ where the thne error term dominates, the error is independent of time

step.

The efficiency of a numerical scheme can be seen by comparing the/2 error against

the linle required to obtain that error level [51]. Both schemes are compared at their

lnaxitllUlTl stable time step. Figure 4.5 shows that the Runge-Kutta scheme is superior

to the GottIie]>Turkel schem_e. Runge-Kutta achieves an equivalent error level in

about an order of magnitude less time than the Gottlieb-Turkel scheme. This result

NASA/TM--2001-210716 34
is due to lhe lower truncation error and larger allowable tinie step of the RuugeI{utta schenw. Reducing the time step of the Gott,liel)-Turkel schemeto maintain

fourth-order accuracy increasesthe accuracy but also increaseslhe coml)utational

cost.

4.2.2 Problem 2

Sample solutions to the smooth nonlinear problem for bolh numerical schemes are

presented in figure 4.6. Because a short sinmlation lime was necessary to maintain a

smooth solution, the differences between the exacl solution and numerical solutions

are not discernible on the plot. Error data for the smooth nonlinear problem are shown

in figures 4.7-4.10. The effects of round-off error, the reduction in accuracy/slope at

small x.t and Ax, make it difficult to report a terminal slope, instead the maxinmm

slope was used. The trends for the nonlinear analysis are similar to the linear case.

However. the benefit for using the l:{unge-Kutta scheme is somewhat less. But, at

large time steps, for maximum efficiency, the Runge-t(utta scheme is superior.

4.2.3 Problem 3

Solutions to the nonlinear problem with a shock (figure 4.11) show that both

schemes produce high frequency oscillations near the discontinuity. This problem does

not provide useful error information and illustrates the limitations of finite difference

schemes. The presence of the discontinuity reduces the accuracy' of all schemes to

first-order (figures 4.12 & 4.13). Results from the Gottlieb-Turkel scheme indicate

that the third term in equation (3.5) provides a damping effect, on the numerical

oscillations that occur near the shock. This damping allowed the scheme to obtain a

converged answer. Figure 4.12 shows that larger time steps yielded lower errors for

NASA/TM--2001-210716 35
this problern. Tile dissipative term in the truncation error is scaled with the time step

and la.rger time steps resulted in greater dissipa, tion of tile non-physical oscillal ions.

To oblain useful solutions to this prol)lem, the }{unge-Kutta scheme required

sonic form of artiticial dissipation. The solution filtering technique of l(ennedv and

(!arpenter [18] was used. Fillering adds to the solution an additional dissipative lernl

lhat removes high frequency numerical oscillations while leaving the low frequency

physical oscillations untouched.

The order of the filter. 2/_ in equation (3.17), determines the magnilude of t,he

dissipation added and the range of frequencies which are damped. If the order of the

filter is la,rger i,han t,he order of the truncation error of lhe numerical scheme, lhen

the filter should have a negligible effec! on the solution error. Increasing the order of

the filler increases the range of fl'equencies which are tlntouched by the fillet'. Figure

4.13 shows that the fillers remove the non-physical oscillations from the solution and

greatly reduce the error of lhe scheme. The insensitivity of the error to the order of the

filter indicates that, the truncation error of the scheme is grea,ter than l,he dissipative

term of the filter and thai, the nulnerical oscillations have frequencies higher than

frequency range of the higltest order filter.

NASA/TM---2001-210716 36
0.6

0.5

0.4

0.3

0.2

0.1

-0.1

385

tl ....."<'1 i 1/_'. I

.................. !.................... r........... ;........ t..........._,t ..........................................

i i , i L ! i

i ', / i ', _ _L

I

, i i I J i . • • i .... i ........ i ....

390 395 400 405 410 415

Figure 4.1: Solution to linear convection problem

NASA/TM--2001-210716 37
w

10 _

1 0 °

10'

10.2

1 0 .3

1 0 .4

1 0 "5

10"

o=2ool i ........ ]

--.-_t,_x=o25oo, p=l_/ __i t

--_ _t,_,=o12so,_=19zI .......__ ............ 1

i i

.............

i i

..... ,,Jl ...... i , ....

10 4 10 ° 1 01

._x

(a) Spatial error

10'

10°

10 4

10 `2

10 ..3

10 `4

....................... I ......

J

10 3 10 .2 10 "1 t0 0 101

_t

(b) Temporal error

Figure 4.2: Error of Gottlieb-Turkel scheme for the linear problem

NASA/TM--2001-210716 38
uJ

1 0 °

10-1

10 2

1 0 3

10 .4`

1 0 -s

10 2

i.......

L

10 4 10 0 101

AX

Figure 4.:3: Improved spatial error of (;ot tlieb-Turkel scheme for the linear problem

NASA/TM--2001-210716 39
10 _

I 0°

10 -_

10 2

10 3

10 .4

10 .5

10 .6

10 `2

F F

--_-- _t/_,=2.00,p=4.00 _ 1

....--_-,,it.,, =1.00,p=4.00!........................... _ ............... !J

F I -'- =_x=0.50, p=4.ooI ._"//"

F-I--*- _'"_"=o25,p. 400!----/-_ .........

i i;iii , i iiiii 

t0 _ 10 0 101

,,Xx

(a) Spatial error

10'

10°

10._

10 .2

1 0 .3

1 04

105

1 04

10 .3

i ........

""!

&idhx=2.00, p=4.001_ _J f._

,&t/Ax = 1.00, p = 4.00

>,t/hx = 0.50, p = 4.00

+ At,'Ax = 0.25, p = 4.00

/

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

10 2 10 4 10 0 10 _

T'1

_t

(b) Temporal error

Figure 4.4: Error of Runge-t(utta scheme for the linear problem

NASA/TM--2001-210716 40
D_

10 _

1 0 °

10 -1

10 2

10 3

10 4

1 0 5

1 0 .6

10 .2

........ = ........ ! ........ ! ........ ! ........ I .......

i

i i i I--'-°T,_X=°._ I

i i i / + GT, ±VAX = f(3.x)I

............................. .........................

:. i _ i

........ t ........ i ........ i ........ i ........ i ........

1 0 * 1 0 ° 1 0 _ 1 02 1 03 1 04

time (s)

Figure 4.5: Efficiency of the numerical schemes for the linear convection t_rol)lem

0.14

0.12

0,1

0.08

= 0.06

0.04

0.02

0

-0.02

-40

: l_analyticall _ i

I-E 7--7

i _ f l i

..... i ........ ?........

-30 -20 - 10 0

.... i .... I ....

17..........................

,. i .... I ....

10 20 30 40

Figure 4.6: Solution to smooth nonlinear convection problem

NASA/TM--2001-210716 41
o

uJ

w_

10 -1

i !

10 "_

10 3

10 .4

1 0 s

10 `6

10 z

I

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

00, p = 1.90

/i./_ [ _ _/'3X = 5.000, p = 1.89

l I 3.000,2.88

i _//Z_.// I _ _w : 0.6s,p=3.40

I 0 _ I 0 ° 0 _

(a) Spatial error

10 "l

1 03

1 0 "4

10 .5

----e--- 3,t/Ax = 6.000, p = 1.90

Atl_.x = 5.000, p = 1.89

_- ±t/._ = 4.000, p = 2.54

At/.3.x = 3.000, p = 2.86

,_t/._x = 2.500, p = 3.20

.MJ&x = 1.250, p = 3.45

At/_.X = 0.625, p = 3.40

10 -6

10-3 10-2 10-_

..................................... :i .__7"........

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

1 0 ° 1 0 _

(b) Temporal error

Figure 4.7: Error of Gottlieb-Turkel schelne for the smooth nonlinear problem

NASA/TM--2001-210716 42
tlJ

10 2

10 .3

10 .4

10 .5

10 6

10-2 10 _

........ I ........ T '

+ Att',x = f(Ax) p =3.34_

./

/

...... I

10 I I 0 °

_x

Figure 4.8: hnproved spatial error of (,ottheb-Turk l scheme for the smooth nonlinear

problem

NASA/TM--2001-210716 43
u3

10 "T

1 0 .4

lO-S

10 6

10 2

....... I ......

....... _-4 I-," ,----_'x- ='12ooo ----o=2,9t__

I--"- _"_=1_0:0_00, p=2.9/

///_._ 5.000,p=3.33 /

___ --.-_,,_=_,00,,_=3.38I__

/// ! I-_-_x- ,_,0p=339/

._ I--'-'_i °82'_=339! _ ...... ! _..I ........

10 _ 100 101

_x

(a) Spatial error

10'

1 0 ,2

10-3

1 0 .4

105

1 0 6

1 04

-_e._ &t,'3.x = 12.000, p = 2.59

+ At/Ax = 10.000, p = 2.91

_ AUAx = 5.000, p = 3.33 __

,.'.'.'.SU.3.x= 2.500, p = 3.38

&t/t_,x = 1.250, p = 3.39 1

_--i i i --

1 0 z 1 0 z 1 0 _ 1 0 ° 10' 1 02

(b) Temporal error

Figure 4.9: Error of Runge-Kutta scheme for the smooth nonlinear problem

NASA/TM--200 I-210716 44
uJ

10 2

10 3

10 "4

10 5

10 .6

10 .2

........ _ ........ ' ........ ' ........ I ........ I .......

RK4, M/Ax 10,0

10 4 10 ° 1 01 1 02 10 3

time (s)

Figure 4.10: Efficiency of tile numerical schelnes for the linear convection problem

NASA/TM--2001-210716 45
0,3

0,25

0.2

0.15

0.1

0.05

..... f -- analytical t .................. _".................. ................ "l"-- .............

! ..... GT,no,lter I ! | i

j ..... G'°'h°rderf"ter/i : _li

--1 ..... RK46thorderf,,tet------7-----_-

-0.05

-20 -10 0 t0 20 30 40 50

(a) Complete wave

0.3

0.25

0.2

0.15

0,1

i / "/

_1 ,2- • ,d

0.05

Imeo.tica, li

- .... GT, no filter _ i

..... GT, 6th order filter ] i

..... RK4, 6th order filte I i

-0.05 .... i .... I .... i .... i .... i ....

3 5 36 37 38 39 40 41

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

i

"i F

(b) Detail of shock

Figure ,I.11' Solution to nonlinear convection problem with shock

NASA/TM--2001-210716 46
I 0°

10 -I

10 2

10 .3

10 `2

I

....... J

10 1

I

_ _1/3.x = 1.25000

---IB-- _t/-_x = 0.62500

__tl',x = 0,31250

_,t/-_x = 0.15625

10 o

L......

10 _

_x

(a) Spatial error

I0 °

10 _

10 .2

10 3

10 .3

_--_ _ ,_t/,&x = 1.25000 i

I .--O-- _t/,_x = 0.62500 1

i _ z_t/,_x = 0.31250 t

1

.... ,,,i ............... l .........

10"2 0 1 100 101

(b) Temporal error

Figure 4.1'2: Error of Gottlieb-Turkel scheme for the nonlinear problem with shock

NASA/TM--2001-210716 47
10'

10°

I0 i

10 2

10 3

10 2

....... T .......

_h orderfilteH 10 "1 10 0 101

(a) Spatial error

10 _

10°

10 -1

10 _

10 3

1 0 2

---,--no,,er'''"l T ........ 1

4th order filter I ) _m

-.--",'-, °rde" fi't°r I _ 1

-'*-''t--. _"o_°rc'rfi't'er' _ J

10 4 10 0 10'

(b) Temporal error

Figure 4.13: Error of Runge-Kutta scheme for the nonlinear problem with shock

NASA/TM--2001-210716 48
CHAPTER 5

LARGE-EDDY SIMULATION

In this chapter the equations for Large-E(ldy Simulation are presented. The modeling of the unresolved (sul)-grid) terlns and implenwntatio_l into the flow sol\'er is

discussed.

The basis of LES is the separation of tile large and small scale turbulent fluctuations. To separate the large (resolved) and small (unresolved) scales the equations

are filtered. A spatial filter G with a filter width _k is used.

f = G(x - g).f(<)d5 (5.1)

• ,5"<,

The overbar represents the resolved, filtered, or large-scale portion of the function.

Commonly used filter functions are a box filter, a Ga.ussia.n fillet', or a spectral cutoff

filter [301. TypicaJly, it is not necessary to know the exact form of the filter function

G, hut only that it. exists. In practice, the solution is not. explicitly filtered. It is

assumed that the numerically computed solution is a filtered representation of the

exact solution. This assumption is justified based on the Fourier analysis results of

sections a.a and 3.4. There. the numerical schemes and solution filtering were shown

to behave as a spectral cutoff filter of the exact solution.

NASA/TM--2001-210716 49
SeveralCOllstraintsare imposedon t].' tiller function.

i) C;(-(,) = (;(_)

2) G(_),/_ = J

4) G(_) is small outside

These constraints are necessary to insure that, the filt.er function will commute with

the derivat ire.

of of

- (s.:_ 0,r 0x

Favre (density) weighting is used in the filtering process. This allows for convenienl

recovery of terms corresponding to the unfiltered equations.

P

5.1 Filtered Equations

The filtering process is apt)lied to the continuity, momentum, and energy equations. Details of this process can be found in Appendix A. The resulting equations

are comprised of resolved and unresolved terms. The resolved terms in the filtered

equations directly correspond, in form, to the unfiltered equations. The additional

unresolved or sub-grid terms are modeled as source terms to the equations.

NASA/TM--2001-210716 50
5.1.1 Continuity Equation

The filtered continuity equation contains only resolved terms. Thus. it corresponds

directly to the unfiltered equation (2.1).

Op O/)h,:

0-7 + &'; - 0 (5.4)

5.1.2 Momentum Equation

Tile moment um equat.ion contains lwo unresolved terms (underbraced) lhal must

be modeled.

O/>gi; ot)a_a i op O&; Or_.i O

57 + o.,., + o_,,,_= 0%-7 _ + _ I_,,,- a_,l

V

tl i t ̧

Three stress tensors occur in (5.5). The filtered stress tensor

(5..-,)

o-,:,_= -- ,{p _i,i.Sxx + -I t,5;.i 15 (_)

The resolved stress tensor

2_ _ -g

&ij -- 5t*(_ij,. _.],. + 2fiSij (5.;)

where/5 = p(T), and

'_ij = _ \ OX i + 0,;t'j ]

The third is the unresolved or sub-grid scale stress tensor, Tij.

provides the effect of the sub-grid scale turbulence on the larger resolved scales.

(5.s)

It is this term that

NASA/TM--2001-210716 51
5.1.3 Energy Equation

The Fax're fillered energy e(luation contains five additional unresolved lernls (underbraced) that must be modeled.

--+ --+

Ot Ox, Ox, Ox, #x, O,c, _ Oxi

i iii

" , =:--- p - (0_- ii,)

+ _ 0,_., " o.,'iI

t I "b

where the tiltered total energy is

(5.9)

e,=_+1⁄2_ (.5.Jo)

and the filtered and resolved heat flux vectors are

OT

q, = -I,,--

OaF i

(5.11a)

where L"= k(T). The sub-grid scale heat flux is given by

(5.1 lb)

(5.12)

The sub-grid scale turbulent dissipation rate is

( = o,j Oxi 6i50xi (5.1_)

and the sub-grid scale turbulent diffusion is

riD,.= _ (5.14)

NASA/TM--2001-210716 52
The eqltalion of state becomes

l, =/)K:_ (5.1.5)

Since the total energy, equation (.5.10). conlains the unknown _, pressure musl be

obtained using the following expression

p (_ - 1) (pgt - _p'u:.u_, 1 ~ ( 16)

5.2 Sub-Grid Scale Modeling

Tile unresolved or sub-grid scale conlribulions to equations (5.9) and (5.5) musl

be modeled. The methodology used here is based on the incompressible sub-grid

scale model of Slnagorinsky [24]. Additional lerms to account for compressibility

were added based on the work of Moin et. a l. [28] and Vreman et. al. [52].

5.2.1 Momentum Equation

The Smagorinsky model is a very popular and widely used model for the subgrid scale stress tensor. It. is an eddy viscosity model where the sub-grid sca.le stress

tensor is modeled as an eddy viscosity multiplying the resolved stress tensor. It was

developed for incompressible flows and has been frequently extended to compressible

flows. The compressible form of the model, given by Moin et. al., was used here.

T,.,.= }c.,n_l,_F_,, .- 2cO_x_lsl(-%- _,_L_,,). (.5._7)

The coefficients C and CI are user defined constants. The second unresolved term in

the equation results from the nonlinearity of the viscous stresses and is neglected.

The final form of the modeled momentum equation is

ODu, 0-7-+ O_i_j o_----f-+ ox_ i)_ 0 [ &ij _ _ii ] (5.18) - oxj (p + _'') _--7

NASAfrM--2001-210716 53
The eddy viscosity is

t,, = ('p-x2l,_l (._.1_))

and ¢ is

(]) _ 2,,, _,2

5.2.2 Energy Equation

There are many different sub-grid inodels for the energy equation and no one

set of models is as popular a,s the Smagorinskv model for the momentunl equation.

This is because there are many forms of the energy equal,ion that resull in different

sub-grid scale terms. It is also because the majority of work in LES has been done

for incompressible flows. The filtered energy equalion (5.9) contains five ullresolved

terlllS.

The first term. the sub-grid scale heat flux is modeled based on Moin's work [28]

O,_ _,, O_

Prt Oxi (5.21)

The model for the second term comes from Vreman [52] and is given by

C3

]'he coefficient (':_ is a user define(/ (onstant. The lasl three terms are typically

considered to be nmch smaller than the other terlns and are neglected.

The final form of the modeled energy equation is

Ot_gt Oi_a/6 O?tip O _,au + _ + pt - --

o_ + o.< + o,r_ - o.,,_ ET,., ..,- _ P,._J + ': (5.2.a)

NASA/TM--2001-210716 54
5.3 Implementation

The' sub-grid scale terms in equations (5.18) and (5.23) and their corresponding models were iml)lemented ill both the two-dimensional/axisymmel tic all([ thl"eedimensional Navier-Stokes flow solvers. It. is important 1o note that lhe modified

definition of tile pressure (equation (.5.16)) must be used.

Because the resolution of scales in a numerical sc[wme is directly related to the

grid resolution, the filter width A was chosen to be a characteristic length of the

COml)utational grid. Since the grid is not uniform, this length varies widely over the

grid. The filter width at each location was defined as the cubed root of the volume

associated with each grid point.

TO maintain a laminar sub-]a._,.'er in wall bounded regions, the effect of the sub-grid

model must be diminished neaz t,he wall. To a.ccomplish this. a Van Driest damping

function [53] was used to scale the effect of the sub-grid terms.

where y+ is the inner variable distance and the constant A + is set to 26.

'The constants for the sub-grid models were chosen based on previous studies.

Erlebacher el. a.1. [54], Moin ['28], and Vreman [52] derived the model constants based

on DNS results. Erlebacher concluded that, C = 0.012, and Moin gave a range o['

va_lues of, 0.008 <_ C ___0.014. Moin a.lso provide a. range of values for Ct, where

0.0025 < ('l _< 0.009. From these data the coefficients were chosen to be C = 0.01'2,

and C_ = 0.00575. \:reman determined that C3 = 0.6.

Results fl'om LES simulations are obta.ined in terms of filtered conservation variables. These variables differ from the exact conservation variables by equation (.5.1).

NASA/TM--2001-210716 55
|tow('ver. since the form o[" th(' filter function. C, is unknown, we cannot r('cover lh( _

exact conservation variables. We must assuine t:hat the difference between tlw two

forms is small. For convenience, when reporting LES simulation results the overbar

( ) and tilde ( _ ) are dropped in the notation and must be asstmwd.

NASA/TM--2001-210716 56
CHAPTER 6

FLOW SOLVER

The numerical methods outlined in Chapter 3 are implemented in two computer

codes 1o solve the unsteady Nax'ier-Sl,okes equations. Separate three-dimensional and

two-dhnensional/axisynnnetric codes were created. The codes were written hi the

FORTRAN 77 computer language and should run on any computer l)latform with a

FORTRAN compiler.

\Vhile the purpose of the codes is to simulate compressible turbulent jets, they were

wrilt, en so that they can be easily adapted to a wide range of flows and geometries.

Large-Eddy simulations can be performed 1)3' solving the filtered unsteady NavierStokes equations with the sub-grid scale model, or direct, solution of the equations

without approximation can be made. By neglecting the viscous terms in the NavierStokes equations the Euler equations can also be solved.

6.1 Generalized Curvilinear Coordinates

To allow easy solution to a wide range of geometries the Navier-Stokes equations

are solved using generalized curvilinear coordinates [,55]. This enables the computational grids to be fitted to complex shapes and allows for grid stretching.

NASA/TM--2001-210716 57
6.1.1 Coordinate Transformation

l';asv and efficient implenlenlation of the numerical scheme is done by working in

a cOtnlmtalional domain. T]ds domain is rectangular and consisls of equally spaced

grid pohlts. The computational domain (_.//, q-) is mal)ped onto the physical domain

(x. !t, z) using the following transformation

,_ = {(.,.,::. :)

< = C(:_'.v- :) ((_.1)

Derivatives in carl esian coordinates are computed using the chain rule

O 0 0 0

O O 0 0

__=: __ + % ____.+ -

O 8 O 0

((_.2)

The terms {_., ¢_y, {:, 71_., %, .r/:, q'_., (u, and ¢': are the grid metrics associated with

the transformation. The grid metrics are computed numerically using the following

NASA/TM--2001-210716 58
method.

(,,, = J(.g,_:< - yc:,_)

_= = ,l (:r,_.q<-:r<.y,,)

q_. = .l(.q<:_ - Y_=c)

q_, = ,l(x_=c - x<=_)

71= = d(.rcy _ -,r_y<)

(_ = J(a'_9, _ -.r,ly_)

The Jacobian of the transformation. J. is defined as

(6.3)

1

d = (6.4)

The terms x_, .r,,, xc., ye, y_, y¢, z{, z,_, and z¢ are computed using fourth-order central

differences in the interior of the grid and one sided differences on the boundaries. The

calculation is easily done in computational space because the grid is evenly spaced in

the _, q. and _ directions.

6.1.2 Chain Rule Form of the Governing Equations

=0

The chain rule form of the governing equations is

0(2 o o _xo

0_T+_ _ (E- E,,)+,_._,_(E- E_,)+ _ _ (E- E,.)

0 , 0

( O(f_ Y,,)+ % (Y- F_,)+ (F- F,,) +

0 0

+ _° (c;- c;_.)+ (G- c_.) + ¢. (c;-(;,.)

NASA/TM--2001-210716 59
This f<)rmof the equations is a weak conservation form..Most ('Fl) solvers use a

strong conservat, ion law fbrm of the equa,t, ions developed by Vinokur [.5(;]. However.

Hixon el. al. [.57] recently showed that in practice, the chain rule form of the equations is more accurale than the strong conservation form when lh(' metric terms are

compuled numerically.

6.2 Time Stepping

The time step used to advance the solution is based on the relation

_SI = ('FL (_) ((i.6)

where .4 is the local waw' propagation speed and ('FL is the ('ore-ant Friedrichs

Lewev number [8], a constant which is based on stability considerations. The term

a,,_£-is a. characteristic time and is computed using the inviscid ('FL condition [58]. A

Xx ( lul I,,l lu, I 1 1

_xt_,_-_- .4 - \_x. + Xyy + _ + _ + --__._+ _ (6.r)

For time accuracy all points in the domain must be advanced at the same time step.

Therefore. --ktceL is computed at every grid point and the minimum value is used to

COml)ute the time step.

6.3

_Xt = Cf'L. AtCFL ....

Treatment of the Viscous Terms

(6.s)

The computation of the derivatives within the viscous fluxes ( equations (2.4),

(2.8), (2.13), & (2.15a)) is <tone differently for both numerical methods. The predictorcorrector schemes require a fairly complex treatment of the viscous terms. Bayliss et..

NASA/TM----2001-210716 60
a.l. [36] proposed the following melhod to maintain the oMer of accuracy of the .MacC.ormack type schemes. If the derivative within the viscous flux is being differenced

in the same direclion as the flux is differenced, then the clerivalive in the viscous flux

is computed using a one-sided difference in the direction opposite of direction the

flux is being differenced. If the derivative wilhin the viscous flux is being differenced

in a direction other than the direction lhe flux is differenced, then the derivative in

the viscous flux is computed using a central difference. (:onsequently, the derivatives

in the viscous terms are COlnputed twice (one-sided and central differenced) for each

step. The Runge-Kutta based schemes shnply use fourth-order central differences for

all viscous derivatives.

6.4 Boundary Conditions

Boundary conditions can be specified on any portion of any grid surface within the

computational grid. This allows for a large range of complex shapes to be modeled

within a single computational domain. Internal portions of the grid that are enclosed

by computational boundaries are specified a.s "holes". The conservation variables

within the hole regions are not updated during the solution process. Figure 6.1

illustrates the use of internal boundaries and holes to model a nozzle wall.

6.4.1 Inflow and Outflow Boundaries

The formulation of these boundary conditions is based on the local one-dimensional

propagation properties of the flow [46]. These properties are obtained from the eigenvalues of the Jacobian matrix in the quasi-linear one-dimensional Euler equations.

NASA/TM_2001-210716 61
Theflow isassumedto beone-dimensionalin th(,direction nor_nalto tho l)oun(lar\'.

The resulting characteristicvariablesare ui/_,, uil_, + a. and Ilill_ -- a. where It; iS lhe

unit normal vector to the boundary and a is the speed of sound. The variable uil_i

occurs twice in two-dimensional flows aim three times in three dimensional flows.

These eigenvalues determine the propagatioll direction of the characteristic variables

al tho boundary.

Subsonic Inflow

For a sut_sonic inflow boundary u,. < a. This results in four (three in twodimensions) positive characteristic variables and one negative variable. Therefore

four variables propagate from the boundary int.o the domain and one variat_le propagates out of the boundary from the domain. To correctly mimic this behavior in

the computation, so that the problem is ma.thematically well posed, four of the five

conservation variables are specified at the boundary aim one is (let.ermined from the

interior.

The flow from the inflow boundary is assumed to be normal to the boundary

so that the transverse velocity components are set. to zero. The total pressure, p0.

and total temperature, 7b are specified. The velocity normal to the boundary is

extrapolated fl'om the interior. For a constant _ surface (constant i) with its normal

NASA/TM--2001-210716 62
vector alignedin Ihex direction, tile conservation variables are are obtained as follow,s.

!

QII,.._.I - HTo u_ _l_ ) --r---7

^,--1 (,+l,jkl

Q2t,,j,_ / = Qli, jj.lu(,+,.j,_. /

(23[r j,k) _-_ 0

4(_,jj,) _--- 0

= ) ?

((_.!))

Supersonic Inflow

For a, supersonic inflow botmdary 'ui > a. This results in five (four in twodimensions) positive eigenvalues. Therefore all information propagates along the

characteristics from the boundary to the interior.

All five conservation variables are specified on the inflow boundary and are held

fixed. For a constant _ boundary the conservation variables are

Subsonic Outflow

Ql(,,_,k) = P_

Q2(,j,k) _-_ p,.,_ Ue¢

Q31,,:,k) = po..v_

Q5(,,:,k) = fl,_:( t.z

For a subsonic outflow boundary u i < a. This results in four (three in twodimensions) negative eigenvalues and one positive one. Here information fiom four

NASA/TM--2001-210716 63
of the characteristicspropagatesfrom tile interior to the boundary and informalion

fronl one characteristic propagates from the boundary to the interior.

The sta.1 ic pressure, p_, is specified on the otllflow boundary, l)ensity and all three

velocity COml)onents are extrapolated from the interior. Total energy is comtmted

using the specified static pressure and the extrapolated densilv and velocilies.

Qlf_max,j,k)

(22(*max,j,kl

Q'3i _ma r ,j.k)

Q4(,*rl.c*x,j,k )

QS(tmax,j,k)

Supersonic Outflow

= Ql(,,,,o___..,.k)

= (22( ......... --1,,,k

Q3(,m.x-l,3.k}

P_ Q22(*ma.r,3,k}

- +

")-1

2

+ Q._( ......... j,k) + Qq( ...... ,_,k_

2Q l i,,,,_.j.k)

(6.11)

For a supersonic outflow boundary ui > a. This results in all five (four in twodimensions) eigenvalues being negative. Information from all characteristics propagates from the interior out through the boundary.

All five conservation variables are extrapolated from the interior onto the boundary.

Q l(_maz,3,k )

Q2(_max,j,kl

Q3(*max,2,k)

Q 4( _m.:,.x,3,k )

Q5( tmax,3,k )

: Q10..ax_l,:.k)

: Q2(,ma__l,j,_)

: Q3(,maz_l,3,k)

Q4(tmax_l,3,k}

Qs[Imax_l,3,k) (6.12)

NASA/TM--2001-210716 64
Non-Reflecting Exit Zone

Wavesthat propagatefrom the interior of the domain t.othe outflow boundaries

ca.n,'effectoff the boundaryand back into the c()ml)utationa]domain. Thesereflected

wavesarenon-physicaland maycontaminateor distort the solution. Thesereflections

occur becausethe extrapolation of interior information to the boundary is donealong

grid linesand not in the direction of the propagatingwaves.

Expanding the computational domain so that the boundariesare fat' away from

the areaof interest reducet,his problem, but adds significantly to the computational

cost. Severalauthors haveattempted to elimhmte these reflections through the rise

of characteristic boundaries I59-61J. }towever, these efforts have met with limited

success for two- and three-dimensional applications.

:\ relatively simple and more effective method to reduce reflections is the us(" of

"exit zones" I18,61+621 . In this inethod a simple outflow boundary condition, such

as those outlined above are used. In the region adjacent to the outflow boundary.

a combination of grid stretching and solution filtering is used to damp the outgoing

and reflected waves. Rapid grid stretching increases ._kx and therefore increases the

associated dissipation in the truncation error of the numerical scheme (equations

(3.5). (3.8) and (3.10)). Solution filtering (section 3.4) provides additional dissipation

that damps a large range of wavelengths.

6.4.2 Solid Surfaces

Two options exist to model solid surfaces. Surfaces in viscous calculations are

specified with a no-slip wall boundary condition. Sin+faces in inviscid calculations

and planes of symmetry are specified using a slip wall boundary condition.

NASA/TM--2001-210716 65
No-Slip Wall

The no-slip wall boundary imposesno relative motion betweenIll(" fluid and the

solid surface. Entbrcing no-slipat the wall is doneby setting all velocity components

on tile boundary to zero. The density and total energy are determined using two

additional conslraints. The first constraint stems t'ronl boundary laver theory. The

pressure gradieni normal to the wall is set to zero. The second constraint assumos no

hea.t tramsfer through the wa.ll, and consequently tile temperature gradienl normal to

the wall is set to zero. For all *1constant wall (constant j). tile conservation varial)les

are written as

Q l(,,j,_,) _-- Qllt,_+I,t,)

Q2(,o,_: ) = 0

Q3(,,j,k,) = 0

Q4i,,,,k ) = 0

Qs(,,:,k) = Qs(i.+,,k) -

2 2

+ + Q4,,,.+,,k,

2Ql(,,)+l,k)

(6.13)

This formulation assumes that tile grid lines are normal to the wall.

Slip Wall

An inviscid wall surface and a symmetry plane share the same constraints and

hence the same boundary condition. The velocity vector of the flow along the surface

nmst be tangent to the surface. To accomplish this, tile velocity vector one point.

off the wall is decomposed into components that are normal and tangential to the

surfa.ce. The normal coml)onent is then removed from the velocity vector and this

vector is imposed on the boundary. The loss in kinetic energy due to the norlnal

NASA/TM--2001-210716 66
velocity is compensated[br by correcling ill(' total energyoll lhe boundary base(1oil

conservationof total entha.lpy. The density gradient a.tlhe wall is assumedlo 1)e

conslant. Vor an // constant wall (constant j), the conservation variables are writ len

a,s

Q I (,,a,k)

Q2(_,j,j,)

Q3( ,,j,kl

Q4(t,j,A)

('2 5{ i,a,k )

(6.14)

where the magnitude of the normal velocity. (':,_. is

7h, %

+ +

F{i,.s+z,k)

The metric terms 71_,, qv, and 'q: are evalua.ted on the l)ounda.ry.

71.=

v6b2.+ '1_+ tl_

(6._5)

6.4.3 Other

The following boundary conditions result from computational necessity, not. physical modeling.

Changes for Fourth-Order Schemes

The fourth-order schemes call for special treatment one grid point away from the

boundaries. Both fourth-order numerical schemes require two points on either side of

a. grid point to compute the deriva.tive a.t that point. One point from a boundary, the

scheme must be modified using a skewed difference stencil to prevent over-running

NASA/TM--2001-210716 67
lhe boundsof lhe storagearraysand to maintain fourth-order accuracy.The "'skewed

forwetrd'"and "skewedbackward" fourth-order differenceslencilsare

Ofo.r.r,:_ = .f_+3 - 6.f;+2 + llS./_+:2 (A.r )- IO.L -:U,- J (6.16a)

i) f bwd

(,).F

Pole Boundary

The collapsing of a grid surface to a line or

3f,+a + ]Of, - 18./)_a -4- 10.1)__ -.L-:3

12(Ax)

(6.161))

"'pole'" is often done to model cylindrical

objects in generalized curvilinear coordinates, t:igure 6.2 illustrates how a volume in

computational space is transforlned to a sector of a right circular cylinder. Vahms for

the variables on the collapsed surface (jmi,_ in figure 6.2) are obtained by averaging the

values over the collapsed index, one point off of the collapsed surface. For example.

on a conslant _l surface where ( lines are collapsed to a point, the values of the

variables at each _ point on the q surface are obtained by taking the mass average of

the va.riables over all _ at ?/+ 1. This bounda.ry condition is used for the axis in a

three-dimensional cylindrical calculation. The average flowfield values are determined

NASA/TM--2001-210716 68
])V

_ f_ _If_

-l

-I

k=l :=

-l

k,k.=lZ .....,k l,k']k\(k_........ [)k"_.]k]

k=l \

-1

_.=1 '1, \_.I,]

(6.17)

The conservation variables Oll the boundary are then

Axis of Symmetry

(6.is)

A boundary condition is needed for axisvmmetric analyses on the bounding grid

line where r = 0, the axis of symmetry. Since the flow is symmetric about this axis

no flow may cross the grid line. Therefore, the velocity vector must be tangent to the

boundary. This restriction is the same as the restriction for the slip wall boundary

and hence, the same method is used to impose the constraints (equation (6,]4)).

NASA/TM--2001-210716 69
Overlap

The grid lOl)ologv commonly ref(,rre(l 1o a,san "'()"-grid ix fr(,(luenlly used to

model round objects an(1domains, such as cylinders, airlbils and pipes. "l'he domain is formed b\ wrapping the grid around upon itself so lhat the n_hiimum and

maximum boundariesof a commonindex are coincident(figure 6.3). The coincident

(Oml)ula,tional l_oundariesare in the interior of lhe l)hysical domain and nltlSt l)e

treat,e(l a.sif they werealso on the interior of the COlnl)utalionaldomain sothat they

(1ouot createa non-physicalbarrier on the interior of the domain. In olher words.

the governinge(luationsshould be solvedon theseboundariesas if they are in lhe

inlerior ()ifthe computational domain and information musl be readily 1)assedacross

the boundaries.This is doneby overlapl)ing the grid so that at every point wherea

skeweddifferenceis performed the result,can be replacedwith a centra, l differenced

value from a coincident poinl in the overlapped region. Table 6.1 lists the skewed

differenced grid points and the coincident central differenced points necessary for the

overlap boundary condition.

boundary point interior point,

(skewed difference) (central difference)

kma z -- 1 3

Table 6.1: Coincident, grid points in overlap region

NASA/TM--2001-210716 70
6.5 Computer Resource Requirements

The menmrv all(:] ('PI" time required are reported for both tile 2I) aud 3D codes.

The resource requirements reported are for Silicon Graphics workstations, which were

nsed for all the large scale calcttlations in this documenl.

Memory requirements for the two-dimensional/axisymmelric and t hree-dimellsional

codes are reporled in table 6.2. Both single and double 1)recisitnl re(luirelllenl[s are

given.

bytes/grid pt.

code single precision double precision

2D/A xi. 2,44 472

3D 360 700

Table 6.2: Memory required

Table 6.3 contains CPIT times for all possible combinations of code (2D, axisymmetric, or 3D), equations (Euler, Navier-Stokes, or LES) and numerical schelne. The

execution times are presented as time per iteration per grid point in seconds. The

CPU times were obta.ined using a. Silicon Graphics Power Indigo 2 workstalion with

a MIPS R8000 processor running at 75MHz. The code was compiled using double

precision data and the compiler options were chosen for maximum performance.

6.6 Parallel hnplementation

Parallel processing is an efficient and inexpensive method t,o speed t,he execution

of computer programs [63-65]. [;NIX workstations with multiple central processing

NASA/TM--2001-210716 71
code equations scheme time/it.er./grid pl. (sec)

2D

Axi.

3D

F, Ulel'

Eu ler

Navier-Slokes

Navier-Stokes

LES

LES

Gottlieb-Turkel

l:{un ge- K u t t a

(k)t tlieb-Turkel

Runge-Kutta

(;ot tlieb-Turkel

R unge- Kut ta

9.092" 10 -'_

t .65:1" 10 -_

1.218-10 -_

2.122.10 -'_

1.655.10 -a

3.196.10 -_

Etller

Euler

Navier-Stokes

Navier-St, okes

LES

LES

(',otllieb-Turkel

11unge- K utta

Got t lieb-Turkel

1Runge-Kutta

C,ott lieb-Turkel

Runge-Kutta

9.415.10 -(_

1.7:32' 10 -a

1.498.10 -_

2.808" 10- _

1.982" 10 -5

:/.945-10 -a

Euler

Euler

Navier-Stokes

Navier-Stokes

LES

LES

Got t lieb-Turkel

Runge-Kutta

Gottlieb-Turkel

Runge-Kutt a

(;ottlieb-Turkel

Runge-Kutta

6.S98.10 -r'

1.069.10 -4

1.182.10 -4

1.555.10 -4

1.:t.94.10 -4

2.074.10 .4

Table 6.3: CPIT time required

units (CPU's) are widely available and cost a fraction of the price of super-computers.

('odes written to take advantage of the parallel architecture of these machines can

execute as fast or faster than on a comparable super-computer. Parallel processing

achieves this performance increase by dividing the work normally done by a single

central processing unit amongst several CPU's.

There are two primary methods used for the parallel processing of CFD solvers.

The first is distributed memory processing. In this method each CPU uses a separate

bank of memory for its calculations. Communication of data between CPU's is done

explicitly through user implemented interfaces. Often the CPU's do not reside in the

N ASA/TM--2001-210716 72
salne conlputer, but operate as separate workstations connected through a u(qwork. A

t,ypica.1 distxibut, ed memory (',FD solver divides the COlnlmt ational domain into several

sub-domains which are each solved on separate processors. The illterfaces between

sub-domains are treated as boundary conditions and are updated by l)assing updated

solution values to adjacent sub-domains. This method requires limited COllanlUlliCation between processors, and is therefore suitable for networked workstat ions, because

data is only passed a.t the end of each iteration or set. of iterations. It works very well

for low order steady-state flow calculalions. High order unsteady ca.lculations are ill

suited for this method because it would be necessary to pass large amounts of data

between interfa.ces a.t every itera.tion to maintain both temporal and spatial accuracy.

The second method, shared lnemory processing, was used in this stlutv. In this

method, a set of' ('P["s share access to a common bank of memory. Communication

between processors is implicit and no special interfaces are required. Shared memory

systems consist of multiple processors housed in the salne computer with a high

speed high volume bus connecting the processors to the memory ba.nk. Because

the processors are closely coupled through the shared Inemory, the division of labor

between processors can also be closely coupled. This division of labor is done at the

loop level. The work performed in a. FORTRAN "DO" loop is divided between (!PU's

with each processor working on a fraction of the total loop. For example, a "'DO" loop

which is indexed from 1 to 100 could be split between 2 processors with one processor

operating from 1 to 50 and the second operating fi'om 51 to 100. hnplementing

this method is done through the insertion of compiler directives into the code. The

directives indicate which loop is to be parallelized and how the variables are shared

in memory. Care must be taken to identify data dependencies so that variables whose

NASA/TM---2001-210716 73
valuedependson lhe current loop inde×are nol overwrilIen I)\"processorsoperalillg

on anot,her portion of the loop. t')xplicit,(',FD methods are well suiled for this type

of paralle]izalion because the data dependencies are limited and easily overcome.

The performance of a parallel code is measured by its speedul_ and efficiency.

Speedup is simply lhe ratio of lhe time required for single processor computation to

lhe {ilne required for a multi-processor COml)uta,tion.

I'[

s(.,,.) = --

rllR

where T_ is lhe time required for an _ processor computation. The eIficiencv of lhe

job is the ratio of the speedup to the number of processors.

S(n) T1 E(,,) - {0.e0)

Ideally a the speedup of a parallel code would vary linearly with the number of

processors. Two primary: factors limit the speedup; processor idle time and communicaliou time between processors. Processor idle time is the time the code spends

outside t,he parallelized loops. During this time only one processor is active and the

potential work of the other processors is wasted. Communication time is the delay during which data is passed between processors and lhe shared memory. I1 is a

property of the comt)uter and varies greatly I)etween machine types.

The performance of the 2D/axisymlnetric {ode was measured on a Silicon Graphics

36 processor Power Challenge computer. The cod<' was run on a representative jet

calculation using a grid containing 38,829 grid points. A calculation consisting of

1.000 iterations using the Runge-Nutta scheme was repeated 10 times using an even

number of processors from 2 to 20. Results were compared to a baseline calculation

performed on one processor. This procedure was repeated four times. Speedup and

NASA/TM---2001-210716 74
efficiencyare])resentedin figures6.4and 6.5. Thevariation in resultsbetweeneachset

of data,is due to the architectureof the computer thai was used. The data (memory)

bus for this machine was insuflicienllv sized to transfer dala for all 36 processors. The

communication delay time varies with lhe number of jobs being ruIl on the svsteln

and the number of processors R)r a given parallel job. Each set of parallel runs were

lnade with the colnl)uter under different loads and hence the tilne required R)r the

runs differed. Regardless. the data do show that a l]nea.r speedup is seen up to about

16 processors. Because the computer's performance varied with its load, the parallel

efficiency cannot be definitively determined. In fa.ct figure 6.5 shows some etficiencies

exceeding one. This is mosl, likely due to the load on the machine being reduced after

the one processor job was run. Consequently. the 1)aseline one processor run is nol

appropriate for computing the efficiency of the current multi-processor run.

6.7 Validation

To insure accurate implementation of the numerical schemes, the codes were valida.ted using three simple test cases in both two and three dimensions. The code's

results were compared to exact solutions I.o insure that they produced acceptable

solutions.

The standard Runge-Nutta scheme (section 3.2.1) was unstable with or without

solution filtering and did not provide a converged solution for any of the two- or threedimensional calculations. The low dispersion Runge-I,/utta scheme (section 3.2.2)

did provide stable solutions for all calculations attempted. All Runge-I,/utta results

reported hereafter were obtained with the low dispersion scheme (equation (3.9)).

NASA/TM--2001-210716 75
6.7.1 Laminar Flat Plate

La,ninar flow over a flal plate was used to check the accuracy of a viscous solul ion.

Skin frictioll along tile plate and tile self-similar velocily profile are con_pared to the

exact solution of Blasius. as givell by Schlicting [661. The skin friction is con_puted

using

and the plate Reynolds nunlber is

(6.2J)

l :%a"

Re,. - (6.22)

The velocity profiles normal to the plale surface can be made self-similar 1)y introducing a dimensionless coordinate.

The Reynolds number based on plate length is 10,000 and the freestream .Mach

number is 0.2.

Two-dimensional solutions were obtained with both numerical schemes. The sixthorder rifler was used for both calculations. The grid that was used is shown in figure

6.7. The grid dimensions were 157 points in the streamwise direction and 117 points

normal to the plate surface. Subsonic inflow and subsonic outflow boundaries were

specified on the left and right boundaries respectively. A no-slip wa.]l wa_sspecified on

the lower boundary and an extrapolation condition (the sarne as supersonic outflow)

was specified on the upper boundary.

Skin friction results for both the [{unge-B[utta and Gottlieb-Turkel scheme are

shown in figure 6.8. The Runge-Nutta scheme slightly over predicts the skin friction

NAS A/TM--2001-210716 76
relative to tile Blasiussolulion. The Gottlieb-Turkel schemeunderpredicts tile skin

frier,ion and is in greatererror than the l{unge-I(utl.a scheme.

The self-similar velocity profilesare shownin figure 6.9. Severalstreamwisestations. /_c_- = 2000, 4000, 6000. and 8000 are plotted i.o verify the code reproduces

the self-similarity. The Rulage-I(ulta results are very close to the Blasius solulion.

The (lottlieb-Turkel profiles show a lhicker boundary layer thickness with a smaller

velocity gradient at the wall. This is consistenl with the skin friction data..

To validate the lhree-dinlensional code. a 3t) grid was created thai consists of

fifteen evenly spaced 2D grid planes. The 3D code was run using both numerical

schemes and the results were compared to those of their 2D counterparts. Skin

fi'iction and velocity profiles a.l /?(x = 4000 are compared fbr both schemes {figures

6.10 ,_,"6. [ 1). The agreement between between the two- and t hree-dimensional codes

is excellent for both schemes.

6.7.2 Supersonic Wedge

Supersonic flow over a two-dimensional wedge was used to test the ability of the

code to predict inviscid flows and shock waves. A 15 degree wedge in a Ma.ch 2

freestream flow was simulated (figure 6.12). Pressure coefficient on the wedge surface

was compared to the exact, solution.

p- p,_:,

(6.24)

The grid dimensions used were 121 points in the axial direction and 81 points in the

vertical direction. The grid is shown in figure 6.13. Uniform supersonic flow was set

at the inflow boundary. Supersonic outflow conditions (extrapolation) were set on

the right, and upper boundaries. A slip wall was specified on the lower boundary.

NASA/TM--2001-210716 77
Two-dinwnsionalresultsare showll in tigure 6.14. ()retail the resultsare in good

agreement wit.h the exact, solution. But, the shock location is smeared over several

grid points and some oscillation in the solution exist both upstreanl and downstream

of the shock wave. A comparison of numerical schemes using the sixth-order filter

shows very little difference. The presence of the shock wave reduces both schelnes to

tirst-order accuracy (the location of the shock is directly related to the grid spacing)

nullifvhlg any accuracy advantage of the 1Runge-Kutta scheme. The affect of fourth,

sixth, and eighth order filters on the shock location and oscillations was examined

wit h t he I{unge-i,_utta scheme (figure 6.14(b)). The lower order fitters reduced t he

oscillation without significant affect on shock local, ion or strength.

Results of the 2D and aD codes are shown in figure 6.15. For both numerical

schemes the 2I) and 3D pressure distributions compare very well.

6.7.3 Supersonic Cone

To test the axisvmmetric terms in the two-dimensional/axisymmetric code and

provide a true three-dimensional flowfield for the three-dimensional code, flow over

a supersonic cone was validated. This test case uses the same flow conditions and

turning angle as the wedge case. But, the body is treated a.s a body of revolution

rather than a planar object. For the axisymmetric simulation, the grid for the wedge

case was used (figure 6.13) and the axisymmetric source terms in the flow solver were

activated.

The three-dimensional grid represents a 15 degree section of the cone and its

ttowfield. It was generated by' rotating the 2D/axisymmetric grid about the centerline.

A 2D/axisymmetric grid plane is located every one degree in the physical domain.

NASA/TM--2001-210716 78
The exact solution for tile pressurecoefficien_on t.hecone surfacewas obtained

from the cha.rt,s in NA(:A [/eport 1135,"li'_qua.tions,Tables, and ('harls for ('ompressibleFlow" [67].

Resultsfor both schemesin 2D avevery good (figure 6.17). [_oth s('he_Hesstneal"

the shock location over several grid poinls, l!nlike the wedge calculation lhe shock

oscillations are weaker and do not over-shoot the pressure behind the shock. The 51I)

results from the (;ott.lieb-Turkel scheme (figure 6.18(a.)) differ fl'om the 2D resull.s

due to a more l)l'onounced pressure oscillation about the shock. The 3D results ft'om

the 1Runge-Kutt.a scheme (figure (i.[8(b)) compare well t.o its 2D counlerpart.

NASA/TM--2001-210716 79
i

I

i

I; i

i

(a) computational grid

(b) internal object in grid (c) grid point specification (filled = aztive, open = hole, both = boundary)

Figure 6.1: Modeling an internal object using hole points

NASA/TM--2001-210716 80
• /

Imin

(a) computational space (b) physical space

Figure 6.2: Pole boundary condition

NASA/TM--2001-210716 81
Jma.t

/

/

Jmm

(a) ('omputational space

/ tma'T

imax

i

i . ,'*

i,,,i,, _ _ "Jm,=

(b) physical space without overlap (c) physical space with overlap

Figure 6.3: Overlap boundary condition

NASA/TM--2001-210716 82
22

20

18

t6

14

12

10

8

6

4

2

0

_ F--",_ea,1 ! i i i/ t

_--+I-+-+e,l ......+.....t--, ........ +............. +++_+_ t

F J'-+-se':i : +.':.......!.......... .i,+._ ........ +...4

F-I --_--Set:l--_--++ "++ i A 3 "I

l--'--+e"l +....L............. L+_++V_-___,2

............ +........... ¢.............. t........... -_........... +--..._-_',_:- ...........

...... t--_---:

+

+i +,+

+

......._--t _ --T-- _

............. i.............. +.......... +_'+ ...........

.._.._+...+..._ +.........+,_ .... + .... , +

+ + ,+ + + :

i + t + t ..... !............ +............;

_.._: + i + + i ! + +

...i,,,i.,.i,..t,,,i,,,i,,.i...i.,,t... ,,,

2 4 6 8 10 12 14 16 18 20 22

Processors

Figure 6.:1: Parallel processing speedup

1.2

0+8

• i

i

0.4 _ _ Set 1 |

I

• .--B_ Set 2 I

• _Set 3 I

0.2

. ._-- Set 4 I

..[ .....

0 iii

0 2 4 6

......_ T"Tj i ...... ''"i''l'''!';! i i

• ' " I " ' '

I

I

I

10 12 14 16 18 20 22

Processors

Figure 6.5: Parallel processing efficiency

NASA/TM--2001-210716 83
X

P

Figure 6.6: _D laminar flal plate

Figure 6.7: 2D grid for laminar fiat plate calculations

NASA/TM--2001-210716 84
o0.05

0.04

0.03

0.02

0,01

J ' ' ' ! " " ' I

_'-'--'_'t _ Blasius

-- - " G°ttlieb-Turkel t

i i

_ i i i i

, ,, i, .. i ,.. i,,, i ..,

2000 4000 6000 8000 10000

Re x

(a) Gottlieb-Turkel scheme

0.05

oo, ........

___ I

0.02

! 1 _ -

0 2000 4000 6000 8000 10000

Re X

(b) Runge-Kutta scheme

Figure 6.8: Skin friction coefficient for laminar fiat plate

NASA/TM--2001-210716 85
1.2

0.8 --

0.8F----L--_ _---"_ ', =2ooo1_..--

i //

0,2 ,', i i _,&h_. _ i __.

0 1 2 3 4 5 6 7

(a) (;ottlieb-Turkel sche]ne

1.2

0.8

0.6

0.4

0.2

.... _........ T.... ! .... ! .... T.... ' ....

4 i if"-- i

...... F .............................................

i_ _" i E -- Btasius I

I J7 i i I ..... Re =20001

i • ! ,1--- _=,,ooo1_

/ i E i l..... .;:,ooo/

i_ 1 F/....."e'='°°°/

V...A .... i_..h_.i.,. /i.... i ........ .

0 1 2 3 4 5 6 8

(b) Runge-Nutta scheme

Figure 6.9: Velocity profiles for laminar flat plate

NASA/TM--2001-210716 86
0.05

0.04

0.03

0.02

0.01

I

\

_ Blasiusl

" i I

.,,t..,I.,,I,..I...

0 2000 _ 6000 8000 10000

Re

x

(a) skin friction coefficient

1.2

0.8

3 o,6

0.4

0,2

0

0

.... I .... ii ........................ i 1i ii ::1 ]!I

....... i .... , .... i ........ , .... , ....

1 2 3 4 6 7

(b) velocity profiles at, Rex = 4000

Figure 6.10: Comparison of 2D and 3D Got.tlieb-Turkel schemes for laminar fiat plate

NASA/TM_2001-210716 87
0.05

0.04

0.03

0.02

0.01

! ...... r i

......................................... _. .... [ -- Blasius .....

• . . i . , , J i

2000 4000 6000 8000 10000

Re X

(a) skin friction coefficient

s

1.2

0.8

0.6

0.4

0.2

.... ! .... ! .... I .... ; .... t ....

i i i _ i

i iy

........ i .... i ........ ]

.... I ....

J

F

- t

_B,a,iusli

..... 2D

-- - - 3D

1 2 3 4 5 6 7 8

(b) velocity profiles at, Re_: = 4000

Figure 6.11: Comparison of 2D and 3D Runge-J<utta schemes for laminar flat plate

NASA/TM--2001-210716 88
Y j"

/"

//

.,.¢"

M _ .//"

P_-. ///J

N\,

\.

\'x

"\'\.

N.

\

x

%,

x

x',N

Figure 6.12: Mach 2 flow over a 1:5 degree wedge

NASA/TM--2001-210716 89
Figure 6.13: 2D grid for wedgeand conecalculations

NASAJTM--2001-210716 90
0.5

0.4

0.3

0.2

0.1

-0.1

....... I .... ! .... I ....

i

I

[

i

................. _L....................

.... i .... I ........

q U-

] i ' i

i

.... i .... i .... i .... , .... , .... J ....

-0.5 0 0.5 1 1.5 2 2.5 3

i

i ....

(a) Comparison of nmnerical schemes with 6th order filter

0.5

0.4

0.3

¢j_" 0.2

0.1

!t:'ti_..[ _ I !

i PI "'1*1] - i

t,

" I--- Exact I

[ ] ._ I I ..... 4thOrderFilterl

_ _ _ - - 6th Order Filter I

i i I ..... ath Order Filterl

i ....

-0.5

-Ol !.. i

-1 0 0.5

i ....

i

1 1.5 2 2.5 3

(b) Comparison of filters with l_unge-I(utta scheine

Figure 6.14: Pressure distribution on wedge surface

NASA/TM--2001-210716 91
0.5

0.4

0.3

0.2

0.1

-0.1

-I

"_L !

I__ I L___Exao,I

". : _ ..... 2D II-

! i !

• / J i i

I ] ! _ i

................... f................... I_................................... i.................. ';.................... _.................... _................

.... i .... I .... I .... I • • • , | i , , , I .... I ....

-0.5 0 0.5 1 1.5 2 2.5 3

(a) Gott,lieb-Turkel scheme

0.5

.... ! ....

J

0.4

0.3

0.2

0.1

0 -'

iO. 1 .... I ....

-1 -O.5 0

f .... ! ........ I .... I ........

I ] I --J_,,_i I !

......... I i i

.... !......t---

[

0.5

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

1.5 2 2,5

(b) Runge-Kut, ta scheme

Figure 6.15: Comparison of 2D and 3D codes fbr wedge calculation

NASA/TM--2001-210716 92
M

P=

......... -1,

Figure 6.16: Mach 2 flow over a 15 degree cone

¢,¢t

0.25

Exact |

.... Runge-Kutta I

- G°ltlieb'Turkel I i

Figure 6.17: Pressure distribution on cone surface

NASA/TM--2001-210716 93
o.

C.)

0.25

0.2

0.15

0.1

0.05

-0.05

-1

J

' [ t :: ] .

|:/ i II--aacll i

Y+'-_-- .............. H ..... 2 D I.- ......... '

I// il --_ I ;

................... l...................

-- i

J

-0.5

]

F

i

I .... I .... i .... I .... I .... i , i i I

0 0.5 1 1.5 2 2.5

(a) Got, tlieb-Turkel scheme

0.25

0.2

0.15

L._: 0.1

0.05

-0.05

-1

hEXaC'I I

--_ - .... 2D 14"_---

-- - - 3D

- i i i

-0,5 0 0.5 1 1.5 2 2.5 3

(b) Runge-Kut, ta scheme

Figure 6.18: Comparisol] of 2D and 3D codes for cone calculation

NASA/TM--2001-210716 94
CHAPTER 7

COMPUTATION OF A NOZZLE FLOWFIELD

Tile flow solver described in Ill(' previous chapter is now a,pplied (o a round supersonic jet. The computat, ional results are compared to high quality detailed data. The

effects of grid resolution and several numerical modeling parameters oll the solution

are studied. Finally. a thorough examination of the flowfield is presented.

7.1 Description of the Nozzle

The nozzle geometry and data of Panda and $easholtz [68-71] were selected for

study. This data set was obtained using the non-intrusive Rayleigh scattering technique. It is based on the measurement of laser light scattered by the air molecules.

This technique eliminates errors due to probe interference in hot, wire and pitot probe

measurements and biasing errors due to seed particles in Laser Doppler Velocimetry

(LDV) and Particle Image \'_locimetry (PIV) measurements. The data consists of

time averaged centerline and radial velocity profiles.

The nozzle is a one inch diameter convergent-divergent round nozzle with a designed exit Mach number of 1.4. The nozzle exhausts into quiescent air. The Reynolds

NASA/TM---2001-210716 95
nunlber basedon nozzlediameler is 1.'2million. Tile nozzlewasoperatedat its lllillimum shock condition which wa,s slightly less t,ha,n ils design Mach munber. The

operating conditions are given in table 7.1.

quantity symbol value milts

ra.tio of specific heals _, 1.4

nozzle plenum pressure P0) 6524.9 _).__z ft 2

nozzle plenum temperature T% 540.0 R

nozzle exit Mach nttnlber :_1: 1.::/95

nozzle exit diameter l)j 0.0833 .fl

jet velocity l,_i 1348.q L

ambient pressure p,_, 2064.8 /T

aml)ieut temperature 71_<. 534.6 R

Reynolds number /?el 1.2- 10 G

Table 7.1" Nozzle operating conditions

7.2 Flowfield Statistics

Analysis of RANS solutions is both simplified and limited by the Reynolds averaging process. The RANS ('odes readily provide a time averaged solution for analysis.

However, these solutions lack any unst, ea.dy flowfield information. The unsteady solutions obtained here contain much more information. But the amount of information

is overwhelming and must be processed and simplified before it can be successfldly

analvzed.

7.2.1 Time Averaging

The instantaneous solutions produced by the analysis code sometimes bear little

resemblance to the time average of the flowfield. To compare the present results

NASA/TM--2001-210716 96
with experimental data, enlpirical correlations and exact solutions a t,ilne average is

computed as

1 ft+At ,f(x. t)_o (r.l) f(.,.)=

and the ]nstai_taneous value of the function can be written as a, sum of the time

averaged and fluctuating quantities

.t = /+ f" (7.2)

In this notation the traditional over-bar (-) and prime {' ) have been replaced with a

double bar (=) and double prime (") to avoid confusion with the filtered quantit, ies

in the LES equations.

For the jet calculations a time average of density, pressure, and the three velocity

components were kept. In addition, the time average of the squares of density and

velocity components are kept to compute turbulent statistics.

7.2.2 Turbulent Statistics

Turbulent statistics can readily be computed from the time averaged flow variables

being stored during the calculation.

f,,---7= f-_7_ _ (7.3)

By subsl_il;uting the velocity components into equation (7.3) we obtain the normal

,,2

components of the Reynolds stress u

simply

-2 It -2

, v , and The root mean square value is

f,.,_, = _ (7.4)

NASA/TM--2001-210716 97
The root meansquarevelocitiesare typically expressedas lul'bulent inlensities.

-_ 7l r m ,s V/f;- 2

Urms _.2

_' - - (7.(i)

t,_i t

U'rn_s _,2

Turbulent kinetic energy is defined as

k=_ 1 (_ + _,,2 + u,,,2 ) (7.s)

For the LES calculations the turbulent kinetic energy is comprised of resolved and

sub-grid scale comi)onenls.

k= _ + : + ) +k,o, (7..q)

where

Ckk

ks_s - - (7.10)

2t_

\Ve nondimensionalize k using the jet velocity.

7.2.3

k

k'- 1 r2 (7.11)

Two Point Correlations

Correlations of velocity signals at two points within the flowfield can lend insight

into the structure of the turbulence [72] & [2]. Chu [73] experimentally applied two

point correlations to jet flows. He obtained a turbulent length scale (eddy size) and

convection velocity from the correlated data. Scott [74] presented preliminary two

NASA/TM--2001-210716 98
I)oint correlation data, obtained with (!FD. for Chu's jet. While the comparison of

tile ('FD and exl)erimental correlations was inconclusive, due to insufticient sinmlation time, this work indicaled (!FD obtained correlations were were possible. Chu's

methodology was followed in this study.

The location where the turbulent information is desired is sl)ecilied by the vector

Oi, whose origin is at, the center of the nozzle exit (figure 7.1). A pair of points evenly

spaced on either side of this location is specified by a separaliol_ c'i. Inslantaneous

velocity data at these three points are saved for the correlations.

The general form of t he two point space-time correlation coefficient is

17;'(0 i -- C'i/_,t)[,T;'(C)i + _'i/O_,t "-_ T)

r) = (7.12)

t 6 (c0 ,t)

where [:e is the component of velocity that makes an angle 0 with the jet axis and r is

a delay time. The correlation coetticient is a measure of similarity of the two signals.

A value of _ near +1 indicates that the two measurements are highly correlated.

The location about which the two point correlations were taken was located six

jet diameters downstream of the nozzle exit on the jet lip line. Eight separations in

the axial direction of up to one jet diameter were used. Table 7.2 summarizes the

location and separation vectors for the two point correlations.

Turbulent Length Scale

Two point space correlations, where the delay time is zero (r = 0), can be used to

determine a turbulent length scale. Several sets of correlation coefficients at different

separations are computed. The separation distance over which the velocities are

highly correlated is used as an estimate of the length scale g.

./f [ = R(' 'g,_,0)dv,_ (7.1:t)

NASA/TM--2001-210716 99
vector coordinates

location oi [6Dj, D i2.0] T

I st sepa.ra.tion (,] [1)_/S, 0, 0] T

2nd separation _'_ [Dj/7, O.O] r

sep ,'a, io,, o,o]

4th separal ion _.,_ [Dj/7). O. oJr

5th separation ¢'_ [D)/4,0,0] r

6th set)aration t/,__ [Dj/3,0, O]"r

7th separation ('; [l)j/2.0. O]T

Sth sel)aral ion _,'2 [Di. O, o]T

Table 7.2: Two point correlation location information

Convection Velocity

Two point space-time correlations can be used to determine the convection velociI 5'

of a t,urbulent eddy. For a given separation, q,i, the correlation coefficient is computed

over a range of delay times. ]'he value of 7- where the correlation coefficient is a

maximum indicates the time necessary for a disturbance at the upstream point (O, -

t/,_/2) to reach the downstream point ((5_ + _b4/2). This time and the separation

distance are then used to compute the convection velocity.

7.3

[ 'il

Axisymlnetric Solutions

Assuming tha! the flow is symmetric about the jet axis reduces the problem to two

dimensions. Using this assumption would significantly decrease the computational

expense of a solution and allow for a greater number of parametric investigations to

be done. However. turbulence is inherently three-dimensional and it is not expected

NASA/TM--2001-210716 100
that a solution obtained with the axisvmmetric versionof the CFD solver will yield

physica.llyrealistic result,s. Here,the a.xisynuuetricanalysesare usedto investigate

munerical issues,grid resolution, numerical scheme,and boundary conditions, and

the knowledge gained is applied to the three-dimeusional simulations.

7.3.1 Grid Generation

The computational grid specifies the nozzle geometry and computational donlain

for the calculaliolls, l+nlike previous LES analyses of jets, the internal nozzle contour

and nozzle lip were modeled. Tile internal nozzle boundary laver and vortex shedding

fi'om the nozzle lip may affect the growth and stability of the shear layer. It was felt,

that it was important to include these effects in 1he simulation. Also. by including

the internal nozzle the inflow boundary is moved away from the region of interest.

the jet shear layer, reducing the influence of tile artificial boundary condilion.

The computational grids used for the analyses were generated using the commercial software package Gridgen [75]. The grid points were clustered near the nozzle

walls using hyperbolic tangent stretching to resolve the boundary layers [76]. Tile

grid was also clustered near tile nozzle exit plane to capture the unsteady vortex

shedding fi'om the nozzle lip and the initial formation of the shear layer. The cell

aspect ratio at the upper and lower corners of the nozzle lip was set to one. (;rid

point clustering in the expected region of the shear layer downstream of the nozzle

lip was also implemented. The computational domain extends 20 nozzle diameters

downstream of the nozzle exit, in the axial direction and 10 nozzle diameters from the

,jet centerline in the radial direction. A representative grid with 301 point in the axial

direction and 129 points in the radial direction is shown in figure 7.2.

NASA/TM--2001-210716 101
7.3.2 Specification of Boundary Conditions

All lhe axisvnnnelricjel computaiions used lhe same boundary con(tilions (figure

7.3). Subsonic inflow using the jel pl('num conditions was specified at the inflow of

tile internal nozzle section. The no-slip wall condition was used for the internal and

external nozzle surfaces. An axis of symmetry/slip wall condilion was used on the,jet

centerline.

Simula, ting the quiescent conditions of the surroundings can cause mmlerica] difficulties. The 1)oundarv conditions are forinulated a,ssunUng a known flow direclion

(inflow or oullqow). In the still air small disturl)ances may cause portions of a boundarv to have an inflow and other portions to have an outflow. This leads 1o over or

under specification of the boundary and nunlerical errors result. To overcome lhis

problenl a small freestream flow is imparted to ensure a coherent flow direction and

properly" specified boundaries. In this case a Mach 0.05 freestreana was used. Subsonic

inflow was specified on the upstream external boundary using the fleestream total

conditions. Conditions on the upper boundary were extrapolated from the fl'eestream.

The subsonic outflow boundary condition imposes a constant pressure over the

entire boundary. In reality' the pressure on the outflow boundary varies in both time

and space, To accommodate this pressure variation while maintaining the correct

fl'eestream flow a combination of subsonic outflow and extrapolation conditions were

used. The subsonic outflow condition was specified on the upper portion of the downstream boundary to maintain the correct pressure level in the freestream. Conditions

on the rest of the downstream boundary; were extrapolated from the interior to accommodate both temporal and spatial pressure variations near the jet centerline.

NASA/TM--2001-210716 102
7.3.3 Effect of Grid Resolution

The effecl of grid resolulion on 1he solution was invesligated. Three grids were

used. where lhe grid spacing in each direcl ion was halved from the previous grid. All

solutions were obtained with the Runge-Kutta schenle using an eighth-order filter.

The sub-grid scale model was not used. Because the grid spacing is not uuiform.

it is ditTicull to quantify. Here, the spacing was represented by a nondimensional

computational grid size.

:xa,"= (7.15)

As the grid resolution was increa.sed, t.he resolution of the unst, eady flowfield iragrid dimensions A.r* u,,,,_, _',,_,_ kT,,_,.,.

coarse 151 x 65 1.0206.10 -_ 0.20006 0.15480 0.047999

medium 301 x 129 5.1031 • t0 -:_ 0.26266 0.20305 0.096392

fine 601 x 257 2.5516- 10 -:_ 0.32909 0.22788 0.12178

Table 7.3: Effect of grid resolution on axisymmetric solution

proved. The change in entropy of the flowfield can be used to visualize the turbulent

structures in the flowfield. The viscous mixing of the jet and ambient air increases the

entropy of the flow. The vortices alter the shape of the mixing region and this altered

shape can be seen in the gradients of entropy in the lnixing layer. Figure 7.4 shows

contours of entropy for the three grids used. The plots show a dramatic increase in the

resolution of vortical structures with grid resolution. A more quantitative measure,

the max]muna of the turbulent statistics in the flowfield, are shown in figure 7..5. The

turbulent intensities and kinetic energy increase as the grid spacing decreases (as the

NASA/TM_2001-210716 103
grid dimensions .._N;r/h

coarse 15l x 65 2081.5

medium 301 x 129 1040.8

fine 601 x 257 520.40

Table 7.4: Comparisonof grid spacingto the lqolmogorovscale

grid resolutionincreases).This result is expectedasthe intensities and kinetic energy

shouhl increaseas the grid is refinedand tile schemeincreasingly cal)tures more of

the smaller energy containing eddies. With even flarther resolution one could expecl

the energy to decrease as the inertial subrange and Kohnogorov scales are reached,

where the energy is dissipated.

TO compare the grid spacing to the Kolmogorov scale, an average cell size was

computed as follows

where L v,,i8 and Hj,.,,_ are the overall length and height of the grid. The Kolmogorov

scale was estimated using equation (1.1) and assuming the integral length scale is

l

a pproxinaately _D 3. Table 7.,1 shows that even for the finest grid the cell size is over

500 times larger than the Kolmogorov scale.

7.3.4 Comparison of Numerical Schemes

A comparison of the two numerical schemes was performed on the mediuln sized

grid (30t x 129). The solutions were run for two characteristic acoustic times (the

time for an acoustic wave to pass through the solution domain) to establish proper

initial conditions for obtaining turbulent statistics. The solutions were then run for an

NASA/TM--2001-210716 104
additional two characteristictimes andthe flowfie[dwasaveraged.Both schemeswere

run usingtheir tn_tximunlst,abletime step. Figure 7.6showslhe entropy contours for

the two schemes.Tile Runge-Kutta scheme exhibits more vortical structure and the

initial vortex roll-uI) of the shear laver occurs earlier. Table 7.5 quantifies the results in

terms of turbulent intensities and kinetic energy. The Runge-Kutta schelne predicts

higher ma.xinmn_ turbulent intensities and kinetic energy than the Gottlieb-Turkel

scheme, indicating better resolution of the flowfield. The simulation time presented

in the lable is the ralio of the time required 1o the time required for lhe GotlliebTurkel scheme. Consistent with the one-dimensional results, the I{ unge-Kutta scheme

is more computationally efficient, using 16.5 percent less ('Pl: time.

scheme =Ilma x =Urea x kT,,,_, , time

Gottlieb-Turkel 0.26033 0.18561 0.080193 1.00000

Runge-I,_utta 0.26266 0.20305 0.096392 0.83523

Table 7.5: Effect of numerical scheme on axisymmetric solution

7.3.5 Evaluation of Exit Zone Boundary Condition

The use of an exit zone outflow boundary has been advocated in the area of

computational a eroacoustics where proper resolution of very' low magnitude acoustic

waves is critical. It is not clear that an exit zone is required when the simulation is

not intended to capture acoustic waves. The need for a.nd effectiveness of the exit

zone boundary condition was tested by comparing solutions with and without this

boundary treatment. The computational grid was modified by adding an exit zone.

NASA/TM--2001-210716 105
a region of highly stretchedgrid cellsdownstreamof the previousoutflow boundary

(figure 7.7). Fifteen a,ddit,ionM grid planeswereaddedusinga ten percentgeometric

stretching factor. The simulations were run using tho Runge-I,_uttaschemewith a

sixth-order filter.

Entropy contours show some differences in the turbulent structures. Bul. it is not

clear that tho differences at'(" due to any reflecting waves. Turbulence levels given in

table 7.6 indicate very little difference between the two solutions.

boundary condition u,,_,,.,. _,,,,_,,. kT,,, _

Exit Zone 0.30181 0.20379 0.096595

No Exit, Zone 0.26266 0.20305 0.096:/.92

Table 7.6: Effect of exit zone boundary condition on axisvmnletric solution

7.3.6 Effect of the Sub-Grid Scale Model

The solutions obtained thus far have not used the sub-grid scale model. The3' are

in effect DNS solutions, t[owever, the grids used are not fine enough to resolve all

the turbulent scales. Scales smaller than the grid size are not computed and their

contribution is lost. Since these scales tend to dissipate the larger eddies, the solutions without the sub-grid model should predict larger more energetic eddies than a

corresponding LES solution. The use of the term DNS to describe these solutions is

inisleading because DNS implies that all the turbulent scales are computed. Therefore, these solutions will be referred to as "coarse grid DNS" solutions.

NASA/TM--2001-210716 106
An LES solution on the medium grid was run and comparedto its "coarsegrid

DNS'" count,erpart. Entropy cent.oursare comparedin figure 7.9. The additional

dissipation flom tile sub-grid model servesto damp tile large scale structures. This

behavior is evident by observit_g the i,itial vortex roll-up it_ t}le mixitlg laver. The tir.st

large eddies occur further downstream in the LES solution due 1o the eddy viscosi)y

fi'om the sub-grid model. Overall, the LES solution produces less turbulen) mixing.

7 percent less turbulent intensi(y and 14 percent less turbulent kinetic energy ().a.1)le

7.7).

small scale modeling u,_,,. v,,,.,. kT,,_.

without sub-grid model 0.26266 0.20305

with sub-grid model 0.24247 0.18797

0.096:}92

0.082770

Table 7.7: Effect of sub-grid scale model on axisymmetric solution

7.3.7 Evaluation of the Axisymmetric LES Solution

To obtain a steady time averaged solution, the LES solution was run for 100,000

iterations, which corresponds to 0.0052 seconds of simulation time. The instantaneous

solution was averaged every 100 iterations.

Analysis of the Flowfield

Contour plots of sample instantaneous and time averaged flowfield quantities are

shown in figures 7.10 through 7.12. The (urbulent structures present in the mixing

layer are clearly seen in the instantaneous contour plots. The density contours (figure

7.10) reveal the shock structure in the potential core. The radial velocity contours

NASA/TM--2001-210716 107
(figure 7.12)providegood visualizationof the vortical structures. The time averaged

plots show smoothly varying gradients ill tile mixing laver and are analogousIo a

ftowfield obtained from a RANS solution or an;dvtic lllel}-lod.

Plots of p, ,. v. and k, are show in figures 7.13 through 7.16. All the turbulent

quantilies exhibit a similar struciure. The detail plots of the nozzle exil show a

small concentrated unsteady region at, the nozzle lip caused by vortex shedding. This

region is followed by a region dominated by very small scale turbulence. Radial

velocity contours in lhis region (figure 7.17(a)) show lhat no large scale structures

exist,. But. a plot of the sub-grid scale turbulent kinetic energy, ra:t:, shows thai the

energy from the small scales peak in this region (figure 7.17(b)). Downstream of the

small scale region, the mixing layer becomes unstable and large vortical structures

begin to form and grow larger wilh increasing distance from the nozzle exil.

Comparison to Experimental Data

The time averaged solution is compared to the experimental data in figures 7.18

and 7.19. Tile experimental data exhibits the behavior of a typical jet [77, 78]. A series

of weak shock waves can be seen near the nozzle exit. The end of ti_e potential core,

as indicated by the start of the velocity decay' on the centertine, is at approximately

7.5 jet diameters. The predicted velocity on the jet: centerline captures the first

few shock waves near the exit. The renlaining shocks are not captured due to grid

stretching away from the nozzle exit,. The predicted centerline velocity does not decay

within the computational domain indicating that the a xisymmetric assumption has

constrained the mixing layer to lie above the jet, axis. Radial profiles of axial velocity

at ,rid a = 2. 4. 6, 8, 10, and 12 are shown in figure 7.19. The experimental data

NASA/TM--2001-210716 108
was taken both al)oveand belowthe jet axis. All thesedata are plotted to show the

a.symmetryin the da.ta..Pot'comparison the LES solution is reflected aboul the axis.

Agreement between the ('FD and experilnent is very good up to 6 .jet diameters.

Beyond this point the prediction departs from the experiment due t.o the lack of

decay of tile 1)otentia] core. The peak in the experimenlal profiles begins _o decrease

and the jet spreads a.t a greater rate lhan the CFD indicates. The good agreement

upstream of the end of tile potential core indicates that the' axisymmetric assumption

is valid in this region. In the region where the potential core breaks down. highly

three-dimensional turbulent struclures may t)e the cause.

7.4 Three-Dimensional Solutions

A large-eddy sinmlation of the full three-dimensional jet flowfield is a large undertaking. The number of grid points and time necessary for a solution severely limit

the size and number of calculations. The results from the axisymmetric solutions.

code validation, and one-dimensional error analysis were used to guide the simulation

process. The Runge-Kutta scheme has been proven superior to the Gottlieb-Turkel

scheme and was used for all 3D computations.

Three computations were performed. The first used the sixth-order filter and

provided time averaged flowfield data. The second calculation used the eighth-order

filter and was used to ascertain the effect of the solution filter and to obtain turbulent

statistics and a series of instantaneous velocity distributions. The third used the

sixth-order filter and a modified grid to obtain two point correlation information.

NASAfFM--2001-210716 109
7.4.1 Grid Generation

The axis\'mnwtric grid study showed an increase in resolution of tile lurbulenl

strttctures with increasing grid resolution. But. a trade between accuracy and computational cost must be made to keep tile simulation time reasonable. The resolution

used in the fine grid would create a prohibitively large 3D grid and lhe coarse grid

showed no visible turbulent structures. The medium grid (301 x 129) showed resolution of large scal_, tm'bulence and provides a. basis for a reasonably sized 3D grid.

:\lso. a reasonable sub-grid model should replace lhe effect of the unresolved terms.

]o form the three-dimensional grid, a grid plane corresponding to the axisymmetric

grid was positioned at every ten degrees around the jet axis resulting in a cylindrical

domain. Four additional grid planes were added overlapping the first four planes

to facilitate use of the overlap boundary condition. The final grid has dimensions

301 x 129 x 40 and contains 1.553.160 points. Figure 7.20 shows st.reamwise and

cross-stream grid planes.

The grid was modified slightly for the two point correlation work. The grid points

were redistributed locally so that a grid point was located at each correlation location

(table 7.2). The overall grid structure and size remained the same.

7.4.2 Specification of Boundary Conditions

Boundary conditions are specified in the same manner as done in the axisymmetric

case (section 7.3.2). No-slip walls are specified on the nozzle surfaces. Subsonic inflow

conditions are used for the nozzle plenum and freestream inflow boundaries. The

combination of extrapolation and subsonic outflow are used on the outflow plane to

allow for pressure variation and extrapolation is used on the upper boundary.

NASAJTM----2001-210716 110
Two additional t)oundavyconditions are necessaryfor the 3D grid. First lhe

overlap condition is usedI,ocreatea continuous doma.il_in the azimuthal direction.

Second.'thejel centerlineconsislsof a grid planethat is collapsedto a line. The pole

l)oundarv is used here,

7.4.3 Presentation of the Time Averaged Flowfield

The LES simulation was run until the time averaged centerline velocity profile

remained unchanged for over 5,000 ilera.liolJs. This required 50.500 i't.eralions using

an average t.ime sl.ep of approximately 48- 10 -_) seconds. [sing 16 processors on the

Silicon Graphics Power Challenge ma.chitle, the simulation required about two months

of calendar lime for completion. The solution was averaged every 100 i`teralions.

(:ontours colnparing instantaneous and time averaged velocil ies in bo`th the streamwise and cross-stream planes are presented in figures 7.21 - 7.26. The difference between the instantaneous and time average is striking. The instantaneous velocities

show large turbulent structures that. average to zero over tithe. Large azimuthal velocit.ies and the appearance of turbulent structures that cross the jet. axis near the

end of the potential core indicate a. highly three-dimensional flowfield.

Radial velocity and sub-grid turbulent kinetic energy near the jet lip are shown

in figure 7.27 and differ from the axisymmetric solution (figure 7.17). In the mixing

layer the vortica.1 structures in the 3D solution appear weaker because of the aD re]ief

effect not accounted for in the axisymmetric solution. Also, the 3D solution does not

show any resoh,ed vortices near the jet. lip. The sub-grid turbulent kinetic is much

higher, in the aD solution, indicating that the sub-grid model has dissipated this

motion. The increased effect of the sub-grid model may be caused by the large grid

NASA/TM--2001-210716 111
spacingin the azimuthal direction. This large grid spacingresults in a larger tilter

width, which increa_sesLheamount of eddy viscosity (equalion (5.19)).

(!ontoursof dilatation. _ are shownin figure 7 98 By rewriting the continuity

Oa"_ ' ....

equation (equation 2.[) we can see thai dilatation is related to the convection of

densiiv waves, which can be related to the conveclion of sound waves.

Ou, l(0p Op )

- e N + (r.1;)

The dilatation contours show that lhe sound for this jet would emanate from the

mixing laver Ileal" and the end of the potential core.

The time averaged velocities closely resemble both experimenta.1 data and RANS

calculations. But. clearly the LES solution is capable of providing nmch more insight

into the flow physics of the jet through the unsteady information.

The time averaged velocity profiles on the jet centerline are presented in figure

7.29. The LES solution predicts that the potential core is shorter than found experiinentallv indicating that the turbulent eddies are too energetic. The sub-grid model

may not be providing adequate dissipation of the large scale eddies. Adjusting the

constants to the model may correct this problem. The rate of velocity decay beyond

the potential core is very close to that of the experimental data.

Radial profiles of axial velocity at x/D:i = 2, 4, 6, 8, 10, and 12 are shown in

figure 7.4.5. At ac/Dj = 2 and 4 the 3D prediction shows less spreading than either

the experimental data or the axisymmetric prediction. This difference may be caused

by the truncation error due to the large grid spacing in the azimuthal direction.

The spacing in the azimuthal direction is much larger than the spacing in the axial

and radial directions and is the reason for the difference between the axisymmetric

NASA/TM--2001-210716 112
and :_Dsolutions. Sil,cethe grid resolution is llOt snffi(ienl t.o capture the turbulenf

structures in this region, the sub-grid model should provide the equivalent effecl

through increased eddy viscosity. Tile fact that the LES solution predicts less jet

sp_wa<lit_g than the experimel_l, i,dicates tl_at the st_b-gri<l model is t_ot adequat<'.

V_,:hile the agreement with t])e experimental data is far t]'om perfecl, lifts calculation is sut)erior to the other high Reynolds number LES calculations in the lit, erature [al aa].

7.4.4 Presentation of the Instantaneous Flowfield and Turbulent Statistics

A second nearly identical sinmlation was run to obtaiil turbulent statistics and instantaneous flowfield data not saved during the first, simulation. The order of tile filter

was changed from sixth to eighth to ascertain its effect. The simulation was again run

until the cenlerline velocity profile did nol change over a period of 5,000 iterations.

The time step was fixed at ,50.10-9 seconds and 70,200 iterations were required. The

solut, ion was averaged every 80 iterations. The increase in the filter's order reduced

the amount of numerical dissipation in the sohltion and increased the resolution of

the high wave number disturbances. This leads to increased turbulent mixing as

shown in the centerline velocity profile (figure 7.31). The increased turbulent structures are most likely the reason why more simulation time was required to obtain a

suitable average. In theory the sub-grid scale model should adjust to the change in

flowfield resolution and filter and grid independence should 1)e possible. The increase

in resolution of the turbulent eddies, with a, corresponding increase in [S[ should be

compensated for by an increase in eddy viscosity (equation ,5.19). It, is evident that

the current model and the choice of coefficients are not adequate for this problem.

NASA/TM---2001-210716 113
Instantaneonsvelocity contours at three times are shown in figures 7.:12- 7.37.

The contours indicate that the jet. flowfield is highly lhree dimensional with very

largescalestructures. Significant variation wilh lime is found not only ill licemixing

laver, bttt also in lhe potential core. Here the coreflow "flaps" and rotates about

theje_ axis. Also. Ihe shockstructure ','ariesconsiderablydue to the variation in the

boundary betweenthe core flow and inixing laver caused by the vortical motion in

the mixing layer.

Turbulent intensities and turbulent kinetic energy are presented in figures 7.:t8

and 7.:{9. The axial turbulent intensity peaks at _ = 0.4123 in the mixing laver

nea.r the nozzle lip. Unsteady motion is also found in the potential core, due to

the unsteady shock structure, and in the freestream, due to acoustic waves. The

radial and azimuthal intensity plots are lower in magnitude and peak at a value of

/i = 0.2002 at. the end of the pot.ential core. They have very similar structures except

for two very low intensity regions. The first is the unsteady radial component of

velocity induced by the shocks in the potential core. The second is a wide region of

unsteadiness in radial velocity at the downstream boundary.

Turbulent kinetic energy is plotted in figure 7.39. The peak values, k" = 0.1760.

occur in the mixing tayer just downstream of the nozzle lip. A second lower peak,

k" = 0.1210. occurs near the end of the potential core.

The anisotropy of the turbulence is easily seen in figure 7.40. The ratio of radial

to axial turbulent intensity, _/f, is shown in three plots with different contour levels

to separate the regions of isotropy and anisotropy. The vast majority of turbulent

structures are anisotropic with a peak ratio of 2.775. The majority of the flowfield

NASA/TM--2001-210716 114
has a ratio in tile range 1.'2 < &/i_"< 2.0 This resull hl<lk'atesihat the isotropic

a,ssuml>tionsof most RANS turbulence models a,re not apl)licable for this flow.

7.4.5 Presentation of the Two Point Correlation Data

iX simulation on the modified 31) grid was run using the sixlh-order filler and axial

and radial velocity data were saved for correlation. These velocilv COml)onenls can

be combined Io form the velocity vector making an angle lhet, a with the jet axis.

[ '_ = u cos 0 -4- _' sin 0

A time history of the velocity signals at t,he points for the fourth separation, 1/21)j,

are shown in figure 7.4.5. Three angles are examined 0, 45, and 90 degrees. There

is very little change in the velocity signal a,s it convecls from the upstream to the

downstream point. The average velocity and t.urbulent intensity for each angle is

given in table 7.8

0 0.65417 0.26473

45 0.46029 0.25601

90 -0.0032162 0.17076

Table 7.8: \_locit, y statistics a,t, the two point correlation location

Two Point Space Correlation

A plot of the two point correlation coefficients versus separation distance is shown

in figure 7.42. The correlation of the velocity signal decreases towards zero with

NASA/TM--2001-210716 115
increashtgseparationdistanceindicating that the averagesizeof tile turbulent structut'esis lessthan the inaxilnuln separationdistance(onejel diameter). The turlmlent

length scah'wasestimated usingequation (7.13). To evaluate tile integral, the data

was fit with a sixth-order polynomial using tile method of least squares [79]. The

polynolnial is of lhe form

(; .5 !J = c6,r + c5.v 4- c4:r 4 4- c3:v:_4- c2.v 2 4- cl.v 4- Co

The coefficients obtained in the curve fitting are found in table 7.!). The t_olynomial

coefficient

angle co C'l "2 c2, c_ c5 c6 12 ]|Orlll

t) 1.OOO 0.04124 -4.230 5.4,58 -1.940 -1.458 1.105 8.564 .10 -(;

.15 1.004 -0.01181 -5.249 6.480 - 1.7:36 -2.002 1.305 9,989 .10 -t:

',-10 1.003 0.03630 -1 t.13 16.4_ -2.633 -8.769 4.687 4.052 .10 -s

Tal)le 7.9: Curve fit coefficients for two point space correlations

was then integrated analytically t,o obtain the length scale. The results are shown in

table 7.10. The data at 0 degrees indicate that the extent of turbulent structures in

0 e/t)j

0 0.50161

45 0.37032

90 0.10422

Table 7.10: Turbulent length scales

this direct, ion are approximately a half of one jet, diameter. The smaller length scales

NASA/TM--2001-210716 116
obtained flom the 45 and 90degreedala arecausedby the relatively large region of

nega, tive correla.{,io_l coefficiel_ts for this da.t,a.. This effect was also noted l,v C'hu [73 I.

Two Point Space-Time Correlations

Plots of the two point sl}a.{'e-till'_e correla,tions for eat}} separation are shown }n

figure 7.4.5. For each set)aration there is a (Ielav time where the signals are clearly

correlated (a distinct peak in the correlation curve). This delay time increases with

increasing separation distance.

The exact location of lhe maxhnum correlation coefficient, cannot be determined

from the discrete data obtained in the analysis. To find this location, the correlation

data was again curve fit with a polynolnial (equation (7.19)). The maximum in

the curve is found by differentiating the polyilonlial equation. The location of the

pea.k in the correlation curves corresponds t,o t,he time required for a signal al lhe

upstream point to convect to the downstream point. The resulting convection speed.

['_., (equation (7.14)) is shown in figure 7.44. The convection speed decreases slightly

with increasing separation distance due to viscous effects. An average of lhe comimted

convection speeds is 0.6260[_ which is consistent with the results of Chu [73].

NASA/TM--2001-210716 117
coefficient

set). F 0 C 1 C2 t's C4 E'5 C(;

0 degrees

l 2 IIOI'lll

1 0.9517 0.5115 -1.585 0.4060 0.5007 -0.:{738 0.07827 9.67:3 .10 -6

2 0.8214 0.84,t4 -1.224 -0.04053 0.3132 -0.06564 -0.004293 1.296 .10 -4

3 0.6650 0.7884 -0.04775 -1.256 0.7415 -0.1015 -0.009149 1.333 .10 -4

4 0.5054 0.6288 0.7048 - 1,583 0.6886 -0.07523 -0.005351 2.419 .10 -4

5 0.3,178 0.,t978 0,9339 -1.208 0.2545 0.06518 -0.01998 2.083 .10 -'i

(i 0.2043 0.3583 1.027 -0.9395 0.1094 0.06236 -0.01272 6.915 .IU -_

7 0,07140 0.4665 0.06638 0,(i771 -0.8747 0.3248 -0.03882 1.119.10 -4

8 0.0052(;,t -0.01440 1.018 -0.3707 -0.1664 0.09207 -0.01089 2.219 .10 -:_

45 degrees

1 0.9332 0.6871 -2.083 0.5811 0.7198 -0.5646 0.1216 2,791 .10 -5

2 0.7569 1.1,t8 -1.619 -0.06277 0.5037 -0.1413 0.003079 3.132 .10 -4

3 (/.5481 1.097 -0.1815 -1.449 0.8416 -0.07731 -0.02293 5.145 .10 -4

4 0.3468 0.8275 0.9012 -1.943 0.7837 -0.04498 -0.01681 1.168 .10 -:_

5 0.1574 0.5952 1,208 -1.316 0.04382 0.2163 -0.04732 5.548 .10 -4

e; -0.003487 0.3590 1.333 -0.9435 -0.1417 ().1981 -0.03372 1.312 .10 -P'

7 -0.1408 -0.4605 0.1026 1.045 -1.271 0.4671 -0.05632 6.7(i2 .10 -4

8 -0.1878 0.2468 -0.3164 1.754 -1.508 0.4642 -0.04896 3.617 .10 -P

90 degrees

1 0.8634 1.365 -4.007 1.197 1.654 -1.23,t 0.2449 1,173 .10 -4

2 0.54,51 1.661 -0.3637 -6.358 7.044 -2.941 0.4,t17 2.954 .10 -4

3 0.2079 1.348 2.144 -6,379 4.328 - 1.096 0.08009 1.018 •10 -s

4 0.05870 0.6330 3.639 -5.640 2..585 -0.3642 -0.007533 3.972 .10 -3

5 -0.2577 0.3414 2.409 -1.468 -0.9120 0.7960 -0.1435 3,272 .10 -3

6 -0.3373 -0.2450 2.282 -0.4172 - 1.284 0.6962 -0.1013 5.808 .10 -3

7 -0.3705 -0.07766 -0.1704 3.000 -2,931 0,9961 -0.1147 6.802 • 10 -:_

8 -0.3132 0.03557 -2.236 5.467 -3.921 1.129 -0.1155 1.882 .lO -2

Table 7.11' Curve fit, coefficients fox" two point space-time correlations

NASA/TM--2001-210716 118
Figure 7.1" Two point correlation schematic

Figure 7.2: Axisymmetric computational grid

NASA/TM--2001-210716 119
Extrapolation

(Supersonic Outflow)

Subsonic Inflow

Freestream Conditions

Subsonic Outflow

Extrapolation

(Supersonic Outflow)

j.No Slip Wall

_-/Plenum Cond_'ons "4z_22::_/ Axis(Slipof Symmetry Wall) 0 (

Figure 7.3: Boundary conditions for axisymmetric calculations

NASAFFM--2001-210716 120
(a) (;,am;s,,_grM

L_j_/

-" _4'

c__i ........ , ] "_ , _r ?

(b) M_'diunl grid

(c) Fine grid

Figure 7.4: Entropy cont, ours for axisymmetric grid study

NASA/TM--2001-210716 121
>,

"E

0.35

0.3

0,25

0.2

0.15

0.002

i

L

0.004 0,006 0.008 0.01 0012

_X

(a) Turbulent intensities

0.14

0.12

0.1

0.08

0.06

0.04

0.002

{

' i

• . . i . . . i . . . i , . . i . . .

0.004 0.006 0.008 0.01 0.012

,3X"

(b) Turbulent kinetic energy

Figure 7.5: Turbulent statistics for axisymmetric grid study

NASA/TM--2001-210716 122
(a) Got _Jieb-Turk_,]

(b) Runge-Kutta

Figure 7.6: Enlropy contours for axisytn_netric _cheme c_m_t)aris(m

Figure 7.7: Axisvmmetric. computa_tional grid_ with exit "zone

NASA/TM--2001-210716 123
(a) wifll exit zone

(b) withoul exit zono

Figure 7.8: Entropy contours for exit zon.e boundary condition, coml:)arison

NASA/TM--2001-210716 124
................................................................................ _iii:_ ................ ...... !ii_:r:_i _

(a)withouts,b-_rklnJod_'l

(b) with sub-grid lnodel

Figure 7.9: EntroPy contours showing effect of l;he sub-grid model

NASA/TM--2001-210716 125
a) lnst, allt, alleOllS

," [/,_ "JL_ :x_F; !_ag"( _, •,_._Y, ..... _ - _ .... __-- ---_

(h) timo averaged

Figure 7.10: Densit, v contours for the a,xisynmlet, ric solution

(a) inst antaneous

(b) time averaged

Figure 7.11: Axial velocity contours for the axisymmetric solution

NASAfrM--2001-210716 126
(a)instantaneous

(b) time averaged

Figure 7.12: Radial velocity contours for the axisymmetric solution

J

J

(a) entire domain

(b) detail of nozzle exit

Figure 7.13: Root mean square density contours for the axisymmetric solution

NASA/TM--2001-210716 127
(a)entiredomain

\

(b)detailofnozzleexit

Figure 7.14: Root mean square axial velocity contours for the axisymmetric solution

(a) entire domain

(b) detail of nozzle exit

Figure 7.15: Root mean square radial velocity contours for the axisymmetric solution

NASA/TM--2001-210716 128
(a) olitire donlain

(t)) d_'i_dl of nozzl(' e,xil

lSigurt_ 7.16: TurbuleTlt ki_leiic er_crgy corltours for the txisviiilllolric solutioll

NASA/TM--2001-210716 129
2i'.

(a) i_stantaneous radial velocity

(b) sub-grid t,urbule_t kit_el,i,'.' energy

Figure 7.17: Nozzle exit. det.a.il R)r axl,. " "'s._,_nH_r_t " ,trl(" soluti(m

NASA/TM--2001-210716 130
1.2

0.8

_- 0_6

0,4

0.2

.... i .... I

o o i i

q J i

............. LES .................... _ .............................................

1 0 15 20

x/Dj

Figure 7.18: Time averaged centerline velocity profile for the axisymmetric LES solution

NASA/TM--2001-210716 131
12

• " ! ! .... ! .... 1 " "I

t

i o ¢o ° il _,

_ _ - _. ; u --....-4.-.J o Expt L-.-]

I ! O_ i I--LES I_

0.6

04

0.2

0

-0.2

o

---7 .......

o...................................... i.............

o

'%

.... , .... i .... i .... I ,

-2 -1 0 1 2

r/D_

(a) ./O_ = "2

"S

1.2

0.8 ........ }

i

J

0.6

O4

• i

i II

0.2 _ t i

I i,

0 i

-02., i

1

.... i .... t .... i .... i , ,

-2 -1 0 1 2

r/Oj

(b) ./r)j = 4

Figure 7.19: Tiine averaged radial profiles of axial velocity ['o1"the axisymmetric I,ES

solution

continued

NASA/TM--2001-210716 132
Figure 7.1.q:conlinued

1.2

1

0.8

0.6

04

0.2

0

-02

• • i .... i ! ' '

............ ............... l.............

• , i , , , , I , • , , t .... I .... I , ,

-2 -1 0 1 2

riD

(c) .r/D.i = 6

1.2

1

0.8

0.6

0.4

0,2

0

-0.2

• • i

i

i

t o Exo,/

-- LES r-

!

.. i .... i .... i, ,,., .... ,..

-2 -1 0 1 2

riD I

(d) ./Dj = 8

continued

NASA/TM--2001-210716 133
Figure 7.19: continued

12

1

0.8

06

04

0.2

0

-0.2

,,v, , ............

_ ..................... _"1--"_ lJ i o i i

....... ....................... i..................... _. _....................................

..I .... i .... i .... i ......

-2 -1 0 1 2

r/Dj

(_) x/Dj = 10

1.2

1

0,8

0.6

04

0.2

0

-02

_....... '= .... _ .... T ......

_- _ o e_, k

o0 =' " O i:/ / f_°o f '

1,,^°o i/ \i -o_ i '

.... , , I .... I , , . . i .... I . .

-2 -1 0 1 2

r/D_

(f) _/o_ = 12

NASA/TM--200 ! -210716 134
(a) streamwise, x- 3. plane

(1:)) crossst, ream, y-z plane at outflow boundary

Figure 7.20: Three-dimensional grid

NASA/TM--2001-210716 135
(a) instantaneous

(b) time averaged

Figure 7.21" Axial velocity contours for 3D LES solution

NASA/TM--2001-210716 136
J

(a) instantaneous

(b) time averaged

Figure 7.22: Radial velocity contours for 3D LES solution

NASA/TM--2001-210716 137
(a) IIISI alll alleOUS

_L

/

.-- _ _-_"-_-_=_'W_. _- _'" _ ',. C_- /_ :'" _, '_-_--

,/

(b) time averaged

Figure 7.23: Azimuthal velocity contours for 3D LES solut, ion

NASA/TM--2001-210716 138
"\_J

\

i

/

./

(a) instantaneous (b) tilne averaged

Figure 7.24: Total velocily contours al :r/Dj = 3 for 3D LES solution

_Ta_ _-:.>.. .....

illl! ,

/'

/.

/

/

/'

/

l

\,

.f

/-"

/

\

\

\

/

/

,/

/"

i//

(a) instant aneous (b) time averaged

Figure 7.25: Total velocity contours at a'/Dj = 6 for 3D LES solution

NASA/TM--2001-210716 139
(a)instantaneous (b)timeaveraged

Figure 7.26: Total velocity contoursat x/D3 = 9 for 3D LES solution

NASA/TM--2001-210716 140
/"

'i' y;

,\t?_._l I \,,"

) /")/!

( .-_._t I'i/-'-

(b) sub-grid iurbuleut kinel, i< energy

t:'igur(' 7.27: Nozzle exit; det:ail for 3I) I..ES sohit:ioii

NASA/TM--2001-210716 141
Figure 7.28: I)ilalation contours for 3I) LES solution

1.2

0.8

0.6

0.4

0.2

o - . \,00% j

oooa_

--LESo Expt.

,,,,I,, i i

5 10 15 2O

xlD

J

Figure 7.29: Time averaged centerline velocity profile for t,he 3D LES solution

NASA/TM--2001-210716 142
12

1

08

06

04

02

0

-0.2

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

I ' i '

o o {

o_ i

) i

i O i O Expt I

............. _ .......................................... 0 ............ "T......................... i--I -- LES I"

/

.......... i........................ ,-o................ _................. _........... !................................ _..........

....

-2 -1 0 1 2

r/Dj

(a) x/D.i = "2

1.2

1

0.8

0.6

04

0.2

0

-02

--4

! o

E'"

__4 .................. ÷.__.,_

io

f i .... i ....

-2 -1

I .... i .... i .... ! .... i • •

i , i

o _°_ ; I o E_p,/

.....Z__le

0 2

o i

o,i

1

r/Dj

(b) ,/D_ = 4

Figure 7.30: Time averaged radial profiles of axial velocity for the 3D LES solution

continued

NASA/TM--2001-210716 143
Figure 7.30: continued

"-3--

12

1

08

0.6

04

02

0

-0.2

.... , .... i " "

'i ........ !

I i :

..... •_,...................... +.--.-_, .......... _.......................... .L ............

i Jo i i

__+ z j_ '_

,,J .... i .... l,ltll .... I.,

-2 -1 0 1 2

r/Dj

(c) x/Dj = 6

._.1 ̧-

1.2

1

0.8

0.6

0.4

0.2

0

-02

''T .... ' .... ! .... ' .... '

_ o_

_ 0

i

- T-T-

: i

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

-2 -1 0 1 2

r/D

J

(d) x/D3 = S

continued

NASA/TM--2001-210716 144
Figure 7.30:continued

Z)--

1.2

1

08

0.6

0.4

0.2

0

-02

• ' 1 ! i , "I

= i

_ _oT

0

i o o

i o i o i

i b oi

__ ..... _ L_-----_-._: ........... 2---

........ _.......................... _................ i............................. +...........................................

I I i I I I I I I I I , i .... [ .... I • •

-2 -1 0 1 2

r/Dj

(e) ajl)j = lO

-0.2

1.2 , , , , t • "

- I o Expt.

't i /--'_Sll 1

o8 i - o"-?o 7 ..... _-

i o i o i i

i o

0.20.4 ,_......................

0 < o :_ i "

td

• , I .... I .... i .... i .... _ . .

-2 -1 0 1 2

r/D]

(f) */>0 = 12

NAS A/TM--2001-210716 145
1.2

0.8

0.6

0.4

0.2

i I o Exot.

A, oO_ o - i o i / -- LES, 6th order filter

!_'_P_ ......... -¢- ........ 1 ..... LES, 8th order filter

-- - o_% i

_',, _'_ _%o

5 10 15 20

xlD

J

Figure 7.31: (!omparison of sixth- and eighth- order filters on 3D LES solution

NASA/TM--2001-210716 146
= 0.0000 ta) r_,

Figu_'e 7.32: h._st,m)t_._eous _xial ve'locity corot.ours-[or the 3D LI'.IS soh_l;km

N ASAJTM_2001-210716

147
(a) _ = 0.0000

(1,) *D,

= 16.149

\\ \,

,

, _ --%--_,,, __

( ) \ ---J -_>

J ,i /

%,

I

i"

\,, i ,,/

"---" J i "jlt/

"i//" ( c,

/

/'

L' p_.

<'-; It

tD_

(c) _ = 3Z297

Figure 7.33: lustantaTieous radial velocity coTltotirs ['()r the 31) I,ES solutio_l

NASA/TM--2001-210716 148
.-:= v._ -

i

i

/

/

i \

a) tD1

= 0.0000

iI

/

!

i

/

/ i

i i

jl / °

f'

(b) _ = 16.149

-.\

\

.,--2-,_%cU,.d' ._ t _ ....

_-._.-¢ " / ,/ _, . '\ _ ,t [

tD,

l:igure 7.34: h:lstal_taneous azimuthal velocity contours for the 3D I,F,S soluiion

NASA/TM--2001-210716 149
tD_

(a) _ = o.000o

(b) _ = 16.149

'3 (c) _ = 32.297

Figure 7.35: Instantaneous total velocity contours at x/Dj = 3 for the 3D LES

solution

NASA/TM--2001-210716 150
(_'1) tl); = i).O000

tD 1

(b) _ = 1(5.149 (c) _ = '32.2.¢ff

Figur( 3 7.:16: Insta.ni,_meous total velocity contours a,t :r/Dj = 6 tbr the 3D LES

solu t ion

NASA/TM--2001-210716 151
tD,

(a) _ = 0,0000

/

/;

r

= 16.149

(c) _ = ",12.2_,_)-I 0

Figure 7..37: hlstantaneot,.s total velocity coilLours aJ :r/Dj = 9 Ibr the 3D LES solut, ioti -'

NASA/TM--2001-210716 152
i

i

(a) axial lurbul¢'nt mlensity

(b) radial turbulent, [ntf:,nsit, y

<:---.

(c) azimuthal turbulent int, ensity

Figure 7.38: Turbulent intensity contours for lhe 3l) I,I_;S solution

NASA/TM--2001-210716 153
Figure 7.39: Turbulent kinetic energy contours for the 3D LES solulion

NASA/TM--2001-210716 154
A

(a) _/f < 0._

A

(b) O.S < _/_ < _.2

A

(_) _/_ > 1.2

Figure 7.40: Ratio of axial to radial turbulent intensity for the 3D LES solution

NASA/TM_2001-210716 155
<-

14

60 80 1O0

t U/D

(a) 0 degrees

12

1

0.8

0,6

0.4

0.2

0

43,2

-0,4

0

/ I _ 1!

20 40 60 80 100

t Ui/D J

(b) 45 degrees

Figure 7.41" Velocity history for two point correlations

continued

NA SA/TM--2001-210716 156
Figure 7.41" continued

D

I I

-0.8 .......... _ ........ h . L__

0 20 40 60 80 100

t u_,

(e) 90 degrees

0.8

0.6

0.4

0.2

0

-0.2

-0,4

1.2 : , ,

' " ............. ................ i............................. 4................ 4..........

!.............. • '........ i..... ............. !...........

0.2 04 0.6 0.8 1 1.2

Ax/D

I

Figure 7.42: Two point space correlation coefficient

NASA/TM--2001-210716 157
,,.,,

1.2 L • ) "-_-'----7 ..... 1

! ,-r'-,--','_

! : .... 0,000 O

,_ : • i -0,2,°: 4

L \" ,_--- .... 01375o

/ \" \/\ '," _ " ', ..... 0.750D -_

:.. .... : ....... ,, ,

/ , /\,_ • ' , ,: , ...... 1.000Dj

, .' ', L i \,

°2F:"::-: ......... .......

0.5 1 1.5 2 2.5 3 3.5 4

tU/DJ }

(a) 0 degrees

o,_V :\\/:_ ....;-,:.:-:',::- ..... °_'_':':.....

• _ - 0.500 D ;

0,6

/ ,\\/_i,_ ', i "\'_\-, - .0875D I

°"I .... t

' t \ , " "', ........ : ....... 1.000D "

I / .' _ ,. , ,

-0.2 ' .'. ..... ; ................ :. --:- -_"-.: .... -'-,,d.,:-:_..-,--_,_,-,<=-.:

-o.__-------_--,---.... _.... i. .. i ......... i..... L--

0 O.S 1 1.5 2 2.5 3 3.5 4

t u_

(b) 45 degrees

Figure 7.43: Two point space-time correlation coefficient

continued

NASA/TM--2001-210716 158
Figure 7.43: continued

1.2 r-_ ......... _ _-_ V-_- " .... , .......... ]

.- 0.000 D

1 _ --- 0.125 D' " "_ "_ - 0.250 O j

o8 L _:,'" ..... 0,375D -

•' _ ',!,, - ' 0,500O_

' ,v.'\ ,;; /\ ': - o,875D

0.2_: . ;'i _ ,\:..',: '".'. _ " " ': ........

.\ • X , / .•

o ............. ..,....... ,.. ..j/1,.:

i , \ \ . / ." .

_,211 .:.....,.........- " "

-04 _ -"" ...._-,-,--;....... " ..... !

0 0,5 1 13 2 25 3 3.5

t U/D

I I

(c) 90 degrees

K 0.65

0.75 .......... _ 1 • e---! 0 degrees I ]

---_- .5_r.,, I .1

_ ! + SOaegr_s I I

0.7 ........ ............ ; ...............

0.55 i , , ,

0 0.2 0.4 O,S 0,8 1 1.2

&xlD.

I

Figure 7.44: Convection velocity

NASA/TM_2001-210716 159

CHAPTER 8

CONCLUSIONS AND RECOMMENDATIONS

A numerical method to simulate high Hevnolds nulnber jel flows was formulated

an([ applied to gain a better understanding of the flow i)hysics. Close attention is

paid to the source_ of error in such calculations and efforts were made to minimize

them whenever possible.

Large-eddy simulation is chosen as the most promising approach to model the

turbulent structures due to its compromise between accuracy and compulational expense. The filtered Xavier-Stokes equations are developed including a total energy

form of the energy equation. Sub-grid scale models are adapted Kom compressible

forms of Smagorinsky's original model.

The effect of using disparate temporal and spatial accuracy in a numerical scileme

was discovered through one-dimensional model problems. The lower order time stepping found in many schemes, such as the Gottlieb-Turkel scheme examined here,

causes the scheme to revert to the lowest order accuracy. A new uniformly fourthorder accurate numerical method was developed based on this work. The scheme

consists of a low-dispersion Runge-Kutta time stepping scheme with a central difference spatial operator. Solution filtering is used to maintain stabilit, y. Results in both

one- ancl two-dimensions showed this new scheme clearly superior to the second-order

NASA/TM--2001-210716 161
i1,tinl_-and fourth-orderin space(-;otllieb-Turkel scheme.The measure used to judge

the schemes was the computational el[iciency, the time requir(,d to reach a given level

of error.

The resulting flow solver was configured to run oi1 a shared memory parallel con>

purer. Poor computer architecture prohibited conclusive results, but 1he evidence

indicates that the code scales well up to 16 processors. Results from validation exercises show thal the code accurately reproduces both viscous (laminar fiat plate)

and inviscid (supersonic wedge and cone flows) flows. The validation exercises also

confirmed the increased accuracy of the new numerical scheme.

Numerous axisymmetric simulations were performed to investigate the effect of

grid resolution, numerical scheme, exit boundary conditions and sub-grid scale model

on lhe solution. While the axisymmetric assumption was not accurate for the jet

flowfield, valuable information was gained for use in the three-dimensiolxal calculations.

The three-dimensional calculations showed that this LES simulation accurately

captures the physics of the turbulent jet. The agreement with experimental data

relatively is good and is much better than results in the current literature. However,

there is still much room for improvement. The improved agreement over previous

work can be attributed to the new numerical scheme and the modeling of the nozzle

lip. A subsequent run using a higher-order solution filter indicated that the modeling

of the unresolved scales needs improvement.

Several techniques were used to gain a better understanding of the underlying

physics. A plot of dilatation was used to provide insight into the sound field by indicating the location of acoustic sources in the jet mixing layer. Turbulent intensities

NASA/TM--2001-210716 162
indicate that the turbulent structures at.this level of modelingare not isotropic and

this information could lend itself to the developmentof improvedsub-grid scalemodels for LEE aT_dtmq)ulencemodels_'brRANg simulations. :\ two point correlation

technique was used to quantify the turbulent, structures. Two point space correlations were used _o obtain a measure of the inlegral length scale, which proved to

1

be approximately :iDj. Two point space-tinw correlations were used to obtain the

convection velocity for tile turbuMlt structures. This velocity ranged from 0.57 to

0.71 l:.i.

There are several recommendations for further work.

The accm'acy of the simula.tions is highly dependent on grid resolution. Accurate

resolution of the shock structure in lhe potential core was f<mnd with the axisvmmetric

calculation. However, the large grid spacing ill tile azimuthal direction in the 3D

calculations diminished this. Further resolution in this direction may improve lhe

prediction of both the shock structure prediction and the mixing layer. A more

systematic study of grid resolution in all three directions is desirable, but may be

computationally prohibitive.

A change in the order of the solution filter drastically changed the centerline velocity decay. The change in turbulent mixing with an increase in resolution of scales

indicates that the sub-grid scale model is not accurately mimicking the effects of the

unresolved scales. Further research into improving the sub-grid scale models is necessary. A promising approach is the dynamic sub-grid model [26] which automatically

adjusts the coefficients based on the filter width.

The time required to run a simulation on this relatively simple geometry severely

limits the useflflness of LES. This limit, may be eased if a more efficient, numerical

NASA/TM--2001-210716 163
melhod was found. As shown, the time slepping scheme is typically tile factor lhal

limits t,he computa.tional elficiency. A high-order accuracy eflicienl lime stepping

scheme would allow faster l urn-around of solutions and more accurate answers on

finer grids.

While there is room for improvement, in accuracy, this research has shown thai

large-eddy simulation can be used, as is. to provide new insight and information about

high Reynolds number jet flows. The characterization of the turbulent strucl, ures, size,

conveclion speed, all(I degree of anisolropy, can be used to develop improved tools

for predicting the fluid mechanics and acoustics of jels.

NASA/TM--2001-210716 164
APPENDIX A

DERIVATION OF THE FILTERED EQUATIONS

A complete derivation of tile Favre filtered Navier-Stokes equalions used in Chapter 5 is presented.

A.1 The Filter

A filtering function (;, is used to separale large and small scale components. The

filtering operation applied to a function f is

,/= [_ G(x - _).f(_)d_ (A.I)

•./--'X

The function can then be decomposed into its resolved/filtered, f, and unresolved, f'

parts

f = f+.f' (A.2)

For most applications, the function is not specified. But., several constraints are

placed on the the function to ensure that the filter colnnmtes with the derivative.

o-7 0.?

- (A.3)

0.r Ox

NASA/TM--2001-210716 165
The constrahils are

I) C;(-() = (,'(_)

2) (;(_)d_ = J

3) (;(_) -+ o _,_ I_l -* c_

4) G(() is small outside

where _.Xis a characteristic width of the filter function. Where it has been necessary

to know the form of lhe filter funclion, researchers have tyl)ica]ly used eilller a box.

(;aussian. or spectral cutoff filter.

Favre (density) weighting is used in the filtering process. This allows for convenient

recovery of terms corresponding to the unfiltered equations.

/-- (A.4)

P

A.2 Continuity Equation

Fi[tering the continuity equation is a. straight forwa, rd process. The spatial filter

is first applied to the conthmity equation (2.1).

Op Op u i

0-7 + 0x_ - 0 (A.,5)

Since the filter commutes with the derivative equation (A.5) is rewritten as

Op O-p-ai

0-7 + Oxi - 0 (A.6)

Then Favre weighting is use(1 to recover an equation of the same form as the unfiltered

equation (2.1).

Or5 Ot_,

-- + - o (A.r)

Ot O.ri

NASA/TM--2001-210716 166
A.3 Momentum Equation

Tile filtering of the nlonlenguln equa.tion and definition of ti_e sub-grid scale sl tess

tensor is ])resented l)elow. The filtering operation is al)l)lied to Ill(' momentunl equat ion (2.2).

Using the the properly in (A.3) the equation is rewritten

O_ O_, Op O/r,j

+ -- + -- = " (A.9)

Ol O.r., O.r,- t)J'.,

Favre weighting is then a.pplied

__O?,a,:+ Or,_-_,i + Of, Oa,, (,.\. 10)

8t /)xj Oa'; O,r.i

where the filtered stress tensor is

_ 2 - -, (A.11) ,5ij - -5Paij,Sa, a-+ 2tt,5'_

The filtered stress tensor and the term _7_,uj are not in useable forms because they are

the filter of a product of two variables. We define a new resolved stress tensor as

2_ q; 9~ 7

0"i2 = --Sffaij , b/,--_-11,_%,_./, '

(A.12)

where _ = y(;F) and

= ( .+ az,.I

(A.13)

We cem then write (A. 11 ) as

(r,:j= 8-,:,+ (a-,:,- _'u)

NASA/TM--2001-210716 167
Similarly. wewrite

(:\.1,3)

\Ve define the sub-grid scale stress tensor as

(.\. 16)

and rewrite (A.15) as

t)_Iitt./ = pttittj + 7"<i (A._7)

Finally. we substitute (A.14) and (A.17)into (A.10

filtered momenl um equ alion

--+--+ _ •

Ot Oa'j Oa'+ O.rj Ox.i

A.4 Energy Equation

and obtain the final form of the

0

+ j0.,.--(% - &.') (A. 1s )

The filt.ering the energy equation is the most involved process. Because the total

energy form of the energy equation is used, several additional manipulations of the

resulting terms are required in order to recover terms for which there are sub-grid

scale models. Most work in compressible large-eddy simulations have used either

static or total ellthalt)y forms of the energy' equation. First. the filtering operation is

applied to the energy equation (2.6)

Ope_ Opuiet Ouip Oujc_i¢ Oqi

0---7-+ Ox,-- + Oz----7= Ox----7 Oxi (A.19)

Commuting the filter operation with the derivative we obtain

Opt t 0_, O_ O_i i 00,

0--7- + Ox,-- + Oxi - Oxi Ox, (A.20)

NASAfrM--2001-210716 168
W|ler_

Fax're weighting equation (A.20) yields

O?_, O?_, aT O_i; Ocl,

-- + _ + -- = - ( :\ "2"2)

Ot O.r, Oa'i O.r, O.r;

where

g, = g + 5"J,."k (:\.2::I)

As was done with the monwntuna equation we rewrite the filter of the product of two

variables to obtain useable forms.

(A.24)

p_-7_= i_._ + (_,: - pgii) (A.25)

¢ = q_+ (¢ - q_) (a._7)

where

OT

0, = -A: .-c- (A.28)

ax,

-o_

_ = -a::r- (A.2.o)

clxi

NASA/TM--2001-210716 169
and _' = /c(F). Substituthlg equations (A.24) -(A.27) into (..\.22) we obtain the

following form of the filtered energy (*quaJ,ion.

i,---i-+ ---7o.r, + o.,., - O:r, Oa,: e):,.:[/' (G,; - g,a:)]

• , ,/

T

0

The underbraced terms can be fitrther manipulated to obtain terms for which sub-grid

scale models have 1)een previously developed.

The argunwnl of the derivative in term (i) is transfornwd as follows.

(A.31)

The first terrll ill (A.31) can be written as the sub-grid scale heat flux.

(A.a2)

_ /:?

The second term in (A.31) is the sub-grid scale turbulent diffusion and is denoted as

_p

'- :D, - = (A.33)

Finally the argument of the derivative of term (i) in (A.30) is simplified to

R

fi (_ - _) - Q_ +/JDi (A.34)

'y-1

lJsing the filtered equation of state

p = pRT (A.35)

NASA/TM--2001-210716 170
the argument of the derivative in term (ii) in (:\.30) can be rewritten a,s the sub-grid

scale hea,t flux a,s wa,s done in equa,tion (A.32).

p_-7,.- pa, = pt_T,, - f)t?7'a,

(,,>- h :\.36)

= tlq_

Terrn (iii) in (A.30) is transformed as follows.

0

07,,_(_J - _a,a,j)_

A.:}7)

where _ is the sub-grid scale turbulent dissipa,tion rate

' (..3_)

{ = O'ij O3"i CriJ O:ri

Substituting (A.34) (A.36) and (A.37) into (A.30) and rearranging terms gives us

the final form of the filtered energy equation.

(A.39)

A.5 Determination of Pressure

Pressure is normally obtained from the total energy as follows

(A,40)

NASA/TM--2001-210716 171
Applyilig l,[le spatial filter and Favre weighting we obt,airl

(:\.41 )

= (_- J)(r,_,- 1⁄2/_,,)

The above equation for pressure contains the filter of a product of variables, iT_-.,a..

This undeterlnin¢'d quantity is eliminated b\' using the sub-grid scale kinetic energy

- (_ _)[_#, ' -- = '-(,;>-,_ a_:a_.)] <A.42)

_-(.. - l)(/_g,- 1⁄2/,a_.a_: - _-,..,.)'

NASA/TM--2001-210716 172
BIBLIOGRAPHY

[1] A. N. Kolmogorov. Local Structure of Turbulence in Incompressible Viscous

Fluid for Very Large Reynolds Number. Doklady Al,'adcmiga Xaul; ,5',5',5'R, 30:2.q.9

303, 1941.

[2] D. C. Wilcox. Y_l.rbuhr_cr :llodeliT_g for CFD. DCW Industries, second edition.

1998.

[:3]

[4]

[5]

[7}

[8]

[9]

W. P. Jones and B. E. Launder. The l}rediction of Laminarization with a TwoEquation Model of Turbulence. lnterlmtional Jourr_al of Heal al_d Mass TraT_,@ P.

15:301-314, 1972.

[10]

S. Sarkar, G. Erlebacher, M. Y. ttussaini, and H. O. Kreiss. The Analysis and

Modeling of Dilatational Terms " ein Compresslbl_ Turbulence. NASA CR 181959.

1989.

O. Zeman. Dilatational Dissipation: The ('.oncept and Application in Modeling

(!ompressible Mixing Layers. Physics of Fhdds A, '2('2):178-188, 1990.

S. B. Pope. An Explanation of the Turbulent Round-Jet/Plane-Jet Anomaly.

AIAA dourrml, 16(3):279-281, 1978.

T. J. Barber, L. M. Chiapetta, J. R. DeBonis, N. J. C,eorgiadis, and D. A. Yoder.

Assesment of the Parameters Influencing the Prediction of Shear Layer Mixing.

Journal of Propulsion a_d Power, 15(1):45-53, 1998.

R. Courant. K. O. Friedrichs, and H. Lewv. On the Partial Difference Equations

of Mathematical Physics. IBM Jourlml of Research a_d Devdopm_l_t. 11:215-

234, 1967.

D. Got tlieb and E. Turkel. Dissipative Two-Four Methods for Time-Dependent

Problems. Mathematics of Computatiol_, 30(136):703-723, 1976.

S. K. Lele. Compact Finite Difference Schemes with Spectral-Like Resolution.

Journal of (!omputational Physics, 103:16-42, 1992.

NASA/TM--2001-210716 173
[11] C. I(. W. Tam and .l. ('. \Vebl_.DispersionRelation-PreservingF'ilti**,l)ifferenc_,

Schem<'sfor (:omputalional Aeroacouslics. Jo_lr,,I of ('ompulaliol_,l P/_t.,ic.,,

107:262 2S1, 1993.

[12] .\I. R. Visbal and D. \7. (',aitonde. tligh-Order Accurate Methods for lrnsteady

Vortical Flows on Curvilinear Meshes. AIAA Paper 95-0131, 1995.

[13] R. Hixon. A new class of compact schemes. AIAA Paper 95-0367. 19!}5.

[14] R. tlixon. Nonlinear (:omparison of High-Order and Optimized Finile-Difference

Schemes. NASA CR 19.98-208670, 1995.

[15] |1. Hixon and E. Turkel. High-Accura.cy ('ompact MacCormack-q'ype Schen_es

for Computational Aeroacouslics. NASA CR 1998-205672. 1995.

[17]

[lS]

[19]

[2o]

.l.B. Freund and P. Moin. Mixing Enhancelnent in .let Exhaust [sing Fluidic Actuators: Direct Numerical Simulations. FEDSM 98-5235, ASME. 19.q5.

Proceedings of FEDSM'98.

J. B. Freund, S. K. Lele, and P. Moin. Direct Sinmlation of a Mach 1.92 .lel and

its Sound Field. AIAA/CEAS Paper 95-2291, 1998.

.1. B. Freund. Direct Numerical Simulation of the Noise fl'om a blach 0.9 .let.

FEDSM 99-7251, ASME, 1999. Proceedings of FEDSM'99.

J. H. Ferziger and D. C. Leslie. Large Eddy Simulation: A Predictive Approach

to Turbulenl Flow Computation. AIAA Paper 79-1471, 1979.

R. S. Rogallo and P. Moin. Numerical Simulation of Turbulent Flows. ,41_1_uaI

Review of Fluid :Uechanics, 16:99-137. 1954.

[21] J. H. Ferziger. Recent Advances in Large-Eddy Simulation. In W. Rodi and

G. Begeles, editors, ET_gineeri_g Turbulence Modeli,n 9 a,d E, vp_r'im_nl._ 3, page

163. Elsevier. 1996.

[22] M. Lesieur and O. Metais. New Trends in Large-Eddy Sinmlations of Turbulence.

Amzual Review of Fluid Mechanics, 28:45, 1996.

[23] lJ. Piomelli. Large-Eddy Simulations: Present State and Future Directions.

AIAA Paper 98-05:34. 1995.

[24] J. Smagorinsky. General Circulation Experiments with the Primitive Equations.

Part I: The Basic Experiment. Monthly 14_ather Review, 91:99-152, 1963.

[25] F. M. White. Viscous Florid Flow. McGraw Hill, 1974.

NASA/TM--2001-210716 174
[26]

[27]

[2s]

[29]

[ao]

[a2]

[aa]

[a4]

[a6]

[3r]

[as]

[a9]

M. Germane. U. Piomelli. P. Moin. and W. (!abet. A 1)vnamic Subgrid-Scale

Eddy Viscosity Model, Pb!].,i('._ ql' b'l,id._ A, 3:1760, 1991,

J. Bardina, J. H. Ferziger. and W. C. I{evnohls. Improved Subgl'id Models for

Large Eddy Simulation. A1AA Pal_er 80-1357, 1980.

P. h[oin. K. Squires, W. (:abet, and S. Lee. A Dynamic Subgrid-Scale Model for

(:ompressible Turbulence and Sca.lar Tra.nsport. Phqsic._ o.f Fluids ,4.3( 11 ):2746

2757, 1991.

Y. Zang. R. I,. St reel, and J. t{. I,[oseff. A Dvuamic Mixed Subgl'id-Scale _lodel

and its Application to Turbu[eut Recirculating Flows. Phgsic., of Fluids .4.

5:3186 3196, 1991.

C. C. Nelson..S'imulalion._ of Spatiall.q E_,oh,in 9 ('ompre._sibl_ Turbuhnce l:si1_ 9

a Local D_amic Sabgrid Model. PhD thesis, Georgia Institute of ]-echnology.

1997.

R. R. Mankl)adi, S. H. Shih. and R. Hixon. Direct Comlmtation of Sound Radiation by ,let Flow Using Large-Scale Equations. NASA TM 106877. 1995.

L. Gamet and J. L. Estivalezes. Application of Large-Eddy Sinmlations and

KirchoffMethod to Jet, Noise Prediction. AIAA Journ(d, 36(12):2170-2178, 19!)8.

D. Choi. T. J. Barber, L. M. Chiapetta, and M. Nishimura. Large Eddy Simulation of High-Reynolds Number ,let Flows. AIAA Paper 99-0230. 1999.

R. W. MacCormack. The Effect of \;iscosity in Hypervelocity Impact Cratering.

AIAA Paper 69-354. 1969.

P. D. Lax and B. Wendroff. Systems of (lonservation Laws. Commu_ications o_

Pure and Applied Mathematics, 13:217 2::17, 1960.

A. Bayliss, P. Parikh, L. Maestrello, and E. Turkel. A Fourth-Order Scheme for

the Unsteady Compressible Navier-Stokes Equations. NASA CR 177994. 1985.

A. Jameson and T. J. Baker. Multigrid Solution of the Euler Equations for

Aircraft Configurations. AIAA Paper 84-0093, 1984.

M. H. Carpenter and C. A. Kennedy. Fourth-Order 2N-Storage Runge-Kutta

Schemes. NASA TM 109112, 1994.

F. Q. Hu, M. 5'. Hussaini. and .1. L. Manthey. Low-Dissipation and LowDispersion Runge-Kutta Schemes for Computational Acoustics. ,]o_rl_al of Computatio_al Phg.qics, 124:177-191, 1996.

NASA/TM--2001-210716 175
' 1

[40] S. (,o llieb and ('. W. ('hu. Total Variation Diminishing Runge-I',:u'_ta Schemes.

NASA ('R 201.'5.9l, 1996.

[41]

[44[

D. Stanescu and W. G. ttabashi. 2N-Slorage Low l)issipalion and Dispersion

Runge-Kutta Schemes for (!omputational Acouslics. Jourlml of ('omputatiomd

Ph!/sics. 1-13(12):6,4- 6,1 1995.

_{. Vichnevet.sky and .I.B. Bowles. bburicr AT_al.vsis of Num_cic, I .4pprozimatio_s of Hgperbolic Equation,s. Soci_'ly for Applied and Industrial Mathematics,

1.:)8'2.

('. K. W. Tam. Computational Aeroacoustics: Issues and Methods. .4IAA

Journal, 33(10), 1995.

T, It. Pulliam. Artificial I)issipation Models for the Euler Equations. A/AA

,]ourlwl, 24(1:2):1932-1940. 1986.

[45] G. Dulikravich and D. Dorney. Artificial Dissipation Sensors for (:omputational

(,asd_namlcs. A1AA Paper 89-06z13, 1989.

[46] ('. tlirsch. Numerical Computaiion of lT_l_rmd (l_d E,rter_al FIo,,s, volume 2,

pages 344-395..John Wiley & Sons, 1990.

[47] E. Turkel and V. N. Vatsa. Effect of Artificial Viscosity on Three-Dimensional

Flow Solutions. AL4A JourT_al. 3_(1):./9-4o. 1994

[48] C. A. Kennedy and M. H. Carpenter. Comparison of Several Numerical Methods

for Simulation of Compressible Shear Layers NASA TP 3484, t997,

[49] R. Hixon. Prefactored Compact Filters for Computational Aeroacoustics. A1AA

Paper .)9-03:_. ( "" _ 1999.

[50] G. B. Whitham. Linear a_d Nonlinear Wave_s. ,lohn Wiley &" Sons. 1974.

[51] R. D. Snyder and J. N. Scott. Comparison of Numerical Schemes for the Analysis

of Aeroacoustics. AIAA Paper 99-0354, 1999.

[52] B. Vreman, B. Geurts. and H. Kuerten. Subgrid-Modelling in LES of Compressible Flows. In P, R. Voke et al., editor, Direct a'nd Large-Eddp Simulation I,

pages 133-144. Kluwer Academic Publishers, 1994.

[53] E. R. Van Driest. On Turbulent Flow Near a Wall. JourTml of the Aeronautical

Sciences, 23:1007-1011, 1956.

[,54] G. Erlebacher, M. Y. Hussaini, C. C. Speziale, and T. A. Zang. Toward the

Large-Eddy Simulation of Compressible Turbulent Flows. NASA CR 187460,

1990.

NASA/TM--2001-210716 176
[55] D. A. Anderson, J. C. Tannehill, and R. H. Pletcher. ('omputalioTml Fluid

:lle:cham:c,_ ol_d Heat Tr_m.@r. Taylor & Francis, 1._)8.:1.

[5(i] M. Vinokur. ( onser_ation Equations of Gas-Dynamics iTl ("urvilinear ('oordinate

Systems..lottrmzl o./" (*ompuhltiolm/ Ph!Isics. 1.1:10:5-12.5. t.cJ74.

[57] R. Hixon, S. H. Shih. T. Dang, and 1R. R. Mankbadi. Evaluation of Generalized Curvilinear Coordinate Transformat.ioxts Applied 1o High-:\ccuracy l:inileDifference Schemes. AIAA Paper 98-0370. 1998.

[88] R. W. _ 'lacCormack. Numerical S'olutio_ of tae I,,t_ ,,oceio,, of a ,s'/,ocZ, ii )t_,( _r,ilh

o Lami_ar' Bou_MaJ' 9 Log_r, volume 8 of Lcctuc_ :Votes i_ Pl_!lsics. pages 151

163. Springer-Verlag, 1971. Proceedings of the Second lnternal ional ('onference

on Numerical _lethods in Fluid Dynamics.

[59] M. Ft. Giles. Nonreflecting Boundary Condilions for Euler Equation (!alculali<ms.

AIAA Jourr_al, 28(i2):9050 2058, 1{190.

[60] R. Hixon. S. H. Shih. and 1R. R. Mankbadi. Evaluation of Boundary ('onditions

tbr Computational Aeroacouslics. NASA TM 10664o, 1995.

[61] T. Colonius, S. K. Lele. and P. Main. Boundary Conditions for Direct ('ompul, a,tion of Aerodynamic Sound Generat, ion. AIAA Journal, 31(9):1574-1582,

1993.

[62] J. B. Freund. Proposed Inflow/Outflow Boundary Condition for Direct Computation of Aerodynamic Sound. AIA.4 Jour_ml, 35(4):740 742, 1997.

[63] M. J. Quinn. Designin9 Efficient Algorithms for Parallel Computers. McGrawHill. 1987.

[64] D. P. Bertsekas and J. N. Tsitsiklis. Parallel and Distributed ('ompulalion. Prentice Hall, 1989.

[65] ,l. H. Ferziger and M. Perid. (?omputalional M_thods for Fluid Dynamics.

Springer-Verlag, 1996.

[66] tI. Schlicting. Boundarg Layer Theo W. McGraw-Hil, seventh edition, 1979.

[67] Ames Research Staff. Equations, Tables and Charts for Compressible Flow.

NACA Report 1135, 195:1.

[68] J. Panda and R. G. Seasholtz. Density measurement in underexpanded supersonic jets using rayleigh scattering. AIAA Paper 98-0281, 1998.

NASA/TM--2001-2 !0716 177
[6.()1

[7o]

ITt]

[74J

[77]

[rs]

.I. Panda and R. G. Seasholtz. Velocity and temperature measurement in supersonic fl'ee jels using spectrally resolved rayleigh scattering. AIAA Paper 99-0296,

1999.

.J. Panda and R. G. Seasholtz. Density fluctuation measurement in supersonic

fully expanded jets using rayleigh sca tlering. AIAA Paper .(I.(I-1870, 1.q99.

J. Panda. and R.. G. Sea.sho[t,z. Invest.iga.tion of density fluctual ions in supersonic

free jets and corrleation with generated noise. AIAA Paper 2000-2099. 2000.

H. Tennekes and J. L. Lumlev. A First ('our.s_ i7_ 7'urb_de_cc. MIT Press, 1972.

W. T. (:hu. Turbulence mea,surelnent, s Relevant to .lel Noise. ITT1AS Report

119. 17niversity of Toronlo Institute for Aerospace Sludies, 1966.

J. N. Scott. A ('omparison of" Numerical and Experimental Results for 17nst, eadv

Flow Associated with ,Jet Noise. AIAA Paper 94-0459, 1994.

Gridgen Version 13. l!ser manual, Pointwise. Inc.. 1998.

K. A. Hofl'mann. Computational Fluid D.tlnamics for E_gi_eers. Engineering

Education System, 1989.

G. N. Abramovich. The The:or, o,l"Turbulent Jets. MIT Press, 1963.

P. O. Witze. Centerline velocity decay of compressible fl'ee jets. AL4A Jourmd.

12(4):417--4186, 1974.

[7t)] E. Kreyszig. Advanced Engirzcering Mathematics. John Wiley & Sons, seventh

edition, 1993.

NASA/TM--2001-210716 178

REPORT DOCUMENTATION PAGE FormApproved

OMB No. 0704-0188

Public reporling burden for this collection of informalion iSestimated to average 1 hour per response, includingthe time for reviewing instructions, searching existing data sources,

gathedng and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of lhis

collection of information,includingsuggestions for reducing this burden, to Washington Headquarters Services, Directoratetot 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

May 2001 Technical Memorandum

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

The Numerical Analysis of a Turbulent Compressible Jet

6. AUTHOR(S)

James R. DeBonis

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

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field

Cleveland, Ohio 44135- 3191

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

National Aeronautics and Space Administration

Washington, DC 20546-0001

WU-708-90-43-00

8. PERFORMING ORGANIZATION

REPORT NUMBER

E-12669

10. SPONSO RING/MO_IITORING

AGENCY REPORT NUMBER

NASA TM--2001-210716

11. SUPPLEMENTARY NOTES

This report was submitted as a dissertation in partial fulfillment of the requirements for the degree Doctor of Philosophy

to the Graduate School of the Ohio State University, Cleveland, Ohio, 2000. Responsible person, James R. Debonis,

organization code 5860. 216--433-6581.

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Categories: 02, 07 and 34 Distribution: Nonstandard

Available electronically at htto://¢ltrs,m'c.nasa.¢ov/GLTRS

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

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

A numerical method to simulate high Reynolds number jet flows was formulated and applied to gain a better understanding of the flow physics. Largeeddy simulation was chosen as the most promising approach to model the turbulent structures due to its compromise between accuracy and computational

expense. The filtered Navier-Stokes equations were developed including a total energy form of the energy equation. Subgrid scale models for the

momentum and energy equations were adapted from compressible forms of Smagorinsky's original model. The effect of using disparate temporal and

spatial accuracy in a numerical scheme was discovered through one-dimensional model problems and a new uniformly fourth-order accurate numerical

method was developed. Results from two- and three-dimensional validation exercises show that the code accurately reproduces both viscous and inviscid

flows. Numerous axisymmetric jet simulations were performed to investigate the effect of grid resolution, numerical scheme, exit boundary conditions

and subgrid scale modeling on the solution and the results were used to guide the three-dimensional calculations. Three-dimensional calculations of a

Mach 1.4 jet showed that this LES simulation accurately captures the physics of the turbulent flow. The agreement with experimental data was relatively

good and is much better than results in the current literature. Turbulent intensities indicate that the turbulent structures at this level of modeling are not

isotropic and this information could lend itself to the development of improved subgrid scale models for LES and turbulence models for RANS simulations. A two point correlation technique was used to quantify the turbulent structures. Two point space correlations were used to obtain a measure of the

integral length scale, which proved to be approximately 1/2 Dj. Two point space-time correlations were used to obtain the convection velocity for the

turbulent structures. This velocity ranged from 0.57 to 0.71 Uj.

14. SUBJECT TERMS

Large eddy simulation; Jet flow; Turbulent jets

17, SECURITY CLASSIFICATION

OF REPORT

Unclassified

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIRCATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

198

16. PRICE CODE

20. LIMITATION OF ABSTRACT

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102