Lectures For ES912, Term 1, 2003
Lectures For ES912, Term 1, 2003
Lectures For ES912, Term 1, 2003
P +
2
u
u = 0 (1)
where the density is assumed to be constant and is the kinetic viscosity due to molecular motion.
The Reynolds number:
R = UL/ (2)
is the ratio of nonlinear term to viscous term (assuming no turbulence) where U is large-scale velocity,
L is a length.
Coherent Structures
Flow behind a cylinder
Laminar transitions
Turbulence: drag independent of R
von Karman vortex street
Flow behind any object
Boundary layer: at plate
Mixing layer
Jets
4
1.5 Properties
Irregularity
All turbulent ows are irregular or contain elements of randomness. It is not purely random (like
an ideal gas), nor is it deterministic.
Diusivity
Fine ow patterns yield strong diusivity
which causes rapid mixing and increased rates of momentum, heat, and mass transport.
GOAL: predicting these transport rates
Large Reynolds Numbers
All of the above assumes high Reynolds numbers, assumptions of randomness and statistical methods
are less useful if R is small.
Three-dimensional vorticity dynamics
Turbulence is inherently 3D and rotational.
All velocity elds can be derived from the vorticity =
u
Vorticity is produced by vortex stretching.
3D vorticity is concentrated in vortex tubes.
Turbulence is the motion of chaotic vortex tubes.
5
Dissipation
Turbulent ows are always dissipative
Viscous shear stresses perform deformation work that converts kinetic energy into heat (internal
energy/temperature)
Turbulent kinetic energy will decay unless there is a continuous supply of energy to sustain it.
(Wave motion might not be turbulent: waves are essentially non-dissipative)
Continuum
Turbulence is described by continuum equations, not particles.
Turbulence is a uid phenomena
Turbulence is a feature of uids: liquids, or gases that are incompressible.
6
1.6 Reynolds averaging
Divide the total velocity u into a mean, or large-scale component, U and a uctuating component u:
u
i
= U
i
+ u
i
(3)
For a nearly steady ow U
i
might be a time average. In large-eddy simulation U
i
is the velocity smoothed
onto a coarser mesh.
DIAGRAM
Properties: average of derivative is derivative of average:
u
i
x
j
=
U
i
x
j
u
i
u
j
= (U
i
+ u
i
)(U
j
+ u
j
) = U
i
U
j
+ u
i
u
j
+ U
i
u
j
+ U
j
u
i
= U
i
U
j
+ u
i
u
j
Correlation coecients: Consider c
ij
= +u
i
u
j
/(u
2
i
u
2
j
)
1/2
Uncorrelated: c
ij
= 0
Perfectly correlated c
ij
= 1
In general in turbulence 0 < |c
ij
| < 1
Reynolds stress:
ij
= u
i
u
j
.
(u
)u
i
= u
j
j
u
i
=
j
(u
i
u
j
) =
j
ij
Neglect time deriviatives to get the Reynolds averaged equations:
(
U
)U
i
=
1
i
P +
2
U
i
ij
(4)
Turbulence modelling is to a large degree the art of modelling
ij
.
LES: Large-eddy simulation: Retain time derivative:
t
U
i
, model
ij
7
2 Classical phenomenological equations
2.1 Mixing length models
Use particle wandering analogy to describe a mixing length model
For xed velocity uctuations and mixing length, equivalence to an eddy viscosity model
Derive log law using mixing length
Application to passive scalar diusion
If more than one velocity or length scale, mixing length models fail
Mixing length models were originally proposed by Prandtl to explain boundary layer scaling. The
objective is to predict the dispersion of a parcel due to a background shear.
Consider a parcel that has been displaced by x
2
in in the y direction in a shear U
1
/y. The
momentum decit will be
M = x
2
U
1
y
The rate at which this momentum decit will be created is obtained by multiplying by u
2
= dx
2
/dt
giving a transport or stress of
12
=
dx
2
dt
x
2
U
1
y
=
1
2
dx
2
2
dt
U
1
y
(5)
so that
12
represents the rate of change of the average displacement of the parcel in y. In the kinetic
theory for an ideal gas, this is determined by the mean free path of a particle and the thermal velocity,
u
, both microscopic quantities determined by the density and temperature. In an analogy to this, in
mixing length theory we will dene a mixing length and a uctuating velocity u
= c
2
U
s
where c
1
and c
2
are ow
8
dependent and must be determined experimentally. Hopefully, we will nd that they are approximately
independent of Reynolds number.
For example in a mixing layer we will use as L the width of layer at a location downstream and U
s
will be the velocity dierence across the layer U
s
= U
2
U
1
.
With these assumptions, (5) becomes
12
= u
U
1
y
= c
1
c
2
LU
s
U
1
y
(6)
Mixing length models depend upon there being only one large-scale length scale and one large-scale
velocity scale in the problem. If there are two or more length or velocity scales, we should not expect
this approximation to be valid. For example, in thermal convection where in addition to the height
above a surface you might need to consider the thickness of the thermal plumes. Where mixing length
theory can give useful relations is in a mechanically driven boundary layer over a smooth surface, shown
below.
2.2 Derivation of a simple eddy viscosity
While the concept of an eddy viscosity goes back to Boussinesq (1877) on purely dimensional grounds,
we will derive the eddy viscosity starting with mixing length theory.
What if boundaries are not an important consideration, for example far enough downstream in a
wake, in a mixing layer, or a jet? Then across the layer, L and U
s
can be assumed to be constant and
we can rewrite (6) as
12
= (c
1
c
2
LU
s
)
U
1
y
=
T
U
1
y
and
ij
=
T
U
i
x
j
(7)
where
T
is the eddy viscosity. For all of the simple cases above, we will be able to derive realistic
velocity proles by assuming that
T
is constant across the layer and dependent upon the large-scale
9
parameters in the manner given, that is
T
= (c
1
c
2
)LU
s
. The underlying principle is that the random
motion of the turbulent eddies mixes momentum or heat in a manner similar to manner of an ideal gas.
This approximation works for an ideal gas because there is a separation in scale between the large-scale
motion and the molecular motion. Applying this approximation to turbulence is limited to regimes
where the turbulence can be considered homogeneous (distribution of uctuations is uniform in space).
The constant c
1
c
2
must still be determined experimentally.
Given
ij
from (7) Show that it has the same form as the Newtonian viscosity.
In large-eddy simulation,
T
is assumed to depend upon some local property of the ow and the
mesh size x . For example,
T
= C
s
(x)
2
|e|, where C
s
is a constant that ultimately must be
determined experimentally (or using direct numerical simulations), but can be estimated using the
statistical theory of turbulence. |e| is an inverse timescale, usually taken to be the modulus of local
strain rate e
ij
= 0.5(u
i,j
+u
j,i
) at the grid scale x. One of the primary objectives of several turbulence
research groups is to determine the best form for the eddy viscosity in terms of these ow conditions
and hopefully will be discussed later in the course.
2.3 Derivation of log-law of a boundary layer using mixing length theory
In the case of a a simple sheared, at plate boundary layer, one can derive a useful relation with mixing
length theory by assuming that will depend solely upon the distance from the wall, so that if L = y,
then = c
1
y. To do this we must start with the Reynolds equation (4), the only term of which is
U
2
U
1
/x
2
giving
U
2
U
1
x
2
=
1
x
2
12
which can be reduced futher because U
2
= 0 giving
12
/x
2
= 0 or
12
=constant.
Therefore, outside the viscous sub-layer, which is the source of the momentum decit, there must be
a constant ux transport, or a constant Reynolds shear stress. This is a constant stress layer. Assume
10
at the wall that
12
(0) = u
2
where u
12
(x
2
) = u
U
1
x
2
= x
2
u
U
1
x
2
where it has been assume that = c
1
L = x
2
and u
= c
2
U
s
= c
2
u
and = c
1
c
2
. Equate this to
12
(0) = u
2
x
2
=
u
x
2
which integrates to
U
1
u
=
1
ln x
2
+ constant (8)
where = 0.4 is the Karman constant. The constant term added is determined by the viscous sublayer.
The eect of the viscous sublayer is ow dependent and depends upon how the boundary layer becomes
unstable and how smooth or rough the walls are. But the
1
0.4
ln x
2
is very robust, that is it is independent
of the ow details.
11
3 Similarity solutions
Start with the Reynolds averaged equations again (4). This time in addition to U
1
= 0 we will assume
that U
2
= 0. This is because the proles are spreading and because U
1
is changing with x
1
, so U
2
= 0
is required by incompressibility:
U = 0.
Assumptions:
U = (U
1
, U
2
, 0)
Wakes : U
1
(x
2
) U
0
U
s
(9)
Jets and mixing layers : U
1
U
s
(10)
/x
1
<< /x
2
or << L
Due to U
1
/x
1
+ U
2
/x
2
= 0
U
2
O(U
s
/L)
Assume a small u (magnitude relative to U
s
to be determined) such that:
uv = u
1
u
2
= O(u
2
), u
2
1
= O(u
2
) u
2
2
= O(u
2
)
3.1 Momentum equations for
U = (U
1
, U
2
, 0)
U
1
U
1
x
1
+ U
2
U
1
x
2
+
x
2
(uv) +
x
1
(u
2
1
) =
1
P
x
1
+
_
_
_
2
U
1
x
2
1
+
2
U
1
x
2
2
_
_
_ (11)
U
1
U
2
x
1
+ U
2
U
2
x
2
+
x
1
(uv) +
x
2
(u
2
2
) =
1
P
x
2
+
_
_
_
2
U
2
x
2
1
+
2
U
2
x
2
2
_
_
_ (12)
12
Cross-stream momentum equation (for U
2
)
U
1
U
2
x
1
+ U
2
U
2
x
2
+
x
1
(uv) +
x
2
(u
2
2
) =
1
P
x
2
+
_
_
_
2
U
2
x
2
1
+
2
U
2
x
2
2
_
_
_
Now we identify the order of the terms. U
1
for the advection of the rst term will be retained and
labeled U
1A
, but express everything else in scaling variables, large scales: U
s
, L, and small scales: u, .
U
1A
U
2
x
1
U
1A
L
U
s
L
U
2
U
2
x
2
_
_
U
s
L
_
_
2
1
x
1
(uv)
u
2
L
x
2
(u
2
2
)
u
2
P
x
2
=?
2
U
2
x
2
1
L
2
U
s
2
U
2
x
2
2
2
U
s
L
13
Write all the terms as a factor times u
2
/ to estimate which terms are signicant.
U
1A
U
2
x
1
U
1A
L
U
s
L
=
_
_
U
1A
u
U
s
u
_
_
L
_
_
2
_
_
u
2
(13)
U
2
U
2
x
2
_
_
U
s
L
_
_
2
1
=
_
_
_
_
U
s
u
_
_
2
_
_
L
_
_
2
_
_
u
2
x
1
(uv)
u
2
L
=
_
_
L
_
_
u
2
x
2
(u
2
2
)
u
2
= [1]
u
2
P
x
2
=?
2
U
2
x
2
1
L
2
U
s
L
=
_
_
U
s
u
1
R
_
_
L
_
_
3
_
_
u
2
2
U
2
x
2
2
2
U
s
L
=
_
_
U
s
u
1
R
_
_
L
_
_
_
_
u
2
where R
= U
s
/. All the terms are multiplied by u
2
/ on the right. Now, compare the terms
in brackets. We expect that u will be the order of U
s
, so U
s
/u will not be very large and will not
compensate /L, so that all but two of the terms have factors of /L multiplying them. They are
(u
2
2
)/x
2
and P/x
2
, which if u
2
= v and x
2
= y gives (v
2
)/y =
1
P/x
2
. Therefore the
mean pressure will be
P + v
2
= P
(14)
Which is just the Bernoulli principle.
14
Streamwise momentum equation (for U
1
)
U
1A
U
1
x
1
+ U
2
U
1
x
2
+
x
2
(uv) +
x
1
(u
2
1
+
1
P) =
_
_
_
2
U
1
x
2
1
+
2
U
1
x
2
2
_
_
_
Now we identify the order of the terms. Retain U
1A
, but express everything else in scaling variables.
U
1A
U
1
x
1
U
1A
U
s
L
U
2
U
1
x
2
_
_
U
s
L
_
_
U
s
x
2
(uv)
u
2
replacing
1
P
x
1
by
v
2
x
1
(14)
x
1
(u
2
1
v
2
)
u
2
L
2
U
1
x
2
1
U
s
L
2
2
U
1
x
2
2
U
s
2
15
Write all the terms as a factor times u
2
/ to estimate which terms are signicant.
U
1A
U
1
x
1
U
1A
U
s
L
=
_
_
U
1A
u
U
s
u
_
_
L
_
_
_
_
u
2
(15)
U
2
U
1
x
2
_
_
U
s
L
_
_
U
s
=
_
_
_
_
U
s
u
_
_
2
_
_
L
_
_
_
_
u
2
x
2
(uv)
u
2
= [1]
u
2
x
1
(u
2
1
v
2
)
u
2
L
=
_
_
L
_
_
u
2
2
U
1
x
2
1
U
s
L
2
=
_
_
U
s
u
1
R
_
_
L
_
_
2
_
_
u
2
2
U
1
x
2
2
U
s
2
=
_
_
U
s
u
1
R
_
_
u
2
L
_
_
O(1) and maybe
_
_
U
s
u
_
_
2
_
_
L
_
_
O(1) (16)
Therefore, the equation we will solve is
U
1A
U
1
x
1
+
x
2
(uv) = 0 (17)
Separate simplications for U
1A
will now be made for wakes and for jets and mixing layers.
16
The question now is how to set U
1A
, U
s
and u. We will assume in all cases that U
s
u.
Wakes: From (9) U
1A
U
0
U
s
. If U
s
u, then by (16),
U
s
U
0
u
U
1A
O(
L
) which is small.
Therefore we can replace U
1A
by U
0
and (17) becomes
U
0
U
1
x
1
+
x
2
(uv) = 0 (18)
Note that in this case U
2
U
s
(/L) U
0
(/L)
2
and is very small.
Jets and mixing layers: Assume U
1A
U
s
. Then the rst two terms in (15) are the same order of
magnitude and (17) is replaced by
U
1
U
1
x
1
+ U
2
U
1
x
2
+
x
2
(uv) = 0 (19)
In this case U
2
U
s
(/L) U
1
(/L) so that U
2
(/x
2
) O(U
1
(/L)/) O(U
1
/L), the same
order as the rst term.
The objective now is to assume self-preservation forms for U
1
and
12
U
1
(y) = U
s
(x
1
)f(y/) and
12
= uv = U
2
s
g(y/) (or u
2
g(y/) ), (20)
make an additional eddy viscosity assumption to relate g and f, and use (18) to predict the f and g
dependence upon y/ and the x
1
dependence of U
s
and . And use (19) for jets and mixing layers. Note
that for (19) incompressibility must be used to predict U
2
, which leads to integral equations.
17
3.2 Wakes: self-preservation
Assume (U
0
U
1
) = U
s
f(y/) and uv = U
2
s
g(y/)
giving
U
x
=
dU
s
dx
f +
U
s
d
dx
f
where = y/ and
refers to derivative with respect to y/.
uv
y
=
U
2
s
d
dx
f
)
U
2
s
= 0
U
0
_
U
2
s
dU
s
dx
_
_f + U
0
_
_
1
U
s
d
dx
_
_
f
= g
=
1
R
T
f
(21)
where the red shows the terms that must be independent of x, the downstream distance. This is possible
if
x
n
and U
s
x
n1
To x n we must introduce two new concepts
The momentum integral
The momentum thickness
and then make an eddy viscosity assumption to determine g
in terms of f
.
18
3.2.1 Momentum integral
Rewrite (19) using incompressibility as (adding derivatives of U
0
, which are 0):
x
[U
1
(U
1
U
0
)] +
y
[U
2
(U
1
U
0
)] +
y
uv = 0
Consider the integral over y of the 1st term, which for the 2nd and 3rd terms yields the values inside
the y-derivatives at large y. Then, because U
2
and uv vanish for large y,
d
dx
_
U
1
(U
1
U
0
)dy = 0
that is the integral is constant or
U
1
(U
1
U
0
)dy = M (22)
To set n for U
s
and for a wake, take the largest term U
1
U
0
and substitute using f from (20) to
get
U
0
U
s
T
U
y
=
T
1
U
s
=
1
R
T
f
(26)
and dene the turbulent Reynolds number as R
T
= U
s
/
T
Now insert (24) U
s
= Ax
1/2
and
= Bx
1/2
and (26) into (21)
U
0
_
U
2
s
dU
s
dx
_
_f + U
0
_
_
1
U
s
d
dx
_
_
f
=
1
R
T
f
to get (f
+ f) + f
= 0 (27)
using
U
2
s
dU
s
dx
=
1
U
s
d
dx
=
1
2
B
A
and dening A, B such that
1
2
R
T
U
0
B
A
= 1 or
A
B
=
R
T
U
0
2
The solution of (27) is
f = exp(
1
2
2
)
20
Summary
U
0
U
1
x
1
+
x
2
(uv) = 0
Use
U
x
=
dU
s
dx
f +
U
s
d
dx
f
Use
uv
y
=
U
2
s
=
T
U
y
=
U
2
s
T
U
s
=
U
2
s
1
R
T
f
U
0
(
dU
s
dx
f +
U
s
d
dx
f
)
U
2
s
= 0
U
0
_
U
2
s
dU
s
dx
_
_f + U
0
_
_
1
U
s
d
dx
_
_
f
= g
=
1
R
T
f
U
0
_
_
B
2A
_
_
f + U
0
_
_
B
2A
_
_
f
=
1
R
T
f
U
s
= const U
s
= Ax
1/2
= Bx
1/2
(f
+ f) + f
= 0
Dening A, B such that
1
2
R
T
U
0
B
A
= 1
The solution is:
f() = exp(
1
2
2
)
21
Find C
U
and C
if
U
s
/U
0
= C
U
(x/)
1/2
/ = C
(x/)
1/2
using R
T
= 12.5 for a plane wake and (22) U
0
U
s
f()d = U
2
0
. Using
_
exp(
1
2
2
)d =
(2)
1/2
we get the condition
=
U
s
U
0
f()d = (2)
1/2
U
s
U
0
= (2)
1/2
AB
U
0
AB = U
0
(2)
1/2
Combining with
A
B
=
R
T
U
0
2
we get A
2
= U
2
0
R
T
()
2(2)
1/2
, A = U
0
()
1/2
_
_
_
12.5
2(2)
1/2
_
_
_
1/2
= 1.58U
0
()
1/2
The large-scale Reynolds number for a plane wake will go as R
= U
s
/ = C
U
C
/U
0
, which will
be constant with x
1
.
Example problem: Find C
= U
s
/ = 0.4U
0
/
Therefore, turbulence will persist arbitrarily far downstream in a plane wake.
This will not be the case for an axisymmetric wake. Experiments show that for an axisymmetric wake
R
T
= 14.1. In that case we still need:
U
2
s
dU
s
dx
=
1
U
s
d
dx
But the momentum integral is cylindrical:
U
0
U
s
2
_
f()(2)d = U
2
0
2
=constant, how will U
s
and depend on x
1
?
When plotted against experimental data, the predicted prole is in very good agreement with the
experimental prole in the middle of the wake. Near the edges the theoretical prediction is somewhat
higher than the experiments. It is believed that this is due to intermittency. Intermittency is the
22
property that vortical motion in turbulence is not uniformly distributed, but tends to be patchy. Near
the centre of the wake the turbulence is close to being homogeneous. But near the edges where there
is constant entrainment of non-turbulent uid, so the degree of intermittency will be higher and the
eective eddy viscosity that should be used should be smaller.
3.3 Jets and mixing layers: self-preservation
In a wake the width gets narrower with respect to x
1
By (16)
u
U
s
O(
L
)
1/2
But in turbulence u U
s
, so /L is constant. That is
Bx
The nasty equation that determines f in this case is
U
s
dU
s
dx
f
2
+ ... + (...)
_
0
fd + ... = g
It is
U
s
dU
s
dx
not U
2
s
in the denominator because U
0
= U
s
in U
0
U
2
s
dU
s
dx
Taking a look at the rst term, we want
U
s
dU
s
dx
= constant
which can be satised by any power law for U
s
= Bx
m
, including m = 0. This is the case for mixing
layers. For a jet the momentum integral is
_
U
2
1
dy = U
2
s
f
2
()d = U
2
j
d
Giving U
s
x
1/2
and R
=
U
s
x
1/2
23
3.4 Energy budget
Let us consider the equation for the large-scale energy in a plane wake, E
L
= (1/2)U
2
1
= (U
0
+ U
1
)
2
.
This can be obtained by multiplying (18) by U
1
= U
0
+ U
1
. The resulting equation is:
U
0
x
1
_
_
1
2
U
2
0
+ U
0
U
1
+
1
2
(U
1
)
2
_
_
+ (U
0
+ U
1
)
x
2
(uv) = 0
The derivative of (1/2)U
2
0
is clearly 0. U
0
(/x
1
)U
1
+ U
0
(/x
2
)(uv) is U
0
multiplying (18), so that
drops out, and we are left with
U
0
x
1
1
2
(U
1
)
2
+ U
x
2
(uv) = 0
or U
0
x
1
1
2
(U
1
)
2
+
x
2
(uvU
1
) = uv
U
1
x
2
(28)
The total energy E
s
=
_
1
2
(U
1
)
2
dy U
2
s
(x
1/2
)
2
x
1/2
x
1/2
which is decaying.
Why does the energy decrease?
The rst term represents the advection of the large-scale energy decit. The term with the Reynolds
stress uv has been split to separate the conservative transport and the drag. The transport is conser-
vative because it integrates to 0. Because uv is negative (for y > 0), the drag term
uv
U
1
x
2
< 0 removes energy from the system.
Where does the energy go?
24
3.5 Turbulent energy budget
The energy budget equation from the rst week of lectures is,
t
(
1
2
u
2
) +
(u
1
2
u
2
) =
1
(uP) +
j
(u
i
j
u
i
) +
ij
i
u
j
(
j
u
i
)
2
(29)
where
ij
= u
i
u
j
is the Reynolds stress and
u = 0 has been used to create the divergences in the
transport due to the nonlinear, pressure, and viscous terms, red is used to highlight the Reynolds stress
terms, which result in the production of turbulent kinetic energy and the viscous or dissipation terms.
If the turbulence at a location is steady (
t
= 0) these will balance.
We now want to consider only those terms appropriate for a wake, using order of magnitudes estimates
as was done to get the prole equations.
Figure 1: Budgets for energy of mean velocity from (28) and for the uctuating turbulent velocity from (30).
Rewrite the turbulent kinetic energy as
1
2
q
2
=
1
2
(u
2
1
+ u
2
2
+ u
2
3
) =
1
2
(u
2
+ v
2
+ w
2
)
25
Its budget for a wake is
0 = U
0
x
1
1
2
q
2
uv
U
1
x
2
x
2
v
_
_
1
2
q
2
+
P
_
_
(30)
The terms are advection, production, transport, and dissipation. Note that the production here is the
opposite of the drag in (28). So globally, excluding the dissipation in (30), energy is conserved.
The dierences with the full budget equation (29) are that (30) uses
1
2
q
2
for
1
2
u
2
, the viscous transport
term is neglected, the time derivative is neglected, and the advection due to the driving velocity U
0
has
been pulled out of
(u
1
2
q
2
). Except for the large scale advection, all derivatives with respect to x
1
are
neglected because /
x1
1/L.
26
Figure 2: Budget for energy of mean velocity from (28).
U
0
x
1
1
2
(U
1
)
2
+
x
2
(uvU
1
) = uv
U
1
x
2
27
Figure 3: Budget for the uctuating turbulent velocity from (30).
0 = U
0
x
1
1
2
q
2
uv
U
1
x
2
x
2
v
_
_
1
2
q
2
+
P
_
_
28
4 Statistical theory of turbulence
4.1 Moments
When Reynolds averaging was introduced we divided the total velocity u into its mean, or large-scale
component, U and a uctuating component u = u
2
= u
2
=
_
dV u
2
B(u)
where is the standard deviation.
Even though the 1st moment is 0, the 3rd moment might not be. The 3rd moment represents any
asymmetry in the distribution. Its non-dimensionalised measure is called the
skewness S = u
3
/
3
Its sign can tell you something about the direction of transport because it is the 3rd moments that enter
into the kinetic energy transport equation (29).
The non-dimensionalised 4th moment is the
kurtosis or atness K = u
4
/
4
29
Closure problem
The equations for the 2nd moments, that is the equations for the energies and Reynolds stresses
u
i
u
j
, depend upon the 3rd moments, the transport terms.
The equations for the 3rd moments depend upon the 4th moments
etc.
That it is impossible to close the higher-order moment equations in terms of combinations of the lower-
order moments is the closure problem. Many turbulence models depend upon assumptions about
replacing these higher-order moments with products of lower order moments. For example replacing
u
4
by A(u
2
)
2
This is the quasi-normal approximation and is infamous for making the equations non-realisable. That
is, it will predict negative energies. All practical closure schemes, such as k , have strategies for
addressing this problem in simple ways.
30
4.2 Correlations and structure functions
The two-point velocity correlation R
ij
and the second-order structure function are dened as follows:
R
ij
(r)
u
i
(x)u
j
(x + r)
(u
2
i
u
2
j
)
1/2
S
ij
(r) =
(u
i
(x + r) u
i
(x))(u
j
(x + r) u
j
(x))
(u
2
i
u
2
j
)
1/2
(31)
Isotropy implies that any component R
ij
of the 2-point velocity correlation tensor can be written in
terms of f(r) and g(r) as
R
ij
= (f(r) g(r))
r
i
r
j
r
2
+ g(r)
ij
Moreover the continuity equation gives
u
i
x
i
= 0 =
R
ij
r
j
= 0 or f +
1
2
r
df
dr
= g (32)
So there is only one independent function f(r) required to completely dene R
ij
. Substitution of (3.11)
into (3.10) shows that
In isotropic ows, i = j, R
ij
(r) = S
ij
(r) = 0
And taking the limit as r 0, the Reynolds stress
ij
= uv = R
ij
(0) = 0
If i = j and both are in the x direction, these are longitudinal correlations and structure functions.
If i = j and both are perpendicular to the x direction, these are transverse correlations and structure
functions.
By (32), the transverse correlations and structure functions are related. An important test for
isotropy is if R
11
(r) and R
22
(r), or S
11
(r) and S
22
(r), are related in this way, where r is in the 1
direction.
31
There are three length scales that matter:
Integral or forcing scale =
0
R
11
(r)dr
Expanding f(r) for small r : f(r) = f(0)
. .
=1
+r f
(0)
. .
=0
+(r
2
/2)f
(0) + ... = 1
r
2
/2
2
where is the Taylor microscale
=
_
u
2
1
(u
1
/r)
2
is a measure of the ne scales because it is based on a representation of the velocity correlation
valid for small separations. For isotropic turbulence, the dissipation rate
= 15(
u
1
r
)
2
= 15
u
2
1
2
or
2
=
u
2
1
/15
(33)
Kolmogorov microscale : The only quantities that matter are:
the dissipation rate
the viscosity .
The only length scale that can be formed from and is
=
_
_
_
_
_
_
1/4
(34)
32
4.3 Measurement
The measurement of correlations, structure functions, and the length scales above is usually by
A single hot-wire probe measuring a single component of velocity at a single point.
What is measured are time signals.
Taylor frozen turbulence assumption:
If the mean ow U
0
>> u, then the assumption is made that
u(r) = u(t) where r = U
0
t
4.4 Fourier transforms
To aid us in developing closure models and as a tool for understanding the dynamics of turbulence, much
eort has been spent on investigating the properties of the Fourier transformed velocity, the Fourier
transformed equations and the kinetic energy spectrum.
A 1D Fourier transform for a periodic domain 0 x 2 where u(2) = u(0) is a real number is
dened as
u(k) =
1
2
_
2
0
u(x)e
ikx
dx where u(x) =
_
u(k)e
ikx
dk
is the forwards Fourier transform equation. The rst equation is the inverse transform equation. Prac-
tically, the integrals are replace by a nite sums of n periodic points x
i
or wavenumbers k
j
.
u(k
j
) =
1
2
i=0,n1
u(x
i
)e
ik
j
x
i
u(x
i
) =
j=n/2+1,n/2
u(k
j
)e
ik
j
x
i
where x
i
= 2i/n and k
j
= j .
u(k) is a complex number where in 1D, u(k) = u
u(q)e
iqx
dq)e
ikx
dx =
1
2
_
2
0
(q k)u(q)dx = u(k)
where
_
e
i(qk)x
dq = (q k) is the Kronecker function
Extensions
By normalising x we can apply this to any periodic domain.
What has been dened above is the real-half complex FFT
To extend to 2 and 3D we need to add full FFTs that transform from complex to complex numbers
4.5 Fourier transformed Navier-Stokes equations
u
t
+ (u
)u =
1
P +
2
u
Solution of P: Apply incompressibility to above equations
{...}
1
2
P = (
i
u
j
)(
j
u
i
) + boundary conditions
In general this is dicult to solve. Except for idealised boundary conditions where we can Fourier
transform the Navier-Stokes equations
u(
k)
t
+
_
d
3
pi p u(
k p)u( p) =
1
(i
kP) k
2
u(
k) (35)
where P =
k
j
k
2
_
d
3
p p u(
k p)u
j
( p).
34
4.6 Solving P in Fourier space
The diculty of solving the incompressible equations in physical space is the pressure term. In physical
space, by applying incompressibility to Navier-Stokes
_
_
_
u
t
+ (u
)u =
1
P +
2
u
_
_
_
where the terms in red drop out due to incompressibility and we are left with
(
i
u
j
)(
j
u
i
) =
1
2
P (36)
In Fourier space (
i
u
j
)(
j
u
i
) =
1
2
P becomes( using p
j
u
j
( p) = 0)
_
d
3
p (k
j
p
j
)p
i
u
i
(
k p)u
j
( p) = k
j
_
d
3
p p
i
u
i
(
k p)u
j
( p) = k
2
P(
k) giving
P =
k
j
k
2
_
d
3
p p u(
k p)u
j
( p). Then to get
P multiply P by ik
i
to get
ik
i
P = i
k
i
k
j
k
2
_
d
3
p p u(
k p)u
j
( p)
and the pressure can be folded into the non-linear term converting (35) into
u
i
(
k)
t
+
_
_
ij
k
i
k
j
k
2
_
_
_
d
3
p i p u(
k p)u
j
( p) = k
2
u
i
(
k) (37)
35
4.7 Fourier energy equation
To get an equation for the energy spectrum E(k) as a function of a 3D wavenumber k = |k| rst
multiply by u, replace the cubic non-linear term that includes the convolution
N with the projection
factor by T
i
(
k),
let the energy E(
k) =
1
2
u
2
(
k) = 2k
2
E(
k) to get
u
_
_
u(
k)
t
+
_
_
ij
k
i
k
j
k
2
_
_
N(
k) = k
2
u(
k)
_
_
which becomes
E(
k)
t
+ T(
k) = (
k)
Finally, to get an energy equation that depends only on |k|, integrate in shells over the solid angle
E(k) =
_
|k|=k
dE(
k) to get
E(k)
t
+ T(k) + F(k) = (k) (38)
where F(k) is some prescribed forcing. The phenomenology of the energy cascade and the formation of
the energy spectrum is based upon this equation.
Energy cascade:
Assume that the forcing and dissipation balance, F(k)dk = dk = .
Assume that F(k) is concentrated at small k
f
and (k) is concentrated at large k
.
Assume that between these regimes, T(k) is constant such that for k
f
<< k << k
, T(k) = ,
where is the ow of energy through wavenumber.
For any such k, assume that E(k) depends only on and k.
then by dimensional analysis, where is a dimensionless constant determined by experiments
E(k) =
2/3
k
5/3
(39)
36
4.8 Energy spectrum
Regimes : E(k) k
4
E(k) =
2/3
k
5/3
E(k) exp(2.1k)
Figure 4: k
5/3
energy spectrum. = k
2
E(k). Note that the Kolmogorov wavenumber is slightly higher
than the peak of the dissipation spectrum.
37
Spectral regimes and structure function scaling:
E(k) k
4
. Energy-containing regime: For k < 1/, the integral scale.
Either at if the forcing is at the largest scale of the ow, or a positive power law.
There exists a theoretical argument that it can be at most k
4
.
E(k) =
2/3
k
5/3
. Inertial-subrange. 1/ < k < 1/
Scaling theory from Kolmogorov (1941), no rigourous theory.
Supported by experiments and DNS for isotropic, homogeneous three-dimensional turbulence
The exponent is very close to -5/3
1.75 if k is the 3D wavenumber.
0.5 if k is the 1D wavenumber.
E(k) k
n
exp(2.1k) where n = 5/3. Dissipation regime. k > 1/
The mathematics of singularities of the Navier-Stokes equations in complex time suggests an
exponential dependence.
Supported by direct numerical simulations, from which the factor 2.1 comes.
Experiments; It is dicult to measure these small scales with hot-wire probes.
Structure functions: S
2
(r) =
2/3
r
2/3
in the inertial subrange > r > .
Viscous regime. S
2
(r) = r
2
/
2
Energy-containing regime S
2
(r) 2u
2
0
as r
In the hypothetical construct of 2D turbulence, a -5/3 spectrum is also predicted. But the sign of T(k)
will be reversed, energy ows up the spectrum. This might be important in geophysical ows.
38
5 Engineering applications and k models
More modelling, less computing
Zero-equation models:
ij
is set by mixing length
Most important advance: Van Driest (1956) damping function
as y 0 mixing length
m
= y[1 exp(y
+
/A
+
0
)] : A
+
0
= 26
One-equation model: energy equation with a mixing length
Two-equation models: Reynolds-averaged Navier-Stokes (RANS) with equations for:
k-. k: energy. : energy dissipation. Combined they give a mixing length = k
3/2
/.
Most widely used industrial model. Variation: k where = /k
DES: Detached Eddy Simulation:
Use log-law, or Van Driest, or your favourite model for boundary layer
Use LES in the interior
LES: Large-eddy simulation.
The eddy viscosity depends on ow variables and changes with time.
DNS: Direct-numerical simulation
More computing, no modelling
39
5.0.1 Numerics
u
t
u(t
0
+ t) u(t
0
)
t
= RHS
where the right-hand side depends upon all the nonlinear, pressure and viscous terms determined at
t = t
0
, that is
RHS = (u(t
0
)
)u(t
0
)
1
P(t
0
) +
2
u(t
0
)
This is a gross oversimplication. This time advancement scheme is know as forwards Euler. Any
attempt to calculate time-development with forwards Euler is known to yield unstable behaviour. A
course in numerics is needed to explain how these problems are resolved.
The reason for introducing this equation is to make you aware of what is actually done with the models
presented. They are used to calculate the RHS for the velocity, energy and dissipation equations so
that the state at a new time t
0
+ t can be determined. Besides the equations, we need additional
conditions. In, particular boundary conditions. The objective of these last lectures is to acquaint
you with the models and approximations that let you do this.
u
i
x
j
=
Non-linear terms
Pressure
Complicated boundaries
Uses
+ Finite-dierence vs.
+
u
i+1
u
i1
2x
+ Easy (u
)u
+ Dicult
+ Possible (u
1
u
0
)/x
+ All practical problems
* Fourier methods
Use FFTS and
= ik
j
u
i
(
k)
* Use FFTs on (u
)u
* Easy
* Impossible
* Idealised: LES, DNS
40
5.0.2 Large-eddy simulations (LES) versus Direct numerical simulations (DNS)
Complicated boundaries
- Re
Simple boundary conditions
- Re
Viscous terms
DNS
Finite dierence, mapping
- very low Reynolds numbers
Fourier or related methods
- high Reynolds number
Newtonian
2
u
LES
Models for wall stresses,
- moderate Re (DES)
same as DNS
- high Re
Smagorinsky
T
(x)
u
5.1 LES: Large-eddy simulation.
The eddy viscosity depends on ow variables and changes with time.
Most common parametrization: Smagorinsky:
T
= C
s
|e
ij
|
2
Usually C
s
= 0.2, and xed due to the property that:
If averaged over a given length scale, on average a turbulent ow dissipates energy
Backscatter is a property of turbulent ows where in some locations, energy ows from the unseen,
modelled subgrid scales to the resolved scales. (Explain averaging procedures)
For xed C
s
, boundary layer instabilities are suppressed and backscatter is underpredicted
All atmospheric mesoscale prediction models (called non-hydrostatic models) are LES models
with C
s
0.1 near the boundaries to compensate for backscatter in an ad hoc manner
Attempts to dynamically predict backscatter: Dynamic models.
At Warwick we are experimenting with adding extra terms. Non-linear models.
41
5.2 k- models
k equation :
k
t
+
U
k
x
i
_
_
k
k
x
i
_
_
=
equation :
t
+
U
x
i
_
_
x
i
_
_
= C
1
k
C
2
2
k
where transport (advective and diusive) is indicated by the terms in red. is the production term.
For the proles we have already considered, then more generally:
special = uv
U
y
=
T
_
_
U
y
_
_
2
generally =
1
2
T
_
_
_
U
i
x
j
+
U
j
x
i
_
_
_
2
where
T
= C
k
2
The underlying problem that these equations address is preventing k < 0. This will happen in
quasi-normal theories if is too large with respect to k.
The way this is addressed is if k is small, /t will be large and negative. That is will automatically
become small.
You still need the RANS equation:
U
t
+ (
U
)
U =
j
_
_
_
T
_
_
_
U
i
x
j
+
U
j
x
i
_
_
_
_
_
_ ,
U = 0, =
P
+
2
3
k
k- depends upon 5 experimentally determined coecients, whose typical values are:
C
= 0.09, C
1
= 1.44, C
2
= 1.92,
k
= 1.0,
= 1.3.
Boundary conditions can be dicult
Most complicated: stress equation models. Too many coecients for industrial applications.
42
5.3 Boundary layer regions
Flat-plate assumed.
Viscous sub-layer: U
x
y
+
Buer layer
Inner turbulent region: Law of the Wall
Outer wake-like region
Unsteady irrotational ow (not turbulent)
Only DNS can handle all of these regions.
Turbulence modelling ignores the inner two regions and starts with the Law of the Wall (with
modied coecients).
U
x
= u
(
1
log y
+
+ a) where y
+
=
u
is the dimensionless distance from the wall and a represents the physics from the boundary condition
whose eects are transmitted through the viscous sub-layer and the buer layers.
43
5.4 Boundary conditions
The k- model as introduced above does not address boundary layers.
Boundary layers must be addressed separately, then used as boundary conditions on a simulation
based upon the k- model.
No-slip boundary condition.
The velocity boundary condition at a surface is u = 0,
For special geometries such as a at plate, DNS can exactly represent this as accurately as
desired.
Almost everything we know about the details in the near-wall region comes from DNS.
In k- models, k = 0 and can be large.
Within the k- model, where y is the direction perpendicular to the wall, we could set
(wall) =
2
k
y
2
= 2
_
_
_
k
1/2
y
_
_
_
2
Modelled boundary conditions
We use models and scaling laws to create less exact boundary conditions, applied at the wall.
Law of the wall
U
x
= u
(
1
log
u
+ a) U
y
0 (40)
k
u
2
C
1/2
u
3
y
2
= (C
2
C
1
)
C
1/2
= 0.43
2
consistent with the experimental value of = 0.39 because these k- coecients are taken from the
same experiments.
44
5.5 Modelling with k- plus log-layer boundary conditions
Advantages of k- modelling over LES due to a smoother larger ow (U):
k- requires enhanced resolution near walls only in the direction normal to the wall.
Time-steps can be relatively long.
You need to choose y and u
for the boundary conditions. This must be done empirically. For example,
the skin-friction coecient is dened as
c
f
w
1
2
U
2
where
w
= uv = u
2
(41)
For a at plate, we know empirically that c
f
= (2 log
10
Re
x
0.65)
2.3
k- solutions are dependent upon
the physical Reynolds number and viscosity only through the boundary conditions given by (40,41).
The choice of y is a balancing act.
Choosing y large puts more of the dissipation into the prescribed boundary layer, the least reliable
part of the model.
Choosing y small makes at the wall large and makes the numerical modelling more dicult due
to three considerations:
1. More grid points are needed in the direction perpendicular to the wall.
2. Stronger velocity gradients in the boundary layer imply that the numerical time-step, determined
by following condition over the ow: t = 1/ max(U/y), is a small value determined in the
boundary layer. This, despite the relatively long time-steps that the interior ow would allow.
3. This is known as the stiness problem and exists in many physical systems where the timescale of
interest is very large, but the physics occurs on very short timescales.
45
5.6 Pitfalls of k-
Separation: k- assumes a at-plat boundary layer which thickens and increases drag as Re increases.
For a curved surface there can be separation (a vortex is shed) and drag decreases dramatically.
Wakes: k- predicts a sharp interface between the turbulent and irrotational uid. In reality, the
interface uctuates and on the average proles are smooth.
Singular behaviour of U and at the boundary implied by (40).
k- is something of an art.
6 Summary of assumptions, their limitations, and review
Reynolds averaging
1. Large-scale time variations are small. This approximation is valid only for simple ows like a
wake, jet or mixing layer. For real applications, unsteadiness should be considered.
2. Downstream scaling obeys simple power laws. Even for simple ows recent work indications
that minor modications are necessary. For complicated ows these assumptions are not true
and are only useful for matching regimes.
3. Know the terms in the RANS equations: advection, production, transport, dissipation.
Statistical theory of turbulence
1. Homogeneity: All statistics can be represented by simple separations or wavenumbers in each
direction. Multi-component hot-wire data is applicable using the Taylor frozen turbulence as-
sumption. This assumption is approximately true in the centre of simple ows, far from any
boundaries. It is not true at all for boundary layers.
46
2. Isotropy: All statistics can be represented by a single direction. Single-component hot-wire
statistical data can be converted into other components. Experiments show that this is never
true.
3. Gaussian statistics: One can use the quasi-normal approximation, or that approximation mod-
ied to preserve realisibility with E > 0 for all times, to close the energy energy. Experiments
show this is never true.
4. Periodic boundary conditions: Simulations in a periodic box are applicable. What is learned
can only be applied where homogeneity is approximately true.
5. Energy cascade depends only upon the energy cascade rate, which is equivalent to the energy
dissipation rate . If this were approximately true it would imply that even if isotropy and
homogeneity were not true on the largest scales, they would be true on smaller scales. This is
an important assumption in large-eddy simulation.
6. Know the spectral regimes: energy-containing, inertial, viscous, and their characteristic length
scales: Large-scale L, Taylor microscale , Kolmogorov . You do not need to know details
about moments, correlations and structure functions, these would be provided on a test or exam.
7. No-slip assumption. True for most real surfaces. Rough surfaces need modications. But for
anything other than direct numerical simulations, using in practical calculations is not feasible.
k- modelling
1. That the dissipation equation can be closed using factors of /k.
2. That there are sharp boundaries between the turbulent ows and the surrounding irrotational
ow. This is actually a consequence of k- modelling, that is the proles predicted have this
property. This property is never true due to intermittency, the tendency for there to be al-
ternating regions of highly rotational, turbulent ow separated by irrotational, nearly laminar
uid.
47
3. Flat-plate assumption: That there is no curvature on surfaces and the wake-like region down-
stream grows with Reynolds number. In general, surfaces are curved and at large enough
Reynolds numbers there is separation. That is vortices are shed, the wake-like region becomes
smaller, and the drag less.
4. No-slip assumption: U = k = 0 at Very dicult to apply in practice.
5. Log-layer assumption: that the scaling of the log-layer can be used to provide the boundary
condition. This is the weakest part of the complete RANS/k- modelling framework. The
thicker the layer, the more the modelled dissipation can lead to errors. The thinner the layer,
the more numerical errors can accumulate. Stiness problem.
6. Know the layers of a at-plate boundary layer: viscous sublayer, (buer layer), log-law, turbulent
wake-like, irrotational outer ow.
Modelling: pros and cons:
1. Know the advantages and disadvantages of dierent types of modelling of turbulent ows.
2. Know situations where each type of modelling is most eective.
3. You are not expected to know the details of either DNS or the latest advances in LES (only
what the Smagorinsky model is).
4. Know the balancing between more modelling/less computing and less modelling/more comput-
ing. Both between methods in general,
5. and specics about k-.
First lectures on compressible Navier-Stokes and related material.
1. Know how to do the example problems for which answers have been provided. Nothing more.
48
7 Coherent structures
7.1 Laminar-turbulence transition (Kelvin-Helmholtz)
7.2 Mixing layer
7.3 Boundary layer
8 Vortex dynamics
49