Turbulence Modelling

(CFD course)

Sławomir Kubacki
[Link]@[Link]

14.11.2016

Copyright  2016, Sławomir Kubacki

Turbulence Modelling Sławomir Kubacki

Contents

1. Reynolds-averaged Navier-Stokes equations ...................................................................................... 3

2. Closure of the modelled terms ............................................................................................................ 9

2.1. Exact Reynolds-stress transport equations .................................................................................. 9

2.2 Transport equation for turbulent kinetic energy ........................................................................ 10

2.3. The standard k- model .............................................................................................................. 11

References ............................................................................................................................................. 17

Appendix A ............................................................................................................................................ 19

2
Turbulence Modelling Sławomir Kubacki

1. Reynolds-averaged Navier-Stokes equations

The discussion of the Reynolds or time-averaged Navier-Stokes (RANS) equations and

turbulence transport equations is limited to the constant-density (incompressible)
fluids. A linear relationship is assumed between the components of the stress and
deformation tensors. An extension to compressible fluids is straightforward and can be
found in many textbooks (see for instance Wilcox, 2006, 2008). In case of compressible
fluid one has to apply the Favre-averaging (density-based averaging) instead of
mentioned Reynolds-averaging to allow for some simplifications in the mass and
momentum conservation equations.
For incompressible fluid, the instantaneous velocity component u i (x, t ) can be
written as the sum of a mean u i ( x) and a fluctuating part ui (x, t ) (fig. 1):

u i x, t   u i x   ui x, t  (1.1)

where the mean velocity is defined as the time-averaged value

Fig. 1 Decomposition of the signal in mean and fluctuating components.

t T
u i x   lim  u x, t dt
1
T  T
i (1.2)
t

T is the averaging time interval. We assume that T. This corresponds to the steady
RANS model. For time-accurate RANS (called URANS) this time interval has to be
sufficiently large with respect to the time scale T1 of the turbulent fluctuations (fig. 2)
and small with respect to the time scale T2 of large scale unsteadiness. So for time-
accurate problems Eq. (1.2) takes the form:
3
Turbulence Modelling Sławomir Kubacki

t T
u i x, t    u x, t dt
1
i
T t

Fig. 2 Time-averaging windows. T1 – time scale of turbulent fluctuations, T2 – time scale

of unsteady motion. The time-averaging window T should be: T1 < T < T2

Two useful properties are (steady RANS):

- the time-average of a time-averaged value is again the same value (the averaged
value
u i ( x) is not a function of time)

u i x 
t T t T
u i x   lim u i x dt  lim
1
T  T 
t
T  T t dt  u i (x) (1.3)

- and that the time-averaged value of the fluctuating part is zero

t T
ui x   lim [u i (x, t )  u i x ]dt  u i (x)  u i x   0
1
T  T 
(1.4)
t

The other averaging properties are listed below for arbitrary quantities ,  and :

     (time averaging is linear) (1.5)

4
Turbulence Modelling Sławomir Kubacki

 
 (1.6)
x i x i

   (1.7)

       (1.8)

   0 (1.9)

    0 (1.10)

    0 (1.11)

The Reynolds-averaging is linear (1.5) and it commutes with the space derivatives (1.6).
Relations (1.7) and (1.8) come from the property (1.3). With (1.9) and (1.10) we take
advantage of the observation that the product of mean and fluctuating quantity is zero
because the mean of the latter is zero. The quantities  and  are correlated. It means
the average of their product is not zero (Eq. 1.11).
For incompressible fluid the conservation of mass and momentum equations (with
the constitutive relation introduced already, see first lecture) are

u i
0
x i (1.12)

u i  u i u j  1 p 
t

x j
 
 x i x j
2Sji  (1.13)

where  is the fluid density,  is the dynamic molecular viscosity, p is the pressure and
Sij is the strain-rate tensor Sij=1/2(ui/xj+uj/xi).

First, we average the continuity equation (1.12). Taking the property (1.6) we obtain:

u i  u i
 .
x i x i
(1.14)

Next, we average the momentum equations (1.13).

u i u i u j  1 p 
t

x j
 
 x i x j
2Sji 
Using the property (1.5) we obtain
5
Turbulence Modelling Sławomir Kubacki

u i u i u j  1 p 
t

x j
 
 x i x j
2Sji 
(1.15)
Here we assume that the time-averaging commutes with the local time derivative. We
can see from figure 1 that for T   the time-derivative of the mean velocity is zero

u i  u i
 0
t t (1.16)

The convective derivative deserves attention. We average the product of ui and uj:

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 i u j (1.17)

The second and third term on r.h.s. of eq. (1.17) cancels out (property 1.9). The last term
is nonzero because there is some correlation between ui and uj fluctuating velocity
components. So using property (1.7) we obtain:

u i u j  u i u j  ui uj (1.18)

Figure 3 shows the momentum exchange process in the shear layers of the plane jet at
Re=20000. This example is used to illustrate a correlation between u’ and v’ fluctuating
velocity components. Let us assume that the red fluid element (located on line R) is
shifted to the left. This shift is visible by negative u’ fluctuating velocity component (see
the coordinate system). Since the red fluid element is moving from the high mean
velocity zone VR to the small mean velocity zone VL this transfer generates the positive v’
fluctuating velocity
(v’=VR-VL, v’>0). The product of these two fluctuations is negative u’v’<0. The time-
averaged product of u’ and v’ is shown in fig. 3 (bottom). This is the turbulent shear
stress profile ii  ui uj . We clearly see that the negative u’ fluctuation generates the
positive v’ fluctuation. The momentum exchange process can also be realized the other
way around, so from left to right (see blue fluid element located on line L). This is due to
the fact that formation of the shear layer is typically related to evolution (in space and
time) of the coherent vortex structures which subsequently breakdown into smaller
forms. It means that at one time instant there is a momentum exchange from right to
left. But at the other time instance there might be the momentum exchange from left to
right. Let us, therefore, assume that the fluid element is shifted from left to right (blue
contour). This shift results in positive u’ fluctuating velocity. Since the blue fluid element
is moving from the low mean velocity zone VL to the high mean velocity zone VR this
transfer generates the negative v’ fluctuating velocity (v’ =VL-VR, v’<0). The product of
these two fluctuating velocities is, again, negative u’v’<0. It means that it leads to
generation of the negative shear stress after time averaging (fig. 3 bottom).
The Reader can perform a similar analysis for the right part of fig. 3. The result of this
analysis will be a formation of the positive shear stress profile (fig. 3, bottom). The sign
6
Turbulence Modelling Sławomir Kubacki

depends on the coordinate system. In both cases the net effect is the momentum transfer
from the mean flow to the fluctuating flow (resulting in increased turbulent shear
stress). But the most important remark is that in turbulent flow the momentum
exchange in one flow direction causes the momentum exchange in the other flow
directions (flow is three dimensional). So there is some correlation between both ui’ and
uj’ fluctuating components. This results in the shear stress tensor (rightmost term in Eq.
1.18).
Going back to the time-averaging. The terms on r.h.s. of equation (1.15) can be easily
averaged:

1 p 1 p  
 
x j
2S ji  
x j
2S ji 
 x i  x i
(1.19)

Finally, we obtain the Reynolds-Averaged Navier-Stokes (RANS) equations (eq. 1.14,

1.18, 1.19):

ui
0
x i (1.20)

 (u j u i ) 1 p 
x j
 
 x i x j
2S ji  ujui  (1.21)

where
1   u i  u j 
Sij   (1.22)
2  x j x i 

The last term on r.h.s of Eq. (1.21) is the divergence of the Reynolds-stress tensor. The
components of the Reynolds –stress tensor are denoted by  ij  ui uj .This term
requires closure model (see below).

7
Turbulence Modelling Sławomir Kubacki

L R
v’>0
u’<0

u’>0

v’<0
Y

Fig. 3. Momentum exchange in the shear layers of the plane jet. Mean V velocity
component (top) and the shear stress profile (bottom).

8
Turbulence Modelling Sławomir Kubacki

2. Closure of the modelled terms

2.1. Exact Reynolds-stress transport equations

If we denote by N(ui) the Navier-Stokes operator:

u i u i 1 p  2u i
N( u i )   uk   0 (2.1)
t x j  x i x k x k

We can obtain the exact Reynolds-stress transport equation, by multiplying the

momentum equation (2.1) by fluctuating quantity uj and adding to this term a similar
term with the indexes interchanged and averaging. The following transport equation is
obtained (Wilcox, 2006):

 ij  ij u j ui u uj 1  ui uj 

 uk   ik   jk  2 i  p  
t x k x k x k x k x k   x j x i 
(2.2)
   ij 1 1 
   ui ujuk  pui jk  puj ik .
x k  x k   

Derivation of Eq. (2.2) is show in Appendix A.

The physical meaning of the terms on r.h.s. of Eq. (2.2) can be described as follows.

- Production (first and second term). Turbulent stresses are generated at the expense of
mean flow energy by mean flow deformation. This term does not need any closure.

- The third term represents the stress dissipation which mainly occurs at the smallest
scales (it should be modelled).

- Pressure fluctuations (fourth term) redistribute the turbulent stress among

components to make turbulence more isotropic.

- Transport term (last term in Eq. 2.2). This term consists of several parts (in square
bracket):

First term is the transport term and it is called the molecular diffusion.

Second, there is the turbulent diffusion term (transport through velocity fluctuations). A
closure model is necessary for this term.

A third and fourth part of the transport term is the pressure transport. This term also
needs closure model.
9
Turbulence Modelling Sławomir Kubacki

2.2 Transport equation for turbulent kinetic energy

An exact equation for the turbulent kinetic energy follows from the equations for the
Reynolds stress components (2.2) by contraction trough putting j = i and making the
sum over i = 1,2,3:

(2.3)

We take advantage of the continuity equation for the fluctuating velocity components.

(2.4)

which means that for an incompressible fluid studied here the pressure strain-rate term
cancels out in Eq. (2.3). We define the turbulent kinetic energy as:

(2.5)

Next, by inserting Eq. (2.5) to Eq. (2.2) we obtain an exact equation for the turbulent
kinetic energy

(2.6)

In Eq. (2.6) the first term on r.h.s describes production of the turbulent kinetic energy.
This term is modeled using the Boussinesq hypothesis:

2
 ik  2 t Sik  k ik (2.7)
3

We later write the production term as Pk  2 t Siku i / x k  2 t Sik Sik . t in Eq. (2.7)
denotes the turbulent viscosity (see discussion below).

The second term on r.h.s. of eq. (2.6) represents dissipation of the turbulent kinetic
10
Turbulence Modelling Sławomir Kubacki

energy

(2.8)

In the frame of the two-equation models discussed here the term (2.8) is obtained by
solving an additional transport equation (discussion below).

The second and third term in brackets describe the turbulent transport process by
velocity fluctuations and the fluctuating pressure–velocity induced diffusion. Both terms
are modeled using the gradient hypothesis:

(2.9)

where k is a certain constant. Finally, the transport equation for the turbulent kinetic
energy reads:

(2.10)

2.3. The standard k- model

The model which is called the standard k- turbulence model was developed by Jones
and Launder (1972). In this model, the RANS-equations are used together with the k-
equation (2.10), the -equation (introduced below) and the eddy-viscosity expression
k2 k2
 t  C . The eddy viscosity expression  t  C  requires the knowledge of the
 
turbulent kinetic energy, which we can derive from equation (2.10). The second
ingredient is the dissipation rate . Remark that we are not restricted to an equation of .
Any quantity which is a combination of k and  may be used. The most obvious is to use
an equation for . The expression for  is known, Eq. (2.11). So, we could use an
approach similar to the one as used for the turbulent kinetic energy to derive a transport
equation. Technically, this is possible, but the equation is much more complicated than
the equation for k and contains many more terms that need modelling. So, looking at the
k-equation (2.10) and taking into account that the parameter that we can form with k
and  is the time scale =k/, the k-equation can be transformed into an -equation by

  c P  c 2       
 uk  1    t  (2.11)
t x k  x k      x k 

where we introduce the constants cε1 and cε2 in the transformation of the production
11
Turbulence Modelling Sławomir Kubacki

term and the dissipation term. We further introduce the diffusion coefficient  (see
discussion latter).
The standard values for the model parameters are:

c  0.09, k  1.0,   1.3, c1  1.44, c 2  1.92

The value for c comes from the observation of thin shear flows with approximate
balance between production and dissipation: 2D free jet mixing layers and the so-called
inertial region or logarithmic layer in a boundary layer flow. Figure 4 shows the terms in
the turbulent kinetic energy equation (Eq. 2.6) along normal to the wall direction in the
turbulent boundary layer at Re=1410. The results have been obtained using DNS by
Laurent et al. (2012). The molecular diffusion is important at the wall. Both production
and dissipation become dominant at y+>5. Note a sign change on the profiles of the
turbulent and molecular diffusion inside the boundary layer (diffusion partly acts as a
sink and partly as a source of the turbulent kinetic energy). This makes the production
and dissipation far more important than the diffusion processes at sufficiently large
distance from the wall. Moreover, we can assume that production almost equals
dissipation for y+>10: Pk=ε. With y as cross-stream coordinate this results in:

u 2
t ( )
 y u ' v ' 2
c   t 2   t 2
( )
k k k

For thin shear flows, the experimental observation is ( u ' v ' )  0.3 , resulting in c = 0.09.
k

production

dissipation

Figure 4. Terms in the turbulent kinetic energy equation inside the turbulent boundary
layer at Re=1410. DNS results by Laurent et al., (2012).

12
Turbulence Modelling Sławomir Kubacki

The constant c2 is determined by considering decaying homogeneous turbulence.

Figure 5 shows the vortex structures obtained for DNS simulation of decaying isotropic
and homogeneous turbulence in the periodic box (Dubief i Delcayre, 2000). In this
interesting flow the statistics of all quantities in brackets on r.h.s. of Eq. (2.6) are
constant. This means that there is no diffusion. Moreover, the statistics of mean velocity
derivatives are zero. It means that there is no production either. The k-equation reduces
to

Dk D 2
  ,  c  2 (2.12)
Dt Dt k

These equations are satisfied for an observer moving with the flow by

1
k  t  n ,   t  m , dla n  , m  n 1 (2.13)
c 2  1

Experiments of Comte-Bellot and Corrsin (1966) indicate values n = 1.2 to 1.3 which
implies c2 = 1.83 to 1.77. The value c2 = 1.92, selected by Jones and Launder differs
somewhat because they did some numerical optimisation of this constant over a range
of flows.

Figure 5. Homogeneous and isotropic turbulence in a periodic box. (Dubief i Delcayre,

2000)

13
Turbulence Modelling Sławomir Kubacki

The constant c1 is determined from a homogeneous shear flow experiment. Figure 6
shows a sketch of such flow. Homogeneous shear flow appears to reach an equilibrium
state with k and  growing in such a manner that the turbulent time-scale k/
approaches an approximately constant value. The governing equations are

Dk D c1Pk  c 2
 Pk   ,  (2.14)
Dt Dt k/

which can be combined with the assumption of a constant turbulent time-scale to yield
the following relations. Assuming D/Dt=0 we arrive at:

D(k / ) Dk 1 k D
 
Dt Dt  2 Dt
Pk k c P  c 2
  1  2 ( 1 k   2 )
  k k (2.15)
P
 c 2  1  (c1  1) k  0


which leads to

c 2  1
c1  1
Pk / 

Using c2 =1.83 and the shear flow data of Tavoularis and Corrsin (1981), where
Pk /   1.8 and k /   constant, the value c1  1.46 results. Again, Jones and Launder
take the somewhat different value c1  1.44 by numerical optimisation.

Figure 6. Flow in 2D mixing shear layer

14
Turbulence Modelling Sławomir Kubacki

The diffusion coefficient σk is a priori taken to unity, since it governs the turbulent
diffusion of the turbulent kinetic energy by the turbulent motion itself.
The second diffusion coefficient σε comes from consideration of the diffusion of ε in
the logarithmic zone of the boundary layer. The following equations can be written in a
zero pressure gradient boundary layer flow:

 u u
0 ( t )   t  cons tan t  w  u 2 (2.16)
y y y

u 2   k
0  t ( )  ( t ) (2.17)
y y  k y

 u 2  2   t 
0  c1  t ( )  c 2  ( ) (2.18)
k y k y  y

In the log-layer,

du u 
 (2.19)
dy y

From Eq. (2.16) and (2.19) we obtain

 t  u y (2.20)

If we assume Pk= and the constant k-profile in the logarithmic part of the turbulent
boundary layer (no diffusion of k) Eq. (2.10) reduces to:

u u
Pk    u ' v '   t ( )2
y y
u u
    t ( )2  u  y(  )2
y y

u 3
 (2.21)
y

k2
Taking the definition of the dynamic turbulent viscosity  t  C and introducing it

into Eq. (2.16) we obtain:

15
Turbulence Modelling Sławomir Kubacki

u k2 u
w  u 2  u ' v '   t  c
y  y
2
k2 u k2 u  c k
 u 2  c  c 3 y  2
 y u y u

u 2
k (2.22)
c

k2
Introducing Eq. (2.19), (2.21), (2.22) and  t  C  into Eq. (2.18) we obtain:


u 3 c u 2 u 6 c   u  y u 3 
0  c1 u  y  c 2   ( 2 ) 
y u 2 ( y)2 ( y)2 u 2 y   y 
u 4 c u 4 c  u4
 c1  c 2  (  )
( y)2 ( y)2 y  y
 c 1 1  4
  (c1  c 2 )  u
2 
 ( y) 2 
 y 


This results in

2
  (2.23)
(c 2  c1) c

In standard k-, the values c1 = 1.44, c2 = 1.92, c = 0.09 and  = 1.3 are used, which
following (2.23) implies  = 0.43, which is a somewhat larger value than what is usually
assumed for the Von Karman constant ( = 0.40-0.41). Again, this difference comes from
the numerical optimisation by Jones and Launder.
One of the shortcoming of the standard k- model is its inability to reproduce a
correct level of the turbulent shear stress in the near-wall region. Figure 7 shows the
predicted by the k- model (dashed line) and computed using DNS (solid line) profiles of
the turbulent viscosity inside the boundary layer for simulation of the channel flow. As
shown the k- model overestimates the turbulent viscosity in the near-wall region. In
order to limit this shortcoming some damping terms are employed in front of the eddy-
viscosity formula. The other approach consists in using separate turbulence model,
which shows better performance than the classic k- model in the near-wall region, and
its blending with standard k- further away from walls. This obviously complicates the
modeling approach using the standard k- model.

16
Turbulence Modelling Sławomir Kubacki

Figure 7. Predicted using the k- model (dashed line) and computed using DNS (solid
line) profiles of the turbulent viscosity inside the boundary layer for simulation of the
channel flow at Re=590. Durbin and Pettersson Reif, (2003).

References

J.B. Cazalbou, P.R. Spalart and P. Bradshaw. On the behaviour of two-equation models at
the edge of a turbulent region. Physics of fluids, 6(5):1797-1804, 1994.

G. Comte-Bellot and S. Corrsin. The use of a contraction to improve the isotropy of grid-
generated turbulence. J. Fluid. Mech., 25:657-682, 1966.

Y. Dubief and F. Delcayre, On coherent-vortex identification in turbulence, Journal of

Turbulence, 1, N11, 2000.

P.A. Durbin and B.A. Pettersson Reif. Statistical Theory and Modeling for Turbulent
Flows. John Wiley, 2001.

W.P. Jones and B.E. Launder. The prediction of laminarization with a two-equation
model of turbulence. AIAA J., 15:301-314,1972.

B.E. Launder and B.I. Sharma. Application of the Energy Dissipation Model of Turbulence
to the Calculation of Flow Near a Spinning Disc. Letters in Heat and Mass Transfer,
1(2):131-138, 1974.

17
Turbulence Modelling Sławomir Kubacki

C. Laurent, I. Mary, V. Gleize, A. Lerat, D. Arnal, DNS database of a transitional separation

bubble on a flat plate and application to RANS modeling validation, Computers & Fluids
61: 21–30, 2012.

Tavoularis and Corrsin., Experiments in Nearly Homogeneous Turbulent Shear Flow

with a Uniform Mean Temperature Gradient. J. of Fluid Mechanics, 104:311-347, 1981.

D.C. Wilcox. Reassessment of the Scale Determining Equation for Advanced Turbulence
Models. AIAA J., 26(11):1299-1310, 1988.

D.C. Wilcox. Comparison of Two-Equation Turbulence Models for Boundary Layers with
Pressure Gradient. AIAA J., 31(8):1414-2031, 1993.

D.C. Wilcox. Turbulence Modeling for CFD. Griffin Printing, Glendale, California, 1993.

D.C. Wilcox. Turbulence Modeling for CFD. DCW Industries, 2006.

D.C. Wilcox, Formulation of the k- turbulence model revisited. AIAA J., 46: 2823-2837,
2008.

18
Turbulence Modelling Sławomir Kubacki

Appendix A

Multiplying Eq. (2.1) by the fluctuating velocity uj and adding this term to a similar
term with the indexes interchanged and averaging we obtain:

(A.1)

After some algebra we obtain:

(A.2)

Assuming the Reynolds decomposition (mean plus fluctuation):

(A.3)

and putting Eq. (A.3) to (A.2) we obtain several terms which are discussed below:

1. Local time derivative

(A.4)

19
Turbulence Modelling Sławomir Kubacki

The property (1.10) was used in order to simplify Eq. (A.4) (second line). The following
relation was used to derive the last term in Eq. (A.4)

(A.5)

2. Convective term

(A.6)

3. Pressure gradient term

(A.7)

4. Viscous term

(A.8)

Finally, putting Eq. (A.4), (A.6), (A.7) and (A.8) to Eq. (A.2) we obtain the exact form of
the Reynolds-stress equation:
20
Turbulence Modelling Sławomir Kubacki

(A.9)

Note that Eq. (A.9) can be rewritten in the form given by Eq. (19) using the following
relations:

(A.10)

21