Fluid mechanics,turbulent ow and turbulence
modeling
Lars Davidson
Division of Fluid Dynamics
Department of Applied Mechanics
Chalmers University of Technology
SE412 96 G¨oteborg,Sweden
http://www.tfd.chalmers.se/lada,lada@chalmers.se
August 25,2013
Abstract
This course material is used in two courses in the International Master's pro
gramme Applied Mechanics at Chalmers.The two courses are TME225 Mechan
ics of uids,and MTF270 Turbulence Modeling.MSc students who follow these
courses are supposed to have taken one basic course in uid me chanics.
This document can be downloaded at
http://www.tfd.chalmers.se/lada/MoF/lecture
notes.html
and
http://www.tfd.chalmers.se/lada/comp
turb
model/lecture
notes.html
The Fluid courses in the MSc programme are presented at
http://www.tfd.chalmers.se/lada/msc/mscprogramme.html
The MSc programme is presented at
http://www.chalmers.se/en/education/programmes/mastersinfo/Pages/AppliedMechanics.aspx
1
Contents
1 Motion,ow
9
1.1 Eulerian,Lagrangian,material derivative
...............9
1.2 Viscous stress,pressure
........................10
1.3 Strain rate tensor,vorticity
.......................11
1.4 Product of a symmetric and antisymmetric tensor
..........13
1.5 Deformation,rotation
.........................14
1.6 Irrotational and rotational ow
....................16
1.6.1 Ideal vortex line
........................17
1.6.2 Shear ow
...........................18
1.7 Eigenvalues and and eigenvectors:physical interpretation
......19
2 Governing ow equations
21
2.1 The NavierStokes equation
......................21
2.1.1 The continuity equation
....................21
2.1.2 The momentumequation
...................21
2.2 The energy equation
..........................22
2.3 Transformation of energy
.......................23
2.4 Left side of the transport equations
..................24
2.5 Material particle vs.control volume (Reynolds Transport Theorem)
.25
3 Exact solutions to the NavierStokes equation:two examples
27
3.1 The Rayleigh problem
.........................27
3.2 Flow between two plates
........................30
3.2.1 Curved plates
.........................30
3.2.2 Flat plates
...........................31
3.2.3 Force balance
.........................33
3.2.4 Balance equation for the kinetic energy
............34
4 Vorticity equation and potential ow
36
4.1 Vorticity and rotation
.........................36
4.2 The vorticity transport equation in three dimensions
.........37
4.3 The vorticity transport equation in two dimensions
..........40
4.3.1 Boundary layer thickness fromthe Rayleigh problem
....41
5 Turbulence
43
5.1 Introduction
..............................43
5.2 Turbulent scales
............................44
5.3 Energy spectrum
............................45
5.4 The cascade process created by vorticity
...............49
6 Turbulent mean ow
53
6.1 Time averaged NavierStokes
.....................53
6.1.1 Boundarylayer approximation
................54
6.2 Wall region in fully developed channel ow
.............54
6.3 Reynolds stresses in fully developed channel ow
..........58
6.4 Boundary layer
.............................60
7 Probability density functions
62
2
3
8 Transport equations for kinetic energy
65
8.1 The Exact k Equation
.........................65
8.1.1 Spectral transfer dissipation ε
κ
vs.true viscous dissipation,ε
68
8.2 The Exact k Equation:2D Boundary Layers
.............69
8.3 Spatial vs.spectral energy transfer
..................70
8.4 The overall effect of the transport terms
................71
8.5 The transport equation for ¯v
i
¯v
i
/2
...................72
9 Transport equations for Reynolds stresses
74
9.1 Reynolds shear stress vs.the velocity gradient
............79
10 Correlations
81
10.1 Twopoint correlations
.........................81
10.2 Auto correlation
............................83
11 Reynolds stress models and twoequation models
85
11.1 Mean ow equations
..........................85
11.1.1 Flow equations
........................85
11.1.2 Temperature equation
.....................86
11.2 The exact
v
′
i
v
′
j
equation
........................86
11.3 The exact
v
′
i
θ
′
equation
........................88
11.4 The k equation
.............................90
11.5 The ε equation
.............................91
11.6 The Boussinesq assumption
......................92
11.7 Modelling assumptions
........................93
11.7.1 Production terms
.......................93
11.7.2 Diffusion terms
........................94
11.7.3 Dissipation term,ε
ij
.....................95
11.7.4 Slow pressurestrain term
...................96
11.7.5 Rapid pressurestrain term
..................99
11.7.6 Wall model of the pressurestrain term
............104
11.8 The k −ε model
............................106
11.9 The modeled
v
′
i
v
′
j
equation with IP model
..............107
11.10 Algebraic Reynolds Stress Model (ASM)
...............107
11.11 Explicit ASM(EASMor EARSM)
..................108
11.12 Boundary layer ow
..........................109
12 Reynolds stress models vs.eddyviscosity models
111
12.1 Stable and unstable stratication
...................111
12.2 Curvature effects
............................112
12.3 Stagnation ow
............................115
12.4 RSM/ASMversus k −ε models
...................116
13 Realizability
117
13.1 Twocomponent limit
.........................118
14 Nonlinear Eddyviscosity Models
120
15 The V2F Model
123
15.1 Modied V2F model
..........................126
4
15.2 Realizable V2F model
.........................127
15.3 To ensure that v
2
≤ 2k/3 [1]
.....................127
16 The SST Model
128
17 Overviewof RANS models
132
18 Large Eddy Simulations
133
18.1 Time averaging and ltering
......................133
18.2 Differences between timeaveraging (RANS) and space ltering (LES)
134
18.3 Resolved &SGS scales
........................135
18.4 The boxlter and the cutoff lter
..................136
18.5 Highest resolved wavenumbers
....................137
18.6 Subgrid model
.............................138
18.7 Smagorinsky model vs.mixinglength model
.............138
18.8 Energy path
..............................139
18.9 SGS kinetic energy
..........................139
18.10 LES vs.RANS
.............................140
18.11 The dynamic model
..........................141
18.12 The test lter
..............................141
18.13 Stresses on grid,test and intermediate level
..............142
18.14 Numerical dissipation
.........................144
18.15 Scalesimilarity Models
........................145
18.16 The Bardina Model
..........................145
18.17 Redened terms in the Bardina Model
................145
18.18 A dissipative scalesimilarity model.
.................146
18.19 Forcing
.................................148
18.20 Numerical method
...........................148
18.20.1 RANS vs.LES
........................149
18.21 Oneequation k
sgs
model
.......................150
18.22 Smagorinsky model derived fromthe k
sgs
equation
.........150
18.23 A dynamic oneequation model
....................151
18.24 A Mixed Model Based on a OneEq.Model
.............152
18.25 Applied LES
..............................152
18.26 Resolution requirements
........................152
19 Unsteady RANS
154
19.1 Turbulence Modelling
.........................157
19.2 Discretization
.............................157
20 DES
160
20.1 DES based on twoequation models
..................160
20.2 DES based on the k −ω SST model
.................161
21 Hybrid LESRANS
163
21.1 Momentumequations in hybrid LESRANS
.............167
21.2 The equation for turbulent kinetic energy in hybrid LESRANS
...167
21.3 Results
.................................167
22 The SAS model
170
22.1 Resolved motions in unsteady
....................170
5
22.2 The von K´arm´an length scale
.....................170
22.3 The second derivative of the velocity
.................172
22.4 Evaluation of the von K´arm´an length scale in channel ow
.....172
23 The PANS Model
175
24 The PITMmodel
179
24.1 RANS mode
..............................179
24.2 LES mode
...............................179
25 Hybrid LES/RANS for Dummies
183
25.1 Introduction
..............................183
25.1.1 ReynoldsAveraging NavierStokes equations:RANS
....183
25.1.2 Large Eddy Simulations:LES
................183
25.1.3 Zonal LES/RANS
.......................184
25.2 The PANS k −ε turbulence model
..................184
25.3 Zonal LES/RANS:wall modeling
...................185
25.3.1 The interface conditions
....................185
25.3.2 Results
.............................186
25.4 Zonal LES/RANS:embedded LES
..................186
25.4.1 The interface conditions
....................186
25.4.2 Results
.............................188
26 Inlet boundary conditions
189
26.1 Synthesized turbulence
........................189
26.2 Randomangles
.............................189
26.3 Highest wave number
.........................190
26.4 Smallest wave number
.........................190
26.5 Divide the wave number range
....................190
26.6 von K´arm´an spectrum
.........................191
26.7 Computing the uctuations
......................191
26.8 Introducing time correlation
......................191
27 Overviewof LES,hybrid LESRANS and URANS models
194
28 Best practice guidelines (BPG)
197
28.1 EU projects
..............................197
28.2 Ercoftac workshops
..........................197
28.3 Ercoftac Classical Database
......................198
28.4 ERCOFTAC QNET Knowledge Base Wiki
..............198
A TME225:ǫ −δ identity
199
B TME225 Assignment 1:laminar ow
200
B.1 Fully developed region
........................201
B.2 Wall shear stress
............................201
B.3 Inlet region
...............................201
B.4 Wallnormal velocity in the developing region
............202
B.5 Vorticity
................................202
B.6 Deformation
..............................202
B.7 Dissipation
...............................202
6
B.8 Eigenvalues
..............................202
B.9 Eigenvectors
..............................203
C TME225:Fourier series
204
C.1 Orthogonal functions
.........................204
C.2 Trigonometric functions
........................205
C.3 Fourier series of a function
......................207
C.4 Derivation of Parseval's formula
...................207
C.5 Complex Fourier series
........................209
D TME225:Why has the energy spectrum,E,such strange dimensions?
210
D.1 An example
..............................210
D.2 An example:Matlab code
.......................212
E TME225 Assignment 2:turbulent ow
216
E.1 Time history
..............................216
E.2 Time averaging
............................216
E.3 Mean ow
...............................217
E.4 The timeaveraged momentumequation
...............217
E.5 Wall shear stress
............................218
E.6 Resolved stresses
...........................218
E.7 Fluctuating wall shear stress
......................218
E.8 Production terms
............................218
E.9 Pressurestrain terms
..........................218
E.10 Dissipation
...............................219
E.11 Do something fun!
...........................219
F TME225 Learning outcomes 2012
220
G MTF270:Some properties of the pressurestrain term
231
H MTF270:Galilean invariance
232
I MTF270:Computation of wavenumber vector and angles
234
I.1 The wavenumber vector,κ
n
j
......................234
I.2 Unit vector σ
n
i
.............................235
J MTF270:1D and 3D energy spectra
236
J.1 Energy spectra fromtwopoint correlations
..............237
K MTF270,Assignment 1:Reynolds averaged NavierStokes
239
K.1 Twodimensional ow
.........................239
K.2 Analysis
................................239
K.2.1 The momentumequations
...................240
K.2.2 The turbulent kinetic energy equation
............241
K.2.3 The Reynolds stress equations
................241
K.3 Compute derivatives on a curvilinear mesh
..............243
K.3.1 Geometrical quantities
....................244
L MTF270,Assignment 2:LES
245
L.1 Task 2.1
................................245
7
L.2 Task 2.2
................................247
L.3 Task 2.3
................................247
L.4 Task 2.4
................................248
L.5 Task 2.5
................................249
L.6 Task 2.6
................................249
L.7 Task 2.7
................................250
L.8 Task 2.8
................................251
L.9 Task 2.9
................................251
L.10 Task 2.10
................................251
L.11 Task 2.11
................................252
M MTF270:Compute energy spectra fromLES/DNS data using Matlab
253
M.1 Introduction
..............................253
M.2 An example of using FFT
.......................253
M.3 Energy spectrumfromthe twopoint correlation
...........255
M.4 Energy spectra fromthe autocorrelation
................257
N MTF270,Assignment 4:Hybrid LESRANS
258
N.1 Time history
..............................258
N.2 Mean velocity prole
.........................259
N.3 Resolved stresses
...........................259
N.4 Turbulent kinetic energy
........................260
N.5 The modelled turbulent shear stress
..................260
N.6 Turbulent length scales
........................260
N.7 SAS turbulent length scales
......................260
O MTF270,Assignment 5:Embedded LES with PANS
262
O.1 Time history
..............................262
O.2 Resolved stresses
...........................263
O.3 Turbulent viscosity
...........................263
O.4 Modelled stresses
...........................264
O.5 Turbulent SGS dissipation
.......................264
P MTF270:Transformation of a tensor
266
P.1 Rotation to principal directions
....................267
P.2 Transformation of a velocity gradient
.................268
Q MTF270:Green's formulas
269
Q.1 Green's rst formula
..........................269
Q.2 Green's second formula
........................269
Q.3 Green's third formula
.........................269
Q.4 Analytical solution to Poisson's equation
...............272
R MTF270:Learning outcomes for 2013
273
S References
279
8
TME225 Mechanics of uids
L.Davidson
Division of Fluid Dynamics,Department of Applied Mechanics
Chalmers University of Technology,G¨oteborg,Sweden
http://www.tfd.chalmers.se/lada,lada@chalmers.se
This report can be downloaded at
http://www.tfd.chalmers.se/lada/MoF/
1.Motion,ow
9
X
i
T(x
1
i
,t
1
)
T(x
2
i
,t
2
)
T(X
i
,t
1
)
T(X
i
,t
2
)
Figure 1.1:The temperature of a uid particle described in L agrangian,T(X
i
,t),or
Eulerian,T(x
i
,t),approach.
1 Motion,ow
1.1 Eulerian,Lagrangian,material derivative
See also [
2
],Chapt.3.2.
Assume a uid particle is moving along the line in Fig.
1.1
.We can choose to study
its motion in two ways:Lagrangian or Eulerian.
In the Lagrangian approach we keep track of its original position (X
i
) and follow
its path which is described by x
i
(X
i
,t).For example,at time t
1
the temperature of
the particle is T(X
i
,t
1
),and at time t
2
its temperature is T(X
i
,t
2
),see Fig.
1.1
.This
approach is not used for uids because it is very tricky to de ne and follow a uid
particle.It is however used when simulating movement of particles in uids (for ex
ample soot particles in gasolineair mixtures in combustion applications).The speed
of the particle is then expressed as a function of time and its position at time zero,i.e.
v
i
= v
i
(X
i
,t).
In the Eulerian approach we pick a position,e.g.x
1
i
,and watch the particle pass
by.This approach is used for uids.The temperature of the u id,T,for example,is
expressed as a function of the position,i.e.T = T(x
i
),see Fig.
1.1
.It may be that the
temperature at position x
i
,for example,varies in time,t,and then T = T(x
i
,t).
Now we want to express how the temperature of a uid particle v aries.In the
Lagrangian approach we rst pick the particle (this gives it s starting position,X
i
).
Once we have chosen a particle its starting position is xed,and temperature varies
only with time,i.e.T(t) and the temperature gradient can be written dT/dt.
In the Eulerian approach it is a little bit more difcult.We a re looking for the
temperature gradient,dT/dt,but since we are looking at xed points in space we
need to express the temperature as a function of both time and space.From classical
mechanics,we know that the velocity of a uid particle is the time derivative of its
space location,i.e.v
i
= dx
i
/dt.The chainrule nowgives
dT
dt
=
∂T
∂t
+
dx
j
dt
∂T
∂x
j
=
∂T
∂t
+v
j
∂T
∂x
j
(1.1)
Note that we have to use partial derivative on T since it is a function of more than one
(independent) variable.The rst termon the right side is th e local rate of change;by local rate
of changethis we mean that it describes the variation of T in time at position x
i
.The second term
on the right side is called the convective rate of change,which means that it describes
Conv.rate
of change
1.2.Viscousstress,pressure
10
x
1
x
2
σ
11
σ
12
σ
13
(a) Stress components on a surface.
x
1
x
2
f
i
(b) Volume force,f
i
= (0,−g,0),acting in
the middle of the uid element.
Figure 1.2:Stress tensor and volume (gravitation) force.
the variation of T in space when is passes the point x
i
.The left side in Eq.
1.1
is called
the material derivative and is in this text denoted by dT/dt.Material
derivativeEquation
1.1
can be illustrated as follows.Put your nger out in the blowi ng wind.
The temperature gradient you're nger experiences is ∂T/∂t.Imagine that you're a
uid particle and that you ride on a bike.The temperature gra dient you experience is
the material derivative,dT/dt.
Exercise 1 Write out Eq.
1.1
,termbyterm.
1.2 Viscous stress,pressure
See also [
2
],Chapts.6.3 and 8.1.
We have in Part I [
3
] derived the balance equation for linear momentum which
reads
ρ˙v
i
−σ
ji,j
−ρf
i
= 0
(1.2)
Switch notation for the material derivative and derivatives so that
ρ
dv
i
dt
=
∂σ
ji
∂x
j
+ρf
i
(1.3)
where the rst and the second term on the right side represent s,respectively,the net
force due to surface and volume forces (σ
ij
denotes the stress tensor).Stress is force
per unit area.The rst termincludes the viscous stress tens or,τ
ij
.As you have learnt
earlier,the rst index relates to the surface at which the st ress acts and the second
index is related to the stress component.For example,on a surface whose normal is
n
i
= (1,0,0) act the three stress components σ
11
,σ
12
and σ
13
,see Fig.
1.2
a;the
volume force acts in the middle of the uid element,see Fig.
1.2
b.
In the present notation we denote the velocity vector by v = v
i
= (v
1
,v
2
,v
3
)
and the coordinate by x = x
i
= (x
1
,x
2
,x
3
).In the literature,you may nd other
notations of the velocity vector such as u
i
= (u
1
,u
2
,u
3
).If no tensor notation is used
the velocity vector is usually denoted as (u,v,w) and the coordinates as (x,y,z).
1.3.Strainratetensor,vorticity
11
The diagonal components of σ
ij
represent normal stresses and the offdiagonal
components of σ
ij
represent the shear stresses.In Part I [
3
] you learnt that the pressure
is dened as minus the sumof the normal stress,i.e.
p = −σ
kk
/3 (1.4)
The pressure,p,acts as a normal stress.In general,pressure is a thermodynamic prop
erty,p
t
,which can be obtained for example fromthe ideal gas law.I n that case the
thermodynamics pressure,p
t
,and the mechanical pressure,p,may not be the same but
Eq.
1.4
is nevertheless used.The viscous stress tensor,τ
ij
,is obtained by subtracting
the trace,σ
kk
/3 = −p,fromσ
ij
;the stress tensor can then be written as
σ
ij
= −pδ
ij
+τ
ij
(1.5)
τ
ij
is the deviator of σ
ij
.The expression for the viscous stress tensor is found in Eq.
2.4
at p.
21
.The minussign in front of p appears because the pressure acts into the surface.
When there's no movement,the viscous stresses are zero and t hen of course the normal
stresses are the same as the pressure.In general,however,the normal stresses are the
sumof the pressure and the viscous stresses,i.e.
σ
11
= −p +τ
11
,σ
22
= −p +τ
22
,σ
33
= −p +τ
33
,(1.6)
Exercise 2 Consider Fig.
1.2
.Show how σ
21
,σ
22
,σ
23
act on a surface with normal
vector n
i
= (0,1,0).Show also how σ
31
,σ
32
,σ
33
act on a surface with normal vector
n
i
= (0,0,1).
Exercise 3 Write out Eq.
1.5
on matrix form.
1.3 Strain rate tensor,vorticity
See also [
2
],Chapt.3.5.3,3.6.
We need an expression for the viscous stresses,τ
ij
.They will be expressed in the
velocity gradients,
∂v
i
∂x
j
.Hence we will now discuss the velocity gradients.
The velocity gradient tensor can be split into two parts as
∂v
i
∂x
j
=
1
2
∂v
i
∂x
j
+
∂v
i
∂x
j
2∂v
i
/∂x
j
+
∂v
j
∂x
i
−
∂v
j
∂x
i
=0
=
1
2
∂v
i
∂x
j
+
∂v
j
∂x
i
+
1
2
∂v
i
∂x
j
−
∂v
j
∂x
i
= S
ij
+Ω
ij
(1.7)
where
S
ij
is a symmetric tensor called the strainrate tensor Strainrate
tensor
Ω
ij
is a antisymmetric tensor called the vorticity tensor
vorticity ten
sor
The vorticity tensor is related to the familiar vorticity vector which is the curl of
the velocity vector,i.e.ω = ∇×v,or in tensor notation
ω
i
= ǫ
ijk
∂v
k
∂x
j
(1.8)
1.3.Strainratetensor,vorticity
12
If we set,for example,i = 3 we get
ω
3
= ∂v
2
/∂x
1
−∂v
1
/∂x
2
.(1.9)
The vorticity represents rotation of a uid particle.Inser ting Eq.
1.7
into Eq.
1.8
gives
ω
i
= ǫ
ijk
(S
kj
+Ω
kj
) = ǫ
ijk
Ω
kj
(1.10)
since ǫ
ijk
S
kj
= 0 because the product of a symmetric tensor (S
kj
) and an anti
symmetric tensor (ε
ijk
) is zero.Let us show this for i = 1 by writing out the full
equation.Recall that S
ij
= S
ji
(i.e.S
12
= S
21
,S
13
= S
31
,S
23
= S
32
) and
ǫ
ijk
= −ǫ
ikj
= ǫ
jki
etc (i.e.ε
123
= −ε
132
= ε
231
...,ε
113
= ε
221
=...ε
331
= 0)
ε
1jk
S
kj
= ε
111
S
11
+ε
112
S
21
+ε
113
S
31
+ε
121
S
12
+ε
122
S
22
+ε
123
S
32
+ε
131
S
13
+ε
132
S
23
+ε
133
S
33
= 0 S
11
+0 S
21
+0 S
31
+0 S
12
+0 S
22
+1 S
32
+0 S
13
−1 S
23
+0 S
33
= S
32
−S
23
= 0
(1.11)
Now let us invert Eq.
1.10
.We start by multiplying it with ε
iℓm
so that
ε
iℓm
ω
i
= ε
iℓm
ǫ
ijk
Ω
kj
(1.12)
The εδidentity gives (see Table
A.1
at p.
A.1
)
ε
iℓm
ǫ
ijk
Ω
kj
= (δ
ℓj
δ
mk
−δ
ℓk
δ
mj
)Ω
kj
= Ω
mℓ
−Ω
ℓm
= 2Ω
mℓ
(1.13)
This can easily be proved by writing all the components,see Table
A.1
at p.
A.1
.Hence
we get with Eq.
1.8
Ω
mℓ
=
1
2
ε
iℓm
ω
i
=
1
2
ε
ℓmi
ω
i
= −
1
2
ε
mℓi
ω
i
(1.14)
or,switching indices
Ω
ij
= −
1
2
ε
ijk
ω
k
(1.15)
A much easier way to go from Eq.
1.10
to Eq.
1.15
is to write out the components of
Eq.
1.10
.Here we do it for i = 1
ω
1
= ε
123
Ω
32
+ε
132
Ω
23
= Ω
32
−Ω
23
= −2Ω
23
(1.16)
and we get
Ω
23
= −
1
2
ω
1
(1.17)
which indeed is identical to Eq.
1.15
.
Exercise 4 Write out the second and third component of the vorticity vector given in
Eq.
1.8
(i.e.ω
2
and ω
3
).
Exercise 5 Complete the proof of Eq.
1.11
for i = 2 and i = 3.
1.4.Productofasymmetricandantisymmetrictensor
13
Exercise 6 Write out Eq.
1.16
also for i = 2 and i = 3 and nd an expression for Ω
12
and Ω
13
(cf.Eq.
1.17
).Show that you get the same result as in Eq.
1.15
.
Exercise 7 In Eq.
1.17
we proved the relation between Ω
ij
and ω
i
for the offdiagonal
components.What about the diagonal components of Ω
ij
?What do you get from
Eq.
1.7
?
Exercise 8 From you course in linear algebra,you should remember how to compute
a vector product using Sarrus'rule.Use it to compute the vector product
ω = ∇×v =
ˆe
1
ˆe
2
ˆe
3
∂
∂x
1
∂
∂x
2
∂
∂x
3
v
1
v
2
v
3
Verify that this agrees with the expression in tensor notation in Eq.
1.8
.
1.4 Product of a symmetric and antisymmetric tensor
In this section we show the proof that the product of a symmetric and antisymmetric
tensor is zero.First,we have the denitions:
• A tensor a
ij
is symmetric if a
ij
= a
ji
;
• A tensor b
ij
is antisymmetric if b
ij
= −b
ji
.
It follows that for an antisymmetric tensor that all diagonal components must be
zero;for example,b
11
= −b
11
can only be satised if b
11
= 0.
The (inner) product of a symmetric and antisymmetric tensor is always zero.This
can be shown as follows
a
ij
b
ij
= a
ji
b
ij
= −a
ij
b
ji
,
where we rst used the fact that a
ij
= a
ji
(symmetric),and then that b
ij
= −b
ji
(antisymmetric).Since the indices i and j are both dummy indices we can interchange
them,so that
a
ij
b
ij
= −a
ji
b
ij
= −a
ij
b
ij
,
and thus the product must be zero.
This can of course also be shown be writing out a
ij
b
ij
on component form,i.e.
a
ij
b
ij
= a
11
b
11
+a
12
b
12
I
+a
13
b
13
II
+a
21
b
21
I
+a
22
b
22
+a
23
b
23
III
+a
31
b
31
II
+a
32
b
32
III
+a
33
b
33
= 0
The underlined terms are zero;terms I cancel each other as do terms II and III.
1.5.Deformation,rotation
14
x
1
x
2
−
∂v
1
∂x
2
Δx
2
Δt
∂v
2
∂x
1
Δx
1
Δt
Δx
1
Δx
2
α
α
Figure 1.3:Rotation of a uid particle during time Δt.Here ∂v
1
/∂x
2
= −∂v
2
/∂x
1
so that −Ω
12
= ω
3
/2 = ∂v
2
/∂x
1
> 0.
1.5 Deformation,rotation
See also [
2
],Chapt.3.3.
The velocity gradient can,as shown above,be divided into two parts:S
ij
and Ω
ij
.
We have shown that the latter is connected to rotation of a uid particle.During rotation rotation
the uid particle is not deformed.This movement can be illus trated by Fig.
1.3
.The
vertical movement (v
2
) of the lowerright corner (x
1
+Δx
1
) of the particle in Fig.
1.3
is estimated as follows.The velocity at the lowerleft corner is v
2
(x
1
).Now we need
the velocity at the lowerright corner which is located at x
1
+ Δx
1
.It is computed
using the rst termin the Taylor series as
1
v
2
(x
1
+Δx
1
) = v
2
(x
1
) +Δx
1
∂v
2
∂x
1
It is assumed that the uid particle in Fig.
1.3
is rotated the angle α during the time
Δt.The vorticity during this rotation is ω
3
= ∂v
2
/∂x
1
− ∂v
1
/∂x
2
= −2Ω
12
.The
vorticity ω
3
should be interpreted as twice the average rotation of the horizontal edge
(∂v
2
/∂x
1
) and vertical edge (−∂v
1
/∂x
2
).
Next let us have a look at the deformation caused by S
ij
.It can be divided into two
parts,namely shear and elongation (also called extension or dilatation).The deforma
tion due to shear is caused by the offdiagonal terms of S
ij
.In Fig.
1.4
,a pure shear de
formation by S
12
= (∂v
1
/∂x
2
+∂v
2
/∂x
1
)/2 is shown.The deformation due to elon
gation is caused by the diagonal terms of S
ij
.Elongation caused by S
11
= ∂v
1
/∂x
1
is
illustrated in Fig.
1.5
.
1
this corresponds to the equation for a straight line y = kx+ℓ where k is the slope which is equal to the
derivative of y,i.e.dy/dx
1.5.Deformation,rotation
15
x
1
x
2
∂v
1
∂x
2
Δx
2
Δt
∂v
2
∂x
1
Δx
1
Δt
Δx
1
Δx
2
α
α
Figure 1.4:Deformation of a uid particle by shear during ti me Δt.Here ∂v
1
/∂x
2
=
∂v
2
/∂x
1
so that S
12
= ∂v
1
/∂x
2
> 0.
x
1
x
2
∂v
1
∂x
1
Δx
1
Δt
Δx
1
Δx
2
Figure 1.5:Deformation of a uid particle by elongation dur ing time Δt.
In general,a uid particle experiences a combination of rot ation,deformation and
elongation as indeed is given by Eq.
1.7
.
Exercise 9 Consider Fig.
1.3
.Show and formulate the rotation by ω
1
.
1.6.Irrotationalandrotationalow
16
x
1
S
x
2
t
i
dℓ
Figure 1.6:The surface,S,is enclosing by the line ℓ.The vector,t
i
,denotes the unit
tangential vector of the enclosing line,ℓ.
Exercise 10 Consider Fig.
1.4
.Show and formulate the deformation by S
23
.
Exercise 11 Consider Fig.
1.5
.Show and formulate the elongation by S
22
.
1.6 Irrotational and rotational ow
In the previous subsection we introduced different types of movement of a uid parti
cle.One type of movement was rotation,see Fig.
1.3
.Flows are often classied based
on rotation:they are rotational (ω
i
6= 0) or irrotational (ω
i
= 0);the latter type is also
called inviscid ow or potential ow.We'll talk more about t hat later on.In this sub
section we will give examples of one irrotational and one rotational ow.In potential
ow,there exists a potential,Φ,from which the velocity components can be obtained
as
v
k
=
∂Φ
∂x
k
(1.18)
Before we talk about the ideal vortex line in the next section,we need to introduce
the concept circulation.Consider a closed line on a surface in the x
1
−x
2
plane,see
Fig.
1.6
.When the velocity is integrated along this line and projected onto the line we
obtain the circulation
Γ =
I
v
m
t
m
dℓ (1.19)
Using Stokes's theoremwe can relate the circulation to the vorticity as
Γ =
Z
ℓ
v
m
t
m
dℓ =
Z
S
ε
ijk
∂v
k
∂x
j
n
i
dS =
Z
S
ω
3
dS (1.20)
where n
i
= (0,0,1) is the unit normal vector of the surface S.Equation
1.20
reads in
vector notation
Γ =
Z
ℓ
v tdℓ =
Z
S
(∇×v) ndS =
Z
S
ω
3
dS (1.21)
The circulation is useful in aeronautics and windpower engineering where the lift
of an airfoil or a rotorblade is expressed in the circulation for a 2D section.The lift
force is computed as
L = ρV Γ (1.22)
1.6.Irrotationalandrotationalow
17
a
b
Figure 1.7:Ideal vortex.The uid particle (i.e.its diagon al,see Fig.
1.3
) does not
rotate.
where V is the velocity around the airfoil (for a rotorblade it is the relative velocity,
since the rotorblade is rotating).In a recent MSc thesis project,an inviscid simula
tion method (based on the circulation and vorticity sources) was used to compute the
aerodynamic loads for windturbines [
4
].
Exercise 12 In potential ow ω
i
= ε
ijk
∂v
k
/∂x
j
= 0.Multiply Eq.
1.18
by ε
ijk
and
derivate with respect to x
k
(i.e.take the curl of) and show that the right side becomes
zero as it should,i.e.ε
ijk
∂
2
Φ/(∂x
k
∂x
j
) = 0.
1.6.1 Ideal vortex line
The ideal vortex line is an irrotational (potential) ow whe re the uid moves along
circular paths,see Fig.
1.7
.The velocity eld in polar coordinates reads
v
θ
=
Γ
2πr
,v
r
= 0 (1.23)
where Γ is the circulation.Its potential reads
Φ =
Γθ
2π
(1.24)
The velocity,v
θ
,is then obtained as
v
θ
=
1
r
∂Φ
∂θ
=
Γ
2πr
To transform Eq.
1.23
into Cartesian velocity components,consider Fig.
1.8
.The
Cartesian velocity vectors are expressed as
v
1
= −v
θ
sin(θ) = −v
θ
x
2
r
= −v
θ
x
2
(x
2
1
+x
2
2
)
1/2
v
2
= v
θ
cos(θ) = v
θ
x
1
r
= v
θ
x
1
(x
2
1
+x
2
2
)
1/2
(1.25)
Inserting Eq.
1.25
into Eq.
1.23
we get
v
1
= −
Γx
2
2π(x
2
1
+x
2
2
)
,v
2
=
Γx
1
2π(x
2
1
+x
2
2
)
.(1.26)
1.6.Irrotationalandrotationalow
18
x
1
x
2
r
θ
θ
v
θ
Figure 1.8:Transformation of v
θ
into Cartesian components.
To verify that this ow is a potential ow,we need to show that the vorticity,ω
i
=
ε
ijk
∂v
k
/∂x
j
is zero.Since it is a twodimensional ow ( v
3
= ∂/∂x
3
= 0),ω
1
=
ω
2
= 0,we only need to compute ω
3
= ∂v
2
/∂x
1
−∂v
1
/∂x
2
.The velocity derivatives
are obtained as
∂v
1
∂x
2
= −
Γ
2π
x
2
1
−x
2
2
(x
2
1
+x
2
2
)
2
,
∂v
2
∂x
1
=
Γ
2π
x
2
2
−x
2
1
(x
2
1
+x
2
2
)
2
(1.27)
and we get
ω
3
=
Γ
2π
1
(x
2
1
+x
2
2
)
2
(x
2
2
−x
2
1
+x
2
1
−x
2
2
) = 0 (1.28)
which shows that the ow is indeed a potential ow,i.e.irrotational (ω
i
≡ 0).Note
that the deformation is not zero,i.e.
S
12
=
1
2
∂v
1
∂x
2
+
∂v
2
∂x
1
=
Γ
2π
x
2
2
(x
2
1
+x
2
2
)
2
(1.29)
Hence a uid particle in an ideal vortex does deform but it doe s not rotate (i.e.its
diagonal does not rotate,see Fig.
1.7
).
It may be little confusing that the owpath forms a vortex but the owitself has no
vorticity.Thus one must be very careful when using the words vortex a nd vorticity.vortex vs.
vorticityBy vortex we usually mean a recirculation region of the mean ow.That the ow has
no vorticity (i.e.no rotation) means that a uid particle mo ves as illustrated in Fig.
1.7
.
As a uid particle moves fromposition a to b on its counterclockwiserotating path
the particle itself is not rotating.This is true for the who le ow eld,except at the
center where the uid particle does rotate.This is a singula r point as is seen from
Eq.
1.23
for which ω
3
→∞.
Note that generally a vortex has vorticity,see Section
4.2
.The ideal vortex is a very
special ow case.
1.6.2 Shear ow
Another example which is rotational is a shear ow in which
v
1
= cx
2
2
,v
2
= 0 (1.30)
1.7.Eigenvaluesandandeigenvectors:physicalinterpretation
19
a
a
b
c
v
1
v
1
x
1
x
2
x
3
Figure 1.9:A shear ow.The uid particle rotates.v
1
= cx
2
2
.
σ
11
σ
12
σ
21
σ
23
x
1
x
2
x
1
′
x
2
′
α
ˆv
1
ˆv
2
λ
1
λ
2
Figure 1.10:A twodimensional uid element.Left:in origi nal state;right:rotated to
principal coordinate directions.λ
1
and λ
2
denote eigenvalues;ˆv
1
and ˆv
2
denote unit
eigenvectors.
with c,x
2
> 0,see Fig.
1.9
.The vorticity vector for this ow reads
ω
1
= ω
2
= 0,ω
3
=
∂v
2
∂x
1
−
∂v
1
∂x
2
= −2cx
2
(1.31)
When the uid particle is moving from position a,via b to position c it is indeed
rotating.It is rotating in clockwise direction.Note that the positive rotating direction
is dened as the counterclockwise direction,indicated by α in Fig.
1.9
.This is why
the vorticity,ω
3
,is negative (= −2cx
2
).
1.7 Eigenvalues and and eigenvectors:physical interpretation
See also [
2
],Chapt.2.5.5.
Consider a twodimensional uid (or solid) element,see Fig.
1.10
.In the left gure
it is oriented along the x
1
−x
2
coordinate system.On the surfaces act normal stresses
(σ
11
,σ
22
) and shear stresses (σ
12
,σ
21
).The stresses forma tensor,σ
ij
.Any tensor has
eigenvectors and eigenvalues (also called principal vectors and principal values).Since
σ
ij
is symmetric,the eigenvalues are real (i.e.not imaginary).The eigenvalues are
obtained from the characteristic equation,see [
2
],Chapt.2.5.5 or Eq.
13.5
at p.
117
.
When the eigenvalues have been obtained,the eigenvectors can be computed.Given
the eigenvectors,the uid element is rotated α degrees so that its edges are aligned
with the eigenvectors,ˆv
1
= ˆx
1
′
and ˆv
2
= ˆx
2
′
,see right part of Fig.
1.10
.Note that the
1.7.Eigenvaluesandandeigenvectors:physicalinterpretation
20
sign of the eigenvectors is not dened,which means that the e igenvectors can equally
well be chosen as −ˆv
1
and/or −ˆv
2
.In the principal coordinates x
1
′
−x
2
′
(right part
of Fig.
1.10
),there are no shear stresses on the surfaces of the uid elem ent.There
are only normal stresses.This is the very denition of eigen vectors.Furthermore,the
eigenvalues are the normal stresses in the principal coordinates,i.e.λ
1
= σ
1
′
1
′
and
λ
2
= σ
2
′
2
′
.
2.Governingowequations
21
2 Governing ow equations
See also [
2
],Chapts.5 and 8.1.
2.1 The NavierStokes equation
2.1.1 The continuity equation
The rst equation is the continuity equation (the balance eq uation for mass) which
reads [
3
]
˙ρ +ρv
i,i
= 0
(2.1)
Change of notation gives
dρ
dt
+ρ
∂v
i
∂x
i
= 0 (2.2)
For incompressible ow ( ρ = const) we get
∂v
i
∂x
i
= 0 (2.3)
2.1.2 The momentumequation
The next equation is the momentumequation.We have formulated the constitutive law
for Newtonian viscous uids [
3
]
σ
ij
= −pδ
ij
+2S
ij
−
2
3
S
kk
δ
ij
(2.4)
Inserting Eq.
2.4
into the balance equations,Eq.
1.3
,we get
ρ
dv
i
dt
= −
∂p
∂x
i
+
∂τ
ji
∂x
j
+ρf
i
= −
∂p
∂x
i
+
∂
∂x
j
2S
ij
−
2
3
∂v
k
∂x
k
δ
ij
+ρf
i
(2.5)
where denotes the dynamic viscosity.This is the NavierStokes equations (sometimes
the continuity equation is also included in the name Navier Stokes).It is also called
the transport equation for momentum.If the viscosity,,is constant it can be moved
outside the derivative.Furthermore,if the ow is incompre ssible the second term in
the parenthesis on the right side is zero because of the continuity equation.If these two
requirements are satised we can also rewrite the rst term in the parenthesis as
∂
∂x
j
(2S
ij
) =
∂
∂x
j
∂v
i
∂x
j
+
∂v
j
∂x
i
=
∂
2
v
i
∂x
j
∂x
j
(2.6)
because of the continuity equation.Equation
2.5
can now for constant and incom
pressible ow be written
ρ
dv
i
dt
= −
∂p
∂x
i
+
∂
2
v
i
∂x
j
∂x
j
+ρf
i
(2.7)
In inviscid (potential) ow,there are no viscous (friction ) forces.In this case,the
NavierStokes equation reduces to the Euler equations Euler
equations
ρ
dv
i
dt
= −
∂p
∂x
i
+ρf
i
(2.8)
2.2.Theenergyequation
22
Exercise 13 Equation
1.3
states that mass times acceleration is equal to the sum of
forces forces (per unit volume).Write out the momentum equation (without using the
summation rule) for the x
1
direction and show the surface forces and the volume force
on a small,square uid element (see lecture notes of Toll &Ek h [
3
]).Now repeat it for
the x
2
direction.
Exercise 14 Formulate the NavierStokes equation for incompressible ow but non
constant viscosity.
2.2 The energy equation
See also [
2
],Chapts.6.4 and 8.1.
We have in Part I [
3
] derived the energy equation which reads
ρ ˙u −v
i,j
σ
ji
+q
i,i
= ρz
(2.9)
where u denotes internal energy.q
i
denotes the conductive heat ux and z the net
radiative heat source.The latter can also be seen as a vector,z
i,rad
;for simplicity,we
neglect the radiation fromhere on.Change of notation gives
ρ
du
dt
= σ
ji
∂v
i
∂x
j
−
∂q
i
∂x
i
(2.10)
In Part I [
3
] we formulated the constitutive law for the heat ux vector ( Fourier's
law)
q
i
= −k
∂T
∂x
i
(2.11)
Inserting the constitutive laws,Eqs.
2.4
and
2.11
,into Eq.
2.10
gives
ρ
du
dt
= −p
∂v
i
∂x
i
+2S
ij
S
ij
−
2
3
S
kk
S
ii
Φ
+
∂
∂x
i
k
∂T
∂x
i
(2.12)
where we have used S
ij
∂v
i
/∂x
j
= S
ij
(S
ij
+Ω
ij
) = S
ij
S
ij
because the product of a
symmetric tensor,S
ij
,and an antisymmetric tensor,Ω
ij
,is zero.Two of the viscous
terms (denoted by Φ) represent irreversible viscous heating (i.e.transformation of
kinetic energy into thermal energy);these terms are important at highspeed ow
2
(for
example reentry fromouter space) and for highly viscous o ws (lubricants).The rst
termon the right side represents reversible heating and cooling due to compression and
expansion of the uid.Equation
2.12
is the transport equation for (internal) energy,u.
Nowwe assume that the ow is incompressible (i.e.the veloci ty should be smaller
than approximately 1/3 of the speed of sound) for which
du = c
p
dT (2.13)
where c
p
is the heat capacity (see Part I) [
3
] so that Eq.
2.12
gives (c
p
is assumed to be
constant)
ρc
p
dT
dt
= Φ+
∂
∂x
i
k
∂T
∂x
i
(2.14)
2
Highspeed ows relevant for aeronautics will be treated in detail in the course Compressible ow in
the MSc programme.
2.3.Transformationofenergy
23
The dissipation termis simplied to Φ = 2S
ij
S
ij
because S
ii
= ∂v
i
/∂x
i
= 0.If we
furthermore assume that the heat conductivity coefcient i s constant and that the uid
is a gas or a common liquid (i.e.not an lubricant oil),we get
dT
dt
= α
∂
2
T
∂x
i
∂x
i
(2.15)
where α = k/(ρc
p
) is the thermal diffusivity.thermal
diffusivity
Pr =
ν
α
(2.16)
is dened where ν = /ρ is the kinematic viscosity.The physical meaning of the
Prandtl number is the ratio of howwell the uid diffuses mome ntumto the howwell it
diffuses internal energy (i.e.temperature).
The dissipation term,Φ,is neglected in Eq.
2.15
because one of two assumptions
are valid:
1.The uid is a gas with low velocity (lower than 1/3 of the speed of sound);this
assumption was made when we assumed that the uid is incompre ssible
2.The uid is a common liquid (i.e.not an lubricant oil).In l ubricant oils the
viscous heating (i.e.the dissipation,Φ) is large.One example is the oil ow in a
gearbox in a car where the temperature usually is more than 100
o
C higher when
the car is running compared to when it is idle.
Exercise 15 Write out and simplify the dissipation term,Φ,in Eq.
2.12
.The rst term
is positive and the second term is negative;are you sure that Φ > 0?
2.3 Transformation of energy
Now we will derive the equation for the kinetic energy,k = v
i
v
i
/2.Multiply Eq.
1.3
with v
i
ρv
i
dv
i
dt
−v
i
∂σ
ji
∂x
j
−v
i
ρf
i
= 0 (2.17)
Using the product rule backwards (Trick 2,see Eq.
8.4
),the rst term on the left side
can be rewritten
ρv
i
dv
i
dt
=
1
2
ρ
d(v
i
v
i
)
dt
= ρ
dk
dt
(2.18)
(v
i
v
i
/2 = k) so that
ρ
dk
dt
= v
i
∂σ
ji
∂x
j
+ρv
i
f
i
(2.19)
Rewrite the stressvelocity termso that
ρ
dk
dt
=
∂v
i
σ
ji
∂x
j
−σ
ji
∂v
i
∂x
j
+ρv
i
f
i
(2.20)
This is the transport equation for kinetic energy,k.Adding Eq.
2.20
to Eq.
2.10
gives
ρ
d(u +k)
dt
=
∂σ
ji
v
i
∂x
j
−
∂q
i
∂x
i
+ρv
i
f
i
(2.21)
This is an equation for the sum of internal and kinetic energy,u + k.This is the
transport equation for total energy,u +k.
2.4.Leftsideofthetransportequations
24
Let us take a closer look at Eqs.
2.10
,
2.20
and
2.21
.First we separate the term
σ
ji
∂v
i
/∂x
j
in Eqs.
2.10
and
2.20
into work related to the pressure and viscous stresses
respectively (see Eq.
1.5
),i.e.
σ
ji
∂v
i
∂x
j
= −p
∂v
i
∂x
i
a
+τ
ji
∂v
i
∂x
j
b=Φ
(2.22)
The following things should be noted.
• The physical meaning of the atermin Eq.
2.22
which includes the pressure,p
is heating/cooling by compression/expansion.This is a re versible process,i.e.
no loss of energy but only transformation of energy.
• The physical meaning of the bterm in Eq.
2.22
which includes the viscous
stress tensor,τ
ij
is a dissipation,which means that kinetic energy is trans
formed to thermal energy.It is denoted Φ,see Eq.
2.12
,and is called viscous
dissipation.It is always positive and represents irreversible heating.
• The dissipation,Φ,appears as a sink termin the equation for the kinetic energy,k
(Eq.
2.20
) and it appears as a source termin the equation for the internal energy,
u (Eq.
2.10
).The transformation of kinetic energy into internal energy takes
place through this source term.
• Φdoes not appear in the equation for the total energy u+k (Eq.
2.21
);this makes
sense since Φ represents a energy transfer between u and k and does not affect
their sum,u +k.
Dissipation is very important in turbulence where transfer of energy takes place at
several levels.First energy is transferred from the mean o w to the turbulent uctua
tions.The physical process is called production of turbulent kinetic energy.Then we
have transformation of kinetic energy from turbulence kinetic energy to thermal en
ergy;this is turbulence dissipation (or heating).At the same time we have the usual
viscous dissipation fromthe mean owto thermal energy,but this is much smaller than
that fromthe turbulence kinetic energy.For more detail,see section 2.4 in [
5
]
3
.
2.4 Left side of the transport equations
So far,the left side in transport equations have been formulated using the material
derivative,d/dt.Let Ψ denote a transported quantity (i.e.Ψ = v
i
,u,T...);the left
side of the equation for momentum,thermal energy,total energy,temperature etc reads
ρ
dΨ
dt
= ρ
∂Ψ
∂t
+ρv
j
∂Ψ
∂x
j
(2.23)
This is often called the nonconservative form.Using the continuity equation,Eq.
2.2
,non
conser
vative
it can be rewritten as
ρ
dΨ
dt
= ρ
∂Ψ
∂t
+ρv
j
∂Ψ
∂x
j
+Ψ
dρ
dt
+ρ
∂v
j
∂x
j
=0
=
ρ
∂Ψ
∂t
+ρv
j
∂Ψ
∂x
j
+Ψ
∂ρ
∂t
+v
j
∂ρ
∂x
j
+ρ
∂v
j
∂x
j
(2.24)
3
can be downloaded from http://www.tfd.chalmers.se/lada
2.5.Materialparticlevs.controlvolume(ReynoldsTransportTheorem)
25
The two underlined terms will form a time derivative term,and the other three terms
can be collected into a convective term,i.e.
ρ
dΨ
dt
=
∂ρΨ
∂t
+
∂ρv
j
Ψ
∂x
j
(2.25)
Thus,left sided of the temperature equation and the NavierStokes,for example,can
be written in three different ways (by use of the chainrule and the continuity equation)
ρ
dv
i
dt
= ρ
∂v
i
∂t
+ρv
j
∂v
i
∂x
j
=
∂ρv
i
∂t
+
∂ρv
j
v
i
∂x
j
ρ
dT
dt
= ρ
∂T
∂t
+ρv
j
∂T
∂x
j
=
∂ρv
i
∂t
+
∂ρv
j
T
∂x
j
(2.26)
The continuity equation can also be written in three ways (by use of the chainrule)
dρ
dt
+ρ
∂v
i
∂x
i
=
∂ρ
∂t
+v
i
∂ρ
∂x
i
+ρ
∂v
i
∂x
i
=
∂ρ
∂t
+
∂ρv
i
∂x
i
(2.27)
The forms on the right sides of Eqs.
2.26
and
2.27
are called the conservative form.conser
vativeWhen solving transport equations (such as the NavierStokes) numerically using nite
volume methods,the left sides in the transport equation are always written as the ex
pressions on the right side of Eqs.
2.26
and
2.27
;in this way Gauss law can be used
to transformthe equations froma volume integral to a surface integral and thus ensur
ing that the transported quantities are conserved.The results may be inaccurate due
to too coarse a numerical grid,but no mass,momentum,energy etc is lost (provided a
transport equation for the quantity is solved):what comes in goes out.
2.5 Material particle vs.control volume (Reynolds Transport The
orem)
See also lecture notes of Toll &Ekh [
3
] and [
2
],Chapt.5.2.
In Part I [
3
] we initially derived all balance equations (mass,momentumand en
ergy) for a collection of material particles.The conservation of mass,d/dt
R
ρdV = 0,
Newton's second law,d/dt
R
ρv
i
= F
i
etc were derived for a collection of particles in
the volume V
part
,where V
part
is a volume that includes the same uid particles all the
time.This means that the volume,V
part
,must be moving and it may expand or contract
(if the density is nonconstant),otherwise particles would move across its boundaries.
The equations we have looked at so far (the continuity equation
2.3
,the NavierStokes
equation
2.7
,the energy equations
2.12
and
2.20
) are all given for a xed control vol
ume.How come?The answer is the Reynolds transport theorem,which converts the
equations frombeing valid for a moving volume with a collection,V
part
,to being valid
for a xed volume,V.The Reynolds transport theoremreads
d
dt
Z
V
part
ΦdV =
Z
V
dΦ
dt
+Φ
∂v
i
∂x
i
dV
=
Z
V
∂Φ
∂t
+v
i
∂Φ
∂x
i
+Φ
∂v
i
∂x
i
dV =
Z
V
∂Φ
∂t
+
∂v
i
Φ
∂x
i
dV
=
Z
V
∂Φ
∂t
dV +
Z
S
v
i
n
i
ΦdS
(2.28)
2.5.Materialparticlevs.controlvolume(ReynoldsTransportTheorem)
26
where V denotes a xed nondeformable volume in space.The divergen ce theorem
was used to obtain the last line and S denotes the bounding surface of volume V.The
last term on the last line represents the net ow of Φ across the xed nondeformable
volume,V.Φ in the equation above can be ρ (mass),ρv
i
(momentum) or ρu (energy).
This equation applies to any volume at every instant and the restriction to a collection
of a material particles is no longer necessary.Hence,in ui d mechanics the transport
equations (Eqs.
2.2
,
2.5
,
2.10
,...) are valid both for a material collection of particles
as well as for a volume;the latter is usually xed (this is not necessary).
3.ExactsolutionstotheNavierStokesequation:twoexamples
27
V
0
x
1
x
2
Figure 3.1:The plate moves to the right with speed V
0
for t > 0.
0
0.2
0.4
0.6
0.8
1
t
1
t
2
t
3
v
1
/V
0
x
2
Figure 3.2:The v
1
velocity at three different times.t
3
> t
2
> t
1
.
3 Exact solutions to the NavierStokes equation:two
examples
3.1 The Rayleigh problem
Imagine the sudden motion of an innitely long at plate.For time greater than zero
the plate is moving with the speed V
0
,see Fig.
3.1
.
Because the plate is innitely long,there is no x
1
dependency.Hence the ow
depends only on x
2
and t,i.e.v
1
= v
1
(x
2
,t) and p = p(x
2
,t).Furthermore,
∂v
1
/∂x
1
= ∂v
3
/∂x
3
= 0 so that the continuity equation gives ∂v
2
/∂x
2
= 0.At
the lower boundary (x
2
= 0) and at the upper boundary (x
2
→ ∞) the velocity com
ponent v
2
= 0,which means that v
2
= 0 in the entire domain.So,Eq.
2.7
gives (no
body forces,i.e.f
1
= 0) for the v
1
velocity component
ρ
∂v
1
∂t
=
∂
2
v
1
∂x
2
2
(3.1)
We will nd that the diffusion process depends on the kinemat ic viscosity,ν = /ρ,
rather than the dynamic one,.The boundary conditions for Eq.
3.1
are
v
1
(x
2
,t = 0) = 0,v
1
(x
2
= 0,t) = V
0
,v
1
(x
2
→∞,t) = 0 (3.2)
The solution to Eq.
3.1
is shown in Fig.
3.2
.For increasing time (t
3
> t
2
> t
1
),the
moving plate affects the uid further and further away fromt he plate.
It turns out that the solution to Eq.
3.1
is a similarity solution;this means that the similarity
solutionnumber of independent variables is reduced by one,in this case fromtwo (x
2
and t) to
one (η).The similarity variable,η,is related to x
2
and t as
η =
x
2
2
√
νt
(3.3)
3.1.TheRayleighproblem
28
If the solution of Eq.
3.1
depends only on η,it means that the solution for a given uid
will be the same (similar) for many (innite) values of x
2
and t as long as the ratio
x
2
/
√
νt is constant.Now we need to transform the derivatives in Eq.
3.1
from ∂/∂t
and ∂/∂x
2
to d/dη so that it becomes a function of η only.We get
∂v
1
∂t
=
dv
1
dη
∂η
∂t
= −
x
2
t
−3/2
4
√
ν
dv
1
dη
= −
1
2
η
t
dv
1
dη
∂v
1
∂x
2
=
dv
1
dη
∂η
∂x
2
=
1
2
√
νt
dv
1
dη
∂
2
v
1
∂x
2
2
=
∂
∂x
2
∂v
1
∂x
2
=
∂
∂x
2
1
2
√
νt
dv
1
dη
=
1
2
√
νt
∂
∂x
2
dv
1
dη
=
1
4νt
d
2
v
1
dη
2
(3.4)
We introduce a nondimensional velocity
f =
v
1
V
0
(3.5)
Inserting Eqs.
3.4
and
3.5
in Eq.
3.1
gives
d
2
f
dη
2
+2η
df
dη
= 0 (3.6)
We have nowsuccessfully transformed Eq.
3.1
and reduced the number of independent
variables fromtwo to one.Nowlet us nd out if the boundary co nditions,Eq.
3.2
,also
can be transformed in a physically meaningful way;we get
v
1
(x
2
,t = 0) = 0 ⇒f(η →∞) = 0
v
1
(x
2
= 0,t) = V
0
⇒f(η = 0) = 1
v
1
(x
2
→∞,t) = 0 ⇒f(η →∞) = 0
(3.7)
Since we managed to transformboth the equation (Eq.
3.1
) and the boundaryconditions
(Eq.
3.7
) we conclude that the transformation is suitable.
Now let us solve Eq.
3.6
.Integration once gives
df
dη
= C
1
exp(−η
2
) (3.8)
Integration a second time gives
f = C
1
Z
η
0
exp(−η
′2
)dη
′
+C
2
(3.9)
The integral above is the error function
erf(η) ≡
2
√
π
Z
η
0
exp(−η
2
) (3.10)
At the limits,the error function takes the values 0 and 1,i.e.erf(0) = 0 and erf(η →
∞) = 1.Taking into account the boundary conditions,Eq.
3.7
,the nal solution to
Eq.
3.9
is (with C
2
= 1 and C
1
= −2/
√
π)
f(η) = 1 −erf(η) (3.11)
3.1.TheRayleighproblem
29
0
0.2
0.4
0.6
0.8
1
0
0.5
1
1.5
2
2.5
3
f
η
Figure 3.3:The velocity,f = v
1
/V
0
,given by Eq.
3.11
.
6
5
4
3
2
1
0
0
0.5
1
1.5
2
2.5
3
τ
21
/(ρV
0
)
η
Figure 3.4:The shear stress for water (ν = 10
−6
) obtained from Eq.
3.12
at time
t = 100 000.
The solution is presented in Fig.
3.3
.Compare this gure with Fig.
3.2
at p.
27
;all
graphs in that gure collapse into one graph in Fig.
3.3
.To compute the velocity,v
1
,
we pick a time t and insert x
2
and t in Eq.
3.3
.Then f is obtained fromEq.
3.11
and
the velocity,v
1
,is computed from Eq.
3.5
.This is how the graphs in Fig.
3.2
were
obtained.
Fromthe velocity prole we can get the shear stress as
τ
21
=
∂v
1
∂x
2
=
V
0
2
√
νt
df
dη
= −
V
0
√
πνt
exp
−η
2
(3.12)
where we used ν = /ρ.Figure
3.4
presents the shear stress,τ
21
.The solid line is
obtained from Eq.
3.12
and circles are obtained by evaluating the derivative,df/dη,
numerically using central differences (f
j+1
−f
j−1
)/(η
j+1
−η
j−1
).
As can be seen from Fig.
3.4
,the magnitude of the shear stress increases for de
creasing η and it is largest at the wall,τ
w
= −ρV
0
/
√
πt
The vorticity,ω
3
,across the boundarylayer is computed fromits denition (E q.
1.31
)
ω
3
= −
∂v
1
∂x
2
= −
V
0
2
√
νt
df
dη
=
V
0
√
πνt
exp(−η
2
) (3.13)
FromFig.
3.2
at p.
27
it is seen that for large times,the moving plate is felt further
and further out in the ow,i.e.the thickness of the boundary layer,δ,increases.Often
3.2.Flowbetweentwoplates
30
x
1
x
2
V
V
P
1
P
1
P
2
h
Figure 3.5:Flow in a horizontal channel.The inlet part of the channel is shown.
the boundary layer thickness is dened by the position where the local velocity,v
1
(x
2
),
reaches 99%of the freestreamvelocity.In our case,this corresponds to the point where
v
1
= 0.01V
0
.FromFig.
3.3
and Eq.
3.11
we nd that this occurs at
η = 1.8 =
δ
2
√
νt
⇒δ = 3.6
√
νt (3.14)
It can be seen that the boundary layer thickness increases with t
1/2
.Equation
3.14
can
also be used to estimate the diffusion length.After,say,10 minutes the diffusion length diffusion
lengthfor air and water,respectively,are
δ
air
= 10.8cm
δ
water
= 2.8cm
(3.15)
As mentioned in the beginning of this section,note that the diffusion length is deter
mined by the kinematic viscosity,ν = /ρ rather than by dynamic one,.
The diffusion length can also be used to estimate the thickness of a developing
boundary layer,see Section
4.3.1
.
Exercise 16 Consider the graphs in Fig.
3.3
.Create this graph with Matlab.
Exercise 17 Consider the graphs in Fig.
3.2
.Note that no scale is used on the x
2
axis
and that no numbers are given for t
1
,t
2
and t
3
.Create this graph with Matlab for both
air and engine oil.Choose suitable values on t
1
,t
2
and t
3
.
Exercise 18 Repeat the exercise above for the shear stress,τ
21
,see Fig.
3.4
.
3.2 Flow between two plates
Consider steady,incompressible ow in a twodimensional c hannel,see Fig.
3.5
,with
constant physical properties (i.e. = const).
3.2.1 Curved plates
Provided that the walls at the inlet are well curved,the velocity near the walls is larger
than in the center,see Fig.
3.5
.The reason is that the ow (with velocity V ) following
the curved wall must change its direction.The physical agent which accomplish this
is the pressure gradient which forces the ow to follow the wa ll as closely as possible
3.2.Flowbetweentwoplates
31
x
1
x
2
P
1
P
2
V
Figure 3.6:Flow in a channel bend.
(if the wall is not sufciently curved a separation will take place).Hence the pressure
in the center of the channel,P
2
,is higher than the pressure near the wall,P
1
.It is thus
easier (i.e.less opposing pressure) for the uid to enter th e channel near the walls than
in the center.This explains the high velocity near the walls.
The same phenomenon occurs in a channel bend,see Fig.
3.6
.The ow V ap
proaches the bend and the owfeels that it is approaching a be nd through an increased
pressure.The pressure near the outer wall,P
2
,must be higher than that near the inner
wall,P
1
,in order to force the ow to turn.Hence,it is easier for the ow to sneak
along the inner wall where the opposing pressure is smaller than near the outer wall:
the result is a higher velocity near the inner wall than near the outer wall.In a three
dimensional duct or in a pipe,the pressure difference P
2
−P
1
creates secondary ow
downstreamthe bend (i.e.a swirling motion in the x
2
−x
3
plane).
3.2.2 Flat plates
The owin the inlet section (Fig.
3.5
) is two dimensional.Near the inlet the velocity is
largest near the wall and further downstreamthe velocity is retarded near the walls due
to the large viscous shear stresses there.The ow is acceler ated in the center because
the mass ow at each x
1
must be constant because of continuity.The acceleration and
retardation of the owin the inlet region is paid for by a pr essure loss which is rather
high in the inlet region;if a separation occurs because of sharp corners at the inlet,the
pressure loss will be even higher.For large x
1
the ow will be fully developed;the
region until this occurs is called the entrance region,and the entrance length can,for
moderately disturbed inow,be estimated as [
6
]
x
1,e
D
h
= 0.016Re
D
h
≡ 0.016
V D
h
ν
(3.16)
where V denotes the bulk (i.e.the mean) velocity,and D
h
= 4A/S
p
where D
h
,
A and S
p
denote the hydraulic diameter,the crosssectional area and the perimeter,
respectively.For ow between two plates we get D
h
= 2h.
Let us nd the governing equations for the fully developed o w region;in this
region the ow does not change with respect to the streamwise coordinate,x
1
(i.e.
∂v
1
/∂x
1
= ∂v
2
/∂x
1
= 0).Since the ow is twodimensional,it does not depend
on the third coordinate direction,x
3
(i.e.∂/∂x
3
),and the velocity in this direction is
zero,i.e.v
3
= 0.Taking these restrictions into account the continuity equation can be
3.2.Flowbetweentwoplates
32
simplied as (see Eq.
2.3
)
∂v
2
∂x
2
= 0 (3.17)
Integration gives v
2
= C
1
and since v
2
= 0 at the walls,it means that
v
2
= 0 (3.18)
across the entire channel (recall that we are dealing with the part of the channel where
the ow is fully developed;in the inlet section v
2
6= 0,see Fig.
3.5
).
Nowlet us turn our attention to the momentumequation for v
2
.This is the vertical
direction (x
2
is positive upwards,see Fig.
3.5
).The gravity acts in the negative x
2
direction,i.e.f
i
= (0,−ρg,0).The momentumequation can be written (see Eq.
2.7
at p.
21
)
ρ
dv
2
dt
≡ ρv
1
∂v
2
∂x
1
+ρv
2
∂v
2
∂x
2
= −
∂p
∂x
2
+
∂
2
v
2
∂x
2
2
−ρg (3.19)
Since v
2
= 0 we get
∂p
∂x
2
= −ρg (3.20)
Integration gives
p = −ρgx
2
+C
1
(x
1
) (3.21)
where the integration constant C
1
may be a function of x
1
but not of x
2
.If we denote
the pressure at the lower wall (i.e.at x
2
= 0) as P we get
p = −ρgx
2
+P(x
1
) (3.22)
Hence the pressure,p,decreases with vertical height.This agrees with our experience
that the pressure decreases at high altitudes in the atmosphere and increases the deeper
we dive into the sea.Usually the hydrostatic pressure,P,is used in incompressible hydrostatic
pressureow.This pressure is zero when the ow is static,i.e.when the velocity eld is zero.
However,when you want the physical pressure,the ρgx
2
as well as the surrounding
atmospheric pressure must be added.
We can now formulate the momentumequation in the streamwise direction
ρ
dv
1
dt
≡ ρv
1
∂v
1
∂x
1
+ρv
2
∂v
1
∂x
2
= −
dP
dx
1
+
∂
2
v
1
∂x
2
2
(3.23)
where p was replaced by P using Eq.
3.22
.Since v
2
= ∂v
1
/∂x
1
= 0 the left side is
zero so
∂
2
v
1
∂x
2
2
=
dP
dx
1
(3.24)
Since the left side is a function of x
2
and the right side is a function of x
1
,we conclude
that they both are equal to a constant.The velocity,v
1
,is zero at the walls,i.e.
v
1
(0) = v
1
(h) = 0 (3.25)
where h denotes the height of the channel,see Eq.
3.5
.Integrating Eq.
3.24
twice and
using Eq.
3.25
gives
v
1
= −
h
2
dP
dx
1
x
2
1 −
x
2
h
(3.26)
3.2.Flowbetweentwoplates
33
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
v
1
/v
1,max
x
2
/h
Figure 3.7:The velocity prole in fully developed channel ow,Eq.
3.28
.
The minus sign on the right side appears because the pressure gradient is decreasing
for increasing x
1
;the pressure is driving the ow.The negative pressure gradient is
constant (see Eq.
3.24
) and can be written as −dP/dx
1
= ΔP/L.
The velocity takes its maximumin the center,i.e.for x
2
= h/2,and reads
v
1,max
=
h
2
ΔP
L
h
2
1 −
1
2
=
h
2
8
ΔP
L
(3.27)
We often write Eq.
3.26
on the form
v
1
v
1,max
=
4x
2
h
1 −
x
2
h
(3.28)
The mean velocity (often called the bulk velocity) is obtained by integrating Eq.
3.28
across the channel,i.e.
v
1,mean
=
v
1,max
h
Z
h
0
4x
2
1 −
x
2
h
dx
2
=
2
3
v
1,max
(3.29)
The velocity prole is shown in Fig.
3.7
Since we know the velocity prole,we can compute the wall she ar stress.Equa
tion
3.26
gives
τ
w
=
∂v
1
∂x
2
= −
h
2
dP
dx
1
=
h
2
ΔP
L
(3.30)
Actually,this result could have been obtained by simply taking a force balance of a
slice of the ow far downstream.
3.2.3 Force balance
To formulate a force balance in the x
1
direction,we start with Eq.
1.3
which reads for
i = 1
ρ
dv
1
dt
=
∂σ
j1
∂x
j
(3.31)
The left hand side is zero since the owis fully developed.Fo rces act on a volume and
its bounding surface.Hence we integrate Eq.
3.31
over the volume of a slice (length
3.2.Flowbetweentwoplates
34
x
1
x
2
τ
w,U
τ
w,L
P
2
P
1
V
h
L
walls
Figure 3.8:Force balance of the ow between two plates.
L),see Fig.
3.8
0 =
Z
V
∂σ
j1
∂x
j
dV (3.32)
Recall that this is the form on which we originally derived the momentum balance
(Newton's second law) in Part I.[
3
] Now use Gauss divergence theorem
0 =
Z
V
∂σ
j1
∂x
j
dV =
Z
S
σ
j1
n
j
dS (3.33)
The bounding surface consists in our case of four surfaces (lower,upper,left and right)
so that
0 =
Z
S
left
σ
j1
n
j
dS+
Z
S
right
σ
j1
n
j
dS+
Z
S
lower
σ
j1
n
j
dS+
Z
S
upper
σ
j1
n
j
dS (3.34)
The normal vector on the lower,upper,left and right are n
i,lower
= (0,−1,0),n
i,upper
=
(0,1,0),n
i,left
= (−1,0,0),n
i,right
= (1,0,0).Inserting the normal vectors and us
ing Eq.
1.5
give
0 = −
Z
S
left
(−p +τ
11
)dS +
Z
S
right
(−p +τ
11
)dS −
Z
S
lower
τ
21
dS +
Z
Supper
τ
21
dS
(3.35)
τ
11
= 0 because ∂v
1
/∂x
1
= 0 (fully developed ow).The shear stress at the upper and
lower surfaces have opposite sign because τ
w
= (∂v
1
/∂x
2
)
lower
= −(∂v
1
/∂x
2
)
upper
.
Using this and Eq.
3.22
give (the gravitation term on the left and right surface cancels
and P and τ
w
are constants and can thus be taken out in front of the integration)
0 = P
1
Wh −P
2
Wh −2τ
w
LW (3.36)
where W is the width (in x
3
direction) of the two plates (for convenience we set W =
1).With ΔP = P
1
−P
2
we get Eq.
3.30
.
3.2.4 Balance equation for the kinetic energy
In this subsection we will use the equation for kinetic energy,Eq.
2.20
.Let us integrate
this equation in the same way as we did for the force balance.The left side of Eq.
2.20
3.2.Flowbetweentwoplates
35
is zero because we assume that the ow is fully developed;usi ng Eq.
1.5
gives
0 =
∂v
i
σ
ji
∂x
j
−σ
ji
∂v
i
∂x
j
+ρv
i
f
i
=0
= −
∂v
j
p
∂x
j
+
∂v
i
τ
ji
∂x
j
+pδ
ij
∂v
i
∂x
j
−τ
ji
∂v
i
∂x
j
Φ
(3.37)
On the rst line v
i
f
i
= v
1
f
1
+ v
2
f
2
= 0 because v
2
= f
1
= 0.The third term on
the second line pδ
ij
∂v
i
/∂x
j
= p∂v
i
/∂x
i
= 0 because of continuity.The last term
corresponds to the viscous dissipation term,Φ (i.e.loss due to friction),see Eq.
2.22
(termb).Now we integrate the equation over a volume
0 =
Z
V
−
∂pv
j
∂x
j
+
∂τ
ji
v
i
∂x
j
−Φ
dV (3.38)
Gauss divergence theoremon the two rst terms gives
0 =
Z
S
(−pv
j
+τ
ji
v
i
)n
j
dS −
Z
V
ΦdV (3.39)
where S is the surface bounding the volume.The unit normal vector is denoted by n
j
which points out from the volume.For example,on the right surface in Fig.
3.8
it is
n
j
= (1,0,0) and on the lower surface it is n
j
= (0,−1,0).Now we apply Eq.
3.39
to the uid enclosed by the at plates in Fig.
3.8
.The second term is zero on all
four surfaces and the rst term is zero on the lower and upper s urfaces (see Exercises
below).We replace the pressure p with P using Eq.
3.22
so that
Z
S
left
&S
right
(−Pv
1
+ρgx
2
v
1
)n
1
dS = −(P
2
−P
1
)
Z
S
left
&S
right
v
1
n
1
dS
= ΔPv
1,mean
Wh
because ρgx
2
n
1
v
1
on the left and right surfaces cancels;P can be taken out of the
integral as it does not depend on x
2
.Finally we get
ΔP =
1
Whv
1,mean
Z
V
ΦdV (3.40)
Exercise 19 For the fully developed ow,compute the vorticity,ω
i
,using the exact
solution (Eq.
3.28
).
Exercise 20 Showthat the rst and second terms in Eq.
3.39
are zero on the upper and
the lower surfaces in Fig.
3.8
.
Exercise 21 Show that the second term in Eq.
3.39
is zero also on the left and right
surfaces in Fig.
3.8
(assume fully developed ow).
Exercise 22 Using the exact solution,compute the dissipation,Φ,for the fully devel
oped ow.
Exercise 23 From the dissipation,compute the pressure drop.Is it the same as that
obtained from the force balance (if not,nd the error;it sho uld be!).
4.Vorticityequationandpotentialow
36
v
1
(x
2
)
x
1
x
2
(x
1
,x
2
)
τ
12
(x
1
−0.5Δx
1
)
τ
12
(x
1
+0.5Δx
1
)
τ
21
(x
2
−0.5Δx
2
)
τ
21
(x
2
+0.5Δx
2
)
p(x
1
−0.5Δx
1
) p(x
1
+0.5Δx
1
)
Figure 4.1:Surface forces in the x
1
direction acting on a uid particle (assuming τ
11
=
τ
22
= 0).v
1
= cx
2
2
and v
2
= 0.∂τ
12
/∂x
1
= 0,∂τ
21
/∂x
2
> 0.v
1
velocity eld
indicated by dashed vectors.
4 Vorticity equation and potential ow
4.1 Vorticity and rotation
Vorticity,ω
i
,was introduced in Eq.
1.8
at p.
11
.As shown in Fig.
1.3
at p.
14
,vorticity
is connected to rotation of a uid particle.Figure
4.1
shows the surface forces in the
x
1
momentumequation acting on a uid particle in a shear ow.Lo oking at Fig.
4.1
it
is obvious that only the shear stresses are able to rotate the uid particle;the pressure
acts through the center of the uid particle and is thus not ab le to affect rotation of the
uid particle.
Let us have a look at the momentum equations in order to show that the viscous
terms indeed can be formulated with the vorticity vector,ω
i
.In incompressible ow
the viscous terms read (see Eqs.
2.4
,
2.5
and
2.6
)
∂τ
ji
∂x
j
=
∂
2
v
i
∂x
j
∂x
j
(4.1)
The right side can be rewritten using the tensor identity
∂
2
v
i
∂x
j
∂x
j
=
∂
2
v
j
∂x
j
∂x
i
−
∂
2
v
j
∂x
j
∂x
i
−
∂
2
v
i
∂x
j
∂x
j
=
∂
2
v
j
∂x
j
∂x
i
−ε
inm
ε
mjk
∂
2
v
k
∂x
j
∂x
n
(4.2)
Let's verify that
∂
2
v
j
∂x
j
∂x
i
−
∂
2
v
i
∂x
j
∂x
j
= ε
inm
ε
mjk
∂
2
v
k
∂x
j
∂x
n
(4.3)
Use the ε −δidentity (see Table
A.1
at p.38
ε
inm
ε
mjk
∂
2
v
k
∂x
j
∂x
n
= (δ
ij
δ
nk
−δ
ik
δ
nj
)
∂
2
v
k
∂x
j
∂x
n
=
∂
2
v
k
∂x
i
∂x
k
−
∂
2
v
i
∂x
j
∂x
j
(4.4)
4.2.Thevorticitytransportequationinthreedimensions
37
The rst term on the right side is zero because of continuity a nd hence we nd that
Eq.
4.2
can indeed be written as
∂
2
v
i
∂x
j
∂x
j
=
∂
2
v
j
∂x
j
∂x
i
−ε
inm
ε
mjk
∂
2
v
k
∂x
j
∂x
n
(4.5)
At the right side we recognize the vorticity,ω
m
= ε
mjk
∂v
k
/∂x
j
,so that
∂
2
v
i
∂x
j
∂x
j
=
∂
2
v
j
∂x
j
∂x
i
−ε
inm
∂ω
m
∂x
n
(4.6)
where the rst on the right side is zero because of continuity,so that
∂
2
v
i
∂x
j
∂x
j
= −ε
inm
∂ω
m
∂x
n
(4.7)
In vector notation the identity Eq.
4.6
reads
∇
2
v = ∇(∇ v) −∇×∇×v = −∇×ω (4.8)
Using Eq.
4.7
,Eq.
4.1
reads
∂τ
ji
∂x
j
= −ε
inm
∂ω
m
∂x
n
(4.9)
Thus,there is a onetoone relation between the viscous termand vorticity:no viscous
terms means no vorticity and vice versa.An imbalance in shear stresses (left side of
Eq.
4.9
) causes a change in vorticity,i.e.generates vorticity (right side of Eq.
4.9
).
Hence,inviscid ow (i.e.frictionless ow) has no rotatio n.(The exception is when
vorticity is transported into an inviscid region,but also in that case no vorticity is
generated or destroyed:it stays constant,unaffected.) Inviscid ow is often called
irrotational ow (i.e.no rotation) or potential ow.The vorticity is always created at potential
boundaries,see Section
4.3.1
.
The main points that we have learnt in this section are:
1.The viscous terms are responsible for creating vorticity;this means that the vor
ticity can't be created or destroyed in inviscid (friction less) ow
2.The viscous terms in the momentumequations can be expressed in ω
i
;consider
ing Item1 this was to be expected.
Exercise 24 Prove the rst equality of Eq.
4.7
using the εδidentity.
Exercise 25 Write out Eq.
4.9
for i = 1 and verify that it is satised.
4.2 The vorticity transport equation in three dimensions
Up to now we have talked quite a lot about vorticity.We have learnt that physically
it means rotation of a uid particle and that it is only the vis cous terms that can cause
rotation of a uid particle.The terms inviscid,irrotation al and potential owall denote
frictionless ow which is equivalent to zero vorticity.There is a small difference be friction
lesstween the three terms because there may be vorticity in inviscid ow that is convected
into the ow at the inlet(s);but also in this case the vortici ty is not affected once it has
4.2.Thevorticitytransportequationinthreedimensions
38
entered the inviscid ow region.However,mostly no distinc tion is made between the
three terms.
In this section we will derive the transport equation for vorticity in incompressible
ow.As usual we start with the NavierStokes equation,Eq.
2.7
at p.
21
.First,we
rewrite the convective termof the incompressible momentumequation (Eq.
2.7
) as
v
j
∂v
i
∂x
j
= v
j
(S
ij
+Ω
ij
) = v
j
S
ij
−
1
2
ε
ijk
ω
k
(4.10)
where Eq.
1.15
on p.
12
was used.Inserting S
ij
= (∂v
i
/∂x
j
+ ∂v
j
/∂x
i
)/2 and
multiplying by two gives
2v
j
∂v
i
∂x
j
= v
j
∂v
i
∂x
j
+
∂v
j
∂x
i
−ε
ijk
v
j
ω
k
(4.11)
The second termon the right side can be written as
v
j
∂v
j
∂x
i
=
1
2
∂(v
j
v
j
)
∂x
i
=
∂k
∂x
i
(4.12)
where k = v
j
v
j
/2.Equation
4.11
can nowbe written as
v
j
∂v
i
∂x
j
=
∂k
∂x
i
no rotation
−ε
ijk
v
j
ω
k
rotation
(4.13)
The last termon the right side is the vector product of v and ω,i.e.v ×ω.
The trick we have achieved is to split the convective term into one term without
rotation (rst termon the right side of Eq.
4.13
) and one termincludingrotation (second
termon the right side).Inserting Eq.
4.13
into the incompressible momentumequation
(Eq.
2.7
) yields
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Comments 0
Log in to post a comment