TableOfContents Turbulence Guide V4.01
User Manual:
Open the PDF directly: View PDF .
Page Count: 214 [warning: Documents this large are best viewed by clicking the View PDF Link!]
Turbulence Models and Their
Application to Complex Flows
R. H. Nichols
University of Alabama at
Birmingham
Revision 4.01
CONTENTS Page
1.0 Introduction
1.1 An Introduction to Turbulent Flow 1-1
1.2 Transition to Turbulent Flow 1-4
1.3 Statistical Concepts for Turbulent Flow Analysis 1-9
1.4 Turbulent Length and Time Scales 1-13
1.5 Vorticity and Strain Tensors 1-14
1.6 Classification of Turbulence Models 1-15
References 1-17
2.0 Classical Turbulence Closure Methodologies
2.1 Reynolds Averaging 2-1
2.2 Favre Averaging 2-4
2.3 Boussinesq Approximation 2-5
References 2-7
3.0 Turbulent Boundary Layer and Shear Layer Theory
3.1 Boundary Layer Theory 3-1
3.2 Shear Layer Theory 3-7
References 3-10
4.0 Algebraic Turbulence Models
4.1 Inner Eddy Viscosity Model 4-1
4.2 Cebeci-Smith Model 4-1
4.3 Baldwin-Lomax Model 4-2
4.4 Wake and Jet Model 4-4
4.5 Algebraic Model Shortcomings 4-4
4.6 Grid Sensitivity for a Flat Plate with Adiabatic Walls 4-6
4.7 Grid Sensitivity for an Axisymmetric Bump 4-10
4.8 Grid Sensitivity for a Flat Plate with Heat Transfer 4-14
4.9 Grid Sensitivity for a Nozzle with Heat Transfer 4-18
4.10 Summary 4-21
Baldwin-Lomax Application Hints 4-22
References 4-22
5.0 One-Equation Turbulence Models
5.1 Theory 5-1
5.2 Spalart-Allmaras Model 5-1
5.3 Rotation/Streamline Curvature Corrections 5-3
5.4 Grid Sensitivity for a Flat Plate with Adiabatic Walls 5-6
5.5 Grid Sensitivity for an Axisymmetric Bump 5-10
5.6 Grid Sensitivity for a Flat Plate with Heat Transfer 5-13
5.7 Grid Sensitivity for a Nozzle with Heat Transfer 5-17
5.8 Summary 5-20
Spalart-Allmaras Application Hints 5-20
References 5-20
6.0 Two-Equation Turbulence Models
6.1 Theory 6-1
6.2 Traditional k-
ε
Models 6-3
6.3 k-
ω
Models 6-6
6.4 SST Model 6-7
6.5 RNG Model 6-9
6.6 Numerical Implementation 6-10
6.7 Compressibility Correction for Shear Layers 6-10
6.8 Initializing Turbulence Values for a Given Profile Shape 6-13
6.9 Rotation and Curvature Correction 6-15
6.10 Grid Sensitivity for a Flat Plate with Adiabatic Walls 6-16
6.11 Grid Sensitivity for an Axisymmetric Bump 6-19
6.12 Grid Sensitivity for a Flat Plate with Heat Transfer 6-22
6.13 Grid Sensitivity for a Nozzle with Heat Transfer 6-26
6.14 Summary 6-29
Two-Equation Model Application Hints 6-30
References 6-30
7.0 Reynolds and Algebraic Stress Models
7.1 Reynolds Stress Models 7-1
7.2 Algebraic Stress Models 7-5
7.3 Grid Sensitivity for a Flat Plate with Adiabatic Walls 7-7
7.4 Grid Sensitivity for an Axisymmetric Bump 7-12
RSM and ASM Model Application Hints 7-18
References 7-18
8.0 Large Eddy Simulation
8.1 The Filtering Operation 8-1
8.2 Derivation of the LES Equations 8-4
8.3 Smagorinski Model 8-7
8.4 Dynamic Smagorinski Model 8-8
8.5 k-Equation Model 8-8
8.6 Inflow Turbulence Boundary Condition 8-9
8.7 Other LES References 8-10
8.8 Spatial Mixing Layer Example 8-10
References 8-13
9.0 Hybrid RANS/LES Turbulence Models
9.1 Theory 9-1
9.2 Circular Cylinder 9-4
9.3 WICS Bay 9-13
9.4 Delayed Detached Eddy Simulation (DDES) 9-25
9.5 Summary 9-30
Hybrid RANS/LES Application Hints 9-31
References 9-31
10.0 Wall Function Boundary Conditions
10.1 Theory 10-1
10.2 Grid Sensitivity for a Flat Plate with Adiabatic Walls 10-7
10.3 Grid Sensitivity for an Axisymmetric Bump 10-10
10.4 Grid Sensitivity for a Flat Plate with Heat Transfer 10-12
10.5 Grid Sensitivity for a Nozzle with Heat Transfer 10-17
Wall Function Application Hints 10-20
References 10-20
11.0 Boundary Layer Transition Simulation
11.1 Transition Models Based on Stability Theory 11-2
11.2 Transition Models with Specified Transition Onset 11-2
11.3 Transition Models with Onset Prediction Capability 11-6
11.4 Transition Models with Onset Prediction Capability 11-10
References 11-14
![](asset-5.png)
1-1
1.0 Introduction
Some basic knowledge of turbulence and an understanding of how turbulence
models are developed can help provide insight into choosing and applying these
models to obtain reasonable engineering simulations of turbulent flows. This
effort is directed at production users of Computational Fluid Dynamics (CFD) and
attempts to provide the basic information required to choose and use currently
available turbulence modeling techniques. The limitations of each modeling
technique and application tips are also provided.
It is believed that the Navier-Stokes equations can be used to fully describe
turbulent flows, but current limitations in computational horsepower have made
the direct solution of the Navier-Stokes equations impractical for all but very
simple flows at low Reynolds numbers. This is because current computers do
not allow for the resolution of the wide range of length and time scales
associated with turbulence. Many complex fluid dynamic applications are
directed at determining time-averaged quantities, and hence it is desirable to find
a means to obtain these mean quantities short of solving the full unsteady
Navier-Stokes equations for all of the length and time scales associated with the
turbulence.
The quest for the ultimate turbulence model has been ongoing for nearly a
century now. Early turbulence models were empirically derived algebraic
relations. As computers developed and numerical simulation evolved differential
equation based transport type turbulence models became the turbulence
simulation methodology of choice. It should always be remembered that
transport models are empirically calibrated. The use of transport type turbulence
models has become standard practice for most engineering applications. Many
current researchers are now solving the unsteady Navier-Stokes equations for
large-scale, or grid realized, turbulence and modeling the smaller, or subgrid,
turbulent scales that cannot be captured on the computational grid. Although the
ultimate general-purpose turbulence model has yet to be developed, turbulence
modeling has matured to the point that reasonably accurate results can be
obtained for a wide range of engineering applications with the current class of
computers. As computer technology continues to improve both the role and the
form of turbulence models will continue to evolve.
1.1 An introduction to Turbulent Flow
An understanding of what constitutes turbulent flow is required before proceeding
to discuss turbulence modeling. Turbulence can be parameterized by several
nondimensional quantities. The most often used is Reynolds number. Reynolds
number represents the ratio of inertial forces to viscous forces. The viscous
forces dominate at low Reynolds numbers and disturbances are damped rapidly.
These disturbances begin to amplify as Reynolds number is increased and
eventually transition into fully turbulent flows. Launder1 gives the following
![](asset-6.png)
1-2
definition for turbulent flow: “At moderate Reynolds numbers the restraining
effects of viscosity are too weak to prevent small, random disturbances in a
shear flow from amplifying. The disturbances grow, become non-linear and
interact with neighboring disturbances. This mutual interaction leads to a
tangling of vorticity filaments. Eventually the flow reaches a chaotic, non-
repeating form describable only in statistical terms. This is turbulent flow.”
Turbulence is a three-dimensional unsteady viscous phenomenon that occurs at
high Reynolds number. Turbulence is not a fluid property, but is a property of the
flow itself. Turbulent flow can be highly nonlinear and is random in nature.
Turbulent disturbances can be thought of as a series of three-dimensional eddies
of different sizes that are in constant interaction with each other. This model of
turbulence is Lagrangian in nature, with these turbulent flow structures being
transported downstream by the mean flow. These structures exist for a limited
amount of time before they are dissipated away by molecular viscosity. An
example of the fluid structures that appear in a turbulent flow can be seen in the
wake of a circular cylinder in Fig. 1.1a and the boundary layer in Fig. 1.1b.
Figure 1.1a Turbulent structures in the wake of a circular cylinder.
![](asset-7.png)
1-3
Figure 1.1b Turbulent structures in a boundary layer.
Turbulence is diffusive in nature. The rates of mass and momentum transfer are
much higher in a turbulent flow than in a laminar flow. This diffusiveness can be
a highly desirable property for enhancing the mixing of physical and/or chemical
properties within a flow. The diffusive nature of turbulence also causes boundary
layers and shear layers to become much thicker than their laminar counterparts,
which is often an undesirable result of turbulent flow. A generic turbulent kinetic
energy spectrum that indicates how the turbulent energy is partitioned among the
various size eddies is shown in Fig 1.2.
The wave number is inversely proportional to the turbulent length scale. It can
be seen in Fig 1.2 that most of the turbulent kinetic energy is in the large
turbulent length scales (or low wave numbers). The large scale eddies take
kinetic energy from the mean flow in the energy production region of the turbulent
spectrum. Smaller eddies feed off of the larger energy producing eddies in the
inertial range. Here the turbulence is essentially in equilibrium and the transfer of
Log(Wave number)
Log(Energy)
Energy
Production Inertial
Range Dissipation
Figure 1.2 Turbulent energy spectrum.
![](asset-8.png)
1-4
energy by inertial forces is the dominant process. At the small-scale end of the
turbulent spectrum the eddies become so small that viscous dissipation converts
the kinetic energy into heat.
As mentioned earlier, it is generally accepted that the Navier-Stokes equations
can be used to directly simulate turbulent flows. This type of fluid modeling is
called Direct Numerical Simulation (DNS). Direct numerical simulation of the
turbulent energy cascade using the Navier-Stokes equations requires that the
grid spacing be smaller than the smallest turbulent length scale. For a typical
full-scale aircraft, the smallest length scale of the turbulence may be of the order
of 10-6 times smaller than the aircraft reference length. More than 1018 grid
points would be required to discretize this example. Since engineering problems
of this size are impractical at this time, some model of the turbulence is required
so that engineering answers may be obtained for large-scale fluid problems.
For many engineering applications the length and time scales of the turbulence
are much smaller than the length and time scales of the problem of interest. This
tends to simplify the modeling of turbulence and remove the need to treat the
unsteady aspects of the turbulent flow since they occur at time scales much
smaller than those of interest. When the turbulent and problem scales become
of the same size the modeling becomes more difficult.
1.2 Transition to Turbulent Flow
Laminar-turbulent transition is an extraordinarily complicated process. The
following is a quick summary of the process. More detail can be found in
Schlichting2. The process begins by transforming external disturbances into
internal instability oscillations in the boundary layer or shear layer. A laminar
boundary layer or shear layer is susceptible to disturbances from both the free
stream and the body surface. The ability of external disturbances to penetrate
the boundary layer or shear layer and then be amplified is defined as the
receptivity of the boundary layer. The free stream disturbances may include
acoustic waves, particles in the flow, transported vorticity from upstream in the
flow, and pressure, density, or temperature fluctuations. Surface disturbances
include the roughness of the body or motion of the body. Many of these
disturbances eventually get damped by viscosity. The transition process starts
when a disturbance in the boundary layer or shear layer is no longer damped but
gets amplified. This results in the formation Tollmien-Schlichting waves which
are the first mode instability of the flow. While these waves travel downstream,
three dimensional waves and vortices begin to develop. At certain points in the
boundary layer, small irregularly shaped turbulent spots will occur and will be
convected in a wedge shaped region as shown in the Figure.1.3. These spots,
known as Emmons spots, appear at random locations on the plates at irregular
time intervals. As these spots move downstream, they grow and eventually fuse
with each other to encompass the entire boundary layer or shear layer. This
process eventually leads to fully developed turbulent flow.
![](asset-9.png)
1-5
Figure 1.3 Growth of turbulent spots in a flat plate boundary layer.
The transition mechanism mentioned above is known as natural transition.
Flows with strong free stream disturbances can “by-pass” the mechanism of
formation of Tollmien-Schlichting waves and form turbulent spots directly from
the influence of the free stream disturbances. This is called as by-pass
turbulence3. For compressible flows, most notably flows with a Mach number of
2.2 or greater, Mack4,5 showed that multiple modes of instability exist. For flows
between M=2.2 and M=4.5, the first mode of the instability, known as the viscous
instability, is the most unstable. As the Mach number increases beyond 4.5, the
second mode, known as the inviscid instability, becomes the most unstable
mode. This second mode includes acoustical disturbances that are characterized
by very large fluctuations in pressure and temperature.
Boundary layer stability theory seeks to predict the receptivity of the boundary
layer in terms the critical Reynolds number for a flow. The critical Reynolds
number is defined as the Reynolds number at which the disturbances in a
laminar boundary layer begin to amplify and the flow becomes turbulent.
According to this theory, the fluid motion is decomposed into a time-averaged
mean flow component and a fluctuating component.
uUU ′
+= (1.1)
Here the overbar signifies a time averaged quantity and the prime is used for
perturbation quantities. If the quantities in the Navier-Stokes equations are
decomposed as in Eq. 1.1 and the mean flow quantities are subtracted out, the
two-dimensional incompressible Navier-Stokes equations become
Laminar Transitional Turbulent
![](asset-a.png)
1-6
v
y
p
x
v
U
t
v
u
x
p
dy
Ud
v
x
u
U
t
u
y
v
x
u
′
∇=
∂
′
∂
+
∂
′
∂
+
∂
′
∂
′
∇=
∂
′
∂
+
′
+
∂
′
∂
+
∂
′
∂
∂
′
∂
+
∂
′
∂
2
2
1
)(
1
ν
ρ
ν
ρ
(1.2)
Eq. 1.2 has three equations and three unknown disturbances ,, vu ′′ and p′. If a
periodic disturbance is assumed then the solution can be expressed as a Fourier
series. This solution is linear, so solutions may be summed up. Note this is only
valid for very small amplitudes for the disturbances. For two-dimensional flow,
the stream-function can be defined as
Ψ(x,y,t) = Φ(y)ei(
α
x-
β
t) (1.3)
Where
λ
π
α
2
= is the wave number and
λ
is the wavelength of the disturbance
and i
ir
β
β
β
+= , where f
r
π
β
2= ( f is the disturbance frequency) and i
i
β
is
the amplification factor for the disturbance. The ratio of
β
/
α
is
icc
i
cir
ir +=
+
==
α
β
β
α
β
(1.4)
Here r
c is the propagation velocity in the x-direction and is known as the phase
velocity. The imaginary part, i
c, is the degree of damping if i
c is positive and is
the degree of amplification if i
c is negative.
Substituting the expressions from Eq. 1.3 into Eq. 1.2
)(
)(
)(
)(
txi
txi
eyi
x
v
ey
y
u
βα
βα
φα
φ
φ
φ
−
−
′
−=
∂
∂
=
′
′
=
∂
∂
=
′
(1.5)
)2(
Re
))(( 422
φαφαφ
α
φφαφ
+
′′
−
′′′′
−=
′′
−−
′′
−i
UcU (1.6)
Eq. 1.6 is called the Orr-Sommerfeld equation. The equation is shown in
nondimensional form. The terms on the left side represent the inertia terms. The
terms on the right side represent the viscous terms from the Navier-Stokes
![](asset-b.png)
1-7
equations. Here the Reynolds number (Re) is define as Re=
ν
δ
U where
δ
is the
boundary layer thickness.
Solving the Orr-Sommerfeld equations first requires a steady state fully laminar
solution of the Navier-Stokes equations. The pressure and velocity from the
laminar solution are used as inputs to the Orr-Sommerfeld equations. An
example solution of the Orr-Somerfield is shown in Fig. 1.4. The vertical axis is
the ratio of the shear layer thickness to the wavelength of the disturbance. The
horizontal axis is the Reynolds number based on the shear layer thickness.
Figure 1.4 Example solution to the Orr-Sommerfield equation.
The curve shown in Fig. 1.4 is the neutral stability curve. Disturbances are
amplified in the crosshatched region, while they are damped in the region outside
the curve. The Reynolds number at which disturbances are first amplified is
called the critical Reynolds number (Rcrit). For Reynolds numbers below Rcrit,
disturbances are damped. As the Reynolds number increases above Rcrit,
almost all disturbances are amplified and turbulent flow occurs.
The effect of pressure gradient is shown in Fig. 1.5. The “0” curve represents no
pressure gradient, the “-5” curve represents an adverse pressure gradient, and
the “+4” curve represents a favorable pressure gradient. An adverse pressure
gradient causes transition to occur more rapidly, while a favorable pressure
gradient delays the onset of transition.
Re
![](asset-c.png)
1-8
Figure 1.5 Effect of pressure gradient on boundary layer stability.
Transition locations in boundary layer flows are generally characterized by two
parameters: Rex and Re
θ
. Rex is the Reynolds number based on distance from
the leading edge. Re
θ
is the Reynolds number based on the momentum
thickness defined as
dy
u
u
u
u
y
yeee
∫
=
=⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−=
δ
ρ
ρ
θ
0
1 (1.7)
Here
δ
is the boundary layer height and the subscript e denotes the edge of the
boundary layer. Fig. 1.6 shows the transition process on the skin friction on a flat
plate. The skin friction changes by almost an order of magnitude during the
transition process.
Re
![](asset-d.png)
1-9
Figure 1.6 Skin friction along a flat plate.
1.3 Statistical Concepts for Turbulent Flow Analysis
Several statistical moments are used to describe turbulent flows. This section
will review those most often associated with modeling turbulence. These
statistics are usually performed on ensembles of discrete observations in both
experimental and computational analysis of turbulent flows. The first moment is
the mean or average and is defined as
∑
=
=n
ii
u
n
u
1
1 (1.8)
Here the average is taken of n discrete observations. The higher moments are
defined relative to the mean. The second moment, or variance, is defined as
()
2
1
22 1∑
=
−=
′
=n
iiuu
n
u
σ
(1.9)
The variance provides information on the spread of the data away from the
mean. The third moment, or skewness, is given by
()
3
1
33
31∑
=
−=
′
=n
iiuu
n
u
S
σσ
(1.10)
![](asset-e.png)
1-10
Skewness gives the amount of time the signal is above or below the mean. The
fourth moment, or kurtosis, is defined by
()
4
1
44
41∑
=
−=
′
=n
iiuu
n
u
K
σσ
(1.11)
Kurtosis gives the amount of time the signal is away from the mean. A simple
example of these moments can be constructed using sine waves. Four cases
are defined as
⎟
⎠
⎞
⎜
⎝
⎛
1000
2
sin:1 t
Case
π
⎟
⎠
⎞
⎜
⎝
⎛
1000
2
sin:2 2t
Case
π
(1.12)
⎟
⎠
⎞
⎜
⎝
⎛
+
⎟
⎠
⎞
⎜
⎝
⎛
000,10
2
sin
1000
2
sin:3 tt
Case
ππ
⎟
⎠
⎞
⎜
⎝
⎛
+
⎟
⎠
⎞
⎜
⎝
⎛
1000
2
sin
1000
2
sin:4 2tt
Case
ππ
The four signals are shown in Fig. 1.7.
Figure 1.7 Example signals for statistical analysis.
The statistical moments for these four cases are shown in Table 1.1.
![](asset-f.png)
1-11
Case Mean
σ
2 S K
1 0.0 0.5 0.0 1.5
2 0.5 0.125 0.0 1.5
3 0.0 1.0 0.0 2.25
4 0.5 0.625 0.759 1.98
Table 1.1 Statistical moments for four test signals.
Signals that are symmetric about the mean such as Case 1 and 3 have a zero
mean and skewness. Case 4 is the only signal that produces a nonzero
skewness.
Time ensembles must be statistically stationary for the moments to have any
meaning. A statistically stationary flow is defined as a flow in which statistical
parameters do not depend on the interval in time used to evaluate them. An
example pressure signal is given in Figure 1.8.
Figure 1.8 Unsteady pressure measurement.
It is very difficult to determine if the signal is stationary by visual inspection. Four
overlapping data windows of 4096 samples each are used to assess whether this
flow is stationary. The overlapped windowing technique is effective if at least two
Window 2 Window 4
Window 3 Window 1
![](asset-10.png)
1-12
of the windows are completely independent as is the case for Window 1 and
Window 4. The number of windows and the number of samples in each window
is generally determined by the number of total samples available in the signal.
Fourier transforms are useful for transforming the signal between the time
domain and the frequency domain. The Fourier transform of the four data
windows and the average of the four transforms are shown in Figure 1.9.
Figure 1.9 Fourier transform of pressure signal.
Sound pressure level is proportional to the log of the variance of the pressure
signal. The Fourier transforms for the individual windows show a large amount of
disagreement. This is due to the discreet nature of the transform and the limited
number of samples in each window. The statistical moments and the frequency
of the maximum spectral peak for each of the data windows and the average of
the windows are shown in Table 1.2. The maximum error in Table 1.2 is defined
as the error between the individual data windows and the window average. This
signal is reasonably stationary for the moments examined. It should be noted
that demonstrating that a signal is stationary with respect to a given moment
does not assure that the signal is stationary for higher moments. The windowing
technique can be applied to both experimental and computational data. This
technique can be used to assess “convergence” for unsteady computations6.
![](asset-11.png)
1-13
Average
Pressure
σ
S K Frequency of
Peak SPL (hz)
Window 1 1.2980 0.20586 0.33207 2.3974 433.63
Window 2 1.2957 0.20787 0.36906 2.4891 461.61
Window 3 1.3003 0.20163 0.29161 2.3725 447.62
Window 4 1.2865 0.20218 0.29161 2.5564 461.61
Window
Average 1.2952 0.02044 0.32309 2.4539 451.12
Max Error 0.66% 1.71% 14.23% 4.18% 3.88%
Table 1.2 Statistics from pressure signal data windows.
1.4 Turbulent Length and Time Scales
It is useful to quantify the range of length and time scales associated with a
turbulent flow. These estimates of the maximum and minimum eddy length and
time scales can be used in selecting the appropriate computational grid spacing
and time step for a given problem.
Turbulent scales are defined in terms of the turbulent kinetic energy (k), turbulent
dissipation (
ε
), and the kinematic viscosity (
ν
). The turbulent kinetic energy and
turbulent dissipation will be derived in Chap. 2. These two quantities are usually
associated with the larger scales of turbulence. From dimensional analysis, the
large scale, or turbulence producing, eddies (see Fig. 1.2) have length scales
characterized by
ε
2
3
k
L= (1.13)
The large scale eddies have time scales of the order of
ε
k
=Τ (1.14)
The smallest scale, or dissipative, eddies (the Kolmogorov scale) have length
and time scales given by
![](asset-12.png)
1-14
4
1
3
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
=
ε
ν
λ
(1.15)
2
1
⎟
⎠
⎞
⎜
⎝
⎛
=
ε
ν
τ
(1.16)
The ratio of the smallest to largest eddy length and time scales is then
4
3
−
=t
R
L
λ
(1.17)
2
1
−
=
Τt
R
τ
(1.18)
where Rt is the turbulent Reynolds number (k2/(
νε
)). These relationships indicate
that the range of turbulent scales may span orders of magnitude for high
Reynolds number flows. Since these turbulent length scales are much smaller
than the physical scales (i.e. wing chord, channel height) associated with a flow
of interest, it is easy to see that large amounts of grid points would be required to
fully simulate a high Reynolds number turbulent flow.
1.5 Vorticity and Strain Tensors
Two tensors that are often used in the development of turbulence models are the
vorticity and strain tensors. The vorticity tensor is defined as
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
=Ω
i
j
j
i
ij x
u
x
u
2
1 (1.19)
The trace of the vorticity tensor (the sum of the diagonal elements) is equal to
zero since Ω11=Ω22=Ω33=0. The largest magnitudes of vorticity typically occur
near the wall in a boundary layer and in the core of a vortex. The strain tensor is
defined as
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
u
x
u
S2
1 (1.20)
The trace of the strain is U
r
•∇ , which is zero for steady incompressible flows.
The largest magnitudes of strain also occur near the wall in a boundary layer.
The magnitude of the strain is almost zero in the core of a vortex. For typical
boundary layers and shear layers the magnitude of the vorticity is equal to the
![](asset-13.png)
1-15
magnitude of the strain since only one of the cross derivatives will dominate both
tensors.
1.6 Classification of Turbulence Models
Turbulence models can be classified by what turbulent scales they choose to
model and what scales they choose to simulate by solving the unsteady Navier-
Stokes equations on a computational grid. The traditional Reynolds Averaged
Navier-Stokes (RANS) use a time averaging process to remove the necessity of
simulating all of the scales of the turbulence spectrum. The RANS approach
uses one length scale to characterize the entire turbulent spectrum, as shown in
Fig 1.10.
The use of a single length scale places a tremendous burden on the turbulence
modeler since it can be difficult to find one length scale that is appropriate for all
cases. When this can be accomplished the flow can be treated as a steady flow
since all the unsteadiness is assumed to occur at scales below the computational
grid size and are handled by the turbulence model. This allows for the use of time
marching numerical algorithms with large numerical dissipation since the object
of the calculation is to dissipate all wavelengths (whether they are of physical or
numerical origin) of the unsteady flow during the convergence process.
Numerical algorithms that converge to a steady state solution rapidly may not be
suited for unsteady flow simulations because these algorithms contain large
amounts of numerical dissipation and may overdamp the actual flow. The RANS
approach is discussed in more detail in Chapter 2.
RANS
LES
Figure 1.10 Region of the turbulent energy spectrum
modeled by different turbulence model approaches.
Log(Wave number)
Log(Energy)
Hybrid RANS/LES
DNS
1-16
Direct numerical simulation (DNS)7 attempts to simulate all of the scales of
turbulence and model nothing. The grid resolution and the maximum allowable
time step for a DNS calculation must be small enough to capture the Kolmogorov
scales (Eq. 1.15 and 1.16) of the turbulent flow. The inverse relationship of these
scales with Reynolds number shown in Eq. 1.17 and 1.18 indicate that the grids
must become increasingly fine as Reynolds number is increased. A DNS
solution is inherently unsteady, so it must be run for long periods to assure that
the resulting solution is statistically stationary and free from the initial conditions
provided to start the computation. The numerical algorithm used in the DNS
solution process must have very low numerical dissipation in order to allow all the
wavelengths of the turbulent flow to naturally persist. The numerical dissipation
must be much smaller than the molecular viscosity or it will manifest itself as an
increase in the effective molecular viscosity, and subsequently as a decrease in
the effective Reynolds number8. Most of the higher order numerical algorithms
currently used in DNS are only valid on computational grids that contain uniform
spacing. Grid stretching tends to reduce the order of the numerical algorithm and
increase the numerical dissipation of the algorithm. The requirement of nearly
uniform grids with small cell sizes leads to extremely large numerical grids even
for small Reynolds numbers. DNS computations also require a model for the
inflow free stream disturbances. For channel flows, the inflow turbulence can be
obtained by recycling the outflow turbulence. To date DNS calculations have
only been performed for low Reynolds numbers and simple geometries.
Large Eddy Simulation (LES) is an attempt to move beyond DNS by modeling
only the smallest turbulent scales in a problem9,10,11 as can be seen in Fig. 1.10.
LES will be discussed in detail in Chapter 7. The smaller turbulent scales are
more nearly isotropic, so they may be modeled with fairly simple turbulence
models. Hence a simple turbulence model is used to simulate the “sub-grid”
turbulence (turbulent scales that can not be realized on the computational grid)
and the Navier-Stokes equations are solved for the remaining scales. This
requires a spatial filtering of the turbulence spectrum. Like DNS, LES solutions
are unsteady, and care must be taken to choose a numerical algorithm with low
numerical dissipation or a decrease in effective Reynolds number will be seen in
the calculation. LES solutions must also be run a large number of time steps to
eliminate starting transients and to allow the solution to become statistically
stationary. While LES allows for simulating much higher Reynolds numbers than
DNS for the same computational resources, it is still not capable of simulating
flight Reynolds numbers with reasonable engineering turn-around times. LES
has been found useful for many lower Reynolds number applications such as
combustion calculations where the chemical reactions are driven by the unsteady
turbulent mixing of the fuel and the oxidizer.
Recently a new class of turbulence models has arisen for unsteady high
Reynolds flows. These turbulence models are called hybrid RANS/LES models
and will be discussed in detail in Chapter 8. These hybrid turbulence models are
extensions of the LES models in which modified RANS turbulence models are
![](asset-15.png)
1-17
used as sub-grid models for the Navier-Stokes equations. The goal for most
applications is to use RANS to calculate the boundary layer where extremely
small turbulent scales are present and use an LES like model to simulate the
smaller turbulent scales away from the body. Hence the large scale turbulent
structures away from the body are simulated by the unsteady Navier-Stokes
equations. These models serve as a bridge between traditional RANS models
and LES models as shown in Fig. 1.9. This allows for coarser grids than LES
since the RANS-type sub-grid model is valid for non-isotropic scales of
turbulence. These models require some sort of spatial filtering to determine the
local value for the sub-grid turbulent viscosity. Hybrid RANS/LES models do not
need the extremely low numerical dissipation required by DNS and LES and
hence can be robust and easy to implement.
Chapter 1 References:
1. Launder, B. E., “An Introduction to the Modeling of Turbulence,” VKI Lecture
Series 1991-02, March 18-21, 1991.
2. Schlichting, H., Boundary-Layer Theory, McGraw-Hill, Inc., 1979.
3. Morkovin, M., "Bypass - Transition Research Issues and Philosophy," in
Instabilities and Turbulence in Engineering Flows, Netherlands: Kluwer
Academic Publishers, 1993, pp. 3-30.
4. Mack, L., "Linear Stability Theory and the Problem of Supersonic Boundary-
Layer Transition," AIAA Journal, Vol. 13, Mar. 1975, pp. 278-289.
5. Mack, L., "Stability of Axisymmetric Boundary Layers on Sharp Cones at
Hypersonic Mach Numbers," AIAA-87-1413, 1987.
6. Nichols, R., “Comparison of Hybrid RANS/LES Turbulence Models for a
Circular Cylinder and a Cavity,” AIAA Journal, Vol. 44, No. 6, June 2006, pp.
1207-1219.
7. Mansour, N., “The Use of Direct Numerical Simulation Data in Turbulence
Modeling,” AIAA-91-0221, Jan. 1991.
8. Moin, P., “Numerical and Physical Issues in Large Eddy Simulation of
Turbulent Flow,” JSME International Journal, 1998, Series B, Vol. 41, No. 2,
pp. 454-463.
9. Breuer, M., “Large Eddy Simulation of the Subcritical Flow past a Circular
Cylinder: Numerical and Modeling Aspects,” International Journal of
Numerical Methods in Fluids, 1998, Vol. 28, pp. 1281-1302.
10. Rodi, W., “Large Eddy Simulation of the Flows past Bluff Bodies,” in:
Launder, B. and Sandham, N. (eds) Closure Strategies for Turbulent and
Transitional Flows, Cambridge University Press, Cambridge, 2000.
11. Menon, S. “A Numerical Study of Secondary Fuel Injection Techniques for
Active Control of Combustion Instability in a Ramjet,” AIAA-92-0777, January
1992.
![](asset-16.png)
2-1
2.0 Classical Turbulence Closure Methodologies
The direct numerical solution of the Navier-Stokes equations for high Reynolds
numbers flow requires a grid fine enough to capture all of the turbulent length
scales and a time step small enough to capture the entire turbulent spectrum.
The computational requirements for such calculations are well beyond the
current state of computers. Obviously, some approximations to the Navier-
Stokes equations are required to allow practical engineering solutions to high
Reynolds number flow problems. Two averaging techniques used to simplify the
Navier-Stokes equations are Reynolds and Favre averaging. Both of these
approaches will be discussed in this chapter. Both of these approaches results
in time-averaged correlation terms that must be modeled to close the equation
set.
Two other methods of closure are large eddy simulation (LES) and hybrid
RANS/LES modeling. Both of these techniques use spatial filtering to decide
what turbulent scales should be modeled and what scales can be simulated.
LES will be discussed in detail in Chapter 8. The hybrid approach will be
discussed in Chapter 9.
2.1 Reynolds Averaging
Osborne Reynolds1 suggested that for a statistically stationary flow (see Chapter
1) the variables in the Navier-Stokes equations could be decomposed into time-
averaged and turbulent-fluctuation terms. The velocity can then be written as
),()(),( txuxUtxU ′
+= (2.1)
where the average velocity U is defined as
∫−
=T
T
Udt
T
U2
1 (2.2)
and the time average of the fluctuating velocity u
′
is
()
0
2
1
2
1=−=
′
=
′∫∫ −− dtUU
T
dtu
T
uT
T
T
T (2.3)
Note that the time scale of the integration T must be much larger than the
turbulent time scales discussed in Chapter 1 or the flow will not be statistically
stationary.
In the Reynolds averaging process the unsteady behavior of the turbulent flow is
replaced by steady-state correlations terms that must then be calibrated with
turbulent flow measurements. This has the advantage of significantly reducing
![](asset-17.png)
2-2
the computer time requirements for obtaining engineering quality values of
steady state variables. It is often difficult to obtain empirical correlations that are
valid for a wide range of turbulent flows, and hence the Reynolds Averaged
Navier-Stokes (RANS) turbulence models must often be tuned to a particular
class of flows.
Assuming that fluctuations in density (
ρ
′
), viscosity (
μ
′
), and thermal
conductivity ( T
κ
′) are negligible (which is often not the case for high speed and
reacting flows) and applying Reynolds averaging to the Navier-Stokes equations
yields
()
0=
∂
∂
+
∂
∂
j
j
U
xt
ρ
ρ
(2.4)
()
()
j
ij
i
ji
j
ji
j
ixx
P
uu
x
UU
xt
U
∂
∂
+
∂
∂
−=
∂
∂
+
∂
∂
+
∂
∂
σ
ρρ
ρ
'' (2.5)
()
()
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
−
∂
∂
+
∂
∂
=
⎟
⎠
⎞
⎜
⎝
⎛++
∂
∂
+
∂
∂
+
∂
∂
j
T
j
iji
j
iji
j
jiijiij
j
j
j
x
T
x
u
x
U
x
uuuuuUeu
x
HU
xt
E
κσσ
ρρρρ
ρ
''
'''''''
2
1
(2.6)
where
⎟
⎠
⎞
⎜
⎝
⎛−= ijkkijij SS
δμσ
3
1
2 (2.7)
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
U
x
U
S2
1 (2.8)
kUUeE ii ++= 2
1 (2.9)
ρ
P
EH += (2.10)
''
2
1
iiuuk = (2.11)
![](asset-18.png)
2-3
Note that the turbulent kinetic energy (k) has been included in the definition of
total energy (
E
). The Reynolds averaging has introduced new unknown
correlations k, '' jiuu , ''eui, ''' jii uuu , and '' iji
u
σ
.
When the velocity correlation '' jiuu is multiplied by the density
ρ
it represents the
transport of momentum due to the fluctuating (i.e. turbulent) motion of the fluid.
The term '' jiuu
ρ
is the transport of xi momentum in the direction of xj (or vice
versa). This term acts as a stress on the fluid and is therefore called turbulent or
Reynolds stress. In most regions of a turbulent flow, the turbulent stresses are
much larger than their laminar counter parts ⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
j
i
x
U
ν
. The Reynolds stress is a
symmetric tensor. The trace of the stress tensor is two times the turbulent kinetic
energy, and is coordinate system independent. The off diagonal terms of the
Reynolds stresses are coordinate system dependent. The velocity-energy
correlation
ρ
'' eui can be thought of as a turbulent heat flux. The ''' jii uuu
ρ
triple
correlation represents turbulent transport or diffusion. The '' iji
u
σ
term is the
turbulent dissipation.
Closure of the equation set requires correlations be developed for these new
terms. These terms are often difficult to measure experimentally. The turbulent
dissipation, '' iji
u
σ
, includes correlations of the perturbed velocity and spatial
derivatives of the perturbed velocity. This quantity cannot be measured
experimentally. Transport equations can be developed for the turbulent stresses
and turbulent dissipation to close the system. This approach to closure is called
Reynolds Stress Modeling (RSM) and will be discussed in Chapter 6. A
simplified approach to modeling the turbulent stresses in which the RSM
equations are replaced with algebraic relationships is called Algebraic Stress
Modeling (ASM) and will also be discussed in Chapter 6.
The Reynolds stresses must meet certain constraints to be physically plausible.
Schumann2 introduced the realizability constraint. This constraint specifies that
all the component energies of the turbulent kinetic energy (the diagonal terms of
the Reynolds stress tensor) remain non-negative and all off-diagonal
components of the Reynolds stress tensor satisfy Schwartz’s inequality. This
can be written as
0
'' ≥
iiuu (2.12)
0
2
'''''' ≥− jijjii uuuuuu (2.13)
![](asset-19.png)
2-4
()
()
0
'
3
'
1
'
2
'
2
'
3
'
2
'
2
'
1
'
3
'
1
'
3
'
1
'
3
'
2
'
3
'
3
'
2
'
1
'
2
'
1
2
'
3
'
2
'
3
'
3
'
2
'
2
'
1
'
1
≥−
+−−
⎟
⎠
⎞
⎜
⎝
⎛−
uuuuuuuuuu
uuuuuuuuuuuuuuuuuu (2.14)
2.2 Favre Averaging
Reynolds averaging for compressible flows produces a very complicated set of
equations because the fluctuations of density, pressure, and temperature must
also be accounted for. Favre3 averaging or Favre3 filtering has been used to
simplify the equations set. Favre averaging is a density weighted averaging
process in space. The Favre filter can be defined for any variable as
),()(
~
),( txuxUtxU ′′
+= (2.15)
where
() ()
()
ρ
ρ
ρ
ρ
U
dtt
dttUt
UT
T
T
T== ∫
∫
−
−
~ (2.16)
Note that the over bar signifies a time averaged quantity, and the tilde represents
a Favre filtered quantity. Other auxiliary relations include
UUU
u
ρρρ
ρ
==
=
′′
~~ 0 (2.17)
Thus, using Favre averaging,
UuUtxUtx
~
~
),(),(
ρρρρ
=
′′
+= (2.18)
The equivalent expression using Reynolds averaging would yield
uUtxUtx ′′
+=
ρρρ
),(),( (2.19)
Note that under these definitions the Reynolds average perturbation 0=
′
u, but
the Favre averaged perturbation 0≠
′′
u.
Applying this filtering operation to the Navier-Stokes equations, and assuming
that the filtering commutes with the derivative operation, yields
![](asset-1a.png)
2-5
()
0
~=
∂
∂
+
∂
∂
i
i
U
xt
ρ
ρ
(2.20)
()
() ()
jiji
ji
ji
j
ji
j
i
xx
P
uu
x
UU
xt
U
σσρρ
ρ
′′
+
∂
∂
+
∂
∂
−=
′′′′
∂
∂
+
∂
∂
+
∂
∂~
~~
~
(2.21)
()
()
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
′′
∂
+
∂
∂
∂
∂
−
′′
+
′′
∂
∂
+
∂
∂
=
⎟
⎠
⎞
⎜
⎝
⎛′′′′′′
+
′′′′
+
′′
∂
∂
+
∂
+
∂
∂
j
T
j
T
j
iji
ij
i
j
iji
j
iij
j
iij
j
j
j
x
T
x
T
x
Uu
x
U
x
uuuuuUeu
x
HU
xt
E
κκσσσ
ρρρρ
ρ
~
~
~
~
2
1
~~~
~
(2.22)
For almost all flows ijij
σσ
′′
>>
~
, and hence the ij
σ
′′ terms are usually neglected.
Similarly, for most flows 2
2
2
2
~
jj x
T
x
T
∂
′′
∂
>>
∂
∂, and hence the 2
2
j
x
T
∂
′′
∂ terms are usually
neglected. The Favre averaged equations are similar in appearance to the
incompressible Reynolds averaged Navier-Stokes equations, but the correlations
are not the same. Favre averaging should only be used in flows with large
density fluctuations such as in combustion or hypersonic flow. In other cases
the simpler Reynolds averaging should be used.
2.3 Boussinesq Approximation
Boussinesq4 suggested that the turbulent stresses can be treated in an
analogous form to the viscous stresses in laminar flows. The Reynolds stresses
are modeled as
ij
i
j
j
i
tji k
x
U
x
U
uu
δρμρ
3
2
'' −
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=− (2.23)
The Boussinesq approximation states that the Reynolds stresses are
proportional to the local mean flow strain rate ⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
i
j
j
i
x
U
x
U. The eddy viscosity
(
μ
t) is the proportionality factor for the Boussinesq approximation. This is
sometimes called a local equilibrium assumption since the turbulent stresses are
assumed to be proportional to the mean flow strain rates. The last term in Eq.
2.23 was added so that the normal stresses would sum to the turbulent kinetic
energy (k). The Boussinesq approximation reduces the turbulence modeling
process from finding the six turbulent stress components as discussed above to
determining an appropriate value for the eddy viscosity.
![](asset-1b.png)
2-6
The eddy viscosity is a property of the flow field and not a property of the fluid.
The eddy viscosity can vary significantly form point to point within the flow, and
also from flow to flow. The eddy viscosity concept was developed assuming a
relationship or analogy exists between molecular motion and turbulent motion.
Although turbulent eddies can be thought of as discrete lumps of fluid that collide
and exchange energy, they are not rigid bodies and their mean free paths are not
small compared to the eddy size, and hence turbulence does not really satisfy
the constraints of kinetic gas theory. In spite of the theoretical weakness of the
eddy viscosity concept, it does produce reasonable results for a large number of
flows.
The eddy viscosity concept has been successful in the prediction of flows in
which the shear stress ( '' jiuu
ρτ
=) is the turbulent stress of greatest importance.
This class of flows includes simple shear layers and attached boundary layers.
The eddy viscosity is introduced as a scalar in Eq. 2.23 so that the eddy viscosity
is the same for all stress components. The assumption of isotropic eddy
viscosity is a simplification that limits the performance of this class of turbulence
models in complex flows.
In a direct analogy to the turbulent momentum transport defined by the
Boussinesq approximation (Eq. 2.23), the turbulent heat or mass transport is
often assumed to be proportional to the gradient of the transported quantity
i
ix
u∂
∂
Γ=−
ϑ
ϑ
(2.24)
where Γ is the turbulent diffusivity of heat or mass. Like the eddy viscosity, Γ is
not a property of the fluid but is a function of the flow. Reynolds analogy
between heat or mass transport and momentum transport suggests that Γ is
proportional to the eddy viscosity
t
t
σ
ν
=Γ (2.25)
where
σ
t is the turbulent Prandtl or Schmidt number. Experiments have shown
that the
σ
t varies little throughout a flow field and hence it is usually taken as a
constant.
Using the Boussinesq approximation and the assumption of a constant turbulent
Prandtl number, the Reynolds averaged Navier-Stokes equations (Eq. 2.4, 2.5,
and 2.6) can be written as
()
0=
∂
∂
+
∂
∂
j
j
U
xt
ρ
ρ
(2.26)
![](asset-1c.png)
2-7
()
j
ij
i
ji
j
i
xx
P
UU
xt
U
∂
∂
+
∂
∂
−=
∂
∂
+
∂
∂
σ
ρ
ρ
(2.27)
()()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
−
∂
∂
=
∂
∂
+
∂
∂
jt
t
j
iji
j
j
jx
T
x
U
x
HU
xt
E
PrPr
μ
μ
σρ
ρ
(2.28)
where
()
⎟
⎠
⎞
⎜
⎝
⎛−+= ijkkijtij SS
δμμσ
3
1
2 (2.29)
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
U
x
U
S2
1 (2.30)
kUUeE ii ++= 2
1 (2.31)
ρ
P
EH += (2.32)
where k is the turbulent kinetic energy. The turbulent kinetic energy is often
ignored in the definition of total energy (Eq. 2.31). Under the assumption of a
constant turbulent Prandtl number (Prt), the equation set is closed once the eddy
viscosity is defined. Some common approaches to defining eddy viscosity
models are discussed in Chapters 3, 4, and 5.
Chapter 2 References:
1. Reynolds, O., “On the Dynamical Theory of Incompressible Viscous Fluids
and the Determination of the Criterion,” Philosophical Transactions of the
Royal Society of London, Series A, Vol. 186, pp. 123-164, 1895.
2. Schumann, U., “Realizability of Reynolds Stress Turbulence Models,” Physics
of Fluids, Vol. 20, pp. 721-725, 1977.
3. Favre, A., “Equations Fondanmentales des Fluids a Masse Variable en
Ecoulements Turbulents,” Pages 24-78 of Favre, A., Kovasznay, L., Dumas,
R., Gaviglio, J., and Coantic, M. (eds.), La Turbulence en Micanique de
Fluids, CNRS, Paris, 1962. English edition: The Mechanics of Turbulence,
Gordon and Breach, New York, 1964.
4. Boussinesq, J., Theorie de l’ecoulement Turbillonnant et Tumultueux des
Liquides Dans les lits Rectilignes a Grande Section, Gaulthier-Villars, 1897.
![](asset-1d.png)
3-1
3.0 Turbulent Boundary Layer and Shear Layer Theory
A brief review of simple theory for turbulent boundary layers and shear layers will
aid in the understanding of the development and calibration of turbulence
models. Most of these relationships are based on empirical correlations of the
shape of the velocity profile across the boundary layer or shear layer. This
shape can be used in conjunction with the Navier-Stokes equations to back out
relationships for the turbulent stresses or the eddy viscosity.
3.1 Boundary Layer Theory
Turbulent boundary layers are usually described in terms of several
nondimensional parameters. The boundary layer thickness,
δ
, is the distance
from the wall at which viscous effects become negligible and represents the edge
of the boundary layer. Two integral parameters across the velocity profile are the
displacement thickness *
δ
and the momentum thickness
θ
.
dy
u
u
ee
∫
∞
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−=
0
*1
ρ
ρ
δ
(3.1)
dy
u
u
u
u
eee
∫
∞
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−=
0
1
ρ
ρ
θ
(3.2)
The integration is performed normal to the wall and the subscript e is used to
denote the edge of the boundary layer at y=
δ
. The displacement thickness is a
measure of the increased thickness of a body due to the velocity defect of the
boundary layer. The momentum thickness is the distance that, when multiplied
by the square of the free-stream velocity, equals the integral of the momentum
defect across the boundary layer.
If a simple power law velocity profile as shown in Eq. 3.3 is assumed and the flow
is incompressible, the relationships in Eqs. 3.4-3.7 can be obtained from
approximations of the Navier-Stokes equations.
7/1
⎟
⎠
⎞
⎜
⎝
⎛
=
δ
y
u
u
e
(3.3)
5/1
Re
371.0
x
x=
δ
(3.4)
5/1
*
Re
046.0
x
x=
δ
(3.5)
![](asset-1e.png)
3-2
5/1
Re
036.0
x
x=
θ
(3.6)
5/12
2
1Re
0577.0
xee
w
fu
C==
ρ
τ
(3.7)
Cf is called the skin friction. The distance from the leading edge in the
streamwise direction is given by x. Both the boundary layer growth rate and the
skin friction decrease as the Reynolds number increases.
Reynolds analogy can be used to develop a relationship for heat transfer.
Reynolds analogy says that the ratio of the shear stress to the heat transfer is a
constant near the wall. Thus the Nusselt number can be defined as
PrRe
2x
f
x
C
k
hx
Nu == (3.8)
Here h is the heat transfer coefficient, k is the thermal conductivity, and Pr is the
Prandtl number. This relationship is independent of the equation used to
determine skin friction.
The power law relationship is not extremely accurate, but is useful for developing
some useful turbulent boundary layer relationships. A more accurate relationship
for skin friction for adiabatic incompressible flow on a flat plate is given by White
and Christoph1
()
x
f
CRe056.0ln
42.0
2
= (3.9)
This relationship is only applicable in the turbulent region of a boundary layer,
and does not apply in the laminar or transitional regions of the boundary layer.
The skin friction is affected by a number of parameters, including pressure
gradient, surface roughness, compressibility, and surface heat transfer. Adverse
pressure gradients cause the skin friction to be reduced as the boundary layer is
pushed toward separation. Boundary layer separation occurs when the skin
friction becomes negative. High values of skin friction are an indication of a very
stable (i.e. difficult to separate) boundary layer. Unfortunately high values of skin
friction also equate to high values of viscous drag. The effect of both
compressibility and heat transfer on the flat plate turbulent skin friction are shown
in Fig. 3.1. Both a hot wall and compressibility tend to reduce the skin friction
over the incompressible adiabatic value.
![](asset-1f.png)
3-3
Figure 3.1 Effect of compressibility and heat transfer on the skin friction on a flat
plate.
As mentioned above, the power law relationship in Eq. 3.3 is not very accurate.
Better approximations of the velocity profile shape are generally written in terms
of the parameters u+ and y+ defined as
τ
u
u
u=
+ (3.10)
w
wyu
y
μ
ρ
τ
=
+ (3.11)
and uτ is the friction velocity
w
w
u
ρ
τ
τ
= (3.12)
The subscript w denotes the value at the wall and
τ
w is the wall shear stress
defined by
![](asset-20.png)
3-4
w
wy
u
∂
∂
=
μτ
(3.13)
The boundary layer velocity profile can be divided into four regions. The
incompressible velocity profile in each of the subregions regions of the inner
region shown in Fig. 3.2 is given by
Laminar sublayer 0 < y+ < 5 u+ = y+ (3.14)
Buffer layer 5 < y+ < 30 u+ = 5 ln y+ - 3.05 (3.15)
Log layer 30 < y+ < 1000 u+ =
κ
1 ln y+ + B (3.16)
The values of
κ
, the von Karmen constant, and B are often debated, but are
generally accepted to be 0.4 and 5.5 respectively. The log layer is also called
the law of the wall.
Figure 3.2 Boundary layer regions.
The y+ value where the profile transitions from the inner to the outer profile varies
with the Reynolds number and the pressure gradient. The outer region is much
more sensitive to pressure gradient. Clauser’s2 equilibrium parameter
β
is often
used to characterize the pressure gradient.
Sublayer Log Layer
Wake
Buffer
Layer
Inner Region Outer
Region
![](asset-21.png)
3-5
x
pe
w∂
∂
=
τ
δ
β
* (3.17)
Coles3 introduced the wake function W given by
⎟
⎠
⎞
⎜
⎝
⎛
=
⎟
⎠
⎞
⎜
⎝
⎛
δ
π
δ
yy
W2
sin2 2 (3.18)
The velocity profile in the outer region is given by
⎟
⎠
⎞
⎜
⎝
⎛
Π
++= ++
δκκ
y
WByu ln
1 (3.19)
where Π is given by
()
75.0
5.08.0 +=Π
β
(3.20)
Spalding4 proposed a composite form for the incompressible velocity profile
given by
(
)
(
)
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡−−−−+=
++
+Π−−++ +
62
1
32 uu
ueeeuy uWB
κκ
κ
κκ
(3.21)
White and Christoph1 give a law of the wall velocity profile that includes the
effects of compressibility, heat transfer, and pressure gradient.
()
⎥
⎦
⎤
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
+
+
−
+−
Γ
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛Θ−Γ
=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛Θ−Γ +
−
+
−
1
1
1
1
ln2
2
sin
2
sin
0
0
0
2/1
0
11
φ
φ
φ
φ
φφ
κ
Q
u
Q
u (3.22)
Where
()
2/1
24Γ+Θ=Q and
(
)
2/1
1+
+= y
αφ
. The values of +
0
y and +
0
u are taken
as 6 and 10 respectively. The parameters Θ, Γ, and
α
represent the effects of
heat transfer, compressibility, and pressure gradient respectively.
τ
ρ
μ
ukT
q
www
ww
=Θ (3.23)
wpTc
ru
2
2
τ
=Γ (3.24)
![](asset-22.png)
3-6
x
p
ue
ww
w
∂
∂
=
τ
τρ
μ
α
(3.25)
Here qw is the wall heat transfer (
w
wy
T
k∂
∂),
μ
w is the wall viscosity, Tw is the wall
temperature,
ρ
w is the wall density, r is the recovery factor (normally taken as the
Prandtl number to the one third power), kw is the wall thermal conductivity, and cp
is the specific heat at constant pressure.
It is not obvious what effect each of the parameters defined in Eqs. 3.23-25 has
on the velocity profile. Fig. 3.3 shows the effect of an adverse pressure gradient
and compressibility has on the velocity profile.
Figure 3.3 Effect of adverse pressure gradient and compressibility on boundary
layer profile shape.
The boundary layer thickens and the skin friction decreases as the pressure
gradient is increased. Compressibility also causes the skin friction to decrease.
Heat transfer effects are shown in Fig. 3.4.
Increasing
α
Increasing
Γ
![](asset-23.png)
3-7
Figure 3.4 Effect of wall heat transfer on boundary layer profile shape.
A cold wall (wall temperature less than the adiabatic wall temperature) creates a
thinner boundary layer and increases the skin friction. A hot wall (wall
temperature greater than the adiabatic wall temperature) thickens the boundary
layer and decreases the skin friction.
The temperature distribution within the inner part of the boundary layer boundary
layer is given by the Crocco-Busemann equation
()
(
)
2
1++ Γ−Θ+= uuTT w (3.26)
where Θ is defined in Eq. 3.23 and Γ is defined in Eq. 3.24. For adiabatic wall
cases,
β
=0 and the Crocco-Busemann equation reduces to
()
⎟
⎠
⎞
⎜
⎝
⎛−
+= 2
wu
2
r1
1TT
γ
(3.27)
The pressure is assumed to be constant in the normal direction from the wall in
the inner part of the boundary layer. Density distributions can be defined based
on the temperature distribution and the equation of state.
3.2 Shear Layer Theory
A free shear layer is always initiated from a surface of some kind. The boundary
layer profile remains for a short period. If no external pressure gradient is
present, the shear layer will eventually become self-similar. A self-similar profile
is one in which the profile shape remains unchanged as you move downstream if
Increasing Θ
(
Cold Wall
)
Decreasing Θ
(
Hot Wall
)
![](asset-24.png)
3-8
the profile is defined in terms of similarity variables. Chapman and Korst5
suggested the similarity variable for two-dimensional shear layers
x
y
σ
η
= (3.28)
where x is the downstream distance from the origin of the shear layer, y is the
normal distance across the shear layer (y=0 denotes the center of the shear
layer, and
σ
is the shear layer spread rate parameter. Brown and Roshko6
suggested using
ω
δ
5.0
*yy
y−
= (3.29)
where the vorticity thickness
δ
ω
is defined as
max
21
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
=
y
u
UU
ω
δ
(3.29)
Here U1 is the velocity at the high speed edge of the shear layer and U2 is the
velocity at the low speed edge of the shear layer. Chapman and Korst5
suggested that the velocity profile was given by
[]
)(15.0
21
*
η
erf
UU
u
u+=
−
= (3.30)
where erf is the error function. This profile shape is valid for both laminar and
turbulent shear layers. Samimy and Elliot7 obtained Laser Doppler Velocimeter
(LDV) data on a shear layer. Measurements were made at several downstream
locations between the trailing edge of the splitter plate and a station 210 mm
downstream of the splitter plate. The flow parameters are given in Table 3.1.
Eq. 3.30 is plotted with data in Fig. 3.5.
T0, K P01,
kPa M1 M
2 M
c U
1,
m/sec U2/U1 ρ2/ρ1 δ1, mm
291.0 314.0 1.80 0.51 0.52 479.5 0.355 0.638 8.0
Table 3.1 Flow parameters for the spatial mixing-layer case
![](asset-25.png)
3-9
Figure 3.5 Shear layer velocity profile.
The shear layer thickness is given by
σ
x
b= (3.31)
The shear layer spread parameter (
σ
) is affected by compressibility. The
accepted value for subsonic flow issuing into quiescent air is
σ
=11. This value
will increase as compressibility effects become greater, and will cause the shear
layer to become thinner. Experimental values8 for
σ
for two-dimensional jets
issuing into quiescent air are shown in Fig. 3.6.
Figure 3.6 Shear layer spread parameter for jet issuing into quiescent air.
![](asset-26.png)
3-10
Subsonic spread parameter experimental results8 as a function of velocity ratio
are shown in Fig. 3.7.
Figure 3.7 Shear layer spread parameter for subsonic jets.
This indicates that the shear layer will become thinner as the ratio of the
velocities of the two streams increases.
Chapter 3 References:
1. White, F. M. and Christoph, G. H., “A Simple New Analysis of
Compressible Turbulent Two-Dimensional Skin Friction Under Arbitrary
Conditions, “ AFFDL-TR-70-133, Feb. 1971.
2. Clauser, F. J., “Turbulent Boundary Layers in Adverse Pressure
Gradients,” Journal of the Aeronautical Sciences, Vol. 21, 1954, pp. 91-
108.
3. Coles, D. “The Law of the Wake in the Turbulent Boundary Layer,” Journal
of Fluid Mechanics, Vol. 1, Part 2, July 1956, pp. 191-226.
4. Spalding, D. B., “A Single Formula for the Law of the Wall,” Journal of
Applied Mechanics, Vol. 28, No. 3, 1961, pp. 444-458.
5. Chapman, A. and Korst, H., “Free Jet Boundary with Consideration of
Initial Boundary Layer,” Proceedings of the Second U. S. National
Congress of Applied Mechanics, The American Society of Mechanical
Engineers, New York, 1954, pp. 723-731.
6. Brown. G. and Roshko, A., “On Density Effects and Large Structure in
Turbulent Mixing Layers,” Journal of Fluid Mechanics, Vol. 64, No. 4,
1974, pp. 775-816.
7. Samimy, M. and Elliot, G., “Effects of Compressibility on the
Characteristics of Free Shear Layers,” AIAA Journal, Vol. 28, No. 3, 1990,
pp. 439-445.
3-11
8. Birch, S. and Eggers, J., “A Critical Review of Experimental Data for
Developed Free Turbulent Shear Layers,” in Free Turbulent Shear Flows,
NASA-SP-321, 1972, pp. 11-37.
![](asset-28.png)
4-1
4.0 Algebraic Turbulence Models
4.1 Inner Eddy Viscosity Model
Most of the earliest turbulence models were based on Prandtl’s mixing length
hypothesis. Prandtl1 suggested that the eddy viscosity could be represented by
y
u
Lmt ∂
∂
=2
ρμ
(4.1)
The Prandtl mixing length relates the eddy viscosity to the local mean velocity
gradient. The secret to successful applications of the mixing length hypothesis is
to find some general method of defining the mixing length. Most algebraic
models divide the boundary layer into an inner and outer region as described in
Chapter 3. The inner layer includes the viscous sublayer, the buffer layer, and
part of the fully turbulent log region. The outer layer includes the remaining part
of the log layer and the wake region. The eddy viscosity in the inner layer follows
Prandtl's form and is given by
()
Ω= 2
m
inner
tL
ρμ
(4.2)
where the mixing length Lm is given by
⎥
⎦
⎤
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
−= +
+
A
y
yLmexp1
κ
(4.3)
where y is the distance normal to the wall,
κ
is the von Karmen constant, and the
vorticity Ω is defined as
x
v
y
u
∂
∂
−
∂
∂
=Ω (4.4)
4.2 Cebeci-Smith Model
Cebeci-Smith2 suggested that the outer eddy viscosity be expressed as
()
*
δραμ
e
outer
tu= (4.5)
Here
α
is usually assigned a value of 0.0168 for flows where the Reynolds
number based on momentum thickness (Reθ) is greater than 5000,
δ
* is the
displacement thickness, and ue is the velocity at the edge of the boundary layer.
The final eddy viscosity is
![](asset-29.png)
4-2
() ()
[
]
outer
t
inner
tt
μ
μ
μ
,min= (4.6)
This model is fairly simple, but it requires knowledge of the conditions at the edge
of the boundary layer and the boundary layer thickness. These quantities are not
always easy to calculate in complicated flows with a Navier-Stokes code since it
is often difficult to define where the boundary layer edge actually occurs.
4.3 Baldwin-Lomax Model
Baldwin-Lomax3 developed a form of the outer eddy viscosity that did not require
knowledge of the conditions at the edge of the boundary layer. This model has
become quite popular for CFD applications. The eddy viscosity in the outer layer
is defined as
()
klebwakecp
outer
tFFKC
ρ
μ
= (4.7)
where Fwake contains the mixing length term and Fkleb is the Klebanoff
intermittency factor. These terms are defined as
{
}
max
2
maxmaxmax /,min FUyCFyF diffwkwake = (4.8)
The quantities Fmax and ymax are taken from the maximum of the F function
defined as
⎥
⎦
⎤
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
−Ω= +
+
A
y
yyF exp1)( (4.9)
and Udiff is given by
(
)
(
)
min
22
max
22 VUVUUdiff +−+= (4.10)
The F function is calculated along a line normal to the wall. The F function for a
typical attached boundary layer is shown in Fig. 4.1.
![](asset-2a.png)
4-3
Figure 4.1 The F function (Eq. 4.9) for an attached boundary layer.
The Klebanoff intermittency factor4 is used to reduct the eddy viscosity to zero at
the outer edge of the boundary layer. The Klebanoff intermittency factor is given
by
1
6
max
5.51)(
−
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+= y
yC
yF kleb
kleb (4.11)
The final eddy viscosity is
() ()
[
]
outer
t
inner
tt
μ
μ
μ
,min= (4.12)
The constants used in the Baldwin-Lomax model are
068.0,4.0,0.1,3.0,6.1,26 ======
+KCCCA wkklebcp
κ
(4.13)
Algebraic models work well for flows that can be characterized by a single length
scale such as attached boundary layers or simple shear layers. The Baldwin-
Lomax model determines the appropriate mixing length from the location of the
peak in the F function. For simple flows there will only be one peak in the F
function. In more complicated flows with multiple shear layers such as separated
boundary layers or wall jets the F function will contain multiple peaks. When
multiple peaks are present it becomes difficult to automatically choose the proper
peak to use, if a single peak can be used to model the turbulent flow. Degani
and Schiff5 recommended a procedure to automatically select the first significant
peak in the F function. This modification to the search procedure for the peak
has been shown to improve the performance of the Baldwin-Lomax for high
angle-of-attack flows with cross-flow separation.
F
y (Fmax,ymax)
![](asset-2b.png)
4-4
4.4 Wake and Jet Model
The Baldwin-Lomax model has been modified for use in wakes and jets (Ref. 6).
First, Udiff in Eq. 4.11 is redefined as
(
)
(
)
max
22
max
22 ydiff VUVUU +−+= (4.14)
The exponential term in the definition of F(y) (Eq. 4.9) is set to zero yielding
Ω= yyF )( (4.15)
The Fwake function (Eq. 4.8) is redefined to be
max
2
max
maxmax Ω
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
Ω
== diff
wake
U
FyF (4.16)
The ymax location in Eq. 4.16 is defined to be the location where max
Ωoccurs.
Finally the Klebanoff intermittency factor Fkleb (Eq. 4.11) is rewritten as
1
6
max
max
5.51
−
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎛
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎛
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
Ω
−
+=
diff
klebkleb U
yy
CF (4.17)
This wake formulation will work for individual simple shear layers. If multiple
shear layers are present, the user must partition the flow and apply the
turbulence model to each shear layer individually. This can easily become
impractical for complicated geometries.
4.5 Algebraic Model Shortcomings
The flow in a two-dimensional channel with a circular arc bump contraction can
be used to demonstrate two shortcomings of the Baldwin-Lomax turbulence
models. These shortcomings are the model’s tendency to switch Fmax peaks and
the lack of any transport terms in the turbulence model. Fig. 4.2 shows the
geometry and the eddy viscosity contours for both the Baldwin-Lomax and
Spalart-Allmaras turbulence models. The eddy viscosity predicted by the
Baldwin-Lomax is seen to reduce itself almost to zero at the start of the bump.
The Spalart-Allmaras one-equation transport turbulence model (discussed in the
next section) predicts a much smoother distribution of eddy viscosity on the
bump. The reason for the anomaly in the Baldwin-Lomax model can be seen in
![](asset-2c.png)
4-5
Fig. 4.3. The location of the Fmax peak moves very close to the wall at the start of
the bump (x=0) and rises back to its original level at the bump trailing edge (x=1).
The relative magnitude of the Fmax peak is also shown if Fig. 4.3, and is seen to
increase in the region of the bump. The eddy viscosity is proportional to the
product of the distance of the peak off the wall and the magnitude of the peak,
and is seen to decrease in the region above the bump because the Fmax peak has
moved very near the wall. The sudden switching of the Fmax peak causes the
turbulence level to be greatly reduced in the favorable pressure gradient region
near the bump leading edge. If the Baldwin-Lomax model included transport
terms, then the higher upstream values of eddy viscosity would be transported
downstream at the beginning of the bump, and the Baldwin-Lomax predicted
eddy viscosity distribution would not demonstrate the discontinuous behavior it
shows here.
Figure 4.2 Eddy viscosity contours for a circular arc bump in a two-dimensional
channel for two different turbulence models.
![](asset-2d.png)
4-6
Figure 4.3 Distance above the wall of the Fmax peak and the Fmax value for a two-
dimensional channel with a circular arc bump.
4.6 Grid Sensitivity for a Flat Plate with Adiabatic Wall
The initial wall spacing of the computational grid and the grid-stretching ratio (the
ratio of the change in grid spacing normal to the wall) can affect the accuracy of
the Baldwin-Lomax model. Figure 4.4 shows the sensitivity of the skin friction to
initial wall spacing for a flat plate. The grid-stretching ratio was 1.2 for all these
cases. The plots include the theoretical skin friction curves of White and of
Spalding.
![](asset-2e.png)
4-7
Figure 4.4 Flat plate skin friction predictions for the Baldwin-Lomax turbulence
model for varying initial wall grid point spacings.
The calculated boundary layer is seen to become fully turbulent around a length
Reynolds number (Rex) of 1x106. The results for y+=0.2 and y+=1 are virtually
identical indicating a grid independent solution. The y+=5 solution shows some
small divergence from the y+=1 solution at the lower length Reynolds numbers
while the y+=10 solution shows large differences from the other solutions.
Predicted velocity profiles for the flat plate boundary layer for various initial wall
grid point spacings are shown in Fig. 4.5. The velocity profile shows little effect
of the initial spacing for all but the y+=10 profile. All of the profiles but the y+=10
profile are in good agreement with the theoretical profile from Spalding. Note
that the theoretical profile does not include the law-of-the wake (see Chapter 3),
and hence the predicted profiles diverge from the theoretical profile in the wake
region of the boundary layer. The predicted eddy viscosity for various initial wall
spacings is shown in Fig. 4.6.
![](asset-2f.png)
4-8
Figure 4.5 Flat plate boundary layer profiles predicted by the Baldwin-Lomax
turbulence model for varying initial wall grid point spacings.
Figure 4.6 Eddy viscosity distribution predicted by the Baldwin-Lomax turbulence
model for varying grid initial wall spacings.
Here again it is seen that the y+=0.2 and the y+=1.0 results are almost identical.
The y+=5 and y+=10 results show the solutions are no longer grid independent at
larger wall spacings.
![](asset-30.png)
4-9
The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig.
4.7. All of these solutions used an initial wall spacing of y+=1.
Figure 4.7 The effect of grid stretching ratio on the skin friction for a flat plate
boundary layer using the Baldwin-Lomax turbulence model.
There seems to be very little effect of grid stretching for this case indicating that
the initial wall spacing is the more critical parameter for skin friction predictions
for flat plates with the Baldwin-Lomax turbulence model. This is also the case for
the velocity profile, as seen in Fig. 4.8. The eddy viscosity does change as the
stretching ratio increases as shown in Fig. 4.9. It is interesting to note that a
wide range of eddy viscosity distributions have little effect on skin friction and the
velocity profile for a flat plate boundary layer. This insensitivity to the absolute
eddy viscosity level is probably a major reason why the eddy viscosity concept
has worked so well in practice.
![](asset-31.png)
4-10
Figure 4.8 The effect of grid stretching ratio on the velocity profile for a flat plate
boundary layer using the Baldwin-Lomax turbulence model.
Figure 4.9 The effect of grid stretching ratio on the eddy viscosity distribution for
a flat plate boundary layer using the Baldwin-Lomax turbulence model.
4.7 Grid Sensitivity for Axisymmetric Bump
A second example of the grid sensitivity of the Baldwin-Lomax turbulence model
that includes a pressure gradient is the NASA Ames transonic axisymmetric
![](asset-32.png)
4-11
bump experiment (Ref. 7). The geometry is shown in Fig. 4.10. The model
consisted of a sharp-lipped hollow cylinder with a 15.2 cm outer surface
diameter. The bump was a circular arc 20.3 cm long and 1.9 cm high that begins
60.3 cm downstream of the cylinder leading edge. The upstream intersection of
the bump and cylinder was faired with a circular arc. The test was run at a Mach
number of 0.875 and a chord Reynolds number of 2.67x106.
c=20.3 cm
h=1.9 cm
D=15.2 cm
M=0.875 Shock
Recirculation
Figure 4.10 Geometry for the transonic axisymmetric bump.
The effect of initial grid spacing on the pressure coefficient distribution along the
bump is shown in Fig. 4.11. The stretching ratio was 1.2 for these cases. The
pressure coefficient seems to be relatively insensitive to the initial grid spacing,
with the y+=10 and y+=20 curves being slightly displaced from the other curves.
The velocity distribution at the aft junction of the bump and the cylinder (x/c=1) is
shown in Fig. 4.12. The y+=20 solution predict a larger velocity in the reverse
flow region than the other solutions. Grid stretching effects on the pressure
coefficient distribution along the bump is shown in Fig. 4.13. The initial grid
spacing was 1.2 for these cases. The pressure distribution coefficient changes
slightly as the grid-stretching ratio is increased to 1.5. The solution in the
separated region differs greatly for a grid-stretching ration of 2.0. The effect on
the velocity distribution at x/c=1 is shown in Fig. 4.14. As with increasing initial
grid spacing, increasing the grid spacing increases the size and the magnitude of
the separated flow region.
![](asset-33.png)
4-12
Figure 4.11 The effect of initial grid wall spacing on the pressure coefficient for
the axisymmetric bump using the Baldwin-Lomax turbulence model.
Figure 4.12 The effect of initial grid wall spacing on the velocity distribution at
x/c=1 for the axisymmetric bump using the Baldwin-Lomax turbulence model.
![](asset-34.png)
4-13
Figure 4.13 The effect of grid stretching on the pressure coefficient for the
axisymmetric bump using the Baldwin-Lomax turbulence model.
Figure 4.14 The effect of grid stretching on the velocity profile at x/c=1 for the
axisymmetric bump using the Baldwin-Lomax turbulence model.
The results for the Ames axisymmetric bump indicate that the grid-stretching ratio
is a critical parameter when pressure gradients are present in the flow. The
stretching ratio should probably be kept between 1.2 and 1.3 to assure that the
grid can capture the pressure gradient effects. This is true in both structured and
unstructured grids.
![](asset-35.png)
4-14
4.8 Grid Sensitivity for a Flat Plate with Heat Transfer
Calculating heat transfer accurately can be more difficult than predicting skin
friction. This can be seen in the subsonic flat plate example when the wall
temperature is specified to be 1.5 times the free-stream temperature. The
sensitivity of the skin friction and heat transfer result with varying initial grid wall
spacing is shown in Fig. 4.15 and Fig. 4.16. The grid stretching ratio was fixed at
1.2 for these results. Both the skin friction and heat transfer seem to be relatively
insensitive to the wall spacing for wall spacings less than y+=5 when no pressure
gradient is present.
Figure 4.15 The effect of wall spacing on the skin friction on a flat plate with heat
transfer using the Baldwin-Lomax turbulence model.
![](asset-36.png)
4-15
Figure 4.16 The effect of wall spacing on the heat transfer (Stanton number) on
a flat plate using the Baldwin-Lomax turbulence model.
The velocity and temperature profiles for a length Reynolds number (Rex) of
1.0x107 are shown in Fig. 4.17 and 4.18 respectively. Both the velocity and
temperature profiles are relatively insensitive to wall spacing for this model when
no pressure gradient is present.
Fig. 4.17 The effect of wall spacing on the velocity profile on a flat plate with heat
transfer using the Baldwin-Lomax turbulence model.
![](asset-37.png)
4-16
Fig. 4.18 The effect of wall spacing on the temperature profile on a flat plate with
heat transfer using the Baldwin-Lomax turbulence model.
Grid stretching effects on skin friction and heat transfer predictions are shown in
Fig. 4.19 and 4.20. The effect of grid stretching on the velocity and temperature
profiles is shown in Fig. 3.22 and 3.23 respectively. The initial wall spacing was
held at y+=0.1 for these calculations. The results reach a grid independent result
for a stretching ratio of less than 1.3.
Figure 4.19 The effect of grid stretching on the skin friction on a flat plate with
heat transfer using the Baldwin-Lomax turbulence model.
![](asset-38.png)
4-17
Figure 4.20 The effect of grid stretching on heat transfer (Stanton number) on a
flat plate using the Baldwin-Lomax turbulence model.
Figure 4.21 The effect of grid stretching on the velocity profile on a flat plate with
heat transfer using the Baldwin-Lomax turbulence model.
![](asset-39.png)
4-18
Figure 4.22 The effect of grid stretching on the temperature profile on a flat plate
with heat transfer using the Baldwin-Lomax turbulence model.
4.9 Grid Sensitivity for a Nozzle with Heat Transfer
Flow through a supersonic nozzle with a constant temperature wall can serve as
a test case for evaluating the performance of the turbulence model in the
presence of strong pressure gradients. Back, Massier, and Gier8 measured the
wall pressure distribution and heat transfer for a converging-diverging nozzle with
a throat diameter of 0.0458 meters and an exit diameter of 0.1227 meters. High-
pressure air was heated by the internal combustion of methanol and flowed along
a cooled constant area duct with a length of 0.4572 meters and a diameter of
0.355 meters before entering the nozzle. The nozzle geometry and boundary
conditions are shown in Fig. 4.23. The gas could be treated as a calorically
perfect gas with a ratio-of-specific heats (γ) of 1.345. The nozzle exit Mach
number was approximately 2.5. The molecular viscosity and thermal conductivity
were assumed to vary according to Sutherland’s law.
Figure 4.23 Nozzle geometry.
Pt=150psia
Tt=1418oR Twall=0.5Tt Supersonic
Outlet
![](asset-3a.png)
4-19
The initial wall spacing was varied with a grid stretching ratio of 1.2 in the
boundary layer. The grid spacing was held constant in the core of the nozzle.
Comparisons of the predicted and measured pressure along the nozzle are
shown in Fig. 4.24 for varying initial wall spacings. The throat is located at
x=0.091 meters. The pressure distribution is seen to be insensitive to the initial
wall spacing. The heat transfer at the wall is shown in Fig. 4.25. The results are
in poor agreement with the data and are quite sensitive to the initial wall spacing,
especially for wall spacing greater than y+=1.
Figure 4.24 The effect of wall spacing on the pressure distribution for a
supersonic nozzle with heat transfer using the Baldwin-Lomax turbulence model.
Figure 4.25 The effect of wall spacing on the wall heat transfer distribution for a
supersonic nozzle using the Baldwin-Lomax turbulence model.
![](asset-3b.png)
4-20
The poor performance of the Baldwin-Lomax model for this case is due to the
difficulty in choosing the proper Fmax peak as was described in Section 4.5. The
F function at the nozzle throat is shown in Fig. 4.26. Notice that multiple peaks
are present in the function. It is extremely difficult to predict which, if any, of
these three peaks is the proper value for this case. The eddy viscosity
distribution at the nozzle throat is shown in Fig. 4.27. The Spallart-Allmaras eddy
viscosity distribution is included in Fig. 4.26 for comparison. The Baldwin-Lomax
eddy viscosity is much lower than the Spalart-Allmaras distribution and is seen to
cut off prematurely. This low value of eddy viscosity results in a low prediction of
the wall heat transfer.
Figure 4.26 The F function at the nozzle throat.
R/Rthroat
![](asset-3c.png)
4-21
Fig. 4.27 The eddy viscosity distribution at the nozzle throat.
4.10 Summary
Algebraic turbulence models were popular in the 1970’s and 1980’s because of
their simplicity and robustness. As Navier-Stokes CFD applications became
more complex in the 1990’s these models began to lose popularity because of
accuracy limitations for flows that contain multiple shear layers or boundary layer
separation. The difficulties encountered when multiple peaks occur in the F
function have been demonstrated. These models have also lost favor in
unstructured grid applications since they require a velocity profile over multiple
grid points aligned with the flow for successful application.
The eddy viscosity predicted by an algebraic turbulence model is only a function
of the local velocity profile used to generate the F function. Thus the eddy
viscosity relates directly to the local instantaneous vorticity field of that profile and
cannot model the transport of turbulence by the flow. This makes these models
a poor choice for unsteady flows and for flows where turbulent transport is
important.
![](asset-3d.png)
4-22
Chapter 4 References:
1. Prandtl, L., “Bericht uber Untersuchugen zur Ausgebildeten Turbulenz,” Z.
Angew Math., Meth., Vol. 5, pp. 136-139.
2. Cebeci, T, and Smith, A. M. O., Analysis of Turbulent Boundary Layers,
Academic Press, New York, 1974.
3. Baldwin, B. S. and Lomax, H., “Thin Layer Approximation and Algebraic
Model for Separated Turbulent Flows,” AIAA-78-0257, Jan. 1978.
4. Klabenoff, P. S., “Characteristics of Turbulence in a Boundary Layer with
Zero Pressure Gradient,” NACA Rep. 1247, 1955.
5. Degani, D. and Schiff, L. B., “Computation of Supersonic Viscous Flows
Around Pointed Bodies at Large Incidence,” AIAA-83-0034, Jan. 1983.
6. Buning, P. G., Jesperson, D. C., Pulliam, T. H., Klopfer, G. H., Chan, W.
C., Slotnik, J. P., Krist, S. E., and Renze, K. J., “OVERFLOW User’s
Manual,” Version 1.8l, July 1999.
7. Johnson, D. A., “Predictions of Transonic Separated Flow with an Eddy-
Viscosity/Reynolds Shear Stress Closure Model,” AIAA-85-1683, July
1985.
8. Back, L. H., Massier, P. F., and Gier, H. L., “Convective Heat Transfer in a
Convergent-Divergent Nozzle,” International Journal of Heat Mass
Transfer, Vol. 7, 1964, pp. 549-568.
Baldwin-Lomax Application Hints
1. The Baldwin-Lomax model requires that the F function be well
defined. This normally requires that at least three points be located
within the sublayer (y+<10). The first point off the wall should be
located about y+<5 for pressure distributions, y+<2 to obtain
reasonable skin friction values, and y+<0.5 for heat transfer. The grid
stretching normal to the wall should not exceed 1.3. Improved heat
transfer results can be obtained by using a constant spacing for the
first three cells off the wall.
2. In order to reduce the probability of finding a second peak well off the
wall, it is usually good to limit the number of points over which the F
function is calculated.
3. Care should be taken not to divide viscous regions such as boundary
layers when dividing the computational domain for blocked or chimera
applications since the entire velocity profile is required to properly
define the Fmax and Udiff quantities.
![](asset-3e.png)
5-1
5.0 One-Equation Turbulence Models
5.1 Theory
The general form for all transport turbulence models is similar to the form of a
species equation in reacting flows
()
)()(
1ZDZP
x
Z
xx
Z
U
t
Z
i
t
ji
i−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
+
∂
∂
=
∂
∂
+
∂
∂
νν
σ
(5.1)
where Z is a general turbulence variable. The terms on the left hand side of Eq.
5.1 represent the convective transport of Z. The first term on the right hand side
of Eq. 5.1 is the diffusion of Z assuming a gradient-diffusion model. The
functions P(Z) and D(Z) on the right hand side of Eq. 5.1 represent the production
and destruction of Z. The production function for turbulence models are usually
functions of the eddy viscosity and the fluid strain or vorticity. A functional
relationship between Z and the eddy viscosity is required to complete the model.
Several one-equation transport models for turbulence have appeared over the
years. The most popular one-equation models have been the Baldwin-Barth1
model and the Spalart-Allmaras model2. The Spalart-Allmaras model will be
discussed in some detail here.
5.2 Spalart-Allmaras Model
The Spalart-Allmaras turbulence model was derived using empirical
relationships, dimensional analysis, and Galilean invariance. The goal was to
produce a turbulent transport model that was fast, numerically stable, and
reasonably accurate for both shear layers and boundary layers. The model uses
a turbulence variable
ν
~
that has the dimensions of viscosity. The model can be
written as
()()()
[]
() ()
νννννν
σ
ν
ν
~~~~~
1
~
~
2
2DPC
x
U
tb
i
i−+∇+∇+⋅∇=
∂
∂
+
∂
∂
(5.2)
where the production term is given by
()
⎟
⎠
⎞
⎜
⎝
⎛+Ω= 2
22
1
~
~~ vb f
d
CP
κ
ν
νν
(5.3)
and the dissipation is given by
()
2
1
~
~⎟
⎠
⎞
⎜
⎝
⎛
=d
fCD ww
ν
ν
(5.4)
![](asset-3f.png)
5-2
Here Ω is the magnitude of the vorticity, d is the distance to the nearest wall, and
fv2 and fw are given by
1
21
1
v
vf
f
χ
χ
+
−= (5.5)
6
1
63
6
63
1⎥
⎦
⎤
⎢
⎣
⎡
+
+
=
w
w
wCg
C
gf (5.6)
and the remaining functions are given by
3
1
3
3
1
v
vC
f+
=
χ
χ
(5.7)
ν
ν
χ
~
= (5.8)
(
)
rrCrg w−+= 6
2 (5.9)
()
2
22 ~
~
v
fd
r
νκ
ν
+Ω
= (5.10)
The remaining constants are given in Table 5.1.
Cb1 C
b2
σ κ
Cw1 C
w2 C
ν
1
0.1355 0.622 2/3 0.41
σκ
/)1(/ 2
2
1bb CC ++ 0.3 7.1
Table 5.1 Coefficients for the Spalart-Allmaras model.
The eddy viscosity is defined as
1
~
vt f
ν
ν
= (5.11)
The original model also includes a boundary layer transition model. The trip
function can be used for simple geometries, but is difficult to apply for complex
geometries and is usually omitted. The turbulent transition location can also be
fixed by setting the production term to zero in the region of the flow upstream of
the desired transition location.
The wall boundary condition is 0
~
=
ν
. At free stream boundaries,
ν
~ is set to a
small number, generally
ν
/10. The turbulence variable
ν
~ is usually limited to be
greater than the free stream value in the field so that it will not go negative. The
model will self-initialize from this small value of
ν
~
in the freestream. The
![](asset-40.png)
5-3
numerical implementation for implicit algorithms is discussed in detail in Ref. 2.
Generally the convective terms are treated with upwind differencing. The
diffusion term is treated with central differencing. The full source term jacobians
(
ν
ν
ν
ν
~)
~
(
,
~)
~
(
∂
∂
∂
∂DP ) should be used in implicit formulations. When the r function
becomes large, the fw function becomes a constant, so r can be clipped at a
value of 10 to avoid numerical overflow problems.
The Spalart-Allmaras model suffers from a deficiency common to transport type
models in that it fails to predict the reduction in shear layer growth rate with
increasing jet Mach number. The model will predict the incompressible shear
layer growth rate regardless of the jet Mach number. Paciorri and Sabetta3
suggest a compressibility correction for this model.
5.3 Rotation/Streamline Curvature Corrections
Conventional linear eddy viscosity models have difficulty predicting flows with
large system rotation or streamline curvature. The Spalart-Allmaras model falls
into this category. The vorticity reaches a local maximum in the core of a vortex
causing the eddy viscosity to increase rapidly. This leads to excessive
dissipation of the vortex core. To remedy this problem the production term (Eq.
5.3) be modified by multiplying by a correction factor4.
() ()
*
~
~~ 12
22
1rFf
d
CP rvb ⎟
⎠
⎞
⎜
⎝
⎛+Ω=
κ
ν
νν
(5.12)
where Fr1(r*) is a correction for rotational flow given
()
[]
1*,0min1*)(
1−+= rCrF vorr (5.13)
and
ε
ε
+Ω
+
=S
r* (5.14)
where |S| is the magnitude of the strain. The strain is defined as
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
u
x
u
S5.0
(5.15)
The magnitude of the strain is given by
ijij SSS 2
2=
(5.16)
![](asset-41.png)
5-4
The
ε
in Eq. (8) is a small number used as a threshold value so that the
turbulence model returns to its baseline form in areas of low vorticity and strain
such as near free stream boundaries. The
ε
value can be defined as
ref
ref
L
U
5.0=
ε
(5.17)
Where Uref and Lref are the reference velocity and length scale for the calculation
respectively.
This correction form implies that curvature effects are a function of the ratio of the
local strain magnitude and vorticity magnitude. The strain goes to zero in the
core of a vortex while the vorticity reaches a local maximum. Thus Eq. (5.13)
converts the production term of the SA model into a dissipation term in the vortex
core region if Cvor is greater than one. Ref. 4 recommended Cvor=2. The
modified form of the production term has almost no effect in boundary layers and
shear layers since the magnitude of the strain and the magnitude of the vorticity
are almost equal for these flows and the correction term goes to zero. This
correction is generally referred to as the approximate Spalart-Allmaras rotational
correction (ASARC).
Shur, et al.5 introduced a correction to the Spalart-Allmaras turbulence model for
rotating and curved flows. The correction is applied to the production term in a
manner similar to the ASARC correction. The correction has the form
()( ) ()
[]
12
1
311 ~
tan1
*1
*2
1)
~
*, rrrrr crcc
r
r
crrF −−
+
+= − (5.18)
where
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛Ω++
Ω
=c
minjmnjnimn
ijjkik SS
Dt
DS
D
S
r
εε
4
2
~ (5.19)
The remaining term D is define as
(
)
222 5.0 Ω+= SD (5.20)
The recommended values for the constants are cr1=1.0, cr2=12.0 and cr3=1.0.
The Lagrangian derivative of the strain tensor in Eq. (5.19) is included to obtain
the proper value for the curvature correction. This term makes the correction
expensive and cumbersome to implement. Some implementations of this
correction ignore the contribution of the time derivative of the Lagrangian
derivative and the reference frame rotation to simplify the coding. This is referred
as the SARC correction.
![](asset-42.png)
5-5
An example of the standard SA, SARC, and ASARC models applied to flow in a
u-duct6 is shown in Fig. 5.1. The distance along the wall is denoted by s. The
standard SA model under-predicts the pressure coefficient and the skin friction at
the exit of the u-bend. Both the SARC and ASARC models improve the
comparison with data.
Figure 5.1. Pressure coefficient and skin friction on the outer wall of a u-duct..
The standard SA model also overdamps the flow in the core of a vortex. This
tends to cause vortices to be damped out prematurely in vortex dominated flows
such as wing or rotor tip vortices. The production term for the SA model (Eq. 5.3)
is based on the vorticity magnitude. The vorticity magnitude reaches a local
maximum in a vortex core, hence the model tends to produce much more eddy
viscosity than is required to properly simulate the flow in the vortex core.
The SARC and ASARC correction terms are negligible in simple shear layers
and boundary layers and do not affect the model since the magnitude of the
strain is approximately equal to the magnitude of the vorticity. The correction
term transforms the production term into a dissipation term in the core of a vortex
where the magnitude of the strain is zero. This additional dissipation tends to
reduce the eddy viscosity in the core of the vortex and significantly improves the
simulation of vortex dominated flows. The wing tip vortex on a NACA 0012 wing4
tip at an angle-of-attack of 10o is shown in Fig. 5.2 for both the original and
modified SA production terms. The tangential velocity predicted by the modified
SA models is in much better agreement with the data. Fig. 5.2b also includes
results for the hybrid RANS/LES Detached Eddy Simulation (DES) models that
are discussed in Chapter 9.
![](asset-43.png)
5-6
Figure 5.2a Axial velocity for a NACA 0012 wing tip vortex.
Fig. 5.2b Vortex tangential velocity.
5.4 Grid Sensitivity for a Flat Plate with Adiabatic Walls
The initial wall spacing of the computational grid and the grid-stretching ratio can
affect the accuracy of the Spalart-Allmaras model. Figure 5.3 shows the
sensitivity of the skin friction to initial wall spacing for a flat plate. The grid-
stretching ratio was 1.2 for all these cases. The plots include the theoretical skin
friction curves of White and of Spalding.
Tangential Velocity
Vortex Radius
![](asset-44.png)
5-7
The boundary layer is seen to become fully turbulent around a length Reynolds
number (Rex) of 1x106. The results for y+=0.2 and y+=1 are virtually identical
indicating a grid independent solution. The y+=5 solution show some small
divergence from the y+=1 solution at the lower length Reynolds numbers while
the y+=10 solution shows large differences from the other solutions.
Predicted velocity profiles for the flat plate boundary layer for various initial wall
grid point spacings are shown in Fig. 5.4. The velocity profile shows little effect
of the initial spacing for all but the y+=10 profile. All of the profiles but the y+=10
profile are in good agreement with the theoretical profile from Spalding. Note
that the theoretical profile does not include the law-of-the wake (see Chapter 3),
and hence the predicted profiles diverge from the theoretical profile in the wake
region of the boundary layer. The predicted eddy viscosity for various initial wall
spacings is shown in Fig. 5.5.
Figure 5.3 Flat plate skin friction predictions for the Spalart-Allmaras turbulence
model for varying initial wall grid point spacings.
![](asset-45.png)
5-8
Figure 5.4 Flat plate boundary layer profiles predicted by the Spalart-Allmaras
turbulence model for varying initial wall grid point spacings.
Figure 5.5 Eddy viscosity distribution predicted by the Spalart-Allmaras
turbulence model for varying grid initial wall spacings.
Here again it is seen that the y+=0.2 and the y+=1.0 results are almost identical.
The y+=5 and y+=10 results show the solutions are no longer grid independent at
larger wall spacings.
The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig.
5.6. All of these solutions used an initial wall spacing of y+=1.
![](asset-46.png)
5-9
There seems to be very little effect of grid stretching for these cases indicating
that the initial wall spacing is the more critical parameter for skin friction
predictions for flat plates with the Spalart-Allmaras turbulence model. This is
also the case for the velocity profile, as seen in Fig. 5.7. The eddy viscosity
does change as the stretching ratio increases as shown in Fig. 5.8. It is
interesting to note that a wide range of eddy viscosity distributions have little
effect on skin friction and the velocity profile for a flat plate boundary layer.
Figure 5.6 The effect of grid stretching ratio on the skin friction for a flat plate
boundary layer using the Spalart-Allmaras turbulence model.
Figure 5.7 The effect of grid stretching ratio on the velocity profile for a flat plate
boundary layer using the Spalart-Allmaras turbulence model.
![](asset-47.png)
5-10
Figure 5.8 The effect of grid stretching ratio on the eddy viscosity distribution for
a flat plate boundary layer using the Spalart-Allmaras turbulence model.
5.5 Grid Sensitivity for and Axisymmetric Bump
A second example of the grid sensitivity of the Spalart-Allmaras turbulence model
that includes a pressure gradient is the NASA Ames transonic axisymmetric
bump experiment described in the previous chapter. The effect of initial grid
spacing on the pressure coefficient distribution along the bump is shown in Fig.
5.9. The stretching ratio was 1.2 for these cases. The pressure coefficient
seems to be relatively insensitive to the initial grid spacing for y+<10. The y+=20
result shows a significant difference from the other results in regions where a
pressure gradient is present. The velocity distribution at the aft junction of the
bump and the cylinder (x/c=1) is shown in Fig. 5.10. The y+=20 solution predict
a larger velocity in the reverse flow region than the other solutions. Grid
stretching effects on the pressure coefficient distribution along the bump is
shown in Fig. 5.11. The initial grid spacing was 1.2 for these cases. The
pressure distribution coefficient changes slightly as the grid-stretching ratio is
increased to 1.5. The solution in the separated region differs greatly for a grid-
stretching ration of 2.0. The effect on the velocity distribution at x/c=1 is shown in
Fig. 5.12. As with increasing initial grid spacing, increasing the grid spacing
increases the size and the magnitude of the separated flow region.
![](asset-48.png)
5-11
Figure 5.9 The effect of initial grid wall spacing on the pressure coefficient for the
axisymmetric bump using the Spalart-Allmaras turbulence model.
Figure 5.10 The effect of initial grid wall spacing on the velocity distribution at
x/c=1 for the axisymmetric bump using the Spalart-Allmaras turbulence model.
![](asset-49.png)
5-12
Figure 5.11 The effect of grid stretching on the pressure coefficient for the
axisymmetric bump using the Spalart-Allmaras turbulence model.
Figure 5.12 The effect of grid stretching on the velocity profile at x/c=1 for the
axisymmetric bump using the Spalart-Allmaras turbulence model.
As was seen with the Baldwin-Lomax model, grid stretching is a critical
parameter for this turbulence model when adverse pressure gradients are
present. Care should be taken to keep the grid-stretching ratio between 1.2 and
1.3.
![](asset-4a.png)
5-13
5.6 Grid Sensitivity for a Flat Plate with Heat Transfer
Calculating heat transfer accurately can be more difficult than predicting skin
friction. This can be seen in the subsonic flat plate example when the wall
temperature is specified to be 1.5 times the free-stream temperature. The
sensitivity of the skin friction and heat transfer result with varying initial grid wall
spacing is shown in Fig. 5.13 and Fig. 5.14. The grid stretching ratio was fixed at
1.2 for these results. Both the skin friction and heat transfer seem to be relatively
insensitive to the wall spacing for initial wall spacings less than y+=1.
Figure 5.13 The effect of wall spacing on the skin friction on a flat plate with heat
transfer using the Spalart-Allmaras turbulence model.
The effect of wall spacing on the velocity and temperature profiles for a length
Reynolds number (Rex) of 1.0x107 is shown in Fig. 5.15 and 5.16 respectively.
The profiles are relatively insensitive to wall spacing for initial wall spacings less
than y+=1.0.
![](asset-4b.png)
5-14
Figure 5.14 The effect of wall spacing on the heat transfer (Stanton number) on
a flat plate using the Spalart-Allmaras turbulence model.
Fig. 5.15 The effect of wall spacing on the velocity profile for a flat plate with
heat transfer using the Spalart-Allmaras turbulence model.
![](asset-4c.png)
5-15
Fig. 5.16 The effect of wall spacing on the temperature profile for a flat plate with
heat transfer using the Spalart-Allmaras turbulence model.
Grid stretching effects on skin friction and heat transfer predictions are shown in
Fig. 5.17 and 5.18. The effect on velocity and temperature profiles for a length
Reynolds number (Rex) of 1.0x107 are shown in Fig. 5.19 and 5.20. The initial
wall spacing was held at y+=0.1 for these calculations. The results reach a grid
independent result for stretching ratios of less than 1.3.
Figure 5.17 The effect of grid stretching on the skin friction on a flat plate with
heat transfer using the Spalart-Allmaras turbulence model.
![](asset-4d.png)
5-16
Figure 5.18 The effect of grid stretching on heat transfer (Stanton number) on a
flat plate using the Spalart-Allmaras turbulence model.
Figure 5.19 The effect of grid stretching on the velocity profile on a flat plate with
heat transfer using the Spalart-Allmaras turbulence model.
![](asset-4e.png)
5-17
Figure 5.20 The effect of grid stretching on the temperature profile on a flat plate
with heat transfer using the Spalart-Allmaras turbulence model.
5.7 Grid Sensitivity for a Nozzle with Heat Transfer
Flow through a supersonic nozzle with a constant temperature wall can serve as
a test case for evaluating the performance of the turbulence model in the
presence of strong pressure gradients. Details of the geometry and boundary
conditions for the converging-diverging supersonic nozzle are given in Chapter 4.
High-pressure air was heated by the internal combustion of methanol and flowed
along a cooled constant area duct before entering the nozzle. The gas could be
treated as a calorically perfect gas with a ratio-of-specific heats (γ) of 1.345. The
nozzle exit Mach number was 2.5. The molecular viscosity and thermal
conductivity were assumed to vary according to Sutherland’s law.
The grid initial wall spacing was varied and the grid stretching ratio was held at
1.2 in the boundary layer. A uniform grid was used in the nozzle core.
Predicted wall pressure distribution results for varying initial wall spacings are
shown in Fig. 5.21. The pressure distribution is seen to be insensitive to the
initial wall spacing. Predicted wall heat transfer is shown in Fig. 5.22. The
results are somewhat sensitive to the initial wall spacing for values of y+ less than
one. The predicted results diverge rapidly from the data for wall spacing greater
than y+=1.
![](asset-4f.png)
5-18
Figure 5.21 The effect of wall spacing on the pressure distribution for a
supersonic nozzle with heat transfer using the Spalart-Allmaras turbulence
model.
Figure 5.22 The effect of wall spacing on the wall heat transfer for a supersonic
nozzle using the Spalart-Allmaras turbulence model.
The grid stretching ratio was varied in the boundary layer while the initial grid
spacing held at y+=0.5. A uniform grid was used in the nozzle core. Predicted
wall pressure distribution results for varying initial wall spacings are shown in Fig.
5.23. The pressure distribution is seen to be insensitive to the grid stretching
![](asset-50.png)
5-19
ratio. Predicted wall heat transfer is shown in Fig. 5.24. The results are
insensitive to the grid stretching ratios less than 1.3. The predicted results
diverge slightly from the data for grid stretching ratios greater than 1.3.
Figure 5.23 The effect of grid stretching ratio on the pressure distribution for a
supersonic nozzle with heat transfer using the Spalart-Allmaras turbulence
model.
Figure 5.24 The effect of grid stretching ratio on the wall heat transfer for a
supersonic nozzle using the Spalart-Allmaras turbulence model.
![](asset-51.png)
5-20
5.8 Summary
The Spalart-Allmaras turbulence model is a very stable and generally reasonably
accurate model for a wide range of turbulent flows. The model is relatively easy
to implement in both structured and unstructured Navier-Stokes codes. One
drawback is the model requires the calculation of the distance to the nearest wall
for all field points. This can be an expensive computation, especially for
unstructured grid codes. The model has been used extensively for three-
dimensional geometries and is well documented. The model has been used with
some success for some unsteady flows. The example given here are for
subsonic and supersonic flows. An excellent validation source for hypersonic
flows can be found in Ref. 7.
Chapter 5 References:
1. Baldwin B. S. and Barth, T. J.,”A One-Equation Turbulence Transport Model
for High Reynolds Number Wall-Bounded Flows,” AIAA-91-0610, Jan. 1991.
2. Spalart, P. R. and Allmaras, S. R., “A One-Equation Turbulence Model for
Aerodynamic Flows,” AIAA-92-0439, Jan. 1992.
3. Paciorri, R. and Sabetta, F., “Compressibility Correction for the Spalart-
Allmaras Model in Free-Shear Flows,” Journal of Spacecraft and Rockets, Vol.
40, No. 3, May-June 2003.
4. Dacles-Mariani, J., Zilliac, G., Chow, J., and Bradshaw, P.,
“Numerical/Experimental Study of a Wingtip Vortex in the Near Field,” AIAA
Journal, Vol. 33, No. 9, September 1995, pp. 1561-1568.
Spalart-Allmaras Application Hints
1. The first point off the wall should be located about y+=1 to obtain
reasonable skin friction values and about y+=0.5 for heat transfer.
The grid stretching normal to the wall should not exceed 1.3.
2. The eddy viscosity should be limited so that it will not run away in
some complex applications. Generally a limit of νt/ν=200,000 is
acceptable.
3. Care should be taken not to divide viscous regions such as boundary
layers when dividing the computational domain for blocked or overset
applications since the model requires the distance from the nearest
wall.
4. This model tends to smear out three-dimensional vortical flows.
Rotation and curvature corrections (SARC and ASARC) can
significantly improve the results.
5. The model can overdamp some unsteady flows.
6. The model contains no corrections for compressibility and will
overpredict the growth rate of high speed shear layers.
5-21
5. Shur, M. L., Strelets, M. K., Travin, A. K., and Spalart, P. R., “Turbulence
Modeling in Rotating and Curved Channels: Assessing the Spalart-Shur
Correction,” AIAA Journal, Vol. 38, No. 5, May 2000, pp.784-792.
6. Nichols, R., “Algorithm and Turbulence Model Requirements for Simulating
Vortical Flows,” AIAA-2008-0337, Jan. 2008.
7. Roy, J. and Blottner, F., “Methodology for Turbulence Model Validation:
Application to Hypersonic Flows,” Journal of Spacecraft and Rockets, Vol. 40,
No. 3, May-June 2003.
![](asset-53.png)
6-1
6.0 Two-Equation Turbulence Models
Two-equation turbulence models have been used for almost fifty years. There
are a tremendous number of these models in the literature. Most of these
models solve a transport equation for turbulent kinetic energy (k) and a second
transport equation that allows a turbulent length scale to be defined. The most
common forms of the second transport equation solve for turbulent dissipation (
ε
)
or turbulent specific dissipation (
ω
). Some of these two-equation models are
valid down to the wall (low Reynolds number models) and some are only valid
outside the inner region of the boundary layer (high Reynolds number models).
All of these models have their own strengths and weaknesses. Several of the
currently popular models are described in this chapter.
6.1 Theory
A transport equation for turbulent kinetic energy can be derived from Navier-
Stokes equations and has the form
j
i
j
i
j
i
ji
jj
i
ji
ix
u
x
u
x
U
uu
p
uu
u
xx
k
U
t
k
∂
∂
∂
∂
−
∂
∂
−
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
−=
∂
∂
+
∂
∂''
''
'
''
'
2
ν
ρ
(6.1)
The left hand side of Eq. 6.1 contains the temporal change and convective terms
for the turbulent kinetic energy. The first term on the right hand side of Eq. 6.1 is
the diffusive transport of the turbulent stresses (Dk). The second term on the
right hand side of Eq. 6.1 is the production of turbulent kinetic energy (Pk). The
last term on the right hand side of Eq. 6.1 represents the dissipation of turbulent
kinetic energy (
ε
). All of the terms on the right hand side must be modeled.
Applying the Boussinesq approximation (Eq. 2.22) the production term becomes
i
i
ijijt
i
i
j
i
ijt
i
i
j
i
i
j
j
i
tk x
U
kSS
x
U
k
x
U
S
x
U
k
x
U
x
U
x
U
P∂
∂
−=
∂
∂
−
∂
∂
=
∂
∂
−
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=3
2
3
2
3
2
ννν
(6.2)
The last term in Eq. 6.2 is zero for incompressible flow since the divergence of
the velocity (
i
i
x
U
∂
∂) is zero. This term is also neglected in many compressible flow
applications since it can cause numerical difficulties near strong shocks and
expansions. Sij is the strain rate defined as
2
2
1⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
u
x
u
S (6.3)
Kato and Launder1 suggested redefining the production term as
![](asset-54.png)
6-2
ijijtk SP Ω=
υ
(6.4)
Where Ωij is the vorticity defined as
2
2
1⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
=Ω
i
j
j
i
ij x
u
x
u (6.5)
The purpose of this modification is to reduce the tendency of many two-
equation models to overpredict the turbulent production in regions with large
normal strain, i.e. regions with strong acceleration or deceleration. This
correction is often used in stagnation regions such as at the leading edge of a
wing or blade.
The diffusion term can be modeled using the gradient-diffusion approximation as
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
=
ik
t
i
kx
k
x
D
σ
ν
ν
(6.6)
where
σ
k can be thought of as a turbulent Schmidt number.
Relationships for the eddy viscosity (
ν
t) and the turbulent dissipation (
ε
) are
required for closure. The square root of the turbulent kinetic energy can be used
to represent a velocity scale for the large-scale turbulent motion. Using this in an
eddy viscosity relationship yields
Lk
t≈
ν
(6.7)
where L is a turbulent length scale. A relationship for the turbulent length scale
of the large scale eddies is given in Eq. 1.13. Using this definition of the
turbulent length scale in Eq. 6.7 produces
ε
ν
μ
2
k
C
t= (6.8)
where C
μ
is a constant. The system of equations would be closed with a
relationship for the turbulent dissipation (
ε
).
The two features that distinguish different two equation models are the treatment
of the wall and the form of the turbulent dissipation equation. It is possible to
derive a transport equation for turbulent dissipation from the Navier-Stokes
equations, but the equation contains fluctuating correlation terms that cannot be
![](asset-55.png)
6-3
easily modeled. Turbulent dissipation cannot be measured directly
experimentally, but must be inferred from measurements of other correlations.
Because of the difficulty in deriving a turbulent dissipation equation based on
physics, turbulence modelers have traditionally resorted to using a transport
equation with empirically determined production and dissipation terms. The two
most popular equations used to close the turbulence model are based on the
dissipation
ε
2 and a turbulence variable
ω
defined by Wilcox3 as
kC
μ
ε
ω
= (6.9)
The eddy viscosity for the Wilcox formulation is
ω
ν
k
t= (6.10)
6.2 Traditional k-
ε
Models
The standard form of the
ε
transport equation is
R
k
fCP
k
fC
xxx
U
tk
i
t
ii
i+−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
=
∂
∂
+
∂
∂2
2211
εεε
σ
ν
ν
εε
εεεε
ε
(6.11)
The production and dissipation terms of the dissipation equation are formed from
the production and dissipation terms of the turbulent kinetic energy equation
scaled by
ε
/k and multiplied by empirically determined constants and wall
damping functions (C
ε
1f
ε
1 and C
ε
2 f
ε
2). An additional damping term must be
included for the eddy viscosity in the k-
ε
equations for applications near walls so
that k and
ε
will have the proper behavior in the near wall region.
The constants C
ε
1, C
ε
2, and C
μ
are empirically determined from comparisons with
experimental data. In grid-generated turbulence the diffusion term and the
production term in Eq. 6.1 and Eq. 6.11 are negligible and C
ε
2 can be determined
from the measured rate of decay of the turbulent kinetic energy k behind the grid.
The value of C
ε
2 was found to lie between 1.8 and 2.0. For shear-layers in local
equilibrium (production=dissipation), Eqs. 6.1 and 6.8 can be combined to
produce
2
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛′′
=k
vu
C
μ
(6.12)
Measurements in equilibrium flows suggest that C
μ
~0.09. In the near-wall region
of the logarithmic region of a boundary layer turbulent production and dissipation
![](asset-56.png)
6-4
are approximately equal and the convection of
ε
is negligible. Under these
assumptions Eq. 6.11 reduces to
με
εε σ
κ
C
CC 2
21 −= (6.13)
Note that all of these empirical constants are defined for near equilibrium
conditions, and thus are not truly valid when production and dissipation are not of
the same size. This is often a source of numerical difficulties when solving these
models and leads to the use of ad hoc relationships to control the behavior of the
turbulence model for cases when production is much greater than dissipation.
The general form of the eddy viscosity relationship for low Reynolds number
models is
ε
νμμ
2
k
fC
t= (6.14)
Many k-
ε
turbulence transport models exist. The coefficients of some of the more
popular are shown in Table 6.1. The wall damping functions are shown in Table
6.2.
Model C
ε
1 C
ε
2
σ
k
σ
ε
C
μ
κ
Standard 1.44 1.92 1.0 1.3 0.09 0.43
Jones-
Launder2 1.45 1.92 1.0 1.3 0.09 0.43
Chien
4
1.35 1.80 1.0 1.3 0.09 0.42
Speziale
5
1.47 1.83 1.0 1.56 0.09 0.41
Table 6.1 Coefficients for low Reynolds number k-
ε
models.
Model f
ε
1 f
ε
2
f
μ
Standard 1 1 1
Jones-
Launder2
2
2
1⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+y
k
ε
ν
()
2
2
2
2
22
exp3.1 ⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+−− y
u
Rt
t
ε
νν
()
⎥
⎦
⎤
⎢
⎣
⎡
+
−2
02.01
4.3
exp
t
R
Chien
4
ε
ν
2
2
1y
k
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−+
⎥
⎦
⎤
⎢
⎣
⎡⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−−
+
2
exp
2
36
exp222.1
2
2
2
y
y
k
R
Ct
ε
ν
ε
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−−
+
87
exp1 y
![](asset-57.png)
6-5
Speziale
5
1
2
2
9.4
exp1*
36
exp222.1
⎥
⎦
⎤
⎢
⎣
⎡⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−−
⎥
⎦
⎤
⎢
⎣
⎡⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−−
+
y
Rt
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
+
70
tanh
45.3
1y
Rt
Table 6.2 Wall damping functions for low Reynolds number k-
ε
models.
The last column of Table 6.1 also contains the von Karman constant (
κ
)
calculated from Eq. 6.10. All of the models shown here produce
κ
values slightly
above the usually accepted value of 0.4-0.41 with the exception of the Speziale
model.
The wall boundary condition for the turbulent kinetic energy is
0=
wall
k (6.15)
Since the turbulent kinetic energy is zero at the wall, then the eddy viscosity is
also zero at the wall. A wall boundary condition for the dissipation (
ε
) can be
found from the turbulent kinetic energy equation
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
=
ii
wall x
k
x
νε
(6.16)
This boundary condition can produce numerical stability problems, so Speziale4
suggested the boundary condition be replaced with the form
2
2⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
=
i
wall x
k
νε
(6.17)
The turbulent kinetic energy is usually assigned to the free stream turbulence
level at free stream boundaries if it is known. If the free stream turbulence level
is not known, then the turbulent kinetic energy is set to a small number at the free
stream boundary. The eddy viscosity is also set to a small number at free stream
boundaries.
1.0<
∞
∞
ν
ν
t (6.18)
![](asset-58.png)
6-6
Once the free stream turbulent kinetic energy and eddy viscosity are known, the
free stream value of turbulent dissipation may be calculated from the definition of
eddy viscosity.
6.3 k-
ω
Model
Using Wilcox’s formulation, a transport equation for
ω
can be written as
2
2
βω
ω
σ
κ
β
ω
σ
ν
ν
ωω
μϖ
ω
−
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
=
∂
∂
+
∂
∂
k
i
t
ii
iP
k
C
xxx
U
t (6.19)
The
ω
transport equation is essentially the
ε
transport equation (Eq. 6.11) with
the addition of a cross diffusion term of the form
Cross Diffusion =
ii
t
x
k
x∂
∂
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
ε
σ
ν
ν
ε
(6.20)
The k-
ω
models are derived for wall bounded flows and require no additional wall
damping terms when used in boundary layer flows. The coefficients for Eq. 6.19
are given in Table 6.3.
β σ
ω
C
μ
κ
0.075 2 0.09 0.41
Table 6.3 Coefficients for the k-ω model.
The wall boundary condition for the turbulent kinetic energy is also given by Eq.
6.15. The definition of
ω
(Eq. 6.9) suggests that
ω
should go to infinity at the wall
since the turbulent kinetic energy (k) is going to zero. Wilcox, based on
asymptotic arguments, suggested that
ω
be given the value
()
2
60
yC
wall Δ
=
μ
ν
ω
(6.20)
at the wall. Here
Δ
y is the normal distance to the first point off the wall. The
following free stream boundary conditions are recommended
∞∞∞
∞
∞
∞
∞===
ωνω
ν
ν
t
tk
L
U,10,001.0 (6.21)
The k-
ω
model has been shown to be sensitive to the free stream value of
ω
chosen.
![](asset-59.png)
6-7
6.4 SST Model
In practice, the k-
ε
models are generally more accurate in shear type flows and
are well behaved in the far field. The k-
ω
models are more accurate and much
more numerically stable in the near wall region. Recognizing that each model
has its strength and weakness and that the forms of the equations are similar,
Menter6 suggested a blended model that is a k-
ω
model near the wall and
transitions to a k-
ε
model away from the wall. The model has the form
()
kCP
x
k
xx
k
U
t
kk
i
tk
ii
i
ωνσν
μ
−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
+
∂
∂
=
∂
∂
+
∂
∂ (6.22)
()
()
ii
k
i
t
ii
i
xx
k
F
P
kxxx
U
t
∂
∂
∂
∂
−
+−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
+
∂
∂
=
∂
∂
+
∂
∂
ω
ω
σ
βω
ω
γ
ω
νσν
ωω
ω
ω
2
11
2
(6.23)
The coefficients are blended forms of the two baseline models. The blending
function F1 is defined as
(
)
4
11 argtanh=F (6.24)
where
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
=2
2
2
1
4
,
500
,maxminarg yCD
k
y
yC
k
k
ω
ω
μ
ρσ
ω
ν
ω
(6.25)
Here y is the normal distance to the wall and CDk
ω
is the positive portion of the
cross-diffusion term
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
=−20
210,
1
2max
ii
kxx
k
CD
ω
ω
ρσ
ωω
(6.26)
The coefficients
σ
k,
σ
ω
,
γ
, and
β
in Eq. 6.23 are computed with the general form
()
2111 1
φ
φ
φ
FF −+= (6.27)
where the
φ
1 corresponds to coefficients from the k-
ω
model and
φ
2 corresponds
to coefficients of the k-
ε
model. The eddy viscosity is calculated from
()
21
1
,max Fa
ka
tΩ
=
ω
ν
(6.28)
![](asset-5a.png)
6-8
where F2 is given by
(
)
2
22 argtanh=F (6.29)
and
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
=
ω
ν
ω
μ
2
2
500
,
2
maxarg y
C
k (6.30)
Here Ω is the magnitude of the vorticity vector. The first term in the denominator
of Eq. 6.28 comes from the traditional definition of eddy viscosity for a k-ω model
(Eq. 6.10). The second term in the denominator represents an attempt to
improve the performance of the model for adverse pressure gradients and comes
from a one equation turbulence modeled developed by Bradshaw7 based on a
relationship for shear stress in a boundary layer similar to Eq. 6.12
kavu 1
'' =− (6.31)
The second term in the denominator of Eq. 6.28 reduces the eddy viscosity in the
logarithmic and wake region of the boundary layer when adverse pressure
gradients are present. The use of the shear stress relationship of Eq. 6.31 in the
definition of the eddy viscosity (Eq. 6.28) has given rise to calling this model the
shear stress transport (SST) model.
The coefficients used with the SST model are given Tables 6.4 and 6.5.
β σ
k
σ
ω
γ
Set 1
k-
ω
0.075 0.85 0.5
μ
ω
μ
κσ
β
C
C
2
1
1−
Set 2
k-
ε
()
1
2
−
εμ
CC 1 0.857 1
1
−
ε
C
Table 6.4 SST model functions.
C
ε
1 C
ε
2 C
μ
κ
a1
1.44 1.92 0.09 0.41 0.31
Table 6.5 SST model coefficients.
The boundary conditions for SST model are the same as the k-
ω
model. The
SST model is relatively insensitive to the free stream value of
ω
.
![](asset-5b.png)
6-9
6.5 RNG Model
One of the major criticisms of the k-
ε
model is that it is not derived from the
Navier-Stokes equations in any systematic fashion. Yakhot and Orszag8 applied
Renormalization Group (RNG) methods to derive the k-
ε
equations. In this
approach, an expansion is made about an equilibrium state with known Gaussian
statistics using the correspondence principle that the effects of the mean strains
can be represented by a random force. Bands of high wave numbers (small
scales) are systematically removed and space is rescaled in a manner analogous
to that used in phase transitions. The successive removal of larger scales
ultimately leads to a k-
ε
model of the form
ε
σ
ν
ν
−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
=
∂
∂
+
∂
∂
k
ik
t
ii
iP
x
k
xx
k
U
t
k (6.32)
R
k
CP
k
C
xxx
U
tk
i
t
ii
i+−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
=
∂
∂
+
∂
∂2
21
εεε
σ
ν
ν
εε
εε
ε
(6.33)
where R is defined as
()
()
k
P
Rk
ε
ζ
ζ
ζ
3
012.01
38.4/1
+
−
−= (6.34)
()
ε
ζ
k
SS ijij
2/1
2= (6.35)
The parameter
ζ
is the ratio of the turbulent to mean strain time scale. The
remaining constants are given in Table 6.6.
C
ε
1 C
ε
2 C
μ
σ
k
σ
e
1.42 1.68 0.085 0.719 0.719
Table 6.6 RNG model coefficients.
The RNG model is similar to the traditional k-
ε
model with the exception of the
additional term R and the lower value of C
ε
2. These two differences decrease
both the rate of production of turbulent kinetic energy and the rate of dissipation
of
ε
, leading to lower values of the eddy viscosity. In free shear flows,
ζ
(and
thus R) is zero. The model presented is the high Reynolds number form. The
additional terms required for the low Reynolds form can be found in Ref. 9.
![](asset-5c.png)
6-10
6.6 Numerical Implementation
The flux terms on the right hand side of Eqs. 6.1, 6.11, and 6.19 are usually
treated with upwind differences. The diffusion terms are generally treated with
central differences. The turbulent kinetic energy production term Pk is limited to
be less than 20
ε
in the turbulent kinetic energy source term to keep the model
from “running away” in complex flows. The maximum value of eddy viscosity is
also capped (normally < 2x105
ν
) to limit non-physical behavior of the model. The
turbulent quantities (k,
ε
, and
ω
) are limited to be positive in the field. The exact
source term Jacobian term can be used for the
ε
or
ω
equations for implicit
computations. Use of the exact source term for the Jacobian of the turbulent
kinetic energy equation usually results in numerical instability. It is
recommended that the turbulent kinetic energy source term Jacobian be replaced
with
TKE source term Jacobian = k
Pk
ε
+
−2 (6.36)
The two equations need not be solved fully coupled, and can be lagged from the
mean flow equations. For time accurate implicit solutions the turbulence
equations should be solved within a Newton iteration loop with the mean flow
equations to assure that the solution is locally converged at each time step.
6.7 Compressibility Correction for Shear Layers
The standard k-
ε
and k-
ω
turbulence models will over predict the mixing for high-
speed compressible shear layers. The models tend to predict the incompressible
growth rate for the shear layer rather than decreasing the growth rate as
compressible effects increase. A common method of correcting the k-
ε
model for
this effect is given in Ref. 10. The compressible turbulent kinetic energy equation
is written as
()
""dpP
x
k
xx
kU
t
k
ck
ik
t
ii
i++−+
⎥
⎦
⎤
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
=
∂
∂
+
∂
∂
εερρ
σ
μ
μ
ρ
ρ
(6.37)
where
εαε
2
1tc M= (6.38)
Mt is the turbulent Mach number defined as
RT
k
Mt
γ
2
= (6.39)
![](asset-5d.png)
6-11
The last term in Eq. 6.37 represents the pressure dilatation. This term is
modeled as
2
3
2
2
"" ttk MMPdp
ρεαρα
+−= (6.40)
Based on direct numerical simulation, the recommended values for the constants
in Eqs. 6.38 and 6.40 are
2.0,4.0,0.1 321 ===
α
α
α
(6.41)
The compressibility correction must also be included in the SST model since the
SST model reduces to the standard k-
ε
model away from the wall. A method for
incorporating the compressibility correction on the SST model is given in Ref. 11.
Further corrections for high temperature jet flows are described in Ref. 12.
The need for a compressibility correction can be seen in predictions of the
supersonic axisymmetric jet shear layer of Eggers13. The jet exit Mach number
was 2.22 and the jet exit static pressure was matched to the quiescent outer air.
Fig. 6.1 shows the predicted and measured axial velocity on the jet centerline.
Fig. 6.2 contains velocity comparisons at various downstream locations in the jet.
Figure 6.1. Centerline axial velocity for the Ref. 12 jet.
![](asset-5e.png)
6-12
Figure 6.2a. Velocity profile at X/R=25 for the Ref. 12 jet.
Figure 6.2b. Velocity profile at X/R=50 for the Ref. 12 jet.
![](asset-5f.png)
6-13
Figure 6.2c. Velocity profile at X/R=100 for the Ref. 12 jet.
Without the compressibility correction all of the turbulence models predict a faster
growth rate than is seen in the experimental data. The SST model with the
compressibility correction is seen to do a good job for accounting for the
compressibility effect.
6.8 Initializing Turbulence Values for a Given Profile Shape
There is really no unique way to construct turbulence variables for a given
velocity profile shape for an inflow boundary condition for transport type
turbulence models. Here are two simple approaches. Often a combination of
these approaches is required in real world applications.
One approach is to specify a turbulent intensity and a turbulent length scale for a
given profile. The turbulent intensity is specified as a percentage (a) of the local
time-averaged fluid velocity
Uau =
′ (6.42)
The turbulent intensity is normally less than twenty per cent of the mean flow
velocity. A turbulent length scale (L) must then be specified. Some typical
choices are shown in Table 6.7.
Flow Boundary
Layer Plane Jet Round Jet Plane
Wake Pipe
L 0.41δ 0.53δ 0.44δ 0.95δ 0.25D
![](asset-60.png)
6-14
Table 6.7 Typical values for the turbulent length scale L.
In Table 6.7, δ is defined as the distance between points where the velocity
differs from the free stream velocity by one percent of the maximum velocity
difference across the layer. For symmetrical flows and boundary layers, δ is the
distance from the symmetry plane or wall to the one percent point at the outer
edge. D is the internal pipe diameter.
The turbulent kinetic energy can be found from
()
2
2
2
3
2
3aUuk =
′
= (6.43)
The turbulent length scale is related to the eddy viscosity by the Kolmogorov-
Prandtl expression
LkC
t
ρμ
μ
= (6.44)
The traditional definition of eddy viscosity for two-equation turbulence models is
given by
ε
ρμ
μ
2
k
C
t= (6.45)
Thus the turbulent dissipation is defined as
L
k2/3
=
ε
(6.46)
Although this method is relatively simple to implement, it is not very accurate for
boundary layers and looses meaning for complex velocity profile shapes.
The second method takes advantage of the turbulence equilibrium assumption
(turbulent production = turbulent dissipation). This requires an estimate of the
eddy viscosity along the profile. The eddy viscosity profile can be obtained from
a simple algebraic turbulence model such as Baldwin-Lomax or a one-equation
turbulence model such as Spalart-Allmaras.
The turbulent dissipation is then given by
k
P=
ε
(6.47)
Finally, the turbulent kinetic energy is obtained from the definition of eddy
viscosity in Eq. (6.45)
![](asset-61.png)
6-15
ρ
εμ
μ
C
kt
= (6.48)
This method works well for general applications and can be useful for starting k-
ε
models. These models often have difficulties generating eddy viscosity when a
solution is started from a uniform free stream condition.
6.9 Rotation and curvature correction
A rotational and curvature correction for the SST model can also be formulated in
a manner similar to the SARC and ASARC corrections used in the SA one-
equation model. The correction has a much smaller effect because the SST
uses strain rather than vorticity in the definition of the production term. The
production term used in both the k and
ω
equations (Eq. 6.22 and 6.23) is
redefined as
()
rrFSSP rijijtk
~
*,
1
μ
= (6.49)
The corrections defined in Eq. (5.13) and Eq. (5.18) are used with Eq. (6.49) to
provide a general three-dimensional rotation and curvature correction for the SST
model in this study. The Eq. (5.13) and Eq. (5.18) corrections are limited
0.0<Fr1<1.25 for the two equation model application. The SST model with the
Eq. (5.13) correction is called the ASSTRC model and the SST model with the
Eq. (5.18) correction is named the SSTRC model.
An example of the standard SST, SSTRC, and ASSTRC models applied to flow
in a u-duct13 is shown in Fig. 6.3. The distance along the wall is denoted by s.
The standard SST model slightly under-predicts the pressure coefficient and the
skin friction at the exit of the u-bend. Both the SARC and ASARC models
improve the comparison with data.
Figure 6.3. Pressure coefficient and skin friction on the outer wall of a u-duct..
![](asset-62.png)
6-16
6.10 Grid Sensitivity for a Flat Plate with Adiabatic Walls
The initial wall spacing of the computational grid and the grid-stretching ratio can
affect the accuracy of the SST model. Figure 6.4 shows the sensitivity of the skin
friction to initial wall spacing for a flat plate. The grid-stretching ratio was 1.2 for
all these cases. The plots include the theoretical skin friction curves of White and
of Spalding.
Figure 6.4 Flat plate skin friction predictions for the SST turbulence model for
varying initial wall grid point spacings.
The boundary layer is seen to become fully turbulent around a length Reynolds
number (Rex) of 1x106. The results for y+=0.2 and y+=1 are virtually identical
indicating a grid independent solution. The y+=5 solution show some small
divergence from the y+=1 solution at the lower length Reynolds numbers while
the y+=10 solution shows large differences from the other solutions. The SST
model tends to be more sensitive to initial grid wall spacing than the Spalart-
Allmaras or Baldwin-Lomax models.
Predicted velocity profiles for the flat plate boundary layer for various initial wall
grid point spacings are shown in Fig. 6.5. The velocity profile shows little effect
of the initial spacing for all but the y+=10 profile. All of the profiles but the y+=10
profile are in good agreement with the theoretical profile from Spalding. Note
that the theoretical profile does not include the law-of-the wake (see Fig. 3.2),
and hence the predicted profiles diverge from the theoretical profile in the wake
region of the boundary layer. The predicted eddy viscosity for various initial wall
spacings is shown in Fig. 6.6.
![](asset-63.png)
6-17
Figure 6.5 Flat plate boundary layer profiles predicted by the SST turbulence
model for varying initial wall grid point spacings.
Figure 6.6 Eddy viscosity distribution predicted by the SST turbulence model for
varying grid initial wall spacings.
Here again it is seen that the y+=0.2 and the y+=1.0 results are almost identical.
The y+=5 and y+=10 results show the solutions are no longer grid independent at
larger wall spacings.
![](asset-64.png)
6-18
The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig.
6.7. All of these solutions used an initial wall spacing of y+=1.
Figure 6.7 The effect of grid stretching ratio on the skin friction for a flat plate
boundary layer using the SST turbulence model.
There seems to be very little effect of grid stretching for this case indicating that
the initial wall spacing is the more critical parameter for skin friction predictions
for flat plates with the SST turbulence model. This is also the case for the
velocity profile, as seen in Fig. 6.8. The eddy viscosity does change as the
stretching ratio increases as shown in Fig. 6.9. It is interesting to note that a
wide range of eddy viscosity distributions have little effect on skin friction and the
velocity profile for a flat plate boundary layer.
![](asset-65.png)
6-19
Figure 6.8 The effect of grid stretching ratio on the velocity profile for a flat plate
boundary layer using the SST turbulence model.
Figure 6.9 The effect of grid stretching ratio on the eddy viscosity distribution for
a flat plate boundary layer using the SST turbulence model.
6.10 Grid Sensitivity for an Axisymmetric Bump
A second example of the grid sensitivity of the SST turbulence model that
includes a pressure gradient is the NASA Ames transonic axisymmetric bump
experiment described in the Baldwin-Lomax Chapter 4. The effect of initial grid
![](asset-66.png)
6-20
spacing on the pressure coefficient distribution along the bump is shown in Fig.
6.10. The stretching ratio was 1.2 for these cases. The pressure coefficient
seems to be relatively insensitive to the initial grid spacing, with the y+=10 and
y+=20 curves being slightly displaced from the other curves. The velocity
distribution at the aft junction of the bump and the cylinder (x/c=1) is shown in
Fig. 6.11. The y+=20 solution predict a larger velocity in the reverse flow region
than the other solutions. Grid stretching effects on the pressure coefficient
distribution along the bump is shown in Fig. 6.12. The initial grid spacing was 1.2
for these cases. The pressure distribution coefficient changes slightly as the
grid-stretching ratio is increased to 1.5. The solution in the separated region
differs greatly for a grid-stretching ration of 2.0. The effect on the velocity
distribution at x/c=1 is shown in Fig. 6.13. As with increasing initial grid spacing,
increasing the grid spacing increases the size and the magnitude of the
separated flow region.
Figure 6.10 The effect of initial grid wall spacing on the pressure coefficient for
the axisymmetric bump using the SST turbulence model.
![](asset-67.png)
6-21
Figure 6.11 The effect of initial grid wall spacing on the velocity distribution at
x/c=1 for the axisymmetric bump using the SST turbulence model.
Figure 6.12 The effect of grid stretching on the pressure coefficient for the
axisymmetric bump using the SST turbulence model.
![](asset-68.png)
6-22
Figure 6.13 The effect of grid stretching on the velocity profile at x/c=1 for the
axisymmetric bump using the SST turbulence model.
As was seen with the Baldwin-Lomax model and the Spalart-Allmaras model,
grid stretching is a critical parameter for this turbulence model when adverse
pressure gradients are present. Care should be taken to keep the grid-stretching
ratio between 1.2 and 1.3.
6.12 Grid Sensitivity for a Flat Plate with Heat Transfer
Calculating heat transfer accurately can be more difficult than predicting skin
friction. This can be seen in the subsonic flat plate example when the wall
temperature is specified to be 1.5 times the free-stream temperature. The
sensitivity of the skin friction and heat transfer result with varying initial grid wall
spacing is shown in Fig. 6.14 and Fig. 6.15. The grid stretching ratio was fixed at
1.2 for these results. Both the skin friction and heat transfer seem to be relatively
insensitive to the wall spacing for wall spacings less than y+=0.1.
![](asset-69.png)
6-23
Figure 6.14 The effect of wall spacing on the skin friction on a flat plate with heat
transfer using the SST turbulence model.
Figure 6.15 The effect of wall spacing on the heat transfer (Stanton number) on
a flat plate using the SST turbulence model.
Profiles of velocity and temperature at a length Reynolds number (Rex) of 1x107
are shown in Fig. 6.16 and 6.17 respectively. The results are relatively
insensitive to the wall spacing for initial wall spacings of y+>5.
![](asset-6a.png)
6-24
Figure 6.16 The effect of wall spacing on the velocity profile on a flat plate with
heat transfer using the SST turbulence model.
Figure 6.17 The effect of wall spacing on the temperature profile on a flat plate
with heat transfer using the SST turbulence model.
Grid stretching effects on skin friction and heat transfer predictions are shown in
Fig. 6.18 and 6.19. The velocity and temperature profiles for a length Reynolds
number (Rex) of 1x107 are shown in Fig. 6.20 and 6.21. The initial wall spacing
was held at y+=0.1 for these calculations. The results reach a grid independent
result for stretching ratios less than 1.3.
![](asset-6b.png)
6-25
Figure 6.18 The effect of grid stretching on the skin friction on a flat plate with
heat transfer using the SST turbulence model.
Figure 6.19 The effect of grid stretching on heat transfer (Stanton number) on a
flat plate using the SST turbulence model.
![](asset-6c.png)
6-26
Figure 6.20 The effect of grid stretching on the velocity profile on a flat plate with
heat transfer using the SST turbulence model.
Figure 6.21 The effect of grid stretching on the temperature profile on a flat plate
with heat transfer using the SST turbulence model.
6.13 Grid Sensitivity for a Nozzle with Heat Transfer
Flow through a supersonic nozzle with a constant temperature wall can serve as
a test case for evaluating the performance of the turbulence model in the
presence of strong pressure gradients. Details of the geometry and boundary
![](asset-6d.png)
6-27
conditions for the converging-diverging supersonic nozzle are given in Chapter 4.
High-pressure air was heated by the internal combustion of methanol and flowed
along a cooled constant area duct before entering the nozzle. The gas could be
treated as a calorically perfect gas with a ratio-of-specific heats (γ) of 1.345. The
nozzle exit Mach number was 2.5. The molecular viscosity and thermal
conductivity were assumed to vary according to Sutherland’s law.
The grid initial wall spacing was varied and the grid stretching ratio was held at
1.2 in the boundary layer. A uniform grid was used in the nozzle core.
Predicted wall pressure distribution results for varying initial wall spacings are
shown in Fig. 6.22. The pressure distribution is seen to be insensitive to the
initial wall spacing. Predicted wall heat transfer is shown in Fig. 6.23. The
results are somewhat sensitive to the initial wall spacing for values of y+ less than
one. The predicted results diverge rapidly from the data for wall spacing greater
than y+=1.
Figure 6.22 The effect of wall spacing on the pressure distribution for a
supersonic nozzle with heat transfer using the SST turbulence model.
![](asset-6e.png)
6-28
Figure 6.23 The effect of wall spacing on the wall heat transfer for a supersonic
nozzle using the SST turbulence model.
The grid stretching ratio was varied in the boundary layer while the initial grid
spacing held at y+=0.5. A uniform grid was used in the nozzle core. Predicted
wall pressure distribution results for varying initial wall spacings are shown in Fig.
6.24. The pressure distribution is seen to be insensitive to the grid stretching
ratio. Predicted wall heat transfer is shown in Fig. 6.25. The results are
insensitive to the grid stretching ratios less than 1.3. The predicted results
diverge slightly from the data for grid stretching ratios greater than 1.3.
![](asset-6f.png)
6-29
Figure 6.24 The effect of grid stretching ratio on the pressure distribution for a
supersonic nozzle with heat transfer using the SST turbulence model.
Figure 6.25 The effect of grid stretching ratio on the wall heat transfer for a
supersonic nozzle using the SST turbulence model.
6.14 Summary
Two-equation models have been the industry standard for many years. These
models have been used for a wide variety of internal and external flows. Many
![](asset-70.png)
6-30
variants are available that have been tuned for specific classes of flows. These
models are more difficult to implement than are the one-equation models and can
often be less numerically stable. For these reasons many code developers have
avoided two-equation models. The examples given here are mainly for subsonic
and supersonic applications. An excellent summary of hypersonic validation of
these models can be found in Ref. 15.
Chapter 6 References:
1. Kato, M. and Launder, B. E., "The Modeling of Turbulent Flow Around
Stationary and Vibrating Square Cylinders", Proc. 9th Symposium on
Turbulent Shear Flows, Kyoto, August 1993, pp. 10.4.1-10.4.6.
2. Jones, W. P., and Launder, B. E., “The Calculation of Low-Reynolds
Number Phenomena with a Two-Equation Model of Turbulence,” Int. J. of
Heat and Mass Transfer, Vol. 16, 1973, pp. 1119-1130.
Two-Equation Model Application Hints
1. The first point off the wall should be located about y+=1 to obtain
reasonable skin friction values and about y+=0.5 for heat transfer.
The grid stretching normal to the wall should not exceed 1.3. A
constant spacing should be used for the first three cells off the wall
for heat transfer calculations.
2. The eddy viscosity should be limited so that it will not run away in
some complex applications. Generally a limit of
ν
t/
ν
=200,000 is
acceptable.
3. Care should be taken not to divide viscous regions such as
boundary layers when dividing the computational domain for
blocked or overset applications since the model requires the
distance from the nearest wall.
4. k-
ε
models are often difficult to start. It is usually easier to come off
the ground with another turbulence model such as Baldwin-Lomax
or Spalart-Allmaras and switch to the k-
ε
model after
ν
t/
ν
>1.0 in the
computational domain. This can be done by first calculating the
turbulent kinetic energy production Pk with Eq. 6.2. Assuming the
turbulence is in local equilibrium, then the dissipation can be found
from
ε
=Pk. The turbulent kinetic energy k can then be found from
the definition of eddy viscosity in Eq. 6.8.
5. This model can overdamp some unsteady flows.
6. Compressibility corrections should be included for high-speed shear
layer flows.
6-31
3. Wilcox, D. C., “Reassessment of the Scale-Determining Equation for
advanced Turbulence Models, AIAA Journal, Vol. 26, No. 11, 1988, pp.
1299-1310.
4. Chien, J. Y., “Predictions of Channel Boundary-Layer Flows with a Low-
Reynolds-Number Turbulence Model,” AIAA Journal, Vol. 20, Jan. 1982,
pp. 33-38.
5. Speziale, C., Abid, R., and Anderson, E., “A Critical Evaluation of Two-
Equation Models for Near Wall Turbulence,” ICASE Report No. 90-46,
June 1990.
6. Menter, F. R. and Rumsey, C. L., “Assessment of Two-Equation
Turbulence Models for Transonic Flows,” AIAA-94-2343, June 1994.
7. Bradshaw, P., Ferriss, D. H., and Atwell, N. P., “Calculation of Boundary
Layer Development Using the Turbulent Energy Equation,” Journal of
Fluid Mechanics, Vol. 28, 1967, pp. 593-616.
8. Yakhot, V., Orszag, S., Thangam, S., Gatski, T., and Speziale, C. G.,
“Development of Turbulence Models for Shear Flows by a Double
Expansion Technique,” Physics of Fluids A, Vol. 4, No. 7, 1992, pp. 1510-
1520.
9. Barber, T., Choi, D., Nedungadi, A., Orszag, S., Konstantinov, A., and
Staroselsky, I., “Advanced Turbulence Models for High-Mach Number
Wall Bounded Flows,” AIAA-98-0536, Jan. 1998.
10. Speziale, C., Abid, R., and Anderson, E., “A Critical Evaluation of Two-
Equation Models for Near Wall Turbulence,” ICASE Report No. 90-46,
June 1990.
11. Suzen, Y. B. and Hoffmann, K. A., “Investigation of Supersonic Jet
Exhaust Flow by One- and Two-Equation Turbulence Models,” AIAA-98-
0322, Jan. 1998.
12. Abdol-Hamid, K. Pao, S. Massey, S., and Elmiligui, A., “Temperature
Corrected Turbulence Model for High Temperature Jet Flow,” AIAA-2003-
4070, Jun. 2003.
13. Eggers, J. “Velocity Profiles and Eddy Viscosity Distributions Downstream
of a Mach 2.22 Nozzle Exhausting to Quiescent Air,” NASA-TN D-3601,
September 1966.
14. Nichols, R., “Algorithm and Turbulence Model Requirements for
Simulating Vortical Flows,” AIAA-2008-0337, Jan. 2008.
15. Roy, J. and Blottner, F., “Methodology for Turbulence Model Validation:
Application to Hypersonic Flows,” Journal of Spacecraft and Rockets, Vol.
40, No. 3, May-June 2003.
![](asset-72.png)
7-1
7.0 Reynolds and Algebraic Stress Models
7.1 Reynolds Stress Models
One means of obtaining closure for the Reynolds or Favre averaged Navier-
Stokes equations is to investigate higher moments of the equations themselves.
Transport equations for each of the Reynolds stresses can be derived from the
momentum equations. These equations are called the second moment
equations. The second moment equations can be obtained by multiplying the
momentum equations by Uj and then Reynolds averaging the resulting
equations. The incompressible form of the second moment equation is
()
()
k
j
k
i
i
j
j
i
ikjjkiji
k
kji
k
k
j
ki
k
i
kj
k
ji
k
ji
x
u
x
u
x
u
x
u
p
uu
p
uu
x
uuu
x
x
U
uu
x
U
uu
x
uu
U
t
uu
∂
∂
∂
∂
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
+
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡++
∂
∂
−
∂
∂
−
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
−=
∂
∂
+
∂
∂
'
'
'
'
'
''
'
'''''
''''
''''
2
ν
ρ
δδ
ρ
ν
(7.1)
The first term on the right hand side of Eq. 7.1 is the turbulent production term.
The second term on the right hand side is the diffusion of the Reynolds stresses.
The third term on the right side is called the pressure strain. The pressure strain
term redistributes energy among the turbulent stress components. The last term
on the right hand side is the turbulent dissipation. The balance of these source
terms is critical to the proper simulation of a turbulent flow. The second moment
equations produce a host of new unknowns and higher order correlations that
must be modeled in order to close the system. These higher order correlations
are difficult, and in many cases impossible, to determine from experimental
measurements. Direct Numerical Simulation (DNS) has been used to determine
these quantities in many cases. Unfortunately DNS is restricted to low Reynolds
numbers and simple geometries at the present time, so the data base for
determining these higher order correlations, whether experimentally or
numerically generated, is limited. Eq. 7.1 may be rewritten as
() ()
ijij
t
ijk
k
v
ijk
k
ij
k
ji
k
ji D
x
D
x
P
x
uu
U
t
uu
ε
−Π+
∂
∂
+
∂
∂
+=
∂
∂
+
∂
∂''''
(7.2)
The production term is defined using the exact form
![](asset-73.png)
7-2
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
−=
k
j
ki
k
i
kjij x
U
uu
x
U
uuP '''' (7.3)
Molecular diffusion also follows the exact definition
()
'' ji
k
v
ijk uu
x
D∂
∂
=
ν
(7.4)
There are several models suggested for the turbulent diffusion. Daly and
Harlow1 suggested
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
=
l
ji
lks
t
ijk x
uu
uu
k
CD
''
''
ε
(7.5)
Hanjalic and Launder2 used the form
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
+
∂
∂
=
l
ji
lk
l
ik
lj
l
kj
lis
t
ijk x
uu
uu
x
uu
uu
x
uu
uu
k
CD
''
''
''
''
''
''
ε
(7.6)
Mellor and Herring3 suggests
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
+
∂
∂
=
k
ji
j
ik
i
kj
s
t
ijk x
uu
x
uu
x
uu
k
CD
''
''
''
2
ε
(7.7)
Many models for the pressure-strain term have been developed since the
pioneering work of Launder, Reece, and Rodi4. A general form is given by
()
ikjkjkik
ijmnmnikjkjkik
ijijmnmnkjikijij
WbWbkC
SbSbSbkC
kSCbbbbCbC
++
⎟
⎠
⎞
⎜
⎝
⎛−++
+
⎟
⎠
⎞
⎜
⎝
⎛−
′
+−=Π
4
3
211
3
2
3
1
δ
δε
(7.8)
where the stress anisotropy tensor is given by
ij
ji
ij k
uu
b
δ
3
1
2−
′′
= (7.9)
![](asset-74.png)
7-3
The strain rate tensor is defined as
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
u
x
u
S2
1 (7.10)
and the rotation tensor is given by
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
=
i
j
j
i
ij x
u
x
u
W2
1 (7.11)
If the turbulence is assumed to be locally isotropic (i.e. the same amount of
turbulent energy is dissipated by each of the turbulent energy components 2'
i
u)
then the dissipation term in Eq. 7.1 can be reduced to
ijij
εδε
3
2
= (7.12)
where
ε
is the total rate of energy dissipation.
A transport equation for the dissipation must also be solved to close the system
of equations. The dissipation equation can be written as
()
tj
ikkk
kD
xxxk
CP
k
C
x
U
t
εεεε
ε
ν
εεεε
∂
∂
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+−=
∂
∂
+
∂
∂2
21 (7.13)
The first term on the right hand side of Eq. 7.13 is the production of dissipation.
The second term represents the dissipation of the dissipation. The last two terms
are the molecular and turbulent diffusion of dissipation respectively. The
production of dissipation is given by
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−=
j
i
ji x
U
uuP ''
ε
(7.14)
The dissipation equation turbulent diffusion term is modeled by
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
=
j
lk
tjx
uu
k
CD
ε
ε
εε
'' (7.15)
Typical values for the constants in the above equations for the Launder-Reece-
Rodi4 (LRR), Gibson-Launder5 (GL) and Speziale-Sarker-Gatski6 (SSG) RSM
models are
![](asset-75.png)
7-4
Model 1
C 1
C′ 2
C 3
C 4
C s
C 1
ε
C 2
ε
C
ε
C
LRR 1.5 0.0 0.4 0.6 0.0 0.25 1.44 1.92 0.15
GL 3.6 0.0 0.8 1.2 1.2 0.11 1.44 1.92 0.11
SSG 3.4+1.8P/
ε
4.2 0.8-1.3(bijbij)
1/2
1.25 0.4 0.11 1.44 1.83 0.11
Table 7.1 Coefficients for three Reynolds stress models.
A system of seven transport equations (six Reynolds stress components and one
turbulent dissipation equation) must be solved at each time step to provide the
Reynolds stresses for the Navier-Stokes equations. The above correlations
must be further modified when applied in the presence of a wall, and hence many
of these models utilize wall functions (see Chapter 10) to avoid the complications
and the additional grid points introduced by the presence of the wall.
The Reynolds stresses produced by these models must meet certain constraints
for the models to be applicable to a large range of problems without significant
tuning. Two of the more important constraints are tensor invariance and
realizability. Tensor invariance requires the replaced terms to have the same
tensor form as the original terms. This will insure that the modeled terms
transform properly in different coordinate systems. Schumann7 introduced the
realizability constraint. This constraint requires the equation for the turbulent
stresses to have the property that all the component energies of the turbulent
kinetic energy (the diagonal terms of the Reynolds stress tensor) remain non-
negative and all off-diagonal components of the Reynolds stress tensor satisfy
Schwartz’s inequality. This can be written as
0
'' ≥
iiuu (7.16)
0
2
'''''' ≥− jijjii uuuuuu (7.17)
()
()
0
'
3
'
1
'
2
'
2
'
3
'
2
'
2
'
1
'
3
'
1
'
3
'
1
'
3
'
2
'
3
'
3
'
2
'
1
'
2
'
1
2
'
3
'
2
'
3
'
3
'
2
'
2
'
1
'
1
≥−
+−−
⎟
⎠
⎞
⎜
⎝
⎛−
uuuuuuuuuu
uuuuuuuuuuuuuuuuuu (7.18)
The modeling of the right hand side source terms in Eq. 7.2 can often cause
numerical stability problems with the second moment closure models. The
difficulty in deriving a turbulent dissipation equation was discussed in Chap. 6
and is a source of weakness in this approach of modeling turbulence as well a
weakness in the application of two-equation turbulence models. Because of the
numerical stability problems, the extra computational time requirements, and the
limited improvement in solution accuracy over lower order models (Ref. 8),
second moment closure is not commonly used for many complex CFD problems
![](asset-76.png)
7-5
today. Some cases where second moment closures have shown improvement
over lower order models can be found in Ref. 9.
7.2 Algebraic Stress Models
Under certain assumptions the Reynolds stress transport equations (Eq. 7.1, 7.2)
can be reduced to a system of algebraic equations that require knowledge of the
turbulent kinetic energy (k) and the turbulent dissipation (
ε
). This class of
turbulence models is called an Algebraic Stress Model (ASM). ASM models
reduce the closure problem to solving two transport equations and a system of
algebraic equations. This is significantly faster than solving the full Reynolds
stress transport equations set. The algebraic equations include most of the
models and assumptions that were used to solve the full equation set. Rodi10
suggested that the transport of '' jiuu is proportional to the transport of the
turbulent kinetic energy k, the proportionality factor being the ratio '' jiuu /k (which
is not a constant). This yields
()
ε
−=
⎟
⎠
⎞
⎜
⎝
⎛−=− oductionPr
k
uu
)k(Diff
Dt
Dk
k
uu
)uu(Diff
Dt
uuD 'j
'
i
'j
'
i
'j
'
i
'j
'
i (7.19)
Equation 7.19 is valid when the temporal and spatial change in '' jiuu /k is small
compared with the change of '' jiuu itself. The algebraic stress model of Rodi is
given by
()
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
∂
∂
−
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
−
∂
∂
−
−
+=
1
3
2
1
3
2
''
1
''''''
3
''
ε
δ
ε
δ
j
i
ji
j
i
jiij
l
i
lj
l
i
li
ijji
x
U
uu
C
x
U
uu
x
U
uu
x
U
uu
C
kuu (7.20)
A more recent algebraic stress model derived by Abid, Morrison, Gatski, and
Speziale11 (AMGS) is given by
()
ij
ijklklkjikkijkkjikijkkijtji
k
SSSS
k
WSWS
k
SSuu
δρ
δ
ε
α
ε
αδμρ
3
2
3
1
3
1
254
*
−
⎥
⎦
⎤
⎢
⎣
⎡⎟
⎠
⎞
⎜
⎝
⎛−−++
⎟
⎠
⎞
⎜
⎝
⎛−=
′′
−
(7.21)
where
![](asset-77.png)
7-6
()( )
ε
ρ
ξηξξηη
ξηη
αμ
2
662222
662
1
*
663
2.013 k
t+++++
+++
= (7.22)
()
2/1
ijij SS
k
ε
η
= (7.23)
()
2/1
ijijWW
k
ε
ξ
= (7.24)
() ()
() () ()
12/
1
2
2
2
4
2
4
2
23
4
51
3544
2
2
43
2
2
3221
−+
=−=−=
−=−=
⎟
⎠
⎞
⎜
⎝
⎛−=
CC
ggC
g
C
g
C
g
C
g
C
αα
ααα
(7.25)
The remaining coefficients for the model of Rodi and AMGS are given in Table
7.2.
Model C1 C2 C3 C4 C5
Rodi 1.5 0.0 0.6 0.0 0.0
AMGS 6.8 0.36 1.25 0.40 1.88
Table 7.2 Coefficients for the ASM models of Rodi and AMGS
Rodi’s model in Eq. 7.20 does not include the necessary terms to be valid near
walls and requires wall functions (described in Chapter 10). The AMGS model is
valid down to the wall.
The solution obtained from an ASM model is not independent of the choice of
model used to provide k and
ε
. The pressure coefficient on the Ames
axisymmetric bump for a free stream Mach number of 0.875 is shown for the
AMGS ASM model using a k-
ω
and k-
ε
two equation models for closure in Fig.
7.1. The velocity distribution at the trailing edge of the bump is shown in Fig. 7.2.
Note the solutions are quite different for the same ASM model.
![](asset-78.png)
7-7
Figure 7.1 Pressure coefficient distribution for the axisymmetric bump at
M=0.875.
Figure 7.2 Velocity Distribution at the trailing edge of the bump for the
axisymmetric bump at M=0.875.
7.3 Grid Sensitivity for a Flat Plate with Adiabatic Walls
The initial wall spacing of the computational grid and the grid-stretching ratio can
affect the accuracy any turbulence model. Figure 7.2 shows the sensitivity of the
![](asset-79.png)
7-8
skin friction to initial wall spacing for a flat plate for the SSG RSM model. Flat
plate results for the AMGS ASM model are shown in Fig. 7.3. The grid-stretching
ratio was 1.2 for all these cases. The plots include the theoretical skin friction
curve of White. The boundary layer is seen to become fully turbulent around a
length Reynolds number (Rex) of 1x106 for both models. These models are
relatively insensitive to wall spacing below a y+=10.
Figure 7.3 Flat plate skin friction predictions for the SSG RSM turbulence model
for varying initial wall grid point spacings.
![](asset-7a.png)
7-9
Figure 7.4 Flat plate skin friction predictions for the AMGS ASM turbulence
model for varying initial wall grid point spacings.
Predicted velocity profiles for the flat plate boundary layer for various initial wall
grid point spacings are shown in Fig. 7.5 and 7.6. The velocity profile shows little
effect of the initial spacing for all but the y+=10 profile for the ASM model.
Figure 7.5 Flat plate boundary layer profiles predicted by the SSG RSM
turbulence model for varying initial wall grid point spacings.
![](asset-7b.png)
7-10
Figure 7.6 Flat plate boundary layer profiles predicted by the AMGS ASM
turbulence model for varying initial wall grid point spacings.
The effect of grid stretching ratio on skin friction for a flat plate is shown in Fig.
7.7 and 7.8 All of these solutions used an initial wall spacing of y+=1.
Figure 7.7 The effect of grid stretching ratio on the skin friction for a flat plate
boundary layer using the SSG RSM turbulence model.
![](asset-7c.png)
7-11
Figure 7.8 The effect of grid stretching ratio on the skin friction for a flat plate
boundary layer using the AMGS ASM turbulence model.
The RSM model shows little sensitivity to grid stretching below a ration of 2. The
model has difficulty converging as can be seen by the spikes in the solution. The
ASM solution shows sensitivity beginning at a ration of 1.5. The boundary layer
profile predicted by both of these models show the same sensitivities to skin
friction as the skin friction. This can be seen in Fig.7.9 and Fig. 7.10.
Figure 7.9 The effect of grid stretching ratio on the velocity profile for a flat plate
boundary layer using the SSG RSM turbulence model.
![](asset-7d.png)
7-12
Figure 7.10 The effect of grid stretching ratio on the velocity profile for a flat plate
boundary layer using the AMGS ASM turbulence model.
7.4 Grid Sensitivity for an Axisymmetric Bump
A second example of the grid sensitivity of these turbulence models that includes
a pressure gradient is the NASA Ames transonic axisymmetric bump experiment
described in the Baldwin-Lomax Chapter 4. The effect of initial grid spacing on
the pressure coefficient distribution along the bump for the RSM and ASM
models is shown in Fig. 7.11 and Fig. 7.12. The stretching ratio was 1.2 for
these cases. The pressure coefficient seems to be relatively insensitive to the
initial grid spacing for the ASM model, with the y+=10 and y+=20 curves being
slightly displaced from the other curves. The velocity distribution at the aft
junction of the bump and the cylinder (x/c=1) are shown in Fig. 7.13 and 7.14.
The trends are similar to those seen with the pressure coefficient .
Grid stretching effects on the pressure coefficient distribution along the bump is
shown in Fig. 7.15 and Fig. 7.16. The initial grid spacing was 1.2 for these
cases. Both models show sensitivity to grid stretching ration in the separated
flow region. The velocity profiles shown in Fig. 7.17 and 7.18 show a sensitivity
to stretching ratio. The sensitivity of the Reynolds shear stress predicted by the
SSG RSM model to both initial wall spacing and grid stretching ratio is shown in
Fig. 7.19 and Fig. 7.20. The shear stress prediction is shown to be extremely
sensitive to stretching ratio.
![](asset-7e.png)
7-13
Figure 7.11 The effect of initial grid wall spacing on the pressure coefficient for
the axisymmetric bump using the SSG RSM turbulence model.
Figure 7.12 The effect of initial grid wall spacing on the pressure coefficient for
the axisymmetric bump using AMGS ASM turbulence model.
![](asset-7f.png)
7-14
Figure 7.13 The effect of initial grid wall spacing on the velocity distribution at
x/c=1 for the axisymmetric bump using the SSG RSM turbulence model.
Figure 7.14 The effect of initial grid wall spacing on the velocity distribution at
x/c=1 for the axisymmetric bump using the AMGS ASM turbulence model.
![](asset-80.png)
7-15
Figure 7.15 The effect of grid stretching on the pressure coefficient for the
axisymmetric bump using the SSG RMS turbulence model.
Figure 7.16 The effect of grid stretching on the pressure coefficient for the
axisymmetric bump using the AMGS ASM turbulence model.
![](asset-81.png)
7-16
Figure 7.17 The effect of grid stretching ratio on the velocity distribution at x/c=1
for the axisymmetric bump using the SSG RSM turbulence model.
Figure 7.18 The effect of grid stretching ratio on the velocity distribution at x/c=1
for the axisymmetric bump using the AMGS ASM turbulence model.
![](asset-82.png)
7-17
Figure 7.19 The effect of initial wall spacing on the shear stress distribution at
x/c=1 for the axisymmetric bump using the SSG RSM turbulence model.
Figure 7.20 The effect of grid stretching ratio on the shear stress distribution at
x/c=1 for the axisymmetric bump using the SSG RSM turbulence model.
![](asset-83.png)
7-18
Chapter 7 References:
1. Daly, B. and Harlow, F., “Transport Equations of Turbulence,” Physics of
Fluids, 1970, pp. 2634-2649.
2. Hanjalic, K. and Launder, B., “Contribution Towards a Reynolds-Stress
Closure for Low-Reynolds-Number Turbulence,” Journal of Fluid Mechanics,
Vol. 74, 1976, pp. 593-610.
3. Mellor, G. and Herring, H., “A Survey of the Mean Turbulent Field Closure
Models,” AIAA Journal, Vol. 11, 1973, pp. 590-599.
4. Launder, B., Reece, G., and Rodi, W., “Progress in the development of a
Reynolds-Stress Turbulence Closure,” Journal of Fluid Mechanics, 1975, Vol.
68, Part 3, pp. 537-566.
5. Gibson, M. and Launder, B., “On the Calculation of Horizontal, Turbulent,
Free Shear Flows under Gravitational Influence,” Journal of Heat Transfer,
Vol. 98, 1976, pp. 81-87.
6. Speziale, C., Sarker, S., and Gibson, T., “Modeling the Pressure-Strain
Correlation of Turbulence: An Invariant Dynamical Systems Approach,”
Journal of Fluid Mechanics, Vol. 227, 1991, pp. 245-272.
7. Schumann, U., “Realizability of Reynolds Stress Turbulence Models,” Physics
of Fluids, Vol. 20, pp. 721-725, 1977.
8. Leschziner, M., Batten, P. and Craft, T., “Reynolds-Stress Modeling of
Transonic Afterbody Flows,” The Aeronautical Journal, June 2001, pp. 297-
306.
9. Nallasamy, M., “Turbulence Models and Their Applications to the Prediction
of Internal Flows: A Review,” Computers and Fluids, Vol. 15, pp. 151-194,
1987.
10. Rodi, W., “Turbulence Models and Their Application in Hydraulics – A State of
the Art Review,” International Association for Hydraulic Research, 2nd Ed.,
Feb. 1984.
11. Abid, R., Morrison, J., Gatski, T., and Speziale, C., “Prediction of
Aerodynamic Flows with a New Explicit Algebraic Stress Model,” AIAA
Journal, Vol. 34, No. 12, Dec. 1996.
RSM and ASM Model Application Hints
1. RSM models generally suffer from numerical stability issues. This is
due to the complexity of the modeled terms. This is the primary
reason that these models are not regularly used in large scale
applications.
2. ASM models are sensitive to the transport equations used to obtain
k and
ε
. These models are generally less numerical stable than
traditional two-equation models.
3. The choice of the pressure strain term is critical to the performance
of both RSM and ASM models. More complex pressure strain
models are generally more accurate but less numerically stable.
![](asset-84.png)
8-1
8.0 Large Eddy Simulation
Theoretically the Navier-Stokes equations can be used to simulate turbulent
flows. The computational grid used in such a simulation would have to be fine
enough to allow the smallest turbulent length scales to be realized and the
computational time step would have to be small enough to simulate the highest
frequencies of the turbulent spectrum. In practice such computations are
prohibitively expensive for high Reynolds numbers. The Reynolds averaged
Navier-Stokes (RANS) equations were derived assuming that all of the
unsteadiness due to the turbulent nature of the flow could be modeled with
empirically derived correlations. This reduces the time and length scales that
must be simulated, but also limits the applicability of the simulation for unsteady
flows. Large eddy simulations (LES) were developed to extend the simulation of
unsteady flows beyond DNS. The desired result of an LES computation is to
obtain a DNS equivalent solution for the large-scale turbulence on a much
coarser grid than is required for DNS. An LES simulation requires:
1. A grid fine enough to discretize the small nearly isotropic scales of the
turbulence
2. A low dissipation numerical scheme
3. A filter function to determine the division of the turbulent spectrum into
grid realized and subgrid regions
4. A subgrid turbulence model
A true LES simulation is more than a high Reynolds number computation run
without a turbulence model. Although the resulting solution from such a
simulation may resemble turbulent flow, the resulting solution will most likely not
represent an equivalent DNS solution.
8.1 The Filtering Operation
The filtering operation is critical to LES. Consider a filtering operation with a
uniform characteristic filter width Δ (which implies isotropic grid elements).
Leonard1 defined the following filter in physical space
∫∞
∞− −=
ξξφξφ
dxGx )()()( (8.1)
Note that in this chapter the over bar represents a filtered quantity, not a time
averaged quantity as in the previous chapters. The filtering operation is a spatial
operation as opposed to the Reynolds averaging operation discussed in Chapter
2 that is a temporal operation. The original function
φ
is then decomposed into a
filtered field (or grid resolved) and a “subgrid” term
φφφ
′
+= (8.2)
![](asset-85.png)
8-2
The function G defined in Eq. 8.1 is the “filter function”. The filter function may be
any function defined on an infinite domain that satisfies the following
requirements (Ref. 1):
1. )()(
ξ
ξ
GG =−
2. ∫
∞
∞−
=1)(
ξξ
dG
3. 0)( →
ξ
G as 0→
ξ
such that all moments ∫
∞
∞−
ξξξ
dG n
)(
(
)
0≥n exist
4. )(
ξ
G is “small” outside ⎟
⎠
⎞
⎜
⎝
⎛ΔΔ
−2
,
2
One important and useful feature of this choice of filter is
xx ∂
∂
=
∂
∂
φφ
(8.3)
Three common filters are
Spectral cutoff filter:
Δ
⎟
⎠
⎞
⎜
⎝
⎛
Δ
=
πξ
πξ
ξ
sin
)(
G (8.4)
Gaussian filter: ⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
ΔΔ
=2
2
exp
1
)(
ξ
π
ξ
CC
G (8.5)
Box filter:
⎪
⎪
⎩
⎪
⎪
⎨
⎧Δ
<
Δ
=
otherwise
if
G
0
2
1
)(
ξ
ξ
(8.6)
The spectral cutoff filter is normally applied in spectral space as
)(
ˆ
)(
ˆ
)(
ˆkkGk
ϕϕ
= (8.8)
![](asset-86.png)
8-3
⎪
⎪
⎩
⎪
⎪
⎨
⎧
Δ
<
=
otherwise
kif
G
0
2
1
)(
ˆ
π
ξ
(8.9)
The constant C in the Gaussian filter (Eq. 8.5) is somewhat arbitrary, and values
from 2 to 6 have been used in practice. Note that the box filter (Eq. 8.6) is valid
on both finite and infinite domains. For this reason box filters are often used to
relate DNS to LES for physical space numerical schemes.
There are advantages to using spectral and pseudo-spectral solution
methodologies for LES. The first advantage is that the approximate field, which
is discrete in spectral space, is a finite sum of continuous functions in physical
space. Thus the approximate field, the spectral filter, and their derivatives are
continuous functions in physical space. The spectral filters defined in Eqs. 8.8
and 8.9 have the following additional useful properties
φφ
= (8.10)
0=
′
φ
(8.11)
This says that filtering a grid resolved quantity yields the original grid resolved
quantity and that the spectral filter of a subgrid quantity is equal to zero.
Another advantage of the spectral approach is that derivatives yield an exact
value for the approximate field. In other words, if one is approximating a function
using N Fourier modes, then as long as one has at least 2(N+1) points in
physical space, it does not matter if one has 20 or 20 million points: the value of
the derivative at any given point is the same. Therefore one can take a flow field
generated with a spectral DNS code, filter and coarsen the solution, and produce
a flow field which satisfies the governing equations for LES.
There are drawbacks to spectral filters. First, they have limited applicability in
realistic flow situations because of boundary constraints. Also, spectral filters are
non-positive (i.e. the filter function is negative at some points in space). The
latter condition results in a subgrid stress tensor that does not satisfy the
Reynolds stress realizability conditions outlined in Chapter 2.
Both Gaussian filters (Eq. 8.5) and box filters (Eq. 8.6) are positive functions, and
thus the Reynolds stress realizability conditions outlined in Chapter 2 will be
satisfied on the resulting subgrid stress tensor. Unfortunately Eq. 8.10 and Eq.
8.11 are not valid for these filters. Furthermore, if a spatial numerical solver is to
be used the filter must also be discretized. Therefore it is necessary to examine
the discretized system that is solved in the actual code because discrete systems
![](asset-87.png)
8-4
sometimes have different properties than the original continuous systems from
which they were derived.
8.2 Derivation of the LES Equations
The Favre averaged or filtered Navier-Stokes equations were derived in Chapter
2. The Favre filter can be defined for any variable as
ρ
ρφ
φ
=
~ (8.12)
Note that the over bar signifies a spatially filtered quantity, and the tilde
represents a Favre filtered, or grid resolved, quantity. Thus the Favre filter may
be thought of as a density weighted filter in space. Applying this filtering
operation to the Navier-Stokes equations, and assuming that the filtering
commutes with the derivative operation, the LES equations are
0
~
=
∂
∂
+
∂
∂
i
i
x
u
t
ρ
ρ
(8.13)
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
−−
∂
∂
+
∂
∂
+
∂
∂
−=
∂
∂
+
∂
∂
jiji
j
ijij
jj
ij
ij
ji
iuuuu
xxxx
p
x
uu
t
u~~
~
~
~~~
~
ρττ
τρ
ρ
(8.14)
()
()
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
∂
∂
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+−
∂
∂
+
∂
∂
+
−
∂
∂
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
−=+
∂
∂
+
∂
∂
iiiii
ijjijj
ii
ijj
ii
i
ii
i
i
i
x
T
x
T
xx
T
x
uu
xx
u
uppu
x
uEEu
x
upE
xt
E
~
~
~
~
~~
~~
~~
~
~
~
~
κκκττ
τ
ρρ
ρ
(8.15)
The resolved viscous stress tensor in Eqn. 8.14 takes the form
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
∂
∂
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=ij
k
k
i
j
j
i
ij x
u
x
u
x
u
δμτ
~
3
1
~
~
2
1
~
2
~ (8.16)
The viscosity and thermal conductivity is assumed to be calculated from the
Favre averaged temperature (
T
~
). The total energy (
E
~
) is defined as
sgs
ji kuueE ++= ~~
2
1
~
~ (8.17)
The subgrid kinetic energy, ksgs, is defined as the effect of the subgrid scales on
the kinetic energy of the resolved field
![](asset-88.png)
8-5
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−= iiii
sgs uuuuk ~~
~
2
1 (8.18)
Finally, the LES equation of state is
TRp
~
ρ
= (8.19)
As with Reynolds averaging, the filtering operation produces terms that must be
modeled in order to close the equation set. The LES momentum equation, Eqn.
8.14, contains two such terms. The first,
()
ijij
j
x
ττ
~
−
∂
∂
, represents the difference
in the viscous terms between Favre and “straight” filtering. The second term,
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
jiji
j
uuuu
x~~
~
ρ
, arises from the convection terms and has been observed to
behave much like a viscous stress. Thus it is frequently redefined in terms of the
“subgrid stress tensor” sgs
ij
τ
defined as
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−= jiji
sgs
ij uuuu ~~
~
ρτ
(8.20)
The LES energy equation, Eq. 8.15, contains four subgrid terms. The first term,
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
ii
i
uEEu
x~
~
~
ρ
, arises from the convection term. The second term is a
pressure-velocity correlation term,
(
)
ii
i
uppu
x~
−
∂
∂
, which is also related to the
convection term. The third term,
(
)
ijjijj
i
uu
x
ττ
~~
−
∂
∂
, is a viscous subgrid term
similar to the term in the LES momentum equation. This term represents the
transfer of energy due to subgrid viscous forces. The fourth term,
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
∂
∂
iii x
T
x
T
x
~
~
κκ
, is a heat flux subgrid term. The subgrid kinetic energy, ksgs,
used in the definition of the Favre filtered total energy in Eq. 8.7 must also be
defined for closure.
The subgrid stress tensor, sgs
ij
τ
is usually treated in a similar manner as the
Boussinesq approximation for the RANS equations (Chapter 2)
![](asset-89.png)
8-6
ij
sgs
ijkkijt
sgs
ij kSS
δρδνρτ
3
2
~
3
1
~
2+
⎟
⎠
⎞
⎜
⎝
⎛−−= (8.21)
In the above equation, ij
S
~
represents the resolved rate of strain tensor
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
+
∂
∂
=
i
j
j
i
ij x
u
x
u
S
~
~
2
1
~ (8.22)
Note that a subgrid eddy viscosity,
ν
t, has been introduced in Eq. 8.21. The
subgrid eddy viscosity accounts for the turbulence that cannot be resolved on the
computational grid.
Introducing the subgrid stress tensor into Eq. 8.14 and assuming that the straight
filtered stress tensor is equal to the Favre filtered stress tensor
(
)
ijij
τ
τ
~
= the LES
momentum equation becomes
j
sgs
ij
j
ij
ij
ji
i
xxx
p
x
uu
t
u
∂
∂
−
∂
∂
+
∂
∂
−=
∂
∂
+
∂
∂
ττρ
ρ
~
~
~
~
(8.23)
The LES energy equation (Eq. 8.15) can be further simplified as follows. Using
the definition of total enthalpy
ρ
p
EH += (8.24)
two of the subgrid terms can be combined to form
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
=−
∂
∂
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
−ii
i
ii
i
ii
i
uHHu
x
uppu
x
uEEu
x~
~
~
~
~~
~
ρρρρ
(8.25)
This term can then be modeled with an eddy diffusion model as
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
Δ−
∂
∂
≈
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
∂
∂
i
sgs
e
i
ii
ix
H
kc
x
uHHu
x
~
~
~
~
ρρρ
(8.26)
where Δ is the local grid spacing and ce is an empirical coefficient. The thermal
conductivity term, ⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
∂
∂
iii x
T
x
T
x
~
~
κκ
, is neglected based on the assumption that
the thermal conductivity is locally constant and that the Favre filtering is
equivalent to the straight filtering. The viscous transport term is modeled as
![](asset-8a.png)
8-7
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
≈−
∂
∂
i
sgs
i
ijjijj
ix
k
x
uu
x
μττ
~~ (8.27)
This term is generally small and is often neglected in practice.
Thus the modeled LES energy equation may now be written as
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+
∂
∂
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
Δ
∂
∂
=+
∂
∂
+
∂
∂
iii
sgs
i
i
ijj
i
sgs
e
i
i
i
x
T
xx
k
x
x
u
x
H
kc
x
upE
xt
E
~
~
~
~
~
~
~
~
κμ
τ
ρρ
ρ
(8.28)
Inspection of the modeled LES equations (Eqs. 8.13, 8.23, and 8.28) shows that
these equations are quite similar to the RANS equations incorporating the
Boussinesq approximation derived in Chapter 2. It should be noted that the time
averaging process used to derive the RANS equations and the spatial filtering
operations used to derive the LES equations are quite different, and hence the
terms in the two equation sets are not the same.
Three things are required to close the modeled LES equations:
1. The subgrid turbulent kinetic energy ksgs
2. The subgrid eddy viscosity νt
3. The empirical coefficient in the energy equation, ce
8.3 Smagorinsky Model
As with the RANS equations, there are numerous models for the LES eddy
viscosity. One of the earliest models for the LES eddy viscosity was proposed by
Smagorinski2 and is given by
()
()
2
1
2~~
2ijijgst SSC Δ=
ν
(8.29)
where Cs is the Smagorinsky coefficient and Δg is the local grid spacing. This
model ignores the subgrid kinetic energy ksgs. The Smagorinsky model provides
a simple closure for the LES equations and has been used effectively for a
number of applications. But, like its algebraic counterpart in the RANS regime, it
has been shown to be lacking in simulating complex turbulent flows. The
Smagorinsky model is not valid near walls, and a wall damping term (Ref. 3) is
often added to the model. The wall-damped form is given by
![](asset-8b.png)
8-8
()()
()
2
1
2
25/ ~~
21 ijij
y
gst SSeC +
−
−Δ=
ν
(8.30)
The coefficient Cs must be defined for this model. In practice no universal value
of Cs exist.
8.4 Dynamic Smagorinsky Model
Several investigators have attempted to dynamically calculate Cs based on
equating the highest wave number resolved turbulent stresses with the subgrid
stress. One dynamic LES closure model outlined in Ref. 4 can be described as
follows. The subgrid stress tensor is written as
()
)
~
()
~
(, 2uSuStxC ij
sgs
ij Δ−=
τ
(8.31)
where C(x,t) is the Smagorinsky coefficient to be determined dynamically and Δ is
the grid spacing. For this purpose a second spatial filter, called the test-filter, of
width larger than the grid filter is introduced. We choose the test-filter scale
Δ=Δ 2. This filter generates a second set of resolvable-scale fields (denoted by
()
). Then, the dynamic SGS model of the ij
τ
at the grid filter and of the ij
T at the
test-filter are written as:
⎪
⎪
⎪
⎩
⎪
⎪
⎪
⎨
⎧
−=
Δ−=−
−=
Δ−=−
ij
ij
kk
ij
ij
ij
ij
kk
ij
ij
C
SSCTT
C
SSC
α
δ
β
τ
δ
τ
2
2
3
2
2
3
2
2
(8.32)
A least squares method is developed to predict C(x,t) as:
(
)
(
)
(
)
mnmnmnmn
ijijij
txC
βαβα
βα
−−
−
−= l
2
1
),( (8.33)
with (/3) 2 2
a
ij ij ij kk ij ij
CC=− =− +ll l
δ
αβ
.
8.5 k-Equation Model
Another approach to closing the LES equations is to introduce a transport
equation for the subgrid kinetic energy, ksgs. The exact subgrid kinetic energy
equation may be written
![](asset-8c.png)
8-9
()
()
i
j
sgs
ij
i
sgs
ijj
ii
i
i
ij
i
i
ij
i
i
i
i
ii
sgs
i
sgs
x
u
x
u
ukku
x
x
u
x
u
x
p
u
x
p
uuk
xt
k
∂
∂
+
∂
∂
−−
∂
∂
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
−
∂
∂
−=
∂
∂
+
∂
∂
~~
~
~
~~~
τ
τ
ρρ
ττ
ρ
ρ
(8.34)
This equation contains several subgrid terms. The equation is modeled in a
similar manner as the RANS turbulent kinetic energy equation. The modeled
equation is
()
(
)
Δ
−
∂
∂
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
=
∂
∂
+
∂
∂2
3
~
~sgs
i
j
sgs
ij
i
sgs
t
i
i
sgs
i
sgs k
c
x
u
x
k
x
uk
xt
k
ε
ρτνρρ
ρ
(8.35)
where the terms on the right-hand-side are the diffusion, production, and
dissipation of the subgrid kinetic energy respectively. An algebraic expression
for the dissipation is generally used instead of a transport equation as in two-
equation RANS turbulence models. The turbulent viscosity is given by
Δ= sgs
tkc
ν
ν
(8.36)
The LES energy equation dissipation coefficient c
ε
and the eddy viscosity
coefficient c
ν
are determined dynamically by equating the highest wave number
resolved turbulent properties with the subgrid properties. The dynamic process
for determining these coefficients is described in detail in Ref. 5.
8.6 Inflow Turbulence Boundary Condition
Flows that include a developing boundary layer at the inflow plane of a simulation
require a boundary condition that includes the specification of the inflow velocity
fluctuations. Soteriou and Ghoniem6 show the development of an
incompressible mixing layer with and without perturbation at the inflow plane.
Compared to the perturbed case, the development of the mixing layer without
perturbation is delayed significantly. This is because disturbances within the
boundary layer are amplified and lead to the Kelvin-Helmholtz instability that
eventually causes the shear layer to roll up. If these disturbances are not
present in the boundary layer, then the shear layer will develop more slowly until
disturbances are generated through minute numerical errors.
Several approaches have been taken to add the perturbations to the inflow. The
simplest approach is to add “white noise” to the velocity at the inflow plane as
was done by Comte et al7. This approach is somewhat unphysical in that the
perturbations have no correlation in space or time. A second approach is to add
perturbations at discreet frequencies as was done by Soteriou and Ghoniem6.
This can be effective for studying forced flows, but it must be remembered that
turbulence is broadband and does not confine itself to discrete frequencies. A
![](asset-8d.png)
8-10
third approach to this problem is to attempt to reconstruct the unsteady flow at
the inflow using correlations for the perturbations. One example of this can be
found in Ref. 5, in which a box of frozen turbulence is generated and saved. This
box is scaled and applied as perturbations at the inflow of the simulation.
8.6 Other LES References
LES methods are still evolving as researches investigate methods for simulating
the near-wall flow physics and methods to improve the performance of the
subgrid models used in LES simulations. Fureby, et al8, provide a good
summary of the LES methods currently being employed and investigated. Ref. 9
includes a wide range of experimental test cases for validation of LES
simulations.
8.7 Spatial Mixing-Layer Example
The mixing layer experiment of Samimy and Elliot8 can be used to demonstrate
the behavior of LES models. The data is also included in Ref. 10. The
experimental setup had an upper supersonic stream and a lower subsonic
stream mixing at a matched static pressure. Laser Doppler Velocimeter (LDV)
measurements were made at several downstream locations between the trailing
edge of the splitter plate and a station 210 mm downstream of the splitter plate.
The flow parameters are given in Table 8.1.
T0, K P01,
kPa M1 M
2 M
c U
1,
m/sec U2/U1 ρ2/ρ1 δ1, mm
291.0 314.0 1.80 0.51 0.52 479.5 0.355 0.638 8.0
Table 8.1 Flow parameters for the spatial mixing-layer case
The solutions were run for 30,000 time steps to wash out initial transients and
then statistics were taken over the next 60,000 time steps. The time step used
was 5.0x10-8 seconds. The flow solver used for these calculations was 4th order
in space and 2nd order in time. The computational grid dimensions were 181 x
121 x 61. The inflow plane perturbations were specified using the “box of
turbulence” approach from Ref. 5. Further details of the experiment and the
computations can be found in Ref. 10.
All of the turbulence models investigated predict roughly the same mixing-layer
thickness for this case as shown in Fig 8.1. The Hybrid k-ε model in Fig. 8.1 is
the multi-scale hybrid RANS/LES turbulence model described in Chapter 9. The
DES model is the Spalart-Allmaras DES hybrid model that is also described in
Chapter 9. The traditional RANS models (k-ε and Spalart-Allmaras) have the
worst agreement with the experimental data.
![](asset-8e.png)
8-11
Figure 8.1 Mixing-layer thickness
The streamwise velocity 210 mm downstream of the splitter plate is shown in Fig.
8.2. Again all of the models are in reasonable agreement with the experimental
data.
Figure 8.2 Nondimensional streamwise velocity at x=210 mm
![](asset-8f.png)
8-12
The streamwise turbulence intensity, Reynolds stress, streamwise velocity
fluctuation skewness, and streamwise velocity fluctuation flatness at x=210 mm
are shown in Fig. 8.3-8.6. The performance of the traditional RANS models is
seen to deteriorate as the order of the turbulent velocity fluctuation correlation
increases. The LES and hybrid RANS/LES models are all in reasonable
agreement with the experimental data.
Figure 8.3 Streamwise turbulence intensity at x=210 mm
Figure 8.4 Reynolds stress profiles at x=210 mm
![](asset-90.png)
8-13
Figure 8.5 Stream wise velocity fluctuation skewness at x=210 mm
Figure 8.6 Streamwise velocity fluctuation flatness at x=210 mm.
Chapter 8 References:
1. Leonard, A., “Energy Cascade in Large-Eddy Simulations of Turbulent
Fluid Flows,” Advances in Geophysics, Vol. 18A, pp. 237-348, 1974.
2. Smagorinsky, J. S., “General Circulation Experiments with the Primitive
Equations I: The Basic Experiment,” Monthly Weather Review, Vol. 91,
pp. 99-164, 1963.
3. Piomelli, U., Zang, T. A., Speziale, C. G., and Hussaini, M. Y., Physics of
Fluids, Vol. 2, pp. 257, 1990.
4. Wasistho, B., Haselbacher, A., Najjar, F.M., Tafti, D., Balachandar, S. and
Moser, R.D., “Direct and Large Eddy Simulations of Compressible Wall-
8-14
Injection Flows in Laminar, Transitional and Turbulent Regimes,” AIAA
2002-4344, 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and
Exhibit, July 7-10, 2002, Indianapolis, IN.
5. Nelson, C. C., “Simulations of Spatially Resolved Compressible
Turbulence Using a Local Dynamic Subgrid Model,” Ph.D. Thesis,
Georgia Institute of Technology, 1997.
6. Soteriou, M. and Ghoniem, A., “Effects of the Free-Stream Density Ratio
on Free and Forced Spatially Developing Shear Layers,” Physics of
Fluids, Vol. 7, No. 8, 1995, pp. 2036-2051.
7. Comte, P., Lesieur, M., Laroche, H., and Normand, X., “Numerical
Simulations of Turbulent Plane Shear Layers,” in Turbulent Shear Flows
6, Springer-Velag, 1989, pp. 360-380.
8. Fureby, C., Alin,N., Wikstrom, N., Menon, S., Svanstedt, N. and Persson,
L., “Large-Eddy Simulation of High-Reynolds-Number Wall-Bounded
Flows,” AIAA Journal, Vol. 42, No. 3, March 2004, pp 457-468.
9. “A Selection of Test Cases for the Validation of Large-Eddy Simulations of
Turbulent Flows,” AGARD-AR-345, April 1998.
10. Samimy, M. and Elliot, G., “Effects of Compressibility on the
Characteristics of Free Shear Layers,” AIAA Journal, Vol. 28, No. 3, 1990,
pp.439-445.
11. Nelson, C. and Nichols, R., “Evaluation of Hybrid RANS/LES Turbulence
Models Using an LES Code,” AIAA-2003-3552, June 2003.
![](asset-92.png)
9-1
9.0 Hybrid RANS/LES Models
9.1 Theory
Several investigators have noted a limitation with RANS turbulence models when
applied to unsteady flows. The turbulence models produce too much eddy
viscosity and over-damp the unsteady motion of the fluid. The problem is
inherent in the construction of the turbulence models because of the assumption
that all scales of the unsteady motion of the fluid are to be captured and modeled
by the turbulence model.
One approach to overcoming the shortcomings of RANS models applied to
unsteady flows is to spatially filter the RANS turbulence models such that the
eddy viscosity does not include the energy of grid resolved turbulent scales. The
spatially filtered RANS turbulence model may be thought of as a subgrid model
for very large turbulent eddies. This class of turbulence models has been called
hybrid RANS/LES models because they incorporate aspects of both forms of
turbulence modeling. It is desirable that the spatial filter function chosen not
degrade the performance of the turbulence model when the largest turbulent
scales present are below the resolution of the grid as is often the case in current
aircraft CFD applications.
There are two possibilities for developing hybrid models from existing RANS
models. A first approach would be to modify the production or dissipation source
terms for existing turbulence model differential equations to include additional
terms to adjust the local eddy viscosity so that it does not include the grid
realized contribution. This means that the turbulence quantities normally
transported by the RANS turbulence model will be created or destroyed based on
the local grid resolution.
A second approach would be to solve the existing RANS turbulence model in the
normal manner and then filter the results to determine the level of eddy viscosity
that will be used in the solution of the Navier-Stokes equations. This approach is
somewhat simpler to develop since it requires no tuning of the differential
transport equations and can be easily extended to different RANS turbulence
models.
Both approaches have been applied to unsteady flows. One example of the first
form of a hybrid RANS/LES turbulence model based on the Spalart-Allmaras
one-equation turbulence model can be found in Ref. 1. The standard Spalart-
Allmaras turbulence model contains a destruction term for eddy viscosity that is
inversely proportional to the square of the distance from the wall (d). Ref. 1
suggests replacing the wall distance (d) in the destruction term with
()
gDES LCdd ,min
~
= (9.1)
![](asset-93.png)
9-2
where CDES is a constant of O(1) and Lg is a grid length scale defined by
()
zyxLgΔΔΔ= ,,max (9.2)
where Δx, Δy, and Δz are the local grid lengths. For two-dimensional or
axisymmetric problems, the out of plane direction is simply ignored in defining the
local grid length. The modified destruction term has the effect of decreasing the
eddy viscosity in regions of tight grid spacing. This modification causes the
Spalart-Allmaras RANS turbulence model to behave like a Smagorinsky LES
turbulence model (see Chapter 8) when the grid spacing (Lg) is less than the
distance from the wall, which is generally the case outside of the boundary layer.
Note that the transition from RANS to LES does not include any turbulent length
scale dependence, but is solely a function of the local grid spacing. Ref. 1
introduces the term Detached-Eddy Simulation (DES) to describe this model.
DES has been applied to a number of unsteady flow problems including flow over
a sphere2, flow over a delta wing3 and flow over an aircraft4.
A similar modification to an SST two-equation model is given by Strelets5. In this
hybrid model the dissipation term (
ε
) in the turbulent kinetic energy equation (k) is
replaced by
()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
===
t
g
DES
gDESt
L
L
C
LCL
k
k,0.1min
,min
2
3
*
ε
β
ω
ε
(9.3)
where
ε
ωβ
2
3
*
kk
Lt== (9.4)
()
DESKWDESKEDES CFCFC 11
1+−= (9.5)
here
ω
is the specific dissipation,
β
*=0.09, CDESKE=0.61, and CDESKW=0.78. The
turbulent dissipation
ε
defined in Eq. 9.3 is effectively increased when the grid
size length scale Lg (Eq. 9.2) is less than the turbulent length scale Lt (Eq. 9.4).
This causes the eddy viscosity to be reduced in these regions. Unlike the one-
equation SA-DES model, this model does include a dependency on the local
turbulent length scale. The turbulent length scale used in this model is a function
of the filtered (subgrid) turbulent quantities. This model behaves like a k-
equation LES subgrid model when the turbulent length scale is greater than the
grid length scale, and the dissipation equation is decoupled from the kinetic
energy equation in this region.
![](asset-94.png)
9-3
Nichols and Nelson6 give an example of the second approach for developing a
hybrid RANS/LES turbulence model which they have designated a multi-scale
model. The method was implemented in conjunction with the SST two-equation
turbulence model and is termed a multi-scale model. The turbulent length scale
used in this effort is defined by
(
)
/,/0.6max 2/3 RANSRANStRANSt kL
εν
Ω= (9.6)
where
ν
tRANS is the unfiltered (large scale) eddy viscosity and Ω is the local mean
flow vorticity. This length scale is a mixture of the traditional turbulent scale
definition for two-equation turbulence models ( 2/3
RANS
k/
ε
RANS) and the definition
usually associated with algebraic turbulence models ((
ν
tRANS/Ω)1/2). The turbulent
length scale definition could be easily adapted to other types of turbulence
models. The subgrid turbulent kinetic energy is defined as
dRANSLES fkk = (9.7)
The damping function is defined as
()()()
2/5.02tanh1 −Λ+=
π
d
f (9.8)
where
3/4
g
t
L
L
0.1
0.1
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+
=Λ (9.9)
where Lg is defined in Eq. 9.2. The eddy viscosity is then calculated from
()
1 tLESddtRANSt ff
ν
ν
ν
−+= (9.10)
where the LES based subgrid eddy viscosity is given by
),0854.0min( tRANS
νν
LESgtLES kL= (9.11)
The multi-scale hybrid model behaves like a traditional SST model on the RANS
end of the spectrum and transitions to a nonlinear k-equation model on the LES
end of the spectrum. The transition function in Eq. 9.10 was chosen to allow a
smooth transition from the standard RANS turbulence model to the LES subgrid
model. This hybrid RANS/LES approach can easily be extended to other RANS
turbulence models with little to no alteration providing that a length scale can be
derived for that model and a value for the turbulent kinetic energy can be
approximated by that model. The multi-scale hybrid model, like the SST-DES
![](asset-95.png)
9-4
model, transitions from RANS to LES as a function of the ratio of the local
turbulent length scale predicted by the RANS model and the local grid spacing
rather than being a function of the grid spacing alone as is the case for the
Spalart-Allmaras DES model.
Nelson and Nichols evaluated the Spalart-Allmaras DES and the multi-scale
hybrid turbulence models as subgrid turbulence models for LES applications in
Ref. 7. The hybrid models were applied to a high-speed shear layer using a true
LES flow solver. The results indicate that the hybrid models perform as well as
more complicated LES subgrid models for this application. Nichols8 investigated
these three hybrid RANS/LES models applied to a circular cylinder and a generic
weapons bay. These results are shown in the following sections.
9.2 Circular Cylinder
Unsteady three-dimensional calculations were performed for the vortex shedding
from a circular cylinder for M=0.2 and Red=8x106. Three computational grids
were used in the simulation:
Fine – 401x201x201
Mid – 201x101x101
Coarse – 101x51x51
The mid and coarse level grids were constructed by removing every other point
from the fine and mid level grids respectively. The grids have a span of ten
cylinder diameters. Periodic boundary conditions were applied at the side planes
of the grid. The initial wall spacing was 2x10-4 diameters for the fine grid, which
corresponds to a y+ of about 20. Wall function boundary conditions were used in
the calculations.
Calculations with three different time steps were run on the mid grid with the
multi-scale hybrid turbulence model to determine an acceptable time step for this
simulation. The results are shown in Fig. 9.1. The primary shedding frequency
for this example is just less than 60 Hz (Strouhal number=0.28). These solutions
correspond to 50, 200, and 800 time steps per shedding cycle. The results from
the largest time scale are quite different from the smaller time step results. The
spectral peak predicted with the large time step occurs at a much lower
frequency than the peak predicted using the smaller time steps. Based on these
results all subsequent simulations were performed using the 9.0E-4 second time
step with 200 time steps per shedding cycle.
![](asset-96.png)
9-5
Figure 9.1. Normal force power spectral density for 3 different time steps using
the SST-MS model on the medium grid.
The data in Ref. 9 was obtained using both air (γ=1.4) and Freon (γ=1.13) as a
test medium. Calculations were performed on the medium grid with the SST-MS
hybrid model to assess the affect of test medium. The normal force spectral
results are shown in Fig. 9.2. This indicates that the test medium is not an
important factor at the conditions of this study. All subsequent solutions were
performed using air as the test medium.
Figure 9.2. Normal force power spectral density for air and Freon using the SST-
MS hybrid model.
![](asset-97.png)
9-6
The calculations were run 10,000 iterations and the final 4096 time steps were
statistically analyzed. The centerline plane of the three grids is shown in Fig. 9.3.
Fine Grid Medium Grid Coarse Grid
Figure 9.3. Cylinder grid centerline.
A comparison of the normal force power spectral density for the three hybrid
models, the Spalart-Allmaras model, and the SST model for the medium grid are
shown in Fig. 9.4. The hybrid models show similar trends in that they all have a
rather broad spectral peak at a Strouhal number of about 0.28. The two RANS
models have a much narrower spectral peak at a similar Strouhal number. The
RANS models have much lower energy away from the peak. This is an
indication that the RANS models are providing too much damping of the
unsteady solutions.
Figure 9.4. Normal force power spectral density on the medium grid.
![](asset-98.png)
9-7
Instantaneous x-vorticity isosurfaces computed with the SST-DES hybrid
turbulence model are shown in Fig. 9.5 for all three grids. The turbulent structure
in the wake of the cylinder is clearly evident in the fine grid solution. The
turbulent structure is reduced in the mid grid solution. There is very little
structure present in the coarse grid solution.
Fine Grid Medium Grid Coarse Grid
Figure 9.5. Instantaneous vorticity isosurfaces colored by Mach number for the
SST-DES hybrid turbulence model.
Instantaneous contours of Mach number and eddy viscosity for the SA-DES,
SST-DES, and SST-MS hybrid turbulence models computed on the three grid
levels are shown in Figs. 9.6-9.11.
a. Fine b. Mid c. Coarse
Figure 9.6. Instantaneous Mach number contours for the SA-DES hybrid
turbulence model.
![](asset-99.png)
9-8
a. Fine b. Mid c. Coarse
Figure 9.7. Instantaneous eddy viscosity contours for the SA-DES hybrid
turbulence model.
a. Fine b. Mid c. Coarse
Figure 9.8. Instantaneous Mach number contours for the SST-DES hybrid
turbulence model.
a. Fine b. Mid c. Coarse
Figure 9.9. Instantaneous eddy viscosity contours for the SST-DES hybrid
turbulence model.
![](asset-9a.png)
9-9
a. Fine b. Mid c. Coarse
Figure 9.10. Instantaneous Mach number contours for the SST-MS hybrid
turbulence model.
a. Fine b. Mid c. Coarse
Figure 9.11. Instantaneous eddy viscosity contours for the SST-MS hybrid
turbulence model.
All three hybrid turbulence models have similar trends on the three grid levels.
Significant turbulent structure can be seen in the fine grid solutions with both
large-scale and small-scale turbulent structure present. The mid level solutions
have large-scale structure present, but the smaller turbulent scales are absent.
The mid level grids also have structure that appears to be more periodic than
does the fine grid solution. The coarse level solutions show almost no turbulent
structure and produce an almost steady state wake away from the cylinder. The
level of eddy viscosity in the wake differs for the three turbulent models. The
eddy viscosity for the SA-DES model shown in Fig. 9.7 indicates that this model
effectively shuts off the eddy viscosity outside the boundary layer. The eddy
viscosity for the SST-DES and the SST-MS models is reduced as the grid is
refined as would be expected for these models. The SST-DES model tends to
predict higher levels of eddy viscosity in the wake region than does the SST-MS
model. The SST-DES model also predicts higher eddy viscosities along the
![](asset-9b.png)
9-10
edges of the wake than either the SA-DES or SST-MS models. Both the SST-
DES and SST-MS models are tending toward a RANS type solution in the far
wake of the coarse grid solution.
The ratio of the turbulent length scale to the local grid length scale for the SST-
MS model is shown in Fig. 9.12 for the three grid levels.
a. Fine b. Mid c. Coarse
Figure 9.12. Instantaneous ratio of turbulent length scale to the grid length scale
for the SST-MS hybrid turbulence model.
The length scale ratio seems to scale with the grid refinement for the mid to fine
grid results, which indicates that the turbulent length scales predicted on the mid
and fine grids are similar. The mid level grid results in a length scale ratio of
greater than two for most of the wake region. This is the desired performance
trend with grid refinement for hybrid turbulence models. The coarse grid results
do not seem to continue the scaling trend.
The power spectral densities (PSD) of the axial normal force coefficient are
shown in Figs. 9.13-9.15 for the three hybrid turbulence models. The mid and
fine level grid results are in general agreement for the SA-DES hybrid model.
The location and intensity of the primary spectral peak is still varying with the
SST-DES and SST-MS hybrid models. This may be due to the fact that the
SST-DES and SST-MS hybrid models use a ratio of the local turbulent length
scale to the local grid length scale to determine the subgrid eddy viscosity, while
the SA-DES model only uses the local grid length scale. The coarse grid
solutions are seen to be significantly different and the energy is contained in
narrower peaks for all of the hybrid models.
![](asset-9c.png)
9-11
Figure 9.13. Power spectral density of the normal force coefficient for the SA-
DES hybrid model.
Figure 9.14. Power spectral density of the normal force coefficient for the SST-
DES hybrid model.
![](asset-9d.png)
9-12
Figure 9.15. Power spectral density of the normal force coefficient for the SST-
MS hybrid model.
The average integrated drag coefficient (Cd), standard deviation of the lift
coefficient (σ(Cl)), and the peak lift coefficient Strouhal number (St) for the 3
hybrid models computed on the three grid levels are shown in Table 9.1. The
mid and fine level grid results are in general agreement with each other and the
experimental data9,10,11. The results are also in reasonable agreement with the
Ref. 12 results obtained at a Reynolds number of 3.0x106.
Model Grid Average Cd σ(Cl) St
SA Coarse 0.266 0.163 0.232
SST Coarse 0.192 0.103 0.244
SA-DES Coarse 0.298 0.253 0.232
SST-DES Coarse 0.216 0.171 0.250
SST-MS Coarse 0.196 0.115 0.250
SA Mid 0.511 0.428 0.281
SST Mid 0.517 0.520 0.281
SA-DES Mid 0.626 0.186 0.256
SST-DES Mid 0.589 0.274 0.293
SST-MS Mid 0.556 0.135 0.281
SA-DES Fine 0.518 0.182 0.281
SST-DES Fine 0.590 0.212 0.256
SST-MS Fine 0.585 0.238 0.268
DATA Ref. 9 0.505-0.540 0.0575-0.0770 0.306-0.308
DATA Ref. 10 0.79 - 0.27
DATA Ref. 11 - - 0.29
SA-DES Ref. 12 0.41-0.51 0.06-0.13 0.33-0.35
Table 9.1 Force coefficient and Strouhal number predictions on the 3D circular
cylinder.
![](asset-9e.png)
9-13
9.3 WICS Bay
Unsteady Navier-Stokes calculations were also performed for the WICS
(Weapons Internal Carriage and Separation) L/D=4.5 weapons bay13 for M=0.95
and Re=2.5x106/ft. The wind tunnel data were obtained in the AEDC four-foot
transonic wind tunnel. The weapons bay was 18 in. long, 4 in. wide, and 4 in.
deep. The computational geometry was a flat plate that extended 15 in.
upstream of the bay to match the experiment and 25 in. downstream of the bay.
The sides of the computational grid extended 50 in. on either side of the bay
centerline. The full bay geometry was modeled using wall functions. The wall
spacing was chosen as 0.0075 in., which corresponds to a y+ of 50 on the
upstream plate. This wall spacing was used for all the grids in the grid
refinement study.
Grid resolution effects were investigated by modifying the bay grid and the grids
in the vicinity of the bay grid. The grids upstream and to the sides of the bay
were not modified. Details of the grid systems and grid spacings used in this
study are shown in Table 9.2. The centerline plane of the three bay grids is
shown in Figure 9.16.
Grid Total
Points Bay Grid
Dimensions Bay Grid
Δxmax
Bay Grid
Δymax Bay Grid
Δzmax
Fine 1.8x10
6
121x61x61 0.3 in. 0.1 in. 0.1 in.
Medium 1.1x10
6
71x41x41 0.6 in. 0.2 in. 0.2 in.
Coarse 7.9x10
5
61x31x31 0.75 in. 0.3 in. 0.3 in.
Table 9.2. Parameters for the grid refinement study.
A time step study was performed using the SST-MS hybrid model and the
medium grid set. Time steps of 8.0x10-6 seconds to 8.0x10-5 seconds were
evaluated. A time step of 8.0x10-6 seconds corresponds to 260 time steps per
cycle for the primary spectral mode of the cavity at 480 hz. The calculations
were run 2000 steps to remove the initial transients and all unsteady results were
processed over the last 8192 time steps. The results are compared to data (Ref.
12) for the K16 and K18 transducer locations in the bay in Fig. 9.17. K16 is
located on the bay ceiling centerline 0.275 inches from the back wall. K18 is
located on the bay back wall centerline 0.725 inches from the bay opening.
The sound pressure levels predicted by the hybrid turbulence models includes
only the contribution of the grid resolved turbulent scales. Subgrid turbulent
energy, which includes most of the energy in the boundary layer, is not included.
Thus the spectrum predicted using the hybrid models should have lower energy
at higher frequencies (i.e. should have a more rapid roll off) than the
experimental data.
![](asset-9f.png)
9-14
Figure 9.16. Centerline plane of the bay grids.
Figure 9.17a. Sound pressure level spectra at the K16 location for varying time
step.
![](asset-a0.png)
9-15
Figure 9.17b. Sound pressure level spectra at the K18 location for varying time
step.
The agreement with the measured spectra is seen to improve at the higher
frequencies as the time step is reduced. The solution with the largest time step
predicts the primary spectral peaks at too low a frequency. The solution with the
smallest time step provides reasonable agreement for the first four spectral
peaks. The solution using 1.6x10-5 second time step has reasonable agreement
with the first three spectral peaks. All subsequent solutions were performed using
the 1.6x10-5 second time step.
A comparison of the sound pressure level for the three hybrid models, the
Spalart-Allmaras model, and the SST model for the medium grid are shown in
Fig. 9.18. The hybrid models show similar trends throughout the spectral range.
The two RANS models predict the first two spectral peaks, but are well below the
data away from the peaks. This is an indication that the RANS models are
providing too much damping of the unsteady solutions.
![](asset-a1.png)
9-16
Figure 9.18a. Sound pressure level spectra at the K16 location for the medium
grid.
Figure 18b. Sound pressure level spectra at the K18 location for the medium
grid.
Instantaneous Mach number and eddy viscosity contours on the bay centerline
are shown in Figs. 9.19-9.21.
![](asset-a2.png)
9-17
a. Fine b. Medium c. Coarse
Figure 9.19a. Instantaneous Mach number contours on the bay centerline using
the SA-DES hybrid model.
a. Fine b. Medium c. Coarse
Figure 9.19b. Instantaneous eddy viscosity contours on the bay centerline using
the SA-DES hybrid model.
a. Fine b. Medium c. Coarse
Figure 9.20a. Instantaneous Mach number contours on the bay centerline using
the SST-DES hybrid model.
![](asset-a3.png)
9-18
a. Fine b. Medium c. Coarse
Figure 9.20b. Instantaneous eddy viscosity contours on the bay centerline using
the SST-DES hybrid model.
a. Fine b. Medium c. Coarse
Figure 9.21a. Instantaneous Mach number contours on the bay centerline using
the SST-MS hybrid model.
a. Fine b. Medium c. Coarse
Figure 9.21b. Instantaneous eddy viscosity contours on the bay centerline using
the SST-MS hybrid model.
The turbulent structure shown in the Mach number contours is seen to
dramatically increase with grid resolution for all three hybrid models. The shear
layer is seen to grow and then burst as it traverses the bay. The eddy viscosity
decreases with increasing grid resolution as is expected for all three hybrid
models. The SST-DES hybrid model has the largest levels of eddy viscosity, and
the maximum eddy viscosity is located in the shear layer. The ratio of the
![](asset-a4.png)
9-19
turbulence length scale to the grid length scale (Lt/Lg) for the SST-MS model is
shown in Fig. 9.22. The contours scale with the grid refinement as is expected
for this model.
The time averaged pressure coefficient on the bay ceiling is shown in Fig. 9.23
for the three hybrid models. The sound pressure level on the bay ceiling is
shown in Fig. 9.24. The solutions are similar for the three models for both of
these properties. There is some sensitivity to grid density for the average
pressure coefficient on the downstream end of the ceiling. The medium and fine
grids produce similar results for the sound pressure. The coarse grid sound
pressure is only slightly different from the finer grids. This indicates that average
pressure coefficient and sound pressure level are not overly sensitive to grid
density.
Fine
Medium
Coarse
Figure 9.22. Instantaneous ratio of turbulent length scale to the grid length scale
contours on the bay centerline using the SST-MS hybrid model.
![](asset-a5.png)
9-20
Figure 9.23a. Average pressure coefficient on the WICS ceiling using the SA-
DES model.
Figure 9.23b. Average pressure coefficient on the WICS ceiling using the SST-
DES model.
![](asset-a6.png)
9-21
Figure 9.23c. Average pressure coefficient on the WICS ceiling using the SST-
MS model.
Figure 9.24a. Sound pressure level on the WICS ceiling using the SA-DES
model.
![](asset-a7.png)
9-22
Figure 9.24b. Sound pressure level on the WICS ceiling using the SST-DES
model.
Figure 9.24c. Sound pressure level on the WICS ceiling using the SST-MS
model.
The sound pressure level spectra at the K16 and K18 locations for the three
hybrid models are shown in Figs. 9.25-9.27. The medium and fine grid spectra
are in reasonable agreement with the first three spectral peaks for all three hybrid
![](asset-a8.png)
9-23
models. The coarse grid spectra are not quite as good, but are not totally
unacceptable. The fine grid solutions underpredict the primary spectral peak,
while the coarse and medium grids generally overpredict the primary spectral
peak. Spectral results seem to be more sensitive to grid density than the time
averaged quantities.
Figure 9.25a. Sound pressure level spectrum at the K16 location using the SA-
DES model.
Figure 9.25b. Sound pressure level spectrum at the K18 location using the SA-
DES model.
![](asset-a9.png)
9-24
Figure 9.26a. Sound pressure level spectrum at the K16 location using the SST-
DES model.
Figure 9.26b. Sound pressure level spectrum at the K18 location using the SST-
DES model.
![](asset-aa.png)
9-25
Figure 9.27a. Sound pressure level spectrum at the K16 location using the SST-
MS model.
Figure 9.27b. Sound pressure level spectrum at the K18 location using the SST-
MS model.
9.4 Delayed Detached Eddy Simulation (DDES)
Hybrid RANS/LES models were developed to provide a RANS solution in the
boundary layer along a body and to transition to LES away from the body. The
![](asset-ab.png)
9-26
functions used to transition from RANS to LES in these hybrid models can
transition too early and lead to solutions that are neither LES nor RANS. This
usually occurs in small pockets of separated flow such as a shock induced
separation on an airfoil. If the grid is fine enough and the turbulent intensities are
high enough in this region the hybrid model will begin to reduce the eddy
viscosity level. The reduced level may not be low enough to allow the flow to
become unsteady. Recently a correction has been proposed for the DES models
to force the models to remain turbulent in the boundary layer. The new
formulation is called Delayed Detached Eddy Simulation (DDES).
The DDES14 version of the SA model modifies Eq. (9.1) to force the modified wall
distance parameter to equal the distance from the wall further into the boundary
layer. This delays the transition from RANS to LES. The new modified distance
function is given by
()
Δ−−= DESd Cdfdd ,0max
~
(9.12)
where
[]
(
)
3
8tanh1 dd rf −≡ (9.13)
and
22 dS
rt
d
κ
νν
+
≡ (9.14)
fd is designed to be 1 in the LES region where rd << 1, and 0 elsewhere. rd can
be thought of as a ratio of the turbulent length scale defined by S
t/
ν
and the
wall distance d. This length scale ratio will work for simple turbulent flows, but
may have difficulty in more complicated flows since the eddy viscosity is a
transported quantity while the strain is a locally derived quantity. It is not difficult
to imagine flows where a large value of eddy viscosity is transported into a region
of relatively small strain. This could cause the DDES model to transition from
LES to RANS mode away from the body.
A DDES version of the SST model15 can be derived by replacing Lg in Eq. (9.3)
with a wall corrected grid length scale Lgcor defined as
(
)
gtDtgcor LLfLL −−= (9.15)
where
![](asset-ac.png)
9-27
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡⎟
⎠
⎞
⎜
⎝
⎛
−=
3
5.1tanh1 d
L
ft
D
(9.16)
Here d is the distance to the wall. The turbulent dissipation in the turbulent
kinetic energy equation is effectively increased when the grid size length scale Lg
is less than the turbulent length scale Lt. The wall correction increases the size
of Lgcor near the wall and causes the turbulence model to remain in RANS mode
in the near wall region without regard to the local grid size in the boundary layer.
Note that the turbulent length scale used here is based on the ratio of the
transported turbulent quantities k and
ω
only.
The pressure distribution for the transonic flow over a NACA 001216 wing test can
be used to show the benefit of the DDES correction. The free-stream Mach
number was 0.799, the angle-of-attack was 2.26 degrees, and the Reynolds
number based on chord was 9x106. A 381x81 two dimensional grid was used for
this study with 301 points on the airfoil.
The pressure coefficient on the airfoil is shown in Fig. 9.28. All of the solutions
were run in an unsteady mode, but all of the models produced steady state
solutions for this case. The SA and SST RANS model produce similar results for
this case. The SA DES model moves the shock on the upper surface upstream
and produces a weak shock on the lower surface. This is due to a decrease in
eddy viscosity in the boundary layer. The SA DDES model produces a solution
that is closer to the SA model result. The SST DES model moves the shock well
upstream on the upper surface and produces a shock on the lower surface. The
SST DDES model produces a solution similar to the SA DDES model.
![](asset-ad.png)
9-28
Figure 9.28a Pressure coefficient for a NACA 0012 airfoil for SA models.
Figure 9.28b Pressure coefficient for a NACA 0012 airfoil for SST models.
The SA DES and SA DDES model modified wall distance ( d
~
) is shown in Fig.
9.29 for the separated flow region behind the shock at x/c=0.75 on the wing
upper surface. The boundary layer edge is approximately at z/c=0.045. Note
that dd =
~
for the SA model. The SA DDES model is seen to delay the transition
to LES to a much larger distance from the wall than does the SA DES model.
The effective turbulent length scale used for transitioning the SST DES and SST
DDES models at the same location is shown in Fig. 9.30. The SST DES is seen
to transition very near the wall. The SST DDES modified turbulent length scale
moves the RANS to LES transition much further away from the wall.
![](asset-ae.png)
9-29
Figure 9.29 Wall distance parameter for the SA, SA DES, and SA DDES
turbulence models.
L
Figure 9.30 Turbulent length scale for the SST, SST DES and SST DDES
turbulence models.
![](asset-af.png)
9-30
9.5 Summary
These simulations indicate a fundamental difficulty in verification and validation
for unsteady flows. As the grid is refined, smaller scale turbulent structures are
resolved in the solution. This process will continue until the grid is refined below
the Komologrov scale and all of the turbulent scales are grid resolved. The
limiting grid refinement case cannot be approached with current computing
hardware for most real world high Reynolds problems of interest. Comparison of
statistical quantities derived from the unsteady solutions can be made to assure
that the important features of the unsteady flow have been captured with a given
grid, but the nature of the statistical analysis makes it difficult to assess whether
a true “grid convergence” has been achieved.
The time step study indicated that about 200 time steps per primary shedding
cycle were required for temporal accuracy with the SST-MS hybrid model. Using
a larger time step causes the primary spectral peak to occur at a lower frequency
than predicted with the smaller time step and than is seen in the data.
The simulations on the mid and fine grid systems are in good general agreement
indicating that a level of grid convergence can be achieved for the large turbulent
scales. The turbulent length scale to grid length scale ratio was greater than two
in the wake region of the cylinder and in the bay. This may serve as a rule of
thumb for grid resolution for hybrid model applications. Although the primary
shedding frequency for the cylinder is approximately right for the coarse grid, the
integrated forces and the spectral shape are quite different from the mid and fine
grid solutions. This indicates that these turbulence models can produce solutions
that appear to capture the relevant physics, but not capture the physical details of
the flow. The coarse grid solutions were totally non-physical.
The hybrid RANS/LES turbulence models are relatively new and will need to be
exercised for a wide variety of problems to determine their accuracy before they
become an accepted tool for fluid modelers. They seem to offer much for
unsteady flow applications, but issues such as grid sensitivity need to be further
addressed. Hopefully more effort will go into these models in the near future so
that they can be matured for use in everyday applications.
![](asset-b0.png)
9-31
Chapter 9 References:
1. Spalart, P., Jou, W-H. Strelets, M., and Allmaras, S. "Comments on the
Feasibility of LES for Wings and on a Hybrid RANS/LES Approach," First
AFOSR Conference on DNS/LES, August 1997, C. Liu and Z. Liu editors,
Greyden Press, Columbus, Ohio.
2. Constatinescu, G. and Squires, K., “LES and DES Investigations of the
Turbulent Flow over a Sphere,” AIAA-2000-0540, Jan. 2000.
3. Mitchell, A., Morton, S., and Forsythe, J., “Analysis of Delta Wing Vortical
Substructures using Detached-Eddy Simulation,” AIAA-2002-2968, Jun.
2002.
4. Forsythe, J., Squires, K., Wurtzler, K., and Spalart, P, “Detached-Eddy
Simulation of Fighter Aircraft at High Alpha,” AIAA-2002-0591.
5. Strelets, M., "Detached Eddy Simulation of Massively Separated Flows,"
AIAA-2001-0879, Jan. 2001.
6. Nichols, R. and Nelson, C., “Applications of Hybrid RANS/LES Turbulence
Models,” AIAA-2003-0083, Jan. 2003.
7. Nelson, C. and Nichols, R., “Evaluation of Hybrid RANS/LES Turbulence
Models using an LES Code,” with C. C. Nelson, AIAA-2003-3552, Jun. 2003.
Hybrid RANS/LES Application Hints
1. These models are very new and have not been used extensively. They
should be used with care and only when validation calculations can be
performed.
2. These models are for unsteady applications, and should not be used with local
time stepping or with other non-time accurate algorithms.
3. Turbulent flows are three dimensional, and hence these models should be
used only in 3D (see Ref. 17).
4. Because of the unsteady nature of these models, they may require a large
number of time-steps to obtain a statistically stationary solution for analysis.
5. These models may be sensitive to the computational mesh because the filter
function is inversely proportional to the grid spacing. A rule-of-thumb is that
the ratio of the turbulent scale to grid length scale should be greater than two
in the region of interest when using hybrid models.
6. As with all unsteady applications, care should be taken to be sure the time
step is small enough to temporally capture the unsteady phenomena of
interest.
7. The models can produce solutions that are neither RANS or LES when
applied to flows with small pockets of flow separation such as shock induced
separation on an airfoil.
9-32
8. Nichols, R., “Comparison of Hybrid RANS/LES Turbulence Models for a
Circular Cylinder and a Cavity,” AIAA Journal, Vol. 44, No. 6, June 2006, pp.
1207-1219.
9. Jones, G. W., Cincotta, J. J., and Walker, R. W., ”Aerodynamic Forces on a
Stationary and Oscillating Circular Cylinder at High Reynolds Number,”
NASA-TR-R-300, October 1968.
10. Roshko, A., “Experiments on the Flow Past a Circular Cylinder at Very High
Reynolds Number,” Journal of Fluid Mechanics, Vol. 10, Part 3, May 1961,
pp. 345-356.
11. Schlichting, H., Boundary Layer Theory, 7th Edition, McGraw-Hill, New York,
1979.
12. Travin, A., Shur, M., Strelets, M. and Spalart, P., “Detached-Eddy
Simulations Past a Circular Cylinder,” Flow, Turbulence and Combustion,
Vol. 63, pp 293-313, 1999.
13. Dix, R. E. and Bauer, R. C., “Experimental and Theoretical Study of Cavity
Acoustics,” AEDC-TR-99-4, May 2000.
14. Spalart, P. R., Deck, S., Shur, M. L., Squires, K. D., Strelets, M. Kh., and
Travin, A., “A New Version of Detached-Eddy Simulation, Resistant to
Ambiguous Grid Densities,” Theor. Comput. Fluid Dyn., Vol. 20, 2006, pp.
181-195.
15. Nichols, R. H., “A Comparison of Hybrid RANS/LES Turbulence Models for a
Generic Weapons Bay With and Without a Spoiler,” AIAA-2008-6629, Aug.
2008.
16. Maksymiuk, C. M. and Pulliam, T. H., “Viscous Transonic Airfoil Workshop
Results Using ARC2D,” AIAA Paper 87-0415, Jan. 1987.
17. Shur, M., Spalart, P., Squires, K., Strelets, M., and Travin, A., “Three
Dimensionality in Reynolds-Averaged Navier-Stokes Solutions around Two-
Dimensional Geometries,” AIAA Journal, Vol. 43, No. 6, June 2005, pp.
1230-1242.
![](asset-b2.png)
10-1
10.0 Wall Function Boundary Conditions
Turbulent transport models that are applicable all the way to the wall are called
low Reynolds number turbulence models since they are valid at low turbulent
Reynolds numbers (
ρ
k2/(
με
)). Low Reynolds number models require tight grid
spacing near the wall to resolve the large gradients in velocity, temperature, and
turbulent quantities near the wall. This can lead to large grid requirements for 3D
Navier-Stokes applications. The stability of a numerical algorithm is limited by
the smallest cell size in a grid, so the small cells near the wall can severely
restrict the maximum allowable time step for a problem. The near-wall
turbulence damping terms are generally expressed in exponential functions. This
can lead to difficulties in converging the numerical solution in the near wall
region.
Turbulence models that do not include wall correction terms in the differential
equations are called high Reynolds number turbulence models. High Reynolds
number turbulence models rely on empirically derived algebraic models of the
near-wall region of the boundary layer to provide boundary condition information
to the mean flow Navier-Stokes equations at the first point off the wall. These
empirically derived relationships are called wall functions. Since the high
gradient region near the wall is modeled with these empirical relationships, the
first point off the all may be placed much farther away than with low Reynolds
number models. This reduces the number of points required to discretize a flow
field and increases the maximum allowable time step.
Wall functions have been used by turbulence modelers for a number of years.
Wall functions were a matter of necessity during the years before computers
began to mature, and are still a convenient and efficient means of modeling
many turbulent flows. Most wall function implementations are based on empirical
relations for incompressible adiabatic flows. Some implementations include a
decoupled relation for the temperature distribution in the lower portion of the
boundary layer that allows the wall functions to be applied to flows with heat
transfer at the wall. The wall function boundary condition formulation developed
here is based on a coupled velocity and temperature boundary layer profile that
is applicable to incompressible flows, compressible flows, and flows with heat
transfer. The wall functions described here are compatible with most transport
type turbulence models.
10.1 Theory
There are six fundamental assumptions used in the development of wall function
boundary conditions for compressible flows:
1. Analytical expressions are available for the velocity and temperature
profiles in the lower part of the boundary layer.
![](asset-b3.png)
10-2
2. Analytical expressions are available for the turbulent transport
variables at the first point off the wall.
3. The pressure is constant in the lower part of the boundary layer.
4. The shear stress (
()
y
u
tw ∂
∂
+==
μμττ
) is constant in the lower part of
the boundary layer.
5. The heat transfer (
()
y
T
kkqq tw ∂
∂
+== ) is constant in the lower part of
the boundary layer.
6. There are no chemical reactions (i.e. the chemistry is frozen) in the
lower part of the boundary layer.
Several empirical relationships for the velocity in a boundary layer are available
for use in wall function boundary conditions. A typical boundary layer velocity
profile is shown in Fig. 3.2. The simplest profile form is the adiabatic
incompressible universal law-of-the-wall
()
Byu += ++ log
1
κ
(10.1)
where u+ and y+ are defined in Eqs. 3.10 and 3.11. The universal law-of-the-
wall is valid for the y+ range of about 30<y+<1000. The velocity in the sublayer
(0<y+<10) is given by
++ =yu (10.2)
The sublayer must be blended with the law-of-the-wall to cover the interim
region. Spalding1 suggested a unified form valid for the log layer and the
sublayer as well as the transition region. This form is given by
(
)
(
)
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡−−−−+=
++
+−++ +
6
u
2
u
u1eeuy
32
uB
κκ
κ
κκ
(10.3)
Note that the +
−uBee
κκ
term in Eq. 10.3 is a restatement of the incompressible
adiabatic law-of-the-wall (Eq. 10.1). The constants
κ
and B are generally taken
as 0.4 and 5.5 respectively. The wall function boundary condition can be tuned
for such effects as surface roughness by tuning the
κ
and B constants (Ref. 2).
Compressibility and heat transfer effects can be included by replacing the
incompressible adiabatic law-of-the-wall term in Spalding’s equation (Eq. 10.3)
with the outer velocity form of White and Christoph3. White and Christoph’s outer
velocity is given by
![](asset-b4.png)
10-3
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
Γ
+Φ+
Γ
=+
+
+
0
lnsin
2
1
y
y
Qu
κ
β
(10.4)
where
()
()
3/1
0
2
1
2
1
www
ww
wp
2
P
r
r
Bexpy
4Q
Q
sin
ukT
q
Tc2
ru
=
−=
Γ+=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
=
=
=Γ
+
−
κ
β
β
ϕ
ρ
μ
β
τ
τ
(10.5)
The nondimensional parameter
Γ
models compressibility effects and the
parameter
β
models heat transfer effects. The outer velocity form of White and
Christoph reduces to the standard law-of-the-wall (Eq. 10.1) for incompressible
adiabatic flow. Eq. 10.3 can then be written as
(
)
(
)
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡+++−+=
++
+−+++
6
u
2
u
u1eyuy
32
B
White
κκ
κ
κ
(10.6)
where
()
B
Q
u
yWhite
κφ
βκ
−
⎪
⎭
⎪
⎬
⎫
⎪
⎩
⎪
⎨
⎧⎥
⎦
⎤
⎢
⎣
⎡−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−Γ
Γ
=
+
−+ exp
2
sinexp 1 (10.7)
The temperature distribution within the boundary layer is given by the Crocco-
Busemann equation
()
(
)
2
wuu1TT ++ Γ−+=
β
(10.8)
For adiabatic wall cases,
β
=0 and the Crocco-Busemann equation reduces to
()
⎟
⎠
⎞
⎜
⎝
⎛−
+= 2
wu
2
r1
1TT
γ
(10.9)
![](asset-b5.png)
10-4
The wall shear stress for adiabatic flows is determined as follows:
1. Set the velocity at the wall to zero for non-moving body problems or to
the grid velocity for moving body problems.
2. Given the velocity and temperature at the first point off the wall, solve
Eq. 10.9 for the wall temperature (Tw).
3. Extrapolate the pressure from the first point off the wall and solve for
the wall density (ρw) using the equation of state.
4. Iteratively solve Eq. 10.6 for the wall shear stress in the wall coordinate
system.
5. Rotate the stress tensor into the computational coordinate system and
replace the viscous flux calculation at the wall for cell-centered
algorithms or at the half node for node-centered algorithms.
The wall shear stress and the heat transfer for constant temperature wall flows
are determined by:
1. Set the velocity at the wall to zero for non-moving body problems or to
the grid velocity for moving body problems.
2. Extrapolate the pressure from the first point off the wall and solve for
the wall density using the equation of state with the given wall
temperature.
3. Iteratively solve Eq. 10.6 and Eq. 10.8 for the wall shear stress and
heat transfer in the wall coordinate system.
4. Rotate the stress tensor and the heat transfer vector into the
computational coordinate system and replace the viscous flux
calculation at the wall for cell-centered algorithms or at the half node
for node-centered algorithms.
One approach to introducing the wall function corrected wall shear stress and
heat transfer into the calculation of the viscous fluxes has been to calculate an
“effective wall eddy viscosity” so that the discrete shear stress at the first half
point off the wall given by
12
12
2t1t
21
5.1 yy
uu
22 −
−
⎟
⎠
⎞
⎜
⎝
⎛+
+
+
=
μμ
μμ
τ
(10.10)
yields the correct value of the wall shear stress. Unfortunately this approach
introduces errors into the energy equation since the slope of the temperature
derivative will not be calculated correctly (i.e. a separate effective wall eddy
viscosity would be required for the temperature).
Sondak and Pletcher4 describe a procedure to perform a transformation of the
stresses for generalized coordinate systems. Although this transformation is
much more general than the effective wall eddy viscosity approach and
introduces no error in the energy equation, it is also much more complex. In their
![](asset-b6.png)
10-5
transformation the wall function generated wall shear stress is assumed to be in
the wall coordinate system, i.e. a system aligned with velocity vector at the first
point off the wall and the normal vector to the wall. This is not always the case
for separated flows, but the assumption seems to produce reasonable results
even when the flow is separated.
Another method of correcting the computational shear stress at the wall with the
wall function value is to calculate the ratio of the magnitude of the computed wall
shear stress to the wall function value and then scaling the computed values by
this ratio. This is equivalent to performing the transformation of Sondak and
Pletcher4, and is much simpler to apply. The temperature derivative at the wall
may be treated in a similar manner.
Once the wall shear stress has been determined the turbulence transport
variables must be determined at the first point off the wall. The eddy viscosity at
the first point off the wall for incompressible adiabatic flow can be found using the
constant stress assumption to be
()
()
(
)
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−−−−=
+
++
2
u
u1uexpBexp
2
w
t
κ
κκκκ
μ
μ
(10.11)
The equivalent form of the eddy viscosity for compressible flows with heat
transfer can also be found using the constant stress assumption. The eddy
viscosity is given by
(
)
w
1w
2
white
w
t
2
u
u1)Bexp(
y
y
1
μ
μ
κ
κκκ
μ
μ
+
+
+
+
+
−
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛++−−
∂
∂
+= (10.12)
where 1w+
μ
is the molecular viscosity at the first point off the wall and +
+
∂
∂
y
ywhite is
given by
()
[]
2/1
2
2
white
white Q/u21
Q
y2
y
y
β
κ
−Γ−
Γ
=
∂
∂++
+
+
(10.13)
Eq. 10.12 reduces to Eq. 10.11 for adiabatic incompressible flows. The values of
the transport model turbulence variables at the first point off the wall must also be
defined. The values for the Spalart variable
ν
~
is given by
3
1
34
~
~
vtt C
μνμν
+= (10.14)
The turbulent kinetic energy and turbulent dissipation at the first point off the wall
for k-
ε
models are given by
![](asset-b7.png)
10-6
μ
τ
C
u
k2
= (10.15)
t
kC
μ
ρ
εμ
2
= (10.16)
The turbulent kinetic energy and specific turbulent dissipation at the first point off
the wall for k-
ω
models are given by
2
w
w
iy075.0
6
ρ
μ
ω
= (10.17)
yC
u
o
κ
ϖ
μ
τ
= (10.18)
2
o
2
i
ϖϖϖ
+= (10.19)
ρ
ϖμ
t
k= (10.20)
The two level model for
ω
in Eqs. 8.17-8.20 was suggested by Veiser, Esch, and
Menter5.
Most upwind flow solvers require no further modifications than the inclusion of
the corrected wall stress and heat transfer in the viscous flux calculations. The
smoothers used with central difference algorithms require some additional
modifications. The fourth-order smoother used to maintain numerical stability will
produce excessive dissipation at the wall when wall functions are employed
because of the large velocity difference between the wall and the first point off
the wall. To prevent this from occurring, the smoother should be modified as
follows. The standard form of the explicit fourth-order smoother in the
ξ
computational direction is
[
]
ξξξ ψ
ΛΔ∇Δ )4(
t (10.21)
where
()
1n
j
qJ +
Δ∇=Λ r
ξξξ
(10.22)
![](asset-b8.png)
10-7
where Δt is the time step, Δ
ξ
is the central difference operator, and
ψ
(4) is the
explicit fourth-order smoothing coefficient. The second derivative term Λ
ξ
is
normally set to zero at the wall for viscous calculations. This term is also set to
zero at the first point off the wall when wall functions are used. This effectively
turns the fourth-order smoothing off at the first point off the wall and reduces it
significantly at the second point off the wall.
Once the wall shear stress is calculated using wall functions, the y+ value of the
first point off the wall is known. It is a simple matter to automatically switch
between the wall functions and integrating to the wall based on the local value of
y+. The turbulence model chosen should contain the low Reynolds number terms
required for the model to be valid in the near wall region if automatic switching is
implemented. Automatic switching offers advantages in complex configurations
since computational grid points can be saved by using the larger wall function
wall spacings on non-essential portions of the geometry and integrating to the
wall where high quality skin friction or heat transfer is required. Automatic
switching is also useful in grid sequencing or multigrid algorithms where the wall
functions can provide an improved estimate of the wall shear stresses and heat
transfer on the coarse mesh solutions.
Two examples of wall function applications for the Spalart-Allmaras and the SST
turbulence models are included here. Other example of wall function
applications for structured and unstructured grids are given in Ref. 6 and Ref. 7
respectively.
10.2 Grid Sensitivity for a Flat Plate with Adiabatic Walls
The predicted and theoretical values for the skin friction distribution on an
adiabatic flat with varying initial grid wall spacings using are shown in Fig. 10.1
for the Spalart-Allmaras turbulence model and in Fig. 10.2 for the SST turbulence
model. The y+=1 results did not use wall functions and are included for
reference. All the results shown here used a grid-stretching ratio of 1.2. The wall
function predictions of skin friction are very similar for both turbulence models.
The models are in reasonable agreement for length Reynolds numbers (Rex)
above 1.0x106. Velocity profiles for the two turbulence models are shown in Fig.
10.3 and 10.4. Again the y+=1.0 results were not run with wall functions and are
included for reference. The velocity profiles are in reasonable agreement for all
of the wall spacings run here.
![](asset-b9.png)
10-8
Figure 10.1 Flat plate skin friction predictions for the Spalart-Allmaras turbulence
model with wall functions for varying initial wall grid point spacings.
Figure 10.2 Flat plate skin friction predictions for the SST turbulence model with
wall functions for varying initial wall grid point spacings.
![](asset-ba.png)
10-9
Figure 10.3 Flat plate velocity profile predictions for the Spalart-Allmaras
turbulence model with wall functions for varying initial wall grid point spacings.
Figure 10.4 Flat plate velocity profile predictions for the SST turbulence model
with wall functions for varying initial wall grid point spacings.
![](asset-bb.png)
10-10
10.3 Grid Sensitivity for an Axisymmetric Bump
A second application of wall functions for adiabatic flows with pressure gradients
is the Ames axisymmetric bump described in Chapter 5. Both the Spalart-
Allmaras and the SST models were run with varying initial grid wall spacings.
Results for the pressure coefficient for the Spalart-Allmaras model and the SST
model are shown in Fig. 10.5 and 10.6 respectively. Results for the velocity
distribution at x/c=1.0 (the trailing edge of the bump) are shown for the Spalart-
Allmaras and the SST models in Fig. 10.7 and 10.8 respectively. The y+=1.0
results were not run with wall functions and are included for reference. For both
the pressure coefficient distribution and the velocity profile at x/c=1.0 the results
for each model are similar. The results tend to start to diverge for y+=200. The
wall functions provide reasonable results in the shock induce separated flow
region.
Figure 10.5 Pressure coefficient predictions for the Ames axisymmetric bump for
the Spalart-Allmaras turbulence model with wall functions for varying initial wall
grid point spacings.
![](asset-bc.png)
10-11
Figure 10.6 Pressure coefficient predictions for the Ames axisymmetric bump for
the SST turbulence model with wall functions for varying initial wall grid point
spacings.
Figure 10.7 Velocity distribution at x/c=1.0 for the Ames axisymmetric bump for
the Spalart-Allmaras turbulence model with wall functions for varying initial wall
grid point spacings.
![](asset-bd.png)
10-12
Figure 10.8 Velocity distribution at x/c=1.0 for the Ames axisymmetric bump for
the SST turbulence model with wall functions for varying initial wall grid point
spacings.
10.4 Grid Sensitivity for a Flat Plate with Heat Transfer
Calculating heat transfer accurately can be more difficult than predicting skin
friction. This can be seen in the subsonic flat plate example when the wall
temperature is specified to be 1.5 times the free-stream temperature. The
sensitivity of the skin friction and heat transfer result with varying initial grid wall
spacing is shown in Fig. 10.9 and Fig. 10.10 for the Spalart-Allmaras model and
in Fig. 10.11 and 10.12 for the SST model. The y+=0.1 results were not run with
wall functions and are included for comparison. The grid stretching ratio was
fixed at 1.2 for these results. Both the skin friction and heat transfer seem to be
relatively insensitive to the wall spacing for the wall spacings presented here.
![](asset-be.png)
10-13
Figure 10.9 The effect of wall spacing on the skin friction on a flat plate with heat
transfer using the Spalart-Allmaras turbulence model.
Figure 10.10 The effect of wall spacing on the heat transfer (Stanton number) on
a flat plate using the Spalart-Allmaras turbulence model.
![](asset-bf.png)
10-14
Figure 10.11 The effect of wall spacing on the skin friction on a flat plate with
heat transfer using the SST turbulence model.
Figure 10.12 The effect of wall spacing on the heat transfer (Stanton number) on
a flat plate using the SST turbulence model.
![](asset-c0.png)
10-15
Profiles of velocity and temperature at a length Reynolds number (Rex) of 1x107
are shown in Fig. 10.13 and 10.14 respectively for the Spalart-Allmaras model
and in Fig. 10.15 and 10.16 for the SST model. Again the y+=0.1 results were
not run with wall functions and are included for comparison purposes. The wall
function results are in reasonable agreement for all wall spacings investigated for
both turbulence models.
Figure 10.13 The effect of wall spacing on the velocity profile on a flat plate with
heat transfer using the Spalart-Allmaras turbulence model and wall functions.
![](asset-c1.png)
10-16
Figure 10.14 The effect of wall spacing on the temperature profile on a flat plate
with heat transfer using the Spalart-Allmaras turbulence model and wall
functions.
Figure 10.15 The effect of wall spacing on the velocity profile on a flat plate with
heat transfer using the SST turbulence model and wall functions.
![](asset-c2.png)
10-17
Figure 10.16 The effect of wall spacing on the temperature profile on a flat plate
with heat transfer using the SST turbulence model and wall functions.
10.5 Grid Sensitivity for a Nozzle with Heat Transfer
Flow through a supersonic nozzle with a constant temperature wall can serve as
a test case for evaluating the performance of the wall function formulation in the
presence of strong pressure gradients. Details of the geometry and boundary
conditions for the converging-diverging supersonic nozzle are given in Chapter 5.
High-pressure air was heated by the internal combustion of methanol and flowed
along a cooled constant area duct before entering the nozzle. The gas could be
treated as a calorically perfect gas with a ratio-of-specific heats (γ) of 1.345. The
nozzle exit Mach number was 2.5. The molecular viscosity and thermal
conductivity were assumed to vary according to Sutherland’s law.
The grid stretching ratio used in this study was 1.2. Comparisons of the
predicted and measured pressure along the nozzle are shown in Fig. 10.17 for
the Spalart-Allmaras model and in Fig. 10.18 for the SST model. The results for
a wall spacing of y+=0.1 were not run with wall functions and are included for
reference. The pressure distribution is seen to be insensitive to the wall spacing.
The wall heat transfer is shown in Fig. 10.19 for the Spalart-Allmaras model and
in Fig. 10.20 for the SST model. The SST results for heat transfer are much less
sensitive to the initial wall spacing than are the Spalart-Allmaras results. The
heat transfer predictions are adequate for many applications even for an initial
wall spacing of y+=100.
![](asset-c3.png)
10-18
Figure 10.17 The effect of wall spacing on the pressure distribution for a
supersonic nozzle with heat transfer using the Spalart-Allmaras turbulence
model and wall functions.
Figure 10.18 The effect of wall spacing on the pressure distribution for a
supersonic nozzle with heat transfer using the SST turbulence model and wall
functions.
![](asset-c4.png)
10-19
Figure 10.19 The effect of wall spacing on the heat transfer distribution for a
supersonic nozzle with heat transfer using the Spalart-Allmaras turbulence
model and wall functions.
Figure 10.20 The effect of wall spacing on the heat transfer distribution for a
supersonic nozzle with heat transfer using the SST turbulence model and wall
functions.
![](asset-c5.png)
10-20
Chapter 10 References:
1. White, F. M., Viscous Fluid Flow, McGraw-Hill, Inc., ISBN 0-07-069710-8,
1974.
2. Rodi, W., “Turbulence Models and Their Application in Hydraulics – A State
of the Art Review,” International Association for Hydraulic Research, 2nd
Ed., Feb. 1984.
3. White, F. M. and Christoph, G. H., “A Simple Analysis of Compressible
Turbulent Two-dimensional Skin Friction Under Arbitrary Conditions,”
AFFDL-TR-70-133, Wright Patterson AFB, OH, Feb. 1971.
4. Sondak, D. L. and Pletcher, R. H., “Application of Wall Functions to
Generalized Nonorthogonal Curvilinear Coordinate Systems,” AIAA Journal,
Vol. 33, No. 1, pp. 33-41, January 1995.
5. Veiser, W., Esch, T., and Menter, F., “Heat Transfer Predictions using
Advanced Two-Equation Turbulence Models,” CFX Technological
Memorandum CFX-VAL10/0602, Nov. 2002.
6. Nichols, R. H., and Nelson, C. C., “Wall Function Boundary Conditions
Including Heat Transfer and Compressibility,” AIAA Journal, Vol. 42, No. 6,
pp 1107-1114, June 2004.
7. Frink, N. T., “Assessment of an Unstructured-Grid Method for Predicting 3-D
Turbulent Viscous Flows,” AIAA-96-0292, Jan. 1996.
Wall Function Application Hints
1. The first point off the wall should not exceed y+=100. A good rule of thumb
is to place the first point off the wall at y+=50.
2. Wall functions should not generally be used for calculations requiring
highly accurate friction drag or heat transfer.
3. Note that the computational domain for the turbulence variables does not
include the first point off the wall when using wall functions since it is
prescribed by the boundary condition.
4. Post processing of terms that depend on the velocity or temperature
gradient at the wall must include the wall function in order to match the wall
velocity or temperature gradient used in the calculation.
![](asset-c6.png)
11-1
11.0 Boundary Layer Transition Simulation
The transition of a boundary layer from laminar to turbulent impacts the
characteristics of a flow field, but its underlying physics has yet to be well
understood. The lack of a validated boundary layer transition simulation
capability is a major source of error for many applications, such as Reynolds
number scaling of wind tunnel results to flight, hypersonic flight, and high altitude
turbine engine propulsion. In spite of the large amount of turbulence transition
data available, no empirical method has been formulated that can reliably predict
transition for a variety of flight conditions and geometries. Though current efforts
using Direct Numerical Simulations are interesting and have shown promising
results, their use in production applications on real geometries is still years away.
One method for predicting the onset of transition is the use of the boundary layer
stability theory described in Chapter 1. This methodology has a number of
shortcomings:
1) The method cannot be used for simulating bypass transition.
2) The method can only be used to predict the onset of transition and not for
modeling the transition region or the fully developed turbulent flow region.
There is no way to directly use the results of this method to obtain a
solution for the entire transitional flow field.
3) The method requires a highly converged steady-state laminar flow solution
to the Navier-Stokes equations. This can be difficult to obtain for many
configurations, especially if large scale flow separation is present.
Probably the most popular method for simulating transition is the use of
Reynolds-Averaged Navier-Stokes (RANS) CFD codes with a low Reynolds
number turbulence model. The RANS equations are solved either directly or in
conjunction with some empirical correlations. The use of low Reynolds number
RANS models has proven unreliable in predicting the change in skin friction and
heat transfer within the transition region. No model of this type performs
satisfactorily under the influence of free-stream turbulence intensity and pressure
gradients. It is extremely difficult to obtain the correct location of the onset of
transition with this class of models. Recently methods have been developed for
simulating boundary layer transition that utilize additional transport equations for
intermittency or for a disturbance kinetic energy. These models rely heavily on
empirical functions to predict the onset and extent of the transition region.
Transport based transition models are of current interest because they can be
easily coupled to existing transport equation turbulence models and can be used
on complex configurations.
Turbulence transition models have been broadly categorized in two groups:
models based on stability theory and models not based on stability theory.
Models not based on stability theory are further divided into two groups: models
with specified transition onset and models with onset prediction capability.
![](asset-c7.png)
11-2
11.1 Transition Models Based on Stability Theory
The en method proposed by Smith and Gamberoni1 and Van Ingen2, based on
linear stability theory, is one of the most popular methods available for transition
prediction. There are three steps in the application of the en method. The first
step involves the computing of the laminar velocity and temperature profiles at
different stream wise locations for the given flow. In the second step, the
amplification rates of the most unstable waves are calculated for each profile by
using the en method. In the third step, these amplification rates are used to
calculate the transition location. There are several problems with the en method.
The major criticisms that the en method has received is that it was developed
based on the linear stability theory with an assumption that the flow is locally
parallel. The value of the ‘n’ factor for transition is not universal and needs to be
determined based on experimental data. This value varies from one wind tunnel
to the next.
The linear Parabolized Stability Equations (PSE) method addresses the non-
parallel effects neglected in the linear stability theory and assumes that the mean
flow, amplitude functions and wave number (α) are dependent upon the stream
wise distance (x). A further development of the linear PSE, known as the
nonlinear PSE, incorporates the nonlinear effects that have been neglected in the
linear stability theory. Methods based on the stability theory have one major
drawback - they need to track the growth of the disturbance amplitude along a
streamline. This limitation poses a significant problem for three-dimensional flow
simulations where the streamline direction is not aligned with the grid. Coupling
of such methods with CFD codes requires an unrealistically high grid density to
yield the boundary layer data with the required level of accuracy. These methods
also require a well-converged steady-state solution, which may not be obtainable
for real-world problems involving local flow separation. The main advantage of
these methods is that they give the correct treatment of the surface curvature.
Some different techniques have been employed to use these stability-based
methods more efficiently3. One method is to generate a database of the solution
of the linear stability equation for different velocity profiles in advance. The local
flow stability can then be determined quickly based on the local velocity
calculated from CFD codes. The validity of these models is limited to the range of
velocity profiles available in the database.
11.2 Transition Models With Specified Transition Onset
The transition region models in this section are unable to predict the location of
the transition. The transition location is determined from empirical data or results
from an en computation. The transition region is modeled by modifying existing
turbulence models. In Ref. 4, six transition models were implemented into the
commercial Navier-Stokes code GASP. These six models were the Baldwin-
![](asset-c8.png)
11-3
Lomax model, the Wilcox k-ω model, the Schmidt and Patankar low-Re k-ε
model which had a production term modified for modeling transition, the Warren,
Harris and Hasan one-equation model, the algebraic transition model developed
at ONERA/CERT, and the linear combination transition model developed by Dey
and Narasimha. These models were used to simulate hypersonic experimental
cases that included transition on a cone at Mach 65, a compression ramp at
Mach 10.086, and five flared cone test cases at Mach 7.937,8. Out of the five
flared cones used there were two with favorable pressure gradients, two with
adverse gradients, and one with a zero pressure gradient.
1) Baldwin-Lomax algebraic turbulence model:
The Baldwin-Lomax model was used to predict the transition region by turning off
the turbulence model for the laminar region by setting the eddy viscosity equal to
zero and then just turning it on at the transition point. It was found that in most of
the cases this model adequately predicted the peak heat transfer, but under
predicted the transition length.
2) Warren, Harris and Hassan (WHH) one-equation model:
This model attempts to include the effect of second mode disturbances in
addition to the first mode9. The transitional stress, incorporating both modes, is
calculated by using the following formula for the eddy viscosity length scale (
μ
l)
()
[]
t
SMTS llll
μμ
Γ+−Γ−= 1 (11.1)
Where TS
l and SM
l are the length scale contributions from the first and second
modes, respectively, which are calculated from experimental correlations. t
l
μ
is
the turbulent length scale, and is the intermittency factor. The intermittency of
the flow is defined as the fraction of the time that the flow is turbulent. The
expression for intermittency used here was developed by Dhawan and
Narasimha10 and is given by
()
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡−
−=Γ 2
2
411.0
1
λ
t
xx
e (11.2)
Where xt is the location of transition onset, is the stream wise distance
between points where = 0.25 and = 0.75. The applicability of this equation
has been confirmed for hypersonic flows. The model was used to simulate cases
in which the first mode disturbances dominate the transition process (M < 4) and
cases in which the second modes are dominant (M > 4) in Ref. 9. In all cases,
the model performed satisfactorily. This model was later implemented in GASP4
and again was found to be quite accurate.
3) Wilcox k-ω turbulence model11,12:
11-4
The low Re k-ω model developed by Wilcox was used to predict the transition
region. The prediction of the transition region was obtained by tripping the
boundary layer at a given point by decreasing the value of the dissipation so as
to destabilize the boundary layer and cause transition. The application of this
model in Ref. 4 showed that it was not very easy to trip the boundary layer at the
desired location due to the sensitivity to the initial conditions. This model
predicted a short transition length and over-predicted the peak heat transfer for
some cases.
4) Schmidt and Patankar production term modifications13:
Schmidt and Patankar have developed modifications to the production term in
the turbulent kinetic energy (TKE) equation of the Lam and Bremhorst k-ε
model14. These modifications limited the production of the kinetic energy. For the
use of this model, a trial and error method was needed to make transition occur
at the desired position by varying the inlet conditions. The results did not
compare well to the experimental cases in Ref. 4, and the method was found to
be very sensitive to the grid spacing near the wall. Due to the defects in this
model a few modifications were suggested in Ref. 4. Since it was found that the
model was difficult to trigger turbulence transition, a spot with high TKE was
introduced into the boundary layer. This spot then grew and caused transition to
take place. In order to improve the prediction of the length of the transition
region, an exponential function was used for the maximum allowable production
of TKE. These modifications improved the results for some cases but gave worse
results for other cases.
5) Algebraic transition model15,16:
The algebraic transition region model was developed at ONERA/CERT and is
described in Arnal15,16. The form of the model of Singer et al.17,18 was
implemented in Ref. 4. This model predicts transition by multiplying the eddy
viscosity by a transition function before adding it to the fluid viscosity. This
function was found to be related to the momentum thickness growth. As a result,
in test cases with severe adverse pressure gradients, where the momentum
thickness decreases, the model did not produce transition. Theoretically this
model should be compatible with any turbulence model. However, it was found
that this model did not perform well with two equation models. In Ref. 4 the
model was used with the Baldwin-Lomax model. Corrections to the calibration of
the transition function for high speed flows were also suggested in Ref. 4. The
new model predicted the cases tested better than the original model.
6) Linear combination transition model10:
This model was developed by Dey and Narasimha19 and is based on the concept
that the transition flow is a combination of the laminar and turbulent flow fields.
The contribution from laminar and turbulent values is proportioned based on the
intermittency factor developed by Dhawan and Narasimha10 mentioned above.
This model requires that a complete laminar flow simulation be run first. This
simulation is followed by a turbulent one, with the turbulent boundary layer
11-5
starting at the point of transition. The model then uses these two solutions to
generate the transitional solution. For example, the mean velocity (U) and skin
friction coefficient (Cf) are given by
()
TL UUU Γ+Γ−= 1 (11.3)
()
fTfLf CCC Γ+Γ−= 1 (11.4)
In the above equations the subscripts L and T stand for values in the laminar and
turbulent boundary layers, respectively. The peak heat transfer was not predicted
in the test cases simulated in Ref. 4. The main difficulty in getting accurate
results with this model was that one of the modeling constants needs to be
modified from case to case to obtain good results. Many researchers such as
Abid20 have used the intermittency function from the linear combination model as
an algebraic transition region function to proportion the amount of the eddy
viscosity added to the fluid viscosity. The results using this method were found to
be very similar to the linear combination model mentioned above, but there are
some noticeable differences4. The transition length was always under-predicted.
For the cases with no pressure gradient and adverse pressure gradients, the
heat transfer predicted at the end of transition and through the turbulent region
was significantly high.
In addition to the above models, efforts have been conducted to modify existing
turbulence models for turbulence transition. In Ref. 21, the performance of the
Spalart Allmaras (S-A)22 and the Baldwin-Barth (B-B)23 one-equation models and
three two-equation models for simulating hypersonic transition were evaluated.
The two equation models assessed in Ref. 21 included a low Re k-ε model with
the modifications of Nagano and Hishida24, the hybrid k-ω model of Menter
(SST)25, and the Wilcox k-ω model26. The Sandia Advanced Code for
Compressible Aerothermodynamics Research and Analysis (SACCARA) was
used to evaluate these models in Ref. 21 using two flow cases. The first case
was the flow over a flat plate at Mach 8 with flow conditions corresponding to an
altitude of 15 km, where a perfect gas was assumed. The second flow case
considered was the flow over a re-entry flight vehicle at Mach 20 and an altitude
of 24.4 km, where real gas effects need to be taken into account. The method
employed in Ref. 21 to specify transition from laminar to turbulent flow was as
follows. The turbulence transport equations were solved over the entire domain,
with a transition plane specified by the user. Upstream of this plane, the effective
viscosity was simply the laminar value whereas at downstream the effective
viscosity was the sum of the laminar and turbulent viscosities. An advantage of
this approach was that the turbulence transport equations were solved over the
whole domain, thus promoting turbulent behavior downstream of the transition
plane. On the other hand, if the turbulence source terms were simply turned on
after the transition plane, the turbulence model might not transition to turbulent
flow until farther downstream, depending on the free stream turbulence level.
However, a disadvantage of the approach was that a discontinuity in the total
![](asset-cb.png)
11-6
viscosity (laminar plus turbulent) could occur at the transition plane. All models
except the S-A model and the low Re k-ε model predicted the transition at the
correct location for the flat plate case at Mach 8. For these two models, the free
stream turbulence values needed to be increased. All the models tested in Ref.
21 provided the correct skin friction levels for this case. For the reentry flight
vehicle, the wall heat flux predicted by the S-A model, the SST model, and the
Wilcox k-ω model were found to be in reasonable agreement with the
experimental data. The B-B and the Low Re k-ε model greatly over-predicted the
heat flux in the turbulent region.
11.3 Transition Models with Onset Prediction Capability
These models not only simulate the characteristic of the transition region, but
also predict the onset of transition.
1) k-ζ turbulence/transition model:
The k-ζ turbulence model27 was used to study the effect of high disturbance
environments (HIDE) on the transitional simulations carried out in conventional
hypersonic facilities. Since HIDE cannot be described by linear stability theory, a
minimum heat flux criterion was used to determine onset of transition. This is
done by assuming initial transition onset points and employing linear interpolation
for interior points. After running a few iterations, the minimum heat flux criterion is
employed to find the locations where the wall heat flux is a minimum. The
solution is independent of the initial guess as long as the initial transition points
are ahead of the actual locations. This approach is similar to the WHH model
mentioned above. In this case, the eddy viscosity was modified using the formula
()
tnte
μ
μ
μ
Γ+Γ−= 1 (11.5)
nt
μ
is the eddy viscosity due to the non-turbulent fluctuations, and can be
calculated as
ntnt k
τ
ρ
μ
09.0= (11.6)
where
τ
nt is the non-turbulent time scale. In addition, the dissipation time scale in
the turbulent kinetic energy equation was also chosen as the combination of time
scales of turbulent and non-turbulent fluctuations. The time scale for calculating
the eddy viscosity and the dissipation time scale were derived for three different
transition mechanisms: cross flow instabilities, second mode instabilities, and
HIDE. These three mechanisms were selected because they were believed to be
responsible for transition over 3D bodies in conventional hypersonic wind
tunnels28. The simulations were performed on an elliptic cone at the Mach
number of 7.93 and were compared with experimental results. It was concluded
from the results that HIDE had a higher impact on the transition mechanism than
the other two mechanisms. The main disadvantage of this model is that it does
![](asset-cc.png)
11-7
not solve the non-turbulent fluctuations using transport equations, which limits
the flexibility of this method. Also, this model requires an initial guess for the
transition location.
2) Papp and Dash model29:
Papp and Dash proposed a concept analogous to that of the WHH model was
used. The SSGZ-J k-ε model developed by So et al.30 was implemented with
compressibility corrections for hypersonic flows. An additional transport equation
was solved for the non-turbulent fluctuations. The non-turbulent fluctuations
included the first- and second-mode mechanisms. The location of the onset of
transition31 was said to be the minimum distance along the surface for
1
1
Re ≥≡
ν
ν
μ
lt
tC, where C
μ
is the turbulent viscosity coefficient (0.09),
ν
is the fluid
kinematic viscosity, and
ν
lt is the eddy viscosity due to the non-turbulent
fluctuations which can be calculated as
ntllt kC
τ
ν
μ
= (11.6)
where
τ
nt is the viscosity time scale obtained for different transition mechanisms,
and kl is the laminar turbulent kinetic energy, which is obtained from a transport
equation in the Papp and Dash model29. This model was incorporated into a
Reynolds-Averaged Navier-Stokes (RANS) flow solver by multiplying the
turbulent eddy viscosity by the intermittency before adding to the fluid viscosity.
In all cases simulated, the transition onset was properly obtained. However, in
some cases the peak in heat transfer was not reproduced correctly. This has
been attributed to the algebraic nature of the intermittency function used. This is
the biggest disadvantage of the model.
3) Suzen and Huang model [32]:
This model uses a transport equation for the intermittency factor. This equation
not only reproduces the intermittency distribution of Dhawan and Narasimha10,
but also gives a realistic variation of the intermittency in the cross-stream
direction. The intermittency transport equation includes source terms from two
different models: the Steelant and Dick model33 and the Cho and Chung model34.
The model is incorporated into the Navier-Stokes solvers by simply multiplying
the eddy viscosity obtained from the turbulence calculations with the
intermittency factor. The Menter’s shear stress transport (SST25) model was used
to calculate the turbulent quantities. The onset of transition was determined by
comparing the local Reynolds number with a transition onset Reynolds number
(Reθ) calculated using the correlation of Huang and Xiong35, where Reθ is a
function of the free stream turbulent intensity and the acceleration parameter.
This model was tested for zero- and variable-pressure gradient flows with
different free stream turbulence intensities. The numerical result showed good
agreement with the experimental data of Savill36,37. This model is not a single
11-8
point model since it uses the free stream turbulence intensity value to calculate
the transition onset Reynolds number, which requires global parameters.
4) Walters and Leylek model38:
This model is based on the concept that bypass transition is caused by very high
amplitude stream wise fluctuations. These fluctuations are very different from
turbulent fluctuations. Mayle and Schulz39 proposed a second kinetic energy
equation to describe these fluctuations. This kinetic energy was called laminar
kinetic energy kL. In the near-wall region, the turbulent kinetic energy (TKE, kT)
was split into small-scale energy and large scale energy. The small-scale energy
(kT,s) contributes directly to the turbulence production, and the large-scale energy
(kT,l) contributes to the production of laminar kinetic energy. These two energies
can be calculated from the kT based on the turbulent length scale. The eddy
viscosities based on both scales are calculated from the respective-scale kinetic
energies. For the onset of transition, a parameter is calculated from kT, the
kinematic viscosity and the wall distance. When this parameter exceeds a certain
threshold, transition is assumed to start. The onset of transition is associated with
the reduction of kL and the consequent increase of kT (indicating the breakdown
of laminar fluctuations into turbulence). This model was incorporated into a
RANS flow solver for the calculation of the total eddy viscosity and eddy thermal
diffusivity to account for contributions from the small-scale as well as large-scale
turbulent kinetic energies. For all test cases simulated, the model responded
correctly to increases in the free stream intensity. It yielded reasonable results for
cases with high pressure gradients and streamline curvatures. Advantages of this
method are that it is very simple to implement it into the existing CFD codes
since it is based on a RANS framework. This is a single point transition model
meaning that it requires only local information, which makes this method easily
applicable to unstructured and parallel computations. The low-Re k-ε models are
typically not calibrated for transition prediction, but provide the transition location
as a by-product of their viscous sublayer formulation. Since this transition model
is developed based on the low-Re k-ε model, the embedded viscous sublayer
formulation coupled with the added transition prediction capability cannot be
calibrated independently. Hence, a change in the transition formulation would
affect the solution in the fully turbulent region. In addition, it is generally observed
that these models are not flexible enough to sufficiently cover the wide range of
transition mechanisms observed in reality40,41.
5) Local Correlation Based Transition Model (γ − Reθ model or LCTM) 40,41:
This model is based on the vorticity Reynolds number. The vorticity Reynolds
number is an extremely important parameter and is a local property and can be
easily calculated in CFD codes. The maximum value of the vorticity Reynolds
number in a boundary layer profile is directly proportional to the momentum
thickness Reynolds number. The vorticity Reynolds number is used in triggering
transition instead of directly using the momentum thickness Reynolds number.
This model solves a transport equation for intermittency and also a transport
equation for the Reynolds number based on the transition onset momentum
![](asset-ce.png)
11-9
thickness. The first transport equation includes two terms that control production.
These are Flength, a parameter which controls the length of transition zone, and
Reθc which is the momentum thickness Reynolds number at the point where the
intermittency starts to increase in the boundary layer. These two variables are
calculated from empirical functions of the transition momentum thickness
Reynolds number ( t
θ
Re ). A second transport equation is required to solve for
t
θ
Re and to include the non-local influence of the turbulence intensity, which
varies with the free stream turbulence kinetic energy and the free stream velocity.
In the case of flows with boundary layer separation, this transition model is
modified so that the intermittency is allowed to exceed the unity when the
boundary layer separates. This event results in larger production of kinetic
energy leading to correct prediction of reattachment41. This model is applied by
modifying the production and destruction terms of the original SST model using
the intermittency. The model was validated with many complicated 2D and 3D
configurations. In all cases, good agreement with the experimental data was
obtained. This model offers two main advantages: 1) it is based on local
variables; 2) it is very flexible and can be used for any mechanism as long as the
empirical correlation can be formulated. However, the empirical correlations used
with this model are proprietary.
6) The model of Lodefier et al.42:
This model is also based on the concept of pre-transitional fluctuations similar to
the Walters and Leylek model. However, this model uses the concept of
intermittency to describe the transition region. The intermittency equation used
was proposed by Steelant and Dick43. The production term of this intermittency
equation was modified by multiplying it with a new factor which is used to locate
the start of transition. This factor is zero before the start of transition and rapidly
goes to unity after the onset point. Similar to LCTM, the vorticity Reynolds
number is used in triggering transition instead of directly using the momentum
thickness Reynolds number ( t
θ
Re ). Unlike LCTM, the equation used to calculate
the critical value of t
θ
Re for transition is calculated from the local free stream
turbulence intensity and not from a transport equation. The empirical correlation
used for t
θ
Re does not include a pressure gradient term. The model is
incorporated into the SST model both by multiplying the eddy viscosity with the
intermittency and by modifying the production terms of the k and ω equations.
These modifications are used to ensure that the turbulence quantities have small
non-zero values at the start of transition as in the concept of pre-transitional
fluctuations. The main disadvantage of this model is that it uses the free stream
intensity to determine the onset of transition, which makes the model non-local
unlike LCTM.
7) Model of Lian and Shyy44:
This model was developed for simulation of flow around the wing of a micro air
vehicle (MAV). The approach used in this model was to couple an
![](asset-cf.png)
11-10
incompressible RANS solver with the en method. The k-ω model of Wilcox [11]
was selected for modeling turbulence in the RANS solver. This coupling is
accomplished as follows. The computation is started with the solution of the
RANS equations; however, the eddy viscosity is not added to the effective
viscosity. The boundary layer parameters required for the solution of the en
method are extracted from the Navier-Stokes solutions to evaluate the
amplification factor. Once the threshold value of the n-factor is reached, the flow
is allowed to become turbulent by multiplying the eddy viscosity with the
intermittency factor and adding it to the effective viscosity. The intermittency in
this case is calculated from an empirical formula. The en method employed in this
case is based on the assumptions that the initial disturbance is small and that the
boundary layer is thin.
8) Model of Arthur and Atkin3:
This method is based on linear stability theory (en method) applied within a RANS
framework. The overall process is as follows. The viscous flow over the
configuration of interest is first calculated with an initial guess of the transition
onset location. A series of pressure distributions is extracted from the RANS
solution at different “line- of-sight” positions across the span. These pressure
distributions are fed into a boundary layer code to predict the boundary layer
parameters with great accuracy and fidelity. The stability analysis, together with
some “n” factor criterion is conducted to yield the transition location. This
information is then passed onto the RANS solver for further solution. This
process is continued until the transition location and the pressure distribution are
converged. For flows with high pressure gradients, it was found that the predicted
transition location can move upstream more easily than downstream during the
iteration process. Thus, for this method it is essential that the initial guess is
downstream of the final, predicted transition location. The method does not have
any intermittency model to predict the nature of the region of transition.
11.4 LCTM Applied to a Flat Plate
The LCTM method solves two additional transport equations for the intermittency
(
γ
) and the transition momentum thickness Reynolds number ( t
θ
Re ). The
intermittency equation can be written as
()
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛+
∂
∂
+−=
jf
t
jxx
EP
Dt
D
γ
σ
μ
μ
ργ
γγ
(11.7)
The intermittency equation source term is given by
[]
()
γγρ
γ
1
5.0
11 1eonsetalength cFScFP −= (11.8)
![](asset-d0.png)
11-11
where S is the strain rate. This term is designed to be zero in the laminar
boundary layer and active everywhere the local vorticity Reynolds number
exceeds the local transition onset criteria. The vorticity Reynolds number is
defined as
Ω=
μ
ρ
ν
2
Re y (11.9)
where Ω is the vorticity magnitude. The vorticity Reynolds number is assumed to
be proportional to the momentum thickness Reynolds number
()
193.2
Remax
Re
ν
θ
= (11.10)
Fonset is defined as follows
c
onset
F
θ
ν
Re193.2
Re
1= (11.11)
(
)
(
)
2,,maxmin 4112 onsetonsetonset FFF = (11.12)
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛⎟
⎠
⎞
⎜
⎝
⎛
−= 0,
5.2
1max
3
3T
onset R
F (11.13)
(
)
0,max 32 onsetonsetonset FFF −= (11.14)
where RT is defined as (
ρ
k)/(
μω
). The two functions, Flength and Re
θ
c must still be
defined. Re
θ
c is the critical Reynolds number where the intermittency first starts to
increase in the boundary layer. This occurs upstream of the transition Reynolds
number, Re
θ
t, because the turbulence must grow to a large enough level to
trigger transition before any change in the laminar profile can be seen. Hence
Re
θ
c is the location where turbulence starts to grow and Re
θ
t is the location where
the velocity profile first departs from a laminar profile. Flength is an empirical
correlation that controls the length of the transition region. Both Re
θ
c and Flength
are defined as functions of a second transport equation variable t
θ
Re .
The destruction/relaminarzation source term is defined as
()
1
22 −Ω=
γ
γ
ρ
γ
eturba cFcE (11.15)
where
![](asset-d1.png)
11-12
4
4⎟
⎠
⎞
⎜
⎝
⎛
−
=
T
R
turb eF (11.16)
The constants for this equation are ce1=1.0, ce2=50.0, ca1=2.0, ca2=0.06, and
σ
f=1.0.
The transport equation for t
θ
Re is
()
()
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
∂
∂
+
∂
∂
+=
j
t
tt
j
t
t
xx
P
Dt
D
θ
θθ
θμμσ
ρ
ReRe (11.17)
Outside the boundary layer, the source term P
θ
t is designed to force the
transported scalar t
θ
Re to match the local value of Re
θ
t calculated from an
empirical correlation. The source term is defined as
()
()
ttttt F
t
cP
θθθθθ
ρ
−−= 0.1ReRe (11.18)
Here t is a time scale defined as
2
500
U
t
ρ
μ
= (11.19)
The blending function F
θ
t is defined as
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎝
⎛
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
−
−= ⎟
⎠
⎞
⎜
⎝
⎛
−0.1,
/10.1
/1
0.1,maxmin
2
2
2
4
e
e
y
waket c
c
eFF
γ
δ
θ
(11.20)
where
U
t
BL
ρ
μ
θ
θ
Re
= (11.21)
BLBL
θδ
2
15
= (11.22)
BL
U
y
δδ
Ω
=50 (11.23)
2
5
101
Re ⎟
⎠
⎞
⎜
⎝
⎛
−
=x
wake eF
ω
(11.24)
![](asset-d2.png)
11-13
μ
ρω
ω
2
Re y
= (11.25)
The Fwake function ensures that the wake function is not active in wake regions
downstream of an airfoil.
The system of equations is closed except for the empirical relationships for Flength,
Re
θ
c, and Re
θ
t. The first two functions are considered proprietary by the
developers. The relationship for Re
θ
t is given by
()
()()
⎪
⎩
⎪
⎨
⎧
>−
≤
⎟
⎠
⎞
⎜
⎝
⎛+−
=
−3.1,5658.05.331
3.1,
2196.0
428.58951.1173
Re
671.0
2
TuifFTu
TuifF
Tu
Tu
t
θ
θ
θ
λ
λ
(11.26)
Tu is the turbulence intensity defined as
100
2x
U
k
Tu
∞
= (11.27)
where k is the turbulent kinetic energy. F(
λ
θ
) is an empirical relationship for the
effect of pressure gradient given by
()
()
()
⎪
⎩
⎪
⎨
⎧
−+
≤−++
=
>
⎟
⎠
⎞
⎜
⎝
⎛−
−
⎟
⎠
⎞
⎜
⎝
⎛
−
0
5.0
35
5.1
32
,1275.01
0,689.40566.123986.121
5.1
θ
λ
θθθθ
θ
λ
λλλλ
λ
θ
Tu
Tu
ee
e
F (11.28)
where
λ
θ
is defined as
ds
dU
U
2
ρθ
λ
θ
= (11.29)
Note this is only valid for -0.1<
λ
θ
<0.1. For a zero pressure gradient case,
λ
θ
=0
and F(
λ
θ
)=1. The remaining two correlation functions are assumed to be
(
)
92.0Re1045.4ReRe 4+−= −ttc x
θθθ
(11.30)
3
7
Re
105.8
t
length x
F
θ
= (11.31)
The LCTM is used to predict the skin friction on a low speed flat plate for free
stream turbulence intensities of 0.3, 0.9, 3.3, and 6.5 percent in Fig. 11.1. The
![](asset-d3.png)
11-14
Model does an excellent job of predicting the onset and length of the transition
region for all but the highest free stream turbulence level. The Tu=6.5 case may
need more grid refinement near the leading edge of the plate to adequately
resolve the transition process.
Figure 11.1 Skin friction for a flat plate predicted using LCTM.
Chapter 11 References:
1. Smith, A. M. O., and Gamberoni, N., “Transition, Pressure Gradient and
Stability Theory,” Douglas Aircraft Co. Rept. ES26388, El Segundo,
California, 1956.
2. Van Ingen, J. L., “A Suggested Semi-empirical Method for the Calculation of
the Boundary Layer Transition Region,” University of Technology, Department
of Aero. Eng., Rep. UTH-74, Delft, 1956.
3. Arthur, M. T., and Atkin, C. J., “Transition Modeling for Viscous Flow
Prediction,” AIAA 2006-3052, 2006.
4. McKeel, S. A., “Numerical Simulation of the Transition Region in Hypersonic
Flow,” PhD dissertation, Blacksburg, Virginia, 1996.
5. Stainback, P., Fisher, M., and Wagner, R., “Effects of Wind-Tunnel
Disturbances on Hypersonic Boundary Transition,” AIAA Paper 72-181, 1972.
6. Holden, M. and Chadwick, K., “Studies of Laminar, Transitional, and
Turbulent Hypersonic Flows Over Curved Compression Surfaces,” Tech.
Rep. 2610-5, Calspan-UB Research Center, 1994.
7. Kimmel, R., 1993, “Experimental Transition Zone Lengths in Pressure
Gradient in Hypersonic Flow,” Symposium on Transitional and Turbulent
Compressible Flows (Kral, L. and Zang, T., eds.), Vol. FED 151, pp. 117-127,
ASME Fluids Engineering Conference, 1993.
11-15
8. Kimmel, R., “The Effect of Pressure Gradients on Transition Zone Length in
Hypersonic Boundary Layers,” Tech. Rep. WL-TR-94-3012, Wright
Laboratory, Wright Patterson AFB.
9. Warren, E. S., Harris, J. E. and Hassan, H.A., “Transition Model for High-
Speed flow,” AIAA Journal, Vol. 33, No. 8, 1995, pp. 1391-1397.
10. Dhawan, S., and Narasimha, R., “Some Properties of Boundary Layer Flow
During Transition from Laminar to Turbulent Motion,” Journal of Fluid
Mechanics, Vol. 3, No. 4, 1958, pp. 418-436.
11. Wilcox, D., “Turbulence Modeling for CFD,” La Canada, CA, DCW Industries,
1992.
12. Wilcox, D., “Turbulence and Transition Modeling for High-Speed Flows,”
NASA CR 191473, 1993.
13. Schmidt, R. and Patankar, S., “Simulating Boundary Layer Transition with
Low-Reynolds Number k-e Turbulence Models: Part 2 - An Approach to
Improving the Predictions,” Journal of Turbomachinery, Vol. 113, 1991, pp.
18-26.
14. Lam, C. and Bremhorst, K., “A Modified Form of the k-e Model for Predicting
Wall Turbulence,” ASME Journal of Fluids Engineering, Vol. 103, pp. 456-
460, 1981.
15. Arnal, D., “Three-Dimensional Boundary Layers: Laminar-Turbulent
Transition,” AGARD-FDP-VKI Special Course on Calcul Des Limites
Tridimensionelles Avec Ou Sans Decollement, 1986.
16. Arnal, D., “Laminar-Turbulent Transition Problems in Supersonic and
Hypersonic Flows,” AGARD-FDP-VKI Special Course on
Aerothermodynamics of Hypersonic Vehicles, 1988.
17. Singer, B., Dinavahi, S., and Iyer, V., “Testing of Transition-Region Models:
Test Cases and Data,” NASA CR 4371, 1991.
18. Singer, B., “Modeling the Transition Region,” NASA CR 4371, 1993.
19. Dey, J. and Narasimha, R., “An Integral Method for the Calculation of 2D
Transitional Boundary Layers,” Journal of Aircraft, Vol.27, no.10, 1990,
pp.859-865.
20. Abid, R., “A Study of Turbulence Models for Prediction of Transitional
Boundary Layers,” Instability, Transition, and Turbulence (Hussaini, Kumar,
and Streett, eds.), Springer-Verlag, 1992.
21. Roy, C. J., and Blottner, F. G., “Assessment of One- and Two- Equation
Turbulence Models for Hypersonic Transitional Flows,” Journal of Spacecrafts
and Rockets, Vol. 38, No.5, 2001, pp. 313-325.
22. Spalart, P. R., and Allmaras, S. R., “A One-Equation Turbulence Model for
Aerodynamic Flows,” AIAA 92-0439, 1992.
23. Baldwin, B. S., and Barth, T. J., “A One-Equation Transport Model for High
Reynolds Number Wall-Bounded Flows,” NASA TM-102847, 1990.
24. Nagano, Y., and Hishida, M., “Improved Form of the k-ε Model for Wall
Turbulent Shear Flows,” J. of Fluid Eng., 109(2), 1987, pp. 156-160.
25. Menter, F. R., “Two-Equation Eddy Viscosity Turbulence Models for
Engineering Applications,” AIAA J., Vol. 32, No. 8, 1994, pp. 1598-1605.
11-16
26. Wilcox, D. C, Turbulence Modeling for CFD, 2nd ed., La Canada, CA, DCW
Industries, 1998, Chap.5, pp. 119-122 and 227-271.
27. Xiao, X., Edwards, J. R. and Hassan, H. A., “Transitional flow over an Elliptic
Cone at Mach 8,” Journal of Spacecraft and Rockets, Vol. 38, No. 6, 2001.
28. Mokovin, M. V., “Bypass Transition to Turbulence and Research Desiderata,”
Transition in Turbines, edited by R. W. Graham, NASA CP2386, 1885, pp.
161-204.
29. Papp, J. L. and Dash, S. M., “Rapid Engineering Approach to Modeling
Hypersonic Laminar-to-Turbulent Transitional Flows,” Journal of Spacecraft
and Rockets, Vol. 42, No. 3, 2005, pp. 467-475.
30. So, R. M. C., Sarkar, S., Gerodimos, G., and Zhang, J., “A Dissipation Rate
Equation for Low-Reynolds-Number and Near-Wall Turbulence,” Theoretical
Computational Fluid Dynamics, Vol. 9, 1997, pp. 47-63.
31. Warren, E. W. and Hassan, H. A., “Alternative to the en Method for
Determining Onset of Transition,” AIAA Journal, Vol. 36, No. 1, 1998, pp. 111-
113.
32. Suzen, Y. B. and Huang, P. G., “An Intermittency Transport Equation for
Modeling Flow Transition,” AIAA 2000-0287, 2000.
33. Steelant, J. and Dick, E., “Modeling of Bypass Transition with Conditioned
Navier-Stokes Equations Coupled to an Intermittency Transport Equation,”
International Journal for Numerical Methods in Fluids, Vol. 23, 1996, pp. 193-
220.
34. Cho, J. R., and Chung, M. K., “A k-e-v Equation Turbulence Model,” Journal
of Fluid Mechanics, Vol. 237, 1992, pp. 301-322.
35. Huang, P. G., and Xiong, G., “Transition and Turbulence Modeling of Low
Pressure Turbine Flows,” AIAA 98-0339, Reno, NV, 1998.
36. Savill, A. M., “Some Recent Progress in the Turbulence Modeling of By-pass
Transition,” Near-Wall Turbulent Flows, edited by R. M. C. So, C. G. Speziale
and B. E. Launder, Elsevier Science Publishers B.V., 1993, pp. 829-848.
37. Savill, A. M., “Further Progress in the Turbulence Modeling of By-pass
Transition,” Engineering Turbulence Modeling and Experiments 2, edited by
W. Rodi and F. Martelli, Elsevier Science Publishers B.V., 1993, pp. 583-592.
38. Walters, D. K. and Leylek, J. H., “A New Model for Boundary Layer Transition
Using a Single-Point RANS Approach,” Journal of Turbomachinery, Vol. 126,
Issue 1, 2004, pp. 193-202.
39. Mayle, R. E., and Schulz, A., ‘‘The Path to Predicting Bypass Transition,’’
Journal of Turbomachinery, Vol. 119, 1997, pp. 405–411.
40. Langtry, R. B. and Menter, F. R., “Transition Modeling for General CFD
Applications in Aeronautics,” AIAA 2005-522, 2005.
41. Menter, F. R., Langtry, R. B., Likki, S. R., Suzen, Y. B., Huang, P. G., Volker,
S., “A Correlation-based Transition Model Using Local Variables Part I -
Model Formulation,” Journal of Turbomachinery, Vol. 128, 2006, pp. 413-422.
42. Lodefier, K., Merci, B., De Langhe, C., and Dick, E., “Intermittency Based
RANS Bypass Transition Modeling,” Progress in Computational Fluid
Dynamics, Vol. 6, Nos. 1/2/3, 2006.
11-17
43. Steelant, J., and Dick, E., “Modeling of Laminar- Turbulent Transition for High
Freestream Turbulence,” Journal of Fluids Engineering, Vol. 123, No. 1, 2001,
pp. 22-30.
44. Lian, Y., and Shyy, W., “Laminar-Turbulent Transition of a Low Reynolds
Number Rigid or Flexible Airfoil,” AIAA Journal, Vol. 45, No .7, 2007, pp.
1501-1513.