B5.

3 Viscous Flow: Lecture Notes

Course synopsis

Overview

Viscous fluids are important in so many facets of everyday life that everyone has some intuition about the
diverse flow phenomena that occur in practice. This course is distinctive in that it shows how quite advanced
mathematical ideas such as asymptotics and partial differential equation theory can be used to analyse the
underlying differential equations and hence give scientific understanding about flows of practical importance,
such as air flow round wings, oil flow in a journal bearing and the flow of a large raindrop on a windscreen.

Reading list

[1] D.J. Acheson, Elementary Fluid Dynamics (Oxford University Press, 1990), chapters 2, 6, 7, 8. ISBN
0198596790.

[2] H. Ockendon & J.R. Ockendon, Viscous Flow (Cambridge Texts in Applied Mathematics, 1995). ISBN
0521458811.

Further reading

[3] G.K. Batchelor, An Introduction to Fluid Dynamics (Cambridge University Press, 2000). ISBN 0521663962.

[4] C.C. Lin & L.A. Segel, Mathematics Applied to Deterministic Problems in the Natural Sciences (Society for
Industrial and Applied Mathematics, 1998). ISBN 0898712297.

[5] L.A. Segel, Mathematics Applied to Continuum Mechanics (Society for Industrial and Applied Mathematics,
2007). ISBN 0898716209.

Synopsis (16 lectures)

Euler’s identity and Reynolds’ transport theorem. The continuity equation and incompressibility condition.
Cauchy’s stress theorem and properties of the stress tensor. Cauchy’s momentum equation. The incompressible
Navier-Stokes equations. Vorticity. Energy. Exact solutions for unidirectional flows; Couette flow, Poiseuille
flow, Rayleigh layer, Stokes layer. Dimensional analysis, Reynolds number. Derivation of equations for high and
low Reynolds number flows.

Thermal boundary layer on a semi-infinite flat plate. Derivation of Prandtl’s boundary-layer equations and sim-
ilarity solutions for flow past a semi-infinite flat plate. Discussion of separation and application to the theory of
flight.

Slow flow past a circular cylinder and a sphere. Non-uniformity of the two dimensional approximation; Os-
een’s equation. Lubrication theory: bearings, squeeze films, thin films; Hele–Shaw cell and the Saffman-Taylor
instability.

Authorship and acknowledgments

The author of these notes is Jim Oliver except for §1.1 and §1.3–1.6 which are modifications of lecture notes
of Peter Howell. The exposition follows closely the course text books [1,2]. Many thanks to Sandy Patel (MI)
for typing assistance. Please email comments and corrections to the course lecturer. All material in these notes
may be freely used for the purpose of teaching and study by Oxford University faculty and students. Other uses
require the permission of the authors.

1
Contents
1 The Navier-Stokes equations 3
1.1 Motivation for studying viscous fluids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 The summation convention and revision of vector calculus . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5 Newtonian constitutive law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.6 The Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.7 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.8 Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.9 Conservation of energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.10 Unidirectional flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.11 Dimensionless Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2 High Reynolds number flows 29

2.1 Thermal boundary layer on a semi-infinite flat plate . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Viscous boundary layer on a semi-infinite flat plate . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.3 Alternative derivation of the boundary layer equations . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4 Blasius’ similarity solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.5 Variable external flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.6 Breakdown of Prandtl’s theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3 Low Reynolds number flows 45

3.1 Slow flow past a circular cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 Slow flow past a sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.3 Lubrication Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.5 The Saffman-Taylor instability and viscous fingering . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.6 Thin films . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

2
1 The Navier-Stokes equations
1.1 Motivation for studying viscous fluids
• Fluid mechanics is the study of the flow of liquids and gases.

• In many practical situations the fluid can be described effectively as incompressible and inviscid, and
modelled by the Euler equations

∇ · u = 0, (1)
 
∂u
ρ + (u · ∇)u = −∇p + ρF, (2)
∂t

where the velocity u and pressure p are functions of position x and time t, ρ is the constant density and F
is the external body force acting per unit mass (e.g. gravity).

• For example:

(i) aerodynamic flows (e.g. flow past wings);

(ii) free surface flows (e.g. water waves).

• However, there are many fluid flow phenomena where inviscid theory fails, e.g.

(i) D’Alembert’s paradox states that there is no drag on an object moving steadily through a fluid (cf. a
ball bearing falling through oil).
(ii) The ability of a thin layer of fluid to support a large pressure, e.g. a lubricated bearing.
(iii) You can’t clean dust from a car by driving fast. Inviscid flow allows slip between car and air:

The tenacious dust suggests that fluid adjacent to the car is dragged along with it:

• The problem with inviscid fluid mechanics which gives rise to these failings, is that it takes no account of
the friction caused by one layer of fluid sliding over another or over a solid object.

• This friction is related to the stickiness, or viscosity, of real fluids.

• No fluid is completely inviscid (except liquid helium below around 1 Kelvin).

• Even for low-viscosity fluids (e.g. air), there will often be regions (e.g. thin boundary layer on a moving
car) where viscous effects are important.

3
 • In this course we will see how the Euler equations (1)–(2) must be modified to obtain the incompressible
Navier-Stokes equations
∇ · u = 0, (3)
 
∂u
ρ + (u · ∇)u = −∇p + µ∇2 u + ρF, (4)
∂t
for flows in which the viscosity µ is important.
• The incompressible Navier-Stokes equations (3)–(4) are nonlinear and in general extremely difficult to solve.
• We will use a combination of asymptotics and partial differential equation theory to analyse (3)–(4) and
hence give scientific understanding about flows of practical importance.

1.2 The summation convention and revision of vector calculus

• We will work in Cartesian coordinates Oxyz and let
i = (1, 0, 0) = e1 , j = (0, 1, 0) = e2 , k = (0, 0, 1) = e3 ,
denote the standard orthonormal basis vectors, so that a position vector may be written
x = x1 e1 + x2 e2 + x3 e3 ,
where (x1 , x2 , x3 ) are the Cartesian coordinates.
• We employ the summation convention of summing over all possible repeated indices in an expression.
• An index which is summed in this way is called a dummy index.
• The summation convention should only be used if it is clear from the context over what ranges the dummy
indices should be summed.
• Example 1: Denote by (u, v, w) = (u1 , u2 , u3 ) the components of the fluid velocity, so that
3
X
u = (u, v, w) = ui + vj + wk = ui ei = ui ei .
i=1

• Example 2: Kronecker’s delta δij is defined by


1 if i = j,
δij =
0 6 j,
if i =
while the Levi-Civita symbol εijk is defined by


 1 if i, j, k in cyclic order,

εijk = ei · (ej ∧ ek ) = −1 if i, j, k in acyclic order,



0 otherwise.
In the identity
εijk εrsk = δir δjs − δis δjr ,
the sum is over k = 1, 2, 3.
• Example 3: The determinant of a matrix A = {aij }3×3 is given by

a11 a12 a13

det(A) = a21 a22 a23 = εijk a1i a2j a3k .

a31 a32 a33

4
• Example 4: The scalar and vector product of two vectors a = ai ei and b = bi ei are given by

a · b = a1 b1 + a2 b2 + a3 b3 = ai bi ,

e1 e2 e3

a∧b = a1 a2 a3 = εijk ei aj bk .

b1 b2 b3

• Example 5: For a differentiable scalar field f (x) and a differentiable vector field G(x) = Gi (x)ei ,
∂f ∂Gj ∂Gk
∇f = ei , ∇·G= , ∇ ∧ G = εijk ei ,
∂xi ∂xj ∂xj
∂f ∂2 f
(u · ∇) f = ul , ∇2 f = .
∂xl ∂xm ∂xm
• Example 6: Since εijk = ei · (ej ∧ ek ),
∂Gk
∇ ∧ G = (ei · (ej ∧ ek )) ei
∂xj
∂Gk
= ej ∧ ek
∂xj
∂
= ej ∧ (Gk ek )
∂xj
∂G
= ej ∧ .
∂xj

• Example 7: The identities of vector calculus may be readily derived using these definitions and vector
identities. E.g. for a differentiable vector field u,
 
∂ ∂u
∇ ∧ (∇ ∧ u) = ei ∧ ej ∧
∂xi ∂xj
∂2
= (ei ∧ (ej ∧ u))
∂xi ∂xj
∂2
= ((ei · u)ej − (ei · ej )u)
∂xi ∂xj
∂2
= (ui ej − δij u)
∂xi ∂xj
 
∂ ∂ui ∂2u
= ej −
∂xj ∂xi ∂xj ∂xj
= ∇ (∇ · u) − ∇2 u.

• Example 8: (The divergence theorem) Let the region V in R3 be bounded by a piecewise smooth surface
∂V with outward pointing unit normal n = nj ej . Let G(x) = Gj (x)ej be a differentiable vector field on
V . Then ZZ ZZZ ZZ ZZZ
∂Gj
G · n dS = ∇ · G dV or Gj nj dS = dV.
∂V V ∂V V ∂xj

• Example 9: The incompressible Navier-Stokes equations (3)–(4) may be written in the form
 
∂uj ∂ui ∂ui ∂p ∂ 2 ui
= 0, ρ + uj =− +µ + ρFi (i = 1, 2, 3).
∂xj ∂t ∂xj ∂xi ∂xj ∂xj

5
1.3 Kinematics
1.3.1 The continuum hypothesis
• Liquids and gases consist of atoms or molecules, which move around and interact with each other and with
obstacles.
• At a macroscopic level, the net result of all these random interactions is that the fluid appears to be a
continuous medium, or continuum.
• The continuum hypothesis is the assumption that the fluid can be characterized by properties (e.g. density
ρ, velocity u, pressure p, absolute temperature T ) which depend continuously on position x and time t
(rather than having to keep track of a large number of individual atoms or molecules).
• The hypothesis holds for the vast majority of practically important flows, but can break down in extreme
conditions (e.g. very low density).

1.3.2 Eulerian and Lagrangian coordinates

• We distinguish two spatial coordinate systems, as follows.
Eulerian Coordinates x = (x1 , x2 , x3 )
– Label points fixed in space.
– Fluid properties at each point x change as different fluid particles pass through that point,
e.g. u(x, t) is fluid velocity at point x at time t.
– The Eulerian time derivative (i.e. holding x fixed) is denoted by
∂ ∂
≡ .
∂t ∂t x
Lagrangian Coordinates X = (X1 , X2 , X3 )
– Label fluid particles and in this sense they “move with the fluid.”
– Fluid properties are described for each fluid particle as it moves through different points in space.
– The convective or material or Lagrangian time derivative (i.e. holding X fixed) is denoted by
D ∂
≡ .
Dt ∂t X
• We choose the label X to be the initial position of a fluid particle at time t = 0, and denote by x(X, t) its
position at time t ≥ 0, i.e.
∂
x(X, 0) = X, x(X, t) = u(x(X, t), t)
∂t X
so that {x(X, t) : t ≥ 0} is the pathline of the fluid particle at X at t = 0.
• The convective derivative D/Dt is related to the Eulerian time derivative ∂/∂t using the chain rule. For a
differentiable scalar field f (x, t),
Df ∂
= f (x(X, t), t)
Dt ∂t X

∂f ∂x1 ∂f ∂x2 ∂f ∂x3 ∂f

= + + +
∂x1 ∂t X ∂x2 ∂t X ∂x3 ∂t X ∂t
∂f ∂x
= + · ∇f
∂t ∂t X

∂f
= + (u · ∇)f
∂t
 
∂
= + u · ∇ f.
∂t

6
 • The convective derivative may be written in the form
D ∂ ∂ ∂ ∂ ∂ ∂ ∂
= +u·∇= +u +v +w = + uj
Dt ∂t ∂t ∂x ∂y ∂z ∂t ∂xj

and applied to vector quantities, e.g. the acceleration of a fluid particle is

 
Du ∂u ∂ui ∂ui
= + (u · ∇)u = ei + uj .
Dt ∂t ∂t ∂xj

• Both Eulerian and Langrangian coordinates can be useful when describing fluid motion:

Eulerian coordinates ⇒ working in a fixed laboratory frame

⇒ convenient for calculations;

Langrangian coordinates ⇒ working in the frame-of-reference of a moving fluid particle

⇒ covenient for the application of conservation principles
(e.g. mass, momentum, energy).

1.3.3 The Jacobian and Euler’s identity

• The continuum hypothesis implies that there is a one-to-one relation between X and x(X, t), i.e. fluid can
never appear from nowhere or disappear.

• We assume in addition that the map from X to x(X, t) is continuous, so that the Jacobian

∂x1 ∂x1 ∂x1

∂X1 ∂X2 ∂X3
∂(x1 , x2 , x3 ) ∂x2 ∂x2 ∂x2
J(X, t) = ≡
∂(X1 , X2 , X3 ) ∂X1 ∂X2 ∂X3
∂x3 ∂x3 ∂x3
∂X1 ∂X2 ∂X3
is positive and bounded.

• The Jacobian J(X, t) measures the change in a small volume compared with its initial volume, i.e.

dx1 dx2 dx3 = J(X, t) dX1 dX2 dX3 ,

with J(X, 0) = 1 because x(X, 0) = X by definition.

• The Jacobian may be written in the form (using the summation convention here and hereafter)

∂x1 ∂x2 ∂x3

J(X, t) = εijk ,
∂Xi ∂Xj ∂Xk

where εijk is the Levi-Civita symbol.

• Hence, the rate of change of the Jacobian J following the fluid is given by
 
DJ D ∂x1 ∂x2 ∂x3
= εijk
Dt Dt ∂Xi ∂Xj ∂Xk
      
∂ Dx1 ∂x2 ∂x3 ∂x1 ∂ Dx2 ∂x3 ∂x1 ∂x2 ∂ Dx3
= εijk + +
∂Xi Dt ∂Xj ∂Xk ∂Xi ∂Xj Dt ∂Xk ∂Xi ∂Xj ∂Xk Dt
 
∂u1 ∂x2 ∂x3 ∂x1 ∂u2 ∂x3 ∂x1 ∂x2 ∂u3
= εijk + +
∂Xi ∂Xj ∂Xk ∂Xi ∂Xj ∂Xk ∂Xi ∂Xj ∂Xk

7
  
∂u1 ∂xm ∂x2 ∂x3 ∂x1 ∂u2 ∂xm ∂x3 ∂x1 ∂x2 ∂u3 ∂xm
= εijk + +
∂xm ∂Xi ∂Xj ∂Xk ∂Xi ∂xm ∂Xj ∂Xk ∂Xi ∂Xj ∂xm ∂Xk
∂u1 ∂(xm , x2 , x3 ) ∂u2 ∂(x1 , xm , x3 ) ∂u3 ∂(x1 , x2 , xm )
= + +
∂xm ∂(X1 , X2 , X3 ) ∂xm ∂(X1 , X2 , X3 ) ∂xm ∂(X1 , X2 , X3 )
 
∂u1 ∂u2 ∂u3 ∂(x1 , x2 , x3 )
= + + ,
∂x1 ∂x2 ∂x3 ∂(X1 , X2 , X3 )
where in the second line we used the fact εijk is constant and the product rule; in the third line we used
the definition
Dxl
= ul ;
Dt
in the fourth line we used the chain rule to write
∂ul ∂ul ∂xm
= ;
∂Xn ∂xm ∂Xn
and in the sixth line we used the fact a determinant is zero if it has repeated rows.
• We have therefore derived Euler’s identity
DJ
= J∇ · u. (5)
Dt
1.3.4 Reynolds’ transport theorem
• Consider a material volume V (t) which is transported along by the fluid and bounded by the surface ∂V (t)
with outward unit normal n.

V (0) n
Pathline: V (t)
Dx
= u(x, t)
Dt x(X, t)
X
∂V (t)
∂V (0)

• Let f (x, t) be any continuously differentiable property of the fluid, e.g. density, kinetic energy per unit
volume. The total amount of f inside V (t) is given by the volume integral
ZZZ
I(t) = f (x, t) dx1 dx2 dx3 .
V (t)

• Transforming from Eulerian coordinates x to Langrangian coordinates X, the integral becomes

ZZZ
I(t) = f (x(X, t), t)J(X, t) dX1 dX2 dX3 .
V (0)

• We can now “differentiate under the integral sign” to find that the time rate of change of I(t) is given by
ZZZ
dI ∂
= (f J) dX1 dX2 dX3
dt V (0) ∂t X
ZZZ
D
= (f J) dX1 dX2 dX3
V (0) Dt
ZZZ  
Df DJ
= J +f dX1 dX2 dX3 .
V (0) Dt Dt
ZZZ  
Df
= + f ∇ · u J dX1 dX2 dX3
V (0) Dt
ZZZ  
∂f
= + ∇ · (f u) dx1 dx2 dx3 ,
V (t) ∂t

8
 where on the fourth line we used Euler’s identity (5) and on the last line we set
Df ∂f ∂f
+ f∇ · u = + (u · ∇)f + f ∇ · u = + ∇ · (f u)
Dt ∂t ∂t
while transforming back to Eulerian coordinates.

• Thus, we have proven Reynolds’ transport theorem: If V (t) is a material volume convected with velocity
u(x, t) and f (x, t) is a continuously differentiable function, then
ZZZ ZZZ
d ∂f
f dV = + ∇ · (f u) dV. (6)
dt V (t) V (t) ∂t

1.3.5 Visualization of Reynolds’ transport theorem

• Let δt be a small time increment, then
ZZZ ZZZ !
I(t + δt) − I(t) 1
= f (x, t + δt) dV − f (x, t) dV
δt δt V (t+δt) V (t)

ZZZ ZZZ
f (x, t + δt) − f (x, t) 1
= dV + f (x, t)dV
V (t) δt δt V (t+δt)\V (t)
| {z } | {z }
Change inside V (t) Change due to moving boundary

• The volume V (t + δt)\V (t) is a thin shell around ∂V (t). The amount of f swept through a surface element
δS of ∂V (t) in the time increment δt is f (u · n δt)δS:

n
∂V (t + δt)
u · n δt
δS
∂V (t)

f × parcel volume = f (u · n δt) δS

• Hence, as δt → 0, ZZZ ZZ
I(t + δt) − I(t) ∂f
→ dV + f u · n dS.
δt V (t) ∂t ∂V (t)

• Finally, apply the divergence theorem

ZZ ZZZ
f u · n dS = ∇ · (f u) dV
∂V (t) V (t)

to recover Reynolds’ transport theorem (6).

1.3.6 Conservation of mass

• Since a material volume V (t) always consists of the same fluid particles, its mass must be preserved, i.e.
ZZZ
d
ρ dV = 0.
dt V (t)

• Apply Reynolds’ transport theorem (6) with f = ρ to obtain

ZZ
∂ρ
+ ∇ · (ρu) dV = 0.
V (t) ∂t

9
 • Since the volume V (t) is arbitrary, the integrand must be zero (assuming it is continuous), i.e.

∂ρ
+ ∇ · (ρu) = 0. (7)
∂t

• This equation is called the continuity equation and represents pointwise conservation of mass.

• Note that an application of Reynolds’ transport theorem (6) with f = ρF implies that
ZZZ ZZZ
d ∂
ρF dV = (ρF ) + ∇ · (ρF u) dV
dt V (t) V (t) ∂t
ZZZ    
∂ρ ∂F
= F + ∇ · (ρu) +ρ + (u · ∇)F dV,
V (t) ∂t ∂t
| {z } | {z }
=0 DF
=
Dt

which yields the following useful corollary

ZZZ ZZZ
d DF
ρF dV = ρ dV (8)
dt V (t) V (t) Dt

for continuously differentiable ρ and F .

• For incompressible fluid, Dρ/Dt = 0, and hence by (7) we obtain the incompressibility condition

∇ · u = 0. (9)

1.4 Dynamics
1.4.1 The stress vector

• Consider a surface element δS with unit normal n drawn through x in the fluid:

n = ej n j
x2 δS
x
e2
O
e3 e1 x1
x3

• The stress vector t(x,t,n) is the force per unit area (i.e. stress) exerted on the surface element by the fluid
toward which n points.

• Example: Fluid flows inside a rectangular box R = {x : 0 < xj < j for j = 1, 2, 3} whose boundary ∂R is
rigid and solid. The force per unit area exerted by the fluid at a point x = x1 e1 + x2 e2 on the bottom face
x3 = 0 is t(x1 e1 + x2 e2 , t, n = e3 ), so the total force exerted by the fluid on the bottom face is given by
Z 2Z 1
t(x1 e1 + x2 e2 , t, e3 ) dx1 dx2 .
0 0

10
1.4.2 Conservation of momentum
• Consider a material volume V (t) with boundary ∂V (t) whose outward unit normal is n:

n
V (t)

∂V (t)

ZZZ
• Its linear momentum is ρudV .
V (t)

• Forces acting on V (t):

(i) Internal forces represented by the stress t(x, t, n) exerted by the fluid outside V (t) on the fluid inside
V (t) via the boundary ∂V (t).
(ii) External forces (e.g. gravity, EM) represented by a body force F(x, t) acting per unit mass.

• Newton’s second law for the material volume V (t) states that the time rate of change of its linear momentum
is equal to the net force applied. Thus,
ZZZ ZZ ZZZ
d
ρu dV = t dS + ρF dV. (10)
dt V (t) ∂V (t) V (t)

Example: Derivation of Euler’s momentum equation

• For an inviscid fluid the stress vector
t = −p n,
where p is the pressure.

• Note the implications:

(i) stress is purely in the normal direction (i.e. no friction);

(ii) the magnitude of the stress (i .e. p) is independent of the orientation of the surface element (i.e. of n).

• For viscous fluid neither of these is correct: we must allow for stress which is not necessarily in the normal
direction, and whose magnitude depends on n.

• Note that the corollary to Reynolds’ transport theorem (8) implies that
ZZZ ZZZ
d Du
ρu dV = ρ dV,
dt V (t) V (t) Dt

while the divergence theorem implies that

ZZ ZZ ZZ ZZ ZZZ
∂
t dS = −pn dS = −ei pδij nj dS = −ei (pδij ) dS = −∇p dV.
∂V (t) ∂V (t) ∂V (t) ∂V (t) ∂xj V (t)

• Hence, by (10), ZZZ

Du
ρ + ∇p − ρF dV = 0;
V (t) Dt
since V (t) is arbitrary, the integrand must be zero (assuming it is continuous), and we recover Euler’s
equation
Du
ρ = −∇p + ρF.
Dt

11
 • In order to generalize this methodology to any continuous medium (and, in particular, a viscous fluid), it
is necessary to convert the surface integral
ZZ
t dS
∂V (t)

into a volume integral. This is accomplished via Cauchy’s stress theorem, which recasts the stress vector
in a form amenable to the divergence theorem.

1.4.3 The stress tensor

• The stress tensor σij (x, t) is the component of the stress in the xi -direction exerted on a surface element
with normal in the xj -direction by the fluid toward which ej points.

• Note that the subscript i corrsponds to the direction of the stress, while the subscript j corresponds to the
direction of the normal. Moreover, by definition,

σij (x, t) = ei · t(x, t, ej ),

so
t(x, t, ej ) = ei σij (x, t).

• Example 1: Fluid flows in the upper half-space x3 > 0 above a rigid solid plate at x3 = 0. The force per
unit area exerted by the fluid on the plate at a point x = x1 e1 + x2 e2 is given by

t(x1 e1 + x2 e2 , t, n = e3 ) = ei σi3 (x1 e1 + x2 e2 , t);

σ13 and σ23 are shear stresses, while σ33 is a normal stress.

• Example 2: Since t = −pn for an inviscid fluid,

σij (x, t) = ei · t(x, t, ej ) = ei · (−p(x, t)ej ) = −p(x, t) δij ,

where δij is Kronecker’s delta.

1.4.4 Action and reaction

• Consider a material volume V (t) having at time t the configuration of a right circular cylinder, with radius
R, height R, centre x and outward unit normals as shown.

∂V1 (t)
n1

!R
R
x ∂V3 (t)

∂V2 (t) n3
n2

• Newton’s second law for the material volume (10) may be written in the form
ZZZ ZZ
Du
ρ − ρF dV = t dS.
V (t) Dt ∂V (t)

12
 • Assuming the integrand is continuous (so that, in particular, the acceleration and body force are finite),
the integral mean value theorem implies that
ZZZ
Du
ρ − ρF dV = O(R3 ) as R → 0.
V (t) Dt

• Moreover, as , R → 0,
ZZ 3 ZZ
X
t(x, t, n) dS = t(xj , t, nj ) dS = t(x, t, n1 )πR2 + t(x, t, n2 )πR2 + O(R2 , R3 ).
x∈∂V (t) j=1 xj ∈∂Vj (t)

• Combining these expressions gives

(t(x, t, n1 ) + t(x, t, n2 )) R2 = O(R2 , R3 ) as , R → 0.

• Since this expression pertains for arbitrarily small  and R, we deduce (setting n = n1 = −n2 )

t(x, t, −n) = −t(x, t, n), (11)

which is Newton’s third law (action-reaction) for a continuous medium.

1.4.5 Cauchy’s stress theorem

• Consider a material volume V (t) having at time t the configuration of a small tetrahedron as shown. Let
the slanting face have area A = L2 and outward unit normal n = ej nj , with nj > 0.

n
e3
A = L2

A2 A1
x
O e2 A3

e1

• Newton’s second law for the material volume (10)

ZZZ ZZ
Du
ρ − ρF dV = t dS.
V (t) Dt ∂V (t)

• Assuming the integrand is continuous, the integral mean value theorem implies that
ZZZ
Du
ρ − ρF dV = O(L3 ) as L → 0.
V (t) Dt

• Since the face with area Aj = nj A = nj L2 (by Q1(b)) has outward unit normal −ej and the slanted face
with area A = L2 has outward unit normal n,
ZZ
t dS = (t(x, t, n) + t(x, t, −ej )nj ) L2 + O(L3 ) as L → 0.
∂V (t)

• Combining these expressions and using Newton’s third law (11) gives

(t(x, t, n) − t(x, t, ej )nj ) L2 = O(L3 ) as L → 0.

13
 • This expression pertains for arbitrarily small L, so there is a local equilibrium of the surface stresses, with

t(x, t, n) = t(x, t, ej )nj .

• Note that this expression holds for an arbitrarily oriented unit normal n by a straightforward generalization
of the above argument.

• Finally, since t(x, t, ej ) = ei σij (x, t) by definition, we deduce Cauchy’s stress theorem

t(n) = ei σij nj (12)

where we have suppressed the dependence of t and σij on x and t.

• Thus, knowing the nine quantities σij we can compute the stress in any direction.

1.4.6 Cauchy’s momentum equation

• We return to the conservation of momentum of a material volume (10).

• Recall that the corollary to Reynolds’ transport theorem (8) implies that
ZZZ ZZZ
d Du
ρu dV = ρ dV.
dt V (t) V (t) Dt

• Using Cauchy’s Stress theorem (12), the net surface force is

ZZ ZZ ZZZ
∂σij
t(n) dS = ei σij nj dS = ei dV
∂V (t) ∂V (t) V (t) ∂xj

after an application of the divergence theorem.

• Combining these expressions we find that (10) may be written in the form
ZZZ
Du ∂σij
ρ − ei − ρF dV = 0.
V (t) Dt ∂xj

• Since V (t) is arbitrary, the integrand must be zero (if it is continuous), and we deduce Cauchy’s momentum
equation
Du ∂σij
ρ = ei + ρF, (13)
Dt ∂xj
which holds for any continuum, not just a fluid.

1.4.7 Symmetry of the stress tensor

• For a material volume V (t), conservation of angular momentum about the origin O is given by
ZZZ ZZ ZZZ
d
x ∧ ρu dV = x ∧ t dS + x ∧ ρF dV. (14)
dt V (t) ∂V (t) V (t)

• Applying Reynolds’ transport theorem and the divergence theorem to this expression gives
ZZZ   ZZZ
Du ∂σij
x∧ ρ − ei − ρF dV = ej ∧ ei σij dV.
V (t) Dt ∂xj V (t)

• We can then deduce from Cauchy’s momentum equation (13) that

ZZZ
ej ∧ ei σij dV = 0.
V (t)

14
 • Since V (t) is arbitrary, the integrand must be zero (if it is continuous), i.e.
0 = ej ∧ ei σij = e1 (σ32 − σ23 ) + e2 (σ13 − σ31 ) + e3 (σ21 − σ12 ) ,
which implies that the stress tensor is symmetric, i.e.
σij = σji .

• The stress tensor may also be shown to be symmetric by taking V (t) to be instantaneously a vanishingly
small cube in (14) and estimating the various terms as in the derivations of Newton’s third law (11) and
Cauchy’s stress theorem (12).
• Note that the symmetry of the stress tensor σij implies that it consists of six independent quantities only,
namely σ11 , σ22 , σ33 , σ12 = σ21 , σ13 = σ31 and σ23 = σ32 .

1.4.8 Change of coordinate system

• Suppose we rotate the coordinate system from Ox1 x2 x3 with orthonormal basis vectors {e1 , e2 , e3 } to
Ox01 x02 x03 with orthonormal basis vectors {e01 , e02 , e03 }.
• A position vector r may be written r = xj ej = x0i e0i , so the rotation of the coordinate system transforms
the coordinates of the vector r according to
x0i = r · e0i = (xj ej ) · e0i = lij xj , (15)
where lij = e0i · ej .
• Equally, we can write r = xi ei = x0j e0j and deduce that the inverse transformation is given by
xi = r · ei = (x0j e0j ) · ei = lji x0j . (16)

• Combining (15) and (16) gives

xi = lji x0j = lji ljk xk ≡ δik xk ,
where δik is Kronecker’s delta.
• Since the last expression holds for arbitrary xk , we deduce that
lji ljk = δik ,
i.e. the matrix L = {lij }3×3 is orthogonal, with
LLT = I = LT L,
where I = {δij }3×3 is the 3-by-3 identity matrix.
• By definition the stress tensor in the primed frame is
0
σrs = e0r · t(e0s ).

• Writing Cauchy’s stress theorem (12) in the original frame in the form
t(n) = ei σij nj = ei σij (n · ej ),
we deduce that under the rotation of the coordinate system due to the orthogonal matrix L = {lij }3×3 , the
stress tensor transforms according to
0


σrs = e0r · ei σij (e0s · ej ) = lri lsj σij , (17)
or equivalently
S 0 = LSLT ,
where
S 0 = {σrs
0
}3×3 , S = {σij }3×3 ,

15
 • That σij transforms according to (17) means that it is a second-rank tensor, which are a generalization of
vector fields or first-rank tensors, cf. (15) and (17).
• It is for this reason that σij is called the stress tensor and the upshot of the above analysis is that there is
an invariantly defined stress in the fluid.
• Note that σij is a second-rank tensor if and only if Cauchy’s stress theorem is independent of the choice of
of coordinate system, i.e. (17) holds iff
t0 (n0 ) = e0r σrs
0
n0s ,
where t0 = t0r e0r , n0 = n0s e0s , with t0r = lri ti , n0s = lsj nj .

1.5 Newtonian constitutive law

1.5.1 Recap
• We have derived for a continuous medium expressions (7) and (13) representing conservation of mass and
momentum, viz.
∂ρ Du ∂σij
+ ∇ · (ρu) = 0, ρ = ei + ρF.
∂t Dt ∂xj

• The number of scalar equations (1 + 3 = 4) is less than the number of unknowns (ρ, u, σij = σji ,
i.e. 1 + 3 + 6 = 10), so we need more information to close the system.

1.5.2 Constitutive relations

• To make progress we must decide how the stress tensor σij depends on the pressure p and velocity u.
• This is called a constitutive relation and cannot be deduced, relying instead on some assumptions about
the physical properties of the material under consideration.
• For example, we would expect the constitutive relation for a solid to be quite different from that for a fluid.
• Examples of simple constitutive relations:
(i) Hooke’s law for the extension of a spring;
(ii) Fourier’s law for the flux of heat energy down the temperature gradient;
(iii) the inviscid stress tensor σij = −pδij .
• Note that
(i) a “thought-experiment” suggests these laws are reasonable;
(ii) they could be confirmed experimentally;
(iii) they will almost certainly fail under “extreme” conditions.

Example: Heat conduction in a stationary isotropic continuous medium

• Let T (x, t) be the absolute temperature in a stationary isotropic continuous medium (e.g. a fluid or a rigid
solid at rest), with constant density ρ and specific heat cv .
• Let q(x, t) be the heat flux vector, so that q · n is the rate of transport of heat energy per unit area across
a surface element in the direction of its unit normal n.
• For a fixed region V in the medium with boundary ∂V whose outward unit normal is n, conservation of
heat energy is given by ZZZ ZZ
d
ρcv T dV = q · (−n) dS,
dt V ∂V
where the term on the left-hand side (LHS) is the rate of increase of internal heat energy and the term on
the right-hand side (RHS) is the rate of heat conduction into V across ∂V .

16
 • Differentiating under the integral sign on the LHS and applying the divergence theorem on the RHS gives
ZZZ  
∂T
ρcv + ∇ · q dV = 0.
V ∂t

• Since V is arbitrary, the integrand must be zero (if it is continuous), i.e.

∂T
ρcv + ∇ · q = 0.
∂t

• A closed model for heat conduction is obtained by prescribing a constitutive law relating the heat flux
vector q to the temperature T .

• Fourier’s Law states that heat energy is transported down the temperature gradient, with

q = −k∇T,

where k is the constant thermal conductivity.

• Hence, T satisfies the heat or diffusion equation

∂T
= κ∇2 T,
∂t
where the thermal diffusivity κ = k/ρcv .

• The SI units of the dependent variables and dimensional parameters are summarized in the following table;
note that kelvin K is the SI unit of temperature, joule J is the SI derived unit of energy (1J = 1N m) and
the newton N is the SI derived unit of force (1N = 1 Kg m s−1 ).

Quantity Symbol SI units

Temperature T K
Heat flux vector q J m−2 s−1
Density ρ Kg m−3
Specific heat cv J Kg−1 K−1
Thermal conductivity k J m−1 s−1 K−1
Thermal diffusivity κ m2 s−1

1.5.3 The Couette flow rheometer

• A layer of viscous fluid of height h is sheared between two parallel plates by moving the top plate horizontally
with speed U .

• The force required to maintain the motion of the top plate is proportional to U and inversely proportional
h.

• Thus, the shear stress exerted by the top plate on the fluid must satisfy
U
σ12 |y=h ∝ .
h

• The liquid flows parallel to the plates, with a velocity profile that is linear in y.

• These observations suggest a constitutive law of the form

∂u
σ12 = µ , (18)
∂y
where µ is a constant of proportionality that depends on the liquid only (in fact µ is the dynamic viscosity).

17
 Ui
y=h

Uy
u= i
h
y=0

1.5.4 The Newtonian constitutive law

• To generalize (18), we begin by writing
σij = −pδij + τij , (19)
where −pδij is the inviscid stress tensor and τij is the deviatoric stress tensor, due to the presence of
viscosity.
• Experiments suggest that
(A) τij is a linear function of the velocity gradients ∂uα /∂xβ ;
(B) the relation between τij and the velocity gradients is isotropic, i.e. invariant to rotations of the
coordinate axes (so that there is no preferred direction).
• These conditions define a Newtonian fluid because they are sufficient to determine the form of τij completely:
together with symmetry of σij , (1)–(2) imply that
τij = λ(∇ · u)δij + 2µeij , (20)
where λ is the bulk viscosity, µ is the dynamic (shear) viscosity and
 
1 ∂ui ∂uj
eij = +
2 ∂xj ∂xi
is the rate-of-strain tensor (which is zero for any rigid-body motion, i.e. if there is no deformation of fluid
elements).
• The expression (20) is the constitutive law for a Newtonian fluid.

Outline Proof (not examinable)

• By property (A),
∂uα
τij = Aijαβ ,
∂xβ
where Aijαβ are constants - 81 of them!
• Since τij and ∂uα /∂xβ are rank-2 tensors, tensor theory implies that Aijαβ is a rank-4 tensor.
• Property (B) means that Aijαβ is an isotropic tensor, i.e. if τij0 = A0ijαβ ∂u0α /∂x0β under rotation of the
coordinate system due to the orthogonal matrix L = {lij }3×3 in §1.4.8, then A0ijαβ = Aijαβ .
• Since Aijαβ is an isotropic rank-4 tensor, tensor theory implies that
Aijαβ = λδij δαβ + µ(δiα δjβ + δiβ δjα ) + µ† (δiα δjβ − δiβ δjα ),
where µ, µ† and λ are constants - just 3 of them!
• Since τij is symmetric, Aijαβ = Ajiαβ , which implies µ† = 0, and hence that
∂uα
τij = Aijαβ = λeαα δij + 2µeij .
∂xβ

• Finally, note that eαα = ∂uα /∂xα = ∇ · u.

18
1.5.5 Incompressibility assumption
• Except where stated we will assume that the density ρ is constant, so that

– the flow is incompressible;

– the continuity equation (7) is replaced by the incompressibility condition (9);
– the Newtonian constitutive law (20) becomes
 
∂ui ∂uj
τij = 2µeij = µ + , (21)
∂xj ∂xi

so that the bulk viscosity λ drops out of the model.

• Note that

– most liquids are virtually incompressible except at extremely high pressures;

– most gases are compressible, but the effects of compressibility are negligible at speeds well below the
sound speed.

• In general, the viscosity µ may depend on the state variables, e.g. ρ, u, p or T , but we will take it to be
constant.

1.6 The Navier-Stokes equations

• For an incompressible Newtonian viscous fluid (19) and (21) give
 
∂ui ∂uj
σij = −pδij + µ + .
∂xj ∂xi

• We calculate
  
∂σij ∂ ∂ui ∂uj
= −pδij + µ +
∂xj ∂xj ∂xj ∂xi
2 2
∂ uj
∂p ∂ ui
= − +µ +µ
∂xi ∂xj ∂xj ∂xj ∂xi
 
∂p 2 ∂ ∂uj
= − + µ∇ ui + µ
∂xi ∂xi ∂xj
∂p ∂
= − + µ∇2 ui + µ (∇ · u)
∂xi ∂xj
∂p
= − + µ∇2 ui
∂xi

since for an incompressible fluid,

∇ · u = 0. (22)

• Hence, we have derived from Cauchy’s momentum equation (13) the incompressible Navier-Stokes equation

Du
ρ = −∇p + µ∇2 u + ρF. (23)
Dt

Remarks

(i) Note that the incompressible Navier-Stokes equations (22)–(23) consist of four scalar equations, which is
the same as the number of unknowns (u1 , u2 , u3 , p).

19
 (ii) In Cartesian coordinates Ox1 x2 x3 ,
∂ ∂2
u · ∇ = uj , ∇2 = ,
∂xj ∂xj ∂xj
so (22)–(23) are given in component form by
 
∂uj ∂ui ∂ui ∂p ∂ 2 ui
= 0, ρ + uj =− +µ + ρFi (i = 1, 2, 3).
∂xj ∂t ∂xj ∂xi ∂xj ∂xj
In other coordinate systems (e.g. cylindrical and spherical polar coordinates) the basis vectors themselves
depend on the coordinates. To calculate the components of the momentum equation in the direction of the
basis vectors, use the identities
 
1 2
(u · ∇) u = (∇ ∧ u) ∧ u + ∇ |u| ,
2
∇2 u = ∇(∇ · u) − ∇ ∧ (∇ ∧ u)

to write the momentum equation in the form

 
∂u p 1 2
+ω∧u+∇ + |u| = −ν∇ ∧ ω + F,
∂t ρ 2
where ω = ∇∧u is the vorticity and ν = µ/ρ is the kinematic viscosity; then use the usual expressions for ∇
and ∇∧ in the relevant coordinate system (e.g. see Acheson Appendices A.6 and A.7 for the Navier-Stokes
equations in cylindrical and spherical polar coordinates).
(iii) The force exerted by the fluid on a solid boundary S with unit normal n pointing into the fluid is given by
ZZ
t(n) dS.
S

For an incompressible Newtonian fluid the stress vector t(n) may be written in the form

t(n) = −pn + µ [2(n · ∇)u + n ∧ (∇ ∧ u)] ,

which quantifies the remarks made in the example in §2.2.2. Care must be taken to evaluate the term
(n · ∇)u for coordinate systems that are not Cartesian. In this course we will work almost exclusively in
Cartesian coordinates.
(iv) The Navier-Stokes equations (i.e. (22)–(23) with µ > 0) are of higher order than the Euler equations
(i.e. (22)–(23) with µ = 0) by virtue of the “diffusive” viscous term µ∇2 u, so it is necessary to impose
more boundary conditions for a viscous fluid than for an inviscid fluid.

1.7 Boundary conditions

1.7.1 Boundary conditions at a rigid impermeable boundary
• Suppose the fluid is in contact with a rigid impermeable surface S that has unit normal n pointing out of
the fluid and velocity U.
• Since fluid cannot flow though the impermeable surface, we prescribe the no-flux condition that

u·n=U·n on S,

so that the normal velocity components of the fluid and boundary are equal.
• For a viscous fluid, we also impose the no-slip condition

u − (u · n)n = U − (U · n)n on S,

so that the tangential velocity components of the fluid and boundary are equal.

20
 • Hence, the combined no-flux and no-slip boundary conditions are given by

u = U on S,

so that the velocities of the fluid and boundary are equal.

• Note that we have prescribed a total of three scalar boundary conditions.

Fluid velocity u

S
n
Boundary velocity U

1.7.2 Boundary conditions at a free surface

• Suppose the fluid has a free boundary Γ that has unit normal n pointing out of the fluid and outward
normal velocity V , the free boundary separating the fluid from a vacuum and being unknown a priori.

• Assuming there is no evaporation, we prescribe the no-flux condition that

u·n=V on Γ,

so that the normal velocity components of the fluid and free boundary are equal.

• Instead of the no-slip condition, we prescribe in the absence of surface tension the no-stress condition that

t(n) = 0 on Γ,

since the vacuum exerts no surface traction on the fluid.

• Note that we have prescribed a total of four scalar boundary conditions, one more than for a rigid imper-
meable boundary because we need to prescribe an additional equation to determine the location of the free
boundary.

n Vacuum

Γ
Fluid velocity u

1.8 Vorticity
• Vorticity ω = ∇ ∧ u is a measure of the local rotation of fluid elements.

• Assuming the body force is conservative (so that ∇ ∧ F = 0) and taking curl of the momentum equation
(23), we obtain the vorticity transport equation

∂ω
+ (u · ∇)ω − (ω · ∇)u = ν∇2 ω,
∂t
where the kinematic viscosity ν = µ/ρ.

• The effect of viscosity is to diffuse vorticity.

21
 • In two-dimensions with velocity
u = u(x, y, t)i + v(x, y, t)j,
the vorticity
ω = ∇ ∧ u = ω(x, y, t)k,
where the z-component of vorticity is given by
∂v ∂u
ω= − .
∂x ∂y

Hence, (ω · ∇)u = ω∂u/∂z = 0, and we obtain the two-dimensional vorticity transport equation

Dω ∂w ∂ω ∂ω
= +u +v = ν∇2 ω. (24)
Dt ∂t ∂x ∂y

• For an inviscid fluid (ν = 0), the two-dimensional vorticity transport equation becomes
Dω
= 0.
Dt
Thus, if ω = 0 at time t = 0, then ω = 0 for all t > 0 (Cauchy-Lagrange Theorem from Part A “Fluid
Dynamics and Waves”).

• Since ∇2 ω ≡ 0 if ω ≡ 0, might expect that adding diffusion (ν > 0) doesn’t change the argument. This is
incorrect because vorticity is generated at boundaries.

• To see this, use the fact the flow is incompressible to write (24) in the conservative form

∂ω
+ ∇ · Q = 0,
∂t
where the “vorticity flux” Q = ωu − ν∇ω; the two terms on the RHS of this expression correspond to
convective and diffusive transport of vorticity.

• Consider a stationary rigid boundary S with unit normal n pointing out of the fluid. In general the no-slip
condition (u = 0 on S) does not imply that Q · n = 0 on S. In particular, the boundary acts as an effective
source (sink) of vorticity if Q · n < 0 (Q · n > 0).

1.9 Conservation of energy

• We now consider the transport of energy by a conducting viscous fluid in the absence of external energy
sources (e.g. radiation, chemical reactions); cf. the heat conduction example in §1.5.2.

• Consider a material volume V (t) whose boundary ∂V (t) has outward unit normal n.

• The total internal energy in V (t) due to heat and kinetic energy is
ZZZ
1
E(t) = ρcv T + ρ|u|2 dV,
V (t) 2

where cv is the specific heat and T (x, t) the absolute temperature.

• Conservation of energy states that the time rate of change of the total internal energy increases due to

(i) conduction of heat into V (t) through ∂V (t), with net rate
ZZ
q · (−n)dS,
∂V (t)

where the heat flux vector q = −k∇T according to Fourier’s law, k being the thermal conductivity;

22
 (ii) work done by surface stresses t(n) on ∂V (t), with net rate
ZZ
t(n) · u dS;
∂V (t)

(iii) work done by body forces in V (t), with net rate

ZZZ
ρF · u dV.
V (t)

• Note that Fourier’s law may be written in the form

∂T
q = −k∇T = −kej
∂xj

and that Cauchy’s stress theorem (12) implies

t(n) · u = (ei σij nj ) · (ek uk ) = δik uk σij nj = ui σij nj .

• Hence, conservation of energy for the material volume V (t) may be written in the form
ZZ   ZZZ
dE ∂T
= k + ui σij nj dS + ρF · u dV.
dt ∂V (t) ∂xj V (t)

• Using the corollary to Reynolds’ transport theorem (8) and the divergence theorem implies that
ZZZ    
D 1 2 ∂ ∂T
ρ cv T + |u| − k + ui σij − ρF · u dV = 0.
V (t) Dt 2 ∂xj ∂xj

• Since V (t) is arbitrary, the integrand must be zero (if it is continuous), i.e.
   
D 1 2 ∂ ∂T
ρ cv T + |u| = k + ui σij + ρF · u.
Dt 2 ∂xj ∂xj

• Assuming that cv and k are constants, and using Cauchy’s momentum equation gives

DT ∂ui
ρcv = k∇2 T + Φ, Φ = σij .
Dt ∂xj

• Finally, substituting the constitutive law for an incompressible Newtonian fluid implies that the viscous
dissipation is given by
 
1 ∂ui ∂uj 2
Φ= µ + .
2 ∂xj ∂xi

• Hence, fluid deformation (⇒ Φ > 0) always increases the temperature.

1.10 Unidirectional flows

• There are hardly any explicit solutions of the Navier-Stokes equations. Almost all of them are for unidi-
rectional flows in which there is one flow direction only.

• Choose x-axis in direction of flow, i.e. set u = u(x, y, z, t)i.

23
 • Since
∂
u·∇=u ,
∂x
(22)–(23) become

∂u
= 0,
 ∂x  2 
∂u ∂u ∂p ∂ u ∂2u ∂2u
ρ +u = − +µ + 2 + 2 ,
∂t ∂x ∂x ∂x2 ∂y ∂z
∂p
0 = − ,
∂y
∂p
0 = − .
∂z

• Hence, u is independent of x; p is independent of y and z. It follows that u = u(y, z, t) and p = (x, t) satisfy
 2 
∂u ∂ u ∂2u ∂p
ρ −µ 2
+ 2 = − .
∂t ∂y ∂z ∂x

• LHS independent of x; RHS independent of y and z; hence LHS = RHS independent of x, y, z.

• Hence, the pressure gradient is a function of time t only, say

∂p
= G(t),
∂x
which must be prescribed to find u(y, z, t).

• The x-component of velocity u(y, z, t) satisfies the two-dimensional diffusion equation

 2 
∂u ∂ u ∂2u G(t)
=ν 2
+ 2 − , (25)
∂t ∂y ∂z ρ

the kinematic viscosity ν = µ/ρ being the diffusion coefficient and G(t) being the applied pressure gradient.

Remarks

(i) In unidirectional flows, all nonlinear terms in the Navier-Stokes equations vanish: the convective term

(u · ∇)u = 0.

(ii) The remaining equation (25) is linear and may be solved analytically using standard techniques in several
physically relevant geometries (which must be invariant to translations in the x-direction, as the flow is in
the x-direction).

(iii) In practical applications, the applied pressure gradient G(t) is zero, constant or oscillatory.

(iv) Further simplifications:

(a) in one-dimensional steady unidirectional flow with u = u(y), (25) reduces to the ordinary differential
equation
d2 u G
2
= ;
dy µ
(b) in one-dimensional unsteady unidirectional flow with u = u(y, t), (25) reduces to the one-dimensional
diffusion equation
∂u ∂ 2 u G(t)
=ν 2 − ;
∂t ∂y ρ

24
 (c) in two-dimensional steady unidirectional flow with u = u(y, z), (25) reduces to Poisson’s equation
∂2u ∂2u G
2
+ 2 = .
∂y ∂z µ
(v) The partial differential equations in (b) and (c) are amenable to standard methods from Moderations
or B568a, e.g. separation of variables and Fourier series methods, simularity reduction to an ordinary
differential equation, integral transforms.

1.10.1 Example: Poiseuille/Couette flow in a channel

• Consider flow in a channel 0 < y < h:

U
y=h

u = u(y)i

y=0

• Suppose lower plate at rest, upper plate moves to right with speed U and constant applied pressure gradient
∂p
= G.
∂x
• Assuming one-dimensional steady unidirectional flow with velocity u = u(y)i, the Navier-Stokes equations
reduce to the ordinary differential equation
d2 u
µ = G,
dy 2
which represents a balance of viscous shear forces and the applied pressure gradient.
• The no-flux boundary conditions on the plates are satisfied automatically, while the no-slip boundary
conditions imply that u(0) = 0, u(h) = U .
• Hence,
G Uy
u(y) = − y(h − y) + .
2µ h
• Special cases:
(i) Poiseuille flow (U = 0) driven by pressure gradient G 6= 0 has a quadratic velocity profile:
y=h

G
u(y) = − y(h − y)
2µ
y=0
2
Gh
u=0 u=
8µ

(ii) Couette flow (G = 0) driven by moving plate has a linear velocity profile:
Ui
y=h

Uy
u= i
h
y=0

25
1.10.2 Example: Shear stress in a Couette flow
Uy
• Consider the Couette flow u = u(y)i, with u(y) = :
h

U
y=h
σ12 (H)
y=H

y=0

• Flow occurs in layers y = H (constant).

• Fluid above y = H exerts a shear stress on fluid below y = H (and vice versa).
• The shear stress is given by
du µU
σ12 (H) = µ (H) = ,
dy h
where the subscript 1 indicates the x-component of stress and the subscript 2 that the normal to y = H is
in the y-direction.
• This shear stress arises because fluid above y = H is moving at a different speed than fluid below.
• Note that σ12 (H) > 0 for U > 0, as the fluid above y = H is moving faster than fluid below, i.e. viscosity
causes fluid above y = H to “drag along” fluid below (cf. inviscid fluid µ = 0).
• The shear stress exerted by the fluid on the lower plate is given by
µU
σ12 (0) = ;
h
by Newton’s third law, the shear stress exerted by the fluid on the upper plate is given by
µU
−σ12 (H) = − .
h
• The force per unit area in the x-direction required to sustain the motion of the upper plate is
µU
σ12 (H) = >0 for U > 0.
h
• If we can measure σ12 (H), U and h in our Couette flow rheometer, then we can calculate the viscosity µ
of an incompressible Newtonian fluid.

1.11 Dimensionless Navier-Stokes equations

• As with all mathematical modelling, we will only ever start to understand the mathematical implications
if we write the model in dimensionless variables.
• Consider the flow of an incompressible Newtonian fluid with far-field velocity U i past a stationary obstacle
of typical size L and with boundary ∂D.

u = 0 on ∂D
u → U i as |x| → ∞

26
 • In the absence of body forces, the flow is governed by the incompressible Navier-Stokes equations (22)–(23)
with F = 0, i.e.  
∂u
ρ + (u · ∇)u = −∇p + µ∇2 u, ∇ · u = 0,
∂t
with boundary conditions in the diagram above.

• Here and hereafter we will denote by [x] the typical dimensional size of the quantity x.

• The typical dimensional sizes of the dependent and independent variables are given by

length scale: [x] = L;

velocity scale: [u] = U ;

L
time scale: [t] = ;
U
pressure scale: [p] = to be determined.

• Hence, we nondimensionalize by scaling

L
x = Lx̂, u = U û, t= t̂, p = patm + [p]p̂.
U
where patm is the atmospheric pressure.

• Since xi = Lx̂i ,
∂ 1 ∂ 1 ˆ
∇ = ei = ei = ∇.
∂xi L ∂ x̂i L
• The incompressibility condition ∇ · u = 0 becomes
1 ˆ ˆ · û = 0.
∇ · (U û) = 0 ⇒ ∇
L

• Similarly, the momentum equation becomes

 
ρU 2 ∂ û ˆ [p] µU ˆ 2
+ (û · ∇)û = − ∇p̂ + 2 ∇ û.
L ∂t L L

• The ratio of the inertia term on the LHS to the viscous term on the RHS is given by

[inertia term] ρU 2 /L ρLU LU

= 2
= = = Re,
[viscous term] µU/L µ ν
which is the (dimensionless) Reynolds number.

• Note that two flows are dynamically similar if they satisfy the same dimensionless problem (i.e. same
geometry, governing equations, boundary conditions and dimensionless parameters).

• We will study both high and low Reynolds number flows.

High-Reynolds number flows Re  1

• Choose the inviscid pressure scale [p] = ρU 2 to obtain

∂ û ˆ 1 ˆ2 ˆ · û = 0.
+ (û · ∇)û = −∇p̂ + ∇ û ∇
∂ t̂ Re

• In this regime we hope to ignore the small viscous term and solve the inviscid Euler equations except in
thin layers on boundaries where viscosity is required to satisfy the no-slip boundary condition.

27
Low-Reynolds number flows Re  1

• Choose the viscous pressure scale [p] = µU/L to obtain

 
∂ û ˆ ˆ 2 û ˆ · û = 0.
Re + (û · ∇)û = −∇p̂ + ∇ ∇
∂ t̂

• In this regime we hope to ignore the inertia term and solve the resulting slow-flow equations:

ˆ 2 û = ∇p̂,
∇ ˆ ˆ · û = 0.
∇

When is viscosity important in practice?

• Typical values of L, U , ν and hence Re = LU/ν for a car travelling at 30 mph through air, a fish swimming
in water and for a marble falling through treacle are shown in the following table.

Object L U ν Re
Car 1m 10 m s−1 10−5 m2 s−1 106
Fish 0.1 m 0.1 m s−1 10−6 m2 s−1 104
Marble 1 cm 1 cm s−1 103 cm2 s−1 10−3

Remarks

• The Reynolds number is large for many everyday flows.

• Warning: Solution may generate its own length scale (e.g. tornado).

• If the Reynolds number is of order unity, then the Navier-Stokes equations must be solved numerically.
However, modern computers can’t get much past Re = 104 in realistic geometries.

• We will consider both large and small Reynolds number flows using asymptotics.

28
2 High Reynolds number flows
2.1 Thermal boundary layer on a semi-infinite flat plate
2.1.1 Dimensional problem
• Paradigm for viscous boundary layer on a semi-infinite flat plate.

• The two-dimensional steady heat convection-conduction problem consists of

– inviscid fluid, velocity U i, temperature T∞ upstream;

– plate at y = 0, x > 0, held at temperature Tp .

T = Tp on y = 0, x > 0
O x

u = Ui T → T∞ as x2 + y 2 → ∞

• Energy equation for temperature T , with µ = 0:

 
∂T
ρcv + (u · ∇)T = k∇2 T.
∂t

• We seek a steady solution T = T (x, y), with u = U i, so that

 2 
∂T ∂ T ∂2T
U =κ + ,
∂x ∂x2 ∂y 2

where κ = k/ρcv is the thermal diffusivity (units m2 s−1 ).

2.1.2 Dimensionless problem

• Choose arbitrary length scale L and set

x = Lx̂, y = Lŷ, T = T∞ + (Tp − T∞ )T̂ .

• The energy equation becomes

∂ T̂ ∂ 2 T̂ ∂ 2 T̂
Pe = + ,
∂ x̂ ∂ x̂ 2 ∂ ŷ 2
where the Peclet number P e = LU/κ; cf. Re = LU/ν.

• Boundary condition on plate: T̂ = 1 on ŷ = 0 < x̂.

• Boundary condition at infinity: T̂ → 0 as x̂ 2 + ŷ 2 → ∞.

• The time scale for heat energy to convect a distance L is L/U , while the time scale for heat energy to
diffuse a distance L is L2 /κ. Hence, the Peclet number

L2 /κ diffusion timescale
Pe = = .
L/U convection timescale

• We seek a solution for P e  1, i.e.  = 1/P e  1.

29
 • Dropping hats on the dimensionless variables, the dimensionless problem is given by
 
∂T ∂2T ∂2T
= + , (26)
∂x ∂x2 ∂y 2

where  = 1/P e  1, with boundary conditions

T =1 on y = 0, x > 0 (27)

and
T →0 as x2 + y 2 → ∞. (28)

2.1.3 Exact solution

• The exact solution to (26)–(28) is given by
Z ∞
2 2
T (x, y) = erfc(η) ≡ √ e−s ds,
π η

where erfc is the complementary error function and

" #1/2
(x2 + y 2 )1/2 − x
η(x, y) = .
2

erfc(η)
2
1 e−η
erfc(η) ∼ √
πη
as η → ∞

O 1 η

• Key observation: y 2 = (4η 2 x) + 2 (4η 4 ).

Deductions from T = erfc(η), y 2 = (4η 2 x) + 2 (4η 4 )

(i) Isotherms η = constant are parabolae.

(ii) For T not close to zero we need η = O(1) as  → 0. As  → 0 with x = O(1),

 
2 2 |y| |y|
η = O(1) ⇒ y ∼ (4η x) ⇒ η∼√
⇒ T ∼ erfc √ . (29)
4x 4x
√
Hence, there is a thermal boundary layer on the plate in which |y| = O( x) as  → 0 with x = O(1), as
illustrated.

Thermal boundary
!√ "
layer: |y| = O !x

Isotherms: η = const.

30
2.1.4 Boundary layer analysis
• Instead of solving exactly and then expanding, let us expand first and then solve.
• In outer region away from plate, expand
∂ T̄0
T ∼ T̄0 + T̄1 + · · · ⇒ = 0 ⇒ T̄0 = 0,
∂x
by the upstream boundary condition (28); in fact, T = O(n ) as  → 0 for all integer n, as we know from
the exact solution that T is exponentially small as  → 0 with |y| = O(1).
• To determine the thickness of the thermal boundary layer δ = δ() on the plate as  → 0, we scale y = δY
so that (26) becomes
∂T ∂2T  ∂2T
= 2 + 2
∂x ∂x δ ∂Y 2
• Since the LHS is of O(1), while the RHS is of O(/δ 2 ), it is necessary that /δ 2 = O(1) for a nontrivial
balance involving both convection and diffusion of heat energy.
• Hence, we set (without loss of generality) δ = 1/2 , so that y = 1/2 Y and
∂T ∂2T ∂2T
= 2 + .
∂x ∂x ∂Y 2
• Now expand T ∼ T0 (x, Y ) + T1 (x, Y ) + · · · to obtain the leading-order thermal boundary layer equation
∂T0 ∂ 2 T0
= . (30)
∂x ∂Y 2
• This is the heat equation with x playing the role of time.
• The boundary condition on the plate (27) becomes
T0 = 1 on Y = 0 < x. (31)

• To ensure that the boundary layer and outer expansions match (i.e. that they coincide in some intermediate
overlap region), we impose the matching condition
T0 (x, Y ) → 0 as |Y | → ∞. (32)

• The leading-order thermal-boundary-layer problem (30)–(32) has a similarity solution, viz.

 
|Y |
T0 (x, Y ) = erfc √ ,
4x
which by (29) is the leading-order term in the expansion of the exact solution as  → 0 because Y = y/1/2 .

2.1.5 Conclusions
• For  = 0 (no diffusion) the upstream condition demands that T ≡ 0, which doesn’t satisfy the boundary
condition on the plate.
• For 0 <   1, this solution applies at leading order except in a thin boundary layer near the plate in
which thermal diffusion (via the ∇2 T term) increases the temperature from its upstream value to that on
the plate.
• This is a singular perturbation problem as a uniformly valid approximation cannot be obtained by setting
the small parameter equal to zero (cf. examples of regular and singular perturbation problems in B568a).
• Singular behaviour arises because the small parameter  multiplies the highest derivative in (26).
• The highest derivative can be ignored except in thin regions where it is sufficiently large that it is no longer
annihilated by the premultiplying small parameter.
• Such regions usually occur near the boundary of the domain, and so they are called boundary layers.

31
2.2 Viscous boundary layer on a semi-infinite flat plate
2.2.1 Dimensional problem
• Consider the two-dimensional steady incompressible viscous flow of a uniform stream U i past a semi-infinite
flat plate at y = 0 < x.

• In the absence of body forces, the flow is governed by the incompressible Navier-Stokes equations (22)–(23)
with F = 0, which become
   2 
∂u ∂u ∂p ∂ u ∂2u
ρ u +v = − +µ + 2 , (33)
∂x ∂y ∂x ∂x2 ∂y
   2 
∂v ∂v ∂p ∂ v ∂2u
ρ u +v = − +µ + 2 , (34)
∂x ∂y ∂y ∂x2 ∂y
∂u ∂v
+ = 0, (35)
∂x ∂y

where u = u(x, y)i + v(x, y)j is the velocity, p(x, y) is the pressure, ρ is the constant density and µ is the
constant viscosity.

• The no-flux and no-slip boundary conditions on the plate are given by

u = 0, v = 0 on y = 0, x > 0. (36)

• The far-field boundary conditions are given by

u → U, v → 0 as x2 + y 2 → ∞. (37)

2.2.2 Dimensionless problem

• Choose arbitrary length scale L and set

x = Lx̂, u = U û, p = ρU 2 p̂

to obtain (dropping the hats ˆ):

 
∂u ∂u ∂p 1 ∂2u ∂2u
u +v = − + + 2 , (38)
∂x ∂y ∂x Re ∂x2 ∂y
 
∂v ∂v ∂p 1 ∂2v ∂2v
u +v = − + + , (39)
∂x ∂y ∂y Re ∂x2 ∂y 2
∂u ∂v
+ = 0, (40)
∂x ∂y

where the Reynolds number Re = LU/ν.

• The no-flux and no-slip boundary conditions on the plate (36) are unchanged, while the far-field conditions
(37) pertain with U = 1, as illustrated in the following diagram.

u = 0, v = 0 on y = 0, x > 0
O x

u → 1, v → 0 as x2 + y 2 → ∞

32
 • By (40), there exists a streamfunction ψ(x, y) such that
∂ψ ∂ψ
u= , v=− .
∂y ∂x

• Eliminating p between (38) and (39) by taking

∂ ∂
(38) − (39)
∂y ∂x
and substituting for u and v gives
∂ψ ∂ ∂ψ ∂ 1 2 2 1 4
∇2 ψ − ∇2 ψ = ∇ (∇ ψ) ≡ ∇ ψ,
∂y ∂x ∂x ∂y Re Re

where ∇4 is the biharmonic operator.

• Setting  = 1/Re, we have

∂(ψ, ∇2 ψ)
= ∇4 ψ. (41)
∂(y, x)
• The boundary conditions on plate become
∂ψ ∂ψ
= =0 on y = 0, x > 0,
∂x ∂y
so we may set
∂ψ
ψ= =0 on y = 0, x > 0 (42)
∂y
without loss of generality; that ψ is constant on the plate means that it is a streamline.

• The far-field boundary conditions imply that

ψ∼y as x2 + y 2 → ∞. (43)

2.2.3 High Reynolds number regime

• We seek a solution to (41)–(43) for Re  1, i.e.
1
=  1.
Re

• For  = 0 (inviscid flow) the solution is ψ ≡ y, which doesn’t satisfy the no-slip boundary condition on the
plate:
∂ψ
u= ≡ 1 6= 0 on y = 0, x > 0.
∂y
• For 0 <   1, expect this solution to apply at leading order except in a thin boundary layer near the plate
in which viscosity (via the ∇4 ψ term) reduces the u from its free stream value to zero:

Inviscid fluid (! = 0) Viscous fluid (0 < ! ! 1)

• Since there is no known exact solution, we use boundary layer theory.

33
2.2.4 Boundary layer analysis
• In the outer region away from the plate, we expand
ψ ∼ ψ0 + ψ1 + · · · ⇒ ψ0 = y (as expected).

• To determine the thickness of the boundary layer δ = δ() on the plate as  → 0, we scale y = δY .
• Since
∂ψ
u= ∼1
∂y
as  → 0 in the outer region, we also scale ψ = δΨ.
• The partial differential equation (41) becomes
 3   
δ ∂Ψ ∂ Ψ δ ∂3Ψ ∂Ψ δ ∂ 3 Ψ δ ∂3Ψ ∂ 4 Ψ 2δ ∂ 4 Ψ δ ∂ 4 Ψ
δ 3 + 2 − δ + = δ + +
δ ∂Y ∂x δ ∂Y 2 ∂x ∂x δ ∂x2 ∂Y δ 3 ∂Y 3 ∂x4 δ 2 ∂x2 ∂Y 2 δ 4 ∂Y 4

• Since the LHS is of O(1/δ), while the RHS is of O(/δ 3 ), it is necessary that /δ 3 = O(1/δ) for a nontrivial
balance involving both inertia and viscous terms.
• Hence, we set (without loss of generality) δ = 1/2 , so that y = 1/2 Y , ψ = 1/2 Ψ and
∂Ψ ∂ 3 Ψ ∂Ψ ∂ 3 Ψ ∂Ψ ∂ 3 Ψ ∂Ψ ∂ 3 Ψ 4
2∂ Ψ ∂4Ψ ∂4Ψ
 + −  − =  + 2 + .
∂Y ∂x3 ∂Y ∂Y 2 ∂x ∂x ∂x2 ∂Y ∂x ∂Y 3 ∂x4 ∂x2 ∂Y 2 ∂Y 4
• Expand Ψ ∼ Ψ0 + Ψ1 + · · · to obtain at leading order a version of Prandtl’s boundary layer equation:
∂Ψ0 ∂ 3 Ψ0 ∂Ψ0 ∂ 3 Ψ0 ∂ 4 Ψ0
− = . (44)
∂Y ∂Y 2 ∂x ∂x ∂Y 3 ∂Y 4
• The boundary conditions on plate (42) imply that
∂Ψ0
Ψ0 = =0 on Y = 0, x > 0. (45)
∂Y
• We consider the boundary layer on the top of the plate, with Y > 0.
• Since the outer solution ψ ∼ ψ0 = y as  → 0 and y = 1/2 Y , ψ = 1/2 Ψ, the matching condition is
Ψ0 ∼ Y as Y → ∞. (46)

Justification for matching condition (not examinable)

• To ensure that the inner (boundary layer) and outer expansions match (i.e. that they coincide in some
intermediate overlap region), introduce e.g. the intermediate variable
y
ȳ = α = (1/2−α) Y,

where 0 < α < 1/2, so that y → 0 and Y → ∞ as  → 0 with ȳ = O(1) fixed.
• In the outer region substitute y = α ȳ and expand as  → 0 with ȳ = O(1) fixed:
ψ(x, y) ∼ ψ0 (x, α ȳ) ∼ α ȳ.

• In the inner region substitute Y = (α−1/2) ȳ and expand as  → 0 with ȳ = O(1) fixed:
ψ(x, y) = 1/2 Ψ(x, Y ) ∼ 1/2 Ψ0 (x, (α−1/2) ȳ).

• Hence, the expansions agree in the overlap region in which ȳ = O(1) provided
1/2 Ψ0 (x, (α−1/2) ȳ) ∼ α ȳ as  → 0,
i.e.
Ψ0 (x, (α−1/2) ȳ) ∼ (α−1/2) ȳ as  → 0,
which can only be the case if
Ψ0 (x, Y ) ∼ Y as Y → ∞.

34
2.3 Alternative derivation of the boundary layer equations
2.3.1 Dimensionless problem

• Recall that the dimensionless two-dimensional steady Navier-Stokes equations are given by
 
∂u ∂u ∂p ∂2u ∂2u
u +v = − + + 2 , (47)
∂x ∂y ∂x ∂x2 ∂y
 
∂v ∂v ∂p ∂2v ∂2u
u +v = − + + , (48)
∂x ∂y ∂y ∂x2 ∂y 2
∂u ∂v
+ = 0, (49)
∂x ∂y

where u = u(x, y)i + v(x, y)j is the velocity, p(x, y) is the pressure and  = 1/Re  1.

• The no-flux and no-slip boundary conditions on the plate are given by

u = 0, v = 0 on y = 0, x > 0. (50)

• The far-field boundary conditions are given by

u → 1, v → 0 as x2 + y 2 → ∞. (51)

2.3.2 Boundary layer analysis

• If  = 0, then u = 1, v = 0; but this solution of the Euler equations doesn’t satisfy the no-slip boundary
condition on the plate.

• If 0 <   1, then u ∼ 1, v = o(1) away from the plate on both sides of which there is a thin viscous
boundary layer.

• To determine the thickness of the viscous boundary layer δ = δ() on the plate as  → 0, we scale y = δY .

• By (49),
∂u 1 ∂v
+ = 0,
∂x δ ∂y
so we scale v = δV for a nontrivial balance.

• By (47),
∂u ∂u ∂p ∂2u  ∂2u
u +V =− + 2 + 2 ,
∂x ∂Y ∂x ∂x δ ∂Y 2
so we set δ = 1/2 for a nontrivial balance involving both inertia and viscous terms.

• Hence, the boundary layer scalings are y = 1/2 Y , v = 1/2 V and (47)–(49) become

∂u ∂u ∂p ∂2u ∂2u
u +V = − + 2 + ,
∂x ∂Y ∂x ∂x ∂Y 2
 
∂V ∂V ∂p ∂2V ∂2u
 u +V = − + 2 2 +  2 ,
∂x ∂Y ∂Y ∂x ∂Y
∂u ∂V
+ = 0.
∂x ∂Y

35
 • Expanding
u ∼ u0 + u1 + · · · , V ∼ V0 + V1 + · · · , p ∼ p0 + p1 + · · · ,
we obtain at leading order Prandtl’s boundary layer equations
∂u0 ∂u0 ∂p0 ∂ 2 u0
u0 + V0 = − + , (52)
∂x ∂Y ∂x ∂Y 2
∂p0
0 = − , (53)
∂Y
∂u0 ∂V0
+ = 0. (54)
∂x ∂Y
• The no-flux and no-slip boundary conditions on the plate (50) are unchanged at leading order, i.e.
u0 = 0, v0 = 0 on y = 0, x > 0. (55)

• The far-field boundary conditions (51) are replaced by the matching condition that
u0 → 1 as |Y | → ∞, x > 0, (56)
which ensures that the leading-order solutions in the outer region (away from the plate) and in the inner
boundary-layer region (on the plate) are in agreement in an intermediate ‘overlap’ region between them
(matching procedure not examinable).
• Note that a matching condition is not required for V0 because there are no terms involving the second-order
derivatives of V0 in the boundary layer equations (52)–(54).

2.3.3 Implications of the matching condition

• By (53), the pressure does not vary across the boundary layer at leaching order, i.e. p0 = p0 (x).
• By (56),
∂u0 ∂ 2 u0
→ 0, →0 as |Y | → ∞,
∂Y ∂Y 2
so taking the limit |Y | → ∞ in (52) we deduce that
∂p0
→0 as |Y | → ∞;
∂x
since p0 is independent of Y ,
∂p0
≡ 0.
∂x

2.3.4 Streamfunction formulation

• By (54), there exists a streamfunction Ψ0 (x, Y ) such that
∂Ψ0 ∂Ψ0
, V0 = −
u0 = ,
∂Y ∂x
in terms of which (52), with ∂p0 /∂x = 0, becomes
∂Ψ0 ∂ 2 Ψ0 ∂Ψ0 ∂ 2 Ψ0 ∂ 3 Ψ0
− = . (57)
∂Y ∂x∂Y ∂x ∂Y 2 ∂Y 3
• The boundary conditions (55) become (setting Ψ0 = 0 on the plate, without loss of generality)
∂Ψ0
Ψ0 = 0, =0 on Y = 0 < x. (58)
∂Y
• The matching condition (56) becomes
Ψ0 ∼ Y as |Y | → ∞. (59)

36
2.4 Blasius’ similarity solution
• Note that for all α > 0 the boundary-layer problem (57)–(59) is invariant under the transformation

x 7→ α2 x, Y 7→ αY, Ψ0 7→ αΨ0 . (60)

• Hence, the solution (assuming it exists and is unique) can involve x, Y and Ψ0 in combinations invariant
under this transformation only, e.g.
  x
Ψ0 

 
 Y2
Y
a function of
Ψ0  


 Y .
x1/2 x1/2

• We choose to seek a solution of the form

Y
Ψ0 = x1/2 f (η), η= (61)
x1/2
and use the chain rule to show that (57)–(59) reduces to Blasius’ equation
1
f 000 + f f 00 = 0, (62)
2
with boundary conditions
f (0) = 0, f 0 (0) = 0, f 0 (∞) = 1, (63)
where prime 0 denotes differentiation with respect to η.

Remarks

(i) The existence of the invariant transformation (60) implies the ansatz (61) reduces the partial-differential-
equation problem (57)–(59) to the ordinary-differential-equation problem (62)–(63). The proof uses group
theory and is beyond scope of course (see e.g. Ockendon & Ockendon, Appendix B for details), though we
shall use such invariant transformations below to facilitate other similarity reductions.

(ii) The ansatz (61) is called a similarity solution because knowledge of Ψ0 at x = x0 > 0 is sufficient to
determine Ψ0 for all x > 0 by a suitable scaling of x, i.e. the solution looks geometrically the same at
different values of x.

2.4.1 Numerical solution of (62)–(63)

• Blasius’ equation is a third-order, nonlinear, autonomous ordinary differential equation.

• There are three boundary conditions at different boundaries (two at η = 0 and one at η = ∞), so (62)–(63)
is a two-point boundary value problem.

• Although there is no explicit solution, the boundary value problem may transformed into an initial value
problem and solved numerically, as follows.

• Consider the transformation

f (η) = γF (ξ), η = ξ/γ,
where γ > 0 is a constant.

• Since
f 0 = γ 2 F 0 , f 00 = γ 3 F 00 , f 000 = γ 4 F 000 ,
Blasius’ equation (62) becomes
1
F 000 + F F 00 = 0. (64)
2

37
• Hence, if we solve (64) subject to the initial conditions

F (0) = 0, F 0 (0) = 0, F 00 (0) = 1, (65)

then f (η) satisfies Blasius’ equation (62) and the boundary conditions

f (0) = 0, f 0 (0) = 0, f 0 (∞) = γ 2 F 0 (∞), (66)

i.e. (62)–(63) provided γ = F 0 (∞)−1/2 .

• The initial value problem (64)–(65) may be solved numerically in maple by formulating it as a system of
first-order differential equations. It is found that γ = C 1/3 , where f 00 (0) = γ 3 = C ≈ 0.332 is a constant
used below.
> # SOLVE NUMERICALLY BLASIUS INITIAL VALUE PROBLEM
ODE := diff(F(t),t) = U(t), diff(U(t),t) = V(t), diff(V(t),t) =
-1/2*F(t)*V(t):
ICS := F(0)=0, U(0)=0, V(0)=1:
SOL := dsolve([ODE,ICS],numeric):
# FIND CONSTANTS
gamma0 := (rhs(SOL(10)[3]))^(-1/2);
C := gamma0^3;
gamma0 := (rhs(SOL(100)[3]))^(-1/2);
C := gamma0^3;
# PLOT VELOCITY PROFILE
with(plots):
odeplot(SOL,[gamma0^2*U(t),t/gamma0],0..6,color=black,thickness=
2,view=[0..1,0..8],labels=["f'(t)", "t"]);

!0 := 0.6924754573
C := 0.3320573956
!0 := 0.6924754573
C := 0.3320573956

t 4

0
0 0.2 0.4 0.6 0.8 1
f'(t)

• This trick was spotted in 1942 by Weyl, who also proved that there exists a unique solution to (62)–(63);
35 years later Blasius’ equation (62) was reduced to a first-order ordinary differential equation.

38
2.4.2 Implications
• The dimensional shear stress on the plate is given by
 
µU 1 ∂u 1/2 ∂V
σ12 = + .
L 1/2 ∂Y ∂x Y =0

Since V = 0 on Y = 0,
 1/2
µU ∂u µU ∂ 2 Ψ0 νU 3
σ12 = 1/2 ∼ 1/2 =ρ f 00 (0),
 L ∂Y Y =0  L ∂Y 2 Y =0 Lx

where f 00 (0) = C ≈ 0.332 from above.

• That σ12 → ∞ as x → 0 reflects the fact that Prandtl’s theory is invalid at the leading edge of the plate. To
find the solution locally it is necessary to solve the full problem (as in the paradigm heated plate problem).

• Ignoring edge singularities, the drag on one side of a plate of length L is given by
Z 1
σ12 L dx ∼ 2f 00 (0)ρ(νU 3 L)1/2 ,
0

where f 00 (0) = C ≈ 0.332. This prediction compares well with experiment for 103 ≤ Re ≤ 105 .

2.5 Variable external flow

• Consider the steady two-dimensional flow of a stream of viscous fluid with far-field velocity U i past an
obstacle of typical dimension L, with Re = LU/ν  1.

u = 0 on boundary
u → U i as |x| → ∞

• We might expect a thin boundary layer around the body of thickness of O(L/Re1/2 ) in which viscous effects
are important. This is only partly true, as we shall see.

• On smooth segments of the boundary the boundary layer analysis is the same, but in curvilinear coordinates.

• We denote by Lx, Ly the distance along and normal to the surface, so that x, y are dimensionless:

Boundary layer
thickness O(L/Re1/2 )

Ly Lx

• The envisaged boundary-layer scaling is then y = Y /Re1/2 , where Y = O(1) as Re → ∞, so that locally
Y = 0 looks like a flat plate.

• It can be shown that there are no new terms in the boundary layer equations, so the only way the boundary
layer “knows” it is on a curved surface is through the matching condition.

39
 • Thus, in the boundary layer the dimensional streamfunction
LU
ψ∼ Ψ as Re → ∞,
Re1/2
where the dimensionless streamfunction Ψ(x, Y ) satisfies the boundary layer equation

∂Ψ ∂ 2 Ψ ∂Ψ ∂ 2 Ψ dp0 ∂ 3 Ψ
− = − + , (67)
∂Y ∂x∂Y ∂x ∂Y 2 dx ∂Y 3
with boundary conditions
∂Ψ
Ψ= =0 on Y = 0, (68)
∂Y
and the matching condition
∂Ψ
u0 = → Us (x) as Y → ∞, (69)
∂Y
where Us (x) > 0 is the dimensionless slip velocity predicted by the leading-order-outer inviscid theory.

• As Y → ∞, the matching condition (69) implies that

∂2Ψ dUs ∂2Ψ ∂3Ψ

→ , → 0, → 0.
∂x∂Y dx ∂Y 2 ∂Y 3

• Hence, taking the limit Y → ∞ in (67) gives

dp0 dUs
= −Us , (70)
dx dx
so that
1
p0 (x) + Us (x)2 = constant,
2
which is Bernoulli’s equation for the leading-order-outer inviscid flow on y = 0+ .

• There is a similarity reduction of the partial differential equation problem (67)–(70) to an ordinary differ-
ential equation problem only if
Us (x) ∝ (x − x0 )m or ec(x−x0 ) ,
where x0 , m and c are real constants.

2.5.1 The Falkner-Skan problem: Us (x) = xm

• Since the partial-differential-equation problem (67)–(69) with Us (x) = xm is invariant to the transformation

x 7→ αx, Y 7→ α(1−m)/2 Y, Ψ 7→ α(1+m)/2 Ψ,

for all α > 0, we seek a similarity solution of the form

Ψ Y
= f (η), η= .
x(1+m)/2 x(1−m)/2
to obtain the Falkner-Skan problem
1 + m 00
f 000 + f f + m(1 − (f 0 )2 ) = 0, (71)
2
with
f (0) = 0, f 0 (0) = 0, f 0 (∞) = 1. (72)

Remarks

(i) The Falkner-Skan equation (71) is a third-order nonlinear ordinary differential equation.

40
 (ii) The Falkner-Skan problem (71)–(72) is a two-point boundary value problem that reduces to Blasius’ prob-
lem for m = 0.
(iii) (71) is autonomous, so the order can be lowered by seeking a solution in which f 0 is the dependent variable
and f the independent one.
(iv) (71)–(72) can be transformed to an IVP for m = 0 only, so the numerics are harder than for Blasius’
problem, but amenable to “shooting” methods.

Numerical Solution of the Falkner-Skan problem (71)–(72)

• In the plots below note that in the boundary layer the x-component of velocity is given by
∂Ψ Y
u0 = = xm f 0 (η), η= .
∂Y x(1−m)/2
• There are three cases, as follows.
(i) For m > 0, there is a unique solution with a monotonic velocity profile:
Y
x(1−m)/2

u
0 Us (x)

(ii) For −0.0904 < m < 0, there are two solutions, one with a monotonic velocity profile as in case (i),
the other with flow reversal:
Y
x(1−m)/2

u
0 Us (x)

(iii) For m < −0.0904, there is no solution.

• We return to these results after we have described the physics required to make sense of them.

2.6 Breakdown of Prandtl’s theory

2.6.1 Favourable and adverse pressure gradients
• Assuming Us (x) > 0, so that the slip velocity is in the x-direction, the sign of Us0 = dUs /dx determines
whether the flow in the boundary layer is driven by
(i) an accelerating outer flow (Us0 > 0) due to a favourable pressure gradient (p00 = −Us Us0 < 0);
(ii) a decelerating outer flow (Us0 < 0) due to an adverse pressure gradient (p00 = −Us Us0 > 0).
• Example: For Us (x) = xm , we calculate
dp0
Us0 (x) = mxm−1 , = −mx2m−1 ,
dx
so the outer flow along the positive x-axis is
(i) accelerating and driven by a favourable pressure gradient for m > 0;
(ii) decelerating and driven by an adverse pressure gradient for m < 0.

41
2.6.2 Flow reversal
• Numerical simulations of the boundary layer equations show that an adverse pressure gradient (p00 > 0)
causes flow reversal near the boundary.

• Flow reversal first occurs at the point of zero skin friction (i.e. zero shear stress) on the boundary:

∂u0 xs ∂u0 ∂u0 x

(x < xs , 0) > 0 (xs , 0) = 0 (x > xs , 0) < 0
∂Y ∂Y ∂Y

• It may be shown that flow reversal is unstable, so

∂u0
>0
∂Y Y =0

is a necessary condition for the validity of Prandtl’s asymptotic solution (in which viscous effects are
confined to a thin-boundary layer on the body).

2.6.3 Separation
• If the skin friction vanishes at x = xs , as in the schematic above, then the boundary layer equations break
down near (xs , 0) due to the formation of a (Goldstein) singularity, with

∂ 2 u0
→ ∞ as (x, Y ) → (xs , 0).
∂x∂Y

• This means that the boundary layer solution cannot be extended into x > xs (because the boundary layer
equations are parabolic, with timelike variable x; see B5b).

• Near (xs , 0), the asymptotic structure is described by “triple deck theory” (see Acheson §8.6, Ockendon &
Ockendon §2.2.4), which predicts separation of the boundary layer from the body at the “separation point”
(xs , 0):

Upper deck: O(Re−3/8 )

Main deck: O(Re−1/2 )

xs
Lower deck: O(Re−5/8 )
O(Re−3/8 )

• Physical explanation of boundary layer separation (Prandtl 1904): While the free fluid on the edge of the
boundary layer has sufficient momentum to traverse an adverse pressure gradient (Us > 0, Us0 < 0, p00 > 0),
fluid deep in the boundary layer, having lost a part of its momentum (due to friction), cannot penetrate
far into a field of higher pressure, and instead turns away from it.

• Separation sheds vorticity into the outer region, rendering invalid the inviscid irrotational approximation,
and hence the whole Prandtl picture.

42
2.6.4 Example: Circular cylinder
• Consider the boundary layer on a circular cylinder of unit radius with far-field velocity i.

x = π/2 r cos θ i + r sin θ j

x j

x=0 O i x=π

• In polar coordinates (r, θ), the velocity potential

 
1
φ(r, θ) = r + cos θ
r

if there is no circulation.

• Denoting by x the distance along the upper boundary from the forward stagnation point at −i, the slip
velocity
∂φ
Us (x) = − (1, π − x) = 2 sin x,
∂θ
giving a pressure gradient
dp0 (x) dUs
= −Us = −2 sin 2x.
dx dx
• Hence, there is

(i) a favourable pressure gradient between the forward stagnation point at x = 0 and a maximum of the
slip velocity at x = π/2;
(ii) an adverse pressure gradient between x = π/2 and the rear stagnation point at x = π.

• A numerical simulation of Prandtl’s boundary layer equations shows that flow reversal, and hence sep-
aration, occurs at x = xs ≈ 1.815, just downstream of the onset of the adverse pressure gradient at
x = π/2 ≈ 1.571.

• This completely destroys the Prandtl picture.

• In practice separation occurs for Re ≈ 5 − 10, with two circulating vortex pairs in the wake:

Boundary layer
separation

• As the Reynolds number is increased, the flow in the wake becomes unstable at Re ≈ 102 and turbulent at
Re ≈ 105 .

• The flow up to separation is only slightly affected in the steady laminar regime (Re ≤ 102 ), the prediction
of the separation point xs ≈ 1.815 being quite good.

• Deep in the turbulent regime at Re ≥ 3.5 × 106 , the boundary layer itself becomes turbulent and the
separation point moves toward the rear stagnation point. The result is a sudden reduction in the size of
the wake and hence the drag on the cylinder - this is called the “drag crisis.”

43
2.6.5 Example: The theory of flight
• Note that Us (x) = xm corresponds to flow past a corner or a wedge in the outer region.
• In polar coordinates (r, θ), the velocity potential
πθ ∂φ
φ(r, θ) = rπ/α cos ⇒ Us (x) = (x, 0) ∝ xm ,
α ∂r
where m = π/α − 1:

m > 0 (α < π) m < 0 (α > π)

Us (x)
α x
α Us (x)
Us (0) = 0 x Us (0) = ∞

• This means we can consider the flow in the boundary layer near the sharp trailing edge of an aerofoil.
• If there were no circulation, the outer inviscid irrotational flow would “turn the corner” at the trailing edge:

• On the upper surface of the aerofoil, Prandlt’s boundary layer problem near the trailing edge would be the
same as the Falkner-Skan problem with m ≈ −1/2 (as m = π/α − 1 and α is close to 2π).
• Since m < 0, there would be a strong adverse pressure gradient as the flow turns the corner, which would
result in separation of the boundary layer from the aerofoil.
• Such a flow is physically unrealistic, an experimental observation consistent with there being no solution
to the Falkner-Skan problem for m < −0.0904.
• Thus, boundary layer theory provides theoretical justification for the Kutta-Joukowski hypothesis, which
says that the appropriate outer inviscid irrotational solution has a circulation just sufficient for the flow to
separate smoothly at the trailing edge (i.e. with finite velocity).
• In practice it is this solution that describes laminar flow past a blunt nosed aerofoil at a small angle of
attack, since the thin viscous wake left behind the aerofoil has a small effect on the outer inviscid flow:

Viscous effects important

44
3 Low Reynolds number flows
3.1 Slow flow past a circular cylinder
3.1.1 Dimensional problem
• Consider the two-dimensional steady incompressible viscous flow of a uniform stream U i past a rigid circular
cylinder of radius a, centre O.
• In the absence of body forces, the flow is governed by the incompressible Navier-Stokes equations (22)–(23)
with F = 0, which become
ρ(u · ∇)u = −∇p + µ∇2 u, ∇ · u = 0.
• The boundary conditions are in the following diagram.

y eθ
er
r
j θ
O i x

u = 0 on r = a
u → U i as r → ∞

3.1.2 Nondimensionalization
• In the slow-flow regime (see §1.11), we nondimensionalize by scaling
µU ∗
x = ax∗ , u = U u∗ , p = p .
a
to obtain (dropping the stars ∗ )
(u · ∇)u = −∇p + ∇2 u, ∇ · u = 0.
where the Reynolds number
ρU a
 = Re = .
µ

3.1.3 Streamfunction formulation

• In two dimensions the incompressibility condition implies the existence of a streamfunction ψ independent
of z such that
u = ∇ ∧ (ψk) ≡ curl(ψk) = ui + vj = ur er + uθ eθ .
where
∂ψ ∂ψ 1 ∂ψ ∂ψ
u= , v=− ; ur = , uθ = − .
∂y ∂x r ∂θ ∂r
• Hence the vorticity
ω = curl(u)

= curl2 (ψk)

= ∇(∇ · (ψk)) − ∇2 (ψk)

 
∂ψ
= ∇ − (∇2 ψ)k
∂z

= −(∇2 ψ)k.

45
 • We take (in Q6) the curl of the momentum equation to eliminate p, and obtain thereby the two-dimensional
steady vorticity transport equation
(u · ∇)ω = ∇2 ω
for the z-component of vorticity ω = −∇2 ψ.

• In Cartesian coordinates (x, y), with ψ = ψ(x, y),

∂ ∂ ∂ψ ∂ ∂ψ ∂
u·∇=u +v = −
∂x ∂y ∂y ∂x ∂x ∂y
so that
∂(ψ, ∇2 ψ)
 = ∇2 (∇2 ψ) ≡ ∇4 ψ,
∂(y, x)
where
∂2 ∂2
∇2 = + .
∂x2 ∂y 2

• In plane polar coordinates (r, θ), with ψ = ψ(r, θ),

∂ uθ ∂ 1 ∂ψ ∂ 1 ∂ψ ∂
u · ∇ = ur + = − ,
∂r r ∂θ r ∂θ ∂r r ∂r ∂θ
so that
 ∂(ψ, ∇2 ψ)
= ∇4 ψ,
r ∂(θ, r)
where
∂2 1∂ 1 ∂2
∇2 = + + .
∂r2 r ∂r r2 ∂θ2

3.1.4 Dimensionless problem

• We have derived the streamfunction formulation
∂(ψ, ∇2 ψ)  ∂(ψ, ∇2 ψ)
 = ∇4 ψ or = ∇4 ψ (73)
∂(y, x) r ∂(θ, r)

of the dimensionless two-dimensional steady incompressible Navier-Stokes equations, and will consider
below the slow-flow regime in which the Reynolds number  is small.

• The no-flux and no-slip boundary conditions on the cylinder become

∂ψ
ψ= =0 on r = 1, (74)
∂r
where we set ψ = 0 on r = 1 without loss of generality; note that r = 1 is a streamline.

• The far-field condition becomes

ψ ∼ y = r sin θ as r → ∞. (75)

3.1.5 Asymptotic solution for r = O(1) as  → 0

• We seek a regular perturbation solution to (73)–(75) by expanding

ψ ∼ ψ0 + ψ1 + · · · as  → 0.

• We obtain at leading-order the slow-flow approximation in which the inertia terms are neglected, i.e. the
biharmonic equation
∇4 ψ0 = 0.

46
• We would like to solve this fourth-order partial differential equation subject to the boundary conditions
(74)–(75), which become
∂ψ0
ψ0 = = 0 on r = 1,
∂r
ψ0 ∼ r sin θ as r → ∞.

• The far-field condition suggests we seek a separable solution of the form

ψ0 (r, θ) = f (r) sin θ,

which implies that

   
2 ∂2 1∂ 1 ∂2 ∂2 1∂ 1
∇ ψ0 = + + f (r) sin θ = + − 2 f (r) sin θ,
∂r2 r ∂r r2 ∂θ2 ∂r 2 r ∂r r

and hence that  2

4 ∂2 1∂ 1
∇ ψ0 = 2
+ − 2 f (r) sin θ.
∂r r ∂r r

• Thus, f (r) satisfies the fourth-order linear ordinary differential equation

 2
d2 1d 1
+ − f (r) = 0, (76)
dr2 r dr r2

with boundary conditions

f (1) = f 0 (1) = 0, f 0 (∞) = 1. (77)

• The ordinary differential equation (76) is homogeneous because it is invariant under the transformation
r → αr (α 6= 0), so we seek solutions of the form f (r) = rn to obtain
 
d2 1d 1
2
+ − 2 f (r) = [n(n − 1) + n − 1]rn−2 = [n2 − 1]rn−2
dr r dr r

so that  2
d2 1d 1
+ − f (r) = [n2 − 1][(n − 2)2 − 1]rn−4 .
dr2 r dr r2

• Hence n = −1, 1, 1, 3, so (76) has general solution

A
f (r) = + Br + Cr log r + Dr3 ,
r
where A, B, C and D are constants.

• Note that this general solution may also be derived by making s = log r the dependent variable.

• By the boundary conditions on the cylinder in (77),

f (1) = 0 ⇒ A + B + D = 0,

f 0 (1) = 0 ⇒ −A + B + C + 3D = 0,

giving
C = 4A + 2B, D = −(A + B).

• Applying the far-field condition in (77),

f 0 (∞) = 1 ⇒ B = 1, C = 0, D = 0.

47
 • But
C = 0, D = 0 ⇒ A = 0, B = 0,
so it is not possible to satisfy both the boundary conditions on the cylinder and the far-field condition.

• Hence, there is no separable solution of the form ψ(r, θ) = f (r) sin θ.

• Can prove rigorously nonexistence of a solution for slow flow past a cylinder of arbitrary cross-section.

• This is known as the Stokes paradox (1851), which wasn’t resolved until 1957!

3.1.6 Resolution of Stokes’ paradox

• The key observation that facilitates the resolution of the Stokes paradox is that (73)–(75) is a singular
perturbation problem because the slow-flow approximation is only valid for
1
r .


• Hence, instead of imposing the far-field condition

ψ0 ∼ r sin θ as r → ∞,

we need to match with a “boundary layer at infinity” (see B568a, Q5.5).

• To balance inertia and viscous terms we scale

x̂ ŷ r̂ ψ̂
x= , y= , r= , ψ= ,
   
so that at first it appears we recover the full two-dimensional formulation:
2
ˆ ψ̂)
∂(ψ̂, ∇ ˆ 4 ψ̂.
=∇
∂(ŷ, x̂)

• However, far away the cylinder has a small effect on the flow, so that

ψ̂ ∼ ŷ + δ ψ̂1 + · · · ,

where δ()  1 is to be determined.

• Obtain at O(δ) Oseen’s equation

2
ˆ ψ̂1 )
∂(ŷ, ∇ ˆ 4 ψ̂1
=∇
∂(ŷ, x̂)
or
∂ ˆ2 ˆ 4 ψ̂1
∇ ψ̂1 = ∇
∂ x̂
or
∂ û ˆ +∇ ˆ 2 û, ˆ · û = 0
= −∇p̂ ∇ (78)
∂ x̂
upon expanding u ∼ i + δ û + · · ·, p ∼ δ p̂, so that

∂ ψ̂1 ∂ ψ̂1
û = i− j.
∂ ŷ ∂ x̂

• There is no closed form solution to Oseen’s equations (78), but can use Fourier transforms to show that
the relevant solution (with zero flow at infinity) has the local expansion

ψ̂1 ∼ E r̂ log r̂ sin θ as r̂ → 0,

where E is an arbitrary constant.

48
 • It is then necessary to match the asymptotic expansions for ψ in the inner region (near the cylinder, with
r = O(1)) and in the outer region (far from the cylinder, with r = r̂/, r̂ = O(1)) by ensuring the constants
A, B, C, D and E are such that the expansions are in agreement in an intermediate ‘overlap’ region between
them.
Matching (details not examinable)
• Further reading: Perturbation Methods by E.J. Hinch; Perturbation Methods in Fluid Mechanics by M. Van
Dyke.
• To ensure that the inner and outer expansions match (i.e. that they coincide in some overlap region),
introduce the intermediate variable
r̂
r̄ = α r = 1−α (0 < α < 1),

so that r → ∞ and r̂ → 0 as  → 0 with r̄ = O(1) fixed.
• Recall that the leading-order inner solution
 
A 3
ψ0 = + Br + Cr log r + Dr sin θ
r
satisfies the boundary conditions on the cylinder r = 1 provided
C = 4A + 2B, D = −(A + B).

• In the inner region fix r̄ and expand as  → 0:

ψ ∼ ψ0 (−α r̄, θ) + · · ·
 r̄ 3 r̄  r̄ 
∼ D α sin θ + C α log α sin θ + · · ·
  
1
1
∼ 3α Dr̄ 3 sin θ + α ( − αC log ) r̄ sin θ + · · · .
 
• In the outer region fix r̄ and expand as  → 0:
1 
ψ ∼ r̂ sin θ + δ ψ̂1 + · · ·
 r̂=1−α r̄

1−α r̄ δ
∼ sin θ + E1−α r̄ log (1−α r̄) sin θ + · · ·
 
1
∼ (1 + (1 − α)δE log ) r̄ sin θ + · · · .
α
• The expansions agree in the overlap region in which r̄ = O(1) provided
D = 0, −αC log  = 1 + (1 − α)δE log .

• This is true for all 0 < α < 1 provided

1
C = δE =  1.
log(1/)

• Hence, the correct expansion for the streamfunction in the inner region in which r = O(1) is given by
 
1 r 1
ψ∼ r log r − + sin θ
log(1/) 2 2r
as  → 0.
• This asymptotic expansion for ψ proceeds in powers of 1/ log(1/), rather than in powers of  as originally
anticipated, so very slow convergence.

49
3.1.7 Exercise: Drag calculation
(i) Show that the dimensionless slow-flow approximation may be written in the form ∇p = −curl(ωk), where
ω = −∇2 ψ, so that in plane polar coordinates
∂p 1 ∂ω 1 ∂p ∂ω
=− , = .
∂r r ∂θ r ∂θ ∂r
(ii) Hence determine the leading-order terms in the expansions of the stress components
∂ur ∂  uθ  1 ∂ur
σrr = −p + 2 , σrθ = r + as  → 0, with r = O(1).
∂r ∂r r r ∂θ
(iii) Deduce that the dimensionless drag per unit length on the circular cylinder r = 1 is given by
Z 2π
4π
σrr (1, θ) cos θ − σrθ (1, θ) sin θ dθ ∼ .
0 ln (1/)

(iv) Deduce that the dimensional drag per unit length on a circular cylinder is approximately
4πµU
 ν ,
ln
Ua
which depends logarithmically on the radius a.

3.2 Slow flow past a sphere

3.2.1 Dimensional problem
• Consider the three-dimensional steady incompressible viscous flow of a uniform stream U k past a rigid
sphere of radius a, centre O.
• In the absence of body forces, the flow is governed by the incompressible Navier-Stokes equations (22)–(23)
with F = 0, which become
ρ(u · ∇)u = −∇p + µ∇2 u, ∇ · u = 0.

• The boundary conditions are in the following diagram.

z (x, y, z)
r
k θ
O
u = 0 on r = a

u → U k as r → ∞

3.2.2 Nondimensionalization and slow-flow approximation

• In the slow-flow regime (see §2.9), we nondimensionalize by scaling
µU ∗
x = ax∗ , u = U u∗ , p = p .
a
to obtain (dropping the stars ∗ )
Re(u · ∇)u = −∇p + ∇2 u, ∇ · u = 0.
with u = 0 on r = 1 and u → k as r → ∞, where the Reynolds number Re = ρU a/µ.

50
 • For Re  1, neglect the inertia term to obtain the slow-flow approximation:
∇p = ∇2 u, ∇ · u = 0.

• We begin by rewriting the momentum equation as

∇p = ∇(∇ · u) − ∇ ∧ (∇ ∧ u) = −curl2 u.

• Taking the curl we obtain the formulation

curl3 u = 0, ∇ · u = 0.

3.2.3 Streamfunction formulation

• Introduce spherical polar coordinate system (r, θ, φ), with origin at the centre of the sphere (Acheson, A.7):

z er
eφ
eθ
r
θ
y
φ
x

• Seeking a steady axisymmetric solution, with velocity u = ur (r, θ)er + uθ (r, θ)eθ , the incompressibility
condition becomes  
1 ∂ 2

∂
0=∇·u= 2 ur r sin θ + (uθ r sin θ) .
r sin θ ∂r ∂θ
• Satisfy this equation exactly by introducing the Stokes streamfunction ψ(r, θ) such that
1 ∂ψ 1 ∂ψ
ur = , uθ = −
r2 sin θ ∂θ r sin θ ∂r
giving
er reθ r sin θeφ  
1 ∂ ∂ ∂ ψ
u= 2 ∂r ∂θ ∂φ = curl eφ .
r sin θ r sin θ
0 0 ψ
• Hence the slow-flow equations reduce to
 
3 4 ψ
0 = curl u = curl eφ .
r sin θ

• Calculate:
er reθ r sin θeφ
1 ∂ ∂ ∂
curl(u) = 2 ∂r ∂θ ∂φ
r sin θ
ur ruθ 0
 
r sin θeφ ∂ ∂ur
= (ru θ ) −
r2 sin θ ∂r ∂θ
    
eφ ∂ 1 ∂ψ ∂ 1 ∂ψ
= − −
r ∂r sin θ ∂r ∂θ r2 sin θ ∂θ
 2 
eφ ∂ ψ cot θ ∂ψ 1 ∂2ψ
= − − 2 + 2 2 .
r sin θ ∂r2 r ∂θ r ∂θ

51
 • Hence  
2 ψ −D2 ψ
curl eφ = eφ
r sin θ r sin θ
where
∂2 cot θ ∂ 1 ∂2
D2 = − + 6= ∇2 .
∂r2 r2 ∂θ r2 ∂θ2

3.2.4 Dimensionless problem

• Hence    
4 ψ 2 −D2 ψ −D2 (−D2 ψ)
curl eφ = curl eφ = eφ ,
r sin θ r sin θ r sin θ
so we will need to solve
D4 ψ = 0. (79)

• The no-flux and no-slip boundary conditions on the sphere become

∂ψ
ψ= =0 on r = 1, (80)
∂r
where we set ψ = 0 on r = 1 without loss of generality.

• As r → ∞,
u → k = cos θer − sin θeθ ,
giving
1
∂ψ 1 ∂ψ
= ur → cos θ, − = uθ → − sin θ,
r2 sin θ
∂θ r sin θ ∂r
so the far-field condition becomes
1
ψ ∼ r2 sin2 θ as r → ∞. (81)
2

3.2.5 Stokes solution

• The far-field condition suggests we seek a separable solution of the form
   2
2 2 ∂2 2 2 4 ∂2 2
ψ = f (r) sin θ ⇒ D ψ= 2
− 2 f (r) sin θ ⇒ D ψ= 2
− 2 f (r) sin2 θ.
∂r r ∂r r

• Seeking a solution of the form f (r) = rn (as for a cylinder) leads to Stokes solution (1851)
 
1 2 3 1 −1
ψ= r − r+ r sin2 θ,
2 4 4

which satisfies the boundary conditions on the sphere and the far-field condition.

• Hence, there is a (separable) solution for slow flow past a sphere, cf. slow flow past a cylinder for which the
nonexistence of a solution (Stokes paradox) was resolved by matching with a “boundary layer at infinity”.

• For slow flow past a sphere nonuniformities arise in the asymptotic solution at O(Re) as Re → 0: this is
the Whitehead paradox and it is also resolved by matching with the relevant solution of Oseen’s equation
in a boundary layer at infinity.

• That the boundary layer at infinity is much weaker for a sphere than for a cylinder is to be expected, as a
cylinder is infinitely long.

• At leading order the streamfunction and hence the streamlines are symmetric fore and aft of the sphere
(and of the cylinder in the slow flow region), though a wake is observed in practice for moderate Reynolds
number.

52
3.2.6 Drag calculation
• The velocity components are given by
   
3 −1 1 −3 3 −1 1 −3
ur = 1 − r + r cos θ, uθ = −1 + r + r sin θ.
2 2 4 4

• Using ∇p = −curl2 u, we deduce that the pressure

3 cos θ
p = p∞ − ,
2r2
where p∞ is the constant pressure at ∞.
• The dimensionless force per unit area on the sphere is given by the dimensionless stress vector (made
dimensionless by scaling it with the pressure scale µU/a)

t(er )|r=1 = σrr er + σrθ eθ |r=1 .

• On the sphere the (nonzero) components of the dimensionless stress tensor (made dimensionless by scaling
them with the pressure scale µU/a) are given by
   
∂ur 3 ∂  uθ  1 ∂ur 3
σrr |r=1 = −p + 2 = −p∞ + cos θ, σrθ |r=1 = r + = − sin θ.
∂r r=1 2 ∂r r r ∂θ r=1 2

• Hence, the dimensionless force on the sphere is given by

ZZ Z 2πZ π
t(er ) dS = (σrr cos θ − σrθ sin θ) r2 sin θ r=1
dθdφ k
r=1 0 0

= (2π + 4π) k (normal stress + shear stress)

= 6π k.

• Recover famous Stokes formula for the dimensional drag:

 
µU
6π a2 = 6πµU a.
a

• Example: Provided Re = ρaU/µ  1, the terminal velocity U of a solid sphere (uniform density ρs ) falling
under gravity g though fluid (density ρ, viscosity µ) is given by
4 2 a2 (ρs − ρ)g
6πµU a = π(ρs − ρ)a3 g ⇒ U= .
3 9 µ

3.3 Lubrication Theory

3.3.1 Motivation
• How can a thin fluid layer support a large normal load?
• Example: sheet of paper sliding across a table.
• Lubrication theory is the analysis of flows in thin layers whose aspect ratio (depth / length) is small.
• Lubrication theory combines ideas from slow flow and boundary layer theory to simplify the Navier-Stokes
equations.
• Use the same methodology as for unidirectional flows: derive a general version of the simplified equation,
then apply it to several configurations that are of practical importance, namely flow in a slider bearing,
squeeze film and Hele-Shaw cell.

53
3.3.2 Governing equations
• Begin with the incompressible Navier-Stokes equations
Du
ρ = −∇p + µ∇2 u, ∇ · u = 0,
Dt
where ρ is the density, µ is the viscosity, u(x, y, z, t) = ui+vj+wk is the velocity, p(x, y, z, t) is the pressure
and the convective derivative
D ∂ ∂ ∂ ∂ ∂
= +u·∇= +u +v +w .
Dt ∂t ∂t ∂x ∂y ∂z
• Consider channel flow between rigid boundaries at z = 0 and z = h(x, y, t), where h(x, y, t) is the prescribed
channel thickness:
∂h z = h(x, y, t)
k
∂t z=0
z
y
Ui + V j
x

• Lower boundary has velocity U i + V j, so no-slip and no-flux boundary conditions are
u = U, v = V, w = 0 on z = 0.
∂h
• Assuming the upper boundary has velocity k (i.e. that it moves in the vertical direction only), the
∂t
no-slip and no-flux boundary conditions are
∂h
u = 0, v = 0, w = on z = h(x, y, t).
∂t
• We will prescribe appropriate boundary conditions at channel sides later on (after reduction of dimension-
ality).

3.3.3 Nondimensionalization
• Let U0 be a typical horizontal flow speed, i.e.
[u], [v], [U ], [V ] = U0 ,
where [ ] denotes here the typical size of the bracketed dimensional quantity.
• Let L and δL be typical horizontal and vertical length scales, respectively, i.e.
[x], [y] = L; [z], [h] = δL,
where the aspect ratio δ  1.
∂w ∂u ∂v
• Continuity equation: =− −
∂z ∂x ∂y
[w] [u] [w] U0
⇒ = ⇒ = ⇒ [w] = δU0 .
[z] [x] δL L
• Hence nondimensionalize by scaling
x = L x̂, y = L ŷ, z = δL ẑ,
u = U0 û, v = U0 v̂, w = δU0 ŵ,
h = δL ĥ, t = L
U0 t̂, p = patm + [p] p̂,

where δ is small (i.e. δ  1), patm is atmospheric pressure and the scale of pressure variations [p] is to be
determined.

54
3.3.4 The reduced Reynolds number
• The convective derivative becomes
 
D ∂ ∂ ∂ ∂ 1 ∂ ∂ ∂ ∂ 1 D
= +u +v +w = + û + v̂ + ŵ ≡ .
Dt ∂t ∂x ∂y ∂z L/U0 ∂ t̂ ∂ x̂ ∂ ŷ ∂ ẑ L/U0 Dt̂

• Hence the x-momentum equation becomes

 
ρU0 Dû [p] ∂ p̂ µU0 ∂ 2 û ∂ 2 û 1 ∂ 2 û
=− + 2 + + .
L/U0 Dt̂ L ∂ x̂ L ∂ x̂2 ∂ ŷ 2 δ 2 ∂ ẑ 2

• Dividing through by µU0 /δ 2 L2 gives

Dû δ 2 L[p] ∂ p̂ ∂ 2 û ∂ 2 û ∂ 2 û
δ 2 Re =− + δ2 2 + δ2 2 + 2 ,
Dt̂ µU0 ∂ x̂ ∂ x̂ ∂ ŷ ∂ ẑ
where
ρLU0 LU0 [x] [u]
Re = = =
µ ν ν
is the Reynolds number based on the length and velocity scale in the x-direction.

• For δ  1, the key parameter in the x-momentum equation is the reduced Reynolds number
(δL) (δU0 ) [z] [w]
δ 2 Re = = ,
ν ν
i.e. the Reynolds number based on the length and velocity scale in the vertical direction.

• In lubrication theory, we assume the reduced Reynolds number is small, i.e.

δ 2 Re  1,

so that inertial terms are much smaller than viscous terms; note that Re need not be small!

• Compare this with boundary layer theory in §3, where we set the reduced Reynolds number equal to unity
in the boundary layer, giving the boundary layer thickness δ = Re−1/2 for Re  1.

• Finally, we balance viscous and pressure terms (to avoid a triviality) by setting
δ 2 L[p] µU0
=1 ⇒ [p] = .
µU0 δ2L

• Hence the pressure scale in lubrication theory is a factor 1/δ 2 larger than the pressure scale in Stokes flow.

3.3.5 Dimensionless problem

• Substitute the scalings into the Navier-Stokes equations and simplify them to obtain
Dû ∂ p̂ ∂ 2 û ∂ 2 û ∂ 2 û
δ 2 Re = − + δ2 2 + δ2 2 + 2 ,
Dt̂ ∂ x̂ ∂ x̂ ∂ ŷ ∂ ẑ
Dv̂ ∂ p̂ ∂ 2 v̂ ∂ 2 v̂ ∂ 2 v̂
δ 2 Re = − + δ2 2 + δ2 2 + 2 ,
Dt̂ ∂ ŷ ∂ x̂ ∂ ŷ ∂ ẑ
Dŵ ∂ p̂ ∂ 2 ŵ ∂ 2 ŵ ∂ 2 ŵ
δ 4 Re = − + δ4 2 + δ4 2 + δ2 2 ,
Dt̂ ∂ ẑ ∂ x̂ ∂ ŷ ∂ ẑ
∂ û ∂v̂ ∂ ŵ
0 = + + ,
∂ x̂ ∂ ŷ ∂ ẑ
where the Reynolds number Re = ρLU0 /µ.

55
 • Scaling U = U0 Û and V = U0 V̂ , the no-flux and no-slip boundary conditions become

û = Û , v̂ = V̂ , ŵ = 0 on ẑ = 0,
∂ ĥ
û = 0, v̂ = 0, ŵ = on ẑ = ĥ(x̂, ŷ, t̂).
∂ t̂

3.3.6 Lubrication approximation

• Assume that both the aspect ratio δ and the reduced Reynolds number δ 2 Re are small, i.e.

δ, δ 2 Re  1,

to obtain at leading order the quasi-steady lubrication equations (dropping the hats ˆ):

∂p ∂ 2 u
0 = − + , (82)
∂x ∂z 2
∂p ∂ 2 v
0 = − + , (83)
∂y ∂z 2
∂p
0 = − , (84)
∂z
∂u ∂v ∂w
0 = + + , (85)
∂x ∂y ∂z

with boundary conditions

u = U, v = V, w = 0 on z = 0,
∂h (86)
u = 0, v = 0, w = on z = h(x, y, t).
∂t

3.3.7 Cross-layer averaged conservation of mass expression

• Integrate the incompressibility condition (85) from z = 0 to z = h:
Z h Z h  Z h 
∂u ∂v ∂w ∂ ∂h ∂ ∂h z=h
0= + + dz = u dz −u + v dz −v + [w]z=0 ,
0 ∂x ∂y ∂z ∂x 0 ∂x z=h ∂y 0 ∂y z=h

by Leibniz’s integral rule.

• Applying the kinematic boundary conditions in (86), we obtain the cross-layer-averaged conservation of
mass expression
∂h ∂   ∂  
+ hū + hv̄ = 0.
∂t ∂x ∂y
where the mean velocities are given by
Z h Z h
1 1
ū := u dz, v̄ := v dz,
h 0 h 0

3.3.8 Integrating the momentum equations

• To determine the mean velocities ū and v̄, need to integrate the momentum equations to obtain u and v.

• This is possible because the z-momentum equation (84) implies that at leading order the pressure is
independent of z.

56
 • Hence, for fixed (x, y), the x- and y-momentum equations (82)–(83) are ordinary differential equations in
z for u and v, giving after two integrations the expressions
1 ∂p 2 1 ∂p 2
u= z + A1 (x, y, t)z + B1 (x, y, t), v = z + A2 (x, y, t)z + B2 (x, y, t),
2 ∂x 2 ∂y
where Ai , Bi are arbitrary functions of (x, y, t).
• Determine Ai , Bi by applying the no-slip boundary conditions in (86) to obtain
1 ∂p  z 1 ∂p  z
u=− z(h − z) + U 1 − , v=− z(h − z) + V 1 − ,
2 ∂x h 2 ∂y h
i.e. Poiseuille / Couette flow in each direction.
• Integrate again to obtain mean velocities
Z Z
1 h h2 ∂p U 1 h
h2 ∂p V
ū = u dz = − + , v̄ = v dz = − + .
h 0 12 ∂x 2 h 0 12 ∂y 2
3.3.9 Reynolds lubrication equation
• In summary, we have obtained the cross-layer averaged conservation of mass expression
∂h ∂   ∂  
+ hū + hv̄ = 0, (87)
∂t ∂x ∂y
where the mean velocities
h2 ∂p U h2 ∂p V
ū = − + , v̄ = − + . (88)
12 ∂x 2 12 ∂y 2
• Substituting for ū, v̄, we obtain a version of Reynolds lubrication equation
 3   3 
∂ h ∂p ∂ h ∂p ∂h U ∂h V ∂h
+ = + + , (89)
∂x 12 ∂x ∂y 12 ∂y ∂t 2 ∂x 2 ∂y
which is an elliptic equation for the pressure p(x, y, t), since the thickness h(x, y, t) is prescribed (cf. the
thin-film equations in §4.3).
• Use Reynolds lubrication equation (89) to study flow in a slider bearing, squeeze film and Hele-Shaw cell.

3.3.10 Stress tensor

• Nondimensionalize the stress tensor σij and stress vector t by scaling them with the pressure scale [p], i.e.
set
µU0 µU0
σij = 2 σ̂ij , t = 2 t̂
δ L δ L
and then drop immediately the hats ˆ denoting dimensionless variables.
• In the lubrication regime, the dimensionless stress tensor is given by
∂u
σ11 = −p + 2δ 2 ∼ −p,
∂x
∂v
σ22 = −p + 2δ 2 ∼ −p,
∂y
∂w
σ33 = −p + 2δ 2 ∼ −p.
∂z
 
2 ∂u 2 ∂v 2 ∂u ∂v
σ12 = δ +δ ∼δ + ,
∂y ∂x ∂y ∂x
∂u ∂w ∂u
σ13 = δ + δ3 ∼δ ,
∂z ∂x ∂z
∂v ∂w ∂v
σ23 = δ + δ3 ∼δ .
∂z ∂y ∂z

57
 • Dimensionless stress on the upper wall is given by

t(n) = ei σij nj ,

where
δhx e1 + δhy e2 − e3
n = nj ej =
(1 + (δhx ) 2 + (δhy ) 2 )1/2
is the downward unit normal to the upper boundary.

• For δ  1, obtain at leading order

   
∂h ∂u ∂h ∂v
t(n) ∼ −δ p + e1 − δ p + e2 + p e3
∂x ∂z z=h ∂y ∂z z=h
   
∂(ph) 1 ∂p U ∂(ph) 1 ∂p V
= −δ − h − e1 − δ − h − e2 + p e3 . (90)
∂x 2 ∂x h ∂y 2 ∂y h

• Hence, vertical component of stress is O(1/δ) larger than horizontal components, which explains why a thin
viscous layer can support large normal loads, with relatively little horizontal resistance (cf. solid contact,
where friction is proportional to the normal load).

3.4 Examples
3.4.1 Steady flow in a one-dimensional slider bearing
• Set h = h(x) on 0 < x < 1, with U = 1, V = 0:

z = h(x)
U =1
x=0 x=1

• Reynolds lubrication equation (89) for the pressure p = p(x) becomes

 
d h3 dp 1 dh
= .
dx 12 dx 2 dx

• Typically ambient pressure at ends of bearing, so

p=0 at x = 0, 1.

• Integrate to obtain the pressure gradient

dp 6(h − h0 )
= ,
dx h3
where the constant h0 is to be determined.

• Note that p(x) has a stationary point at h = h0 .

• Integrate again and apply boundary conditions to obtain the solution

Z x Z 1
h(x) − h0 h(x) − h0
p(x) = 6 dx, dx = 0.
0 h(x)3 0 h(x)3

58
 • By (90), the dimensionless components of force on the slider in the x- and z-directions are given by
Z 1 Z 1
4h(x) − 3h0
Fx = δ dx, Fz = p(x)dx.
0 h(x)2 0

• Note Fx = O(δFz ), so high loads, low friction.

• Converging channel h0 < 0, gives p > 0 and slider supports load.

• Diverging channel h0 > 0, gives p < 0 and hence cavitation (i.e. gas bubble formation) if dimensional
pressure becomes too small.

• Cavitation may lead to contact and wear: see e.g. Ockendon & Ockendon, Viscous flow, §4.1.

3.4.2 Axisymmetric flow in a squeeze film

• Lower boundary stationary (U = V = 0) and upper boundary a flat circular unit disc at z = h(t), r < 1:

1
h(t)

• Seeking an axisymmetric solution p = p(r, t), Reynolds lubrication equation (89) becomes
 
1∂ rh3 ∂p dh 3ḣ 2

= ⇒ p(r, t) = 3
r + A(t) ln r + B(t) ,
r ∂r 12 ∂r dt h

where A(t), B(t) are arbitrary functions of t.

• Demand p bounded and impose ambient pressure p = 0 at r = 1 to obtain

 
3ḣ 2
p(r, t) = 3 r − 1 .
h

• By (90), the dimensionless upward force on the disc is given by

Z 2πZ 1
3π ḣ
Fz = p(r, t) rdrdθ = − .
0 0 2h3

• If ḣ > 0, this force is negative (suction makes disc adhere to plate) and vice versa.

• If a constant load N is applied to disc, then

Z t  
3π ḣ(τ ) 3π 1 1
Nt = − dτ = − ,
0 2h(τ )3 4 h(t)2 h(0)2

so it takes infinite time to squeeze fluid out of the gap!

• In practice, this prediction is rendered invalid by neglected effects, e.g. surface roughness, fluid compress-
ibility, plate curvature.

• Note that in the nondimensionalization, the typical plate velocity W0 (say) determines the appropriate
horizontal velocity scale U0 = W0 /δ (by the continuity equation) and time scale [t] = δL/W0 = L/U0 (by
the no-flux boundary condition on the plate).

59
3.4.3 Flow in a Hele-Shaw cell
• Flow between two parallel stationary plates (h constant, U = V = 0) driven by an external pressure
gradient:

• The mean velocities (88) are given by

h2 ∂p h2 ∂p
ū = − , v̄ = − ,
12 ∂x 12 ∂y

so Reynolds lubrication equation (89) reduces to Laplace’s equation for the pressure p(x, y, t):

∂2p ∂2p
+ = 0.
∂x2 ∂y 2

• Remarkably, since
∂ ū ∂v̄ ∂ ū ∂v̄
+ = 0, − = 0,
∂x ∂y ∂y ∂x
the mean flow (and in fact the flow at any fixed z) is the same as for Part A two-dimensional incompressible
inviscid irrotational flow with velocity potential −h2 p/12.

• However, there are three important distinctions.

(i) Although the streamlines are the same, the pressure is different (for an incompressible inviscid irrota-
tional fluid flow, the pressure is given by Bernoulli’s equation).
(ii) The circulation Γ around any closed curve C lying in a horizontal plane is zero, whether or not the
fluid domain is simply connected:
Z Z
1 ∂p ∂p 1
Γ= udx + vdy = − z(h − z) dx + dy = − z(h − z)[p]C = 0
C 2 C ∂x ∂y 2
because p is a single valued function of x and y. Prediction: can use Hele-Shaw cell to visualize
inviscid flows with zero circulation!
(iii) Near a smooth obstacle the lubrication approximation breaks down in a boundary layer of horizontal
width of the order of the gap thickness in which the cross-layer flow is governed at leading order by the
quasi-two-dimensional Stokes flow equations (as the local Reynolds number δ 2 Re  1), so boundary
layer structure completely different from Prandtl picture and no separation.

3.4.4 Injection at a point in a Hele-Shaw cell

• A blob of fluid Ω(0) containing the origin lies within a Hele-Shaw cell (h constant, U = V = 0) when at
time t = 0 a source of constant strength Q is introduced at the origin.

Ω(t)

60
• This is a free boundary problem, as fluid domain Ω(t) must be determined as part of solution.
• Convenient to introduce vector notation
∂ ∂
ū(x, y, t) = ūi + v̄j, ∇=i +j ,
∂x ∂y
so that Reynolds lubrication equation (89) and mean velocities (88) become
h2
∇ · (hū) = 0, ū = − ∇p.
12
• Since h is constant, pressure p(x, y, t) satisfies Laplace’s equation
∇2 p = 0 in Ω(t).

• Source strength Q at origin means

Q
hū ∼ er as r = (x2 + y 2 )1/2 → 0,
2πr
so that Z
hū · er ds → Q as  → 0.
r=

• On free boundary ∂Ω(t) need to impose two boundary conditions to determine its location.
• Suppose free boundary is given by f (x, y, t) = 0 and has outward unit normal
∇f
n= .
|∇f |

• The simplest boundary conditions are the zero pressure (neglecting surface tension) and kinematic condi-
tions
Df
p = 0, = 0 on f = 0,
Dt
where
D ∂
= + ū · ∇.
Dt ∂t
• The kinematic condition says the outward normal velocity of the fluid is equal to that of the boundary, i.e.
ū · ∇f
ū · n = = vn on ∂Ω(t),
|∇f |
the outward normal velocity of the boundary being given by
1 ∂f
vn = − .
|∇f | ∂t f =0

• Seek an axisymmetric solution in which

h2 ∂p
p = p(r, t), ū = − er , Ω(t) = {r < R(t)}.
12 ∂r
• Choose f (x, y, t) = x2 + y 2 − R(t)2 , so that
∇f 2xi + 2yj
n= = = er ,
|∇f | 2r
then outward normal velocity of boundary

1 ∂f −2RṘ
vn = − =− = Ṙ,
|∇f | ∂t f =0 2r
r=R
as expected.

61
 • Hence, lubrication problem reduces to
 
1∂ ∂p
r =0 for r < R(t),
r ∂r ∂r
h2 ∂p Q
− ∼ as r → 0,
12 ∂r 2πhr
h2 ∂p
p = 0, − = Ṙ on r = R(t).
12 ∂r

• Integrate and apply the boundary condition at the origin to obtain

h2 ∂p Q
− = .
12 ∂r 2πhr

• Integrate again and apply p(R(t), t) = 0 to obtain

 
6Q r
p(r, t) = − 3 log .
πh R(t)

• Finally, apply kinematic condition to obtain

Q
= Ṙ ⇒ πhR(t)2 = πhR(0)2 + Qt,
2πhr r=R

which represents global conservation of mass.

• Prediction: flow reversible upon changing sign of Q.

3.4.5 Squeeze film with a free boundary

• A blob of fluid lies between two parallel plates whose separation is h(t) (so h = h(t) and U = V = 0).

• Similar to injection problem, except flow driven by plate motion.

• Seeking an axisymmetric solution p = p(r, t) with a circular blob of radius R(t) results in the problem
 
1∂ h3 r ∂p
= ḣ for r < R(t),
r ∂r 12 ∂r
p bounded as r → 0,
h2 ∂p
p = 0, − = Ṙ on r = R(t).
12 ∂r

• The solution is
3ḣ(t) 2
p(r, t) = (r − R(t)2 ),
h(t)3

with πh(t)R(t)2 = πh(0)R(0)2 , representing global conservation of mass.

• Prediction: flow reversible upon reversing path of h(t).

62
3.5 The Saffman-Taylor instability and viscous fingering
• For the examples in §3.4.4 and §3.4.5, experiments show that the flow is not time reversible: when the
interface is advancing small perturbations to the interface shrink, so the circular base state is stable;
however, when the interface is retreating small perturbations grow, so the circular base state is unstable.

• To analyze this phenomena we use linear stability theory to study the local problem near the interface:
consider a Hele-Shaw cell consisting of a region of saturated flow in y < H(x, t) which is separated from a
dry region in y > H(x, t) by an interface at y = H(x, t).

• In the dimensionless variables of §3.4.3, with h2 /12 = 1 for brevity, the flow is governed by

∂ ū ∂v̄ ∂p ∂p
+ = 0, ū = − , v=− .
∂x ∂y ∂x ∂y

• Impose the zero-pressure and kinematic conditions

∂H ∂H
p = 0, v̄ = + ū on y = H(x, t).
∂t ∂x

• Take the base state to be the travelling wave solution moving upward with constant speed V , i.e.

ū = 0, v̄ = V, H = V t, p = p0 (y, t) = −V (y − V t),

so that the flow is driven by the prescribed pressure gradient

∂p
→ −V as y → −∞.
∂y

• The key step is to write the evolution problem relative to the base state by changing to a travelling wave
frame (x, η) moving with the interface by setting

y − V t = η,

ū = û(x, η, t),

v̄ − V = v̂(x, η, t),

p − p0 (y, t) = p̂(x, η, t),

H − V t = Ĥ(x, t),

where  is an arbitrary small parameter.

• Use the chain rule to derive the new problem

∂ û ∂v̂
+ =0 in η < Ĥ(x, t), (91)
∂x ∂η
with
∂ Ĥ ∂ Ĥ
p̂ = V Ĥ, v̂ = + û on η = Ĥ(x, t), (92)
∂t ∂x
and
p̂ = o(−η) as η → −∞, (93)
where
∂ p̂ ∂ p̂
û = − , v̂ = − . (94)
∂x ∂η
• As   1, can now proceed as in part A small amplitude water waves, as follows.

63
 (i) Linearize the boundary conditions (92) and impose them on η = 0, so that at leading order (91)–(92)
become
∂ û ∂v̂
+ = 0 in η < 0, (95)
∂x ∂η

with
∂ Ĥ
p̂ = V Ĥ, v̂ = on η = 0. (96)
∂t
Hence, by (93)–(96), p̂(x, η, t) satisfies Laplace’s equation

∂ 2 p̂ ∂ 2 p̂
+ =0 in η < 0, (97)
∂x2 ∂η 2

with boundary conditions

∂ p̂ ∂ Ĥ
p̂ = V Ĥ, =− on η=0 (98)
∂η ∂t

and
p̂ = o(−η) as η → −∞. (99)

(ii) Seek a separable solution to (97)–(99) of the form

   
p̂(x, η, t) = Re g(η) eikx+λt , Ĥ(x, t) = Re eikx+λt ,

where k is the (real) wavenumber, Re(λ(k)) is the growth rate and Re means ‘real part’. Deduce from
(97)–(99) that g(η) satisfies the second-order linear ordinary differential equation

g 00 − k 2 g = 0, (100)

with boundary conditions

g(0) = V, g 0 (0) = −λ; g(−∞) bounded. (101)

(iii) The problem (100)–(101) is an eigenvalue problem for λ(k): there is a nontrivial solution for g(η) if
and only if
λ = −V |k|. (102)

Determine the region of instability from the dispersion relation (102):

(a) V > 0 6= k ⇒ λ < 0 ⇒ perturbations decay, so interface linearly stable.
(b) V < 0 6= k ⇒ λ > 0 ⇒ perturbations grow, so interface linearly unstable.

• This is an example of the Saffman-Taylor instability: the interface between two immiscible fluids is unstable
when it moves toward the more viscous fluid.

• In the unstable regime experiments show that small perturbations grow into “viscous fingers” that evolve
in a highly nonlinear manner sensitive to the driving force, surface tension and substrate properties.

• This instability causes significant problems in e.g. the oil extraction industry because the flow of oil and
water through a porous media (such as rock) is governed by a similar set of differential equations (the
Darcy flow model: see e.g. Ockendon & Ockendon, chapter 5).

64
3.6 Thin films
3.6.1 Motivation

• Consider gravity-driven flow of a thin two-dimensional viscous layer of a priori unknown thickness h(x, t)
on a stationary horizontal plate at z = 0 below a vacuum:

−gk
z = h(x, t)

z=0

• Such a thin film may model spreading of a raindrop, molten lava or an ice sheet.

• This is a free boundary problem, as fluid domain 0 < z < h(x, t) must be determined as part of solution.

• New features: need to incorporate gravity into the incompressible Navier-Stokes equations and prescribe
appropriate free boundary conditions at the free surface z = h(x, t).

• Key idea: assume aspect ratio (depth / length) is small and apply lubrication theory.

3.6.2 Dimensional problem

• The flow is governed by the incompressible Navier-Stokes equations (22)–(23) with F = −gk, so that

Du
ρ = −∇p + µ∇2 u − ρgk, ∇ · u = 0,
Dt

where ρ is the density, µ is the viscosity, g is the acceleration due to gravity, u(x, z, t) = ui + wk is the
velocity, p(x, z, t) is the pressure and the convective derivative

D ∂ ∂ ∂ ∂
= +u·∇= +u +w .
Dt ∂t ∂t ∂x ∂z

• On the substrate impose the no-flux and no-slip boundary conditions

u = 0, w = 0 on z = 0.

• On the free surface z = h(x, t) need to impose three boundary conditions to determine its location: the
simplest are

∂h ∂h
no-flux (1 scalar boundary condition) : w= +u ,
∂t ∂x

zero-stress (2 scalar boundary conditions) : t(n) = 0,

where t(n) = ei σij nj is the stress exerted on the free surface by fluid and n is the downward pointing unit
normal, i.e.
hx e1 − e3
n= .
(1 + hx2 )1/2

65
3.6.3 Nondimensionalization
• We make the usual lubrication scalings:

x = L x̂ (L = [x] prescribed initially),

h = δL ĥ (δ = [h]/[x] prescribed initially),
z = δL ẑ ([z] = [h])
u = U û (U = [u] to be determined),
w = δU ŵ ([w]/[z] = [u]/[x] by cty. eq.),
L
t = t̂, ([h]/[t] = [w] by kinematic boundary condition),
U
µU
p = p̂ ([p]/[x] = µ[u]/[z]2 by x-mom. eq.),
δ2L
µU
σij = σ̂ij ([σij ] = [p]),
δ2L
µU
t = t̂ ([t] = [p]).
δ2L

3.6.4 Dimensionless problem

• Substitute scalings into the two-dimensional Navier-Stokes equations and simplify to obtain

Dû ∂ p̂ ∂ 2 û ∂ 2 û
δ 2 Re = − + δ2 2 + 2 ,
Dt̂ ∂ x̂ ∂ x̂ ∂ ẑ
Dŵ ∂ p̂ ∂ 2 ŵ ∂ 2 ŵ δ 3 ρgL2
δ 4 Re = − + δ4 2 + δ2 2 − ,
Dt̂ ∂ ẑ ∂ x̂ ∂ ẑ µU
∂ û ∂ ŵ
0 = + ,
∂ x̂ ∂ ẑ
where the Reynolds number Re = ρLU/µ and dimensionless convective derivative

D ∂ ∂ ∂
= + û + ŵ .
Dt̂ ∂ t̂ ∂ x̂ ∂ ẑ

• Key step: assume gravity drives flow by setting

δ 3 ρgL2
U= .
µ

• On the substrate the no-flux and no-slip boundary conditions are unchanged:

û = 0, ŵ = 0 on ẑ = 0.

• On the free surface the kinematic condition is unchanged

∂ ĥ ∂ ĥ
ŵ = + û on ẑ = ĥ(x̂, t̂).
∂ t̂ ∂ x̂

• The zero-stress condition becomes

t̂(n) = 0 on ẑ = ĥ(x̂, t̂),
where
t̂(n) = e1 (σ̂11 n1 + σ̂13 n3 ) + e3 (σ̂31 n1 + σ̂33 n3 ) ,

66
 with
∂ û ∂ û ∂ ŵ ∂ ŵ
σ̂11 = −p̂ + 2δ 2 , σ̂13 = σ̂31 = δ + δ3 , σ̂33 = −p̂ + 2δ 2 ,
∂ x̂ ∂ ẑ ∂ x̂ ∂ ẑ
and
δ ĥx̂ 1
n1 =  1/2 , n3 = −  1/2 .
1 + (δ ĥx̂ ) 2 1 + (δ ĥx̂ ) 2

• Hence, the stress vector

" ! #
∂ ĥ ∂ û
t̂(n) ∼ −δ p̂ + e1 + p̂ e3 as δ → 0.
∂ x̂ ∂ ẑ
ẑ=ĥ

3.6.5 Lubrication approximation

• Assume that both the aspect ratio δ and the reduced Reynolds number δ 2 Re are small, i.e.
δ 5 gL3
δ, δ 2 Re = 1,
ν2
to obtain at leading order the quasi-steady lubrication equations (dropping the hats ˆ):
∂p ∂ 2 u
0 = − + , (103)
∂x ∂z 2
∂p
0 = − − 1, (104)
∂z
∂u ∂w
0 = + , (105)
∂x ∂z
• On the substrate, the no-flux and no-slip boundary conditions are unchanged:
u = 0, w = 0 on z = 0. (106)

• On the free surface, the kinematic and zero-stress boundary conditions become

∂h ∂h 
w = +u (no flux), 

∂t ∂x 


p = 0 (zero normal stress) on z = h(x, t). (107)


∂u 


= 0 (zero shear stress) 
∂z
3.6.6 Cross-layer averaged conservation of mass expression
• Integrate the incompressibility condition (105) from z = 0 to z = h:
Z h Z h  " #z=h
∂u ∂w ∂ ∂h
0= + dz = u dz − u + w ,
0 ∂x ∂z ∂x 0 ∂x z=h z=0
by Leibniz’s integral rule.
• Applying the no-flux boundary conditions on z = 0, h in (105)–(106), we obtain
Z h 
∂ ∂h
0= u dz + .
∂x 0 ∂t
• Hence, defining the mean x-velocity by
Z
1 h
ū := u dz, (108)
h 0
we obtain the cross-layer averaged conservation of mass expression
∂h ∂
+ (hū) = 0. (109)
∂t ∂x

67
3.6.7 Integrating the momentum equations
• To determine the mean velocity ū, need to integrate the momentum equations to obtain u.

• This is possible because z-momentum equation (104) implies that

p = h(x, t) − z

by the normal stress condition in (107).

• At leading order the pressure is hydrostatic, receiving a contribution from weight of fluid above only.

• For fixed (x, t), x-momentum equation (103) is an ordinary differential equation in z for u, subject to the
no-slip condition in (106) and the zero-shear-stress condition in (107), giving after two integrations the
expression  
∂h 1 2
u= z − hz .
∂x 2

• By (108), the mean velocity

Z h  
1 ∂h 1 2 h2 ∂h
ū = z − hz dz = − .
h 0 ∂x 2 3 ∂x

3.6.8 Thin-film equation for h(x, t)

• Hence, the cross-layer averaged conservation of mass expression (109) becomes
!  3 
∂h ∂ ∂ h ∂h
=− hū = ,
∂t ∂x ∂x 3 ∂x

which is a nonlinear version of the heat equation.

• Coefficient of highest derivative (the mobility h3 /3) vanishes when h = 0, so the partial differential equation
is degenerate at h = 0.

• Two important manifestations of this degeneracy are the existence of compactly supported mass preserving
similarity solutions and “waiting time” behaviour (see B5a).

68