NASA/TM–2005–213771

Assessment of an Explicit Algebraic

Reynolds Stress Model

Jan-Rene ́e Carlson

Langley Research Center, Hampton, Virginia

December 2005
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated

to the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a

key part in helping NASA maintain this

important role.

The NASA STI Program Office is operated

by Langley Research Center, the lead center

for NASA’s scientific and technical

information. The NASA STI Program Office

provides access to the NASA STI Database,

the largest collection of aeronautical and

space science STI in the world. The Program

Office is also NASA’s institutional

mechanism for disseminating the results of

its research and development activities.

These results are published by NASA in the

NASA STI Report Series, which includes the

following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results

of NASA programs and include extensive

data or theoretical analysis. Includes

compilations of significant scientific and

technical data and information deemed to

be of continuing reference value. NASA

counterpart of peer-reviewed formal

professional papers, but having less

stringent limitations on manuscript length

and extent of graphic presentations.

• TECHNICAL MEMORANDUM.

Scientific and technical findings that are

preliminary or of specialized interest, e.g.,

quick release reports, working papers, and

bibliographies that contain minimal

annotation. Does not contain extensive

analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION.

Collected papers from scientific and

technical conferences, symposia, seminars,

or other meetings sponsored or

co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to

NASA’s mission.

Specialized services that complement the

STI Program Office’s diverse offerings

include creating custom thesauri, building

customized databases, organizing and

publishing research results . . . even

providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home

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

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA STI Help Desk

at (301) 621–0134

• Phone the NASA STI Help Desk at (301)

621–0390

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076–1320
NASA/TM–2005–213771

Assessment of an Explicit Algebraic

Reynolds Stress Model

Jan-Rene ́e Carlson

Langley Research Center, Hampton, Virginia

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23681-2199

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

constitute an offical endorsement, either expressed or implied, of such products or manufacturers by the

National Aeronautics and Space Administration.

Available from:

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS)

7121 Standard Drive 5285 Port Royal Road

Hanover, MD 21076–1320 Springfield, VA 22161–2171

(301) 621–0390 (703) 605–6000
Assessment of an Explicit Algebraic Reynolds

Stress Model

This study assesses an explicit algebraic Reynolds stress turbulence

model in the in the three-dimensional Reynolds averaged NavierStokes (RANS) solver, ISAAC (Integrated Solution Algorithm for

Arbitrary Configurations). Additionally, it compares solutions for two

select configurations between ISAAC and the RANS solver PAB3D.

This study compares with either direct numerical simulation data,

experimental data, or empirical models for several different geometries

with compressible, separated, and high Reynolds number flows. In

general, the turbulence model matched data or followed experimental

trends well, and for the selected configurations, the computational

results of ISAAC closely matched those of PAB3D using the same

turbulence model.

1 Introduction

The capabilities of application–level three–dimensional Reynolds

averaged Navier-Stokes (RANS) solvers have increased due to the

rapid growth of the speed and the size of computational resources

concurrent with advances in turbulence modeling. Requirements of

more completely resolving wake–boundary layer flows, highly curved

flows and jet shear flow physics, for example, have provided impetus

for the implementation of more, higher order, turbulence models

into RANS solvers. Two equation turbulence transport models,

such as k-ε or k-ω, though a step beyond the previous 0-, 1/2-,

and 1-equation models, still make significant compromises in the

physical modeling in the flow calculations.

Many current three–dimensional Navier-Stokes methods utilize the

higher order turbulence models such as full Reynolds stress model of

Launder et al;

1

the nonlinear model by Shih, Zhu, and Lumley;2 and

1 B. Launder, G. Reece, and W. Rodi.

Progress in the development of a

Reynolds stress turbulence closure.

J. Fluid Mech., 68:537–566, 1975.

2 T.-H. Shih, J. Zhu, and J. Lumley. A

New Reynolds Stress Algebraic Model.

NASA TM–166644, 1994. the explicit algebraic Reynolds stress model (EASM) of Gatski and

Speziale.3These models have had only a minor additional impact to 3 T. Gatski and C. Speziale. On Explicit Algebraic Stress Models for Complex Turbulent Flows. J. Fluid Mech.,

254:59–78, 1993

the overall computational effort and, for the most part, have been

an improvement over linear eddy viscosity turbulence models.

It is important to be able to evaluate the turbulence models within

a numerical code whose numerical characteristics are well-known.

In an attempt to quantify potential numerical biases present in the

code, a code-to-code comparisons is performed. The challenge is to

be sure that the turbulence model is similarly implemented in each

code. This level of effort is time consuming but is becoming more

essential in the current technical environment where an increasing

number of validation studies are undertaken.

1
The purpose of this paper is to describe and assess the implementation of a particular explicit algebraic Reynolds stress model

(EASM) in the RANS solver ISAAC (Integrated Solution Algorithm

for Arbitrary Configurations). Earlier results using this model in

the solver PAB3D are published in papers by Hamid,4 and Carlson.5

4 K. Abdol-Hamid. Implementation

of Algebraic Stress Models in General

3D Navier Stokes Method (PAB3D).

NASA CR–4702, December 1995.

5 J.-R. Carlson. Prediction of High

Reynolds Number Flow Using Algebraic Reynolds Stress Turbulence Models, Part 1: Incompressible Flat Plate.

J. Propulsion and Power, 13:610–619,

1997; and J.-R. Carlson. Prediction

of High Reynolds Number Flow Using Algebraic Reynolds Stress Turbulence Models, Part 2: Transonic ShockSeparated Afterbody. J. Propulsion

and Power, 13:620–628, 1997.

The first two sections briefly describe the ISAAC code and

details of the turbulence model equations used herein. The next

section contains the validation cases of the EASM simulating four

different flows and comparing them with the DNS (Direct Numerical

Simulation), empirical or experimental data where appropriate.

The last section are N-version tests (code-to-code comparisons) of

numerical solutions of ISAAC compared with PAB3D, using two

of the four test cases. The nomenclature is found in Appendix A

followed by a discussion on turbulence transition indicators used for

this study.

2 Numerical approach

Numerical algorithm

ISAAC 6is a multiblock, upwind, finite volume code that solves 6 J.H.Morrison. A Compressible

Navier-Stokes Solver with TwoEquation and Reynolds Stress Turbulence Closure Models. NASA CR–4440,

May 1992.

the three-dimensional Favre averaged compressible Navier-Stokes

equations with an upwind finite volume formulation. The inviscid

terms are upwind-biased spatial differenced and the viscous terms

are centrally differenced. This study also uses the flux-difference

splitting procedure of Roe7 with MUSCL8(Monotone Upstream- 7 P. L. Roe. Approximate Riemann

Solvers, Parameter Vectors, and Difference Schemes. J. Comp. Phys.,

43:357–372, 1981

8 B. van Leer. Towards the Ultimate

Conservative Difference Schemes V.

A second Order Sequel to Godunov’s

Method. J. Comp. Phys., 32:101–136,

1979

center Scheme for Conservative Laws) differencing and it advances

solutions in time with an approximately factored three-factor scheme.

The limiter used in the MUSCL scheme was Venkatakrishnan’s

limiter,9 unless otherwise noted.

9 V. Venkatakrishnan. Convergence to

Steady State Solutions of the Euler

Equations on Unstructured Grids with

Limiters. J. Comp. Phys., 118:120–130,

1995

k-ε equations transport equations

The transport equations for the turbulent kinetic energy, k =

1

2

u′·u′, where u

′

is the fluctuating velocity vector, and the turbulent

dissipation rate, ε, are10

10 Several variations of the dissipation

rate transport equation exist in the literature. Many of the variations involve

the destruction of dissipation term and

with the singularity of the term as k

approaches zero at the wall. The current model uses the model proposed

by Hanjali ́c and Launder. (K. Hanjali ́c and B. Launder. Contribution

towards a Reynolds-stress closure for

low-Reynolds-number turbulence. J.

Fluid Mech., 74:593–610, 1976.)

D

Dt (ρk) = ρPk − ρε + Dk (1)

D

Dt (ρε) = Cε1

f1

ρε

k

Pk − Cε2f2

ρεε ̃

k

+ Dε (2)

with

Cε1 = 1.44, Cε2 = 1.92, σk = 1.0, σε = 1.3 (3)

and where D/Dt = ∂/∂t+ui∂/∂xi. The equation for the production

of turbulent kinetic energy, Pk, uses the Reynolds stresses and the

2
mean flow gradients,

Pk = −τij

∂ui

∂xj

(4)

The term Dk derives its form from the turbulent convection and

pressure transport terms of the Reynolds stress equations. Often, as

is the case presently, the gradient-diffusion model of Daly–Harlow11

11 B. Daly and F. Harlow. Transport

Equations of Turbulence. Phys. Fluids

B, 11:2634–2649, 1970.

is used resulting in

Dk = ∇ · μ +

μt

σk



∇k



(5)

The dissipation rate transport equation uses a fully modeled term,

Dε, for the viscous diffusion and turbulent transport of ε.

Dε = ∇ · μ +

μt

σε



∇ε



(6)

The eddy viscosity, μt, in these two transport equations is the high

Reynolds number form with the equilibrium model coefficient Cμ0

,

μt = Cμ0

ρk2

ε

, Cμ0 = 0.09 (7)

The functions in the dissipation rate equation are

f1 = 1, f2 = 1 − 0.3e

−Rt

2

, Rt =

ρk2

με

(8)

and the singularity in the destruction of dissipation term is removed

through the term  ̃ε

ε ̃ = ε − 2ν







∇

√

k







2

(9)

The boundary conditions enforced on the wall for the turbulence

equations are

kw = 0, εw = 2ν





∇

√

k







2

w

(10)

Freestream levels for k and ε are derived from the ISAAC user

input variables TKEINF 12, 13 and RMUTNF ≡ μt/μ. TKEINF is then

12 It is interesting to note that if the

turbulence intensity is related to the

ratio of the turbulent kinetic energy to

the mean flow kinetic energy, I ∼

p

k/E ̄

(where k = 1

2

u′· u′, E = ̄ 1

2U·U), then,

assuming in the freestream that U2∞ =

U · U, an expression for the freestream

turbulent kinetic energy can be written

as k∞ ∼ 1

2

I

2∞U2∞ = 1

2

I

2∞M2∞a2∞.

13 More often, and as is the case for

this study, the streamwise turbulence

intensity is defined using the streamwise mean and fluctuating velocities

as I ∼

p

u′u′/U. For isotropic flow

k reduces to 3

2

u′u′, and again solving

for k and expressing the velocity in

terms of the Mach number and the

freestream speed of sound, results in

k∞ = 3

2

I

2∞M2∞a2∞. The difference

comes from the definition of I and the

method of introducing k back into the

equations. Typical levels for I∞ are approximately 10−2

for grid turbulence

and 10−1for shear flows. The input

variable TKEINF is equivalent to 3

2

I

2∞.

nondimensionalized by the square of the speed of sound and used

as the freestream value for k. RMUTNF determines the freestream

value for ε via the expression for the turbulent eddy viscosity.14 The

14 μt = Cμ

ρk2

ε

, R ≡ Reynolds number

per unit length, and M∞ ≡ freestream

Mach number.

nondimensional forms designated by the ′ are then,

k

′

∞ ≡ TKEINF × M2∞ =

3

2

I

2

∞M2∞ (11)

ε

′

∞ ≡

Cμρ∞k

′2

∞

RMUTNF 

R

M∞



=

Cμρ∞k

′2

∞

(μt/μ)∞



R

M∞



(12)

3
Algebraic Reynolds stress model

One can express the Reynolds stress transport equations in terms of

the stress anisotropy tensor b (= τ /2k − I/3),

2k

Db

Dt

=

Dτ

Dt

−

τ

k

Dk

Dt

=



P −

τ

k

Pk



+ Π +



D −

τ

k

Dk



(13)

where D represents the combined terms of viscous diffusion and

turbulent transport of the Reynolds stresses; Dk = {D}/2 similarly

denotes the viscous diffusion and turbulent transport of the kinetic

energy equations. The braces, { }, denote a tensor inner product.

Π represents the stress redistribution tensor, here modeled by the

pressure–strain relation developed by Speziale, Sarkar and Gatski

(SSG).15

15 C. Speziale, S. Sarkar, and T. Gatski.

Modeling the pressure-strain correlation of turbulence: an invariant dynamical systems approach. J. Fluid Mech.,

227:245–272, 1991.

Π = −C1εb + C2kS + C3k



bS + Sb −

2

3

{bS}I



−C4k



bW − Wb(14)

with the symmetric and antisymmetric strain rate tensors, respectively,

S =

1

2



∂ui

∂xj

+

∂uj

∂xi



, W =

1

2



∂ui

∂xj

−

∂uj

∂xi



. (15)

The coefficients of the SSG model are denoted C1 − C4:

C1 = C

1

1

Pk

ε

+ C

0

1

, C0

1 = 3.4, C11 = 1.8,

C2 = 0.36, C3 = 1.25, C4 = 0.4 (16)

Similarly, an expression for the production, P, in terms of the

tensors, b, S and W is,

P = −

4

3

kS − 2k (bS + Sb) + 2k (bW − Wb) (17)

Substituting eqs. (17) and (14) into eq. (13) the implicit anisotropy

transport equation is

2k

D

Dtb =



−2Pk − C1ε



b −

4

3

− C2



kS

−



2 − C3



k



b S + S b +

I

3k

Pk



+



2 − C4



k



bW − W b+



D −

τ

k

Dk



(18)

One can form various models for D, both differential and algebraic,

(e.g. the gradient-diffusion model of Daly-Harlow16 or algebraic

16 B. Daly and F. Harlow. Transport

Equations of Turbulence. Phys. Fluids

B, 11:2634–2649, 1970.

4
forms that model D = f(τ, b, Dk

, k)

17). The present work uses the 17 T. Gatski and C. Speziale. On

Explicit Algebraic Stress Models for

Complex Turbulent Flows. J. Fluid

Mech., 254:59–78, 1993; J.-R. Carlson, N. Duquesne, C. Rumsey, and

T. Gatski. Computation of Turbulent

Wake Flows in Variable Pressure Gradient. Computers & Fluids, 30:161–187,

2001.

algebraic model for the turbulent transport and viscous diffusion

by Gatski & Speziale and the solution of the explicit equation

by Girimaji.18 The standard algebraic model for D,

19 eq.(19), is

18 S. Girimaji. Fully Explicit and SelfConsistent Algebraic Reynolds Stress

Model. Theoret. Comp. Fluid Dynamics, 8:387–402, 1996

19 Most explicit algebraic Reynolds

stress models make this assumption

for the viscous diffusion and turbulent

transport terms, D = Dμ + Dt, of the

stress anisotropy transport equations.

A good discussion on algebraic stress

models and explicit representations can

be found in T. Gatski and T. Jongen.

Nonlinear eddy viscosity and algebraic

stress models for solving complex turbulent flows. Progress in Aerospace

Sciences, 36(8):655–682, 2000.

substituted into eq.(18)

D =

τ

k

Dk (19)

so that the trailing term is now identically zero. Next, one must

derive the expression for the production of turbulent kinetic energy

in terms of the anisotropy and symmetric strain rate tensor from

one-half the trace of eq. (17),

Pk =

1

2

{P} = −2k{bS} (20)

Substituting in the expression for Pk, eq. (20), and assuming

an equilibrium condition of Db/Dt = 0, one writes the implicit

equation for the anisotropy tensor,

0 = −



C

0

1 − 2(2 + C

1

1

){bS}T

 b

2T

− a1S

−a3



b S + S b −

2

3

{bS}I



+ a2



bW − W b(21)

a1 =

1

2

4

3

− C2



, a2 =

1

2



2 − C4



, a3 =

1

2



2 − C3



(22)

a4 =

T

γ0

Pk

ε

+ γ1

, T =

k

ε

(23)

One can represent the anisotropy tensor by an expansion of a set of

basis tensors T(n) with scalar coefficients αn as,

b =

X

3

n=1

αnT

(n)

(24)

with

T

(n) =







S , n = 1

SW − WS , n = 2

S

2 −

1

3

{S

2}I , n = 3

(25)

A basis set of three is able to fully describe two-dimensional

flows.20 Substitution of eq. (24) into eq. (21) results in a cubic

20 An evaluation of 3- vs 5-term basis

sets can be found in T. Jongen and

T. Gatski. General Explicit Algebraic

Stress Relations and Best Approximation for Three-Dimensional Flows.

IJES, 36:739–763, 1998. equation in scalar coefficient α1, i.e.,

α1

T

3

+ p

α1

T

2

+ q

α1

T



+ r = 0 (26)

5
where 21

21 The first two invariants of S and W

are usually defined as η1 = {S

2} and

η2 = {W2}, respectively. Some previous articles have used the nomenclature

of η

2

(≡ η1) and R (≡ −η2/η1).

p = −

γ1

η

∗

1

γ0

(27)

q =

1

(2η

∗

1

γ0)

2



γ

2

1 − 2η

∗

1γ0a1 −

2

3

η

∗

1a

2

3 − 2η

∗

2a

2

2



(28)

r =

γ1a1

(2η

∗

1

γ0)

2

(29)

γ0 =

1

2



C

1

1 + 2

, γ1 =

1

2



C

0

1 − 2



(30)

η

∗

1 = T2

{S

2

}, η∗

2 = T2

{W2} (31)

There are three roots to eq. (26) and can be written as,

R1 = −

p

3

+ 2r

−

a

3

cos 

θ

3



(32)

R2 = −

p

3

+ 2r

−

a

3

cos 

θ

3

+

2π

3



(33)

R3 = −

p

3

+ 2r

−

a

3

cos 

θ

3

+

4π

3



(34)

with

a = q −

p

2

3

(35)

b = 2p

3 − 9q2 + 27r (36)

d =



b

2

2

+

a

3

3

(37)

θ = arccos(φ), φ =

−b/2

p

−a

3/27

(38)

The only physically realizable solution comes from the smallest

negative root of the cubic equation. An examination of the

equations when d < 0 and whether b is either positive or negative

reveal that R2 will always be the smallest negative root using the

following logic:

If d < 0, R2 will always be the smallest negative root, therefore

C

∗

μ = −

α1

T



= −R2 (39)

If d > 0, there are two roots,

R1 = −

1

3

p + sign −

b

2

+

√

d

 



 −

b

2

+

√

d







1/3

+sign −

b

2

−

√

d

 



 −

b

2

−

√

d







1/3

(40)

6
R2 = −

1

3

p −

1

2

sign −

b

2

+

√

d

 



 −

b

2

+

√

d







1/3

−

1

2

sign −

b

2

−

√

d

 



 −

b

2

−

√

d







1/3

(41)

where the model coefficient is

C

∗

μ = −

α1

T



= −min (R1, R2) (42)

α2 = a2a4α1 (43)

α3 = −2a3a4α1 (44)

One then writes the Reynolds stress relation as,

τ =

2

3

kI + 2kα1 T

(1) + 2kα2 T(2) + 2kα3 T(3) (45)

=

2

3

kI − 2C∗

μkTS −

2C∗

μkT

2

4γ0C∗

μT2 + γ1



a2



SW − WS− 2a3



S

2 −

1

3

{S

2

}I

 (46)

In order to capture the boundary layer log-layer region, one uses the

expression C∗

μ

fμ in place of C∗

μ

in the linear stress term, making the

final expression for the Reynolds stresses

τ =

2

3

kI − 2C∗

μ

fμkTS −

2C∗

μkT

2

4γ0C∗

μT2 + γ1



a2



SW − WS− 2a3



S

2 −

1

3

{S

2

}I

 (47)

3 Results and Discussion

This section presents four configurations concerning the initial

validation of this EASM in the CFD (computational fluid dynamics)

code ISAAC.

1. Channel flow DNS at Rτ = 590.22

22 R. Moser, J. Kim, and N. Mansour.

DNS of Turbulent Channel Flow up to

Reτ = 590. Phys. Fluids, 11:943–945,

1999.

2. Zero pressure gradient flat plate flow DNS at Rθ = 1410.23

23 P. Spalart. Direct Simulation of

Turbulent Boundary Layers up to Reθ

= 1410. J. Fluid Mech., 187:61–98,

1988.

3. Backward facing step experiment.24

24 D. Driver and H. Seegmiller. Features of a Reattaching Turbulent Shear

Layer in Divergent Channel Flow.

AIAA J., 23:163–171, 1985.

4. Transonic flow over an axisymmetric afterbody experiment.

25

25 D. Reubush. Experimental Investigation to Validate Use of Cryogenic Temperatures to Achieve High Reynolds

Numbers in Boattail Pressure Testing.

NASA TM X–3396, August 1976.

This study will also compare solutions utilizing the ISAAC flow

solver of the the zero pressure gradient flat plate and the transonic

flow over the axisymmetric afterbody with solutions from using the

solver PAB3D.

7
ISAAC validation test cases

Turbulent channel flow at Rτ = 590

DNS data—Computational results from the RANS code will be

compared with the numerical experiment (DNS) of Moser et al.26

26 R. Moser, J. Kim, and N. Mansour.

DNS of Turbulent Channel Flow up to

Reτ = 590. Phys. Fluids, 11:943–945,

1999.

for a fully developed channel flow at a Reynolds number (Reτ ) of 590.

The direct numerical method used a Chebychev-tau formulation

in the wall-normal direction and a Fourier representation in the

two other directions. The DNS calculation used a third-order

Runge-Kutta time discretization for the nonlinear terms, and it

applied periodic boundary conditions in the streamwise and spanwise

directions. The computational space was dimensioned 384×257×384

with a streamwise length of 2πδh, ( δh ≡ channel half-height ).

Computational mesh for ISAAC—The grid consisted of a single

block of an H-type mesh topology with the dimensions of 201×201

with the streamwise direction in i and the normal (viscous) direction

in j, respectively. Figure 1 shows the one-quarter density grid. The

i

j

subsonic

inflow

subsonic

outflow

farfield

farfield

(1/4 grid density)

wall

wall

x=5

y=0

x=0 (leading edge of wall)

y=0.01

Figure 1: Channel flow computational grid with boundary conditions.

inflow boundary was -0.5 units upstream of the channel entrance.

The channel length was 5 units, the channel full height (i.e., from

wall to wall) was 0.01 units, and the first grid point away from

the wall was at 5.0×10−6 units (the first cell center was placed at

y

+ ≈ 0.3 ) and geometrically stretched away from the wall at a rate

of approximately 4 percent. The grid spacing at the channel leading

edge was 0.025 units and stretched downstream from the leading

edge at a rate of approximately 0.5 percent.

This study applied a ‘subsonic inflow’ boundary condition to the

inflow face i = 1, setting the total temperature and total pressure to

the input freestream conditions; additionally it applied a ‘subsonic

outflow’ boundary condition to the outflow face i = imax, with the

back pressure to the reference static pressure, p/p∞, set to 1. Since

the subsonic inflow boundary cannot immediately abut a no-slip

8
wall, the first 40 cells of the upper and lower boundaries had a

farfield (1-D Riemann invariants) boundary condition imposed.

This study also set the initial freestream turbulent eddy viscosity

ratio (RMUTNF) to 10, and turbulence intensity (TKEINF) to 1.×10−4.

The freestream Mach number is 0.20, and the Reynolds number per

unit grid length is 5.8×106in order to achieve a channel Reynolds

number, Rτ , based on the local friction velocity, uτ , of 590 in the

channel. The input unit Reynolds number was manually iterated

until the correct Rτ developed in the channel.

Grid density study—This study calculated solutions at four different

grid density levels. Each coarser grid was generated by removing

every other grid point from the previous level grid.

Figure 2(a) shows the variation with the grid density of the

channel Reynolds number versus local Reynolds number.

Rx

Rτ

106107108

0

200

400

600

800

1000

(a) Rτ vs. Rx ; —–, Baseline grid ; - - - -, 1/2 grid ; −·−·,

1/4 grid ; − · · − ··, 1/8 grid

Ncell

Rτ

103104105

300

400

500

600

700

800

Rτ= 590

(b) Rτ vs. Ncell ; #, Rx = 2.0×107; 2, Rx = 5.0×106

Figure 2: Effect of grid density.

With the exception of the region of developing flow (Rx < 7×106),

the 1/4, 1/2 and baseline grids all result in very similar Rτ

distributions. The channel Reynolds number at two streamwise

stations, Rx = 5.0 × 106 and Rx = 2.0 × 107is plotted against

the total number of computational cells, fig. 2(b). The channel

Reynolds number at each streamwise station in the region of

developing flow, Rx = 5.0 × 106, appears to have some local

grid sensitivity changing roughly 1 percent between the 1/2 and

baseline grids. The channel Reynolds number at the station further

downstream at Rx = 2.0 × 107 varies less than 0.3 percent between

the 1/2 and baseline grids.

9
y

+

u

+

10-1 100101102103

0

5

10

15

20

25

(a) Law-of-the-wall

y

+

k

+

0 200 400 600

0

1

2

3

4

(b) Turbulent kinetic energy

y

+

-τ

+12

0 200 400 600

0

0.2

0.4

0.6

0.8

1

(c) Turbulent shear stress

Figure 3: Effect of grid density at Rτ = 590: —–, Baseline grid (Rτ = 590) ; 2, 1/2 grid (Rτ = 588) ; ,

1/4 grid (Rτ = 586) ; , 1/8 grid (Rτ = 565). .

Mean flow, turbulent kinetic energy, and turbulent shear stress

are plotted with grid density in fig. 3. All of the profiles show good

grid convergence. Very little shift is observed between solutions

generated using the 1/4, 1/2 and baseline grids.

y

+

u

+

10-1 10010110210

0 3

5

10

15

20

25

(a) Mean flow (log-linear scale)

y

+

k

+

0 100 200 300 400 500 600

0

2

4

6

(b) Turbulent kinetic energy

y

+

-τ

+12

0 100 200 300 400 500 600

0

0.2

0.4

0.6

0.8

1

(c) Turbulent shear stress

Figure 4: Comparisons with DNS at Rτ = 590, linear scale. ——, EASM; #, DNS [Moser, Kim & Mansour].

EASM comparisons with DNS—Figures 4–7 show an evaluation of

the EASM for predicting mean and turbulent flow characteristics of

the channel flow. In general, there is good agreement between the

EASM and the DNS data, with the exception of the anisotropy of

the normal stresses approaching the wall and the channel centerline.

The normal Reynolds stresses are equal to 2/3 k with linear

two-equation eddy viscosity turbulence models. The pressure-strain

relation in the Reynolds stress transport equation (the basis for

the EASM) can redistribute turbulent kinetic energy between the

10
y

+

u

+

10-1 10010110210

10 3

-1

100

101

102

O(y)

(a) Mean flow

y

+

k

+

10-1 10010110210

10 3

-2

10-1

100

101

O(y2)

(b) Turbulent kinetic energy

y

+

-τ

+12

10-1 10010110210

10 3

-5

10-4

10-3

10-2

10-1

100

O(y3)

(c) Turbulent shear stress

Figure 5: Comparisons with DNS at Rτ = 590, log-linear scale. ——, EASM; #, DNS [Moser, Kim &

Mansour].

different stress components resulting in anisotropy in the flow as one

may observe in figs. 6a and 6b. Though in the region approaching

y

+

0 100 200 300 400 500 600

0

2

4

6

8

τ11

τ33

τ22

(a) Linear scale

y

+

10-1 100101102103

10-4

10-3

10-2

10-1

100

101

τ τ 11 33

τ22

O(y2)

(b) Log-log scale

Figure 6: Comparisons with DNS at Rτ = 590, turbulent normal

stresses. ——, EASM; #, DNS [Moser, Kim & Mansour].

the wall, the stress anisotropy due to the presence of the solid

boundary in not correctly modeled. The normal Reynolds stresses

predicted by the CFD begin to depart from the DNS around y+ ≈

100 and, approaching the wall, become isotropic by y+ ≈ 0.5 (see

fig. 6b).

Figure 7 shows a comparison of the predicted turbulent kinetic

energy budget with DNS. As a result of the damping of the linear

terms in the expression for the Reynolds stresses (see eq. (47) and

the location of the maximum and general profile of the production

of turbulent kinetic energy, Pk, the EASM fairly closely matches the

11
y

+

0 20 40 -0.4

-0.2

0

0.2

0.4

Pk

-ε

Tk

Dk

Figure 7: Comparisons with DNS at Rτ = 590, turbulent kinetic energy

budget. ——, EASM ; #, DNS [Moser, Kim & Mansour].

DNS calculation. As a consequence the mean velocity profile is also

well matched as seen in fig. 4a. The dissipation rate, ε, balances

the viscous diffusion of k, Dk, at the wall although the wall values

of these two quantities are elevated. This is likely due to elevated

near-wall asymptotic value of k (see fig. 5b).

The absence of a damping function in the turbulent eddy viscosity

in the k and ε equations contributes to the excessively high level of

the turbulent transport of k, Tk, compared to the DNS. Additionally

it contributes to the Tk profile peak occurring too close to the wall

(see EASM at y

+ ≈ 5, fig. 7).

Zero pressure gradient flat plate

DNS data—Spalart27 used a fully spectral method in space based 27 P. Spalart. Direct Simulation of

Turbulent Boundary Layers up to Reθ

= 1410. J. Fluid Mech., 187:61–98,

1988.

on Fourier series in the directions parallel to the plate and an

exponential mapping in the normal direction. The time integration

uses a Runge-Kutta scheme for the transport terms and a CrankNicolson scheme for the viscous terms. The following section uses

mean and turbulent flow data at the Reynolds number based on

momentum thickness of 1410 from this DNS.

Computational mesh for ISAAC—The 5-unit long, flat plate, single

block grid had an H-type mesh topology. The computational domain

included an inflow block extending 2.5 units upstream from the

leading edge of the 5 unit flat plate. The initial streamwise grid

spacing at the leading edge of the plate was 5.0×10−4 units and was

exponentially stretched from the leading edge to the trailing edge

at a rate of 4 percent with a total of 161 grid points. The first cell

12
height of the baseline grid was 1.0 × 10−6 units ( first cell center at

y

+ ≈ 0.1 ) fixed at both ends of the plate and exponentially stretched

from the surface to the outer boundary at a rate of 6 percent with a

total of 201 grid points. The upper boundary was 5 units away. The

baseline grid had the dimensions 201 × 201. Numerical transition

to turbulent flow occurred around Rx ≈ 0.5 × 106 or Rθ ≈ 500,

which corresponds to a physical distance of approximately 0.004

units downstream of the plate leading edge. Grid cell counts were

divisible by eight to allow three levels of mesh sequencing and/or

the use of a multigrid to assist solution convergence.

Ncell

cf

103104105

0

0.002

0.004

0.006

(a) cf vs. total number of computational cells

Ncell

δ

*

103104105

0

0.002

0.004

0.006

(b) δ

∗ vs. total number of computational cells

Figure 8: Effect of grid density on local skin friction and

displacement thickness at fixed

stations. #, Rθ = 1410. ; 2,

Rθ = 1.0 × 104; —–, Baseline

grid reference line.

Grid density study—Figure 8 shows grid density effects for the

zero pressure gradient flat plate. The solution functionals, local

skin friction and boundary layer displacement thickness, are plotted

against total cell count in fig. 8. These solution functionals are not

monotonically grid converged at Rθ = 1410, circle symbols in fig.

8a, because the skin friction of the 1/2 density solution is slightly

higher than both the 1/8 density and baseline grid density solutions.

The variations between the skin friction at Rθ = 1410 for these 3

finest grids are within 1 percent. One sees a similar trend in the

variation of the boundary layer displacement thickness at the same

two stations, fig. 8b, where the 1/2 density grid solution produces

results below those of both the 1/8 density and baseline density

solutions. This is due to the point of transition from laminar to

turbulent flow occurring further downstream for the 1/2 density

grid than any other grid.

Further downstream of the turbulent transition point, the

variation of local skin friction at Rθ = 1.0 × 104, square symbols

in fig. 8a, did not change for the 3 finest grids. Variations in

the turbulent transition point can be masked because the skin

friction correlation ( cf vs. Rθ ) use local functional values, that

is, the local skin friction, cf = 2τw/ρu2

∞, and a Reynolds number

based on the local momentem thickness. Changes in the transition

location globally affect the displacement thickness, hence the same

non-monotonic grid density trend occurs at both Rθ stations, see

square symbols in fig. 8b.

Figures 9 and 10 shows the variation in the solution with the grid

density for the mean flow, turbulent kinetic energy, Reynolds shear

stress at (Rθ = 1410), the skin friction correlation and boundary

layer thickness. The mean flow solution changed little with increase

in grid resolution. Slightly more variance is observed in the turbulent

kinetic energy and the turbulent shear stress at y+ ≈ 200 and above.

This could be in part due to the differing states of transition of the

boundary for each grid resolution solution.

The survey station of Rθ = 1410 was only slightly downstream of

the start of the numerical laminar-to-turbulent transition phenomena

13
y

+

u

+

10-1 10010110210

0 3

5

10

15

20

25

(a) Law-of-the-wall profiles

y

+

k

+

0 100 200 300 400 500 600

0

1

2

3

4

(b) Turbulent kinetic energy

y

+

-τ

+12

0 100 200 300 400 500 600

0

0.2

0.4

0.6

0.8

1

(c) Turbulent shear stress

Figure 9: Effect of grid density at Rθ = 1410. #, Baseline grid ; 2, 1/2 grid ; , 1/4 grid ; , 1/8 grid .

Rθ

cf

10110210310410

10 5

-4

10-3

10-2

10-1

Rθ= 1410

(a) Skin friction correlation

Rx

δ

*

103104105106107108

10-5

10-4

10-3

10-2

Rθ= 1410

(b) δ

∗ vs Rx

Figure 10: Effect of grid density at Rθ = 1410. #, Baseline grid ; 2, 1/2 grid ; , 1/4 grid ; , 1/8 grid .

14
(constrast the location of Rθ = 590 in figures 10a and 10b). Some

distance appears to be required before ‘fully’ turbulent flow is

achieved; therefore, in this transitional region the profiles of various

quantities are subject to some distortion. The difference in the state

of the boundary layer is not as apparent in fig. 10a as it is in fig.

10b, noting where Rθ = 1410 falls in each plot, and that the variance

in δ

∗

is observed. Appendix B briefly discusses various indicators of

laminar or turbulent flow.

EASM comparisons with DNS—Figures 11 through 13 shows

comparisons between the EASM and DNS for various quantities.

Here as well, the EASM closely matched the DNS for the mean

flow, turbulent kinetic energy, and the Reynolds shear stress. The

general characteristics of the normal Reynolds stresses are captured

similar to the channel flow, again with the exception of missing the

nearwall anisotropy.

y

+

u

+

10-1 10010110210

0 3

5

10

15

20

25

(a) Mean flow (log-linear scale)

y

+

k

+

0 100 200 300 400 500 600

0

2

4

6

(b) Turbulent kinetic energy (linear scale)

y

+

-τ

+12

0 100 200 300 400 500 600

0

0.2

0.4

0.6

0.8

1

(c) Turbulent shear stress (linear scale)

Figure 11: A comparison of the mean velocity profiles with DNS at Rθ = 1410. ——, EASM ; #, DNS

[Spalart]. .

The skin friction correlation, fig. 14a, has a distinct laminar

characteristic followed by an apparently abrupt transition to

turbulent like flow starting around Rθ = 400. The growth of the

displacement thickness, δ

∗

, displays a more gradual transition from

the laminar to turbulent like flow spanning the Reynolds number

range, Rx ≈ 4×105to roughly 2.5 ×106, reference fig. 14b.28 The

28 δ∗ = 1.728x/√

Rx (Blasius) ;

δ

∗ = 0.046xR

−1/5

x (turbulent); (F. M.

White. Viscous Fluid Flow. McGrawHill, Inc., 1974)

trend and level of both the skin friction correlation and the growth

of displacement thickness duplicate the classical theory well.

Backward facing step

Experimental data—Driver et al.29 obtained the data for this section

29 D. Driver and H. Seegmiller. Features of a Reattaching Turbulent Shear

Layer in Divergent Channel Flow.

AIAA J., 23:163–171, 1985

from an experimental investigation that collected flow survey data

for a rectangular channel with a backward facing step. The channel

15
y

+

u

+

10-1 10010110210

10 3

-1

100

101

102

O(y)

(a) Mean flow (log-linear scale)

y

+

k

+

10-1 10010110210

10 3

-2

10-1

100

101

O(y2)

(b) Turbulent kinetic energy (log-log scale)

y

+

-τ

+12

10-1 10010110210

10 3

-5

10-4

10-3

10-2

10-1

100

O(y3)

(c) Turbulent shear stress (log-log scale)

Figure 12: A comparison of the mean velocity profiles with DNS at Rθ = 1410. ——, EASM ; #, DNS

[Spalart]. .

y

+

0 100 200 300 400 500 600

0

2

4

6

8

τ11

τ33

τ22

(a) Turbulent normal stresses (linear scale)

y

+

10-1 10010110210

10 3

-4

10-3

10-2

10-1

100

101

τ11 τ33

τ22

O(y2)

(b) Turbulent normal stresses (log-log scale)

Figure 13: A comparison of the mean velocity profiles with DNS at

Rθ = 1410. ——, EASM ; #, DNS [Spalart]. .

16
Rθ

cf

10110210310410

10 5

-4

10-3

10-2

10-1

(a) Skin friction correlation

Rx

δ

*

10310410510610710

10 8

-5

10-4

10-3

10-2

(b) δ

∗ vs Rx

Figure 14: Comparisons with classical theory. #, EASM ;——, Blasius

; - - - -, turbulent.

consisted of a 1-m long, 10.1-cm high by 15.1-cm wide channel,

followed by a 1.27-cm high step down. The top wall of the channel

varied in angle from -2 degrees (inward) to 10 degrees (outward),

creating either an adverse or positive pressure gradient downstream

of the step. This study examines top wall angles of 0 and 6 degrees.

In units of the step height H, the upstream channel was 8 units

high, followed by a constant channel height of 9 units for the 0

degree wall configuration. The channel height was 12.2 units for

the 6 degree wall configuration at the outflow boundary 40 units

downstream of the step face. In absolute units, the reference

conditions were Uref = 44.2 m/s, Mref = 0.128, δbl = 1.9 cm, and

Rθ = 5000 obtained from data at station x/H = -4.

Computational mesh for ISAAC—The computational domain

consisted of two blocks; the upstream block modeled the region

upstream of the backward facing step; and the downstream block

included the step and region downstream of the step. The upstream

and downstream blocks had the dimensions 121×113, and 201×225

respectively, with a 1:1 block interface. The cell height of the first

cell was 1.0 × 10−6 m and was stretched away from the wall at the

rate of approximately 15 percent until reaching the centerline of the

channel. The streamwise grid spacing at the step wall was 5.0×10−4

m and was stretched away from that point at approximately 5

percent in the upstream and downstream directions.

The location of the inflow boundary was placed 1.0 m upstream

of the step with subsonic conditions applied along that boundary.

Turbulent flow developed “naturally” along the 1.0 m wall leading

up to the location of the step. An alternative would be to define

17
j

i

wall

wall

backstep (height H)

subsonic

inflow

subsonic

outflow

(1/4 grid density)

x=-80 H x=0 x=40H

inflow boundary height = 8H

outflow boundary height = 9H

Figure 15: Backward facing step grid .

a set of fixed conditions at the reference station x/H = -4 derived

from the specified experimental quantities. The backpressure along

Table 1: Flow Parameters for the Backward Facing Step

Case (p/p∞)

outf low Mref δbl,ref (cm) Rθ,ref

Reference data 0.128 1.9 5000

EASM–α = 0 deg. 1.00126 0.131 1.70 6160

EASM–α = 6 deg. 1.0065 0.121 1.87 4750

the downstream outflow boundary was varied to approximate the

stated conditions at the reference station of x/H = -4, see table 1.

Figure 15 shows the entire computational domain and fig. 16 shows

the detail of the gridding in the vicinity of the step.

(baseline grid density)

H

j

i

Figure 16: Detail of step for the

backward facing step grid.

Grid convergence—Figure 17 shows the variation of static pressure

coefficient and local skin friction with grid density for the α = 0

configuration. Similar solutions were calculated using the 1/4, 1/2,

and baseline grids, the 1/8 density grid departing significantly from

the other three grids.

Comparisons with experimental data—Figures 18 and 19 show

comparisons of the experimental data with the ISAAC solutions

using the EASM for both α = 0 and α = 6 degree configuration.

The change in the static pressure coefficient with a variation in the

wall angle is quite well captured though the actual levels are slightly

elevated. The change in the local skin friction coefficient with the

wall angle was similarly predicted though not quite matched as well.

These results are similar to what has been published previously for

the Rumsey-Gatski EASM for the α = 0 configuration.30 The earlier

30 C. Rumsey and T. Gatski. Summary

of EASM Turbulence Models in CFL3D

With Validation Test Cases. NASA

TM–2003-212431, June 2003.

results predicted the shift in the boundary layer edge velocity with

the wall angle, fig. 18c, though the absolute levels were low in the

region downstream of the step probably due to a high back pressure.

These results also predicted a change in the flow reattachment point

with the wall angle though the level was shifted upstream of the

18
x/H

Cp

-5 0 5 10 15 20 25 30 35

-0.2

0

0.2

0.4

(a) Static pressure coefficient

x/H

cf

-5 0 5 10 15 20 25 30 35

-0.002

0.000

0.002

0.004

0.006

0.008

(b) Local skin friction coefficient

Figure 17: Variation with grid density, backward facing step, α = 0. —–, Baseline grid ; - - - -, 1/2 grid ;

− · −·, 1/4 grid ; − · · − ··, 1/8 grid.

x/H

Cp

-5 0 5 10 15 20 25 30 35

-0.2

0

0.2

0.4

(a) Static pressure coefficient

x/H

cf

-5 0 5 10 15 20 25 30 35

-0.002

0.000

0.002

0.004

0.006

0.008

(b) Local skin friction coefficient

x/H

ue/uref

0 20 40

0.7

0.8

0.9

1

1.1

(c) Normalized boundary layer edge velocity

Figure 18: Variation with wall divergence angle, backward facing step. #, Data, α = 0 ; 2, Data, α = 6

[Driver and Seegmiller] ; —–, EASM, α = 0 ; - - - -, EASM, α = 6.

19
data, fig. 19. The error bar on the experimental data is from the

published uncertainty in the determination of the point of flow

reattachment. The error bar applied to the CFD data was the local

grid spacing. The local grid spacing was considered to be the finest

resolution possible for determining the reattachment point.

α

xr/H

-4 -2 0 2 4 6 8 10 12

4

6

8

10

12

Figure 19: Variation of reattachment point with wall angle,

backward facing step. •, Data

[Driver and Seegmiller] (Error

bar indicates stated experimental uncertainty) ; #, EASM

(Error bar indicates local grid

spacing) .

Axisymmetric afterbody

Test facility—This test case was an axisymmetric geometry that

was part of a series of models tested in both the Langley 1/3-m

Pilot Transonic Cryogenic Tunnel, fig. 20, and the Langley 16-Foot

Transonic Tunnel. The Pilot Tunnel had an octagonal test section

with slots at the corners of the octagon similar to the Langley 16-Foot

Transonic Tunnel test section.31 The test medium for the cryogenic

31 R. Kilgore, J. Adcock, and J. Edward. Flight Simulation Characteristics of the Langley High Reynolds

Number Cryogenic Transonic Tunnel.

AIAA Paper 74-80, January 1974.

Figure 20: Axisymmetric afterbody model in 1/3-meter Pilot Transonic

Cryogenic Tunnel .

tunnel was liquid nitrogen cooled air. High Reynolds number

data were obtained in the Pilot Tunnel through a combination of

cryogenic freestream temperatures and freestream total pressure

that were independently controllable. Approximately 5 atm of

pressure and 100K total temperature produced a unit Reynolds

number of 260 × 106/m.

This test case conducted the experiment over a range of temperatures from approximately 100K to 300K and pressures from 1

to 5 times the standard atmospheric level. Several variations of

freestream total temperatures or pressures resulted in the same

free stream Reynolds number. Surface pressure coefficients and

nozzle-boattail drag were shown to be similar, regardless of the

temperature and pressure combinations that created the equivalent

Reynolds numbers.

Previously32 high Reynolds number simulations with a CFD

32 J.-R. Carlson. Prediction of High

Reynolds Number Flow Using Algebraic Reynolds Stress Turbulence Models, Part 2: Transonic Shock-Separated

Afterbody. J. Propulsion and Power,

13:620–628, 1997.

solver were obtained through increased total pressure alone rather

than through a combination of freestream total pressure and using

20
cryogenic temperatures. The temperature range of Sutherland’s law

used for calculating bulk viscosity imposes some limits to the input

flow conditions. ISAAC input parameters are Reynolds number per

unit grid length, static temperature and Mach number, so similarly,

the high Reynolds number conditions were obtained through a

change in density rather than cyrogenic temperatures.

x

y

-40 -20 0 20 40

0

20

40

60

80

(a) 1/4 grid.

x

y

15.5 16 16.5 17 17.5

0

0.5

1

1.5

2

(b) Detail of afterbody region-baseline grid

Figure 21: Axisymmetric afterbody computational grid

Axisymmetric Afterbody Model and Conditions—The model validation geometry presented herein is one of six models that were

built for a Reynolds number investigation performed by Reubush.33

33 D. Reubush. The Effect of Reynolds

number on Boattail Drag. AIAA paper 75-63, January 1975; D. Reubush

and L. Putnam. An Experimental and

Analytical Investigation of the Effect

on Isolated Boattail Drag of Varying

Reynolds Number up to 130 × 106.

NASA TND–8210, May 1976; and

D. Reubush. Experimental Investigation to Validate Use of Cryogenic Temperatures to Achieve High Reynolds

Numbers in Boattail Pressure Testing.

NASA TM X–3396, August 1976 .

This model had a characteristic length of 16 in. and the boattail

geometry was a circular-arc, with a length-to-maximum-diameter

ratio (fineness ratio) of 0.8 boattail. The nose of the model was

a 28◦cone, 1.7956 in. long, fairing to the cylindrical body via

a 1.3615-in.-radius circular-arc whose center is 2.125 in. downstream

of the model nose and 0.8615 in. below the model centerline.

The circular-arc fairing is tangent at its endpoints to the conical

nose (1.7956 in. from the nose) and cylindrical body (2.125 in. from

the nose). The model was sting mounted with the diameter of the

sting being equal to the model base diameter. The length of the

constant diameter portion of the sting (6.70 in. measured from the

nozzle connect station) was such that there should be no effect of

the sting flare downstream of the nozzle trailing edge on the boattail

pressure distributions. These numbers will produce an analytically

consistent surface definition. The model was cast aluminum and

from conversations with the original test enginer, it was felt that

the models were possibly not manufactured to tolerances that would

21
have produced an analytically smooth shape.

The grid utilizes an H-type mesh topology, see fig. 21. The block

dimensions are divisible by 4 to facilitate investigating grid density

effects on the flow solution. The mesh is gridded with the streamwise

flow direction oriented along the i index and the j index in the

normal direction. The body is described using 100 cells extending

from the leading edge of the nose to the nozzle connect station

and 80 cells from the nozzle connect station to the nozzle-boattail

trailing edge, as well as, 80 cells extending 80 body radii from the

body surface to the far field. The boundary layer grid expansion

rate away from the body is approximately 16 percent. The inflow

boundary is 40 body radii upstream of the model nose and the

outflow boundary is 35 radii downstream of the nozzle boattail.

Solid walls are treated as no-slip adiabatic surfaces. Riemann

invarinats along characteristics are specified for boundary conditions

along the inflow and the lateral freestream outer boundaries of the

flow domain. The extrapolation boundary condition is applied on

the downstream outflow face. The flux limiter is min-mod and mesh

sequencing is used to enhance numerical convergence of the solution.

x/dm

Cp

-1 0 1 2 3

-0.4

-0.2

0

0.2

(a) R = 7×106

x/dm

Cp

-1 0 1 2 3

-0.6

-0.4

-0.2

0

0.2

0.4

(b) R = 128.5×106

Figure 22: Variation of afterbody pressure coefficient with grid density, M = 0.9, L = 16. #, baseline grid;

2, 1/2 grid, , 1/4 grid .

The freestream conditions for the axisymmetric CFD cases were

M = 0.9, T∞ = 300F, using air at γ = 1.4. The first cell height

of the grid of each configuration was different for each freestream

Reynolds number according to table 2. Input Reynolds numbers of

Table 2: Schedule of Reynolds

Number and Grid Spacing Normal to the Wall

R(×106) R(×106)/in. h1 (in.)

7.0 0.4375 6×10−5

128.3 8.0313 2×10−6

4.375 × 105/in. and 8.0313 × 106/in. were used to develop the low

and high Reynolds number solutions respectively.

22
Figure 22 shows grid sensitivity of the EASM at M = 0.9 at

the lowest and highest Reynolds number for this test case, (a) R

= 7 × 106 and (b) R = 128 × 106respectively. These sensitivities

were relatively consistent for the other turbulence models and

other viscous models investigated.34 The EASM prediction of

34 J.-R. Carlson. Prediction of High

Reynolds Number Flow Using Algebraic Reynolds Stress Turbulence Models, Part 2: Transonic Shock-Separated

Afterbody. J. Propulsion and Power,

13:620–628, 1997. the afterbody shock strength and pressure recovery was Reynolds

number dependent, see fig. 23.

x/dm

Cp

-1 -0.5 0 0.5

-0.4

-0.2

0

0.2

(a) R = 7×106

x/dm

Cp

-1 -0.5 0 0.5

-0.6

-0.4

-0.2

0

0.2

0.4

(b) R = 128.5×106

Figure 23: Comparison of afterbody pressure coefficients with experimental data at M = 0.9, L = 16. ——,

EASM ; #, φ = 0; 2, φ = 120; ⋄, φ = 240 Data [Reubush] .

N–Version Testing

Flat plate

Figures 24 through 30 compare the solutions produced by the

ISAAC code with those produced by the PAB3D code for the zero

pressure gradient flat plate. Several flow variables, such as mean

flow and turbulence quantities, are compared at two stations, Rθ =

1410 and Rθ ≈ 10000, followed by comparisons of the development

of boundary layer displacement thickness, shape factor and the skin

friction correlation.

A very good match has been obtained between the two solvers for

the mean flow at station Rθ = 1410, see fig. 24, but the turbulent

kinetic energy and Reynolds shear stress at station Rθ = 1410 are

not as close, see figs. 25 , 26 and 27. Further downstream of the

numerical transients caused by the transition to turbulent flow at

the station Rθ ≈10000, the two codes appear to match more closely,

see figs. 28 and 29.

The first boundary layer shape factor, H12 and the skin friction

23
y

+

u

+

10-1 10010110210

0 3

5

10

15

20

25

(a) Mean flow, Log-linear scale

y

+

u

+

10110210

10 3

15

20

25

(b) Mean flow, Detail of log-layer

Figure 24: Zero pressure gradient flat plate at Rθ = 1410. •, ISAAC; #, PAB3D .

y

+

k

+

0 100 200 300 400 500 600

0

1

2

3

4

(a) Turbulent kinetic energy.

y

+

-τ

+12

0 100 200 300 400 500 600

0

0.2

0.4

0.6

0.8

1

(b) Reynolds shear stress.

Figure 25: Zero pressure gradient flat plate at Rθ = 1410, linear scale. •, ISAAC; #, PAB3D .

24
y

+

k

+

10-1 10010110210

10 3

-3

10-2

10-1

100

101

O(y2)

(a) Turbulent kinetic energy.

y

+

-τ

+12

10-1 10010110210

10 3

-5

10-4

10-3

10-2

10-1

100

O(y3)

(b) Reynolds shear stress.

Figure 26: Zero pressure gradient flat plate at Rθ = 1410, log-log scale. •, ISAAC; #, PAB3D .

y

+

-α1/T

10-1 100101102103

0

0.2

0.4

0.6

0.8

1

(a) −α1/T

y

+

P

k

/

ε

10-1 100101102103

0

0.5

1

1.5

2

(b) Pk/ε.

Figure 27: Zero pressure gradient flat plate at Rθ = 1410. •, ISAAC; #, PAB3D .

25
y

+

u

+

10-1 100101102103104

0

10

20

30

(a) Log-linear scale.

y

+

u

+

10110210310

10 4

15

20

25

30

(b) Detail of log-layer.

Figure 28: Zero pressure gradient flat plate at Rθ = 104, mean flow. •, ISAAC; #, PAB3D ; · · · , log-layer.

y

+

k

+

0 1000 2000 3000 4000

0

1

2

3

4

(a) Linear scale

y

+

k

+

10-1 10010110210310

10 4

-3

10-2

10-1

100

101

O(y2)

(b) Log-log scale

Figure 29: Zero pressure gradient flat plate at Rθ = 104, turbulent kinetic energy. •, ISAAC; #, PAB3D ;

· · · , log-layer.

26
Rx

δ

*

103104105106107108

10-5

10-4

10-3

10-2

(every 3rd point plotted)

transition region

(a) Momentum thickness with local Reynolds

number

Rθ

H12

101102103104105

1

1.5

2

2.5

3

turbulent

laminar

(b) First shape factor

Rθ

cf

10110210310410

10 5

-4

10-3

10-2

10-1

(c) Skin friction correlation

Figure 30: Zero pressure gradient flat plate. •, ISAAC; #, PAB3D ; —–, Linear theory-laminar ; - - - -, Linear

theory-turbulent.

correlation, cf vs Rθ, are also quite closely matched with only

minor variances seen in the upstream laminar portion of the flat

plate, see figs. 30b and 30c. Similarly, the development of the

the displacement thickness with Reynolds number, fig. 30a, for

the most part is identical between ISAAC and PAB3D. A slight

departure occurs in the the laminar-to-turbulent transition region,

3 × 105 < Rx < 2 × 106, but is relatively local and does not appear

to affect the downstream match significantly.

Axisymmetric afterbody

Using the two codes, fig. 31 compares the static pressure coefficients

on the afterbody. A slightly higher pressure occurs in the ISAAC

solution compared with the PAB3D solution for the lower Reynolds

number case, see fig. 31a just downstream of x/dm = 0.0. The

high Reynolds number comparison show virtually identical static

pressure distributions on this portion of the model.

Initially, the solutions developed using ISAAC had a considerably

weaker shock because of the type of limiter being applied to

the inviscid fluxes. Since the PAB3D solutions were developed

using the min-mod limiter, changing the limiter being used in the

ISAAC calculations from Venkatakrishnan’s smoothing limiter to

the min-mod resolved this difference in the solutions between the

two solvers.

4 Summary

The turbulence equations being solved in the three-dimensional

Reynolds averaged Navier-Stokes (RANS) code PAB3D have

27
x/dm

Cp

-1 -0.5 0 0.5

-0.4

-0.2

0

0.2

(a) R = 7.0×106

x/dm

Cp

-1 -0.5 0 0.5

-0.6

-0.4

-0.2

0

0.2

0.4

(b) R = 128.5×106

Figure 31: Comparison of static pressure coefficients, axisymmetric afterbody, M = 0.9, L = 16. —–, ISAAC

; - - - -, PAB3D ; #, φ = 0; 2, φ = 120; ⋄, φ = 240, Data [Reubush].

been implemented in the three-dimensional RANS code ISAAC.

Initial validation of the implementation was accomplished through

comparisons with data for four test cases: channel flow, zero pressure

gradient flat plate, backward facing step and an axisymmetric

transonic afterbody flow. These test cases covered a wide range of

fluid dynamic situations, including compressible flow, separated flow,

and high Reynolds number flow. The turbulence equations were

able to simulate all of these flow situations fairly well. Additionally,

excellent agreement was achieved in comparisons between solutions

developed using ISAAC and with solutions developed using the

solver PAB3D.

28
A Nomenclature

Symbol Definition

b, bij Reynolds stress anisotropy tensor

C

∗

μ Variable model coefficient, −α1/T

Cp Pressure coefficient, (p - p∞) / q∞

cf Local skin friction coefficient, τw / q∞

dm Body maximum diameter (1.0 in.) [axisymmetric afterbody]

H Backward facing step height

H12 First boundary layer shape factor, (δ

∗/θ)

h, δh Channel half-height

h1 Physical height of first computational grid from a wall

L Reference length [axisymmetric afterbody]

M Mach number

Ncell Total number of computational cells

P Production, −τik(∂uj/∂xk) − τjk(∂ui/∂xk)

Pk Production of turbulent kinetic energy, − 1

2

{P} = τij (∂ui/∂xj )

p Static pressure

q Dynamic pressure

R Reynolds number based on model reference length, (u∞L/ν)

Rt Cell turbulent Reynolds number, k

2/(νε)

Rx Reynolds number based on distance x, (u∞x/ν)

Rθ Reynolds number based on momentum thickness, (u∞θ/ν)

Rτ Reynolds number based on wall shear stress, (uτ h/ν)

S, Sij Symmetric velocity gradient tensor, (∂ui/∂xj + ∂uj/∂xi) /2

T Kinematic time scale, k/ε

u, ui Velocity vector

u+ Law-of-the-wall coordinate, u/uτ

uτ Local friction velocity, p

τw/ρw

W, Wij Antisymmetric velocity gradient tensor, (∂ui/∂xj − ∂uj/∂xi) /2

x Streamwise distance

x, xi Position vector

xr Point of flow reattachment

y Vertical distance

y+ Law-of-the-wall coordinate, y uτ / ν

α Wall divergence angle (backward facing step)

αn scalar coefficient

δ Boundary-layer thickness, value of y at 0.995 umax

δ

∗ Displacement thickness, R δ

0

(1-u/ue)dy

η

∗

1

Flow invariant, T2{S

2}

η

∗

2

Flow invariant, T2{W2}

θ Momentum thickness, R δ

0

(u/ue) (1 - u/ue) dy

μ Laminar viscosity

μt Turbulent eddy viscosity, Cμ0 ρk2/ε

ν Kinematic viscosity, μ / ρ

ρ Density

τ, τij Reynolds stress tensor

τw Shear stress at the wall, ρw du/dn

φ Angular location of pressure orifices, deg. (axisymmetric afterbody)

bl Boundary layer

ref Reference conditions

w Condition at wall surface

∞ Freestream condition

29
B A Comparison of Transition Indicators

There are several parameters that can indicate the presence of

turbulent flow, or as the case may be, developing turbulent flow over

the flat plate. These indicators occur both at a primative level, e.g.,

peak turbulent kinetic energy and dissipation rate level at the wall;

or at a derived level, e.g., local skin friction, the rate of boundary

layer growth with local Reynolds number, or boundary layer shape

factors, to name a few examples.

Typically, a variation of the local skin friction with the local

Reynolds number, fig. 32 is plotted, which clearly deliniates the

Rx

cf

103104105106107108

10-4

10-3

10-2

Rx≈400000

Rθ≈410

Figure 32: Local skin friction

with local Reynolds number. ◦,

EASM.

change from laminar flow by the rather abrupt discontinuity in the

curve around Rx = 400,000.

If the Reynolds number based on the momentum thickness of the

boundary layer is used instead of the local Reynolds number, an

equally abrupt change in the skin friction correlation from laminar

flow to turbulent flow levls is observed, fig. 33. At first glance,

Rθ

cf

101102103104105

10-4

10-3

10-2

Rθ≈410

Figure 33: Skin friction correlation. ◦, EASM ;- - - -, KarmanShoenherr empirical curve.

the local skin friction in these two figures appears to attain a fully

turbulent level almost immediately downstream of the transition

point. (As noted on the figure, the momentum Reynolds number of

410 is approximately coincident with the local Reynolds number of

400,000 and is used as a reference point for the discussion.)

Rx

δ

*

10410510610

10 7

-5

10-4

10-3

10-2

Rθ≈410

Figure 34: Displacement thickness with local Reynolds number (detail).

◦, EASM ;——, Blasius ;- - - -, turbulent (flat plate theory) .

A slightly different conclusion could be drawn (concerning the

transition of the flow from laminar to turbulent flow) from a plot of

the rate of growth in the boundary layer displacement thickness with

local Reynolds number, fig. 34. The displacement thickness only

gradually approaches the turbulent level as the Reynolds number

increases downstream of the transition point marked (Rθ ≈ 410).35

35 M. Rai and P. Moin. Direct Numerical Simulation of Transition and

Turbulence in a Spatially Evolving

Boundary Layer. Journal of Computational Physics, 109:169–192, 1993.

Similarly the first shape factor, H12 = δ

∗/θ, shows a clear

departure from laminar flow, fig. 35, but a rather prolonged region

30
Rθ

H12

10110210310410

1 5

1.5

2

2.5

3

Rθ≈410

laminar (H12=2.6)

turbulent (H12=1.27)

Figure 35: First shape factor with momentum thickness Reynolds

number. ◦, EASM .

of developing turbulent flow before the shape factor approaches the

classically derived level.

If the “simpler” parameter of peak turbulent kinetic energy is

plotted against Reynolds number, a less definitive picture appears

as to the point of transition, fig. 36. For reference, the figures show

the locations of Rx ≈ 400000 and Rθ ≈ 410.

Though still somewhat abrupt at transition, the peak turbulent

kinetic energy has a region of gradual increase in intensity before

increasing another two orders of magnitude greater than the level in

the region of laminar flow, in the current example, starting from a

value around 10−6.

Rx

k

max

10310410510610710

10 8

-6

10-5

10-4

10-3

Rx≈400000

(a) kpeak vs. Rx.

Rθ

k

max

10110210310410

10 5

-6

10-5

10-4

10-3

Rθ≈410

(b) kpeak vs. Rθ.

Figure 36: Streamwise variation of peak turbulent kinetic energy.

31
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704–0188

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection

of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports

(0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be

subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18

1. REPORT DATE (DD-MM-YYYY)

01-12-2005

2. REPORT TYPE

Technical Memorandum

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Assessment of an Explicit Algebraic Reynolds Stress Model

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

26-065-30-35

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

23-065-30-35

6. AUTHOR(S)

Carlson, Jan-Rene ́e

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

NASA Langley Research Center

Hampton, VA 23681-2199

8. PERFORMING ORGANIZATION

REPORT NUMBER

L–19103

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

National Aeronautics and Space Administration

Washington, DC 20546-0001

10. SPONSOR/MONITOR’S ACRONYM(S)

NASA

11. SPONSOR/MONITOR’S REPORT

NUMBER(S)

NASA/TM–2005–213771

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

Subject Category 01

Availability: NASA CASI (301) 621-0390

13. SUPPLEMENTARY NOTES

An electronic version can be found at http://ntrs.nasa.gov

14. ABSTRACT

This study assesses an explicit algebraic Reynolds stress turbulence model in the in the three-dimensional Reynolds averaged

Navier-Stokes (RANS) solver, ISAAC (Integrated Solution Algorithm for Arbitrary Con gurations). Additionally, it compares solutions

for two select con gurations between ISAAC and the RANS solver PAB3D. This study compares with either direct numerical simulation

data, experimental data, or empirical models for several di erent geometries with compressible, separated, and high Reynolds number

ows. In general, the turbulence model matched data or followed experimental trends well, and for the selected con gurations, the

computational results of ISAAC closely matched those of PAB3D using the same turbulence model.

15. SUBJECT TERMS

CFD, turbulence modeling, code assessment, Reynolds stress, explicit algebraic stress model

16. SECURITY CLASSIFICATION OF:

a. REPORT

U

b. ABSTRACT

U

c. THIS PAGE

U

17. LIMITATION OF

ABSTRACT

UU

18. NUMBER

OF

PAGES

36

19a. NAME OF RESPONSIBLE PERSON

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

19b. TELEPHONE NUMBER (Include area code)

(301) 621-0390