LECTURE NOTES ON

INTERMEDIATE FLUID MECHANICS

Mihir Sen

Department of Aerospace and Mechanical Engineering

University of Notre Dame
Notre Dame, IN 46556

December 1996

Contents
1 Mathematical preliminaries

1.1 Vectors . . . . . . . . . . . . . . . . . .
1.1.1 Dot product . . . . . . . . . . . .
1.1.2 Cross product . . . . . . . . . . .
1.1.3 Line integrals . . . . . . . . . . .
1.1.4 Surface integrals . . . . . . . . .
1.1.5 Dierential operators . . . . . . .
1.1.6 Identities . . . . . . . . . . . . .
1.2 Special theorems . . . . . . . . . . . . .
1.2.1 Green's theorem . . . . . . . . .
1.2.2 Gauss's theorem . . . . . . . . .
1.2.3 Green's identities . . . . . . . . .
1.2.4 Stokes's theorem . . . . . . . . .
1.3 Cartesian coordinates . . . . . . . . . .
1.3.1 Gradient of a scalar . . . . . . .
1.3.2 Divergence of a vector . . . . . .
1.3.3 Curl of a vector . . . . . . . . . .
1.3.4 Laplacian . . . . . . . . . . . . .
1.4 Orthogonal curvilinear coordinates . . .
1.4.1 Cylindrical polar coordinates . .
1.4.2 Spherical polar coordinates . . .
1.5 Index notation for Cartesian coordinates
1.6 Cartesian tensors . . . . . . . . . . . . .
1.6.1 Symmetric tensors . . . . . . . .
1.6.2 Anti-symmetric tensors . . . . .
1.7 Complex numbers and functions . . . .
1.7.1 Analytic functions . . . . . . . .
1.7.2 Polar form . . . . . . . . . . . .
1.7.3 Integrals . . . . . . . . . . . . . .
Problems . . . . . . . . . . . . . . . . . . . .

2 Flow, rate of deformation and stress

2.1 Flow lines . . . . .
2.1.1 Streamlines
2.1.2 Pathlines .
2.1.3 Streaklines

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.
1

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

7
7
8
8
9
9
10
10
10
11
11
11
11
11
12
12
13
13
14
14
16
18
21
22
23
23
24
25
26

31

31
31
31
31

2.1.4 Timeline . . . . . . . . . . . . .
2.2 Rate of deformation . . . . . . . . . .
2.3 Vortex lines and tubes . . . . . . . . .
2.4 Eulerian and Lagrangian descriptions .
2.4.1 Material derivative . . . . . . .
2.5 Stress . . . . . . . . . . . . . . . . . .
Problems . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

3 Dierential equations of motion

3.1 Approaches to balance equations . . . . . . . . . .

3.1.1 Dierential element . . . . . . . . . . . . .
3.1.2 Integral element . . . . . . . . . . . . . . .
3.2 Mass balance . . . . . . . . . . . . . . . . . . . . .
3.2.1 Elemental control volume . . . . . . . . . .
3.2.2 Elemental control mass . . . . . . . . . . .
3.2.3 Integral control volume . . . . . . . . . . .
3.2.4 Integral control mass . . . . . . . . . . . . .
3.2.5 Equivalence . . . . . . . . . . . . . . . . . .
3.3 Linear momentum equation . . . . . . . . . . . . .
3.4 Linear momentum equation in a non-inertial frame
3.5 Angular momentum equation . . . . . . . . . . . .
3.6 Energy equation . . . . . . . . . . . . . . . . . . .
3.7 Entropy equation . . . . . . . . . . . . . . . . . . .
3.8 Additional relations . . . . . . . . . . . . . . . . .
3.8.1 Equations of state . . . . . . . . . . . . . .
3.8.2 Constitutive equations . . . . . . . . . . . .
3.9 Governing equations for special cases . . . . . . . .
3.9.1 Mass balance . . . . . . . . . . . . . . . . .
3.9.2 Linear momentum equation . . . . . . . . .
3.9.3 Energy equation . . . . . . . . . . . . . . .
3.9.4 Boundary conditions . . . . . . . . . . . . .
3.10 Nondimensionalization . . . . . . . . . . . . . . . .
3.11 Molecular approach . . . . . . . . . . . . . . . . . .
Problems . . . . . . . . . . . . . . . . . . . . . . . . . .

4 Special theorems
4.1
4.2
4.3
4.4
4.5

Circulation . . . . . . .
Kelvin's theorem . . . .
Vorticity equation . . .
Helmholtz's theorem . .
Bernoulli's theorem . . .
4.5.1 Steady ow . . .
4.5.2 Irrotational ow
Problems . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
2

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

31
33
34
34
34
35
36

39

39
39
40
40
41
42
42
42
43
43
45
45
46
48
49
49
49
50
50
51
52
52
53
54
54

59

59
59
60
61
62
63
63
63

5 Ideal ow
5.1
5.2
5.3
5.4
5.5
5.6

Two-dimensional ows . . . . . . . . .
Properties . . . . . . . . . . . . . . . .
Complex representation . . . . . . . .
Polar form . . . . . . . . . . . . . . . .
Summary of equations . . . . . . . . .
Simple ows . . . . . . . . . . . . . . .
5.6.1 Uniform ow . . . . . . . . . .
5.6.2 Source or sink . . . . . . . . . .
5.6.3 Vortex . . . . . . . . . . . . . .
5.6.4 Sector with angle =n . . . . .
5.7 Combined ows . . . . . . . . . . . . .
5.7.1 Doublet . . . . . . . . . . . . .
5.7.2 Cylinder without circulation . .
5.7.3 Cylinder with circulation . . .
5.8 Forces on a submerged body . . . . . .
5.8.1 Cylinder with circulation . . .
5.9 Conformal transformation . . . . . . .
5.9.1 Joukowski transformation . . .
5.10 Three-dimensional axisymmetric ow .
5.10.1 Cylindrical coordinates . . . .
5.10.2 Spherical coordinates . . . . . .
Problems . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

6.1 Flow between at plates . . . . . . . . . . . . . . . .

6.1.1 Alternative derivation of governing equations
6.1.2 Velocity and temperature proles . . . . . . .
6.1.3 Couette ow . . . . . . . . . . . . . . . . . .
6.1.4 Poiseuille ow . . . . . . . . . . . . . . . . . .
6.1.5 Heat generation . . . . . . . . . . . . . . . . .
6.2 Flow between coaxial rotating cylinders . . . . . . .
6.3 Flow in a circular pipe . . . . . . . . . . . . . . . . .
6.3.1 Isothermal wall . . . . . . . . . . . . . . . . .
6.4 Flow over a porous wall . . . . . . . . . . . . . . . .
6.5 Natural convection between vertical at plates . . .
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

6 Incompressible viscous ow: exact solutions

7 Incompressible viscous ow: negligible inertia

7.1 Stokes's ow . . . . . . . .
7.2 Uniform ow past a sphere
7.3 Lubrication theory . . . . .
7.4 Flow in porous media . . .
Problems . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

65

65
66
67
67
70
70
70
70
71
71
71
71
72
72
72
73
73
73
75
76
76
78

79

79
79
82
83
84
84
85
87
88
88
89
90

93

93
93
94
96
96

8 Incompressible viscous ow: at plate boundary layer

8.1 Prandtl boundary layer equations .
8.2 Blasius solution . . . . . . . . . . .
8.3 Momentum integral method . . . .
8.4 Vertical plate natural convection .
Problems . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

9 Hydrodynamic stability and turbulence

99

99
101
103
104
105

107

9.1 Hydrodynamic stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

9.2 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

10 Compressible ow in gases

10.1 Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10.1.1 Equations of state for real gases . . . . . . . . . . . . . .
10.1.2 Thermodynamics of perfect gases . . . . . . . . . . . . . .
10.2 One-dimensional steady ow equations . . . . . . . . . . . . . . .
10.2.1 Dierential form . . . . . . . . . . . . . . . . . . . . . . .
10.2.2 Constant area without friction, heat or work transfer . . .
10.3 Speed of sound . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.4 Stagnation and critical properties . . . . . . . . . . . . . . . . . .
10.5 Normal shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.6 Oblique shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.7 Flow in ducts . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.7.1 Adiabatic, frictionless ow in variable area ducts . . . . .
10.7.2 Frictionless ow with heat transfer in constant area ducts
10.7.3 Adiabatic ow with friction in constant area ducts . . . .
10.8 Multi-dimensional ow . . . . . . . . . . . . . . . . . . . . . . . .
10.8.1 Stagnation enthalpy . . . . . . . . . . . . . . . . . . . . .
10.8.2 Crocco's theorem . . . . . . . . . . . . . . . . . . . . . . .
10.8.3 Irrotational ow . . . . . . . . . . . . . . . . . . . . . . .
10.8.4 Small-perturbation theory . . . . . . . . . . . . . . . . . .
10.8.5 Subsonic ow . . . . . . . . . . . . . . . . . . . . . . . . .
10.8.6 Supersonic ow . . . . . . . . . . . . . . . . . . . . . . . .
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11 Compressible eects in liquids and two-phase ow

11.1 Waterhammer . . . . . . . . . . . . . . . . .
11.2 Two-phase systems . . . . . . . . . . . . . .
11.2.1 Homogeneous two-phase mixture . .
11.2.2 Two-phase ows in horizontal pipes
Problems . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

12.1 Ordinary dierential equations . . . . . . . .

12.1.1 Fourth-order Runge-Kutta integration
12.2 Partial dierential equations . . . . . . . . . .
12.2.1 Primitive variables . . . . . . . . . . .
12.2.2 Stream function-vorticity . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

12 Numerical methods

111

111
111
112
113
114
114
115
117
117
119
120
121
123
126
129
129
130
131
131
132
133
135

137

137
138
138
140
141

143

143
143
144
144
144

12.3 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

12.3.1 Finite dierences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

A Governing equations
A.1
A.2
A.3
A.4
A.5

Integral form . . . . . .
Dierential form . . . .
Cartesian coordinates .
Cylindrical coordinates .
Spherical coordinates . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

149
149
150
151
151
152

B Use of MATLAB

153

C Use of MAPLE
Bibliography

157
159

B.1 Graphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

B.2 Plotting streamlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Chapter 1

Mathematical preliminaries
1.1 Vectors
This is a short introduction to the algebra and calculus of scalar, vector and tensor elds. The
scalar, vector or tensor quantity is considered to be a function of position r in three-dimensional
space.
A scalar eld is one in which a scalar (r) varies with position. In most areas of mechanics,
including uid mechanics, it is also common to talk of vectors which are quantities which have both
magnitude as well as direction. Examples are velocities, accelerations and forces. Here we will
represent vectors using bold-faced letters such as a. Other commonly used notations are a~, a, and
!a . A unit vector n is sometimes indicated as n^ . If a vector f is a function of where it is evaluated,
i.e. f(r) where r is the position vector of the location, f is referred to as a vector eld.
Vectors in three-dimensional physical space have three components. The Cartesian components
a1 , a2 and a3 of the vector a can be indicated explicitly in one of various ways:
2

a1
4 a2 5
a3
[a1 a2 a3 ]T
(a1 ; a2 ; a3 )
a1 e1 + a2 e2 + a3 e3

=
=
=
=

Unit vectors in the x1 , x2 , and x3 directions are e1 , e2 , and e3 . It is also common to use x; y; z
for x1 ; x2 ; x3 and i,j,k for e1 ,e2 ,e3 .

1.1.1 Dot product

The scalar (also dot or inner) product of a and b is written as a
b. and is dened by

a
b

= a1 b1 + a2 b2 + a3 b3
=

3
X

i=1

ai bi

(1.1)
(1.2)

 a



 a


O
x2

x1
a3

a1

x3

Figure 1.1: Vector a with its Cartesian components a1 ; a2 ; a3 .

The dot product is also equivalent to

a
b = ab cos 
(1.3)
where a and b are the magnitudes of a and b respectively, and  is the angle between them. The
vectors a and b are orthogonal if a
b = 0.
1.1.2 Cross product

The cross (or vector) product c of a and b is written as c = ab. The magnitude of c is ab sin ; it
is normal to both a and b in the direction given by the right-hand rule. The cross product is also
equivalent to

ab=

e e e
1

a1 a 2 a3
b1 b2 b3

(1.4)

1.1.3 Line integrals

A line integral is of the form

f
dl
(1.5)
C
where f is a vector eld, and dl is an element of curve C .
In general the value of a line integral depends on the path. If, however, we take f = r, where
 is a scalar eld (see below for the gradient operator), then the integral I is independent of path.
f is then called a conservative eld, and  is its potential.
I=

b


  

-a

Figure 1.2: Arbitrary vectors a and b

1.1.4 Surface integrals

A surface integral is of the form

f
n dA
(1.6)
A
where f is a vector eld, A is an open or closed surface, dA is an element of this surface, and n is a
unit vector normal to this element.
I=

1.1.5 Dierential operators

Surface integrals can be used for coordinate-independent denitions of the gradient, divergence and
curl operators. Thus
Z
1
r = Vlim
n dA
(1.7)
!0 V A
Z
1
n
f dA
(1.8)
r
f = Vlim
!0 V A
Z
1
r  f = Vlim
n  f dA
(1.9)
!0 V A
where (r) is a scalar eld, and f (r) is a vector eld. V is the region enclosed within a closed surface
S , and n is the unit normal to an element of the surface dA. The operators r, r
and r are to
be read as gradient of, divergence of, and curl of, respectively. The notation grad, div, and curl is
also used for these operators.
A related operator is the Laplacian, r2 , where
r2 = r
r
(1.10)
For a vector, the identity (1.21) below is often used. The Laplacian should not be confused with the
operator rr
which is also possible for a vector, but is a dierent operator.
The directional derivative of (r) along a unit vector a is
@
= r
a
(1.11)
@a
where a is the coordinate in that direction. In the special case of a being any tangent to the surface
 = constant, we have @=@a = 0. Thus r is orthogonal to a and must be normal to the surface.
9

1.1.6 Identities

a  (b  c)
(a  b)
(c  d)

r  (r)
r
(r  f )
r
(f )
r  (f )
r
(f  g)
r  (f  g)
r(f
g)
r  (r  f )
d
(f
g)
dt

=
=
=
=
=
=
=
=
=
=

b(a
c) c(a
b)
(a
c)(b
d) (a
d)(b
c)
0

0
r
f + r
f
r  f + r  f
g
(r  f ) f
(r  g)
(g
r)f (f
r)g + f (r
g) g(r
f )
(f
r)g + (g
r)f + f  (r  g) + g  (r  f )
r(r
f ) r2 f
dg df
= f
+
g
dt dt
d
dg df
(f  g) = f  +  g
dt
dt dt

(1.12)
(1.13)
(1.14)
(1.15)
(1.16)
(1.17)
(1.18)
(1.19)
(1.20)
(1.21)
(1.22)
(1.23)

Example 1.1

Prove by expanding into components that

(a  b)
(c  d) = (a
c)(b
d) (a
d)(b
c)

 

 e1 e2 e3   e1 e2 e3 

 
LHS =  a1 a2 a3 
 c1 c2 c3 
b1

b2

b3

d1 d2 d3
c3 d2 ) + (a3 b1 a1 b3 )(c3 d1 c1 d3 ) + (a1 b2 a2 b1 )(c1 d2 c2 d1 )
a2 b3 c3 d2 a3 b2 c2 d3 + a3 b1 c3 d1 + a1 b3 c1 d3
a3 b1 c1 d3 a1 b3 c3 d1 + a1 b2 c1 d2 + a2 b1 c2 d1 a1 b2 c2 d1 a2 b1 c1 d2
(a1 c1 + a2 c2 + a3 c3 )(b1 d1 + b2 d2 + b3 d3 ) (a1 d1 + a2 d2 + a3 d3 )(b1 c1 + b2 c2 + b3 c3 )
(a1 c1 b1 d1 + a1 c1 b2 d2 + a1 c1 b3 d3 + a2 c2 b1 d1 + a2 c2 b2 d2 + a2 c2 b3 d3
+a3 c3 b1 d1 + a3 c3 b2 d2 + a3 c3 b3 d3 )
(a1 d1 b1 c1 + a1 d1 b2 c2 + a1 d1 b3 c3 + a2 d2 b1 c1 + a2 d2 b2 c2 + a2 d2 b3 c3
+a3 d3 b1 c1 + a3 d3 b2 c2 + a3 d3 b3 c3 )

= (a2 b3 a3 b2 )(c2 d3
= a2 b3 c2 d3 + a3 b2 c3 d2

RHS =
=

= LHS

1.2 Special theorems

1.2.1 Green's theorem

Let f = fx (x; y) i + fy (x; y) j be a vector eld, C a closed curve, and A the region enclosed by C , all
in the x-y plane. Then
I
ZZ
y @fx
) dx dy
(1.24)
f
dl = ( @f
@x
@y
C
A
10

1.2.2 Gauss's theorem

Let S be a closed surface, and V the region enclosed within it, then
Z

f
n dA = r
f dV

(1.25)

where dV an element of volume, dA is an element of the surface, and n is the outward unit normal
to it.

1.2.3 Green's identities

Applying Gauss's theorem to the vector f = r , we get

Z

r

n dA =

r
(r ) dV

From this we get Green's rst identity

Z

A
Interchanging  and

n dA =

r

(r2 + r
r ) dV

(1.26)

in the above and subtracting, we get Green's second identity

Z

r)
n dA =

(r

(r2

r ) dV
2

(1.27)

1.2.4 Stokes's theorem

Let S be an open surface, and the curve C its boundary. Then

Z

(r  f )
n d A =

f
dl

(1.28)

where n is the unit vector normal to the element dA, and dl an element of curve C .

1.3 Cartesian coordinates

Consider the element of volume in Cartesian coordinates shown in Figure 1.3. The dierential
operations in this coordinate system can be deduced from the denitions. The r operator is seen
to have the Cartesian denition
e1 @x@ 1 + e2 @x@ 2 + e3 @x@ 3
(1.29)

1.3.1 Gradient of a scalar

We take the reference value of  to be at the center of the element P . At the center of the two faces
which are a distance  dx1 =2 away from P in the x1 -direction, it is (  @=@x1 )dx1 =2. Writing
11

x2
dx1
P

x1

dx2
x3
dx3

Figure 1.3: Element of volume in Cartesian coordinates

V = dx1 dx2 dx3 , equation (1.7) gives

r

@ dx1
1
( +
)e dx dx
= lim
V !0 V
@x1 2 1 2 3

(

+ similar terms from the x2 and x3 faces

@ dx1
)e dx dx
@x1 2 1 2 3

@
e + @ e + @ e
@x1 1 @x2 2 @x3 3

(1.30)

1.3.2 Divergence of a vector

Equation (1.8) becomes

r
f

1
@f dx
= lim
(f1 + 1 1 ) dx2 dx3
V !0 V
@x1 2

(f1

@f1 dx1
) dx2 dx3
@x1 2


+ similar terms from the x2 and x3 faces

=

@f1 @f2 @f3

+
+
@x1 @x2 @x3

(1.31)

1.3.3 Curl of a vector

From equation (1.9), we have

rf

1
@f dx
= lim
(f2 + 2 1 )e3 dx2 dx3
V !0 V
@x1 2
12

(f3 +

@f3 dx1
)e dx dx
@x1 2 2 2 3

(f2

@f2 dx1
)e dx dx + (f3
@x1 2 3 2 3

@f3 dx1
)e dx dx
@x1 2 2 2 3

+ similar terms from the x2 and x3 faces

=


 1
 @

 @x1
 f1

@
@x2

f2

@
@x3

f3

(1.32)

1.3.4 Laplacian
The Laplacian of a scalar is

r  = @@x + @@x + @@x

(1.33)

and that of a vector

@f
r f = @x

(1.34)

2
1

2
2

2
1

2
3

@2f @2f
+
@x22 @x23

is similar.

1.4 Orthogonal curvilinear coordinates

Let the orthogonal curvilinear coordinates of a point be (q1 ; q2 ; q3 ), where qi = qi (x1 ; x2 ; x3 ), the xi s
being the Cartesian coordinates of the point. We can dene scale factors h1 ; h2 ; h3 such that

hi =

s


@x1
@qi

2

@x2
@qi

2

@x3
@qi

2

(1.35)

The dierential operators can be written using these scale factors.

Let (r) be a scalar eld and f(r) a vector eld, both being functions of the position vector r. It
turns out that:
1 @
e + 1 @ e + 1 @ e
h1 @q1 "q1 h2 @q2 q2 h3 @q3 q3
#
@
@
@
1
r
f = h h h @q (f1 h2 h3 ) + @q (f2 h3 h1 ) + @q (f3 h1 h2 )
1 2 3
1
2
3




1  h1@eq1 h2@eq2 h3@eq3 
r  f = h h h  @q1 @q2 @q3 
1 2 3 
f1 h1 f2 h2 f3 h3 
"
!
!
!#
@ h2 h3 @
@ h3 h1 @
@ h1 h2 @
1
2
r  = h h h @q h @q + @q h @q + @q h @q
1 2 3
1
1
1
2
2
2
3
3
3

r

(1.36)
(1.37)
(1.38)
(1.39)

where f1 ; f2 ; f3 are the components of f in the q1 ; q2 ; q3 directions respectively, and eq1 ; eq2 ; eq3 are
the unit vectors in these directions.
13

x2




R( P

O
(( ((
x

x3

Figure 1.4: Cylindrical polar coordinates (R; ; z ) of point P .

1.4.1 Cylindrical polar coordinates

Coordinates:

q1 = R
q2 = 
q3 = z
Relation to Cartesian coordinates:

x1 = R cos 
x2 = R sin 
x3 = z
Scale factors:

h1 = 1
h2 = R
h3 = 1

1.4.2 Spherical polar coordinates

Coordinates:

q1 = r
q2 = 
q3 = '
14

x2



r ( P

O(
(!((
x


x3

Figure 1.5: Spherical polar coordinates (r; ; ') of point P .

Relation to Cartesian coordinates:

x1 = r sin  cos '

x2 = r sin  sin '
x3 = r cos 
Scale factors:

h1 = 1
h2 = r
h3 = r sin 
Example 1.2

Find expressions for the gradient, divergence, and curl in cylindrical coordinates (r; ; z) where
x1 = r cos 
x2 = r sin 
x3 = z
The 1,2 and 3 directions are associated with r, , and z, respectively. From equation (1.35) the scale factors are
r

1 2
2 2
3 2
( @x
) + ( @x
) + ( @x
)
@r
@r
@r

hr

h

= cos2  + sin2 
= 1r
1 2
2 2
3 2
= ( @x
) + ( @x
) + ( @x
)
@
@
@

15

hz

so that

= r2 sin2  + r2 cos2 
= rr
1 2
2 2
3 2
= ( @x
) + ( @x
) + ( @x
)
@z
@z
@z
= 1

r = @
e + 1 @ e + @ e
@r r r @  @z z
@
r
f = 1r [ @r
(fr r) + @@ (f ) + @z@ (fz r)]


 er
1
r  f = r  @r@


fr

re
@
@

f r




@z 
fz

e@z

1.5 Index notation for Cartesian coordinates

Some of the vector relations can be written in a compact form by using the index notation. Let
x1 ; x2 ; x3 represent the three coordinate directions and e1 ; e2 ; e3 the unit vectors in those directions.
Then a vector a may be written as

a=

3
X

i=1

ai ei

(1.40)

where a1 , a2 , and a3 are the three Cartesian components of a. According to the Einstein summation
convention, repetition of indices indicates summation, so that the  symbol can be left out. One
has to take care that an index does not appear more than twice in a given product. Thus we can
simply write
a = ai ei
(1.41)
It is also common to write a = ai , the single free index on the right side indicating that an ei is
assumed.
Two additional symbols are needed for later use. They are the Kronecker delta

ij =

0 if i 6= j
1 if i = j

(1.42)

and the substitution symbol (or Levi-Civita density)

ijk =

8
<

if indices are in cyclical order 1,2,3,1,2,

if indices are not in cyclical order
if two or more indices are the same

(1.43)

The identity

ijk lmn = il jm kn + im jn kl + in jl km
16

il jn km

im jl kn in jm kl

(1.44)

relates the two. The following identities are also easily shown:

ii
ij
ij jk
ijk lmk
ijk ljk
ijk ijk

=
=
=
=
=
=

3
ji
ik
il jm
2il
6

im jl

(1.45)
(1.46)
(1.47)
(1.48)
(1.49)
(1.50)

In this notation the scalar and vectors products can be written as

and
Furthermore

a
b = aibi

(1.51)

a  b = ijk aibj ek

(1.52)

@
e
@xi i
@fi
r
f = @x
i
@fj
ek
r  f = ijk @x
i
2

r2  = @x@ @x
i i
2
@
fj
r2 fj = @x @x
i i
Gauss's theorem, equation (1.25), can be written as

r

(1.53)
(1.54)
(1.55)
(1.56)
(1.57)

@fi
dV
(1.58)
A
V @xi
The index notation sometimes can greatly simplify the demonstration of vector relations, as in
the example below.
fi ni dA =

Example 1.3

Prove using the index notation that

(a  b)
(c  d) = (a
c)(b
d) (a
d)(b
c)
LHS = (ijk ai bj )(lmk cl dm )
= (il jm im jl )ai bj cl dm
= al bm cl dm am bl cldm
= RHS

17

x02 x2

EE

EE

a


EE 
EE 
(( x0
E( (((((( x

1
1

x3

x03

Figure 1.6: Cartesian coordinate systems S and S 0

1.6 Cartesian tensors

Consider the two right-handed Cartesian coordinate systems, S with coordinates x1 ; x2 ; x3 and S 0
with x01 ; x02 ; x03 , shown in Figure 1.6. A scalar  at the common origin O has the same value for either
of the two systems. But now let us look at a vector a represented by the arrow. It has components
a1 ; a2 ; a3 in S and a01 ; a02 ; a03 in S 0 . It is easy to show that the components are related by
aj = Aij a0
(1.59)
i

where Aij = cos(x0i ; xj ) is a transformation array. We may adopt the denition that a quantity with
components that are transformed in this manner is a vector. The inverse transformation is
The transformation array Aij satises

a0j = Aji ai

(1.60)

Aki Akj = ij
(1.61)
This may shown from xi = Aki x0k and x0k = Akj xj , dierentiating the rst with respect to xj and
then substituting the second. The matrix with elements Aij is orthogonal.
A quantity T with components Tij is called a tensor of second order if the components transform
under the relations
Tij = Aki Alj Tkl0
(1.62)
0
Tij = Aik Ajl Tkl
(1.63)
Matrices are often used to represent vectors and tensors; a vector may be represented by a 3  1
column vector, and a tensor by a 3  3 matrix. For example the Kronecker delta, ij , is the identity
matrix I. If a matrix A is Aij , its transpose AT is Aji .
18

The tensor product ab of two vectors a and b is the tensor ai bj . The multiplication Aa of a
3  3 matrix A and a 3  1 column vector a is a vector that is represented in index notation as
Aij aj . Similerly the product AB of two 3  3 matrices A and B is Aij Bjk .
Equation (1.61) is equivalent to
AT A = I
(1.64)
which means that A is an orthogonal matrix with

A
where A

is the inverse of A. Thus

=A

(1.65)

AAT = I

(1.66)

which can also be written as

Aik Ajk = ij
The transformation equations can also be written as
a = AT a0
a0 = Aa
T = AT T0 A
T0 = ATAT

(1.67)
(1.68)
(1.69)
(1.70)
(1.71)

Example 1.4
Consider two Cartesian coordinate systems: S with unit vectors (i; j; k), and S 0 with (i0 ; j0 ; k0 ), where i0 = i,
j0 = (j k)=p2, k0 = (j + k)=p2. The tensor T has the following components in S :

1 0 0
0 1 0
0 0 3

Find its components in S 0 .

The transformation matrix is
A=
From equation (1.71) we get,
T0 =

1 0p
0
0 1=p2 1=pp2
0 1 = 2 1= 2
1 0 0
0 1 2
0 2 1

Tensors of higher order may be similarly dened. The components of a tensor of order three,
Sijk , transform as
0
Sijk = Ali Amj Ank Slmn
(1.72)
and so on.
In summary, the order of a tensor is dened by the transformation rule of its components if the
coordinate system is rotated.
19

Transformation

=
a0i = Aij aj
Tij0 = Aik Ajl Tkl
0
Sijk = Ail Ajm Akn Slmn

Quantity

 is a tensor of order zero, or scalar

ai are components of a tensor of order one, or vector)
Tij are components of a tensor of order two
Sijk are components of a tensor of order three

If there is no risk of confusion, we will refer to a tensor of order two as simply a tensor. We may
also say that fi is a vector and Tij is a tensor, instead of being components of a vector and a tensor.
It can be shown that a tensor of order two T is an operator which transforms a vector a into
another b. Thus
bi = Tij aj
(1.73)
where ai and bi are the components of a and b respectively.
Example 1.5
If Tij

and aj are the components of a second order tensor and vector, respectively, show that bi = Tij aj
are the components 0of a vector.
In the system S , the components are
b0i = Tij0 a0j



= Aik Ajl Tkl Ajm fm
= Aik lm Tkl am
= Aik Tkl al
= Aik bk
proving that gi are the components of a vector.

Similarly it can be shown that Tij = ai bj are the components of a tensor of order two if ai and
bj are the components of vectors. This
product of vectors a and b is called a dyad and written as a
Pk
b. A linear combination of dyads, n=1 Cn an bn is a dyadic.
Notice that ai bi is a scalar while ai bj is a tensor. The process of equating indices is called a
contraction, also indicated by a dot, so that the contraction of ab is a
b. For two second order
tensors, a single contraction R
S = Rij Sjk will give a tensor and a double contraction R : S = Rij Sji
a scalar.
The gradient operator increases the order of the tensor. Thus, for example, the gradient of
a vector fi is the second-order tensor @fi =@xj . The divergence operator decreases the order of a
tensor (for this reason it cannot be applied to a scalar). The divergence of a tensor Tij is the vector
@Tij =@xi . The tensor version of Gauss's theorem, equation (1.58), is
Z

@Tij
dV
A
V @xi
An isotropic tensor is one that is invariant to rotation of the coordinate system.
Tij ni dA =

Example 1.6

(1.74)

Show that the tensor represented by ij is0 isotropic.

The tensor has components ij in S . In S , the components of this tensor will be Aik Ajl kl = Ail Ajl = ij .
Thus this tensor has the same form in the rotated coordinate system.
20

A symmetric tensor is one for which Tij = Tji , while for an anti-symmetric tensor Tij = Tji .
Since
1
1
(1.75)
Tij = (Tij + Tji ) + (Tij Tji )
2
{z
}
{z
}
|
|2
symmetric

anti symmetric

any tensor can be written as the sum of two tensors, one symmetric and the other anti-symmetric.

1.6.1 Symmetric tensors

For a any tensor or matrix T an eigenvalue  is obtained from the equation

Tz = z

(1.76)

where z is a non-zero vector called eigenvector, and  is an eigenvalue. There can be more than
one eigenvalue and eigenvector for a given tensor. For a symmetric tensor the eigenvalues are real
numbers. Taking z = [a b c]T , the scalar equations corresponding to the vector equation above are
(T11 )a + T12 b + T13 c = 0
T21 a + (T22 )b + T23 c = 0
T31 a + T32 b + (T33 )c = 0
For a nontrivial solution for a, b, and c, we must have







(1.77)

T11  T12
T13 
T21
T22  T23  = 0
T31
T32
T33  

(1.78)

The cubic equation in  thus obtained is called a secular equation, and its three solutions give the
eigenvalues 1 ; 2 ; 3 . For each i equations (1.77) provide a nonunique eigenvector zi .
The three eigenvectors z1 ; z2 ; z3 form an orthogonal coordinate system that is rotated with
respect to the original. With unit vectors along these directions (called principal axes), the original
tensor becomes diagonal
2
3
1 0 0
4 0 2 0 5
(1.79)
0 0 3
There is thus considerable simplication obtained as a result of using principal axes as coordinate
directions. For instance, the equation
Tij xi xj = 1
(1.80)
describing a quadric surface, is really equivalent to

T11 x21 + T12 x1 x2 + T13x1 x3 + T21 x2 x1 + T22 x22 + T23 x2 x3 + T31 x3 x1 + T32 x3 x2 + T33 x23 = 1 (1.81)
with a number of cross terms xi xj (i 6= j ). For a symmetric Tij , one can change to new coordinates
along its principal directions to get an equation for the surface of the form

1 x21 + 2 x22 + 3 x23 = 1

21

(1.82)

which can be interpreted geometrically as an ellipsoid, a paraboloid, or a hyperboloid.

Example 1.7

Find the eigenvalues and eigenvectors of

"

1 0 0
0 1 2
0 2 1

Eigenvalues are solutions of the equation



 1 
0
0 

1  2  = 0
 0
 0
2 1 
which gives
(1 )( 3 + )(1 + ) = 0
so that 1 = 1; 2 = 1; 3 = 3.
Equations (1.77) become
(1 i )ai = 0
(1 i )bi + 2ci = 0
2bi + (1 i )ci = 0
For 1 = 1: 2a1 = 0; 2b1 + 2c1 = 0; 2b1 + 2c1 = 0. This gives a1 = 0; b1 = c1 . Choose b1 = 1
arbitrarily. Thepeigenvector is e2 e3 . To make it a unit vector, and remove some arbitrariness,
we divide by
its magnitude 2. So the unit eigenvector corresponding to this eigenvalue is (e2 e3 )=p2.
For 2 = 1: 2c2 = 0; 2b2 = 0. This gives b2 = c2 = 0, while a1 is arbitrary. Choose a1 = 1 to make it a
unit vector. Then the eigenvector corresponding to this eigenvalue is e1 .
For 3 = 3: 2a3 = 0; 2b3 + 2c3 = 0; 2b3 2c3 = 0. This gives a3 = 0; b3 = c3 . Choose b3 = 1. Then
the eigenvector is e2 + e3 . To make it a unitpvector, we divide by its magnitude p2. So the unit eigenvector
corresponding to this eigenvalue is (e2 + e3 )= 2).
Vectors in directions opposite to the ones chosen here would also be unit eigenvectors.

1.6.2 Anti-symmetric tensors

The diagonal terms of an anti-symmetric tensor Tij are zero and it has only three independent terms
as shown:
2
3
0
T21 T13
T32 5
T = 4 T21 0
(1.83)
T13 T32
0
The three terms can be made to constitute a vector. Thus, we can let

a=T e

32 1

This vector has the property that

where b is any vector.

+ T13e2 + T21 e3

T b=ab
22

(1.84)
(1.85)

Im(z )

Re(z )

Figure 1.7: Complex number z = x + iy in Argand plane.

1.7 Complex numbers and functions

A complex number z has the form

z = x + iy

(1.86)

1. x is the real part of z and y is its imaginary part. Keeping the

where u and v are real, and i =
denition of i in mind, these numbers can be added, subtracted, multiplied, and divided as usual.
Two complex numbers are equal only if their real parts are equal and the imaginary parts are also
equal.
A complex number can be geometrically represented in a plane (called the Argand plane) as
shown in Fig. 1.7. The real part of z is plotted on the abscissa and the imaginary part on the
ordinate.
The modulus or absolute value of the complex number is
p

jz j = + x

+ y2

(1.87)

which is the distance r in the gure.

The conjugate of this number, z, is dened by
so that

z = x iy

(1.88)

jz j

(1.89)

= z z

1.7.1 Analytic functions

A real function f (x) is said to be analytic at x = x0 if f and all its derivatives exist at this point.
An analytic function can be expanded in a Taylor series in the neighborhood of the point at which
it is analytic.
A function of a complex variable, or complex function for short, dened in some domain is of
the form f (z ) where z is a complex number. A derivative of a complex function, f (z ), at z = z0 is
dened as
df
f (z ) f (z0 )
= zlim
(1.90)
!
z
dz
z z0
0
23

In contrast to the case of real functions, in this limiting process the approach z ! z0 can be in
dierent ways and the derivative exists only if all of these give the same value. A complex function
is analytic at a point if its derivative exists at that point.
Let us suppose that the function f (z ) is analytic at z = z0 . Writing z = x + iy, we have f as a
complex function of the real variables x and y. Thus

f0 =
=

df
dz

(1.91)

df
d(x + iy)

(1.92)

As a special case the approach z ! z0 can be by changing x only, which gives

@f
@x
On the other hand we can approach z0 by changing y only, giving
f0 =

(1.93)

@f
@y

(1.94)

f = u + iv

(1.95)

f0 = i
If we separate f into its real and imaginary parts
we get from the above that

@v
@u
=
@x
@y
@u
@v
=
@y
@x

(1.96)
(1.97)

These are the Cauchy-Riemann equations. If the partial derivatives exist, these are necessary and
sucient conditions for a complex function to be analytic at a point. For a complex function it can
be shown that if f 0 (z ) exists, then so do all other derivatives.
From the above relations we can show that

@2u @2u
+
= 0
@x2 @y2
@2v @2v
+
= 0
@x2 @y2

(1.98)
(1.99)

The real and imaginary parts of a complex function are harmonic, i.e. they satisfy Laplace's equation.
It can happen that a function is analytic in a neighborhood of point z0 , but not at z0 itself. A
singular point is one at which the function is not analytic.

1.7.2 Polar form

A complex number may also be written in polar form using Euler's formula

z = rei
24

(1.100)

from which we can see that

z = re i
From Fig. 1.7 the relation between the Cartesian and polar forms is

Thus, we have Euler's relation

(1.101)

x = r cos 
y = r sin 

(1.102)
(1.103)

ei = cos  + i sin 

(1.104)

and de Moivre's theorem

ein = cos n + i sin n

(1.105)
Examples of use of the polar form include the following:
(a) Write ei(A+B) = eiA eiB and expand both sides. The real and imaginary parts give the trigonometric relations for cos(A + B ) and sin(A + B ).
(b) The relation z n = rn ein can be used to nd the roots of complex numbers.
(c) The real and imaginary parts of the logarithm of a complex number may be found by using
ln z = ln r + i.
To nd the Cauchy-Riemann conditions in polar coordinates, we write the derivative of a complex
function as
df
f0 =
(1.106)
d(rei )
Moving in the radial direction, this becomes
@f
f 0 = e i
@r
and in the circumferential direction it is
ie i @f
f0 =
r @
For a complex function of the form f = u + iv, we get
1 @v
@u
=
(1.107)
@r
r @
@v
1 @u
=
(1.108)
r @
@r

1.7.3 Integrals

If the function f (z ) is analytic inside and on a closed contour C , then according to Cauchy-Goursat's
theorem
Z
f (z ) dz = 0
(1.109)
C
Furthermore, by Cauchy's theorem
Z
dn f
n!
f (z )
(z ) =
dz for n  1
(1.110)
dz n 0 2i C (z z0 )n+1
and
Z
1
f (z )
f (z0 ) =
dz
(1.111)
2i C (z z0 )
25

Problems

1. Using the identity

derive

ijk lmn = il jm kn + im jn kl + in jl km
ijk lmk
ijk ljk
ijk ijk

= il jm
= 2il
= 6

il jn km

im jl kn

in jm kl

im jl

2. Prove by expanding into components that

(a) a  (b  c) = b(a
c) c(a
b)
(b) r  r = 0
3. Using a  (b  c) = b(a
c) c(a
b) show that ijk lmk = il jm im jl .
4. Show that
r  (r  f) = r (r
f) r2 f
5. Find the principal axes of the symmetric tensor
"
2 4
6 #
4 2
6
6 6 15
Conrm using MAPLE.
6. Show that equation (1.85) holds for any vector b, where T and a are given by equations (1.83) and (1.84)
respectively.
7. Write the expressions for r, r
f, and r  f in (a) cylindrical polar coordinates, and (b) spherical polar
coordinates.
8. Find the principal axes of the tensor
"
3 0 0#
0 3 4
0 4 3
9. Apply Stokes's theorem to the plane vector eld f(x; y) = fx (x;H y)i + fy (x; y)j and a closed curve enclosing a
plane region. What is the result called? Use this result to nd C f
dl, where f = yi+ xj and the integration
is counterclockwise along the sides C of the parallelogram with corners at (0,0), (1,0), (2,1), and (1,1).
10. Two-dimensional bipolar coordinates (; ) are dened by
 sinh 
x =
cosh  cos 
 sin 
y =
cosh  cos 
Find expressions for the gradient and divergence operators. Also computer generate some constant  and
constant  lines.
11. Using the index notation, show that
r  (f  g) = (g
r)f (f
r)g + f(r
g) g(r
f)
where f and g are vector elds.
12. Consider
two Cartesian coordinate systems: S with unit vectors (i; j; k), and S 0 with (i0 ; j0 ; k0 ), where i0 = i,
j0 = (j k)=p2, k0 = (j + k)=p2. The tensor T has the following components in S :
1 0 0!
0 1 0
0 0 3
0
Find its components in S .
26

13. Find the matrix A that operates on any vector of unit length in the x-y plane and turns it through an angle 
around the z-axis without changing its length.
14. Prove the following identities using index notation:
(a) (a  b)
c = a
(b  c)
(b) a  (b  c) = b(a
c) c(a
b)
(c) (a  b)
(c  d) = (a  b)  c
d
15. The position of a point is given by R = ia cos !t + jb sin !t. Show that the path of the point is an ellipse.
Find its velocity V and show that R  V = constant. Show also that the acceleration of the point is directed
towards the origin and its magnitude
is proportional to the distance from the origin.
H
16. Use Green's theorem to calculate C f
dl, where f = x2 i + 2xyj, and C is the counterclockwise path around a
rectangle with vertices at (0,0), (6,0), (0,4) and (6,4).
17. Derive an expression for the divergence of a vector in orthogonal paraboloidal coordinates
x = uv cos 
y = uv sin 
1 2 2
z =
2 (u v )
Determine the scale factors. Find r, r
f, r  f, and r2  in this coordinate system.
18. Derive an expression for the gradient, divergence, curl and Laplacian operators in orthogonal parabolic cylindrical coordinates (u; v; w) where
x = uv
1 (u2 v2 )
y =
2
z = w
19.

20.
21.

22.
23.
24.

where 0  u < 1, 1 < v < 1, and 1 < w < 1.

Orthogonal elliptic cylindrical coordinates (u; v; z) are related to Cartesian coordinates (x; y; z) by
x = a cosh u cos v
y = a sinh u sin v
z = z
where 0  u < 1, 0  v < 2 and 1 < z < 1: Determine r; r
f ; r  f and r2  in this system , where
 is a scalar eld and f is a vector eld.
Determine the unit vector normal to the surface x3 xyz + z3 = 1 at the point (1,1,1).
Show using indicial notation that
r  r = = 0
r
rf = 0
r(f
g) = (f
r)g + (g
r)f + f  (r  g) + g  (r  f)
1 r(f
f) = (f
r)f + f  (r  f)
2
r
(f  g) = g
r  f f
r  g
r  (r  f) = r(r
f) r2 f
r  (f  g) = (g
r)f (f
r)g + f(r
g) g(r
f)
Show that the Laplacian operator @ 2 =@xi @xi has the same form in S and S 0 .
Show the identities (1.45)-(1.50).
For potential ow the velocity eld is given by u = r, where  is the velocity potential. Consider the velocity
potential
(x; y; z ) = x2 + y2 2z 2
(a) Show that it satises Laplace's equation r2  = 0.
27

25.

26.

27.
28.
29.
30.
31.

(b) Find the velocity eld u.

(c) Show that u has zero divergence.
Show that the following scalar eld in cylindrical coordinates (which represents the stream function for ow
around cylinders)


2
= UR 1 Ra 2 sin 
where U and a are constants, satises Laplace's equation.
Show that the following scalar eld in spherical coordinates (which represents the velocity potential for ow
around spheres)


1 a3 cos 
 = Ur 1
2 r3
where U and a are constants, satises Laplace's equation.
Show that the Cartesian vector eld
u = (y + z)i + (z + x)j + (x + y)k
is divergence free and irrotational. Also nd its gradient.
Find expressions for the gradient, divergence and curl operators in spherical coordinates.
Find expressions for the gradient, divergence and curl operators in cylindrical coordinates.
Show that the sum of two analytic functions is also analytic.
Determine the integral
Z
C

z 2 dz

where C is the path from (0,0) to (1,1) in the three dierent ways:
(a) straight lines from (0,0) to (1,0) to (1,1)
(b) straight line from (0,0) to (1,1)
(c) straight lines from (0,0) to (0,1) to (1,1).
32. If z and z0 are the same complex number but in the unprimed and primed coordinate system, respectively,
show that
z = z 0 ei
33. Prove Cauchy-Goursat's theorem1 .
34. Show that
(r; ) = a0 + a1 ln r + a2  + a3  ln r +

in polar coordinates is a harmonic function.

1 h
X
n=1

Consult your favorite book on complex variable theory.

28


i
B 
D 
An rn + nn cos n + Cn rn + nn sin n
r
r

y0

BB
B
B

BB

 x

BB


B
x
0

Figure 1.8: Coordinate systems S and S 0

29

30

Chapter 2

Flow, rate of deformation and

stress
2.1 Flow lines
2.1.1 Streamlines
These are lines the tangents to which are in the direction of the velocity vector. Thus, in Cartesian
coordinates
dx dy dz
= =
(2.1)
u
v
w
is the equation of the streamlines, where u = uex + vey + wez .

2.1.2 Pathlines
These are lines traced out by the uid particles. Thus

dx
= u
dt
dy
= v
dt
dz
= w
dt
where (x; y; z ) is the position of a uid particle at time t.

(2.2)
(2.3)
(2.4)

2.1.3 Streaklines
These are lines traced out by particles that have all passed through a given point at some previous
time.

2.1.4 Timeline
This is the locus of a set of particles that initially denes a line.
31

Example 2.1

For a ow with velocity eld u = U cos !ti + U sin !tj, nd (a) the streamline passing through the origin
at time t =  , (b) pathline of particle that is at the origin when t = 0, and (c) streakline at time t =  for a
source at the origin.
(a) The velocity eld at time t =  is u = U cos ! i+ U sin ! j. The equation for the streamline is a solution of
dx
= dy
U cos ! U sin !
This gives
Ux sin !t = Uy cos !t + a
Since x = 0 for y = 0, we get a = 0. Thus the required streamline is
y = x tan !
(b) The pathline is a solution of
dx
= U cos !t
dt
dy
= U sin !t
dt
Thus
U
sin !t + x0
x =
!
U
y =
cos !t + y0
!
Since the particle is at (0,0) when t = 0, we have x0 = 0, y0 = U=!. The parametric form of the pathline is
U
x =
sin !t
!
U
y =
(1 cos !t)
!
Eliminating t, we get the equation
2
!2 2  !
x
+
y
1
=1
U2
U
This is a circle of radius U=! and center (0; U=!).
(c) The particle position at time t =  is given by
U
x =
sin ! + x0
!
U
y =
cos ! + y0
!
Since the particles that compose the streakline have been at the origin at a time t = s (say), we have
0 = U! sin !s + x0
0 = U! cos !s + y0
Using these values of x0 and y0 , we get the parametric equation of the streakline at time t = 
U
x =
(sin ! sin !s)
!
U
y =
(cos ! + cos !s)
!
Eliminating s, we have

2 
2
!
!
x sin ! +
y + cos ! = 1
U
U
which is a circle of radius U! and center ( U! sin !; U! cos ! ).
32

C





 0

d 
B
0 d
D0

x2

C
dx2

dx1

x1

Figure 2.1: Deformation of uid material

2.2 Rate of deformation

The gradient of the velocity vector is the deformation rate tensor D, where its components are

Dij =

@ui
@xj

(2.5)

We can split it up into the sum of an antisymmetric and a symmetric tensor, so that

Dij =
ij + "ij
where
1 @ui

ij =
2 @xj

@uj
@xi

1 @ui @uj
"ij =
+
2 @xj @xi

(2.6)
!

(2.7)
!

(2.8)

ij and "ij are called the rotation rate and shear rate tensors respectively.
The vector corresponding to twice the anti-symmetric tensor
ij is the vorticity given by

=ru

(2.9)

A velocity eld for which ! is zero everywhere is said to be irrotational. Physically, the vorticity
vector ! is twice the rotation rate of innitesimal uid elements and the symmetric shear rate tensor
"ij represents the rate at which these uid elements are being deformed. This is shown below.
33

Consider the uid in a material rectangle in Figure 2.1 of size dx1  dx2 moving with the ow.
After an interval dt the points A, B , C , and D move to A0 , B 0 , C 0 , and D0 . The instantaneous rates
of rotation of sides AB and AD are
d
@u2
=
(2.10)
dt
@x1
d
@u1
=
(2.11)
dt
@x2
Thus we have
1 @u2
Rotation rate =
2 @x1

@u1
@x2

1 @u2 @u1
Shear rate =
+
2 @x1 @x2

(2.12)
!

(2.13)

The shear rate in all three coordinate planes are the o-diagonal components of the shear rate
tensor 21 (@ui =@xj + @uj =@xi ). The diagonal components of this tensor are the extensional strain
rates @u1=@x1 , @u2=@x2 , and @u3 =@x3 .

2.3 Vortex lines and tubes

Tangents to vortex lines are in the direction of the vorticity vector. Thus

dx dy dz
=
=
!x !y !z

(2.14)

is the equation of a vortex line, where

! = !x x + !y y + !z z

(2.15)

2.4 Eulerian and Lagrangian descriptions

In the Eulerian approach, we describe the ow eld by its properties as a function of space and
time. In the Lagrangian approach, however, we follow a uid \particle" and observe the change
of properties with time as it is swept through the ow. Newton's laws are usually formulated in
a Lagrangian sense, while the uid equations as well as measurements are easier handled with the
Eulerian description. Thus some relation between the two descriptions is necessary.

2.4.1 Material derivative

In this, as in other sections in this chapter, we will derive the appropriate expressions in Cartesian
coordinates, and then generalize to other coordinate systems. Consider a property (x1 ; x2 ; x3 ; t),
scalar, vector or tensor, which is function of position expressed in Cartesian coordinates x1 ; x2 ; x3
and of time t. If we follow a uid particle, the change in this property after a time dt will be d,
where
@
@
@
@
dt +
dx1 +
dx2 +
dx
(2.16)
d =
@t
@x1
@x2
@x3 3
34

where dx1 ; dx2 ; dx3 are the components of the displacement of the particle in the coordinate
directions. Dividing by dt, we have

d @
@
@
@
=
+ u1
+ u2
+ u3
dt @t
@x1
@x2
@x3

(2.17)

where

dx1
= u1
dt
dx2
= u2
dt
dx3
= u3
dt

(2.18)
(2.19)
(2.20)

Frequently, the symbol D=Dt is used instead of d=dt to emphasize the fact that the time derivative
is calculated following a uid particle. It is called a material, total, or substantial derivative. Thus

D
@
@
= +u
Dt @t i @xi

(2.21)

The rst term on the right hand side is called the local derivative, and the second the convective
derivative. For other coordinates we can use

D
@
= +u
r
Dt @t

(2.22)

where u is the velocity vector.

2.5 Stress
The stress is dened as the force per unit area. It is a tensor since both force and area can be
considered to be vectors (an elemental area has orientation). Consider an elemental volume of uid
as shown in Figure 2.2. There is a force on each one of the six surfaces. Let the stress tensor at the
center of the element P be  , so that those on the two faces normal to the x1 direction are


= 

@  dx1
@x1 2

The components of the force per unit area on the face with normal in the +x1 direction are marked
+
+
+
as 11
, 12
, and 13
, while those on the face with normal in the x1 direction are 11 , 12 13 . The
stresses on the + and faces may, in general, be dierent.  and  + tend to  as the size of
the element goes to zero. The components of the stress tensor at P can thus be represented as ij ,
where the rst index represents the face on which the force is acting, and the second the direction
of the force.
If  is the state of stress at a point in the uid, the force per unit area on an element of surface
with unit normal n, sometimes called the traction, is 
n, or ij nj . If the area of the element is
dA, the force vector is 
n dA, or ij nj dA.
35

x2
dx3

+
12

13

11


dx2


?

x3

-

+
11

x1

+
13

12

dx1

Figure 2.2: Forces per unit area on elemental volume

Problems

1. Find an expression for the acceleration of a particle Du=Dt for the two-dimensional velocity eld
(a) u = x1 e1 x2 e2 in Cartesian coordinates, and
(b) u = U (1 a2 =r2 ) cos er U (1 + a2 =r2 ) sin e in cylindrical coordinates.
2. Show that 21 curl u is the vector corresponding to the anti-symmetric vorticity tensor
ij .
3. For the velocity eld u = x1 e1 x2 e2 , nd the deformation-rate, vorticity, and strain-rate tensors.
4. For an axisymmetric ow in spherical coordinates (i.e. with u! = @=@! = 0) show that the only nonzero
component of the vorticity,  , is in the direction of the axis of symmetry, and that
"
#
1
@ (ru ) @ur
=
r

@r

@

5. Find an expression for the acceleration of a particle in the two-dimensional velocity eld u = x22 e1 + x21 e2 .
6. For the two-dimensional velocity eld u = x1 e1 x2 e2 , nd (a) the shape of a portion of uid that is initially
a rectangle with corners at (a; b), (a + a; b), (a + a; b + b) and (a; b + b), where a, a, b and b are all
non-negative, after an interval of time t. Taking the appropriate limits, nd (b) the rate of change in area,
(c) the rates of rotation of the two diagonals of the rectangle and their average. Repeat for the velocity eld
u = x2 i1 .
7. For the velocity eld u = (3x1 + x2 +2x3 )e1 +( x1 +3x2 +7x3 )e2 +( 2x1 + x2 +3x3 )e3 nd the deformation
rate, vorticity, and shear rate tensors. Find the principal axes of the shear rate tensor. Check using Maple.
8. Show that
Du @ u 1
= + r (u
u) + (r  u)  u
Dt @t 2
9. For the two-dimensional velocity eld u = x1 e1 x2 e2 , nd the time derivative of the velocity determined
by a relative-velocity
sensor (i.e. the sensor measures the uid velocity relative to it) moving with velocity
v = te1 + t2 e2 . The sensor is at the origin at time t = 0.
10. For the velocity eld u = x1 e1 x2 e2 , nd the equations governing the motion of a skywriting airplane that
is required to produce a perfect circle of unit radius centered at (1,1)
at time t = T . The plane moves at a
constant speed V and begins to write at t = 0 from the point (2e T ; eT ).
36

11. Find the force on an elemental area of size dA and unit normal n, if the stress tensor at that location is  .
12. A temperature sensor starts from the origin and moves with velocity us = ti + t2 j, where t is time. Find its
reading, Ts (t), in a time-dependent temperature eld T (x; y; t) = txy.
13. Find the shape of a material line at time t that is initially a circle with unit radius and center (2,2) moving
with the ow eld u = xi yj.
14. The surface ij xi xj = 1, where ij is symmetric, is called the Cauchy stress quadric. Show that the normal
component of the stress on any plane is inversely proportional to the square of the distance from the center of
the quadric to its surface in the direction of the normal to the plane.

37

38

Chapter 3

Dierential equations of motion

3.1 Approaches to balance equations
There are dierent ways in which the balance equations for uid mechanics may be derived. The
approach may be dierential, i.e. applied to a element, or it may be integral, i.e. applied on a
larger scale. Furthermore, the application may be to a xed mass of uid that is moving around
and changing its shape in space, or to a xed region in space through which the uid ows in and
out. So there are at least four ways which must be somehow related.
Consider a (scalar, vector or tensor) property  that is to be balanced.

3.1.1 Dierential element

Though the following may be easily carried out for any coordinate system, we will analyze the
situation in Cartesian coordinates. Consider the box-like element shown in Fig. 3.1.

Control volume
The box represents a region that is xed in space and unchanging. Thus

d
@
( dx1 dx2 dx3 ) = () dx1 dx2 dx3
dt
@t

(3.1)

Control mass
Now the box is an element of uid that is moving and changing its shape and volume. Thus we have

D
D
D
( dx1 dx2 dx3 ) = () dx1 dx2 dx3 +  (dx1 dx2 dx3 )
Dt
Dt
Dt
We also have
D
D
D
D
(dx dx dx ) = dx2 dx3 (dx1 ) + dx3 dx1 (dx2 ) + dx1 dx2 (dx3 )
Dt 1 2 3
Dt
Dt
Dt






@u1
@u2
@u3
= dx2 dx3
dx + dx3 dx1
dx + dx1 dx2
dx
@x1 1
@x2 2
@x3 3
= dx1 dx2 dx3 r
u
39

x2
dx3
P
dx2

x1

x3
dx1

Figure 3.1: Element in Cartesian coordinates

from which

D
D
( dx1 dx2 dx3 ) =
() + r
u dx1 dx2 dx3
Dt
Dt

(3.2)

3.1.2 Integral element

We consider now a volume V , not necessarily small.

Control volume
The volume is a xed region in space. Thus
Z
Z
@
d
 dV =
() dV
dt V
@t
V
is the rate of change of the total property within the volume.

(3.3)

Control mass
If V represents a uid material, its shape, volume and position will change with time. The Reynolds
transport theorem

Z
Z 
D
D
 dV =
() + r
u dV
(3.4)
Dt V
V Dt
describes the rate of change of the total property of this material as it moves around.

3.2 Mass balance

We will derive this equations in its dierent forms and show that they are all equivalent.
40

x2
dx3

-m_

m_

dx2

+
1

x1

x3
dx1

Figure 3.2: Mass ux in elemental volume

3.2.1 Elemental control volume

Consider an elemental control volume within the ow, Figure 3.2, of size dx1  dx2  dx3 centered
at P (x1 ; x2 ; x3 ) and aligned with the coordinate directions. There is mass ux through all six faces
of this volume. The mass ux through each face is equal to the product of the local uid density,
the velocity normal to the face, and the area of the face. The uxes in the two faces with normals
in the  x1 direction are
"

m_

@ (u1 ) dx1
u1 +
+ : : : dx2 dx3
@x1 2

+
1

"

m_ 1

u1

@ (u1 ) dx1
@x1 2

(3.5)

: : : dx2 dx3

(3.6)

where  is the uid density, and u1 is the velocity component at P . The net gain of mass per unit
time by the control volume due to ow in this direction is

@ (u1 )
dx1 dx2 dx3 + : : :
@x1
There are similar gains due to ow in the other two directions. The total is
m_ 1

"

m_ +1 =

(3.7)

@ (u1 ) @ (u2 ) @ (u3)

+
+
dx1 dx2 dx3
@x1
@x2
@x3

On the other hand there is accumulation of mass within the control volume, the rate of which is
@
dx dx dx
@t 1 2 3
41

The rate of accumulation of mass should be equal to the total gain due to ow. Thus,

@
dx dx dx =
@t 1 2 3
from which

"

@ (u1) @ (u2 ) @ (u3 )

+
+
dx1 dx2 dx3
@x1
@x2
@x3

(3.8)

@ @ (u1) @ (u2 ) @ (u3 )

+
+
+
=0
@t
@x1
@x2
@x3

(3.9)

@ @ (ui )
+
=0
@t
@xi

(3.10)

@
+ r
(u) = 0
@t

(3.11)

@u
D
+ i =0
Dt
@xi

(3.12)

In index notation
For other coordinate systems
Another way of writing this equation is

or

D
+ r
u = 0
Dt
The equation of conservation of mass is often called the continuity equation.

(3.13)

3.2.2 Elemental control mass

Following a uid element that is initially dx1  dx2  dx3 , the mass of this element must be constant.
Thus
D
( dx1 dx2 dx3 ) = 0
(3.14)
Dt

3.2.3 Integral control volume

The rate of change of mass with a control volume V plus the net mass ow out of it must be zero.
Thus
Z
Z
d
 dV + u
n dA = 0
(3.15)
dt V
A

3.2.4 Integral control mass

The rate of change of mass following a control mass of initial volume V is zero. So
Z

D
 dV = 0
Dt V
42

(3.16)

3.2.5 Equivalence
(a) Using equation (3.2) in (3.14), we get equation (3.13).
(b) Using the Stokes theorem on the second term of equation (3.15), we get
Z 

Z 

@
+ r
(u) dV = 0
(3.17)
V @t
Since this equation holds for any V , the integrand must be zero everywhere. This gives equation
(3.11).
(c) Using the Reynolds transport theorem, equation (3.16) becomes
D
+ r
u dV = 0
(3.18)
V Dt
Again since this equation holds for any V , the integrand must be zero everywhere. This gives
equation (3.13).

3.3 Linear momentum equation

We will use only an elemental control mass instantaneously of size dx1  dx2  dx3 . According to
Newton's second law in an inertial frame of reference, the product of its mass and its acceleration is
equal to the force on it. Let f be the force on the uid per unit volume. Then, writing this equation
per unit volume
Du

=f
(3.19)
Dt
The force per unit volume on a uid element can be considered to be the sum of body forces fb (such
as gravity or electromagnetic forces) which act over the entire volume of the element, and surface
forces per unit volume fs (such as stress) which act through the surface.
A conservative body force is one for which

fb = rG

(3.20)

where G is a scalar potential. If gravity is the only body force, then

fb = g

(3.21)

where g is the gravity force per unit mass.

The surface forces can be understood by considering a control mass shown in Figure 3.3. If
the components of  on the two faces with normals in the  x1 direction were the same, their net
contribution would be zero since they would be in opposite directions. So a uniform stress eld
produces no force on the uid element. For an inhomogeneous stress eld in which  is a function of
position, there is a net surface force. For the moment let us consider the forces in the  x1 direction
only. Each one of the six faces has a component in this direction. Thus
+
11
= 11 +

11 = 11
43

@11 dx1
@x1 2
@11 dx1
@x1 2

(3.22)
(3.23)
(3.24)

x2

+
22
dx3 dx1
6- 21+ dx3 dx1

11 dx2 dx3

+
12
dx2 dx3

-6


?

+
11

12 dx2 dx3

21 dx3 dx1

?
 dx
22

dx2 dx3

dx1
x1

Figure 3.3: Forces for linear and angular momentum equations

The pair of faces with normals in the  x1 direction contributes a force (@11 =@x1) dx1 dx2 dx3 .
Similarly, on including also the contributions from the other two pairs of faces, we have the total
surface force in this direction. Per unit volume, we get

fs1 =

@11 @21 @31

+
+
@x1 @x2 @x3

(3.25)

The linear momentum in the x1 direction gives

Du1 @11 @21 @31

=
+
+
+ f1
Dt
@x1 @x2 @x3

(3.26)

Generalizing this expression to include components of linear momentum in the other directions
also, we have
@
Du
 j = ij + fj
(3.27)
Dt
@xi
which, for other coordinate systems, is

Du
= r
 + f
(3.28)
Dt
As a reminder, the left hand side represents the product of the mass times the acceleration of an
element of uid, while the two terms on the right are the surface and body forces respectively on
this element, all per unit volume.
Multiplying the continuity equation, equation (3.10), by ui and adding to equation (3.27), we
get
@
@
@
(uj ) +
(ui uj ) = ij + fj
(3.29)
@t
@xi
@xi


44

or

@
(u) + r
(uu) = r
 + f
@t
which is another form of the momentum equation.
The linear momentum equation is also called the Cauchy equation.

(3.30)

3.4 Linear momentum equation in a non-inertial frame

If the coordinate system with respect to which the velocities and positions are measured is noninertial, i.e. rotating and/or linearly accelerating, other terms have to be included in the momentum
equation.
Consider a non-inertial frame characterized by r0 , u0 , and a0 which are the position, velocity,
and acceleration of its origin, and
, and
_ which are its rotation rate, and angular acceleration, all
with respect to an inertial frame. The positions of a particle measured in the two frames are then
related by
r0 = r0 + r
(3.31)
where r0 and r are in the inertial and non-inertial frames, respectively. Similarly velocities are related
by
u0 = u0 + u +
 r
(3.32)
and accelerations by

a0 = a

+ a +
_  r + 2
 u +
 (
 r)

(3.33)

The linear momentum equation is then


or

@
Dui
+ ai + jki
j rk + 2jki
j uk + lmi njm
l
n rj = ij + fj
Dt
@xi


Du
+ a +
_  r + 2
 u +
 (
 r) = r
 + f

Dt

(3.34)
(3.35)

2
_  u is the Coriolis and
 (
 r) the centripetal acceleration term.

3.5 Angular momentum equation

Consider an elemental control mass of instantaneous size dx1  dx2  dx3 . Figure 3.3 shows the
stresses on a section of a uid element normal to the x3 axis. We can use Euler's equation which
govern the angular dynamics of a rigid body1 . For the x3 axis

d
3
+ (I2 I1 )
1
2
(3.36)
dt
where Mi ; Ii ;
i are the total moment, the moment of inertia and angular velocity of rotation
about the principal axis xi . We will assume that the uid is non-polar, i.e. that there are no body
moments, so that the only contribution to the moment comes from the stresses.
M3 = I3

See, for instance, H. Goldstein,

, Addison Wesley, Reading, MA, 1950.

Classical Mechanics

45

The moment of inertia about the x3 axis is

1
I3 =
(mass of uid element)(dx21 + dx22 )
12

=
(dx dx dx )(dx21 + dx22 )
12 1 2 3
Similarly for I1 and I2 . Equation (3.36) becomes

dx
@ dx
dx
@12 dx1
) dx2 dx3 1 + (12 + 12 1 ) dx2 dx3 1
@x1 2
2
@x1 2
2
@21 dx2
dx2
@21 dx2
dx2
(21
) dx1 dx3
(21 +
) dx1 dx3
@x2 2
2
@x2 2
2

2
2 d
3
dx dx dx (dx + dx2 )
=
12 1 2 3 1
dt

2
2
+ dx1 dx2 dx3 [(dx1 + dx3 ) (dx22 + dx23 )]
1
2
12
We divide by dx1 dx2 dx3 and take the limit as dx1 , dx2 , and dx3 go to zero, to obtain

(3.37)

(12

12 = 21

(3.38)
(3.39)

This can be repeated for the other two coordinate axes giving 13 = 31 and 32 = 23 . Thus, the
conservation of angular momentum of a uid element implies that the stress tensor is symmetric:
or

ij = ji

(3.40)

= T

(3.41)

3.6 Energy equation

We consider the elemental control mass shown in Figure 3.4. The total energy per unit mass of the
uid, et , is the sum of its internal and kinetic energy. Thus
1
(3.42)
et = e + ui ui
2
where e is the specic internal energy. The rate of in ow of energy in the form of heat and work is
equal to the rate of increase of total energy of the uid. Thus, per unit volume, we have
De
 t = Q_ + W_
(3.43)
Dt
where Q_ and W_ are the in ow of heat and work respectively per unit volume, the sign convention
being that energy coming into the system is considered positive.
The heat in ow is due to the spatial variation in the heat ux q_ , which is a vector in the direction
of the heat ow with units of energy per unit time per unit area. Let q_1 and q_1+ be the heat ux
components in the two faces which have normals in the  x1 direction. The ow of heat through
these faces are q_1 dx2 dx3 and q_1+ dx2 dx3 , respectively. Since

q_1+ = q_1 +

@ q_1
dx + : : :
@x1 1
46

(3.44)

x2
dx3
P

-q_

dx2

q_

+
1

x1

x3
dx1

Figure 3.4: Energy uxes in elemental volume

the net energy gain due to heat ow in the x1 -direction is

@ q_1
dx dx dx
@x1 1 2 3
Similarly, one can obtain the contributions from the heat ow in the other two directions. Summing
them up we get
@ q_i
Q_ =
(3.45)
@xi
on a unit volume basis.
Now let us consider the rate of work in ow W_ . We know that the rate of work done by a force
is the dot product of the force and velocity of its point of application. Let us look at the work done
by the surface forces rst. The rate of work done on the two face in the x1 -direction is
(11 u1 + 12 u2 + 13 u3 ) dx2 dx3 = 1j uj dx2 dx3
where the negative sign comes from the fact that the force and the velocity are in opposite directions.
The work done on the + face in the same direction is
+ +
+ +
+ +
(11
u1 + 12
u2 + 13
u3 ) dx2 dx3 = 1+j u+j dx2 dx3

Since
the work in ow is

1+j u+j = 1j uj +

@ (ij uj )
dx1
@x1

@ (1j uj )
dx1 dx2 dx3
@x1
47

On including the contributions from the other two pairs of faces also, we have the rate of work in ow
due to the surface forces W_ s , where
@ ( u )
W_ s = ij j
(3.46)
@xi
per unit volume. In addition the work done by the body force W_ b is

W_ b = ui fi

(3.47)

per unit volume. So

W_ = W_ s + W_ b
Using equations (3.45){(3.48), equation (3.43) becomes

or

@ q_i @ (ij uj )
Det
=
+
+ ui fi
Dt
@xi
@xi

(3.48)
(3.49)

Det
= r
q_ + r
(
u) + u
f
(3.50)
Dt
As a reminder, the left hand side represents the rate of change of total energy; the rst term on
the right is the heat in ow, the second is the work done by the surface forces, and the third the
work done by gravity, all per unit volume.
The above is the equation of conservation of total energy. It can be simplied by subtracting
out the mechanical energy part. Take the dot product of the velocity and the momentum equation
(3.27)
@
Du
(3.51)
uj j = uj ij + uj fj
Dt
@xi
This is the mechanical energy equation; subtract this from the total energy equation to get



or

De
@ q_i
@uj
=
+
Dt
@xi ij @xi

De
= r
q_ +  : ru
Dt
This is the thermal energy equation, sometimes simply referred to as the energy equation.


(3.52)
(3.53)

3.7 Entropy equation

In one form the second law of thermodynamics states that the change in entropy of a system, S , is
governed by
Q
S 
(3.54)
T
where Q is the heat input to the system, and T is its absolute temperature. The equality holds for
a reversible process.
Considering an elemental control mass of size dx1  dx2  dx3 , we nd that

Ds

Dt

@ q_i
@xi T
48

(3.55)

or

Ds

Dt
per unit volume, where s is the specic entropy.


 
r
Tq_

(3.56)

3.8 Additional relations

To complete the mathematical formulation of the problem, additional equations which involve properties specic to a particular uid material are needed.

3.8.1 Equations of state

Incompressible uid

In many occasions the density of uid following a particle is taken to be constant. In this case

D
=0
(3.57)
Dt
In general the density may be dierent for dierent particles, though the homogeneous uid is a
special case of this.

Ideal gas
This is a commonly used approximation for the behavior of gases:

p = RT

(3.58)

where R is the particular constant for the gas.

3.8.2 Constitutive equations

Inviscid uid

The constitutive relation for a zero viscosity uid is

ij = pij
or

= pI

(3.59)
(3.60)

The stress tensor is diagonal.

Newtonian uid
A constitutive relation for a uid is that which relates the stress and strain rate tensors. A Newtonian
uid has the following properties:
1. For a uid at rest the stress is hydrostatic, and the pressure is the thermodynamic pressure.

2. The stress tensor  is linearly related to the deformation-rate tensor D, and does not depend
on the rate of rotation tensor.
49

3. There are no preferred directions in the uid properties.

Under these conditions the constitutive relation for a Newtonian uid can be shown to be

ij =
or

@u
@ui @uj
+
p +  k ij + 
@xk
@xj @xi

= ( p + r
u) I + 

ru + (ru)T

(3.61)
(3.62)

where  and  are the rst and second coecients of viscosity.

Non-Newtonian uids
A uid which does not have a shear stress{rate of deformation relation given by equation (3.61) is
non-Newtonian. In a broad sense, the behavior of non-magnetic continua can be divided into several
categories:
(I) Purely viscous uids: The shear rate depends only on the shear stress. The uid can be Newtonian
or non-Newtonian.
(II) Time-dependent uids: The shear rate depends not only on the shear stress but also on the
duration of the stress.
(III) Viscoelastic materials: The shear rate depends on the imposed stress as well as the strain.
(IV) Complex rheological bodies: Materials displaying combinations of the characteristics above.

Fourier's law
The heat ux due to conduction heat transfer in a uid is governed by the relation

@T
@xi

(3.63)

q_ = krT

(3.64)

q_i = k
or

where T is the uid temperature, and k(T ) is the coecient of thermal conductivity.
There are uids in which a heat ux is produced by a concentration gradient also (diusionthermo or Dufour eect).

3.9 Governing equations for special cases

3.9.1 Mass balance
For an incompressible uid, equation (3.10) reduces to

or

@ui
=0
@xi

(3.65)

r
u=0

(3.66)

50

3.9.2 Linear momentum equation

Inviscid uid

The momentum equation becomes


or

Duj
@p
=
+ fj
Dt
@xj


This is called the Euler equation.

Du
=
Dt

(3.67)

rp + f

(3.68)

Newtonian uid
Substituting this into the momentum equation (3.27) gives the so-called Navier-Stokes equation


 

@p
@
@u
@
@ui @uj
Du
+
 k +

+
 j=
Dt
@xj @xj @xk
@xi
@xj @xi
or



+ fj

(3.69)

h n
oi
Du
= rp + r (r
u) + r
 ru + (ru)T + f
(3.70)
Dt
Because of the continuity equation (3.65), the constitutive relation for an incompressible Newtonian uid reduces to


@ui @uj
ij = pij + 
+
(3.71)
@xj @xi
or
h
i
 = pI +  ru + (ru)T
(3.72)

The second coecient of viscosity does not play a role in the mechanics of incompressible uids.
Taking  to be constant, the Navier-Stokes equation, equation (3.69), becomes


or

Duj
@p
@ 2 uj
=
+
+ fj
Dt
@xj
@xi xi

Du
= rp + r2 u + f
Dt
Often the equation is divided out by  and written as

or
where

(3.73)
(3.74)

Duj
1 @p
@ 2 uj
=
+
+f
Dt
 @xj
@xi xi j

(3.75)

Du
1
= rp +  r2 u + f
Dt


(3.76)

=
51




(3.77)

3.9.3 Energy equation

Inviscid uid

With Fourier's law the energy equation is


or

De
@2T
=k
Dt
@xi xi

@ui
@xi

(3.78)

De
= kr2 T pr
u
(3.79)
Dt
where k has been taken to be constant. The last term on the right is the rate of work due to the
pressure.


Newtonian uid
For a Fourier-Newtonian uid with constant properties, we have


or

@2T
De
=k
Dt
@xi @xi
De
= kr2 T
Dt

where
or

@ui
=
@xi

2

@ui
+
@xi

pr
u + 


@ui @uj @uj

+
+
@xj @xi @xi

 =  (r
u)2 + 

ru + (ru)T : ru

(3.80)
(3.81)
(3.82)
(3.83)

 is a positive quantity called the dissipation function; physically it represents the rate of change of
energy per unit volume from mechanical to thermal form.
For an incompressible uid, we can write
DT

= r2 T +
(3.84)
Dt
c
where e = cT and  = k=c, where c is the specic heat and  is the thermal diusivity.

3.9.4 Boundary conditions

Hydrodynamic and thermal boundary conditions are needed for the dierential equations.
 For a viscous uid the velocity at a solid wall is usually taken to be that of the wall itself. For
an inviscid uid only the normal velocities are equated.

At a free surface the stress is continuous across the surface if surface tension is not considered.
Otherwise there is a pressure dierential given by


1
1

p =
+
(3.85)
R1 R2
where R1 and R2 are the principal radii of curvature of the free surface, and is the coecient
of surface tension.
52

The temperature at the boundary is that of the wall itself. Sometimes, however, it is more
convenient to to prescribe the normal derivative of the temperature at the wall; for an adiabatic
wall it is zero, while for a wall with prescribed heat ux, it is a given quantity.

3.10 Nondimensionalization
The purpose of nondimensionalization is to normalize the variables and to bring out the relative
importance of the dierent terms in the equations. Often the procedure depends on the problem
being considered. Careful use of a method based on the dimensions of the physical quantities being
considered will often provide considerable qualitative and quantitative information of the ow.
As an example of nondimensionalization let us look at incompressible Navier-Stokes and energy
equations. Using asterisks for nondimensional quantities, we can choose the following dimensionless
variables
t = tU=L; x = x=L; u = u=U; p = p=U 2; T  = (T T0 )=
T
(3.86)
where U and L are characteristic velocity and length scales; T0 is a reference temperature and
T is
a characteristic temperature dierence. Writing f = geg , where eg is the dimensionless unit vector
in the direction of gravity, the equations become
Du
1
1
= r p + r2 u + eg
(3.87)
Dt
Re
Fr

DT
1
=
r2 T  + Ec

(3.88)
Dt
Re P r
Re
where r and r2 are the nondimensional gradient and Laplacian operators respectively, and 
is the nondimensional dissipation function. In this case the important parameters are Re, P r, and
Ec. These and some other nondimensional groups that commonly occur in uid mechanics and heat
transfer are listed below along with their physical signicance.
Name

Symbol

Denition

Bi
Ec
Fo
Gr
Nu
Pe
Pr
Ra
Re

hL=ks
U 2 =c
T
t=L2
g
T L3= 2
hL=k
cUL=k
c=k
g
T L3=
UL=

Biot number
Eckert number
Fourier number
Grashof number
Nusselt number
Peclet number
Prandtl number
Rayleigh number
Reynolds number

Physical signicance

ratio of solid to uid thermal resistance

ratio of dissipation to internal energy change
dimensionless time
ratio of buoyancy to viscous forces
dimensionless heat transfer coecient
ratio of convection to conduction heat transfer
ratio of momentum to thermal diusivities
= Gr P r
ratio of inertia to viscous forces

Additional nomenclature:

g
h
k
ks




acceleration due to gravity

convective heat transfer coecient
thermal conductivity of uid
thermal conductivity of solid
thermal diusivity
coecient of thermal expansion
kinematic viscosity = =
53

3.11 Molecular approach

See I. Michelson, Molecular Basis of Fluid Mechanics.

Problems

1. Show that  = constant, u = x1 e1 x2 e2 satisfy the continuity equation.

2. Show that the following velocity eld in polar coordinates
2
2
u = U (1 ar2 ) cos er U (1 + ar2 ) sin e :
is incompressible and irrotational.
3. Show that the two-dimensional velocity eld
Kx1
2
u = x2Kx
+ x2 e1 + x2 + x2 e2
1

4.
5.
6.
7.

where K is a constant, is incompressible and irrotational.

From an elemental control volume in cylindrical coordinates, derive the continuity equation in cylindrical
coordinates.
Show that the continuity equation can also be written as
D
+ r
u = 0:
Dt
Using equations (2.21) and (3.10), show that the dierence between equations (3.49) and (3.51) is (3.52).
Show that the dissipation function in equation (3.82) can also be written as
@ui
 =  @x
j

j
+ @u
@xi

@ui
@xj

j
+ @u
@xi

8. Show that the continuity equation for axisymmetric ow in spherical coordinates can be satised by the Stokes
stream function (r; ), where
1 @ ;u = 1 @
ur = 2
r sin  @ 
r sin  @r
9. Coordinates (; ;  ) in an elliptic cylindrical coordinate system are related to the Cartesian coordinates (x; y; z)
through x = a cosh  cos ; y = a sinh  sin ; z =  . If u , u and u are the uid velocities in the directions
of increasing , , and  respectively, write the continuity equation in elliptic cylindrical coordinates for an
incompressible uid.
10. For a Newtonian uid with velocity eld u = x1 e1 x2 e2 , nd the stress tensor and, by substituting into the
governing equations, the pressure distribution p(x; y).
11. Find the stress tensor for a Newtonian uid with the velocity eld u = x22 e1 + x21 e2 .
12. Neglecting viscous forces, nd the diameter of a laminar water jet coming down from a faucet as a function of
the distance downstream.
13. Find the stress tensor for a Newtonian uid with velocity eld u = x1 e1 x2 e2 in Cartesian coordinates.
14. For the toroidal coordinates (r; ; ) in the gure (the cross section has been enlarged to indicate r and ) show
that
x = (R + r cos ) cos ; y = (R + r cos ) sin ; z = r sin 
where (x; y; z) are suitably chosen Cartesian coordinates.
Find expressions for the gradient, divergence and curl in this coordinate system.
15. Given a two-dimensional velocity eld (x1 + x2 )e1 + (3x1 + 4x2 )e2 , show that the ow is compressible. Find
the strain rate and vorticity tensors, and (c) the principal axes of the strain rate tensor.
16. For steady ow and a divergence-free velocity eld, show that the ow lines are also constant density lines.
54

'$
i
&%
'$
n-n````
&%


Figure 3.5: Toroidal coordinate system (r; ; )

17. For the velocity eld u = x1 e1 x2 e2 , nd (a) the acceleration of the elemental volume of size dx1  dx2  dx3
at (x1 ; x2 ; x3 ), and (b) from Newton's second law, the net force on the volume required to produce this
acceleration.
18. Derive the linear momentum equation in cylindrical coordinates.
19. Show that


@
@
(u ) = @x ji uj ui
@t i
i

20. Show the theorem of stress means:

Z

F ij nj dA =

Z 

ij

@F
@xj

+ F

Dui
Dt

fi



dV

for any function F (xi ; t).

21. Starting from an elemental
in Cartesian coordinates, prove the momentum equation.
22. For each one of the N chemical species in a uid mixture, mi is the mass of one mole, ni is the number of
moles per unit volume of uid, ui is the velocity, and r_i is the generation by chemical reaction per unit volume
per unit time, where i = 1; : : : ; N . Show that the mass balance equation for each species is
@i
+ r
(i ui ) = r_i
@t
where i = mi ni (no sum). Show also that
@
+ r
(u) = 0
@t
control volume

where PNi=1 i =  and u = PNi=1 i ui .

23. Use the Reynolds's transport theorem to derive the mass, momentum and total energy equations.
24. The governing equations for mass, momentum and internal energy balance for any uid may be combined and
written in the form
@f @p @q @r
+ + + =s
@t

@x

@y

@z

in Cartesian coordinates. Find the column matrices f, p, q, r, and s in terms of the density and components
of the velocity vector, stress tensor, and body force vector.
25. For a Newtonian uid with velocity eld u = x1 e1 x2 e2 , nd, by substituting into the governing equations,
the pressure distribution p(x1 ; x2 ) and the stress tensor.
26. Show that the tensors
and

ij

ij pq +  ip jq + iq jp

ip jq

iq jp

are isotropic, where , , and are arbitrary scalars.

27. Write down the parameters that might aect the period of a swinging pendulum and nd, by dimensional
analysis, the period as a function of these parameters.
55

28. If r' and r are the position vectors of a particle in an inertial and non-inertial frame of reference, show that
d2 r0 d2 r
= dt2 + a + ddt
 r + 2
 ddtr +
 (
 r)
dt2
where a and
are the linear acceleration and rotation vectors of the non-inertial frame2.
29. The momentum equation for an incompressible, Newtonian uid in a coordinate system rotating at an angular
speed
is
@u
+ u
ru + 2
 u +
 (
 r) = 1 rp +  r2 u
@t
where r is a position vector. Nondimensionalize this equation, nd the nondimensional groups, and from
this determine if the rotation of the earth is important in the study of oceanographic phenomena? Choose
appropriate scales.
30. Show that for a steady ow, the no-slip condition at a xed wall requires that r
u also vanish at the wall.
31. For an incompressible, Newtonian uid with constant properties in a closed container of volume V , show that
the rate of decrease of total kinetic energy of the uid
is
Z


32.
33.
34.
35.

36.
37.
38.

!
!

dV

where ! is the vorticity vector. The body force is conservative. In a similar way, nd the rate of change of
total internal energy of the uid.
For steady, inviscid, barotropic ow show that Bernoulli's equation holds along a vortex line.
Using the vorticity equation for an inviscid uid, show that if the vorticity of a uid particle is once zero it is
always zero.
Find the stress tensor for a Newtonian uid corresponding to the velocityeld u = x1 e1 x2 e2 . Find also the
components of the same tensor in a coordinate system that is rotated 45 about the x3 -axis.
The momentum equation for an incompressible, Newtonian uid in a coordinate system rotating at a constant
angular speed
is
@u
+ u
ru + 2
 u +
 (
 r) = 1 rp +  r2 u
@t
where r is a position vector.
Show that the vorticity equation has an additional term 2(

r)u as compared
to the non-rotating case3 .
Applying the momentum equation Z
Z
d
u dV +
u (u
n) dA = F
dt CV
CS
over an elemental control volume in Cartesian coordinates, derive the momentum equation in dierential form.
Write the vorticity equation for an incompressible, Newtonian uid in component form using polar coordinates
(r; ) and the notation
u = ur er + u e
! = !z ez
(a) Show that the moment of the surface forces  (r) about a point O within an arbitrary volume V with surface
S is
Z
S

ijk rj km nm dA

where r is the position vector relative to O of an elemental surface area dA with normal n.
(b) Using Gauss's theorem write this as
"
#
Z
Z
0
@km
@ (rj km )
0

dV = 
 + rj
dV
V

ijk

@xm

ijk

kj

@xm

where  0 is the stress at O.

(c) As V ! 0, show that this reduces to ijk kj0 = 0, from which ij0 = ji0 .
2
Consult your favorite book on classical mechanics.
3 Hint: You may use Cartesian coordinates for some of the steps.
56

39. Show that the incompressible Navier-Stokes equation with constant properties can be written as
Du
= 1 rP  r  (r  u)
Dt
where
P = p G
f = rG
40. Consider the incompressible Navier-Stokes equation for a uid in coordinates that are rotating at a rate
=

k, where
is a 2constant. Neglect the body force. Use the following scalings: U for velocity, L for length, L=U
for time, and U for pressure to non-dimensionlize
the equation. Show that the Rossby number  = U=
L
and the Ekman number E = =
L2 appear as nondimensional parameters.
41. For ow in a coordinate system xed to the earth, nd numerical values for the Rossby and Ekman numbers
for (a) 1 m/s water velocity on a scale of 30 cm and (b) 100 km/hr wind at a scale of 100 km. For each indicate
if the rotation of the earth is important.
42. Using the Gibbs relation
 
1
T ds = de + p d

and the denition of enthalpy
p
h= e+

show that the energy equation for a Fourier-Newtonian uid can be written as
Ds
T
= r
( k rT ) + 
Dt
or
Dh
= r
(krT ) + Dp
+

Dt
Dt
where


 =  (r
u)2 +  ru + (ru)T : ru
43. Show that the 2nd law of thermodynamics for a Fourier-Newtonian uid can be written as
 + Tk rT
rT  0

57

58

Chapter 4

Special theorems
4.1 Circulation
The circulation

around a closed curve C is dened as

I

=
or

From Stokes's theorem, we have

so that

u
dl =
=

(4.1)

u
dl

(4.2)

r  u
n dA

(4.3)

where A is any surface with the boundary C .

ui dli

!

n dA

(4.4)

4.2 Kelvin's theorem

We consider an inviscid, barotropic uid (one for which the pressure is only a function of density)
with a conservative body force. The circulation around a closed material curve C is
=
so that
But

uj dlj

I 

D
Duj
D(dlj )
=
dlj + uj
Dt
Dt
C Dt
I

D(dlj )
=
uj
Dt
C
59

uj duj

(4.5)


(4.6)

= 0
Also

1
d uj uj
C 2

(4.7)
(4.8)

Duj
1 @p @G
=
+
Dt
 @xj @xj

Since

(4.9)

dp
1 @p
dxj =

@x
j
C 
C

and

(4.10)

@G
dxj =
dG
@x
C
C j
= 0

the only non-zero term is

(4.11)

D
dp
=
Dt
C 
For a barotropic ow for which  and p are uniquely related, this is also zero. Thus
D
=0
Dt

(4.12)

(4.13)

4.3 Vorticity equation

For a Newtonian uid with constant density and properties and conservative body force, the momentum equation is
 
@u
p
+ u
ru = r
+  + r2 u + rG
(4.14)
@t

Because of the identity
r (u
u) = 2u
ru + 2u  r  u
(4.15)
it can be written as


 
1
@u
+r u
u
u  ! = r p +  r2 u
(4.16)
@t
2
Since

rr

p
1
u


u
+
G = 0
2

r  (u  !) = ur
! !r
u
= u
r! + !
ru

(4.17)

u
r! + !
ru

(4.18)
(4.19)

the curl of the momentum equation gives us

@!
+ u
r! = !
ru +  r2 !
@t
60

(4.20)

n

'$
&%
2

I@@
@
  
!

'$
n&%
1

Figure 4.1: Portion of vortex tube.

The term !
ru represents stretching and tilting of the vortex lines.
For a two-dimensional ow in the x-y plane, the vorticity
!

is in the z -direction, so that

=  ez

@
+ u
r =  r2 
@t

(4.21)
(4.22)

4.4 Helmholtz's theorem

Many of the ideas regarding vorticity, vortex lines and vortex tubes are summarized in the following
theorems. A portion of a vortex tube is shown in Fig. 4.1.

Vortex lines move with the uid:

Consider a small area dA on the side of a vortex tube. Since the vorticity vector and the unit
normal to this area are perpendicular, the circulation around dA is zero. After an interval to
time, the uid forming this area has moved to a new position, but its circulation, by Kelvin's
theorem, is still zero. This can be said for all dA lying on the surface of the vortex tube so
that the uid elements forming the new tube also form a vortex tube. Thus vortex tubes move
with the uid. In the limit an innitesimal vortex tube is a line moving with the uid.
Circulation around a vortex tube is constant:

The sides of a vortex tube are vortex lines. Since

!

=ru

61

(4.23)

we have that

r
! =0

(4.24)

We integrate this over the volume of a vortex tube

Z

V
Using Gauss's theorem, we have

r
! dV = 0

(4.25)

!
n dA = 0
(4.26)
A
where n is the unit normal vector to the surface A. Since the integrand vanishes on the sides
of the vortex tube, we have that
Z

A1

!

n dA +
1

from which
1

A2
=

!

n dA = 0
2

(4.27)
(4.28)

A vortex tube cannot end with a uid; it must end at a boundary or form a closed loop:

Since the circulation is constant along a vortex tube, the tube cannot end within a uid.

Circulation around a vortex tube is constant in time:

This follows from Kelvin's theorem since the vortex tubes are made of material lines.

4.5 Bernoulli's theorem

For an inviscid uid with a conservative body force, the momentum equation is

We have the identity

@u
1
+ u
ru = rp + rG
@t


(4.29)

u
ru = 21 r (u
u) u  r  u

(4.30)

Furthermore, if the ow is barotropic

Z

r dp = 1 rp
Thus

@u
+r
@t

Z

(4.31)


dp 1
+ u
u G =u!
 2

There are two special cases.

62

(4.32)

4.5.1 Steady ow
Dening

B=

we see that

dp 1
+ u
u G
 2

(4.33)

rB = u  !

(4.34)

DB
@B
=
+ u
rB
Dt
@t
= 0

(4.35)

Since

B is a constant along a pathline which, for steady ow, is also a streamline. The constant may be
dierent for dierent streamlines.

4.5.2 Irrotational ow
Since ! = 0, we can take u = r. Thus, from equation (4.32), we have


+B =0
r @
@t

Thus,

(4.36)

@
dp 1
+
+ u
u G
@t
 2
is constant everywhere, though it may vary with time.

(4.37)

Problems

1. For an inviscid, barotropic uid with a conservative body force in a coordinate system rotating at a constant
rate
, show that Kelvin's theorem is
D
Dt

I

u
dl + 2

n dA = 0

where A is an area bounded by the closed curve C . [Hint: Start with u0 = u +

 r, where u0 and u are the
uid velocities in non-rotating and rotating frames, and r is the position vector in the rotating frame.]
2. For steady, inviscid, barotropic ow show that Bernoulli's equation holds along a vortex line.
3. Using the vorticity equation for an inviscid uid, show that if the vorticity of a uid particle is once zero it is
always zero.
4. Derive the vorticity equation for a compressible Newtonian uid with constant properties and with a conservative body force.

63

64

Chapter 5

Ideal ow
In this chapter we will discuss the motion of an incompressible, inviscid uid without vorticity.
For an inviscid uid we know that a uid particle, once irrotational, is always irrotational. Thus
the vorticity ! = r  u is zero everywhere for such a ow, so that

ru=0

(5.1)

u = r

(5.2)

r
u=0

(5.3)

everywhere. This is satised by

where  is called the velocity potential.

In addition let us consider the uid to be incompressible, so that the continuity equation is
from which

r =0

(5.4)
For an inviscid uid the normal velocity at a solid or impermeable boundary must be zero. Thus
2

@ 
=0
@n w

(5.5)

is the boundary condition for equation (5.4), where n is the coordinate normal to the wall.

5.1 Two-dimensional ows

Now, we will also, for simplicity, conne ourselves to plane, two-dimensional ow. Thus, we can
write the continuity equation as
@u @v
+ =0
(5.6)
@x @y
This is satised by
@
(5.7)
u =
@y
@
(5.8)
v =
@x
65

where (x; y) is the stream function.

In two dimensions the vorticity ! = !z k =  k, where

=

@v
@x

@u
@y

(5.9)

For an irrotational ow  = 0, so that we can write

@
@x
@
v =
@y

u =

(5.10)
(5.11)

which corresponds to equation (5.2).

The stream function and velocity potential both satisfy Laplace's equations

@2
@2
+ 2 = 0
2
@x
@y
2
@  @2
+
= 0
@x2 @y2

(5.12)
(5.13)

5.2 Properties
(a) Consider a line with

= constant. Along this line

0 = d
@
@
=
dx +
dy
@x
@y
= v dx + u dy

(5.14)
(5.15)
(5.16)

from which

dx dy
=
(5.17)
u
v
Since this is the equation of a streamline, it follows that the stream function is constant along a
stream line.
(b) Let Q be the volume ow rate per unit depth between two stream lines with stream functions
1 and
2 , respectively. Then

dQ = u dy v dx
@
@
=
dy +
dx
@y
@x
= d
so that we have

Q=

66

(5.18)
(5.19)
(5.20)
(5.21)

(c) The

= constant lines have a slope



dy
dx

v
u

(5.22)

On the other hand, the  = constant lines have a slope



from which
The constant

dy
u
=
dx 
v

 

dy
dy
= 1
dx
dx 
and constant  lines are thus orthogonal to each other.

(5.23)
(5.24)

5.3 Complex representation

Comparing equations (5.10) and (5.11) with (5.7) and (5.8), we have

@
@
=
@x
@y
@
@
=
@y
@y
These are the Cauchy-Riemann conditions for the complex function
F (z ) = (x; y) + i (x; y)

(5.25)
(5.26)
(5.27)

to be analytic everywhere. z = x + iy is the position coordinate, and F (z ) is the complex potential.

The derivative of the complex potential is

dF
(5.28)
dz
where W is referred to as the complex velocity. Since F (z ) is an analytic function, we can take its
derivative in any direction. Choosing the x-direction, we get
@ @
W (z ) =
+i
@x @x
= u iv
(5.29)
W (z ) =

Furthermore

W W = u2 + v2

(5.30)

5.4 Polar form

We will frequently have to use polar coordinates (r; ), and their corresponding velocity components
(ur ; u ) as shown in Fig. 5.1. The velocity vector can be written as

u = ur er + u e
67

(5.31)

u

I@@ 

ur

r


Figure 5.1: Velocity components ur and u in polar coordinates.

The incompressibility condition is
1@
1 @u
(ru ) +
(ru ) = 0
r @r r r @ r
so that we can dene the stream function as
1@
ur =
r @
@
u =
@r
The irrotationality condition is

1 @
(ru )
 =
r @r 
= 0
from which

@ur
@

@
@r
1 @
u =
r @
ur =

Laplace's equations for

and  are

@ @
1 @2
r
+
= 0
@r @r
r @2


@
@
1 @2
r
+
= 0
@r @r
r @2

(5.32)
(5.33)
(5.34)
(5.35)
(5.36)
(5.37)
(5.38)
(5.39)
(5.40)
(5.41)
(5.42)

Since F =  + i , we have

W =

dF
dz

(5.43)
68

= e

@ @
+i
@r
@r
i
iu )e

i

= (ur

(5.44)
(5.45)

Referring to Fig.5.1 the relation between the Cartesian and polar components of the velocity is
found to be

u = ur cos  u sin 
v = ur sin  + u cos 

(5.46)
(5.47)

ur = u cos  + v sin 
u = u sin  + v cos 

(5.48)
(5.49)

and

The polar unit vectors can be written in terms of the Cartesian unit vectors as

er
e

= cos i + sin j
=
sin i + cos j

(5.50)
(5.51)

Notice that the polar unit vectors depend on . Thus

@ er
=
@
=
@ e
=
@
=

e

sin i + cos j

(5.52)
(5.53)

cos i sin j

(5.54)
(5.55)

er
69

5.5 Summary of equations

Cartesian

Incompressibility

@u
@x

Irrotationality

@v
@x

Velocity potential
Stream function
Cauchy-Riemann eqns.

@v = 0
+ @y
@u
@y

=0

u = @
@x
v = @
@y

Polar
@
@u 
r @r (rur ) + @ = 0

@
r @r (ru )

@ur  = 0
@

ur = @
@r
u = r1 @
@

u=
v=

@
@y
@
@x

ur =
u =

@ =
@x
@ =
@y

@
@y
@
@x

@ = 1 @
@r
r @
@
1 @
=
r @
@r

@
r @
@
@r
1

2
2
@ (r @ ) + 1 @ 2  = 0
Laplace's eqn. @@x2 + @@y2 = 0 1r @r
@r
r2 @2
@ 2 + @ 2 = 0 1 @ (r @ ) + 1 @ 2 = 0
@x2
@y2
r @r @r
r2 @2

Complex potential

F =+i

Complex velocity

W = u iv

F =+i
W = (ur

iu )e i

5.6 Simple ows

5.6.1 Uniform ow
F =
=
 =
W =
u =
v =

Uz
Uy
Ux
U
U
0

(5.56)
(5.57)
(5.58)
(5.59)
(5.60)
(5.61)

5.6.2 Source or sink

m
ln z
2
m
=

2
m
 =
ln r
2

F =

70

(5.62)
(5.63)
(5.64)

m
2z
m
ur =
2
u = 0

W =

(5.65)
(5.66)
(5.67)

5.6.3 Vortex
F =
=

 =
W =
ur =
u =

ln z
2
ln r
2

2
i
2z
0
2r

(5.68)
(5.69)
(5.70)
(5.71)
(5.72)
(5.73)

5.6.4 Sector with angle =n

F =
=
 =
W =
ur =
u =

Uz n
Urn sin n
Urn cos n
nUz n 1
nUrn 1 cos n
nUrn 1 sin n

(5.74)
(5.75)
(5.76)
(5.77)
(5.78)
(5.79)

5.7 Combined ows

Since the Laplace's equation is linear, solutions may be added together to produce other solutions.

5.7.1 Doublet
F =
=

 =
W =
ur =
u =


z

sin 

r
cos 

r

z2
cos 
 2
r
sin 
 2
r
71

(5.80)
(5.81)
(5.82)
(5.83)
(5.84)
(5.85)

5.7.2 Cylinder without circulation



a2
z2


a2
Ur 1 2 sin 
r


a2
Ur 1 + 2 cos 
r


a2
U 1 2
z


a2
U 1 2 cos 
r


a2
U 1 + 2 sin 
r

F = Uz 1 +
=

 =
W =
ur =
u =

(5.86)
(5.87)
(5.88)
(5.89)
(5.90)
(5.91)

5.7.3 Cylinder with circulation



a2
z
+ i ln
2
z
2 a


2
a
r
Ur 1 2 sin  + ln
r
2 a


a2

Ur 1 + 2 cos 
r
2


a2
U 1 2 +i
z
2z


2
a
U 1 2 cos 
r


a2
U 1 + 2 sin 
r
2a

F = Uz 1 +
=

 =
W =
ur =
u =

(5.92)
(5.93)
(5.94)
(5.95)
(5.96)
(5.97)

5.8 Forces on a submerged body

For steady, incompressible, irrotational ow we can use Bernoulli's equation

p 1
+ u
u=B
 2

(5.98)

to determine the pressure, where B is a constant.

On the surface of a body the only force is due to pressure. The total force on the body is given
by
Z
F = pn dA
(5.99)
A
where dA is an elemental area of the surface A with normal n.
72

5.8.1 Cylinder with circulation

The pressure on the surface of the cylinder is given by


2
1
p = p1 + U 2  2U sin  +
(5.100)
2
2a
where p1 and U are the pressure and velocity of the ow far from the cylinder. The components of
the force on the cylinder are

Fx =

al

= 0

Fy =

al

Z 2
0

Z 2
0

p cos d
(5.101)

p sin d

= U l
where l is the length of the cylinder. This is the Kutta-Joukowsky theorem.

(5.102)

5.9 Conformal transformation

An analytic function  = f (z ) that maps the z -plane to the  -plane is called a conformal transformation. Let z = x + iy and  =  + i. Then it can be shown that if
@2 @2
+
=0
(5.103)
@x2 @y2
then
@2 @2
+
=0
(5.104)
@ 2 @2
also. We can also show that ratios of length and angles are preserved in the transformation. The
complex velocity in the two planes is related by
dF
W (z ) =
dz
d dF
=
dz d
d
W ( )
(5.105)
=
dz

5.9.1 Joukowski transformation

This is

c2

where c is a constant. The function is analytic for  6= 0. The reverse mapping
z=+
r

(5.106)

z2 2
z
= 
c
(5.107)
2
4
is non-unique. As an example the Joukowski transformation can be applied to a uniform ow to
obtain ow around a cylinder.
73

z -plane

 -plane








Figure 5.2: Mapping from a circle to an ellipse

Flow around an ellipse

Let c be real and positive. Then any point on a circle with center at the origin, of radius a where
a < c, in the  -plane can be written as  = aei , or
 = a cos 
(5.108)
 = a sin 
(5.109)
Substituting in the transformation and separating the real and imaginary parts, we have


c2
x = a+
cos 
a


c2
y = a
sin 
a
Eliminating  , we have
2 
2

y
x
+
=1
a + c2 =a
a c2 =a
which is the equation for an ellipse. Figure 5.2 shows the circle and its elliptic image.
The complex potential for a uniform ow at an angle  around a cylinder of radius a is


a2 i
i
F ( ) = U e + e

Substituting equation (5.107), we get the complex function
2

F (z ) = U 4

z

2

z2
4

c2 e i + @

74

a2

z2
4

(5.110)
(5.111)
(5.112)

(5.113)

A ei 5

(5.114)

for ow around an ellipse.

Flow around a at plate

For a body with sharp trailing edge, the rear stagnation point is at the trailing edge. This is the
Kutta condition. Consider a at plate 2a  x  2a, where a = c. Add a circulation of strength
= 4Ua sin 
and complex potential

F=

i
z
ln
2 a

(5.115)
(5.116)

to satisfy the Kutta condition. Thus




a2
F ( ) = U e i + ei + i2Ua sin  ln

a

(5.117)

Substituting equation (5.107), we get the complex potential in the z -plane.

Determing the force on the at plate, we nd that

Fy = U l
= 4U 2 al sin 

(5.118)

Interms of the area A = 4al, the coecient of lift is

Fy
U 2 A
= 2 sin 

CL =

1
2

(5.119)

Flow around a symmetrical Joukowski airfoil

This is obtained by a transformation of a circle that has its center on the real axis but displaced
from the origin, as shown in Fig. 5.3.

Flow around a circular-arc airfoil

The circle is now on the imaginary axis, but not at the origin as shown in Fig. 5.4.

Flow around a Joukowski airfoil

This is shown in Fig. 5.5.

5.10 Three-dimensional axisymmetric ow

Either cylindrical or spherical coordinates may be used to deal with axisymmetric problems.
75

 -plane





Figure 5.3: Mapping from a circle to a symmetrical airfoil

5.10.1 Cylindrical coordinates

5.10.2 Spherical coordinates
For axisymmetric ow there is zero velocity and dependence in the '-drection. The vorticity vector
is thus given by


@ur
1 @
!=
(ru )
e
(5.120)
r @r 
@ '
For irrotational ow this is zero, so that the velocity potential is dened by

@
@r
1 @
u =
r @
ur =

(5.121)
(5.122)

The mass conservation equation is

1 @
1 @ 2
(r ur ) +
(u sin ) = 0
2
r @r
r sin  @ 

(5.123)

This can be satised by the Stokes stream function

1 @ s
r2 sin  @
1 @ s
u =
r sin  @r

ur =

76

(5.124)
(5.125)

 -plane




Figure 5.4: Mapping from a circle to a circular arc

Uniform ow
 = Ur cos 
1 2 2
Ur sin 
s =
2
ur = U cos 
u = U sin 

Source and sink

m
4r
m
(1 + cos )
=
s
2
m
ur =
4r2
u = 0
 =

Doublet
 =
s

ur =
u =


cos 
4r

sin2 
4r

cos 
2r3

sin2 
2r3
77

(5.126)
(5.127)
(5.128)
(5.129)

(5.130)
(5.131)
(5.132)
(5.133)

(5.134)
(5.135)
(5.136)
(5.137)

 -plane





Figure 5.5: Mapping from a circle to a Joukowski airfoil

Flow around a sphere



a3
cos 
2r3


1 2
a3
Ur 1 3 sin2 
s =
2
r


a3
ur = U 1 3 cos 
r


a3
u = U 1 + 3 sin 
2r
 = Ur 1 +

Problems

(5.138)
(5.139)
(5.140)
(5.141)

1. Show that the complex potential, F (z), that satises z = c cosh F represents ow through an aperture.
2. Computer generate plots of the ow around a symmetrical Joukowsky airfoil.
3. If F 2 = U 2 (z2 + c2 ), show that the stream function satises


2 U 2 (x2 + c2 ) + 2
2
y =
U 2 (U 2 x2 + 2 )
and that it represents a stream of velocity U past a thin obstacle of length c projecting perpendicularly from
a straight boundary.
4. State and prove the circle theorem.
5. A source is symmetrically placed near a corner. Find the force on the walls due to potential ow from the
source.
6. Find the velocity and pressure elds for ow over a wedge of angle .
78

Chapter 6

Incompressible viscous ow: exact

solutions
6.1 Flow between at plates
Consider uid between two innite at plates separated by a distance L. The coordinate system is
as shown. The lower plate is stationary while the upper one is moving with a velocity U parallel
to itself. The plates are also maintained at dierent temperatures. Neglect gravity and assume a
steady ow in the x-direction, with u = u(y); v = w = 0; T = T (y). Of the governing equations,
equation (A.1) is identically satised. Equations (A.3) and (A.4) become

@p
(6.1)
@y
@p
(6.2)
0=
@z
from which we know that the pressure is a function of x alone. Let us assume that a pressure
gradient is imposed externally so that dp=dx is a known constant. Equation (A.2) becomes
0=

0=

dp
d2 u
+ 2
dx
dy

Equation (A.5) simplies to

d2 T
du
0=k 2 +
dy
dy

(6.3)

!2

(6.4)

6.1.1 Alternative derivation of governing equations

Consider the forces in the x-direction on an element of uid of dimension dx  dy as shown. They
are due to pressure on each one of the vertical faces and to shear stress at the horizontal faces. These
forces will be calculated per unit length in the z -direction. The net force due to pressure is

Fp = dp dy =
79

dp
dx dy
dx

(6.5)

U-

T = T1

T = T0

Figure 6.1: Flow between at plates

 (y + dy)

p dy

- dy

(p + dp) dy

dx

 (y)

Figure 6.2: Force balance in ow between at plates

The net force due to shear stress is

"

F =  (y + dy)  (y) dx =

d
dy dx
dy

(6.6)

Since the uid element is moving at constant velocity and not accelerating, the sum of the forces on
it must be zero. Thus
Fp + F = 0
(6.7)
from which, on dividing by dx dy, we get the momentum equation
dp d
+ =0
(6.8)
dx dy
The constitutive relation for a Newtonian uid, equation (3.61), gives the stress eld as
2
4

p  du
dy
 du
p
dy
0
0

Thus, with the present notation

 =
80

du
dy

0
0
p

3
5

(6.9)

W (y + dy)

Q(y + dy)

Wp (x)

- dy


dx

Wp (x + dx)

W (y)
Q(y)
Figure 6.3: Energy balance in ow between at plates
so that

d2 u
dp
+ 2 =0
dx
dy

(6.10)

This is the same as equation (6.3).

Let us now consider the rate at which energy is owing into the uid element, again per unit
length in the z -direction.
There is no heat conduction through the vertical faces since the temperature does not vary in
the x-direction. The pressure, however, does work on the uid element.

Wp (x) = pu dy
"

Wp (x + dx) =

p+

(6.11)

dp
dx u dy
dx

(6.12)

remembering that u does not vary in the x-direction and that the pressure force and velocity are in
opposite directions for Wp (x + dx). W (y) and W (y + dy) are the work done by the shear forces.
Referring to Fig. 6.3, we have

W (y) = u dx
"

W (y + dy) =

u +

(6.13)

d(u)
dy dx
dy

(6.14)

The heat uxes are Q(y) and Q(y + dy) where

Q(y) =

k
"

Q(y + dy) =

dT
dx
dy
dT d
k +
dy dy

(6.15)

dT
k
dy

dy dx

(6.16)

Summing up, conservation of energy gives the relation

Q(y) Q(y + dy) + Wp (x) + Wp (x + dx) + W (y) + W (y + dy) = 0

81

(6.17)

Substituting for the dierent terms and dividing by dx dy, we get

d2 T
dy2

dp d(u)
+
=0
dx
dy

(6.18)

This is the total energy equation. If we subtract u times the momentum equation (6.1.1) from this,
we have
d2 T
du
k 2 + =0
(6.19)
dy
dy
which is the thermal energy equation. Using equation (6.9) for the shear stress, we obtain

du
d2 T
k 2 +
dy
dy

!2

=0

(6.20)

which is equation (6.4).

6.1.2 Velocity and temperature proles

The solution to equation (6.3) is
1 dp y2
u(y) =
+ Ay + B
 dx 2

(6.21)

where A and B are the constants of integration. The boundary conditions are

u = 0 at y = 0
u = U at y = L

(6.22)
(6.23)

The rst condition gives B = 0. From the second we have

A=
The velocity prole is

y
u(y) = U
L

U
L dp=dx

L
2

L2 dp y
1
2 dx L

Nondimensionally

(6.24)

y
L

u =  + P (1 )
where the dimensionless velocity, distance and pressure gradient are
u
u =
U
y
 =
L
!
L2
dp
P =
2U
dx
82

(6.25)
(6.26)
(6.27)
(6.28)
(6.29)

respectively. The velocity of the upper plate U can be used to scale the uid velocity only if it
is nonzero. Otherwise, it can be some other characteristic velocity, like the mean ow velocity for
instance.
Substituting the velocity prole from equation (6.25) into the energy equation (6.4), we have

d2 T

=
dy2
k

"

U
L

!2

U dp
y
L2 dp
1 2
+ 2
 dx
L
4 dx

!2

y
1 2
L

!2 #

(6.30)

the solution to which is

"

T = A + By

 1 Uy
k 2 L

!2

Uy2 dp
1
2 dx

2y
L2 y2 dp
+ 2
3L
8 dx

!2

4y 2y 2
+
3L 3L2

!#

(6.31)

where A and B are constants. Introducing the boundary conditions

T = T0 at y = 0
T = T1 at y = L

(6.32)
(6.33)

we can evaluate A and B and get the temperature prole

T = T0 + (T1

U 2 y
y
+
1
T0 )
L
2k L
L4 dp
+
24k dx

y
L

!2

UL2 p y
y
y2
1 3 +2 2
6k x L
L
L

y
y
y2
1 3 +4 2
L
L
L

In terms of nondimensional temperature T  = (T

can be written as

T0)=(T1

y3
2 3
L

(6.34)

T0 ) and distance  = y=L, this

"

1
1
1
T  =  + Br (1 ) P (1 3 + 22 ) + P 2 (1 3 + 42
2
3
6

23 )

(6.35)

where the Brinkman number Br = Ec P r, with

Eckert number Ec =

U2

c(T1 T0 )
c
Prandtl number P r =
k
The temperature dierence T1 T0 can be used as a scale only if it is nonzero.

(6.36)
(6.37)

6.1.3 Couette ow
If no pressure gradient is imposed on the ow, P = 0 and the dimensionless velocity and temperature
proles become
u = 
(6.38)
1
(6.39)
T  =  + Ec P r (1 )
2
83

The maximum temperature occurs within the ow at  = m where

1
1
m = +
2 Ec P r
if 0    1. Then the maximum temperature is
1
1
1
Tm = 1 +
+ Ec P r
2
Ec P r 8

(6.40)

(6.41)

6.1.4 Poiseuille ow
On the other hand, if the velocity of the upper plate is zero, the ow is driven solely by the pressure
gradient, U = 0. The velocity prole is then

L2 dp y
u(y) =
1
2 dx L

y
L

(6.42)

The maximum velocity is at the center y = L=2, so that

L2 dp
(6.43)
8 dx
The mean velocity U is dened by the volume ow rate per unit area. In this case it is given by
Z
1 L
U =
u(y) dy
(6.44)
L 0
L2 dp
=
(6.45)
12 dx
2
= umax
(6.46)
3
Using this to get a nondimensional velocity u , we have
u
u = 
(6.47)
U
= 6(1 )
(6.48)
umax =

6.1.5 Heat generation

The heat ux through either of the plates is k dT=dy per unit area. For simplicity, let us consider
the Couette ow problem with zero dp=dx. We can write a heat balance for the control volume
shown. Since
!
dT T1 T0 U 2
y
=
+
1 2
(6.49)
dy
L
2kL
L
we have

Q0 =
Q1 =

kA
kA

T1
T1
84

L
L

T0

U 2
+
2kL

T0

U 2
2kL

(6.50)
!

(6.51)

Q1

?6
6

Q0

Figure 6.4: Heat generation in ow between at plates

where A is the cross sectional area of the control volume. So the heat generated within the control
volume is

QG = Q1 Q0
(6.52)
U 2 A
(6.53)
=
L
This is due to the work done on the uid by the shear stress. There is no pressure work; also there
is no work done at the lower plate due to shear stress since the uid velocity there is zero. The rate
of work done at the upper plate is
W =
=



AU 


(6.54)

U 2 A
L

(6.55)

y=h

which is converted into heat.

6.2 Flow between coaxial rotating cylinders

Fluid is contained within two innitely long cylinders as shown in Fig. 6.5. The inner and outer
cylinders have radii r0 and r1 , rotate counterclockwise at angular speeds of !0 and !1 , and are kept at
temperatures T0 and T1 , respectively. We use cylindrical coordinates with the velocity components
ur = uz = 0; u = u (r), the pressure p = p(r), and the temperature T = T (r). Gravity is neglected.
The continuity equation (A.6) is identically satised. Equations (A.6) - (A.6) give

u2
r

dp
dr

(6.56)

"

1 @ @u u
0 = 
r
r @r @r r2
"

1 @ @T
r
0 = k
r @r @r
85

!#

(6.57)
"

@ u
+ r
@r r

!#2

(6.58)

!0

r1

!1

r0

Figure 6.5: Flow between rotating cylinders

Equation (6.57) can be simplied to
!

d2 u d u
+
=0
dr2 dr r
On integrating twice, we have

(6.59)

B
r
where A and B are constants. The boundary conditions are
u = Ar +

(6.60)

u = !0 r0 at r = r0
u = !1 r1 at r = r1

(6.61)
(6.62)

The constants can be determined to give the velocity prole

u =

r12

"

r02

(! r

2 2
1 1

! r )r (!1
2 2
0 0

r2 r2
!0 ) 0 1
r

(6.63)

The pressure distribution from equation (6.56) is

p = p0 +

(r

2
1

2(!1

"

r)

2 2
0

(!12 r12

!0 )(! r

!02 r02 )2

r2
2

! r ) ln r (!1

2 2
1 1

2 2
0 0

where p0 is a constant of integration.

The temperature distribution

r02
r2

T T0
r4 (1 ! =! )
= Br 1 4 1 4 0 1
T1 T0
r1 r0
86

!"

(6.64)

rr
!0 )
4r

4 4
2 0 1

(6.65)

ln(r=r0 )
ln(r=r0 )
+
ln(r1 =r0 )
ln(r1 =r0 )

(6.66)

T = T0

Figure 6.6: Flow in a circular pipe

is obtained by solving equation (6.58). The Brinkman number is given by

Br =

r02 !02
k(T1 T0 )

(6.67)

6.3 Flow in a circular pipe

Flow is through a long pipe with a pressure gradient dp=dz as in Fig. 6.6. We neglect gravity and
let uz = uz (r); ur = u = 0. The continuity equation (A.6) is identically satised. Equations
(A.6){(A.6) reduce to

@p
@r
1 @p
r @q

0 =
0 =

(6.68)
(6.69)

so that p = p(z ) only. Equations (B.7) and (B.8) become

"

1 d duz
dp
+
r
0=
dz
r dr dr

!#

(6.70)

The solution of this equation is

R2
uz =


dp
dz

!"

1 r
4 R

!2

r
+ A ln
+B
R

(6.71)

At the centerline the velocity must be nite, so that A = 0. The other condition is that

uz = 0 at r = R

(6.72)

After calculating B , we obtain the velocity prole as

1
uz =
4

dp
dz
87

(6.73)

The maximum velocity Umax is at the centerline r = 0, where

R2
Umax =
4

dp
dz

(6.74)

The mean velocity is

Z

R
1
u 2r dr
R2 0 z
!
dp
R2
=
8
dz
1
= Umax
2

U =

so that

r2
R2

uz = 2U 1

Equation (6.71) can be written as

"

1 @ @T
k
r
r @r @r

!#

(6.75)
(6.76)
(6.77)

(6.78)

16U 2r2
R4

(6.79)

6.3.1 Isothermal wall

The wall is at a constant temperature T0 . Since T = T (r), the energy equation becomes
"

1 d dT
r
0=k
r dr dr

!#

duz
+
dr

!2

(6.80)

Integrating, and using the conditions that T is nite at r = 0, and T = T0 at r = R, we get

U 2
1
T = T0 +
k

r4
R4

(6.81)

6.4 Flow over a porous wall

A uniform ow of velocity U exists over a porous wall as shown in Fig. 6.7. There is suction
through the porous wall resulting in a normal velocity of V . The pressure p is constant. The twodimensional velocity eld is of the form u = u(y). From the continuity equation @v=@y = 0, so that
V is independent of y. Since v = V at y = 0, it is v = V everywhere. The x-momentum equation
is
du
d2 u
V =  2
(6.82)
dy
dy
with boundary conditions

u = 0 at y = 0
u = U at y ! 1
88

(6.83)
(6.84)

-U

6-x
?

Figure 6.7: Flow over a porous wall

The solution is
where the kinematic viscosity  = =.

u = U (1 e V y= )

(6.85)

6.5 Natural convection between vertical at plates

In natural convection, body forces due to changes in density are taken into account. For a small
change in the temperature, the following linear relationship may be used
"

 = 0 1  (T

T0 )

(6.86)

where  is the coecient of volumetric expansion, and 0 is the density at the reference temperature
T0 . Under the Boussinesq approximation we keep the properties of the uid constant, but use a
variable density only in the body force term of the momentum equation. Viscous dissipation may
usually be neglected.
A simple solution is obtained for the case of fully developed ow between at plates at dierent
temperatures as shown in Fig. 6.8. The reference temperature is taken to be the average of the two.
The pressure gradient is hydrostatic, i.e. dp=dx = 0 g.
For u = u(y); v = w = 0; T = T (y), the governing equations are
d2 u
0 =  2 + 0 g (T T0)
(6.87)
dy
d2 T
0 = k 2
(6.88)
dy
with boundary conditions
u = 0; T = T0
T at y = L
(6.89)
u = 0; T = T0 +
T at = L
(6.90)
Solutions are
!
0 g
T L2 y y2
u =
1
(6.91)
6
L L2
y
(6.92)
T = T0 +
T
L
89

gravity
vector

?

T0

x
y

2h

T

T0 +
T

Figure 6.8: Natural convection between parallel walls

Nondimensionally, if

0 Lu

T T0
T =

T
y
 =
L
u =

then

Gr 2
u =
(
6
T = 

(6.93)
(6.94)
(6.95)
1)

(6.96)
(6.97)

where the Grashof number is Gr = g
T L3= 2 .

Problems

1. A 5 cm diameter shaft rotates at2 3000 rpm within a bearing of length 5 cm. The clearance is 0.1 mm and lled
with oil of viscosity 0.01 N s/m . Approximating the ow of oil to be that between at plates, nd the rate of
heat generation within the bearing.
2. For Poiseuille ow between at plates, nd the principal axes of the stress tensor at one of the plates.
3. For Couette ow between at plates, we can dene the Nusselt number as Nu = hc L=k, where qwall =
hc (T1 T0 ). Find the Nusselt numbers at the walls in terms of the Brinkman number.
4. Find the stress and strain rate tensors for the problem of ow between rotating cylinders.
5. Fluid of viscosity  is contained in the small gap of width , between two cones. Find the torque necessary to
rotate the inner cone at an angular velocity
, keeping the outer one stationary. Find also the heat generated
within the uid. Use a at plate approximation for the ow in the gap.
90

@@@
@@@
@@
@@@
@@
@@
@@@

@R

I@

Figure 6.9: Flow between cones

6. An incompressible Newtonian liquid ows down an innite vertical plane in a lm of constant thickness under
the action of gravity. Determine the velocity prole, and write the simplied dierential equation governing
the temperature prole.
7. Using a force balance on the elemental volume shown, derive equation (6.70).
8. Two porous plates are separated by a distance 2L; the ow through each plate is as indicated. The ow is
driven by a constant pressure gradient dp=dx. Find the velocity eld.
9. Show that the volume ow rate Q for laminar ow in a circular pipe is given by
Q=

R4
8

dp
dx

10. The ow rate of glycerine through

the apparatus shown is measured to be 0.66 cm3/min. Given that the
density of glycerine is 1260 kg/m3 , determine its viscosity.
11. Determine the wall shear stress w and the coecient of friction corresponding to this, Cf = w = 21 U 2 , for
ow over an innite porous wall with suction.
12. SAE 30 motor oil (thermal conductivity = 0.145 W/m K) occupies the space between two parallel plates which
are kept at 20 C and 2 mm apart. One plate is stationary, while the other is moving parallel to itself at 20
m/s. What is the maximum temperature in the oil?
13. A 5 cm diameter shaft rotates at 3000 rpm within a bearing of length 5 cm. The clearance is 0.1 mm and
lled with SAE 30 motor oil (thermal conductivity = 0.145 W/m K). Approximating the ow of oil to be that
between at plates, nd the maximum temperature in the oil if the shaft and bearing are both kept at 20 C.
14. Show, by taking the curl of the two-dimensional momentum equation for an incompressible Newtonian uid,
that
D
=  r2 
Dt
where  is the magnitude of the vorticity normal to the plane of the ow.
15. Determine the velocity and temperature proles in the pressure-driven ow in the annular space between two
long cylinders.
16. A viscous, incompressible liquid between two innite, parallel plates falls due to gravity alone. Find the velocity
prole.
91

17. Two vertical, porous plates are separated by a distance 2L; the ow through each plate, U , is to the right.
The left and right plates have temperatures T0
T and T0 +
T , respectively. The ow is driven by natural
convection. Neglect viscous dissipation.
(a) Find the simplied equations for the velocity and temperature elds.
(b) Solve the temperature eld.
(c) Solve the velocity eld in terms of constants of integration. Indicate how you would determine the constants.
18. An incompressible, viscous uid occupies the annular space between two innite concentric cylinders. The
inner cylinder of radius r1 is pulled in an axial direction with constant velocity U , while the outer one of radius
r2 is stationary. Find the velocity prole in the uid.
19. Find the velocity eld in a layer of liquid of constant thickness owing down an inclined plane due to gravity.
20. Find the heat dissipation rate for a journal bearing carrying a 150 mm diameter shaft rotating at 720 rpm.
The bearing is 200 mm long and there is an average20.15 mm clearance between the shaft and the bearing.
Lubrication is by turbine oil of viscosity 0.0012 Ns/m . Assume Couette ow in the clearance.
21. Find the fully developed velocity eld in an axisymmetric lm of liquid owing down due to gravity outside a
vertical, circular rod. The radius of the rod is a, and the lm thickness is .

92

Chapter 7

Incompressible viscous ow:

negligible inertia
7.1 Stokes's ow
At low Reynolds numbers, the inertia forces may be neglected. Under steady ow conditions then,
we have
div u = 0
grad p + r2 u = 0

(7.1)
(7.2)

This is a linear set of equations. On taking the divergence of the momentum equation, and using
the mass conservation equation, we have
r2 p = 0
(7.3)
On taking the curl, however, we have

r !=0

(7.4)

since the curl of any gradient is zero, and ! = curl u.

7.2 Uniform ow past a sphere

The governing equations (7.1) and (7.2) can be solved analytically for steady uniform ow around
a sphere, as shown in Fig. 7.1. In a mixture of spherical and Cartesian coordinates


3 ax2 a2
1a
a2
1
3+ 2 +1
3
2
4 r
r
4r
r
 2

3 axy a
v = U 3
1
4 r
r2


3 axy a2
1
w = U 3
4 r
r2
3 Uax
p p0 =
2 r3

u = U

93

(7.5)
(7.6)
(7.7)
(7.8)

a
U

Figure 7.1: Flow past a sphere

hhhhhhh
y

6


hhhhhhh
hhh
6
h(x)

-x

-U

Figure 7.2: Lubrication ow

where

r 2 = x2 + y 2 + z 2
(7.9)
This velocity eld is zero at r = a, and becomes u = U; v = w = 0 at r ! 1.
The drag due to the pressure is 2aU , and that due to viscous stress is 4aU . The total drag
force Fd on the sphere is thus
Fd = 6aU
(7.10)
This result is valid only for small Reynolds number Re = U (2a)= < 1.

7.3 Lubrication theory

Fluid is contained within the space between two surfaces sliding against each other as in Fig. 7.2.
The lower surface has a horizontal motion of velocity U ; the upper one which is slightly inclined is
stationary. The ow is mostly in the x-direction so that v  u. Since w = 0 from two-dimensionality,
the y- and z -momentum equations are negligible compared to the x-momentum equation. Thus,
p = p(x). The viscous terms in the x-momentum equation are of order @ 2 u=@y2  U=h2 , the
other viscous terms being smaller since @ 2 u=@x2  @ 2 u=@y2. The inertia term is of order u@u=@x 
U 2 =a. The ratio, called a reduced Reynolds number Re = (Ua=)(h2 =a2 ), is normally small in
lubrication problems so that inertia can be neglected altogether.
So the x-momentum equation becomes

@2u
dp
= 2
dx
@y
94

(7.11)

Since the gradient dp=dx is a function of x, this can be integrated twice to give
1 dp 2
y + A(x)y + B (x)
2 dx

u=

(7.12)

The boundary conditions to be applied are

y = 0; u = U
y = h(x); u = 0
so that we get

y
h(x)

u=U 1

(7.13)
(7.14)

1 dp
y (h(x) y)
2 dx

(7.15)

The continuity equation can be written as

Q=

Z h(x)

u dy

(7.16)

where Q, the volume ow rate per unit length, is a constant. From this

Q=

h
U
2

h2 dp
6 dx

(7.17)

At both ends, x = 0 and x = a, the pressure must be atmospheric. Integrating this equation, and
using the rst condition at x = 0, we have

p = p0 + 6U

Z x

Z x
1
1
dx
12
Q
dx
2
3
h(x)
0 h(x)

(7.18)

The condition at the other end x = a gives

1
Q = UH
2

(7.19)

where the characteristic thickness H is given by

Ra

1
2 dx
H = R0a h(x1)
dx
0 h(x)3

(7.20)

Per unit length the total force on the lower surface has the two components

Fx =
Fy =

Z a
Z a
0

@u

@y

p(x) dx
95

y=0

dx

(7.21)
(7.22)

7.4 Flow in porous media

The continuity equation for incompressible ow in a porous medium is
div u = 0

(7.23)

For the momentum equation, the simplest model is that due to Darcy in which the ow velocity is
taken proportional to the imposed pressure gradient. Thus

K
(grad p f )
(7.24)

Here K is called the permeability of the medium. Thus, in the incompressible Navier-Stokes equation
with constant properties, the inertia terms are dropped and the viscous force per unit volume
is represented by (=K )u. The condition on the velocity is that of zero normal velocity at a
boundary, allowing for slip in the tangential direction.
From equations (7.23) and (7.24), we get

u=

r p=0
2

from which the pressure distribution can be determined.

The energy equation is


@T
f cf  + u
grad T = ke r2 T
@t
where ke is the eective thermal conductivity of the porous medium, and

=

f cf + (1 )s cs
f cf

(7.25)

(7.26)

(7.27)

Subscripts f and s refer to the uid and solid respectively, and  is the porosity of the material.

Problems

1. A 2 mm diameter stainless steel ball (density = 8000 kg/m3 ) is dropped in glycerine (viscosity = 0.8 Ns/m2).
Find the terminal velocity of the ball, i.e. the velocity reached after the initial transient has died down.
2. In the lubrication problem, nd the components of the force on the moving lower plate if the gap is h(x) = b cx
for 0  x  a, where a < b=c, b > 0, c > 0.
3. For Stokes
ow around a sphere, show that the drag coecient is given by Cd = 24=Re, where Cd =
Fd =(R2 21 U 2 ); Re = 2UR=.
4. Find the
terminal velocity of fall in air of a 10 mm diameter aluminum particle. Density of aluminum = 2700
kg/m3 , density of air = 1.2 kg/m3 , viscosity of air = 1:9  10 5 Ns/m2. Check that the Reynolds number is
small enough for Stokes ow.
5. Show that the velocity and pressure elds given by equation (7.5){(7.8) are solutions of the governing equations
(7.1) and (7.2).
6. In the lubrication problem, show that the special caseRofa a wedge-like shape with h(x) = b(c x), gives
H = 2bc(c a)=(2c a). Find the force per unit length 0 p(x) dx.
7. If the entrance and exit gap widths for the previous problem are h1 and h2 respectively, show that H =
2h1 h2 =(h1 + h2 ). Find the pressure distribution, the pressure and the viscous forces.
8. A glass-ber porous lter
is placed within a duct with a square section 30 cm  30 cm. If the permeability of
the lter is5 0:5  102 9 m2 , nd the pressure drop for 0.2 m3/s of air ow through the lter. Viscosity of air =
1:9  10 N s/m .
96

9. Find the radial volume ow rate of uid through an annular porous material. The inner and outer pressures are
known. Find also the temperature distribution if the inner and outer surfaces are kept at dierent temperatures.
10. Show that r2 p = 0 for ow of an incompressible uid through a porous medium.
11. A viscous uid occupies the space between an upper plate of shape h(x) = a=(1 + bx) and a at lower
plate
moving with velocity U , where a, b and U are constants. Given that U = 5 m/s,  = 0.01 Ns/m2 and the
distances h1 = 2 mm, h2 = 1 mm, L = 10 cm, nd (a) the characteristic thickness, and (b) the pressure
distribution.
12. A cylindrical, compressed-air tank (inner radius a, outer radius b, length L) leaks air (viscosity , density
) through its cylindrical, porous walls to the atmosphere. If the permeability of the wall is K , nd the air
pressure, p(t), in the tank as a function of time. Assume an ideal gas law for the air and isothermal conditions.
13. Estimate the net pressure force on a shaft rotating at speed ! in a journal bearing. The shaft radius is a, the
bearing radius is a + , and the distance between the shaft and bearing centers is . Assume that =a  1, and
make any related approximations.

97

98

Chapter 8

Incompressible viscous ow: at

plate boundary layer
8.1 Prandtl boundary layer equations
For high Reynolds number, the entire ow region can be divided into two parts for hydrodynamic
analysis: an outer region where the eect of viscosity is not important, and an inner region near
boundaries and solid walls where it is. This is shown in Fig. 8.1 The viscous region is called a
hydrodynamic boundary layer. A similar division can be made for thermal analysis, as shown in
Fig. 8.2. The region near a solid wall where temperature gradients are steep and conduction is
important is the thermal boundary layer. Outside this region conduction eects can be neglected.
We take the leading edge of the plate to be the origin of coordinates, with x being along and
y normal to the plate. The thicknesses of the layers (x) and T (x) have been exaggerated in the
gure. For large Reynolds number, they are usually very thin, so that we can make the following
approximations within the boundary layers: v  u; @=@x  @=@y. Under these conditions, the
governing equations can be simplied.
According to the y-momentum equation, the pressure does not depend on y. The pressure p(x)
is thus determined by the uid mechanics outside the boundary layer where viscous eects are not
important. Using the Bernoulli equation, we have
1 @p
dU
=U
 @x
dx

(8.1)

where U (x) is the free-stream velocity outside the boundary layer.

The boundary layer continuity, x-momentum, and energy equations are

@u @v
+
= 0
@x @y
@u
@u
dU
@ 2u
u +v
= U
+ 2
@x
@y
dx
@y
2
@T
@T
@ T
u +v
=  2
@x
@y
@y
Viscous dissipation has been neglected in the energy equation.
99

(8.2)
(8.3)
(8.4)

U-

inviscid

y
6
aa- x

viscous

6
?

Figure 8.1: Hydrodynamic boundary layer over a at plate

U-

non-conductive

conductive

6
aa- x

Figure 8.2: Thermal boundary layer over a at plate

The boundary conditions for these equations are

u = v = 0; Tw at y = 0
u = U; T = T1 as y ! 1

(8.5)
(8.6)

These equations can be written in terms of a stream function (x; y) where

@
(8.7)
@y
@
v =
(8.8)
@x
satisfying the continuity equation (8.2). The x-momentum and energy equations, (8.3) and (8.4)
respectively, become
u =

@ @2
dU
@3
@ @2
=
U
+

@y @x@y @x @y2
dx
@y3
2
@ @T @ @T
@T
=  2
@y @x @x @y
@y

(8.9)
(8.10)

with the boundary conditions

@
@
=
= 0; T = Tw at y = 0
@x @y
@
= U; T = T1 as y ! 0
@y
100

(8.11)
(8.12)

The thickness of the boundary layer can be dened in several ways. One of the most common is
where u() = 0:99U . Another, the displacement thickness, is dened by
Z 1
u
 =
1
dy
(8.13)
U
0
The momentum thickness is

=

1u
1
U
0

u
dy
U

(8.14)

8.2 Blasius solution

Consider a at plate where the velocity U is a constant. We will have a zero pressure gradient
dp=dx = 0, so that the pressure is constant everywhere. Now let
r

U
p x
= xUf ()

 = y

The momentum equation becomes the Blasius's equation

2f 000 + ff 00 = 0

(8.15)
(8.16)
(8.17)

with

f = f 0 = 0 at  = 0
f 0 = 1 as  ! 1

(8.18)
(8.19)

This problem can be integrated numerically; the results of f () and f 0 () are shown in Fig. 8.3.
Since u = Uf 0 , the f 0 curve also represents u=U .
According to the numerical solution u() = 0:99U at  = 5. Thus from equation (8.15) we have
r

5=

U
x

This can be written as

p5
=
x
Rex
where the local Reynolds number Rex = Ux= . The local shear stress at the wall

(8.20)
(8.21)

can be written as

@u 
wx =  
@y y=0
r

(8.22)

U 3 00
f (0)
(8.23)
x
Since f 00 (0) is numerically found to be 0.332, we can obtain the local coecient of friction dened
by Cfx = wx = 21 U 2 to be
0:664
Cfx = p
(8.24)
Rex
wx = 

101

f ( )
3

f 0 ( )

Figure 8.3: Solution of Blasius's equation

Using the similarity variable , and T = T (), the energy equation becomes
1
(8.25)
T 00 + P r fT 0 = 0
2
To obtain simple boundary conditions, the temperature can be nondimensionalized as
T Tw
T =
(8.26)
T1 Tw
The temperature prole is then given by
1
(8.27)
T 00 + P r fT 0 = 0
2
with the boundary conditions
T  = 0 at  = 0
(8.28)

T = 1 as  ! 1
(8.29)
Again this problem can be solved numerically for a given Prandtl number. For P r > 0:6, it is seen
that a good approximation for the slope of the temperature prole at the wall is
T 0 (0) = 0:332 P r1=3 for P r > 0:6
(8.30)
102

The local convective heat transfer coecient may be expressed as

k(@T=@y)y=0
Tw T1
The local Nusselt number dened as Nux = h(x)x=k comes out to be
h(x) =

(8.31)

Nux = 0:332 Re1x=2 P r1=3 for P r > 0:6

(8.32)

= P r1=3 for P r > 0:6

T

(8.33)

In addition it is found that

8.3 Momentum integral method

Equation (8.3) can be written as

@2u
@ (u2 ) @ (uv)
+
= 2
@x
@y
@y
which, on integrating across the boundary layer from y = 0 to y = , becomes

(8.34)

@ (u2 )
@u 
dy + Uv(x; ) =  
(8.35)
@x
@y y=0
0
since u(x; ) = U; (@u=@y)y= = 0. Integrating the continuity equation (8.2) within the same limits,
we get
Z
@u
dy
(8.36)
v(x; ) =
@x
0
Applying the Leibnitz formula
Z
Z  (x)
d (x)
@f (x; y)
d
f (x; y) dy =
dy + f (x;  )
dx (x)
@x
dx
(x)

we can write

d
dx

Z
0

u(U

u) dy = U

"Z

@u
d
dy + U
@x
dx

"Z

f (x; )

d
dx

@ (u2 )
d
dy + U 2
@x
dx

(8.37)
#

(8.38)

Substituting equations (8.35) and (8.36) into the right hand side, we get the momentum integral
equation

Z
d
@u 
u(U u) dy =  
(8.39)
dx 0
@y y=0
For an approximate solution, assume a cubic polynomial for the velocity prole
 y 2
u
y
= a0 + a1 + a 2
U

Using the conditions


@u 
u(x; 0) = 0; u(x; ) = U and  = 0
@y y=
103

(8.40)
(8.41)

the constants can be determined as

a0 = 0; a1 = 2 and a2 = 1
so that

u
= 2  2
U
where
y
=

Substituting into the momentum integral equation (8.39) gives



d 2 2
2U
U =
dx 15

(8.42)
(8.43)
(8.44)

(8.45)

Using the initial condition (0) = 0, the solution of this equation is

= 30

x
U

which in nondimensional form is

5:48
=p
x
Rex
The local coecient of friction is found to be
0:73
Cfx = p
Rex

(8.46)
(8.47)
(8.48)

These can be compared with the numerical results of equations (8.21) and (8.24).

8.4 Vertical plate natural convection

Let us consider a at plate at temperature Tw which is higher than the temperature T1 of the
surroundings. Fluid ow is then upwards in the vicinity of the plate. The Boussinesq approximation
can be used for the governing equations. We make the usual boundary layer assumptions and take
u  v; @=@y  @=@x, where y is normal to the plate and x is upwards.
The boundary layer equations are

@u @v
+ =0
@x @y

@u
@u
@2u
+ v =  2 + g (T T0 )
@x
@y
@y
@T
@2u
@T
u +v = 2
@x
@y
@y

(8.49)
(8.50)
(8.51)

where  is the coecient of thermal expansion (for example, for a perfect gas,  = 1=T , where T is
the absolute temperature.) The boundary conditions are

u = v = 0; T = Tw at y = 0
u = 0; T = T1 at y ! 1
104

(8.52)
(8.53)

Similarity solutions can be obtained by using the variables

p1 Grx= xy

 =

(8.54)

1 4

 df
x d 
 df
Grx1=4

x d
T1
T1

u = 2 Grx1=2

p1

v =

2
T

T =
Tw

(8.55)
3f

(8.56)
(8.57)

where the local Grashof number is

Grx =

g(Tw

T1 )x3

(8.58)

The continuity equation (8.49) is satised. The y-momentum equation (8.50) and the energy equation (8.51) become


d2 f
df 2
d3 f
+
3
f
2
+ T = 0
d3
d2
d
dT 
d2 T 
+ 3 Pr f
= 0
2
d
d

(8.59)
(8.60)

The boundary conditions are

f = f 0 = 0; T  = 1 at  = 0
f 0 = T  = 0 as  ! 1
This problem can be numerically solved.
The local Nusselt number is

Nux =
=

@T
Tw T1 @y y=0


1
@T 
1= 4
p Grx @
2
=0
x

(8.61)
(8.62)

(8.63)
(8.64)

Numerical computations show that the values approximate



1=4
Grx P r2
Nux =
2:435 + 4:884 P r1=2 + 4:953 P r

Problems

(8.65)

1. Show that the total drag per unit

length FD = 0L wx dx on a at plate of length L can be expressed as a
p
coecient of drag CD = 1:328= ReL .
R
2. Show that the average heat transfer coecient h = L1 0L h(x) dx on a at plate of length L can be expressed
as a Nusselt number NuL = 0:664Re1L=2 P r1=3 for P r > 0:6.
105

3. Using

4.
5.
6.
7.
8.

9.
10.

u(x; y)

=
v(x; y) =

U (x)f 0 ()
Z
@ y
u(x; y) dy
@x 0

where  = (x; y) and U (x) is the inviscid x-velocity outside the boundary layer, show that the boundary layer
equations reduce to a form
f 000 + g1 f 00 = g2 ff 00 + g3 (f 02 ff 00 1)
Find g1 , g2 and g3 .
Use the momentum integral method with the velocity prole u = U sin(y=2) to determine =x and Cfx.
A square
at plate of side 30 cm, maintained at temperature 50 C, is exposed to a ow of air (density = 1.2
kg/m3 , viscosity = 1:9  10 5 Ns/m2 , Prandtl number = 0.71) at 20 C parallel to it. Find the drag force
on one side of the plate as well as the rate of heat transfer from it. Determine also the thicknesses of the
hydrodynamic and thermal boundary layers at the trailing edge of the plate.
Find the extra term Rin the at plate boundary layer problem with viscous dissipation.
Show that exp( 12 P r f () d) is an integrating factor for equation (8.27). Solve for T  in terms of f ().
Given T  = (T T1 )=(Tw T1 ) and the similarity variables
1
y
 = p Grx1=4
2p x
2
f () =
4 Grx1=4
where (x; y) is the stream function, reduce the boundary layer equations (8.50){(8.53) for natural convection
in a vertical at plate to ordinary dierential equations (8.59){(8.60). Transform also the boundary conditions.
Determine the extra term in the equation (8.60) corresponding to viscous dissipation, if it is included in the
vertical at plate natural convection analysis.
Show that the momentum integral equation corresponding to the vertical at plate natural convection problem
is

Z
Z

d 2
@u 
g
(T T1 ) dy dx u dy =  @y 
(8.66)
0

y=0

11. An incompressible uid enters the space between two at plates at velocity U . Estimate the distance downstream at which the two boundary layers meets. State your assumptions.
12. The one-seventh power law for the velocity prole u = U (y=)1=7 is often used for a turbulent at plate
boundary layer. Find =x and Cfx using the momentum integral. Wrong?
13. Using a cubic velocity prole for a laminar boundary layer over a at plate, nd the boundary layer thickness.
14. In the at plate boundary layer problem, use the momentum integral method with a linear velocity prole
to determine the (a) boundary layer thickness as a function of distance downstream x, and (b) the local
coecient of friction.
15. Findthe rate of heat transfer by natural convection from one side of a vertical 10 cm 3 10 cm
at plate which
is 10 C hotter than the ambient air. The coecient of expansion of air  = 3:4  10 K 1 .

106

Chapter 9

Hydrodynamic stability and

turbulence
9.1 Hydrodynamic stability
We will consider the stability of the parallel ow of an incompressible uid with constant properties.
The eld velocity for this ow is u = u(y), v = w = 0, the pressure eld is p = p(x), and the ow
domain is y > 0 with a wall at y = 0. For this the governing equations are
0=

@2u
1 dp
+ 2
 dx
@y

(9.1)

where the bars represent the basic ow. Let us apply small two-dimensional time-dependent perturbations so that

u(x; y; t) = u(y) + u0 (x; y; t)

v(x; y; t) = v0 (x; y; t)
p(x; y; t) = p(x) + p0 (x; y; t)

(9.2)
(9.3)
(9.4)

where the primed quantities are small. Substitute in the governing equations, neglect terms with
products of primed terms, and use equation (9.1) to get
@u0 @v0
+
= 0
(9.5)
@x @y


1 @p0
@ 2u0 @ 2 u0
@u0
@u0 0 du
+u
+v
=
+
+
(9.6)
@t
@x
dy
 @x
@x2 @y2
 2 0

@v0
@v0
1 @p0
@ v @ 2 v0
+u
=
+
+
(9.7)
@t
@x
 @y
@x2 @y2
Boundary conditions on the velocity components are
u0 = v0 = 0 for
u0 = v0 ! 0 for
107

y=0
y!1

(9.8)
(9.9)

Using the perturbation stream function dened by

@
u0 =
@y
@
v0 =
@x

(9.10)
(9.11)
(9.12)

and eliminating p0 , we get



@
@
+u
@t
@x



@2
@2
+
@x2 @y2

d2 u @
@4
@4
@4
=

+
2
+
dy2 @x
@x4
@x2 @y2 @y4

(9.13)

We make a normal mode expansion of the form

(x; y; t) = (y)ei(x ct)

(9.14)

where  is real and c = cr + ici is complex quantity. Thus we can write

(x; y; t) = (y)eci t ei(x cr t)

(9.15)

Physically, is the amplitude of the perturbation wave, ci is its temporal growth rate,  is its wave
number, and cr is its frequency.
With these variables we get the Orr-Sommerfeld equation


d2
(u c)
dy2

d2 u
 d4

=
dy2
i dy4

d2
2 2 + 4
dy
2

(9.16)

Boundary conditions are

d
=0
dy
d
=
!0
dy
=

for

y=0

(9.17)

for

y!1

(9.18)

The problem can be nondimensionalized using

u
(9.19)
u =
U
c
c =
(9.20)
U

 =
(9.21)
UL
 = L
(9.22)
y

y =
(9.23)
L
where the velocity scale is U and the length scale is L. The nondimensional Orr-Sommerfeld equation
is
 2 

 4 

2 
d
d2 u 
1
d



2 

2d
4 
(u c )

=
2
+
(9.24)
dy2
dy 2
iRe dy4
dy2
108

where

Re =

and with

UL


d 
 =  = 0
dy
d

 =  ! 0
dy

(9.25)

for

y = 0

(9.26)

for

y ! 1

(9.27)

From equation (9.15) we see that the imaginary part of c determines the stability of the ow; if
the imaginary part is negative it is stable, but if it is positive it is unstable. For a given ow eld
u (y) and values of  and Re, we have an eigenvalue problem for the determination of c that will
satisfy the boundary conditions. This is usually done numerically.

9.2 Turbulence
At low Reynolds number a ow is usually laminar. As Re is increased, it becomes unstable and
disturbances whicha re naturally present begin to grow. For higher Re, the ow becomes turbulent
with \random" uctuations of the velocity and pressure elds.
For fully turbulent ow we can write the velocty and pressure elds as
(9.28)
ui = ui + u0i
0
p = p+p
(9.29)
where the bars indicate the time-averaged value and the primes the uctuationg components. Note
that the time average of the primed quantities is zero.
Substituting into the governing equations for an incompressible uid with constant properties

@ui
= 0
@xi

@ui
@u
@p
@ 2ui

+ uj i
=
+
@t
@xj
@xi
@xj @xj

(9.30)
(9.31)

and averaging over time, we get

@ui
= 0
(9.32)
@x
i




@u
@p
@
@u
 uj i
=
+
 i u0i u0j
(9.33)
@xj
@xi @xj @xj
which resemble the equations for steady laminar ow in the time-averaged quantities except for the
additional stress term u0iu0j due to turbulence which is called the Reynolds stress. Additional
equations must be introduced for these terms before the governing equations can be solved. One
common method is to model it as


@ui @uj
(9.34)
u0i u0j = 
+
@xj @xi
where  is called the eddy viscosity.
109

There are various turbulence models for determining . For example, one of them is the mixing
length theory of von Karman for ow next to a wall which proposes that
 
 du 

dy 

 = 2 y2 

(9.35)

where  = 0:41.

Problems

1. For turbulent ow of an incompressible uid with constant properties, show that the Navier-Stokes equation
becomes




@u
@p
@
@u
+
 i u0i u0j
 uj i =
@xj
@xi @xj
@xj
where the overbars indicate a time average.
2. Show the following for the problem of two-dimensional stability of
parallel ow between at plates.
inviscid

6y

-x

(a) For the velocity eld u = U (y) + u0 (x; y; t), v = v0 (x; y; t), w = 0, neglect the products of primed quantities
and show from the governing equations that

@
@  2 0 d2 U @v0
+
U
r v dy2 @x = 0
@t
@x
with v0 = 0 at y = 0; 1, where 1 is the distance between the plates.
(b) Substitute
v0 (x; y; t) = v^(y)ei(x ct)
where  is real, and v^ and c are complex. Show that


d2 U
d2 v^
2
 v^
v^ = 0
(U c)
dy2

dy2

with v^ = 0 at y = 0; 1.
(c) Divide by U c, multiply by the complex conjugate v^ , integrate by parts, and use the boundary conditions
to show that

Z 1  2
Z 1
d2 U jv^j2
 dv^ 
2
2
+

j
v
^
j
dy
+
dy = 0
 
dy
dy2 U c
0

(d) Show that the imaginary part of this is

ci

1
0

d2 U
dy2 (U

jv^j2

cr )2 + c2i

dy = 0

where c = cr + ici , which can only be satised if d2 U=dy2 changes sign in the interval y = 0; 1.

110

Chapter 10

Compressible ow in gases
10.1 Basic equations
10.1.1 Equations of state for real gases

Perfect or ideal gas

p = RT

(10.1)

Van der Waals

p=

RT

a2

(10.2)

a2
(1 + b)T 1=2

(10.3)

Redlich-Kwong
p=

RT

Virial equation
p = RT + a(T )2 + b(T )3 + : : :

Beattie-Bridgeman
p = 2 RT

(10.4)

1
+ b (1 ) A2


(10.5)

RT
1 b

(10.6)

Berthelot
p=

111

a2
T

10.1.2 Thermodynamics of perfect gases

We will assume for the most part that the ow is that of a (thermally) perfect gas which satises
the relation
p = RT
(10.7)
where R is the particular gas constant.
The specic internal energy e and specic enthalpy h = e + p= are given by

de = cv dT
dh = cp dT

(10.8)
(10.9)

where cv and cp are the specic heats at constant volume and constant pressure respectively. The
dierence is
cp cv = R
(10.10)
and the ratio is
cp
=
(10.11)
cv
From these
R
cp =
(10.12)
1
R
cv =
(10.13)
1
In this chapter the specic heats will be taken to be constant (i.e. calorically perfect gas). Thus the
integrated versions of equations (10.8) and (10.9) are

e2 e1 = cv (T2 T1)
h2 h1 = cp (T2 T1 )
The thermodynamic properties are related by the Gibbs equations
p
T ds = de 2 d

dp
= dh


(10.14)
(10.15)
(10.16)
(10.17)

Changes in specic entropy are determined from

dT
T
dT
= cp
T

ds = cv

d

dp
R
p
R

(10.18)
(10.19)

which can be integrated to give

T2
T1
T2
= cp ln
T1

s2 s1 = cv ln

112

2
1
p2
R ln
p1
R ln

(10.20)
(10.21)

((
(
(
(
(
(
(
(
(( (
((((((
? ?
 ; u ; p - Control
 ;u ;p
volume
Q_

W_

hhhhhhA hfriction A
h hhhhhh
hhhh
1

Figure 10.1: One-dimensional ow

For an isentropic process s2

s1 = 0, from which
T2
=
T1
T2
=
T1
p2
=
p1

2 1
1
 ( 1)=
p2
p1
 
2
1

(10.22)
(10.23)
(10.24)

10.2 One-dimensional steady ow equations

Consider the ow of a compressible uid through a variable area duct. Take a control volume
as indicated in Fig. 10.1 where the conditions upstream and downstream are shown. The onedimensional governing equations for steady ow conditions in a nite control volume can be written
as follows.
Mass:
Momentum:
Energy:

m_ = 1 u1 A1 = 2 u2 A2
X

F = m_ (u2

(10.25)

u1 )

(10.26)


u22
u2
=0
(10.27)
Q_ + W_ + m_ h1 + 1 h2
2
2
where the
P subscripts 1 and 2 denote quantities at the inlet and outlet of the control volume, respectively; F is the sum of all forces on the control volume including pressure and frictional forces; m_
is the mass ow rate through the duct; Q_ and W_ are the heat ow rate and work ow rate to the
control volume.
113

10.2.1 Dierential form

Dierential relations can also be written for an innitesimal control volume. In this case 1 = ,
u1 = u, p1 = p and A1 = A at section 1, and 2 =  + d, u2 = u + du, p2 = p + dp and A2 = A + dA
at section 2.
The mass conservation equation is

uA d + A du + u dA = 0

(10.28)

where products of dierentials have been neglected. Dividing by uA, we have

d du dA
+ +
=0
(10.29)

u
A
This could also have been obtained by taking the logarithm of uA = constant, and dierentiating.
In the momentum equation the forces are
X
dp
F = pA + (p + ) dA (p + dp)(A + dA) w P dx
(10.30)
2
P is the perimeter of the section and w is the shear stress at the wall. The pressure on the lateral
surface is assumed to be p + dp=2. The wall shear stress is often approximated by
1
w = fu2
(10.31)
8
where f is the Darcy friction factor that is a function of the local Reynolds number and wall
roughness. Thus, the momentum equation is

u2
dx = 0
(10.32)
2Dh
where products of dierentials have been neglected again. Dh is the hydraulic diameter of the duct
dened by Dh = 4A=P .
In the energy equation, we neglect work transfer W_ and write the heat transfer rate as Q_ = q0 dx,
where q0 is the heat in ow per unit length. Thus
q0
dx + dh + u du = 0
(10.33)
m_
where m_ = uA.
dp + u du + f

10.2.2 Constant area without friction, heat or work transfer

If the area is constant, the mass conservation equation (10.25) is

1 u1 = 2 u2

(10.34)

The constant area momentum equation (10.26) without friction at the wall becomes

p1 + 1 u1 = p2 + 2 u2

(10.35)

Without heat and work transfer the energy equation (10.27) reduces to

h1 +

u2
u21
= h2 + 2
2
2
114

(10.36)

wave moving with

velocity
-c

uid velocity u
 + ; p + p; T + T

uid at rest
; p; T

(a) Moving wave

wave stationary

  uid velocity c

uid velocity c u
 + ; p + p; T + T

; p; T

(b) Fixed wave

Figure 10.2: Pressure wave

10.3 Speed of sound

Consider the propagation of a small pressure, temperature, density, etc. disturbance with a velocity
c. The equations of conservation can be easily written if we consider the steady problem in a frame
moving with the wave as shown in Fig. 10.2.
The mass conservation equation (10.25) is

cA = ( + )(c u)A

(10.37)

where A is the cross-sectional area. Thus

u = c


 + 

(10.38)

The momentum equation (10.26) is

pA (p + p)A = Ac [(c u) c]

(10.39)

p = c u

(10.40)

from which
Combining equations (10.38) and (10.40)

c2 =

p

1+



As  ! 0, we have

c=
115

dp
d

(10.41)

(10.42)

motion
 of
source

#

SoS
qk

"!

A B

(a) subsonic case

#

q 
k"!

XXz

C
SoS

Z
A Z B
ZZ
Z

Mach cone
motion
 of
source

(b) supersonic case

Figure 10.3: Propagation of sound
This is the speed of an innitesimal pressure wave, commonly termed an acoustic wave or sound.
For an isentropic process
p = constant 
(10.43)
from which dp=d = constant 

= p=. Thus
r

p

p
=
RT

c =

(10.44)
(10.45)

For air at 20C, c = 343:2 m/s. Equation (10.42) can be directly used for liquids or solids; typical
values are c = 1400 m/s for water, and c = 6000 m/s for steel.
A sound pattern from a moving source depends on whether it is moving faster (supersonic) or
slower than the speed of sound (subsonic). Figures 10.3 (a) and (b) represent the position at an
instant of time of the spherical wave fronts emitted when the source was at dierent points along
its line of travel. A is the current position of the source, B its position a time
t ago, and BC the
radius of the wave front which grew from B . Thus, if U is the speed of the source, AB = U
t and
BC = c
t. In the supersonic case, the sound is kept conned to a Mach cone, the half angle of
which is  = 6 CAB , where

c
U
1
=
M

sin  =

where the Mach number M = U=c.

116

(10.46)
(10.47)

10.4 Stagnation and critical properties

The thermodynamic properties obtained on reducing the velocity of a ow through a frictionless,
adiabatic process are dened as the stagnation properties (indicated by the subscript 0). Thus, in
the case of enthalpy
u2
h0 = h +
(10.48)
2
The stagnation temperature is given by

T0 = T +

u2
2cp

(10.49)

To write this in terms of the Mach number

1 2
T0
=1+
M
T
2
since u2 =cp T = u2 ( 1)= RT = ( 1)u2=c2 = (
In addition, because the process is isentropic

1)M 2.


T0 =( 1)
T


1 2 =(
= 1+
M
2
 1=( 1)
T0
=
T


1 2 1=(
= 1+
M
2

p0
=
p
0


(10.50)

1)

1)

(10.51)

(10.52)

Critical properties (indicated by *) are those obtained when the ow is accelerated or decelerated
to sonic velocity in an isentropic manner. At M = 1, equations (10.50), (10.51) and (10.52) become

T0
+1
=
T
2


+ 1 =(
p0
=
p
2


0
+ 1 1=(
=

2

(10.53)
1)

1)

(10.54)
(10.55)

10.5 Normal shocks

Consider Fig. 10.2 (b) in which the wave is of large amplitude, i.e. p, , T and u are not
necessarily small. In the absence of heat and work transfer, equations (10.27) and (10.48) indicate
that the stagnation enthalpy and hence stagnation temperature are the same on either side of the
shock. Thus
T01 = T02
(10.56)
117

Using equation (10.50)

T2 2 + ( 1)M12
=
T1 2 + ( 1)M22
p
Since T2 =T1 = (1 =2 )(p2 =p1 ) = (u2 =u1)(p2 =p1) = (M2 =M1 )(p2 =p1) (T2 =T1 ), we have
p2 M1
=
p1 M2
Also, since 2 =1 = (p2 =p1)(T1 =T2 ),

2 M1
=
1 M2
For the stagnation pressures
so that

(10.57)

2 + (
2 + (

1)M12
1)M22

(10.58)

2 + (
2 + (

1)M22
1)M12

(10.59)

p =p p
p02
= 02 2 2
p01
p01 =p1 p1


p02 M1 2 + ( 1)M12 ( +1)=[2( 1)]

=
p01 M2 2 + ( 1)M22
Since u2 = (p=RT )u2 = pu2=c2 , the momentum equation (10.35) reduces to
p1 + p1 M12 = p2 + p2M22
from which
p2 1 + M12
=
p1 1 + M22
Combining with equation (10.58)

(10.60)
(10.61)
(10.62)
(10.63)

M1 2 + ( 1)M12 M2 2 + ( 1)M22
=
1 + M12
1 + M22
This can be solved to give either M1 = M2 (no shock) or
s

M2 =

2 + ( 1)M12
2 M12 ( 1)

(10.64)

(10.65)

Equation (10.65) can be substituted into equations (10.57){(10.61) to get T2 =T1, p2 =p1, 2 =1
and p02 =p01 in terms of either M1 or M2 . Thus, we have


p2
2
= 1+
M2 1
(10.66)
p1
+1 1
2
( + 1)M12
=
(10.67)
1
( 1)M12 + 2
The shock is a process with entropy change. Equations (10.20) and (10.21) can be reduced to




s2 s1
1

2 + ( 1)M12
2 M12 ( 1)
+
=
ln
ln
(10.68)
R
1
( + 1)M12
1
+1
Figure 10.4 shows that s2 s1 < 0 if M1 < 1, which would violate the second law of thermodynamics. So a subsonic to supersonic shock is not physically possible.
118

2
s2 s1
R

1.5

0.5

M1

Figure 10.4: Entropy change as function of incoming Mach number

10.6 Oblique shocks

Figure 10.5 shows an oblique shock in which the ow is de ected through an angle . The tangential
velocity is the same on either side of the shock, so that

u1 cos  = u2 cos(

(10.69)

The normal velocity undergoes a shock as described in the previous section. We can use the normal
shock relations as long as we substitute M1 sin  for the Mach number on the upstream side and
M2 sin( ) on the downstream side. Thus

From continuity

( + 1)M12 sin2 
2
=
1 ( 1)M12 sin2  + 2

(10.70)

1 u1 sin  = 2 u2 sin(

(10.71)

and equation (10.69) we get

2
tan 
=
1 tan( )
Equating the two expressions for 2 =1 and simplifying, we get
M12 =

2 cos( )
sin  [sin(2 ) sin ]
119

(10.72)

(10.73)



u1

u

 >
 

-
 





Figure 10.5: Oblique stationary shock

For given M1 and , two values of  are possible. One is called a strong shock and the other a
weak shock.
The downstream Mach number is
s

1
2 + ( 1)M12 sin2 
M2 =
sin( ) 2 M12 sin2  ( 1)

(10.74)

The pressure, density, temperature, stagnation temperature and stagnation pressure ratios are
given by

p2
p1
2
1
T2
T1
T02
T01
p02
p01

2
M 2 sin2  1
+1 1
( + 1)M12 sin2 
=
( 1)M12 sin2  + 2

 2 M12 sin2  ( 1)
= 2 + ( 1)M12 sin2 
( + 1)2 M12 sin2 

= 1+

(10.75)
(10.76)
(10.77)

= 1
=

(10.78)

( + 1)M12 sin2  =(
2 + ( 1)M12 sin2 

1)

1=(
+1
2 M12 sin2  ( 1)

1)

(10.79)

10.7 Flow in ducts

The momentum equation (10.32), divided by p becomes

dp + u du + f
are

u2
dx = 0
2Dh

(10.80)

To summarize, the governing equations for a variable area duct with heat transfer and friction

p = RT
120

(10.81)

d du dA
+ +
= 0

u
A

(10.82)

dp + u du =

(10.83)

u2
2

(10.84)

d h+
In addition we will use

u2
f
dx
2Dh
q0 dx
=
m_

RT
c =
u
M =
c
R
cv =
1
R
cp =
1

(10.85)
(10.86)
(10.87)
(10.88)

10.7.1 Adiabatic, frictionless ow in variable area ducts

Thermodynamically, the process is isentropic. Since
du dA
d
+
=
u
A

du
= u d from equation (10.83)
dp
d du
= u2
dp u
2
u du
= 2
c u
du
= M2
u
Thus
1 dA
du
=
u M2 1 A
Also du=u = dp=u2, so that
dp
1
dA
=
2
2
u
1 M A
2
2
2
Furthermore, since dp=u = (c =u )(d=) = (1=M 2)(d=), we have

Similarly
From the equation of state

(10.89)
(10.90)

d
M 2 dA
=
 1 M2 A

(10.91)

dM 2 + ( 1)M 2
=
dA
M
2(M 2 1)

(10.92)

dT dp
=
T
p
121

d


(10.93)

M <1
M >1
(subsonic) (supersonic)
dA < 0
du > 0
du < 0
(converging) dM > 0
dM < 0
dp < 0
dp > 0
d < 0
d > 0
dT < 0
dT > 0
dp0 = 0
dp0 = 0
d0 = 0
d0 = 0
dT0 = 0
dT0 = 0
ds = 0
ds = 0
dA > 0
du < 0
du > 0
(diverging)
dM < 0
dM > 0
dp > 0
dp < 0
d > 0
d < 0
dT > 0
dT < 0
dp0 = 0
dp0 = 0
d0 = 0
d0 = 0
dT0 = 0
dT0 = 0
ds = 0
ds = 0
Table 10.1: Duct ow with variable area
Since u2=p =  RT M 2=p = M 2, we have

dT
= (
T

1)

M 2 dA
1 M2 A

(10.94)

Table 10.1 shows the changes in velocity, pressure, density and Mach number which take place in
converging (dA < 0) and diverging (dA > 0) ducts. The behavior of subsonic and supersonic ows
is seen to be dierent. A duct that decreases the uid velocity is called a diuser and that which
increases it is a nozzle.
Equating the mass ow rate to that at critical conditions, uA =  u A , so that
A
 u  0 u
=
=
A
 u 0  u
p
 0 RT 
=
0  u
r
r
p
 0 RT  T  T0
=
0  u
T0 T



1
2
1 2 ( +1)=[2( 1)]
=
1+
M
(10.95)
M +1
2
where we have used equations (10.50), (10.52), (10.55), and (10.53). As shown in Fig. 10.6 this has
a minimum at M = 1 where Amin = A .
122

A=A
5
4
3
2
1
0

0.5

1.0

1.5

2.0

2.5

3.0

Figure 10.6: Variation of area with Mach number

From

p up

m_ = uA =
RT A = pM
A
RT c
RT


r

2 ( +1)=[2( 1)] 
A
= p0
RT0 + 1
we get

m_ max =


2 ( +1)=(
+1

1)

A 0 RT0

(10.96)

(10.97)

10.7.2 Frictionless ow with heat transfer in constant area ducts

1 and 2 are two sections along the length of a duct that is frictionless, but with heating. The
momentum equation gives
p1 + 1 u21 = p2 + 2 u22
(10.98)
Since u2 = (p=RT )u2 = pM 2, we have

p2 1 + M12
=
p1 1 + M22
123

(10.99)

Using equation (10.51) for the stagnation pressure, we have



 

p02
1 + M12
=
p01
1 + M22

1 + 2 1 M22
1 + 2 1 M12

! =(

1)

(10.100)

Similarly, it can be shown that

T2
T1
T02
T01
2
1
u2
u1
s2 s1
R

=
=
=
=
=

1 + M12 2 M2 2
1 + M22
M1

 
 

2 2
1 + M1
M2 2 2 + ( 1)M22
1 + M22
M1
2 + ( 1)M12




2
1 + M22
M1
1 + M12
M2

 

1 + M12 2 M2 3
1 + M22
M1
"
#



2 ( +1)=( 1)
1 + M1
M2 2 =( 1)
ln
1 + M22
M1

(10.101)
(10.102)
(10.103)
(10.104)
(10.105)

Table 10.2 shows the changes in ow variables for subsonic or supersonic ow with either heating
or cooling.
The energy equation is
u2
u2
h1 + 1 + Q_ = h2 + 2
(10.106)
2
2
from which
h02 h01 = Q_
(10.107)
and
cp (T02 T01 ) = Q_
(10.108)
Thus


Q_
T
T
= 01 02 1
(10.109)
cp T1
T1 T01
#
 " 2  
 


M2
1 + M12 2 2 + ( 1)M22
1 2
M1
1
(10.110)
= 1+
2
M12
1 + M22
2 + ( 1)M12
The maximum heat is transferred when M2 = 1. Thus
2
(Min
1)2
Q_ max
=
2
cp Tin 2Min( + 1)

(10.111)

where Tin and Min are the inlet temperature and Mach number. respectively. A T (s) graph is called
the Rayleigh line.
Example 10.1

Air ows along a frictionless duct of diameter 5 cm with heat transfer. The mass ow rate and the inlet
temperature are 1 kg/s and 300 K, respectively. The duct is long enough for the ow to reach sonic velocity at
124

M <1
M >1
(subsonic) (supersonic)
Q_ < 0
du < 0
du > 0
(cooling) dM < 0
dM > 0
dp > 0
dp < 0
d > 0
d < 0

dT < 0
dp0 > 0
dp0 > 0
d0 > 0
d0 > 0
dT0 < 0
dT0 < 0
ds < 0
ds < 0
Q_ > 0
du > 0
du < 0
(heating) dM > 0
dM < 0
dp < 0
dp > 0
d < 0
d > 0
**
dT > 0
dp0 < 0
dp0 < 0
d0 < 0
d0 < 0
dT0 > 0
dT0 > 0
ds > 0
ds > 0
Table 10.2: Duct ow with heat transfer; * dT < 0 for M < 1=2 , dT > 0 for M > 1=2 ; **
dT > 0 for M < 1=2 , dT < 0 for M > 1=2

its end. The heat rate in is 1000 W/m. Find the Mach number, velocity, temperature, pressure, density and
entropy at dierent distances along the duct, if (a) Min = 0:1, and (b) Min = 2. Take sin as the reference
entropy. Also draw the Rayleigh line.
We march from M = Min to M = 1 in steps of 0.1, using the following equations:


2
1
+
Min
M 2
T = Tin
1 + M 2 Min

1 2
T0;in = Tin 1 +
2 Min



1 + Min2 2  M 2 2 + ( 1)M 2
T0 = T0;in
1 + M 2
Min
2 + ( 1)Min2
p
u = M RT
4m_
 =
Dh2 u
p = RT
R
x =
( 1)q0 (T0 T0;in )
"

1 + Min2 ( +1)=( 1)  M 2 =( 1)
s sin = R ln
1 + M 2
Min
where Dh = 0:05 m, m_ = 1 kg/s, Tin = 300 K, q0 = 1000 W. For air, R = 287 J/kg s, and = 1:4.
(a) Subsonic case: Min = 0:2. See Table 10.3.
125

(m)
0
303.3285
622.1643
906.2821
1129.551
1286.332
1383.38
1432.606
1446.48

.2
.3
.4
.5
.6
.7
.8
.9
1

(m/s)
69.43774
146.5223
239.6283
339.4734
438.7874
532.7679
618.787
695.8091
763.8152

(K)
300
593.6832
893.1949
1147.259
1331.055
1441.684
1488.992
1487.597
1452

(Pa)
631506.8
592247.9
544829.3
493978.7
443398.4
395534.5
351725.3
312498.2
277863

(kg/m3 )
7.334573
3.475897
2.125359
1.500254
1.16069
.9559439
.8230557
.7319481
.6667794

s sin

(J/kg K)
0
704.0557
1138.304
1417.877
1598.144
1711.128
1777.251
1810.247
1819.632

Table 10.3: Subsonic case

x

(m)
0
14.18464
29.32387
45.34278
62.07742
79.23067
96.30785
112.5285
126.7066
137.1034
141.258

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

(m/s)
694.3774
683.1945
670.5458
656.1908
639.8452
621.1749
599.7907
575.244
547.0294
514.597
477.3844

(K)
300
321.7892
345.3839
370.8109
398.0154
426.8109
456.8078
487.3127
517.1887
544.6774
567.1876

(Pa)
63150.68
68846.13
75288.03
82598.99
90923.75
100432.4
111323.3
123824.9
138194.5
154712.2
173664.4

(kg/m3 )
.7334573
.745463
.7595249
.7761405
.7959679
.8198918
.8491232
.8853568
.9310216
.9896992
1.066847

s sin

(J/kg K)
0
45.64732
91.05396
135.8114
179.3697
220.9881
259.6678
294.0566
322.3155
341.9308
349.4442

Table 10.4: Supersonic case

(b) Supersonic case: Min = 2. See Table 10.4.

Figure 10.7 is the Rayleigh line. Notice that for the supersonic part the curve has been shifted to the right
so that the entropy for M < 1 and M > 1 match at the maximum point.

10.7.3 Adiabatic ow with friction in constant area ducts

Let 1 and 2 be two sections along the length of a duct that with frictionless, but adiabatic. The
energy equation is
u2
u2
h1 + 1 = h2 + 2
(10.112)
2
2
so that
T02
=1
(10.113)
T01
Thus
T2
T =T
= 01 1
T1
T02 =T2
126

1500
T

1250
1000

M <1

M>1

750
500
250
0

250

500

750

1000

1250

1500
s

1750 2000
sin (J/kg K)

Figure 10.7: Rayleigh line

=
on using equation (10.50).
Similarly, we can derive

p2
p1
p02
p01
2
1
s2 s1
R

=
=
=
=

2 + (
2 + (



1)M12
1)M22


M1
2 + ( 1)M12 1=2
M2
2 + ( 1)M22



M1
2 + ( 1)M22 ( +1)=[2( 1))]
M2
2 + ( 1)M12



M1
2 + ( 1)M12 1=2
M2
2 + ( 1)M22
"
#


M2
2 + ( 1)M12 ( +1)=[2( 1)]
ln
M1
2 + ( 1)M22

(10.114)

(10.115)
(10.116)
(10.117)
(10.118)

Table 10.5 shows the changes in ow variables in the ow direction in subsonic or supersonic
ow. The change in Mach number with distance x from the inlet is given by
f dx
4
1 M2
dM
=
(10.119)
2
2
Dh
M 2 + ( 1)M M
which can be integrated to give


M
f x
1
+1
2M 2
=
ln
(10.120)
Dh
M 2
2
2 + ( 1)M 2 Min
127

M <1
M >1
(subsonic) (supersonic)
du > 0
du < 0
dM > 0
dM < 0
dp < 0
dp > 0
d < 0
d > 0
dT < 0
dT > 0
dp0 < 0
dp0 < 0
d0 < 0
d0 < 0
dT0 = 0
dT0 = 0
ds > 0
ds > 0
Table 10.5: Duct ow with friction
where f is the average value of f , and M = Min at x = 0. If Lmax is the length to reach sonic
velocity, then
2
2
 max 1 Min
fL
+1
( + 1)Min
=
+
ln
(10.121)
2
2
Dh
Min
2
2 + ( 1)Min
The T (s) graph is called the Fanno line.
Example 10.2

Air ows along an adiabatic duct of diameter 5 cm. The mass ow rate and the inlet temperature are 1 kg/s
and 300 K, respectively. The duct is long enough for the ow to reach sonic velocity at its end. The average
friction factor is 0.02. Find the Mach number, velocity, temperature, pressure, density and entropy at dierent
distances along the duct, if (a) Min = 0:1, and (b) Min = 2. Take sin as the reference entropy. Also draw the
Fanno line.
We march from M = Min to M = 1 in steps of 0.1, using the following equations:


2 + ( 1)Min2
T = Tin
2 + ( 1)M 2
p
u = M RT
4m_
 =
Dh2 u
p = RT

M

Dh
2M 2
1 + 1 ln
x =
M 2
2
2 + ( 1)M 2
f
"

  2 + (

( +1)=[2(

Min

1)]
1)Min
s sin = R ln
2
2 + ( 1)M
where Dh = 0:05 m, m_ = 1 kg/s, Tin = 300 K, f = 0:02. For air, R = 287 J/kg s, and = 1:4.
(a) Subsonic case: Min = 0:2. See Table 10.6.
(b) Supersonic case: Min = 2. See Table 10.7.
Figure 10.8 is the Fanno line. Notice that for the supersonic part the curve has been shifted to the right so
that the entropy for M < 1 and M > 1 match at the maximum point.

M
Min

128

(m)
0
23.08504
30.56193
33.66051
35.10611
35.81282
36.15244
36.29688
36.33316

.2
.3
.4
.5
.6
.7
.8
.9
1

(m/s)
69.43774
103.6438
137.2511
170.087
201.9992
232.8589
262.5617
291.0284
318.2037

(K)
300
297.053
293.0233
288
282.0895
275.4098
268.0851
260.2409
252

(Pa)
631506.8
418931.6
312060.3
247499.1
204121.9
172877.7
149242.9
130705.1
115757.1

(kg/m3 )
7.334573
4.91391
3.710688
2.994327
2.521278
2.187145
1.93972
1.749988
1.600535

s sin

(J/kg K)
0
107.869
178.6735
227.8277
262.3004
285.9084
301.0229
309.258
311.7904

Table 10.6: Subsonic case

x

(m)
0
.0766702
.1577754
.2429831
.3315991
.422366
.5131462
.6004115
.6783964
.737654
.7624915

2
1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1
1

(m/s)
694.3774
674.4331
653.1241
630.3724
606.1022
580.2419
552.7264
523.5005
492.5214
459.7626
425.2175

(K)
300
313.5889
327.6699
342.2053
357.1429
372.4138
387.9311
403.5875
419.2547
434.7827
450

(Pa)
63150.68
67963.24
73331.94
79349.09
86128.81
93814.3
102588
112686.8
124424.3
138226.4
154687

(kg/m3 )
.7334573
.7551471
.7797847
.8079291
.8402811
.8777308
.9214255
.9728668
1.034059
1.107737
1.197731

s sin

(J/kg K)
0
23.42138
45.72262
66.6891
86.07614
103.6034
118.9502
131.7471
141.5664
147.9067
150.1722

Table 10.7: Supersonic case

10.8 Multi-dimensional ow
10.8.1 Stagnation enthalpy
Taking the dot product of the momentum equation


by the vector u we get

Du
=
Dt

rp

D 1
u
u =
Dt 2
The energy equation without heat conduction is


Dh Dp
=
Dt Dt
129

(10.122)

u
rp

(10.123)

(10.124)

(K)

300

M <1

250
200

M >1

150
100
50
0

50

100

150

200

250

300
350
(J/kg K)

s sin

Figure 10.8: Fanno line

Adding the two, we get

D
1
Dp
h+ u
u =
Dt
2
Dt

which is equivalent to


For steady ow, we have

u
rp

Dh0 @p
=
Dt
@t

(10.125)
(10.126)

Dh0
=0
Dt
which means that h0 is constant along a streamline.

(10.127)

10.8.2 Crocco's theorem

Equation (10.122) is
since

@u
1
+r u
u
@t
2

u! =


u
ru = r 12 u
u

Using the Gibbs relation

T ds = dh
130

1
dp


1
rp


u!

(10.128)
(10.129)
(10.130)

which is equivalent to

we get Crocco's equation

rs = rh 1 rp

u  ! + T rs = r

For steady ow this becomes

(10.131)


@u
1
h+ u
u +
2
@t

(10.132)

u  ! + T rs = 0

(10.133)

Thus irrotationality implies that the entropy is constant everywhere, and vice-versa.

10.8.3 Irrotational ow
If ! = r  u = 0, we have

u = r

(10.134)

It can be shown from the governing equations that

r
u = c12
from which

u
[(u
r)u] + @t@ @

+u
u
@t



(10.135)




r  = c1 r
[(r
r)r] + @t@ @
+ r
r
@t
2

(10.136)

The pressure, obtained from the energy equation, is



1
p = p1 1 + 2 U 2
2c1



=(

u
u

1)

(10.137)

where at innity the pressure and velocity are p1 and U . The pressure coecient dened by

Cp =
becomes

2
Cp =
2
M1

p p1
1
 U2
2 1

(

1
1 + 2 U2
2c1

u
u

(10.138)


=(

1)

(10.139)

10.8.4 Small-perturbation theory

For steady ow over a slender body, we can write

u = U ex + u 0

(10.140)

 = Ux + 
where u0 = r. Equation (10.136) becomes

(10.141)

or

r  = c1 (U ex + r)
[(U ex + r)
r] (U ex + r)

2

131

(10.142)

y = (x)

-x

Figure 10.9: Subsonic ow past wall

We assume that the devation from uniform ow is small, so that ju0 j  U . Under this condition,
we get
2 2
r2  = Uc2 @@x2
(10.143)
This can be further approximated by
2 2
(10.144)
r2  = cU2 @@x2
1
which gives

@ 2
@2 @2
2
1 M1
+ 2 + 2 =0
(10.145)
2
@x
@y
@z
In the small-perturbation approximation, the pressure coecient becomes
u0
Cp = 2
U
2 @
=
(10.146)
U @x
where u0 is the x-component of u0 .

10.8.5 Subsonic ow
Consider two-dimensional ow next to a wall of shape y = (x) as in Fig. 10.9. The smallperturbation approximation gives

@2
@2
2
+ 2 =0
(10.147)
1 M1
2
@x
@y
with boundary conditions
d
@
= U
at wall
@y
dx
@
= nite as y ! 1
@x
132

hy hhhhh hhhhy =hhu (x)

hhh
6- x y =hhh(x)
l

Figure 10.10: Airfoil in supersonic ow

With new variables dened by
0 =
y0 =
we have

p
p

2 
1 M1
2 y
1 M1

@ 2 0 @ 2 0
+
=0
@x2 @y02

(10.148)

with boundary conditions

@ 0
d
= U
at wall
@y0
dx
@ 0
= nite as y0 ! 1
@x
Thus

@ =@x
Cp
=
0
Cp
@ 0 =@x
1
= p
2
1 M1

(10.149)
(10.150)

which is the Prandtl-Glauert rule; Cp is the coecient of pressure in the compressible ow, and Cp0
is the corresponding value for incompressible ow.

10.8.6 Supersonic ow
Consider ow around an airfoil with upper and lower surfaces y = u (x) and y = l (x), respectively,
as shown in Fig. 10.10. Thus
@2
1 @2
=0
(10.151)
2
2
@x
M1 1 @y2
with
@
d
= U u;l
at wall
@y
dx
@
= nite as y ! 1
@x
This is a wave equation with a solution of the type
p

2
 = f (x  M1

133

1 y)

(10.152)

The solutions for the upper and lower parts of the ow eld are dierent, so that
=

2
fu (x p M1
1 y) upper
2
fl (x + M1 1 y) lower

(10.153)

The boundary condition at the wall gives

d
2
M1
1 fu0 = U u
dx
p
dl
0
2
M1 1 fl = U
dx

Thus, from equation (10.146), we get

Cp =

8
<
:

pM 2
1
pM 2
2

lower surface

2 0
f
U u;l

Cp =
which gives

upper surface

du
dx
dl
1 dx

(10.154)

upper surface
lower surface

(10.155)

Lift coecient
The coecient of lift is

L
1
(p p ) dx
(10.156)
CL = 1
 U 2L 0 l u
2 1
where L is the chord and pu and pl are the pressures at the upper and lower surfaces. Substituting

1
p = p1 + 1 U 2 Cp
2

(10.157)

we get

2
[l + u ]xx==0L
2
M1 1
Taking l (L) = u (L) = 0 and l (0) = u (0) = L, where  is the angle of attack

CL =

CL =

4
2
M1
1

(10.158)

(10.159)

Drag coecient
The coecient of drag is

CL =
which becomes

2
CD = p 2
M1

1
1

U 2L
2 1
1
1 L

Z L
0

(pl


Z L "
d 2

dx

134

pu ) dy


du
+
dx

(10.160)
2 #

dx

(10.161)

Problems

1. An airplane is ying at a height of 5 km at a Mach number of 1.5. Find L.

2. Show that for M  1, the equation
p0  1 2  =( 1)
= 1+ 2 M
p
can be expanded to
1
p0 = p + u2 + : : :
2
3. For a normal shock wave, show that (s2 s1 )=R = ln(p02 =p01 ).
4. An air stream with temperature 300 K undergoes a normal shock. The downstream Mach number is half that
upstream. Determine the air velocities upstream and downstream of the shock.
5. Show from equation (10.95) that A=A is a minimum at M = 1.
6. A converging-diverging air nozzle has a throat area of 1 cm2 . Find the inlet and exit areas if the Mach numbers
there are 0.1 and 2, respectively.
7. A large pressurized tank containing air at2 1 MPa pressure, 100 K temperature oating in space develops a
small leak through a hole of area 0.01 mm . Find the mass ow rate through the leak.
8. A stationary normal shock exists in a ow of a gas for which the ratio of specic heats is 1.3. If the upstream
Mach number is very large, nd (a) the downstream Mach number, and (b) the density ratio.
9. Conditions upstream of a shock in a ow of air are M1 = 3, T1 = 300 K, and p1 = 100 kPa. The downstream
pressure is 1 MPa. Find if the shock is normal or oblique and the downstream temperature and Mach number.
10. Air ows at 60 m/s into an adiabatic commercial steel pipe of diameter 5 cm at a temperature and pressure of
100 C and 1 MPa. Determine (a) the Mach number at a distance of 50 m down the pipe, and (b) the length
of the pipe for choked ow. Take the average friction factor to be 0.02.
11. Air at a stagnation pressure and stagnation temperature
of 300 kPa, 400 K ows isentropically through a
converging nozzle from an inlet area of 100 cm2 to an exit area of 50 cm2 . Find the mass ow rate for a back
pressure of 100 kPa.
12. Air enters a heated, frictionless duct of square cross-section at a pressure of 100 kPa, temperature 300 K, and
Mach number 0.2. Find the amount of heat added per unit mass necessary to choke the ow.
13. A supersonic wind tunnel is formed at the divergent section of a CD nozzle. The ow
in a 15 cm  15 cm
square test section of such a tunnel is at a Mach number of 3 at temperature 20 C and pressure 50 kPa.
Calculate the mass ow rate.
14. The inlet to a frictional, adiabatic, 10 cm diameter, 8.2 m long pipe is air at 100 kPa, 300 K, and Mach number
of 2.5. If the friction factor is 0.002, nd the temperature at the outlet.
15. Air at a stagnation pressure of 150 kPa ows isentropically through a converging nozzle. Find the back pressure
for which the nozzle is just choked.
16. Air ows through a diverging duct. Inlet conditions are: T =2 500 K, p = 1000 kPa, A = 0:01 m2, u = 1000 m/s;
outlet conditions are T = 1000 K, p = 100 kPa, A = 0:08 m . There is heat transfer between the surroundings
and the air in the pipe. Determine (a) the velocity at the outlet, (b) the mass ow rate, (c) direction and
magnitude of the heat transfer, and (d) change in specic entropy between inlet and outlet.
17. Find the speed of travel of a shock wave in air given that the pressure and temperature ahead are 0.1 MPa,
300 K, respectively, and the pressure behind is 0.5 MPa.
18. Air ows without friction through a duct. At the inlet the temperature
is 300 K, the pressure is 105 Pa, and
the velocity is 120 m/s. At the outlet the pressure is 0:5  105 Pa. Determine the change in Mach number and
entropy from the inlet to the outlet.
19. Using
uA = m_
u2
h+
2 = h0
p = RT
h = cp T
h   i
p 1
s s1 = cv log
p1 
135

show that the equation of the Fanno line is

s s1
cv

= log h(h0

h)
(

"

=  + log 1 R
p1 cp

1) 2

where the subscript 1 refers to inlet conditions.

20. Using

2A2
m_ 2

(

1) 2

= m_
h+
2 = h0
p = RT
h = cp T
h   i
p 1
s s1 = cv log
p1 
show that the equation of the Fanno line is
"
#
R  2A2 ( 1)=2



s s1
= log h(h0 h)( 1)=2 + log p11 cp m_ 2
cv
where the subscript 1 refers to inlet conditions.
21. Air ow steadily through a round tube of diameter D with both friction and heat transfer. Conditions T1 , p1 ,
and M1 at the inlet and p2 and M2 at the outlet are known. Calculate the heat added per unit mass and the
friction force exerted by the gas on the tube.
22. Find the drag coecient for the double-wedge airfoil in supersonic ow at zero angle of attack using the slender
body approximation.
uA
u2

136

Chapter 11

Compressible eects in liquids and

two-phase ow
11.1 Waterhammer
Compressibility in liquids leads to the eect known as waterhammer which is due to the rapid rapid
opening or closing of valves or other changes in a liquid pipeline. The eect is due to the nite
speed of sound resulting from nite compressibility. The bulk modulus K represents the ratio of the
pressure change on a given material and the relative density change that it produces. Thus we have

K=

dp
d

(11.1)

K


(11.2)

The speed of sound c in any material is given by

s

c=

For an ideal incompressible uid, K ! 1, so that c ! 1. Another factor that is important is that
the material of which the wall of the pipe carrying liquid is not absolutely rigid. As an acoustic
wave travels through it, the high and low pressures cause not only a compression and expansion
of the liquid but also deformation of the pipe wall. This aects the overall compressibility of the
liquid-pipe system, and must be taken into account in the speed of sound. In a thin-walled pipe the
speed of sound is
s
K
c=
(11.3)
(1 + K
E )
where E is the Young's modulus of the material of the wall, and a is a constant which depends on
the kind of restraints on the pipe. As a special case E ! 1 for a rigid pipe material, which gives
equation (11.2). Table 11.1 gives the values of a for specic cases of pipe restraints.
where D and e are the diameter and wall thickness of the pipe, and  is the Poisson's ratio of
the material of the wall.
Consider now a long pipe of length L with a constant pressure liquid inlet at one end and a valve
at the other. Liquid ows through the pipe at velocity U and approximately uniform pressure p0 . At
137

Restraint

Anchored against longitudinal movement

Anchored at upper end
With frequent expansion joints

D (1
e


D
e

D (1:25
e

2 )
)

Table 11.1:  for specic restraints

time t = 0, the valve is suddenly closed leading to the progression of events schematically indicated
in Fig. 11.1; liquid velocity and pressure and pipe diameter are shown at dierent instants of time.
The liquid that is coming towards the valve builds up a pressure which compresses the liquid near
the valve and expands a portion of the pipe. The time for this compression to travel the length of
the pipe as a wave is L=c so that we will look at time intervals of L=2c. At t = L=2c, the valve half
of the pipe has expanded and the liquid has come to a stop there. The liquid in the inlet half is still
coming in. At t = L=c, the wave has reached the inlet; the entire pipe is expanded and the liquid
velocity is zero. The compression wave re ects back from the inlet; at t = 3L=2c the inlet half of
the pipe contracts to its normal size squeezing the liquid towards the inlet. At t = 2L=c, the entire
pipe is at normal size with liquid velocity towards the inlet. This liquid motion reduces the pipe
diameter so that at t = 5L=2c the valve half of the pipe is at less than normal size. In other words
the compression wave coming from the inlet has re ected from the valve as an expansion wave. At
t = 3L=c, this wave reaches the inlet and the entire pipe is at reduced area with zero liquid velocity.
At 7L=2c the expansion wave has re ected from the inlet. At t = 4L=c it reaches the valve. This
completes the cycle, with the whole process repeating itself until damped out.
The frequency of the pressure wave set up by the sudden closing of the valve in this example is
c=4L. Since 2L=c is the time taken for the compression wave to start from the valve and return, any
valve closing time smaller than this would cause waterhammer eects to be felt in the pipe. The
peak pressure rise is due to the conversion of the entire kinetic energy of the uid, i.e.
p = U 2 =2.

11.2 Two-phase systems

Multiphase ows in pipelines involving solids, liquids or gases are fairly common. Here we will look
at some aspects of two-phase liquid-gas ows.

11.2.1 Homogeneous two-phase mixture

Consider a homogeneous dispersion of liquid droplets in a gas, or gas bubbles in a liquid as shown
in Fig. 11.2. The speed of sound is given by equation (11.2) where the bulk modulus K is that of
the mixture. This has to be determined in terms of those of the pure liquid Kl and pure gas Kg . If
 is the volume fraction of gas, and 1  is that of the liquid, then the density  is given by
Since
V = V
p=K = V
p=Kg + (1 )V
p=Kl we have
1
(?)
+ 1Kl
1
1
=  1 
Kg + Kl g + (1 )l

K =
c2


Kg

138

(11.4)
(11.5)

@@

p = p0
V

V =0
p = p0 +
p

p = p0
V =0
p = p0 +
p

p = p0

t=0

p = p0
V

@@

t = L=2c

@@

t = L=c

V =0
p = p0 +
p @@

V =0

p = p0

p = p0
p

V =0
p = p0
p
V

t = 3L=2c

p = p0

V =0
p = p0
p

p = p0

@@

t = 2L=c

@@

t = 5L=2c

@@

t = 3L=c

@@

t = 7L=2c

@@

t = 4L=c

Figure 11.1: Sequence of events in waterhammer

=
where


g c2g

1
1
1 

+
(1
)l
+ l c2 g
l

Kg
g
Kl
c2l =
l
c2g , and for  ! 0, c2 ! c2l . We can take Kl
c2g =

In the limit of  ! 1, c2
approximations, so that

(11.6)

(11.7)
(11.8)

 Kg , l  g

as

K
1
c2 = g
l (1 )

(11.9)

 = g + (1 )l

(11.10)

139

b b bb

bb bb b bb b b b bb b b b b

Figure 11.2: Homogeneous two-phase ow

The minimum value is at  = 1=2

K
c2min = 2 g
l

(11.11)

11.2.2 Two-phase ows in horizontal pipes

The following are patterns in a horizontal pipe:
Bubble of froth ow
Plug ow
Stratied ow
Wavy ow
Slug ow
Annular ow
Spray or dispersed ow

dispersed bubbles of gas throughout liquid

alternate plugs of liquid and gas along upper part of pipe
liquid at bottom, gas on top with smooth interface
stratied with interfacial waves
frothy slug moves at velocity larger than liquid velocity
liquid lm, gas core
liquid as ne droplets

Flow pattern prediction

This can be done from the accompanying diagram (not yet included) where
l = (G =air )(L =water )
y = (water =L)[(L =water )]1=3 (water =L)2=3
G = m=A
_ is a mass velocity
Air and water properties are at 20C, 101.3 kPa.

Lockhart-Martinelli pressure drop prediction

To determine the pressure drop in a horizontal pipe, nd
1. Reynolds numbers ReL, ReG and ow regime (laminar or turbulent) for each single-phase
uid.


2. Pressure drops
pL ,
pG for each phase using
pi = f (L=D) 21 i Ui2 .
3. X =

pL=
pG .

4. YL or YG from graph
5.
pT P = YL
pL = YG
pG
140

Critical condition in ashing ow

As p decreases due to friction, the saturation temperature also decreases and the liquid vaporizes.
Below a certain backpressure pb the mass ow rate m_ does not change, This is the critical mass ow
rate m_ c and depends on the thermodynamics state of the uid and the ow regime. The critical or
maximum mass velocity is given in the graph1 (not yet included), where E = h01 + q = h02 , pc =
pressure p2 inside tube at 2 at critical condition, and x = vapor weight fraction (quality) at 2.

Problems

1. Saturated water ows in a 2.5 cm diameter boiler tube where there is vaporization due to pressure drop and
heat addition. The outlet quality and pressure are 0.8 and 2 atmospheres, respectively. What is the maximum
mass ow rate?
2. Water is owing through a cast iron pipe of length 10 m anchored against longitudinal movement. The wall
thickness is 1 mm and the pipe diameter is 1 cm. A valve at the end of the pipe is suddenly closed. How fast
must the pipe be closed for waterhammer eects to be felt? What is the frequency of the pressure wave? Use
a thin wall approximation.
3. Water ows at 5 m/s through a 30 cm diameter PVC
pipe of length 50 m and thickness 5mm. The modulus
of elasticity
and
Poisson's
ratio
of
PVC
are
2
:8  109 Pa and 0.45, respectively. The bulk modulus of water is
2:19  109 Pa. The pipe can be assumed to be anchored against longitudinal movement. There is a valve at
the end of the pipe that is being closed. (a) How slowly must this valve be closed so that waterhammer eects
are not produced? If the valve were closed suddenly, nd (b) the time it takes for the pressure wave to travel
upstream to the inlet of the pipe, (c) the frequency of oscillation of the pressure within the pipe, and (d) an
estimate of the pressure rise within the pipe.
4. A mixture of gaseous CO2 and liquid gasoline ows along a horizontal2 pipe of 100 m length, inner diameter 202
cm, and roughness 0.2 mm. The mass velocity of the CO2 is 10 kg/sm , and that of the gasoline is 250 kg/sm .
Find (a) the ow pattern of the two-phase ow, and (b) the pressure drop. Use the following property values:
(kg/m3 ) (Ns/m2 ) 4
(N/m)
Gasoline
680
2
:92  10
2
:16  10 2
Water (at 20 C,1 atm) 998
10 3
7:28  10 2
Air (at 20 C, 1 atm)
1.2
1:8  10 55
CO2
1.82 1:48  10
5. Air and water ow in a 50 cm diameter horizontal pipe at mass ow rates of 0.1 and 0.2 kg/s respectively.
determine the ow pattern.
6. Water under saturated conditions enters a long 10 cm diameter tube. Find the critical mass ow rate if the
exit quality and pressure are 0.4 and 700 kPa.
7. Crude oil5 (density2 800 kg/m3 , viscosity 60  10 5 Ns/m2) and natural gas (density 0.8 kg/m3 , viscosity
1:4  10 Ns/m ) ow along a 5 cm diameter, 300 m long pipe with roughness 0.15 mm. Find the pressure
drop for crude oil and natural gas ow rates of 0.64 kg/s and 0.04 kg/s.
Fluid

Density

Viscosity

Surface tension

Hodge, B.K., Analysis and Design of Energy Systems, Prentice-Hall, Englewood Clis, NJ, 1990.
141

142

Chapter 12

Numerical methods
In this chapter we discuss some of the basic ideas regarding numerical methods as they relate to the
solution of the ordinary and partial dierential equations which occur in uid dynamics.

12.1 Ordinary dierential equations

If x(t), y(t) and z (t) are the Cartesian coordinates of a uid particle moving in a velocity eld
ui + vj + wj, then

dx
= u(x; y; t)
(12.1)
dt
dy
= v(x; y; t)
(12.2)
dt
dz
= w(x; y; t)
(12.3)
dt
These equations can be integrated using any numerical procedure to determine the ow lines.

12.1.1 Fourth-order Runge-Kutta integration

This is a popular integration technique for ordinary dierential equations with given initial values.
Choose an appropriately small time step
t, for example
t = 0:01. Let

xi = x(i
t)
yi = y(i
t)
Then x0 ; y0 is the initial condition which is known. Knowing any xn ; yn we can nd xn+1 ; yn+1 in
the following manner. First determine the six numbers

p0 =
t u(xn ; yn ; t)
q0 =
t v(xn ; yn; t)
1
1
1
p1 =
t u(xn + p0 ; yn + q0 ; t +
t)
2
2
2
143

1
1
1
q1 =
t v(xn + p0 ; yn + q0 ; t +
t)
2
2
2
1
1
1
p2 =
t u(xn + p1 ; yn + q1 ; t +
t)
2
2
2
1
1
1
q2 =
t v(xn + p1 ; yn + q1 ; t +
t)
2
2
2
p3 =
t u(xn + p2 ; yn + q2 ; t +
t)
q3 =
t v(xn + p2 ; yn + q2 ; t +
t)
In the right-hand sides, the functions u(x; y; t) and v(x; y; t) are evaluated at the values of x, y and
t indicated in each expression.
Then
1
xn+1 = xn + (p0 + 2p1 + 2p2 + p3 )
6
1
yn+1 = yn + (q0 + 2q1 + 2q2 + q3 )
6
From x0 ; y0 get x1 ; y1 , then x2 ; y2 , and so on as long as necessary.

12.2 Partial dierential equations

There are dierent ways in which the incompressible ow problem in two dimensions with constant
properties may be written. As an example we will consider the problem of incompressible, twodimensional ow without body force.

12.2.1 Primitive variables

In nondimensional terms, we have
@u @v
+
= 0
@x @y#
"
#

@u
@u
@u
@p 1 @ 2 u @ 2 u
+u +v
=
+
+
@t
@x
@y
@x Re @x2 @y2
"

"

where

@v
@v
@v
+u +v
@t
@x
@y

@T
@T
@T
+u +v
@t
@x
@y
" 

1
@u
=
2
Re
@x

@p 1 @ 2 v @ 2 v
+
+
@y Re @x2 @y2

=
#

@v
+2
@y

2

(12.6)

(12.5)

1
@2T @2T
Ec
=
+ 2 + 
2
Re P r @x
@y
Re

2

(12.4)

@u @v
+
+
@y @x

(12.7)

2 #

12.2.2 Stream function-vorticity

The stream function (x; y) dened by

u =

@
@y

144

(12.8)

r
6
?r

n

-s

Figure 12.1: Boundary condition on vorticity

@
(12.9)
@x
satises the continuity equation (12.4). The vorticity, !, has only the z -component,  , which is
dened by
@v @u
(12.10)
=
@x @y
and which satises
@2
@2
+
+ =0
(12.11)
@x2 @y2
The vorticity equation
v =

@ @ @
+
@t @y @x

@ @
1 @2 @2
=
+
@x @y Re @x2 @y2

(12.12)

is obtained by taking the curl of the momentum equation.

Boundary conditions
At a solid boundary the value of may be specied. In Fig. 12.1 the value of  at a solid boundary,
for example w at the point w, can be related to the at an interior point, i , in the following way.
At the boundary, let s and n be the coordinates along and normal to it, respectively. Then, from a
Taylor series expansion, we have


1 @ 2 
@ 
2

n +
 (
n) + : : :
i= w+
@n w
2 @n2 w

(12.13)

where
n is the distance between w and i along the coordinate n normal to the wall. The normal
and tangential velocity components at the wall must vanish. The rst of these gives


@ 
=0
@s w

(12.14)

indicating that w must be a constant along the wall. The other velocity component gives


@ 
 =0
@n w
145

(12.15)

i; j + 1

i 1; j

i; j

i + 1; j

y

i; j

x
Figure 12.2: Computational mesh

Since

@ 2 
w =
@s2 w

@ 2 
@n2 w

(12.16)

and w is constant along the wall, we get from equation (12.13) that

w = 2 w 2 i
(
n)

(12.17)

12.3 Numerical methods

There are many numerical methods that have been developed for the governing equations of uid
mechanics: nite-dierence, nite-element, boundary-element, spectral methods, among others. At
this point we will describe only the rst.

12.3.1 Finite dierences

There are several variations of the nite dierence method, the following being one of them. With
reference to Fig. 12.2 he x-derivatives of a function (x; y) are written in the approximate forms
shown:
@ i;j i 1;j
=
backward dierence
(12.18)
@x

x
@ i+1;j i 1;j
=
central dierence
(12.19)
@x
2
x
2
@  i+1;j 2i;j + i 1;j
=
central dierence
(12.20)
@x2
(
x)2
The derivatives in the y-direction may be similarly written.
The steady state may be obtained by starting from certain initial conditions and integrating
forward in time. One way of doing this with the vorticity-stream function formulation is to rst
146

change equation (12.11) to a ctitious form

@
@2
@2
= 2 + 2 +
@t @x
@y

(12.21)

which reduces to the correct form under steady-state conditions. The time derivatives in this and
equation (12.12) are then approximated by

@i;j ki;j+1 ki;j

=
@t

t

(12.22)

where k + 1 and k refer to instants in time t +
t and t, respectively. The rest of the terms in these
equations are assumed to be at time t and are hence known. The integration process is halted when
the dependent variables are seen not to change very much with time.
Example 12.1

Compute the developing ow in a plane channel with uniform in ow.

?
Figure 12.3: Flow between at plates

The characteristic velocity and length are taken to be the uniform inlet velocity and the channel width,
respectively. In nondimensional terms they are both unity, while the length of the channel is `. Suitable
boundary conditions are
u = 1; v = 0 at x = 0
(12.23)
1
u = v = 0 at y = 
(12.24)
2
@u @v
= = 0 at x = `
(12.25)
@x @x
The rst condition refers to the uniform inlet velocity, the second to no slip at the walls, and the third to the
fully developed nature of the ow eld at the exit section.
Finite dierencing may be applied on a regular mesh that is obtained by dividing the region into N and
M parts in the x and y-directions, respectively, so that
x = `=N and
y = 1=M . Either the primitive or the
vorticity-stream function formulation may be employed. Symmetry around the y = 0 plane may be used, in
which case @u=@y = v = 0 at y = 0, instead of condition (12.24) at y = 1=2.

147

Problems

1. By numerical integration, determine where the uid particles originally within a rectangle of size 0:1  0:1
centered at (1,1) are after 50 time units. The velocity eld is two-dimensional:
u = y
v = x x3 + cos t
Take (a) = 0, and (b) = 1.
2. Solve the following problem numerically. A viscous uid enters a channel with a uniform velocity prole. The
channel is long enough for the ow to become fully developed. Using suitable geometrical and ow parameters,
compute the velocity eld.
Your report should be comprehensive; include details of how you set up the problem, description of the numerical
method that you used, graphs of the results that you obtained, and comparisons with available analytical
solutions.
3. In the problem above, determine the temperature eld if the entering uid and the walls of the channel are all
at dierent temperatures.
4. For the two-dimensional, unsteady velocity eld ui + vj, where
u = y
v = x x3 + cos t
determine, by numerical integration, where the uid particles originally within a square of size 0:1  0:1 centered
at (1,1) are after 50 time units. Consider two cases: (a) = 0, and (b) = 1. Take as initial conditions a
large number of dierent points within the square, say for example 100 evenly-spaced points, and integrate up
to t = 50. Store the nal coordinates in a le. Plot the initial as well as the nal states of all the points. The
particles should end up well-mixed (chaotically) for = 1, but not for = 0.

148

Appendix A

Governing equations
A.1 Integral form
Assumptions: nondeformable, inertial control volume
Conservation of mass

d
 dV +
V
n dA = 0
dt CV
CS
Newton second law for linear momentum
Z
Z
d
V dV +
V (V
n) dA = F
dt CV
CS
Newton second law for angular momentum

d
 (r  V) dV +
(r  V) (V
n) dA = T
dt CV
CS
First law of thermodynamics
Z

p
d
e dV +
e + (V
n) dA = Q_ + W_ s
dt CV

CS
Second law of thermodynamics
Z

d
1
s dV +
s (V
n) dA 
dt CV
CS
CS T
149

Q_

A dA

A.2 Dierential form

Conservation of mass

@
+ r
(u) = 0
@t
@ @ (ui )
+
= 0
@t
@xi
or

D
+ r
u = 0
Dt
D
@u
+ i = 0
Dt
@xi

Newton second law for linear momentum

Du
= r
 + f

Dt
Du
@ij
 j =
+ fj
Dt
@xi
Newton second law for angular momentum
= T
ij = ji


First law of thermodynamics

De
=
Dt
De

=
Dt


r
q_ +  : ru
@ q_i
@u
+ ij j
@xi
@xi

Second law of thermodynamics



Ds
@ q_i
 
Dt
@xi T
 
Ds
q_
  r

Dt
T
150

A.3 Cartesian coordinates

Velocity vector

u = ui + v j + w k
Material derivative

D
@
@
@
@
= +u +v +w
Dt @t
@x
@y
@z

Laplacian

@2
@2
@2
+
+
@x2 @y2 @z 2
Equations of motion (incompressible, Newtonian uid with constant properties)

@u @v @w
+ +
@x
@y @z

@u
@u
@u
@u

+u +v +w
@t
@x
@y
@z


@v
@v
@v
@v

+u +v +w
@t
@x
@y
@z


@w
@w
@w
@w

+u +v +w
@t
@x
@y
@z


@T
@T
@T
@T
+u +v +w
c
@t
@x
@y
@z
where
 = 2

"

@u
@x

2

@v
+
@y

2

@w
+
@z

2 #

+

"

= 0
=
=
=

(A.1)

@p
+ r2 u + fx
@x
@p
+ r2 v + fy
@y
@p
+ r2 w + fz
@z

(A.2)
(A.3)
(A.4)

= k r2 T + 

@u @v
+
@y @x

2

(A.5)

@u @w
+
+
@z @x

2

@v @w
+
+
@z @y

2 #

A.4 Cylindrical coordinates

Velocity vector
Material derivative
Laplacian

u = ur er + u e + uz ez
@
@ u @
@
D
= +u
+
+ uz
Dt @t r @r r @
@z


1@
@
1 @2
@2
r
+ 2 2+ 2
r @r @r
r @
@z
Equations of motion (incompressible, Newtonian uid with constant properties)

1@
1 @u @uz
(rur ) +
+
= 0
r
@r
r @
@z


@ur
@ur u @ur u2
@ur
@p

+ ur
+
+ uz
=
+  r2 ur
@t
@r
r @
r
@z
@r
151

ur
r2

2 @u
+ fr
r2 @

@u
@u u @u u u
@u
1 @p
u 2 @ur

+ ur  +   + r  + uz  =
+  r2 u
+
+ f
@t
@r
r @
r
@z
r @
r2 r2 @


@u u @u
@u
@p
@uz
+ ur z +  z + uz z =
+ r2 uz + fz

@t
@r
r @
@z
@z


@T
@T u @T
@T
c
+ ur +
+ uz
= k r2 T + 
@t
@r r @
@z
where
"


2 
 #
@ur 2 1 @u
@uz 2
 = 2
+ 2
+ ur +
@r
r @
@z
+ 

"

@u 1 @uz
+
@z r @

2

@uz @ur
+
+
@r
@z

2

2 #

1 @ur
@  u 
+
+r
r @
@r r

A.5 Spherical coordinates

Velocity vector

u = ur er + u e + ue

Material derivative

D
@
@ u @
u @
= +u
+
+ 
Dt @t r @r r @ r sin  @

Laplacian

@ 2@
1 @
@
1
@2
r = r12 @r
r
+ 2
sin 
+ 2 2
@r
r sin  @
@
r sin  @2
Equations of motion (incompressible, Newtonian uid with constant properties)
1 @
1 @u
1 @ 2

r ur +
(u sin ) +
= 0
2
r @r
r sin  @
r sin  @#
"
@ur
@u u @u
u @ur u2 + u2
@p
2ur

+ ur r +  r + 
=
+  r2 ur
@t
@r
r @ r sin  @
r
@r
r2
2

2u cot 
2 @u
+ fr
2
2
r
r sin  @


1 @p
u 2 @ur
2
+  r u
+
+ f
r @
r2 r2 @
@p
+ r2 uz + fz
@z

2 @u
r2 @

@u u @u u u
@u
@u
+ ur  +   + r  + uz  =

@t
@r
r @
r
@z


@uz
@uz u @uz
@uz
+ ur
+
+ uz
=

@t
@r
r @
@z


@T
@T u @T
@T
c
+ ur + 
+ uz
= kr2 T + 
@t
@r r @
@z
where
 = 2
+ 

"

"

@ur
@r

2

1 @u
+ 2
+ ur
r @

@u 1 @uz
+
@z r @

2

2

@uz
+
@z

@uz @ur
+
+
@r
@z
152

2

2 #


2 #

@  u 
1 @ur
+
+r
r @
@r r

Appendix B

Use of MATLAB
Matrix manipulation and plotting can be done using MATLAB1 which is available for several machines including Macintosh and UNIX systems. In a SPARCstation, enter MATLAB by typing
matlab

on the command line. The MATLAB prompt

>>

will be displayed.

B.1 Graphing
The following steps will enable you to read and plot a data le of points.

To load data the command is

>> load nameofdatale

where nameofdatale is a le containing your data.

To plot the data as points the command is

>> plot(nameofdatafile(:,1),nameofdatafile(:,2),'c')

The values in column one of your data will be plotted against the corresponding values in
column two. The character c species which character you will use to mark points with, and
must be between single quotes. Choices are
.
o
x
+
*

If this character is not included in the plot command, the points will be connected by a line.
1

See University of Notre Dame document number U4910,

153

Using MATLAB 4.2 on the SPARCstations

The commands to title your graph are

>> title('Your title')
>> xlabel('Your x-axis label')
>> ylabel('Your y-axis label')

To add more data to your graph you must rst type

>> hold

Otherwise a new graph will be made.

B.2 Plotting streamlines

As an example of stremline plotting we will consider ow around a Joukowski airfoil. The Joukowski
tranformation is given by
c2
z=+ ;

where the  -plane corresponds to the tranformed space. Once the solution is obtained in this space,
the inverse mapping is required to obtain the solution in the z -plane. To this end, the above equation
can be solved for  to give
r
 z 2
z
= 
c2 :
2
2
Since it is desired to have  ! z for large values of z , the positive root of the inverse mappling must
be chosen. This conclusion is valid for z expressed in polar form (i.e., z = exp i); however, care
must be taken when expressing z in Cartesian form (i.e., z = x + iy). For instance, for z = b where
b is a large real number, the above inverse mapping (with + root taken) predicts that   0; thus
the negative root should be used to preserve the desired mapping.
Therefore, when plotting the streamlines for the Joukowski airfoil in which z = x + iy, it is
necessary to change from the positive to the negative root in the above inverse transform for x < 0
in order to assure that   z for large z .
The following are hints for plotting the streamlines in MATLAB.

Generate a computational grid in the z-plane:

>> x = xmin : dx : xmax;
>> y = ymin : dy : ymax;
>> [X,Y] = meshgrid(x,y);
>> Z = X + i*Y;

Use the FOR and IF commands to create a loop in which  at each z = x + i  y location can
be computed. Remember to specify the ( ) root of the inverse transform for x < 0. (also,
type \help for" or \help if" at MATLAB prompt for information on how to create a loop).




Compute the complex potential F (z ).

Determine the values of the stream function at each (x; y) location:
>> psi = imag(F);

154

Contour plot the stream function values to obtain the streamlines (type \help contour" to get
online help information concerning this plot procedure):
>> contour(X,Y,psi,50)

155

156

Appendix C

Use of MAPLE
On the Sparcstations, type
to load MAPLE,

maple;
with(linalg);

to include the linear algebra packet, and

to get help on topic and examples.

help(topic);

157

158

Bibliography
Anderson, J.D., Jr., Modern Compressible Flow, McGraw-Hill, 2nd Ed., New York, 1990.
Aris, R., Vectors, Tensors, and the Equations of Fluid Mechanics, Dover, New York, 1962.
Batchelor, G.K., An Introduction to Fluid Dynamics, Cambridge University Press, 1967.
Chevray, R. and Mathieu, J., Topics in Fluid Mechanics, Cambridge University Press, Cambridge,
U.K., 1993.
Currie, I.G., Fundamental Mechanics of Fluids, 2nd Ed., McGraw-Hill, New York, NY, 1993.
Greenspan, H.P., The Theory of Rotating Fluids, Cambridge University Press, Cambridge, 1969.
Homann, K.A., Computational Fluid Dynamics for Engineers, Engineering Education System,
Austin, TX, 1989.
John, J.E.A., Gas Dynamics, Allyn and Bacon, Boston, 2nd Ed., 1984.
Kundu, P.K., Fluid Mechanics, Academic Press, San Diego, 1990.
Lamb, H., Hydrodynamics, Dover, New York, NY, 1945.
Landahl, M.T. and Mollo-Christensen, E., Turbulence and Random Processes in Fluid Mechanics,
Cambridge University Press, Cambridge, U.K., 1986.
Landau, L.D. and Lifshitz, E.M., Fluid Mechanics, Pergamon Press, Oxford, 1959.
Milne-Thompson, L.M., Theoretical Hydrodynamics, 5th Ed., Macmillan, London, 1968.
Saad, M.A., Compressible Flow, Prentice-Hall, Englewood Clis, NJ, 1985.
Schlichting, H., Boundary Layer Theory, 6th Ed., McGraw-Hill, New York, 1968.
Warsi, Z.U.A., Fluid Dynamics, Theoretical and Computational Approaches, CRC Press, Boca Raton, 1993.
White, F.M., Viscous Fluid Flow, McGraw-Hill, New York, NY, 1974.

159