Flnotes
Flnotes
Flnotes
Dr Mark Blyth
Contents
0. Introduction
1. Preliminary ideas
2(a) Stress in fluids
2(b) Strain in fluids
2(c) Stress-strain relationship for a Newtonian fluid
3. The Navier-Stokes equations
4. Vorticity dynamics and the stream function
5. Exact solutions of the Navier-Stokes equations
6. Stokes flow
7. Lubrication theory
8. Boundary layers
9. Hydrodynamic stability
Appendix: Invariants of a matrix under rotation.
1
0. Introduction
The study of fluid mechanics encompasses many varied applications in both natural and
engineering processes. Some examples are
Aerodynamics – flow past wings; computation of drag and lift; D’Alembert’s paradox
Biological fluid dynamics – fish swimming; flow in arteries and veins; blood flow in
vascular tumours; Korotkoff sounds
Geophysical fluid dynamics (GFD) – oceanography and meteorology; climate and
weather prediction
Industrial fluid dynamics – flow past structures and buildings; coating flows (e.g. depo-
sition of photographic films); spraying processes
Magnetohydrodynamics (MHD) – flow of conducting fluids; Maxwell’s equations
Physiochemical fluid dynamics – flow under surface tension in the presence of chemical
contaminants
Turbulence – transition to turbulence from laminar flow; boundary layer receptivity
1. Preliminary ideas
What is a fluid?
A simple fluid cannot withstand an externally applied stress without deforming in some
manner. This definition of a fluid encompasses both liquids and gases. In order to develop
a mathematical basis for the study of fluid motion, we adopt the continuum hypothesis,
2
which states that the material under scrutiny is assumed to behave as if it is perfectly con-
tinuous. In real materials, fluctuations at the molecular level mean that this hypothesis is
inaccurate at very small scales. However, in most situations we are concerned with macro-
scopic observations of a given fluid motion and are not troubled by the irregular variations
apparent on molecular lengthscales. Focussing on a tiny volume of a given fluid, we make
the approximation that physical quantities, such as mass and momentum, are spread uni-
formly over that volume rather than unevenly distributed within it. Under this assumption
we may define, for example, the fluid velocity u(x, t), density ρ(x, t), or temperature T (x, t),
valid at a point x and time t in the flow field and treat them as continuous, differentiable
quantities over the domain of interest.
A viscous fluid is distinguished by its ability to resist shearing motions. A perfect or inviscid
fluid offers no resistance to shear. The viscosity, µ, of a fluid is a measure of its ability to
withstand such shearing action. The kinematic viscosity ν of a fluid is defined to be the
ratio between the viscosity and the density, so ν = µ/ρ.
Compressibility
What role do compressibility effects play in determining the motion of a fluid? Air is clearly
much more easily compressible than water, say, but to what extent is this important in a
particular flow?
Bernoulli’s theorem states that in a steady, inviscid flow with no body forces acting, the
pressure is related to the velocity via
1
p + ρu2 = constant along a streamline.
2
It follows that a change in speed of a fluid element from zero to U requires a pressure change
3
of O(ρU 2 ).
By definition,
dp
c2 = ,
dρ
where c is the speed of sound. Thus a small change in density is related to a small change
in pressure via
dρ dp
= 2,
ρ ρc
and so ρU 2
dρ
2
=O = O M ,
ρ ρc2
on defining the Mach number
U
M= .
c
We can deduce that if M ≪ 1, the effects of compressibility can be safely ignored.
Examples:
1. In air, c ≈ 340ms−1 , so that for U = O(102 ms−1 ), M = 0.29 and compressibility will
not be important.
2. For water, c ≈ 1.5 × 103 ms−1 , so would need U ≫ O(1500ms−1 ) for significant com-
pressibility effects.
3. For an incompressible fluid, changes in pressure occur with no corresponding change in
density. In this case dp/dρ is infinite, so the speed of sound in the incompressible medium is
infinite. A consequence of this is that the fluid feels the effect of a local change everywhere
in the fluid instantaneously.
The forces acting on a fluid in motion are divided into two categories:
1. Long-range forces: Externally applied forces such as gravity, an imposed pressure gra-
dient, or electromagnetic forces in the case of conducting fluids.
2. Short-range forces: These refer to local dynamic stresses developing within the fluid
itself as it moves; specifically, the forces acting on a given element of fluid by the
surrounding fluid.
4
To describe those in the second category, we introduce the stress tensor σij , defined as
the i-component of stress acting on an element of surface with unit normal n
in the j-direction.
where i, j both run over 1,2,3, corresponding to x, y, z in a Cartesian frame.
This makes the stress on a surface easy to find if the surface in question has normal pointing
in one of the coordinate directions. We just pick out the relevant entry of σij . For example,
the stress on a plane wall at y = 0 is found by noting that the normal points in the y
direction. So, for the component of stress pointing in the x direction, for example, we need
σ12 .
Similarly, the component of stress in the z direction on a wall at x = 0 is σ31 .
In general, of course, the unit normal will not point in a coordinate direction. Suppose
instead we have a solid boundary over which flows a fluid. We show below that the force
exerted at a point on the boundary by the passing fluid is
3
X
σij nj
j=1
where the unit vector n is normal to the boundary and points into the fluid.
5
F
n
2
0000000
1111111
1111111
0000000
0000000
1111111
0000000
1111111
0000000
1111111
0000000
1111111
1 1111111
0000000 dS
A3 dS1
n
dS2 dS
P x2
A2
A1 dS3
x1
We let
A1 A2 A3
dS
dS
P A2 A3 dS 1 n1 dS
the area of be = . (2.2)
P A1 A3 dS 2 n dS
2
P A1 A2 dS 3 n3 dS
(Note: To find the areas of each surface of the tetrahedron, use the simple half-base times
height rule.)
Now,
dS
Hi (x, n) dS
dS 1 −Hi (x, N1 ) dS 1
the force exerted by the fluid outside dV on is .(2.3)
dS 2
−H i (x, N2 ) dS 2
dS 3 −Hi (x, N3 ) dS 3
Note: Here we are using the Einstein summation convention of summing over a repeated
index. So σij nj means σi1 n1 + σi2 n2 + σi3 n3 .
6
Newton’s second law states that
If the dimensions of the tetrahedron are of size O(ε), the left hand side of (2.4) is O(dS) =
O(ε2 ), while the right hand side is O(dV ) = O(ε3 ). Therefore, if we let ε tend to zero, (2.4)
produces the result that Force = 0, and we conclude that and so
Example 1: Find the stress in the x direction on a wall placed at y = 0. In this case
n1 = n2 = 0, n3 = 1. So Fi = σij nj = σi2 . The component in the x direction is F1 = σ12 ,
as we had above.
Example 2: Find the stress in the z direction on the upper side of the plane x + y + z = 1
due to the fluid above it. In this case
∇(x + y + z − 1) 1
n=+ = √ (1, 1, 1).
|∇(x + y + z − 1)| 3
Then the stress in the z direction is F3 , with
1 1 1
F3 = σ3j nj = σ31 √ + σ32 √ + σ33 √
3 3 3
1
= √ (σ31 + σ32 + σ33 ).
3
We know that
X d
Moments = {Angular momentum}.
dt
7
The left hand side is made up of moments due to
while the right hand side is of size (cube volume) × (cube length) = O(ε4 ). Therefore, on
the face shown in the figure,
ε
2 ε2 [σ12 − σ21 ] = O(ε4 ) =⇒ σ12 = σ21 .
2
Following similar arguments, we deduce that
σij = σji .
Consider the relative positions of two individial fluid particles under any particular fluid
motion. At time t, the particles are close together, with the first at position x and moving
with velocity u(x, t), and the second at x + dx, moving with velocity u(x + dx, t).
u(x+dx,t)
x+dx
x
u(x,t)
Expanding the second particle’s velocity in a Taylor series, we have the following expression
for the relative velocity of the two particles:
∂ui
ui (x + dx, t) − ui (x, t) = dxj + O(dx2j ). (2.6)
∂xj x
We now re-write the derivative on the right hand side like this
∂ui
= eij + ξij
∂xj
where
! !
1 ∂ui ∂uj 1 ∂uj ∂ui
eij = + , ξij = − − . (2.7)
2 ∂xj ∂xi 2 ∂xi ∂xj
8
Note that eij is symmetric, while ξij is antisymmetric (So eij = eji , ξij = −ξji ).
eij is called the rate of deformation or strain tensor.
ξij is sometimes called the vorticity tensor.
To see why these two tensors are given their names, we examine their meanings in terms of
the kinematics of the fluid particles.
v3
v2
y
v1
x
x’
A set of Cartesian axes (x, y, z) and a rotated set of principal axes (x′ , y ′ , z ′ )
along which are aligned the three orthogonal eigenvectors of e′ij .
We denote it as e′ij in the new frame. Writing both out as a matrix we have
e11 e12 e13 e′11 0 0
eij = e21 e22 e23 and e′ij = 0 e′22 0
e31 e32 e33 0 0 e′33
The three eigenvalues of the new matrix are λ1 = e′11 , λ2 = e′22 and λ3 = e′33 and, since e′ij is
symmetric, they are all real and have mutually orthogonal eigenvectors. These eigenvectors,
v1 , v2 , v3 are aligned along the principal axes and form a basis.
The action of matrix e′ij on a general vector, q = 3i=1 αi vi , for constants αi , is
P
3
X 3
X
′ ′
e ·q = αi e · vi = λi αi vi .
i=1 i=1
9
Local expansion of the fluid
Suppose we try the above ideas out on a small cube of fluid aligned with the principal axes
and with sides of length dx′1 , dx′2 , dx′3 and volume V = dx′1 dx′2 dx′3 .
z’
dx 3
y’
dx 1
dx 2
x’
After a short time dt, the side dx′1 will have become deformed so that its new length is
old length + (rate of stretching in the x′1 direction)×dt
So its new length will be (1 + e′11 dt) dx′1 . Similar arguments apply for the other sides so the
new edges are of lengths
{(1 + e′11 dt)(1 + e′22 dt)(1 + e′33 dt) − 1}V = (e′11 + e′22 + e′33 )V dt = Trace(e′ij )V dt,
neglecting terms of size O(dt2 ). Now, the trace of a matrix is invariant under a rotation to
a new coordinate frame (see Appendix A). So, the change in volume (divided by V dt) is
equal to
∂ui
Trace(e′ij ) = Trace(eij ) = eii = = ∇ · u.
∂xi
A fluid is defined to be incompressible if the local change in volume of a fluid element such
as the box just considered is zero. In other words
10
so, for example, ε123 = ε231 = 1, ε132 = −1, ε133 = 0. The alternating tensor provides a
convenient way of expressing a vector (cross) product in index notation.
Exercise: Show that εijk aj bk is equivalent to a × b by writing out the components.
To understand the kinematics associated with ξij , we reason as follows. Since ξij is antisym-
metric, it has only 3 independent components and may therefore be written quite generally
as
The factor −1/2 is included for convenience. Check for yourself that the right hand side of
(2.8) has only 3 independent components, ωk .
Referring back to the Taylor expansion (2.6), the antisymmetric contribution to the relative
velocity of the two particles is
dx
Local rotation of two fluid particles at angular rate 12 |ω| about an axis pointing
in the direction ω.
What are the components of ω in terms of the fluid velocity? First we note the useful
identity involving the alternating tensor,
11
where the Kronecker delta δij is 1 if i, j are equal and 0 otherwise. Hence
Now, operating with εijl on both sides of (2.8) and using the definition (2.7), we find
∂u ∂ui ∂uj
j
ωl = −εlij ξij = 12 εlij − = εlij i.e. ω = ∇ × u.
∂xi ∂xj ∂xi
The vector ω is called the vorticity. Any flow for which ω = 0 is called irrotational.
To proceed with the modelling of fluid flow, we need to understand the relationship between
the stress in a fluid and the local rate of strain. In other words, how much is the fluid being
deformed by the stresses exerted by the surrounding fluid?
For a fluid at rest, there are no viscous forces, and the stress is simply related to the static
pressure p by
So σii = −p δii = −3p and therefore p = −σii /3. Hence the stress exerted on an element
with unit normal n is σij nj = −pni and has magnitude p in the −n direction. So the
pressure force exterted on a wall by a fluid at rest acts in a direction perpendicular to that
wall.
For a fluid in motion, we write
where p is the dynamic pressure and dij is the deviatoric stress tensor.
We define dii = 0 so that once again σii = −3p and thus p = −σii /3. However, note that
since the fluid is now moving, p in equation (2.10) is not in general the same as p in equation
(2.9).
The deviatoric stress dij is due to the motion of the fluid, which consists of translation,
rigid rotation and stretching or contracting.
Recall that ξij corresponds to rigid rotation.
eij corresponds to stretching/straining.
Newtonian fluid
For a Newtonian fluid, the deviatoric stress depends only on eij and the dependence is
linear.
12
Notes
1. For a Newtonian fluid dij is therefore symmetric.
2. The instantaneous stress at any point in the fluid does not depend on the past history
of the fluid motion. Non-Newtonian, visco-elastic fluids can have a “memory” whereby the
stress depends on the state of the fluid at previous times.
Therefore,
n o
dij = µ[δik δjl + δil δjk ] + λδij δkl ekl
= 2µeij + λδij ekk . (2.12)
Remember that ekk = ∇ · u represents the amount of dilatation or local expansion of the
fluid. Therefore, for an incompressible fluid, ekk = 0 and so
13
This provides a constitutive relation between the stress and strain for an incompressible,
Newtonian fluid. The constant µ is called the viscosity of the fluid. We assume equation
(2.15) throughout the rest of the course.
Our ultimate goal is to derive an equation describing the motion of a fluid. To move towards
this goal, we introduce the important concept of a derivative which tells us how things
change as we move with the fluid. The full usefullness of this idea will become apparent
later on when we apply Newton’s second law to derive the equations of fluid motion.
u
z
y
At rest
Consider a set of axes moving with the local fluid velocity u and another set fixed in space
as shown in the diagram. If the origin of the moving set of axes is at x = (x(t), y(t), z(t))
at time t, we can write its velocity as
dx
u(x, t) =
dt
Suppose now we wish to calculate the acceleration of this set of axes. By definition, the
acceleration is given by
14
on writing dx = dt.(dx/dt). Also by Taylor’s theorem,
∂u
u(x + dx, t + dt) = u(x, t) + dx · ∇u + dt + ··· . (3.3)
∂t
d2 x du 1 n ∂u o
2 = dt = lim dt + dx · ∇u + · · ·
dt dt→0 dt ∂t
∂u dx
= + · ∇u
∂t dt
du ∂u
Thus = + u · ∇u. (3.4)
dt ∂t
The derivative on the left hand side of (3.4) is often written with a large D and referred to
as the convective derivative (sometimes the material or substantial derivative).
We will follow this convention and henceforth write
Du ∂u
= + u · ∇u (3.5)
Dt ∂t
The convective derivative may be applied to any scalar or vector function. For example,
consider the function f (x(t), y(t), z(t)). Proceeding as above, we find
Df ∂f
= + u · ∇f
Dt ∂t
∂f ∂f ∂f ∂f
= +u +v +w .
∂t ∂x ∂y ∂z
This is really just the statement of the chain rule:
Df df ∂f ∂f dx ∂f dy ∂f dz ∂f ∂f ∂f ∂f
≡ = + + + = +u +v +w . (3.6)
Dt dt ∂t ∂x dt ∂y dt ∂z dt ∂t ∂x ∂y ∂z
Df /Dt tells us how the function f changes as we move with the flow.
For example, Dρ/Dt expresses the amount by which the density of a fluid particle changes
as it is convected along with the flow. For an incompressible fluid, we therefore have
Dρ
= 0.
Dt
15
Working in the Lagrangian framework, we consider the motion of individual point particles;
that is, given a fluid particle at a particular point in space at a particular time, we monitor
how that particle travels through space as time progresses. To this end, we label all particles
over the entire fluid domain according to their position at time t = 0. Then, at a later time,
an individual fluid particle is identified by its label, as if it was carrying a flag along with
it as it moves in the flow.
We therefore ascribe to each particle a label a, where a is the position vector of that particle
at t = 0.
Any position in the flow-field x may be regarded as a function of its label and time. Moreover
we can write
where χt (a) is a mapping from the initial particle positions at t = 0, to their current
positions at time t. The subscript on the mapping function indicates that the particles are
being mapped from where they were at t = 0 to their location at time t.
χ 1010
t
x ( a,t)
0011 a
Consider now the differential vector (da1 ,0,0) in labelling space and denote its counterpart
in convected space as dx(1) . Using (3.7) and applying the chain rule we can write:
∂χt
dx(1) = da1 . (3.8)
∂a1
16
At this point, we note that the expression
represents a volume element in the convected space. Convince yourself that this is true.
Therefore, using (3.8) we can write
∂χ ∂χt ∂χt
t
(dx(1) × dx(2) ) · dx(3) = × · da1 da2 da3 .
∂a1 ∂a2 ∂a3
Using the notation χit to denote the ith component of χt , we note that
∂χ ∂χt ∂χt ∂χi
t t
× · = det ,
∂a1 ∂a2 ∂a3 ∂aj
dV c = J dV l , (3.9)
where
∂χ1t ∂χ1t ∂χ1t
∂a1 ∂a2 ∂a3
∂χ2t ∂χ2t ∂χ2t
∂χi
t
J = det = (3.10)
∂aj ∂a1 ∂a2 ∂a3
∂χ3t ∂χ3t ∂χ3t
∂a1 ∂a2 ∂a3
we have immediately managed to change an integral defined over a volume which is changing
in time to one over a volume which is fixed. This will make differentiation a lot easier.
17
To find the rate of change of the parcel’s volume, we form
dV c d dJ DJ
Z Z Z
= J dV l = dV l = dV l , (3.12)
dt dt Vl Vl dt Vl Dt
since J is a function of x, the position in convected space, and so it follows by the definition
of the convective derivative that
dJ ∂J DJ
= + u · ∇J ≡
dt ∂t Dt
Note that we were able to take the d/dt inside the integral in (3.12) since Vl is independent
of time.
Equally, we may argue as follows: Since the particles making up the parcel move with the
00
11 V(t+dt)
00
11
111111111111
000000000000
c
000000000000
111111111111
000000000000
111111111111
00000000000011
00
volume swept
111111111111
out by surface 00
11
patch
V(t)
c
u n
fluid velocity u, the volume swept out by a small element dS c of its surface Sc in time dt
(see figure) must equal u · n dS c , where n is the surface unit normal, and thus
dV c
Z
= u · n dS c .
dt Sc
dV c
Z Z
= ∇ · u dV c = ∇ · u J dV l . (3.13)
dt Vc Vl
DJ
= J ∇ · u. (3.14)
Dt
We will find this useful in the proof of the following important theorem.
18
This result tells us how an integral of any scalar quantity defined over a time-dependent
volume of fluid V (t) changes in time. In fact, it states that
d DG
Z Z
G dV = + G ∇ · u dV . (3.15)
dt V (t) V (t) Dt
Proof
Using (3.9), we may convert the volume integral into one defined over the equivalent volume
in labelling space, Vl . So,
d d
Z Z
G dV = G J dV l
dt V (t) dt Vl
d DG
Z Z
G dV = J + GJ ∇ · u dV l (3.16)
dt V (t) V Dt
Z l
DG
= + G ∇ · u J dV l (3.17)
Dt
ZVl
DG
= + G ∇ · u dV , (3.18)
V Dt
using (3.9). ⋄.
Example
If we identify the scalar G(x, t) with the fluid density ρ(x, t), we can use (3.15) to write
d Dρ
Z Z
ρ dV = + ρ ∇ · u dV .
dt V (t) V (t) Dt
Since mass is conserved, it follows immediately that the integrand on the right hand side is
zero, so
Dρ
+ ρ ∇ · u = 0.
Dt
19
We noted above that if the fluid is incompressible, then
Dρ
= 0.
Dt
It therefore follows that, for an incompressible fluid,
∇ · u = 0,
which is a restatement of the result ekk = 0 found before using the small cube of fluid
argument.
Note: It is important to realise that incompressible does not mean constant density. This
is a common mistake. For example, consider the two-dimensional flow with velocity and
density fields
u = u(x, y)i + v(x, y)j, ρ = ρ(z).
Then
Dρ ∂ρ ∂ρ ∂ρ
= +u +v =0
Dt ∂t ∂x ∂y
so the condition for incompressibility is satisfied, but the density field varies with z.
We will make use of the Reynolds transport theorem when we derive the Navier-Stokes
equations.
Up to now, we have been concerned with the kinematics of a flow due to a predetermined
velocity field u(x, t). But how do we determine the underlying velocity field in the first
place?
To formulate the necessary equations governing the motion of a viscous fluid we will apply
the following fundamental principles:
1. Conservation of mass.
2. Newton’s second law of motion.
We will also discuss how energy is dissipated by viscous forces in a moving fluid and derive
an equation for the conservation of energy in a parcel of fluid.
20
where ρ is the fluid density, S is the volume surface, and n is the normal to V pointing
outwards. Using the fact that V is fixed in space and applying the divergence theorem, this
becomes
∂ρ
Z n o
+ ∇ · (ρu) dV = 0.
V ∂t
Since our choice of volume V was arbitrary, it must be true that
∂ρ
+ ∇ · (ρu) = 0. (3.19)
∂t
Equation (3.19) therefore expresses the principle of mass conservation. It is usually referred
to as the continuity equation, which will frequently be abbreviated to cty equation.
Using the definition of the convective derivative, we can re-express the continuity equation
as
Dρ
+ ρ∇ · u = 0. (3.20)
Dt
This is just the equation we found above using the Reynolds transport theorem.
d DG
Z Z
G dV = + G∇ · u dV .
dt V (t) V (t) Dt
for any scalar quantity G.
We are now in a position to derive the equation of motion of a fluid. Let F represent some
body force per unit mass acting on the fluid (i.e. a long-range force such as gravity). We
apply Newton’s second law of motion to a volume of fluid V (t), with surface area S(t) and
unit outward normal n, moving with the flow. Thus,
d
Z Z Z
ρui dV = ρ Fi dV + σij nj dS.
dt V (t) V (t) S(t)
The left hand side expresses the rate of change of momentum of the fluid; the right hand
side the forces acting on the fluid. The last term represents the stresses acting over the
surface of the chosen volume by the surrounding fluid.
Applying the Reynolds transport theorem to the first term and the divergence theorem to
the last, and using (3.20) we find
Dui ∂σij
Z Z Z
ρ dV = ρ Fi dV + dV . (3.21)
V (t) Dt V (t) V (t) ∂xj
21
Note that we have used the fact that
D(ρui ) Dui Dρ
=ρ + ui
Dt Dt Dt
and that Dρ/Dt = 0 since the fluid is incompressible.
If the fluid is Newtonian, we may use the constitutive relation introduced earlier,
2
σij = −p δij − µ∆δij + 2µeij ,
3
where
∆ = ∇ · u = ekk
represents the local expansion of the fluid.
Accordingly,
∂σij ∂p 2 ∂ ∂
=− − (µ∆) + 2 (µeij ).
∂xj ∂xi 3 ∂xi ∂xj
∆ = 0,
But ∂uj /∂xj ≡ ∇ · u = 0. Also, since our original volume was chosen arbitrarily, the
integrand in (3.23) must vanish, as the statement must be true for all possible volume
choices. Finally, we have
Dui ∂p ∂ 2 ui
ρ = ρ Fi − +µ 2 .
Dt ∂xi ∂xj
22
direct from (3.21).
In vector form, together with the continuity equation,
Du
∇ · u = 0, ρ = ρ F − ∇p + µ∇2 u. (3.25)
Dt
∂2 ∂2 ∂2
∇2 = + + .
∂x2 ∂y 2 ∂z 2
v2
1 2 u 2
ρ ut + uur + vuθ + wuz − = −pr + Fr + µ ∇ u − 2 − 2 vθ ,
r r r r
1 uv 1 2 2 v
ρ vt + uvr + vvθ + wvz + = − pθ + Fθ + µ ∇ v + 2 uθ − 2 ,
r r r r r
1
ρ( wt + uwr + vwθ + wwz ) = −pz + Fz + µ∇2 w,
r
u 1
ur + + vθ + wz = 0.
r r
∂2 1 ∂ 1 ∂2 ∂2
∇2 = + + + .
∂r 2 r ∂r r 2 ∂θ 2 ∂z 2
If the body force, F, is conservative, i.e. we can write it as the gradient of a scalar potential,
F = −∇Ω,
ρ F − ∇p = −∇(ρ Ω + p) = −∇P,
23
where we have defined the modified pressure P to be
P = p + ρ Ω.
Fluid-Solid
The usual inviscid no-penetration condition applies at a solid impermeable surface. This
states that u · n = 0, i.e. fluid cannot flow into the wall.
For a viscous fluid a further boundary condition is required to complete the flow description
due to the extra derivatives imported by the last term in (3.25). Experiments reveal that
this supplementary condition should stipulate that the fluid does not slip over the solid
boundary; that is to say the local fluid velocity at the wall is zero. Mathematically this is
stated as u · t = 0, where t is any vector tangential to the wall.
In summary, we require
Flow
y
x v
u
24
The no-slip condition states that u = v = 0.
Note: We can use the continuity equation to deduce a further piece of information about
the flow at z = 0. We know that
∂u ∂v ∂w
+ + = 0,
∂x ∂y ∂z
and, on z = 0,
∂u ∂u ∂v ∂v
= = = = 0,
∂x ∂y ∂x ∂y
since u and v are both zero everywhere on z = 0. Therefore,
∂w
=0 on z = 0.
∂z
25
The equation of the wall is z = sin t. Set f = z − sin t. Then,
Df ∂f ∂f ∂f ∂f
= + us + vs + ws
Dt wall ∂t ∂x ∂y ∂z
= − cos t + ws = 0,
=⇒ ws = cos t.
Fluid 1 t
n
f(x,y,z,t)=0
These state that the normal and tangential components of the velocities are continuous
across the interface.
Dynamic condition: We also need a condition expressing the balance of forces prevailing at
the interface between the fluids. To obtain this condition, we consider a small rectangular
volume of fluid of length ds as shown in the figure below.
Fluid 2
s
τ τ
t
ds
Fluid 1 n
The height of the rectangular volume is assumed negligible in comparison with its length
ds. If s measures arc length along the interface, the local unit normal and tangent vectors
to the interface are denoted as n(s) and t respectively at position s. Surface tension τ (s)
tugs on the volume at either end in the tangential direction. The force exerted on the lower
26
portion of the rectangular volume by fluid 1 is given by
(σ (1) · n) ds,
where σ (1) is the stress tensor in fluid 1. This follows from the definition of the stress tensor
given above. Similarly, the force exerted on the top by fluid 2 is given by −(σ (2) · n) ds.
The minus sign is needed since the unit normal must point into the fluid exerting the force.
Applying Newton’s second law to the rectangular control volume (and disregarding its
negligible inertia), we obtain the force balance
h i
τ (s + ds)t(s + ds) + τ (s)t(s) + ds σ(1) · n − σ(2) · n = 0.
The first two terms represent surface tension at either ends of the control volume. Since ds
is small, we may expand in Taylor series, so,
h dτ i h dt i h i
τ (s) + ds t(s) + ds + τ (s)t(s) + ds σ (1) · n − σ (2) · n + O(ds2 ) = 0.
ds ds
Expanding the first pair of square brackets and cancelling relevant terms,
dt dτ h i
ds τ + ds t + ds σ (1) · n − σ (2) · n + O(ds2 ) = 0.
ds ds
Dividing by ds, taking the limit as ds → 0, and defining the local surface curvature κ so
that
dt
κn = − ,
ds
we obtain,
dτ h i
−κτ n + t + σ(1) · n − σ (2) · n = 0.
ds
Rearranging, this yields
h i dτ
σ(1) − σ(2) · n = κτ n − t
ds
If the surface tension is constant, as it is generally speaking, the last term disappears. In
index notation, the jump in stress is then written as
(1) (2)
∆Fi ≡ [σij − σij ] nj = 2τ κ ni .
This states that the jump in stress ∆Fi at the interface is balanced by surface tension. A
more general derivation valid for a two-dimensional curved surface can be found in Batch-
elor. Refer to pages 64 and 69.
N.B. Surface tension is a very weak force. For example, at 20◦ C the surface tension between
air and water is
τ = 72.8 dyne/cm.
27
p (1)
p (2)
n
28
where ρF is the body force per unit mass. To proceed, we take the dot product of this
equation with the velocity u:
Du
ρu · = −u · ∇p + ρF · u.
Dt
Integrating over a volume V moving with the fluid, we find
Du
Z Z Z
ρu · dV = − u · ∇p dV + ρF · u dV .
V Dt V V
So,
D 1 2
Z Z Z
ρu dV = − u · ∇p dV + ρF · u dV .
V Dt 2 V V
Note that
∇ · (pu) = p∇ · u + u · ∇p = u · ∇p
since ∇ · u = 0, and so,
d 1 2
Z Z Z
ρu dV = − ∇ · (pu) dV + ρF · u dV .
dt V 2 V V
on using the Reynolds Transport Theorem on the LHS. Using the divergence theorem, we
have finally
d 1 2
Z Z Z
ρu dV = (−pn) · u dS + ρF · u dV . (3.26)
dt V 2 S V
Now, since −pn is the pressure force acting on the blob due to the surrounding fluid, the
first term on the RHS represents the rate of working of the external pressure force on the
blob. Similarly, the second term on the RHS represents the rate of working of the body
force on the blob. Recognising
1 2
Z
ρu dV = K
V 2
as the total kinetic energy in the blob, we can read the equation as meaning that in each
unit of time, the
Increase of K.E. of blob = work done by pressure force + work done by body force
In other words, no energy at all is lost. All work put in by the pressure force and the body
force goes directly into increasing the kinetic energy of the blob. This remark does not hold
true if viscosity is included, as we shall see in a moment.
Firstly, we recall that equation (3.26) applies for a volume V (t) moving with the fluid. If
instead we consider a volume V which is fixed in the fluid, we proceed as follows.
29
As before, we take the dot product of the Euler equation with the velocity u:
Du
ρu · = −u · ∇p + ρF · u.
Dt
Integrating over the fixed volume V , we find
Du
Z Z Z
ρu · dV = − u · ∇p dV + ρF · u dV .
V Dt V V
For a steady flow with a conservative body force such that F = −∇V , this becomes
1 2
Z Z Z
ρu · ∇ |u| dV = − u · ∇p dV + ρu · ∇V dV ,
V 2 V V
where we have used the fact that u · ∇u = ∇(u2 /2) for a divergence free velocity field. Next
we note that
1 2 1 2
ρu · ∇ |u| = ∇ · ρ|u |u , u · ∇V = ∇ · (V u),
2 2
Therefore, Z
1 2
ρ|u | + p + ρgy u · n dS = 0
S 2
is the statement of conservation of energy over a closed volume fixed in space1 .
where ρFi is the body force per unit mass (see equation 3.24).
As above, we take the dot product with u. Skipping steps identical to before (albeit this
time in index notation) we find
dK ∂σij ∂σij
Z Z Z
= ui ρFi + dV = ρ u · F dV + ui dV .
dt V ∂xj V V ∂xj
1
Compare Longuet-Higgins (1975), Proceedings of the Royal Society of London, 342, equation 1.7).
30
Note that the effect of the pressure is here contained within the σij of the final term since
This term therefore also incorporates the new effect due to viscosity.
Now,
∂σij ∂(ui σij ) ∂ui
= − σij ,
∂xj ∂xj ∂xj
and so
∂σij ∂(ui σij ) ∂ui
Z Z Z
ui dV = dV − σij dV
V ∂xj V ∂xj V ∂xj
∂ui
Z Z
(using divergence theorem) = ui σij nj dS − σij dV .
S V ∂xj
But, writing
∂ui 1 ∂ui ∂ui
= + ,
∂xj 2 ∂xj ∂xj
we have
∂ui 1 ∂ui ∂ui 1 ∂ui 1 ∂ui 1 ∂ui 1 ∂uj
σij = σij + = σij + σij = σij + σij
∂xj 2 ∂xj ∂xj 2 ∂xj 2 ∂xj 2 ∂xj 2 ∂xi
on swapping i and j in the second term and using the fact that σij is symmetric. Thus,
To summarize, we have
dK
Z Z Z
= ρ u · F dV + ui σij nj dS − σij eij dV . (3.27)
dt V S V
Also,
σij eij = −peii + 2µeij eij = 2µe2ij ,
since eii = ∇ · u = 0.
Therefore, finally, the rate of change of the volume’s kinetic energy
dK
Z Z Z Z
= ρ u · F dV + (−pn) · u dS + 2µ ui eij nj dS − 2µ e2ij dV . (3.28)
dt V S S V
31
Compare this with equation (3.26) for inviscid flow. If we set µ = 0, equation (3.28) reduces
to (3.26).
The first new term in (3.28), Z
2µ ui eij nj dV
V
is the rate of working of the deviatoric viscous stresses on the surface of the fluid blob V .
Thus the viscous stress makes a contribution to the increase of kinetic energy in the blob,
as we would have expected.
So far so good – all of the first three terms on the RHS of (3.28) represent work done on
the blob to increase its kinetic energy. But...
...the second new term in (3.28), Z
−2µ e2ij dV
V
is definitely negative and marks a reduction in kinetic energy. So some energy is being lost.
We say is it lost due to viscous dissipation and define the dissipation function Φ, where
Φ = 2µeij eij .
This represents the rate at which energy is dissipated per unit volume by viscous action.
For a compressible fluid we find
∆2
Φ = 2µ eij eij − , (∆ = eii )
3
∆δij 2
(non-trivial step) = 2µ eij − .
3
Summary
No energy is lost in an inviscid flow. This does not include such phenomena as hydraulic
jumps and breaking waves - in these examples discontinuities are introduced into the flow,
the preceding arguments do not work, and energy is in general lost.
Energy is lost from a viscous flow. The rate of loss of energy is represented by the dissipation
function, Φ.
32
Consider a fluid moving with no body forces acting. The total energy, E, inside a moving
volume of fluid, V , is given by
E = K + I,
where K is the kinetic energy, and I is the internal energy associated with the movement
of the constituent molecules within the fluid. (Note there is no potential energy because
we’ve left out body forces.) When work is done on the volume by the surrounding fluid,
some of the work goes into increasing the kinetic energy (K) and some into increasing the
internal energy (I).
From equation (3.27), we know that
dK
Z Z Z Z
= ui σij nj dS − 2µe2ij dV = ui σij nj dS − Φ dV ,
dt S V S V
The rate of change of the total energy E in the moving volume must be balanced by the
rate of working of the local stresses on the volume surface, so that
dE
Z
= ui σij nj dS.
dt S
33
The First Law of Thermodynamics dictates that, in a unit of time, the internal energy of
the volume is increased by (a) any work done, W , on the volume, and (b) any heat, Q,
entering the volume2 . In symbols,
∆E = Q + W, (3.29)
where ∆E means the change in internal energy per unit time. From above, we have over
the time interval dt,
Z Z Z
Q = dt q dS = dt κ n · ∇T dS = dt κ ∇2 T dV ,
S S V
Z Z
W = dt 2µe2ij dV = dt Φ dV .
V V
where c is the specific heat capacity of the fluid and T is the temperature of the fluid. So,
according to the First Law of Thermodynamics (3.29),
dI d
Z Z
2
= ρ c T dV = κ∇ T + Φ dV .
dt dt V V
Conclusion: The rate of increase of temperature as we move with the fluid is proportional
to the rate of loss of energy due to viscous dissipation.
34
3.7 The momentum integral equation
It is sometimes convenient to consider the balance of forces acting on a given parcel of fluid.
To this end, it is useful to consider a volume of fluid which is fixed in space, as we did for
the conservation of mass argument above. However, instead of considering the mass flux
into and out of the fixed volume, we now examine the flux of momentum.
Newton’s second law states that the force acting on the fixed volume V is equal to rate of
change of its momentum. Hence3 ,
∂
Z Z Z Z
(ρ u) dV = − ρ u(u · n) dS + ρF dV + σ · n dS. (3.30)
V ∂t S V S
In turn:
The term on the left hand side is the rate of change of momentum inside the fixed volume.
The first term on the right hand side expresses the flux of momentum through the surface
S. Since the volume is fixed, some fluid particles are entering and some are leaving the
volume, carrying momentum with them.
The second term on the right hand side expresses the total body force acting on the fixed
volume.
The third term on the right hand side expresses the total of the surface forces acting on S
due to the surrounding fluid.
For an inviscid fluid, σ = −pI, where I is the identity matrix4 . For a steady flow in the
absence of a body force, we then have
Z
ρ u(u · n) + pn dS = 0.
S
This relation may also be obtained by integrating the Euler equations directly over a volume
V and using appropriate vector identities.
35
Ui
direction at constant pressure P . Let D be the total force on the object. Then
Z
D=− p n dS,
P
Assuming that the only body force acting is gravity, F = −gk, and writing
F = ∇ϕ, ϕ = −gz,
since gravity is a conservative force, the volume integral becomes a surface integral courtesy
of the divergence theorem,
Z Z Z Z
ρF dV = ρ∇ϕ dV = ρϕ n dS = − ρgz n dS,
V V S S
n n
n
P
Ui n
B
H
F
36
Assuming that the bounding box is sufficiently large that the disturbance caused by the
object has decayed and the flow has settled to a uniform stream, then u · n = 0 on U and L,
and u · n = U on F and u · n = −U on B. Also, since the fluid can’t penetrate the object,
u · n = 0 on P . So the first integral in (3.32) gives
Z Z Z
− ρ u(u · n) dS = − ρ U 2 i dS + ρ U 2 i dS = 0.
Ω F F
To determine the second integral, we note that Bernoulli’s equation in the fluid states that,
along a streamline,
1
p + ρu2 + ρgy = C,
2
for constant C. Letting the size of the box approach infinity we see that
pU − pL = −ρgH
pB = pF ,
where A is the area of the upper (or lower) face and VP is the volume occupied by the
object.
So (3.32) gives
0 = D − ρgH k + ρg[HA − VP ] k.
and thus
D = ρgVP k.
37
So the total force just involves the Archimedes upthrust of the solid object in the vertical
direction. Taking the horizontal component in the direction of the flow,
D · i = 0,
and we predict zero drag on the object!
This is D’Alembert’s Paradox. It states that an object of an shape placed in a uniform
stream in an inviscid fluid experiences zero drag.
It is important to realise that D’Alembert’s paradox does not mean than any object placed
in any inviscid flow experiences zero drag. The previous argument applies for an object
placed in a uniform stream. If there is a free surface, for example, an object such as a ship
may experience drag to wave resistance, as is explored in the next section.
U
H z
x P
W
The wall is at z = 0 and the undisturbed height of the water is H. If the Froude number,
F , defined as
U
F = ,
(gH)1/2
is less than 1 then the free surface of the stream, T , may develop waves as shown in the
figure.
Let us calculate the drag on the object P for this flow. The inviscid steady form of the
momentum integral equation (3.30) in two-dimensions is
Z Z Z
0=− ρ u(u · n) dC − p n dl + ρF dA,
C C A
over an area A enclosed by the closed contour C of incremental arc length dl. Assuming
that the only body force acting is gravity, F = −gk, and writing F = −∇(gz) as in the
previous section, the equation becomes
Z Z Z
0=− ρ u(u · n) dl − p n dl − ρgz n dl. (3.33)
C C C
38
Take S to be the dotted contour shown in the figure below. Note that there is no direction
associated with the contour as we are integrating with respect to arc length.
T
d
n
n n
F n B
n P
W
The contour comprises the ends F and B, the wall W , the object P , and the free surface
up to one of the crests T . Without loss of generality, we take p = 0 on the free surface, T .
Noting that u · n = 0 on the object, on the free surface T , and on the wall W , and that nx
is zero on the wall, where n = (nx , ny ), the horizontal component of (3.33) gives
Z Z Z
−Dx = −i · D = −ρ u(u · n) dl − p nx dl − ρg z nx dl,
F,B F,B P,F,B,T
where Z
D=− p n dl
P
is the hydrodynamic force on the object. Now,
1
Z Z
z nx dl = ηη ′ dx = [η 2 ]L
−L = 0
P P 2
where d is the height of a wave crest over and above the undisturbed stream height (see
figure). Moreover,
1 1
Z Z
z nx dl = f f ′ dx = [f 2 ]X 2 2
−∞ = [(H + d) − H ]
T T 2 2
39
So,
1 1
Z Z Z
2 2
−Dx = ρ u dl − ρ u dl − p nx dl + ρg[(H + d)2 − H 2 ] − ρg[(H + d)2 − H 2 ]
B F F,B 2 2
Z Z
2
= (p + ρu ) dl − (p + ρu2 ) dl
B F
which in general is non-zero. Longuet Higgins shows how this expression may be written in
terms of the potential energy of the wave train at infinity downstream (see Vanden-Broeck,
JFM 1980, V.96, p.610).
Hence the drag on an object is non-zero if on one side there is a train of waves extending to
infinity. Under some circumstances, one can obtain waves trapped over the obstacle with
the free surface becoming flat at the same height on either side like this:
U U
In this case, since the conditions are the same a long way upstream or downstream, we get
Z Z
2
Dx = (p + ρu ) dl − (p + ρu2 ) dl = 0
B F
and so there is no drag on the object (see Forbes, 1982, J.Eng.Math V 16, p.171-180).
where C is the bubble contour and s is arc length around the bubble contour. Using the
Laplace-Young equation, this becomes
Z Z
− pB n ds − γ κn ds.
C C
40
The first integral is easily shown to be zero by the divergence theorem. For the second
integral, we note that, by the definition of curvature,
dt
κn = − ,
ds
where t is the unit tangent to the bubble contour pointing in the direction of increasing arc
length. The force becomes
dt
Z Z
−γ ds = −γ dt = [t]C = 0
C ds C
since the bubble is closed. Hence the bubble experiences no net force in any direction.
A spherical bubble placed within an inviscid flow also experiences zero no force whatever
the far-field flow (i.e. we do not need a free stream at infinity as in section 3.8). The force
on the bubble is Z
D=− pn dS,
B
where B is the closed surface of the bubble. By the Laplace-Young equation given above,
we have
p = pB + γκ,
where pB is the constant pressure inside the bubble, and κ is the bubble curvature, and γ
is the surface tension. So,
Z Z Z
pn dS = pB n dS + γ κn dS
B B B
Since the bubble is spherical, then κ is a constant and the second integral also vanishes, with
the result that the bubble experiences no drag regardless of the nature of the surrounding
flow.
In fact, a closed three-dimensional bubble of any shape experiences zero net force in an
inviscid fluid since it can be shown that
Z
κn dS = 0
S
where κ is the mean curvature at a point on the surface S, for any choice of S (see Blackmore
& Ting, SIAM Review, 27(4), p.569-572). This is consistent with physical thinking: since
the bubble interface has zero inertia, the resultant force acting upon it must be zero (see
also the force balance on an interfacial segment in section 3.4.5).
41
So the total force on the bubble is zero. This means that we cannot have a bubble of fixed
shape in a fluid with gravity, since then the total force on the bubble is
Z Z
− p+ n dS + p− n dS
S S
where p+ and p− are the pressures on the inside of the bubble respectively. In a static field,
p = −ρgy + B for constant B. So, on using the divergence theorem, the total force is
Z Z
+
− ∇p dV + ∇p− dV
V V
Z Z
= ρ+ g dV + ρ− g dV = (ρ+ − ρ− )gV.
V V
But this must be zero by the previous argument. So there must be another force to balance
this buoyancy force. This comes, for example, from the visous drag on the bubble, in which
case the total force acting is a combination of the pressure force and the normal viscous
stress, which must sum to zero to give a total overall force of zero on the bubble.
Summary
We have derived equations governing
Vector identities
The following vector identities are useful in fluid mechanics:
∇2 u = ∇(∇ · u) − ∇ × (∇ × u),
u · ∇u = ∇( 12 |u|2 ) + ω × u, where ω = ∇ × u,
∇ · (φu) = φ∇ · u + u · ∇φ,
∇ × (φu) = φ∇ × u + (∇φ) × u.
42
4. Vorticity dynamics and the stream function
∇ · ω = 0. (4.6)
Physical interpretation
43
w1
w2
Consider a thin vortex tube made up of vortex lines pointing in the direction of the vorticity
vector ω 1 as shown in the left of the figure above. Since the tube is thin, the vorticity |ω 1 | is
approximately constant inside the tube. Suppose that a while later the tube has expanded
a little and the vorticity inside is now |ω 2 |. Then the fact that ∇ · ω = 0 means that
|ω 2 | < |ω 1 |, since the “flux of vorticity” through the tube is conserved.
The term Dω/Dt in (4.5) represents the convection of vorticity with the flow.
The ν∇2 ω in (4.5) represents viscous diffusion of vorticity.
The term ω · ∇u in (4.5) represents stretching of the vortex lines.
Example: Snow clearing at telegraph poles
pole
Wind blows in the x direction over flat ground at y = 0 with velocity u = yi. The vorticity
is ω = ∇ × u = −k, so the vortex lines point along the z axis as shown in the diagram. The
flow encounters a telegraph pole at some x. The vortex lines wrap themselves around this
obstacle, becoming very stretched close to the pole. Therefore, since ∇ · ω = 0, thin vortex
tubes experience greatly intensified vorticity at the pole and the rate of local spinning of
fluid particles increases. This mechanism is responsible for the cleared areas you sometimes
see at the foot of telegraph poles when it snows in the winter.
Notes
1. There is no vortex stretching in two-dimensional flows.
44
This is because if u = u(x, y)i + v(x, y)j, then ω = (vx − uy )k = ωk and so
∂
ω · ∇u = w u = 0.
∂z
Burger’s vortex
A famous example of a fluid flow which exhibits all of the features described above is
Burger’s vortex. In terms of cylindrical polar coordinates, this flow is given by
1
w = kz, u = − kr, v = v(r) (4.7)
2
1 d(rv)
ω = ω(r)k, ω(r) = , (4.8)
r dr
for some chosen constant k. Note that none of the components depend on θ. This is an
example of an axisymmetric flow.
Substituting (4.7, 4.8) into the vorticity transport equation (4.4), we obtain
d(ωr 2 ) d2 ω dω
k + 2ν r 2 + = 0. (4.9)
dr dr dr
kΓ −kr2/4ν
ω= e k, (4.10)
4πν
for some p
constant Γ. So most of the vorticity is concentrated in a ‘core’ region of radius of
order O( ν/k).
It is now possible to find v(r) by integrating the last of equations (4.8). We find
Γ 2
v= (1 − e−kr /4ν ).
2πr
The Burger’s vortex described above is a steady and, remarkably, exact solution of the
Navier-Stokes equations. It involves a perfect balance between the three mechanisms of
convection of vorticity, vortex stretching, and viscous diffusion of vorticity.
45
Illustration of Burger’s vortex.
Consider the steady version of the vorticity transport equation for two-dimensional flow,
so that ω = ωk. Recall that in 2-D the term ω · ∇u is identically zero. If the kinematic
viscosity ν is very small, we may write, to a good approximation,
(u · ∇)ω = 0. (4.11)
This means that the vorticity, ω, is constant along streamlines and hence just depends on
the streamfunction ψ, i.e. ω = ω(ψ). Now,
∂ω ∂ω h ∂ψ ∂ψ i
∇ × ω = ∇ × (ωk) = i− j = ω ′ (ψ) i− j = ω ′ (ψ) u. (4.12)
∂y ∂x ∂y ∂x
We can show (see Problem Sheet 2) that, if there exists a closed streamline C in steady
viscous flow with vorticity ω, then the line integral
Z
(∇ × ω) · dx = 0.
C
Inserting (4.12) into this result yields (since ω is constant along the streamline C)
Z Z
′ ′
ω (ψ) u · dx = ω (ψ) u · dx = 0.
C C
So unless by fluke this line integral happens to be zero, we must have ω ′ (ψ) = 0. Therefore
in a region of closed streamlines the vorticity is constant.
46
The equation of conservation of mass (continuity) states that
∇ · u = 0.
u = ∇ × A,
∇ × A′ = ∇ × (A + ∇q) = ∇ × A + ∇ × ∇q = ∇ × A.
ω = −∇2 A. (4.13)
Examples
i) Two-dimensional flow
In this case, u = u(x, y) i + v(x, y) j.
If we pick A = ψ(x, y)k, then
∂ψ
∇ · (ψk) = = 0.
∂z
Thus we can represent the flow velocity as
∂ψ ∂ψ
u = ∇ × (ψk) = i− j = u i + v j.
∂y ∂x
So
∂ψ ∂ψ
u= , v=− ,
∂y ∂x
and the vorticity is given by (4.13).
47
ii) Axisymmetric flow
a) In cylindrical polar coordinates, the velocity is
u = u(r, z) r̂ + w(r, z) k
So
1 ∂ψ 1 ∂ψ
w= , u=− ,
r ∂r r ∂z
with the vorticity given by (4.13).
and 1 ∂ ∂ψ
u = ∇ × (ψ θ̂) = (rψ) k − r̂ = w k + u r̂.
r ∂r ∂z
So
1 ∂ ∂ψ
w= (rψ), u=− .
r ∂r ∂z
The representation of a flow in terms of a stream function can therefore be done in different,
but equivalent, ways.
48
then ψ(r, θ) 1 ∂ψ
∇·A =∇· φ̂ = = 0.
r sin θ r 2 sin2 θ ∂φ
So
1 ∂ψ 1 ∂ψ
u= , v=− ,
r 2 sin θ ∂θ r sin θ ∂r
with the vorticity given by (4.13).
Physical interpretation
For simplicity, we concentrate on 2-D flow. So u = ∇ × ψk. Now, using the identity
∇ × (φA) = φ∇ × A + ∇φ × A,
we can write
dy
dx
ψ ψ+ d ψ
x
49
Start with the vorticity transport equation
∂ω
+ u · ∇ω − ω · ∇u = ν∇2 ω. (4.14)
∂t
Recall that in 2-D, there is no vortex stretching, and so ω · ∇u = 0. Therefore the left hand
side of (4.14) is
∂ ∂ψ ∂ ∂ψ ∂ n ∂∇2 ψ ∂(ψ, ∇2 ψ) o
(−∇2 ψk) + (−∇2 ψk) − (−∇2 ψk) = − − k,
∂t ∂y ∂x ∂x ∂y ∂t ∂(x, y)
where the Jacobian
∂(ψ, ∇2 ψ) ∂ψ ∂∇2 ψ ∂ψ ∂∇2 ψ
= − .
∂(x, y) ∂x ∂y ∂y ∂x
∂∇2 ψ ∂(ψ, ∇2 ψ)
− = ν∇4 ψ.
∂t ∂(x, y)
50
Then, we have
ρ(x) 1
Z Z
∇2x0 φ(x0 ) = ∇2x0 dV (x) = ∇2x0 ρ(x) dV (x).
V r V r
But recall that we can express the delta function δ(x̂), where x̂ = x − x0 , in the form
1 2 1
δ(x̂) = − ∇
4π r
where r = |x̂| = |x − x0 |. We can do this since G = 1/r is the Green’s function for the
singularly forced Laplace’s equation,
∇2 G = −4π δ(x̂).
= −4πρ(x0 ),
ω = ∇ × u, ∇·u=0
for as yet unknown velocity field u. Let’s express u as the curl of its vector potential as
follows,
u=∇×B
where we have chosen ∇ · B = 0. Then,
∇ × u = ∇ × (∇ × B) = −∇2 B = ω.
So,
∇2 B = −ω.
From the above, a solution to this Poisson equation is
1
Z
ω(x)
B(x0 ) = dV (x).
4π V r
51
We need to check ∇ · B = 0. So,
1 1
Z Z
ω(x) ω(x)
∇ x0 · B(x0 ) = ∇x · dV (x) = ∇ x0 · dV (x)
4π 0 V r 4π V r
1 1
Z
= ω(x) · ∇x0 dV (x)
4π V r
1
Z ω(x)
= ∇x · dV (x).
4π V r
Using the divergence theorem, we have
1 1
Z
∇ x0 · B(x0 ) = ω(x) · n(x) dS(x),
4π S r
1
Z
ω(x)
u= ∇ x0 × dV (x).
4π V r
Exact solutions to the Navier-Stokes equations which may be written down in closed form
are very rare. Most problems are far too complicated to expect to be able to write down
the solution in a simple way. Instead approximate solutions to the equations have to be
found numerically on a computer.
However, there are a number of very well-known examples of exact solutions which can
be found fairly easily. Often these are obtained by reducing the full equations using some
symmetry or other special property of the problem. We begin with a discussion of flows
which only move in one direction.
Keep in mind that we are always trying to find solutions to the Navier-Stokes equations,
∂u
∇ · u = 0, + u · ∇u = −∇p/ρ + ν∇2 u, (5.1)
∂t
with boundary conditions appropriate to the particular problem in hand.
Unidirectional flows
52
These are very special types of flow for which there is only one component of velocity. We
can write
u = wk
so that the fluid is moving with speed w in the z direction only.
The continuity equation gives us some useful information about w straightaway. In Carte-
sians,
∂u ∂v ∂w
∇·u = + + = 0,
∂x ∂y ∂z
but for this unidirectional flow, u = v = 0. So we have
∂w
=0
∂z
which states that w does not change as we move along in the z direction. Hence w can only
be a function of x and y, i.e. w(x, y).
This observation leads to the following remarkable property of unidirectional flows. Consider
the nonlinear term in the Navier-Stokes equations (5.1),
u · ∇u = (wk · ∇)(wk)
∂
= w (wk)
∂z
= 0
since neither w nor k depend on z. So the nonlinear term vanishes! This leaves us with a
much simpler, linear problem to solve.
To summarize, for a unidirectional flow, we have to solve
∂u
u = wk, w = w(x, y), = −∇p/ρ + ν∇2 u, (5.2)
∂t
together with appropriate boundary conditions.
In component form,
0 = −px + µ(uxx + uyy ),
0 = −py + µ(vxx + vyy ),
0 = −pz + µ(wxx + wyy ).
53
Since u = v = 0, the first two of these equations tells us that
px = py = 0.
−d
Suppose that fluid is driven along the channel by a pressure gradient of size −G, where
G > 0, i.e.
dp
= −G.
dz
We try looking for a unidirectional solution with w = w(y) in the z direction and v = 0 in
the y direction. Then, from (5.3), we need to solve
d2 w
0=G+µ , w(d) = w(−d) = 0. (5.4)
dy 2
The boundary conditions are those of no-slip at the walls y = d and y = −d. Note that the
normal flow condition v = 0 at y = ±d is satisfied automatically since v ≡ 0 everywhere.
Integrating (5.4) twice with respect to y, we find
1 2
Ay + B = Gy + µw
2
for arbitrary constants A, B. These are set by applying the boundary conditions:
1 2
w(d) = 0 =⇒ Ad + B = Gd ,
2
1 2
w(−d) = 0 =⇒ −Ad + B = Gd
2
Solving for A and B we obtain the solution.
G 2
w= (d − y 2 ). (5.5)
2µ
54
y
It is customary to plot the velocity profile with the velocity on the horizontal axis like this:
The maximum velocity wmax = Gd2 /2µ occurs on the centre line y = 0.
Flux
One feature of the flow we might be interested in is the volume flow rate or flux through
the channel. This follows from integrating the velocity profile across the channel width. So
the flux
Z d
G d 2
Z
Q= w dy = (d − y 2 ) dy
−d 2µ −d
2Gd3
= .
3µ
This represents the volume of fluid passing through a vertical line across the channel in unit
time.
0 z
In this case the fluid flow is driven purely by the upper wall motion. We expect a unidirec-
tional solution of the form w = w(y). In that case, equations (5.3) imply
d2 w
0= .
dy 2
55
So
w = U y/d,
and the velocity profile looks like this:
y
d
w
U
r=a
dp
= −G, G>0
dz
as the pressure gradient driving the fluid flow.
The continuity equation states that
∂u u 1 ∂v ∂w
∇·u = + + + = 0.
∂r r r ∂θ ∂z
But, since w = wk and u = v = 0,
∂w
= 0.
∂z
Therefore we just have w = w(r, θ). Assuming no variation in θ we take w = w(r).
So now (5.3) requires us to solve
d2 w 1 dw
0=G+µ + , w(a) = 0. (5.6)
dr r dr
56
r dA=rdrdθ
The boundary condition at the wall is that of no-slip. We will need to apply another
condition at r = 0 later on.
We can re-write (5.6) as
µ d dw
0=G+ r . (5.7)
r dr dr
Integrating once with respect to r,
1 2 dw
A= Gr + µ r ,
2 dr
for constant A. And again,
1 2
A log r + B = Gr + µw.
4
To prevent w from diverging at r = 0, we need to set A = 0. This is the second boundary
condition alluded to above. The constant B is fixed by requiring w(a) = 0, so
1
B = Ga2 .
4
The maximum velocity wmax = Ga2 /4µ occurs on the centre line r = 0.
Flux
To compute the flux through the circular tube, we once more integrate the velocity over
the cross-section. Recall that an element of area in the cross-section is dA = rdrdθ. So,
Z 2π Z a
G a
Z
Q= w rdrdθ = 2π r(a2 − r 2 ) dr
0 0 4µ 0
πGa4
= .
8µ
57
Furthermore, we define the mean velocity,
Flux Q Ga2
umean = = = .
x-sec area πa2 8µ
Shear stress
To calculate the axial shear stress at the wall, recall that Fi = σij nj . In this case the unit
normal at the wall pointing into the fluid is n = −r̂. So the stress in the z direction at the
wall is
Fz = −σzr = −2µezr .
We know that erz = ezr and we are given that
1 ∂w ∂u
erz = +
2 ∂r ∂z
Gr
(in this case) = − .
4µ
Therefore, at the wall,
Ga
Fz = −2µerz |r=a = , (5.9)
r=a 2
which is the shear stress exerted by the fluid on the wall in the downstream (i.e. z) direction.
(x − a)2 + y 2 = a2 ,
58
where the first factor in brackets may be identified as the solution for Poiseuille flow in a
circular pipe. Note that the first factor is identically zero on the smaller circle of radius b
and the second factor is identically zero on the larger circle of radius a. In this sense the
solution resembles that for undirectional flow in an equilateral triangular pipe where each
factor is zero of one the three walls.
A contour plot showing lines of constant velocity u is shown in the figure below for the
case a = 1 and b = 0.3.
0.8
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
Notes
1. CPF is often used by medical workers as a model of blood flow in arteries. The shear stress
is important because of its link with the onset of atherosclerosis, the disease responsible for
about a third of deaths in Western society.
2. The result (5.9) shows that, for a fixed flux through the tube, the driving pressure
gradient G ∝ a−4 . Applying CPF to blood flow, this states that build-up of fatty deposits
on the arterial walls due to an unhealthy lifestyle rapidly increases the amount of effort
required of the heart to drive the same amount of blood through the body. If the artery’s
diameter is halved, for example, the pressure gradient needs to be 16 times larger to produce
the same flux!
3. Experimentally, one might set up CPF as shown: A circular tube of length L has one
p p
1
2
CPF
Reservoir
59
end in a reservoir of water so that the entry pressure p2 is held fixed. The pressure at the
exit, p1 , is also constant. The pressure gradient is then
p − p ∆p
2 1
G=− =− , say.
L L
This experiment was first performed by Osborne Reynolds in 1883. His original apparatus
is kept in the Engineering department at Manchester University. Reynolds found that for
sufficiently low flow rates Q, the CPF solution (5.8) was indeed observed. He injected a
blob of dye at the entrance to the pipe and it stretched out in a perfectly straight line, as
this solution predicts.
streak of dye
However, as Q was increased, Reynolds found that at some distance along the pipe the
laminar CPF flow eventually broke up and thereafter became turbulent.
As Q was increased, the start of the turbulent region moved further and further upstream.
The reason for the break-up of the fluid flow in a circular pipe into turbulence is still not
fully understood!
f(x,y)=0
With a pressure gradient −G in the z direction driving the flow, the problem to be solved
is
G ∂2 ∂2
∇2 u = − , where ∇2 = +
µ ∂x2 ∂y 2
60
with
u=0 on f (x, y) = 0.
For very special cross-sectional geometries (a circle, an ellipse, a triangle – see Sheet 3) the
solution can be found in closed form. Otherwise a solution can be sought using a complex
variable - conformal mapping technique. Alternatively, we can seek a series solution, as in
the next example.
2b x
2a z
−a ≤ x ≤ a, −b ≤ y ≤ b,
If we look for a separable solution of type u(x, y) = X(x)Y (y), we find that the velocity
can be expressed in the series form
∞ X
∞
X h πx i h πy i
u= Amn cos (2m + 1) cos (2n + 1) , (5.11)
m=0 n=0
2a 2b
61
Exercise 2: Show that the channel solution above reduces to plane Poiseuille flow (PPF) in
the limit b/a → 0.
y g
α
h
x
Fix Cartesian axes as shown. The body force, F, due to gravity is given by
From before, we know that for unidirectional flow u · ∇u = 0. The momentum equations
are therefore,
∂p d2 u
0 = − + ρg cos α + µ 2 ,
∂x dy
∂p
0 = − − ρg sin α.
∂y
62
Therefore, the first momentum equation becomes
d2 u
0 = ρg cos α + µ .
dy 2
Integrating twice,
ρgy 2
u=− cos α + βy + γ,
2µ
for constants β, γ.
No-slip at y = 0 =⇒ γ = 0.
At y = h there is a fluid-fluid interface between the water (say) in the film and the air above
it. The stress condition at y = h is given by (assuming negligible surface tension)
(air) (water)
[σij − σij ] nj = 0,
where
n = (0, 1, 0) = j.
So,
(air)
−pa ni + 2µ(air) eij |y=h nj = −p|y=h ni + 2µ eij |y=h nj ,
i.e.
p|y=h = pa ,
since e22 = ∂v/∂y = 0, and
(air)
2µ(air) e12 |y=h = 2µ e12 |y=h .
µ(air) ≪ µ,
Exercise: Follow a similar argument to above to derive the shear stress condition at the
surface of a liquid film flowing down the outside of a cylindrical rod of radius a.
Applying the condition of zero surface stress, we find
ρgh cos α
0=− + β,
µ
63
and so
ρgy
u= (2h − y) cos α,
2µ
which can be re-written
ρg h 2 i
u= h − (h − y)2 cos α.
2µ
This is the Nusselt solution. It represents ‘half’ of plane Poiseuille flow.
1 Poiseuille flow
2
In other words, given a flux of fluid F and slope angle α, a steady unidirectional motion is
possible for a film of height
3µF 1
3
h= .
ρg cos α
What happens if F is suddenly increased by a small amount dF ?
Since dF is small, we expect a corresponding increase in the layer thickness,
dh
dh = dF .
dF
The jump travels downstream at a finite speed, V , so that, with axes fixed in the jump, we
have the picture:
Continuity requires that
64
h+
dh y
jump
F+dF
h
F
x
h+ V
dh
stationary
F+dF V
h
F
x
i.e.
(F + dF ) − (h + dh)V = F − hV.
Thus, the speed of the jump is twice the maximum speed of the fluid.
∂u
u = wk, w = w(x, y, t), = −∇p/ρ + ν∇2 u.
∂t
As in the steady case, the first and second momentum equations imply that
px = py = 0 =⇒ p = p(z, t).
65
1. Impulsively started flows - Rayleigh layers
a) Unbounded flow
Consider a semi-infinite mass of fluid at rest above a plane, rigid wall.
y
Fluid
νWyy − Wt = 0,
y = αŷ, t = α2 t̂
leaves the equations and boundary conditions unchanged. This motivates us to seek a
solution with the same property. We note that the grouping
y
t1/2
is unchanged by the α-transformation. Accordingly, we try looking for a solution in the
form
y
W = f (η), where η = .
2(νt)1/2
66
This is called a similarity solution, for reasons which will become clear below. It has the
advantage of reducing the partial differential equation to an ordinary one. Applying the
chain rule,
∂ 1η ∂ ∂ 1 ∂ ∂2 1 ∂2
=− , = , = .
∂t 2 t ∂η ∂y 2(νt)1/2 ∂η ∂y 2 4νt ∂η 2
It is known that ∞ √
π
Z
2
e−ξ dξ = .
0 2
Thus η
2
Z
2
f =1− √ e−ξ dξ = 1 − erf(η) = erfc(η),
π 0
where erf is the Error function. Plotting y against f :
y
O( νt 1/2)
t increasing
1
f
Thus there is a layer of thickness O((νt)1/2 ), in which the fluid is most affected by the wall
motion. This is called a Rayleigh layer. The (y, f ) profiles shown in the figure above are
self-similar. This means that they are just re-scaled versions of each other. For example,
67
we could multiply the first profile by a suitably-chosen and get it to lie directly on top of
any of the other profiles.
Considering the vorticity of this impulsively-driven flow,
∂w U dW U 1 2
ω = ∇ × u = ∇ × (wk) = i= 1/2
i=− 1/2
√ e−η i.
∂y 2(νt) dη (νt) π
So,
U 1 2
ω=− 1/2
√ e−y /(4νt) i.
(νt) π
Now, if y 2 ≫ 4νt, the exponential is very small indeed and ω is approximately zero. How-
ever, if y 2 ∼ 4νt, the exponential is not small and ω will have an appreciable value. So, as
time progresses, measurable values of vorticity lie in the approximate y range [0, 2(νt)1/2 ].
This range extends in time like the square root of t. Thus the vorticity (shown as a shaded
area) diffuses away from the wall in time.
y y
O( t 1/2)
O( t 1/2)
z z
later time
Thus, the vorticity behaves like heat diffusing away from an impulsively heated plate at
y = 0.
Similarity solutions
We can now see why we called the previous solution a similarity solution. As time changes,
the velocity profile is simply rescaled, but retains its original shape. Any profile at any
given time may be simply rescaled to lie on top of any other profile.
A similarity solution is possible in a problem like this when there is no natural lengthscale. If
the plate were finite we could choose its length, L, say, as a basis for non-dimensionalization.
In this case, the plate is infinite and so we cannot do this. Consider the governing equation:
νWyy = Wt .
We know that ν has dimensions [ν] = L2 T −1 . Thus the dimensions of the LHS are
[ν] L2 1 1
2
= 2
=
[y ] T L T
and those of the RHS are
1 1
= .
[t] T
68
A dimensional balance therefore gives
[ν][t]
= 1.
[y 2 ]
So the only possible dimensionless grouping of the variables is
y
(νt)1/2
or any power thereof. So the solution must feature y, t and ν grouped only in this way.
We can conclude that the fluid velocity must be expressed in the form
y
w = W f (η), η= ,
(νt)1/2
a similarity solution with similarity variable η.
z
W0 cos(ω t)
Again, it is reasonable to assume that the pressure gradient G = 0 and that the flow is
unidirectional, with w = w(y, t). So continuity is satisfied and the problem is
νwyy = wt ,
with
w → 0 as y → ∞, w = W0 cos(ωt) on y = 0.
Since the frequency of the wall motion is ω and the problem is linear, we expect the frequency
of the fluid velocity w also to be ω.
Thus we seek a solution by writing
n o
w = W0 ℜ W (y)eiωt ,
69
W =1 on y = 0, W →0 as y → ∞.
So
W = Aeσy + Be−σy ,
for constants A, B and where r
iω
σ= .
ν
√ √
Now, i = (1 + i)/ 2 and so, to satisfy the infinity condition, we need to set A = 0. Thus,
h i
W = B exp − (1 + i)y/(2ν/ω)1/2 ,
and so
w √ω r
ω
= e− 2ν y cos ωt − y .
W0 2ν
ω −√ ω y
r r
∂w ω π
ω = ∇×u = i = −W0 e 2ν cos ωt − y+ . (5.12)
∂y ν 2ν 4
According to the exponentially
p decaying term in (5.12), most of the vorticity is confined
in a region of thickness O( ν/ω). This region is called a Stokes layer. The higher the
frequency, the thinner the Stokes layer. Note that the Stokes layer thickness does not grow
in time. This contrasts with the growing Rayleigh layer studied above.
y
Stokes layer
1/2
O( (ν/ω) ) z
Therefore there is a phase difference of π/4 between the velocity w and the wall shear.
70
y
111111111111111111111111111111
000000000000000000000000000000 x
∇2 ψ = 0.
v = 0 on y = 0 or ψ = 0 on y = 0.
ψ = kxy,
Clearly this solution does not satisfy the no slip condition at y = 0. To manage to do this,
we need to incorporate viscous effects.
By analogy with the inviscid solution, we seek a solution to the viscous flow problem of a
similar form. Specifically, we write
u = xf ′ (y), v = −f (y),
where f is a function to be determined. With this choice, we note that the continuity
equation is satisfied straightaway, since
ux + vy = f ′ − f ′ = 0.
71
To determine f , we need to solve the x-momentum equation,
1
uux + vuy = − px + ν[uxx + uyy ].
ρ
However, we don’t yet know the form of the pressure. To resolve this issue, consider the
y-momentum equation,
1
uvx + vvy = − py + ν[vxx + vyy ],
ρ
which becomes,
py = −ρ[f f ′ + f ′′ ] ≡ Q(y), say.
Z y
=⇒ p = p(x) + ρ Q(ξ) dξ,
0
for function of integration p(x).
Rearranging the x momentum equation, we find
px = ρx[νf ′′′ + f f ′′ − f ′2 ].
νf ′′′ + f f ′′ − f ′2 = β, (5.13)
β = −k2 .
νf ′′′ + f f ′′ + k2 − f ′2 = 0,
f = f ′ (0) = 0 on y = 0; f′ → k as y → ∞.
The first two boundary conditions are those of no normal flow and no slip at the wall.
72
70
60
50
40
y
30
20
10
0
-1 -0.5 0 0.5 1
x
Notes
1. The solution for f must be found numerically.
2. The streamlines for this flow look like this:
3. The stagnation point flow on a plane wall is also relevant to that occuring on a curved,
bluff body.
In the neighbourhood of the stagnation point, the body may be considered locally flat. For
small x, the external inviscid flow may be expanded in a power series
U (x) = u0 x + u1 x2 + u2 x3 + · · · ,
and so locally U ∝ x.
73
So, for an axisymmetric flow, the velocity field is given by
u = u(r, z, t) r̂ + w(r, z, t) k.
Ω1 Ω2
u = v(r) θ̂.
v2 1 ∂p d2 v 1 dv v
− =− , 0= 2 + r dr − r 2 .
r ρ ∂r dr
74
Solving the second of these,
B
v = Ar + ,
r
for constants A, B.
Thus the solution is a combination of a rigid-body rotation
u = r θ̂,
Thus,
b2 Ω2 − a2 Ω1 Ω2 − Ω1
A= , B = a2 b2 .
b2 − a2 a2 − b2
The vorticity
b2 Ω2 − a2 Ω1
ω = ∇×u = 2 k,
b2 − a2
and is therefore constant.
What is the couple on the inner cylinder? [Recall that couple (or torque or moment) means
the turning effect of a force acting at a distance from a given point.]
To determine this, we will need σθr at r = a. We know that (see sheet),
1 d v B
eθr = r = − 2.
2 dr r r
Thus,
2µB
σθr = − .
r2
Hence, the couple or torque per unit length in the z direction is given by
Z 2π
τa = (aσθr |r=a ) a dθ = −4πµB.
0
So,
Ω2 − Ω1
τa = 4πµa2 b2 .
b2 − a2
75
Similarly, the couple on the outer cylinder is found to be
Ω2 − Ω1
τb = −4πµa2 b2 .
b2 − a2
Special cases:
but still everything is independent of θ. We continue to assume the flow is also independent
of z.
76
The continuity equation implies,
1 d
(ru) = 0,
r dr
and so
ru = constant.
Suppose that u = −U at r = a, with U > 0, so there is suction of fluid at the inner cylinder.
Then
aU
u=− ,
r
and, at r = b,
aU
u=− ,
b
so there is injection of fluid at the outer cylinder.
The r-momentum equation states that
du v 2 1 dp d2 u 1 du u
u − =− +ν + − ,
dr r ρ dr dr 2 r dr r2
while the θ-momentum equation has
dv uv d2 v 1 dv v
u + =ν + − .
dr r dr 2 r dr r 2
d2 v U a 1 dv Ua v
+ 1 + − 1 − = 0.
dr2 ν r dr ν r2
This is an equation for v(r). Since u is already known, once v is determined, we can in
principle solve the r-momentum equation to obtain the pressure p.
Defining a Reynolds number based on the suction velocity, R = U a/ν, we can re-write the
v equation as
d2 v 1 dv v
2 + (1 + R) r dr − (1 − R) r 2 = 0.
dr
Seeking a solution of the form v = kr n , we get
and so
n = −1, (1 − R).
In the special case of R = 2, n has a repeated root of −1.
77
So there are two possibilities:
C
v = + Dr (1−R) if R 6= 2,
r
C ln r
v = +D if R = 2.
r r
v = aΩ1 on r = a,
v = bΩ2 on r = b.
Then,
a2 Ω1 h b2−R
1−R
i
v = − r
b2−R − a2−R r
a2 Ω1 h 1 − (r/b)2−R i
=
r 1 − (a/b)2−R
a2 Ω1 h 1 − (b/r)R−2 i
= .
r 1 − (b/a)R−2
Consider what happens when the Reynolds number becomes large, i.e. R ≫ 1. Then,
a2 Ω1 a R−2 a R−1
v≈ = aΩ1
r r r
Since r > a and R is very large, this is essentially zero unless r is very close to a. In other
words, the only significant swirl occurs very near to the inner cylinder.
The figure shows that the solution goes from near zero to one over a very short distance.
This becomes more severe as R gets larger.
This is our first example of a boundary layer, a very narrow region over which there is a
very rapid change in the solution. We will deal with these in more detail later in the course.
78
a=1, b=2, Omega1=1
1
R=10
R=100
0.8
0.6
v(r)
0.4
0.2
0
1 1.5 2
r
Make a cup of coffee and stir it round a few times. How long does it take for the coffee to
come to rest?
We consider a simplified version of this problem, that of how long it takes swirling fluid in
an infinite cylinder to stop once the forcing is removed. Suppose viscous fluid occupies the
cylinder r = a and suppose that both the cylinder and the fluid are initially rotating with
angular velocity Ω, so
u = v θ̂, v = Ω r for t < 0.
This is the special case discussed above of rigid body rotation with no inner cylinder.
At t = 0 the cylinder is suddenly stopped. The ensuing unsteady fluid motion is governed
by the θ-momentum equation
∂v ∂ 2 v 1 ∂v v
=ν + − ,
∂t ∂r 2 r ∂r r 2
with
v = 0 on r = a for t > 0.
79
for separation constant λ. Thus,
Ġ = −νλ2 G,
which has solution
2
G = e−νλ t .
This leaves
1 1
H ′′ + H ′ + λ2 − 2 H = 0,
r r
or
r 2 H ′′ + rH ′ + (λ2 r 2 − 1)H = 0.
This is Bessel’s equation of order 1 (see MTH-2C23). It has the general solution
H(a) = 0 =⇒ J1 (λa) = 0,
0.6
0.4
0.2
0 2 4 6 8 10 12 14 16 18 20
x
–0.2
80
Defining λn = βn /a, the general solution is given by
∞
2 2
X
v= An J1 (βn r/a) e−νβn t/a ,
n=1
for constants An .
We can exploit the orthogonality of the Bessel functions to find the An . We know that
Z a
rJ1 (λn r)J1 (λm r) dr = 0 if n =6 m.
0
The initial condition is that v = Ωr at t = 0 and so, using the general solution,
∞
X
Ωr = An J1 (λn r).
n=1
Multiplying both sides by rJ1 (λn ), integrating from 0 to a, and using the orthogonality
property above, we find
Z a Z a
2
Ω r J1 (λn r) dr = An r[J1 (λn r)]2 dr.
0 0
Thus, Ra
Ω 0 r 2 J1 (λn r) dr
An = R a 2
.
0 r[J1 (λn r)] dr
t = te = a2 /(νβ12 ).
81
So by this time, the swirling fluid has significantly slowed.
Going back to the cup of coffee, if we take a = 4cm and ν = 10−6 m2 s−1 (the value for
water), we get
te ≈ 2 minutes,
which seems rather long. Experiment suggests that in fact the coffee slows to 1/e of its initial
value in about 15 seconds. The discrepancy lies in the crucial role played by a boundary
layer (see later lectures) of fluid on the bottom of the cup (see Acheson pp. 45, 165, 284).
We have shown that the steady flow of a Newtonian, incompressible fluid of kinematic
viscosity ν and density ρ satisfies the Navier-Stokes equations (written here without a body
force),
Stokes flow, sometimes referred to as creeping flow, describes the motion of an extremely
viscous or sticky fluid such as treacle or golden syrup. Under these conditions, we might
expect
|u · ∇u| ≪ |∇2 u|,
since the flow is very sluggish and so the magnitude of the velocity at any point is small.
Therefore, since the left hand side is quadratic in the small velocity, we would expect it to
be smaller than the term on the right hand side as a reasonable approximation.
We may derive these equations more formally as follows. Suppose that characteristic length
and velocity scales are L and U respectively. For example, we may wish to model a bacterium
swimming at speed U , where the length of the bacterium is L. Then, we nondimensionalize
(i.e. remove all dimensions from the problem) by writing
µU ∗
u = U u∗ , v = U v∗ , p= p , x = L x∗ , y = L y∗.
L
82
Then all quantities with an asterisk superscript are dimensionless. In component form, the
momentum and continuity equations for a 2D steady flow with no body force acting are:
1
uux + vuy = − px + ν(uxx + uyy ),
ρ
1
uvx + vvy = − py + ν(vxx + vyy ),
ρ
ux + vy = 0.
Substituting in, we obtain (using the chain rule to convert the derivatives):
U2 ∗ ∗ νU νU ∗
u ux∗ + v ∗ u∗y∗ = − 2 p∗x∗ + (u ∗ ∗ + u∗y∗ y∗ ),
L L L2 x x
U2 ∗ ∗ νU νU ∗
u vx∗ + v ∗ vy∗∗ = − 2 p∗y∗ + (v ∗ ∗ + vy∗∗ y∗ ),
L L L2 x x
u∗x∗ + vy∗∗ = 0.
Formally taking the limit R → 0, we derive the Stokes equations of creeping fluid motion:
∇ · u = 0, 0 = −∇p + µ∇2 u.
83
The key step is taking the limit R → 0, which means we are considering flows for which the
Reynolds number is small enough for the nonlinear term u · ∇u to be neglected. By the
definition of R,
UL
R= .
ν
This means that for R ≪ 1, we need U to be small, or L to be small, or for ν to be large.
So Stokes flow arises if
• The characteristic velocity scale is sufficiently small (e.g., a bubble moving through
treacle)
If the scales U , L, and ν combine to make R small, we may use the Stokes equations (6.2)
to approximate the fluid flow. The great advantage of doing this lies in the fact that the
Stokes equations are linear. In general, linear problems are easier to solve than nonlinear
ones.
∇2 u = ∇(∇ · u) − ∇ × (∇ × u),
∇ · u = 0, ∇p = −µ∇ × ω (6.4)
0 = ∇ × (∇ × ω),
and hence we derive the equations of vorticity transport for a Stokes flow,
∇ · ω = 0, 0 = ∇2 ω. (6.5)
Instead, taking the divergence of the second of equations (6.4), we see that, for a Stokes
flow, p satisfies Laplace’s equation,
∇2 p = 0, (6.6)
84
1111111111111111111111111111111111
0000000000000000000000000000000000
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
xo
0000000000000000000000000000000000
1111111111111111111111111111111111
1 V
0
Ω
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
0000000000000000000000000000000000
1111111111111111111111111111111111
∇2 φ = 0
[Sketch region D]
Green’s third identity in two-dimensions states that, at a point x0 in D,
1 1 1
Z Z
φ(x0 ) = ln r ∇φ(x) · n(x) dl(x) + (x0 − x) · n(x) φ(x) dl(x), (6.7)
2π C 2π C r 2
where n(x) is the unit normal at a point x on the curve C, and dl(x) is the increment of
arc length along C at the point x.
Choose C to be a circle of radius a centred about the point x0 . Then (6.7) gives
ln a 1
Z Z
φ(x0 ) = ∇φ(x) · n(x) dl(x) + (x0 − x) · n(x) φ(x) dl(x).
2π C 2πa2 C
85
Hence we are left with
1
Z
φ(x0 ) = φ(x) dl(x). (6.8)
2πa C
This is Gauss’ mean value theorem for harmonic functions in two-dimensions. It states that
φ at a point x0 is the average of the φ values over a surrounding circle of radius a.
Clearly ∇p · n cannot be single-signed over the arbitrary surface S. It follows that p must
be increasing or decreasing arbitrarily close to the point x0 , contradicting the supposition
that it is either a maximum or a minimum. Therefore the points of maximum or minimum
pressure in a closed-region Stokes flow can only occur on the boundary. ⋄
Nifty move: This observation hints at a more profound characteristic of Stokes flows. In
fact it is possible to describe a Stokes flow entirely by the values the various flow variables
assume on the boundary. This leads to a very efficient means of computing Stokes flows,
known as the boundary integral method. For more information on this advanced procedure,
see Pozrikidis, ‘Boundary Integral and Singularity Methods’, Cambridge (1992).
If the flow in question is two-dimensional with velocity field u = ui + vj, we may introduce
a stream function defined by
∂ψ ∂ψ
u= , v=− .
∂y ∂x
86
equation
∇4 ψ = 0. (6.10)
We will make use of this feature of two-dimensional Stokes flows later when we study flow
in corners.
∇ · u = 0, µ∇2 u = ∇p
and
∇ · u′ = 0, µ∇2 u′ = ∇p′
with u = u′ on the boundary S.
Let eij , e′ij be the stress tensors associated with (u, p) and (u′ , p′ ) respectively.
Then,
1 ∂ ∂
Z Z h i
(eij −e′ij )eij
dV = (ui − u′i ) + (uj − u′j ) eij dV
V 2 V ∂xj ∂xi
1h ∂ 1h ∂
Z i i
(since eij symmetric) = (ui − u′i ) eij + (ui − u′i ) eij dV
V 2 ∂xj 2 ∂xj
∂
Z
= eij (ui − u′i ) dV
V ∂x j
∂ ∂eij
Z Z
= [(ui − u′i )eij ] dV − (ui − u′i ) dV
V ∂xj V ∂xj
∂eij
Z Z
(by divergence theorem) = (ui − u′i )eij nj dS − (ui − u′i ) dV
S V ∂xj
∂eij
Z
(since ui = ui on S) = − (ui − u′i )
′ dV .
V ∂xj
87
1 h ∂ 2 ui ∂ ∂uj i
Z
= − (ui − u′i ) + dV
V 2 ∂xj ∂xj ∂xi ∂xj
1 ∂ 2 ui
Z
(by continuity) = − (ui − u′i ) dV
V 2 ∂x2j
1 ∂p
Z
(♣) (since µ∇2 u = ∇p) = − (ui − u′i ) dV
V 2µ ∂xi
1 ∂
Z
(by continuity again) = − [p(ui − u′i )] dV
2µ V ∂xi
1
Z
(by divergence theorem) = − p(ui − u′i )ni dS
2µ S
(since ui = u′i on S) = 0.
Hence,
Z
(eij − e′ij )eij dV = 0. (6.11)
V
µ∇2 u′ 6= ∇p′ ,
88
i.e. (u′ , p) does not represent a Stokes flow.
From above, we know that,
Z
(eij − e′ij )eij dV = 0,
V
but now Z
(e′ij − eij )e′ij dV 6= 0 since µ∇2 u′ 6= ∇p′ .
V
Recall that the viscous dissipation for the flow (u′ , p′ ) is defined as
Z
′
Φ = 2µ e′ij e′ij dV .
V
Clearly
Φ′ ≥ 0.
The viscous dissipation for flow (u, p) is
Z Z
Φ = 2µ eij eij dV = 2µ e2ij dV .
V V
Hence, Z
Φ′ − Φ = 2µ (e′ij − eij )2 dV ≥ 0.
V
QED ⋄
Notes
1. Navier-Stokes flows always dissipate more energy than Stokes flows.
2. Since plane (PPF) and circular Poiseuille flow (CPF) happen to satisfy the Stokes
equations, they are therefore minimum dissipation flows.
89
so if we change u to −u, p to −p and F to −F, we get another equally valid solution. If
the boundary conditions involve u we must also switch the sign in these. This property is
not valid for flows at non-zero Reynolds number, as the nonlinear term u · ∇u scuppers this
sign-reversal process.
It follows that if we reverse the velocity on the boundary of a Stokes flow, we also reverse
the flow velocity throughout the fluid.
Puzzle: There is an interesting consequence of this feature of Stokes flows. A solid sphere
moves slowly through a viscous fluid parallel to a plane wall. Does the sphere feel any force
repelling it away from or attracting it towards the wall?
Answer: No; it can’t. It if did, when we reversed the flow it would naturally feel the same
force in the opposite direction. But the symmetry of the problem makes this ridiculous.
Therefore the force acting on the sphere perpendicular to the wall must be zero.
Fix axes in the sphere and adopt spherical polar coordinates (r, θ, φ). The oncoming fluid
moves at speed U and the entire flow is taken to be steady.
For this problem, it is convenient to work with the Stokes equations in the form
∇ · u = 0, ∇ × (∇ × ω) = 0,
since ∇2 ω = ∇ × (∇ × ω).
The fluid velocity
u = u r̂ + v θ̂ + w φ̂.
Exploiting the axisymmetry of the sphere, we can assume that there are no changes in φ,
∂
≡ 0, and also w ≡ 0.
∂φ
We can therefore introduce a stream function ψ (see earlier section on these), defined by
1 ∂ψ 1 ∂ψ
u= , v=− .
r 2 sin θ ∂θ r sin θ ∂r
90
In spherical polars, for vector function F = Fr r̂ + Fθ θ̂ + Fφ φ̂,
1 ∂ 1 ∂ 1 ∂ 1 ∂Fr
∇×F = (Fφ sin θ) r̂ − (rFφ ) θ̂ + (rFθ ) − φ̂.
r sin θ ∂θ r ∂r r ∂r r ∂θ
So the vorticity,
n 1 ∂ 1 ∂ψ 1 ∂ 1 ∂ψ o
ω =∇×u= − − φ̂,
r ∂r sin θ ∂r r ∂θ r 2 sin θ ∂θ
or,
ω = ωφ φ̂,
where
D2 ψ ∂ 2 ψ sin θ ∂ 1 ∂ψ
ωφ = − , D2 ψ ≡ + 2 .
r sin θ ∂r 2 r ∂θ sin θ ∂θ
Now,
1 ∂ 1 ∂
∇×ω = (ωφ sin θ) r̂ − (rωφ ) θ̂
r sin θ ∂θ r ∂r
1 ∂ 1 ∂
= − 2 (D 2 ψ) r̂ + (D 2 ψ) θ̂.
r sin θ ∂θ r sin θ ∂r
Hence
∇ × (∇ × ω) = 0 =⇒ D 2 (D2 ψ) = 0 =⇒ D 4 ψ = 0.
u → U cos θ, v → −U sin θ as r → ∞.
91
These mean that
1 2 2
ψ→ U r sin θ as r → ∞, (6.14)
2
and
∂ψ ∂ψ
= =0 on r = a.
∂r ∂θ
These boundary conditions suggest that we seek a solution of the form
ψ = f (r) sin2 θ.
f = f ′ = 0 on r = a.
In this case,
n ∂2 sin θ ∂ 1 ∂ o
D2 ψ = + f (r) sin2 θ
∂r 2 r 2 ∂θ sin θ ∂θ
f
= f ′′ − 2 2 sin2 θ
r
f
= H(r) sin2 θ, where H(r) = f ′′ − 2 2 .
r
So,
H
D4 ψ = H ′′ − 2 2 sin2 θ
r
h d2 2 2 i 2
= − f sin θ.
dr 2 r 2
= 0.
92
for constants A, B, C, D.
The infinity condition (6.14) implies that
U
B= , D = 0,
2
while
∂ψ A U a2
= 0 =⇒ + + Ca = 0,
∂θ r=a a 2
∂ψ A
= 0 =⇒ − 2 + U a + C = 0,
∂r r=a a
giving
U a3 3
A= , C = − U a.
4 4
Thus,
1 3a a3
ψ = U r 2 sin2 θ 1 − + 3 ,
2 2r 2r
or
1 a 2 a
ψ = U r 2 sin2 θ 1 − 1+ .
2 r 2r
Now we know ψ, we can compute the velocity and vorticity components:
f
r sin θ ωφ = − f ′′ − 2 sin2 θ (= −D 2 ψ)
r
3 Ua
= − sin2 θ
2 r
3U a
=⇒ ωφ = − 2 sin θ.
2r
Also,
1 ∂ψ 3a a3
u = = U cos θ 1 − + ,
r 2 sin θ ∂θ 2r 2r 3
1 ∂ψ 3a a3
v = − = −U sin θ 1 − − 3 .
r sin θ ∂r 4r 4r
To calculate the fluid pressure, we make use of the alternative form of the momentum
equation
∇p = −µ∇ × ω,
n 1 ∂ 1 ∂ o
= −µ (ωφ sin θ) r̂ − (rωφ ) θ̂ .
r sin θ ∂θ r ∂r
93
So,
1 ∂p 1 3U a
= − − 2 sin θ cos θ
µ ∂r r sin θ r
p 3U a
=⇒ = − 2 cos θ + g(θ), arbitrary g,
µ 2r
and,
1 1 ∂p 1 ∂
= (rωφ )
µ r ∂θ r ∂r
3U a
= sin θ
r3
p 3U a
=⇒ = − 2 cos θ + f (r), arbitrary f .
µ 2r
Combining these,
3µU a
p=− cos θ + const.
2r 2
Letting r → ∞, p → p∞ , its value at infinity. So,
3U a
p=− cos θ + p∞ .
2r 2
Notes:
1. The solution for ψ is symmetric about θ = π/2. Therefore the streamlines are symmetric
fore and aft the sphere.
Exercise: Use Maple to plot the streamlines around the sphere.
2. If we impose a velocity −U on the solution, we can obtain results for a sphere moving
through a fluid at rest.
94
1
So Ψ = ψ − U r 2 sin2 θ
2
3 a2
= − U ar sin2 θ 1 − 2 .
4 3r
∂uj o
n ∂u
i
Fi = σij nj = −pni + +
nj
∂xj ∂xi
= −pni − µ(n × ω)i (see Problem Sheet 2).
To get the desired component, we need to take the dot product of F with i, the unit vector
pointing along the line of motion. Denote i as x̂ = (1, 0, 0) in a Cartesian set of axes. Then
θ r
we get
95
Thus, on r = a,
3µU
F · x̂|r=a = −p∞ cos θ + .
2a
Therefore, the total force in the x̂ direction is:
Z 3µU
− p∞ cos θ dS
sphere 2a
Z 2π Z π
2 3µU
= 4πa − p∞ a2 sin θ cos θ dθ dφ (6.15)
2a 0 0
since the surface element dS = a2 sin θdθdφ in spherical polar coordinates. So the total x
force is Z π
2 3µU
4πa − 2π p∞ a2 sin θ cos θ dθ
2a 0
The integral in the second term is zero, and so the total drag on the sphere,
Ftot = 6µπaU.
Taking the Reynolds number to be R = U a/ν, we can define a dimensionless drag coefficient
Ftot 12
CD = 1 2 2
= .
2 ρU πa
R
The agreement with experiment is good for R < 1 and very good for R < 21 .
Consistency analysis
Is the above solution for the sphere everywhere a consistent approximation to the correct
solution of the full Navier-Stokes equations when R ≪ 1? In other words, are the approxi-
mations leading to the Stokes equations valid for all values of r?
To answer this question, we examine the order of magnitude of the terms.
We have
1 ∂ψ 3a a3
u = = U cos θ 1 − + ,
r 2 sin θ ∂θ 2r 2r 3
1 ∂ψ 3a a3
v = − = −U sin θ 1 − − 3 ,
r sin θ ∂r 4r 4r
and so,
∂u U 2a
u ∼O for large r,
∂r r2
∂2u Ua
ν 2 ∼O ν 3 for large r.
∂r r
96
Recall that we neglected the nonlinear terms, u · ∇u, from the Navier-Stokes equations. So,
terms neglected U 2 a/r 2
∼ O
terms retained νU a/r 3
Ur
= O
ν
Rr
= O .
a
i.e., for any fixed Reynolds number, R, no matter how small, we can still find r sufficiently
large such that the terms neglected are as important, or more important, than the retained
viscous terms.
So, although the Stokes solution is consistent near the sphere for 0 < R ≪ 1, it breaks
down as a valid approximation when
a
r=O .
R
This is a major problem, but at least it’s not as bad as the equivalent flow past a circular
cylinder, where the Stokes solution doesn’t work at all (see Problem Sheet 4)!
The difficulty in the far-field was partially tied up by Oseen in 1910, and later more fully
by Proudman & Pearson in 1957.
U
θ
z
Drop
The axes are fixed in the moving drop, so the outer fluid appears to be flowing towards it
at speed U .
For this problem we have to consider the flows inside and outside of the drop. In both cases,
we assume that the conditions of Stokes flow are satisfied. Thus, the equations of motion
are
∇ · u = 0, ∇ × (∇ × ω) = 0.
97
As for the solid sphere calculation above, we assume the flows are axisymmetric,
∂
≡ 0, and also w ≡ 0.
∂φ
Following the results for the solid sphere, the equation of motion reduces to
∂ 2 ψ sin θ ∂ 1 ∂ψ
D 4 ψ = 0, where D 2 ψ ≡ + 2 .
∂r 2 r ∂θ sin θ ∂θ
ψo = f (r) sin2 θ.
U
B=− , D = 0,
2
98
and so
A 1 2
f= − U r + Cr.
r 2
Before we can find A and C we need to consider the motion within the drop.
ud = uo = 0 on r = a =⇒ ψd = ψo = 0 on r = a. (6.16)
In addition, we need the tangential stress at the surface to be continuous (see previous
discussion on fluid-fluid boundary conditions). Thus we want
∂ 1 ∂ψd ∂ 1 ∂ψo
2µd edrθ = 2µo eorθ on r = a =⇒ κ = , (6.18)
∂r r 2 ∂r ∂r r 2 ∂r
where
κ = µd /µo is the viscosity ratio,
since
r ∂ v 1 ∂u
erθ = +
2 ∂r r 2r ∂θ
in spherical polars, and u is everywhere zero on the drop’s surface.
The form of these conditions at the interface suggest that inside the drop we also look for
a solution of the form
ψd = f (r) sin2 θ.
Thus,
E
f= + F r 2 + Gr + Hr 4 ,
r
for constants E, F, G, H.
Now, using the definition of ψ, we have
2 cos θf (r) 1
u= , v = − sin θf ′ (r).
r2 r
Therefore, to avoid a singular velocity at r = 0, we need E = G = 0.
99
Condition (6.16) implies that
A 1 2
F a2 + Ha4 = − U a + Ca = 0,
a 2
while condition (6.17) demands
A
fd′ (a) = fo′ (a) =⇒ 2aF + 4a3 H = − − U a + C,
a2
and condition (6.18) requires
h 1 i′ h 1 i′
κ 2 fd′ = 2 fo′ on r=a
r r
which implies
2F 4A U 2C
κ − 2 + 4H = 5 + 2 − 3 .
a a a a
Solving for the various constants we find
a3 U κ aU (2 + 3κ) U U
A=− , C= , F = , H=− .
4(1 + κ) 4(1 + κ) 4(1 + κ) 4a2 (1 + κ)
Thus, we have
U r 2 sin2 θ h r2 i
ψd = 1− 2 in the drop.
4(1 + κ) a
Exercise: Verify the above formula for the drag on the moving drop.
As the fluid in the drop becomes more and more viscous, we would expect this result to
approach that for the sphere. Indeed, taking the limit κ → ∞, we obtain
100
Low Reynolds number swimming
The reversibility of Stokes flows has an interesting consequence with regard to the swimming
of very small organisms, such as a spermatozoa. Since the lengthscale of such organisms is
very small, the Reynolds number is likely to be very small and the conditions of Stokes flow
apply.
If this idealised little spermatozoa swishes its tail up and down as shown, since the flow
around it is reversible, it will never get anywhere!
We will look at an idealised model of a swimming motion which overcomes the reversibility
problem.
y1
0.8 λ=2π/k
0.6
0.4
0.2
-0.2
(x s ,ys )
-0.4
sheet
-0.6
-0.8
-1
0 2 4 6 8 10 12
x
A thin flexible sheet ‘swims’ in a viscous fluid under conditions of Stokes flow. A point on
the sheet (xs , ys ) moves according to
xs = x, ys = a sin(kx − ωt),
i.e. the wave moves from left to right with speed c = ω/k. The wavelength λ = 2π/k.
Assume the fluid above the sheet (y ≥ ys ) is in two-dimensional motion, satisfying the
equation
∇4 ψ = 0,
101
where
u = ψy , v = −ψx .
As y → ∞, we require u to be bounded.
Non-dimensionalize by writing
Then ∇4 ψ = 0 becomes
∂2 ∂ 2 2 ∗
k3 a ω + ψ = 0.
∂x∗2 ∂y ∗2
and the boundary conditions,
∂ψ ∗ ∂ψ ∗
u∗ = = 0, −v ∗ = = cos(x∗ − t∗ ) on y ∗ = ka sin(x∗ − t∗ ).
∂y ∗ ∂x∗
Dropping all the asterisks for convenience, and setting X = (x−t), the boundary conditions
are
∂ψ ∂ψ
(x, ε sin X) = 0, (x, ε sin X) = cos X.
∂y ∂x
∂ψ ∂ψ ∂2ψ
(x, ε sin X) = (x, 0) + ε sin X 2 (x, 0) + · · · = 0,
∂y ∂y ∂y
and
∂ψ ∂ψ ∂2ψ
(x, ε sin X) = (x, 0) + ε sin X (x, 0) + · · · = cos X.
∂x ∂x ∂x∂y
In view of these boundary conditions, we seek a solution in the form of the expansion
ψ = ψ0 + εψ1 + ε2 ψ2 + · · · ,
and so we obtain h i h i
∇4 ψ0 + ε ∇4 ψ1 + · · · + εn ∇4 ψn + · · · = 0.
102
The boundary conditions become
∂ n o ∂2 n o
ψ0 + εψ1 + · · · + ε sin X 2 ψ0 + εψ1 + · · · + · · · = 0,
∂y ∂y
and
∂ n o ∂2 n o
ψ0 + εψ1 + · · · + ε sin X ψ0 + εψ1 + · · · + · · · = cos X.
∂x ∂x∂y
Taking the limit ε → 0 yields all the terms in the equations of size O(ε0 ):
∇4 ψ0 = 0,
with conditions
∂ψ0 ∂ψ0
= cos X, =0
∂x ∂y
on y = 0.
At the next ‘order’ we pick out all terms of size O(ε1 ), namely,
∇4 ψ1 = 0,
with conditions
∂ψ1 ∂ 2 ψ0 ∂ψ1 ∂ 2 ψ0
+ sin X = 0, + sin X =0
∂x ∂x∂y ∂y ∂y 2
on y = 0.
We can continue in this vein for all orders O(εn ), n = 2, 3, 4, . . . as far as we like.
We look for a solution to the O(ε0 ) problem by trying
ψ0 = f0 (y) sin X.
Then, since
∂ ∂X ∂ ∂
= = ,
∂x ∂x ∂X ∂X
we get:
∇2 ψ0 = (f0′′ − f0 ) sin X = F (y) sin X, say.
So,
(iv)
= [f0 − 2f0′′ + f0 ] sin X.
103
So
(iv)
f0 − 2f0′′ + f0 = 0.
Writing f0 = emy , we have
=⇒ m = 1, 1, −1, −1.
Accordingly, n o
ψ0 = (A0 + B0 y)ey + (C0 + D0 y)e−y sin X.
But, on y = 0,
∂ψ0
= cos X =⇒ C0 = 1.
∂x
And, again on y = 0,
∂ψ0
= 0 =⇒ −C0 + D0 = 0 =⇒ D0 = 1.
∂y
The boundary conditions suggest that we look for a solution of the form
Then
∇2 ψ1 = f11
′′ ′′
+ (f12 − 4f12 ) cos 2X.
104
So,
(iv) (iv)
∇4 ψ1 = f11 + (f12 − 4f12
′′ ′′
) cos 2X − 4(f12 − 4f12 ) cos 2X
(iv) (iv)
′′
= f11 + (f12 − 8f12 − 16f12 ) cos 2X
= 0.
We therefore require
(iv) (iv)
′′
f11 = 0 and f12 − 8f12 − 16f12 = 0.
Thus,
f11 = A11 y 3 + B11 y 2 + C11 y + D11 .
and so
f12 = (A12 + B12 y) e2y + (C12 + D12 y) e−2y .
Now,
∂ψ1
=0 on y = 0 =⇒ f12 (0) = 0,
∂x
and
∂ψ1 1 ′ 1 ′ 1
= (1 − cos 2X) on y = 0 =⇒ f11 (0) = , f12 (0) = − .
∂y 2 2 2
Hence C11 = 1/2 and, because the velocities are bounded as y → ∞, we need A11 = B11 = 0.
For f12 we shall also require A12 = B12 = 0 in order for the velocities to be bounded at
infinity. Then
f12 (0) = 0 =⇒ C12 = 0.
Also,
′ 1 1
f12 (0) = − =⇒ D12 = − .
2 2
Therefore
1 1
ψ1 = f11 (y) + f12 (y) cos 2X = y + D11 − y e−2y cos 2X.
2 2
105
We can take D11 = 0 without loss of generality (since it is just adding a constant to a
streamfunction). So,
1
ψ1 = y(1 − e−2y cos 2X).
2
Hence
ψ = ψ0 + εψ1 + · · ·
1
= (1 + y)e−y sin X + εy(1 − e−2y cos 2X) + · · · .
2
In dimensional form,
εωa
u = −ωkaye−ky sin(kx − ωt) + 1 + (2ky − 1) e−2ky cos 2[kx − ωt] + · · · .
2
In real life, similar ‘non-symmetric’ motions have to be adopted at low Reynolds numbers
to avoid the reversibility problem.
A more realistic picture of the spermatozoa has it wiggling its tail in a rotary motion,
sending helical waves down to the tip. This overcomes the problems with reversibility and
allows it to move forwards. Bacteria, which often have several tails (or flagellae) adopt
similar non-symmetric strategies for swimming.
106
Now we will look at slow flow of viscous fluid in a sharp corner or crevice of angle 2α. The
fluid satisifes no slip and no normal flow at θ = ±α, and can be thought of as driven by
some agency a distance from the corner r = 0. For example, uniform flow over a sharp
indentation will drive a flow consisting of eddies inside the crevice.
U
θ=α 000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
r 000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
θ=−α 000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
000000000000000000000000000000000000000000000
111111111111111111111111111111111111111111111
The flow is two-dimensional and satisfies the Stokes equations. We use a stream function,
ψ, defined by
1 ∂ψ ∂ψ
u= , v=− ,
r ∂θ ∂r
for velocity field u = u r̂ + v θ̂.
The equation of motion is ∇4 ψ = 0, and so
∂2 1 ∂ 1 ∂ 2 2
+ + ψ = 0. (6.19)
∂r 2 r ∂r r 2 ∂r 2
ψ = r λ f (θ), λ 6= 0
and so
n o n o
∇4 ψ = ∇2 (∇2 ψ) = r β−2 β 2 F + F ′′ = r β−2 β 2 [λ2 f + f ′′ ] + [λ2 f ′′ + f (iv) ] .
β 2 λ2 f + (λ2 + β 2 )f ′′ + f (iv) = 0.
Thus, since β = λ − 2,
107
Seeking a solution of the form f ∝ emθ , we find
Hence,
A cos λα + C cos(λ − 2)α = 0,
Aλ sin λα + C(λ − 2) sin(λ − 2)α = 0,
or,
! !
cos λα cos(λ − 2)α A
=0 (6.20)
λ sin λα (λ − 2) sin(λ − 2)α C
So, for non-trivial solutions, we need the determinant of the matrix to vanish:
108
1
0.8
0.6
sin(x)/x
0.4
0.2
-0.2
0 5 10 15 20
x
If we choose 2α = 2.55, we get B = 0.217, and −0.217 just touches the global minimum of
the solid curve, as shown. So if α < 2.55, there will be no real solutions. In degrees, this
means that there are no real solutions if
2α < 146.3◦ .
is single signed for all r and so the fluid moves in and out of the corner something like this:
109
are clearly solutions to the equation of motion ∇4 ψ = 0, and since this equation is linear
we can add them together, so
1h i
ψ = r λ f (θ) + r λ f (θ) = ℜ{r λ f (θ)}
2
is a solution.
Suppose we fix θ = 0, so we are focussing on the centre line. Then f (0) is just some complex
number, call it c. Since λ = p + iq, the azimuthal velocity component
Writing
r iq = eiq log r = cos(q log r) + i sin(q log r),
we have
Intriguingly, it is clear from (6.22) that there must be an infinite number of these eddies,
as v will repeatedly change sign as r → 0.
These eddies have been visualized in experiments, although the ones very close to the corner
are far too weak to be observed in practice.
7. Lubrication theory
Lubrication theory is the study of thin films of fluid. The key feature of such problems
is that one dimension is much larger than the other. For example, you can get a coin to
“stick” to a surface by putting a small amount of water between the two. A surprisingly
large force is required to pull the coin away from the table. Compared to the diameter of
the coin, the height of the film of water is very small (h ≪ L in the figure below) and as
such is amenable to lubrication analysis. We shall deal with this problem later.
110
z
F F
11111111111
00000000000
00000000000
11111111111
h
x
L
∇ · u = 0, u · ∇u = −∇p/ρ + ν∇2 u
and make some simplifications, using the fact that we are dealing with a thin film of fluid.
Specifically, we work on the assumption that if h and L are typical horizontal (x, y) and
vertical (z) dimensions respectively, and if
h
δ= then δ ≪ 1.
L
Let U be a typical horizontal flow speed. Then, since u = 0 at the solid walls z = 0 and
z = h, u and v must change by an amount of order U as z varies over a distance h. We can
estimate sizes of the terms in the equations using these type of arguments. In the film we
expect
∂u U ∂2u U ∂u U ∂2u U
∼ , ∼ 2, ∼ , ∼ 2,
∂z h ∂z 2 h ∂x L ∂x2 L
with similar estimates for v and its derivatives.
Now, the continuity equation
∂u ∂v ∂w ∂w U Uh
∇·u= + + = 0 =⇒ ∼O =⇒ w ∼ O .
∂x ∂y ∂z ∂z L L
So we can estimate
U2 U2 U 2h
u · ∇u ∼ O i+O j+O k,
L L L2
and
∂2u U U U
∇2 u ≈ ∼ O i + O j + O δ 2 k.
∂z 2 h2 h2 h
These arguments suggest that we can assume
UL
|u · ∇u| ≪ |∇2 u| if δ2 ≪ 1.
ν
111
Defining the Reynolds number R = U L/ν as usual, and introducing the modified Reynolds
number
Rm = δ2 R,
this means that if Rm ≪ 1, we may neglect the nonlinear terms and approximate the thin
film flow with the lubrication equations:
∂p ∂2u ∂p ∂2v ∂p ∂2w
0=− +µ 2, 0=− + µ 2, 0=− +µ 2, (7.1)
∂x ∂z ∂y ∂z ∂z ∂z
∂u ∂v ∂w
+ + = 0. (7.2)
∂x ∂y ∂z
Applications
The applications of the lubrication approximation are manifest. Some examples:
1. Flow in a slider-bearing (see below).
2. Blood flow in slowly-tapering arteries. In some situations, a lubrication model is appro-
priate.
3. Peeling of the retina away from the choroid during retinal detachment, when part of the
vitreous humour in the back of the eye has liquefied. This is a topic of current research.
4. Modelling the thin film of fluid drawn over the eye during a blink. A lubrication model
predicts that this film will break-apart after approximately 15 seconds. (Try not to blink
for longer than that!) This is also the a subject of current research.
V
Solid
U
h
x
L Fluid
Rigid wall
In a typical configuration, such as that shown above, we have two lengthscales, h and L
(and generally speaking we are interested in the case h ≪ L). It is therefore sensible to
112
non-dimensionalize by setting
x = LX and y = hY.
If the solid wall is moving with horizontal speed U , this suggests we take u = U Û . For an
order of magnitude balance of terms in the continuity equation, we expect
v U h
∼ =⇒ v ∼ O U ,
h L L
x = LX, y = hY
2
u = U Û , v = U (h/L)V̂ , p = µ(U L/h )p̂.
Substituting into the two-dimensional Navier-Stokes equations we find
h i
Rm Û ÛX + V̂ ÛY = −p̂X + ÛY Y + δ2 ÛXX ,
h i
Rm Û V̂X + V̂ V̂Y = −p̂Y + δ2 V̂Y Y + δ4 V̂XX ,
ÛX + V̂Y = 0,
0 = −p̂X + ÛY Y ,
0 = −p̂Y
ÛX + V̂Y = 0.
113
Thus, for a lubricating flow, the above analysis shows that we can use an approximate form
of the Navier-Stokes equations, given by
u = v = 0 on y = 0; u = U, v = V, on y = h(x). (7.6)
For the next bit, we will need to use the following result (see Problem Sheet 4),
h(x) h(x)
∂ ∂f dh
Z Z
f (x, y) dy = dy + f (x, h) .
∂x 0 0 ∂x dx
Demanding that v = V on y = h we find
h(x)
∂u
Z
V = − dy,
0 ∂x
h(x)
∂ dh
Z
thus −V = u dy − U . (7.7)
∂x 0 dx
Meanwhile, the flux
h(x) hn
1 h 2 Uy o
Z Z i
u dy = y − yh px + dy,
0 0 2µ h
1 h h3 i U h
= px − +
2µ 6 2
h(x) 3
∂ ∂ 1 h U dh
Z h i
and so u dy = px − + . (7.8)
∂x 0 ∂x 2µ 6 2 dx
114
Combining (7.8) with (7.7), we obtain the Reynolds lubrication equation:
1 ∂ 3 ∂p 1 dh
h =− U + V. (7.9)
12µ ∂x ∂x 2 dx
This equation involves a second order derivative in x and so two boundary conditions are
required. The height function h is normally given. We usually solve (7.9) subject to
p = p1 at x = x1 ,
p = p2 at x = x2 .
Note: The lubrication equation (7.9) also applies when U and V are functions of time, t,
provided that
∂ ∂2
≪ ν 2.
∂t ∂y
We will make use of this fact later.
Assume that the angle α ≪ 1. The pressure on either side of the bearing is atmospheric,
which we are at liberty to set to zero. For convenience we change variable from x to h
and so, since h ≈ αx + const, d/dx = αd/dh. Then the Reynolds Lubrication Equation
becomes:
1 2 d n 3 dp o 1
α h = − U α.
12µ dh dh 2
115
Integrating,
6µU E
p= + F − 2, const E, F.
αh 2h
Choosing E, F such that p = 0 at h = h1 , h2 , we find
12µ U h1 h2 6µU 1
E= , F =− . (7.10)
α h1 + h2 α (h1 + h2 )
Then,
6µU (h2 − h)(h − h1 )
p= > 0 and U > 0 for h1 < h < h2 .
h2 α(h1 + h2 )
If the system is closed and the forces at both ends are negligible, then a force balance shows
that the lift and drag forces on the bearing are equal and opposite to those on the wall, as
shown:
y 1111111111111111111
0000000000000000000
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
L
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
D
0000000000000000000
1111111111111111111
0000000000000000000
1111111111111111111
D x
Now,
∂u h dp U
=− + ,
∂y 2µ dx h
Z h2 n
U 4 6h1 h2 o
=⇒ D = µ − dh.
α h1 h (h1 + h2 )2 h2
And so,
Uh h 6(h − h ) i
2 1 2
D=µ 4 log − .
α h1 h1 + h2
116
Now the lift force, L:
∂v ∂v
Z Z n o
L = −σ22 dx = − − p + µ{ + }y=0 dx
LW ∂y ∂y
h2
p
Z
∂v
(since ∂y = − ∂u
∂x = 0 on y = 0) = dh
h1 α
6µU n h 2(h − h ) o
2 2 1
= log − .
α2 h1 h1 + h2
D αn o
= 1 + O(α)
L 3
L
=⇒ ∼ 1/α and so L ≫ D if α ≪ 1.
D
In other words, the normal force is much larger than the tangential force and the bearing
can withstand heavy loads.
Adhesive problems
We now return to the coin problem mentioned in the introduction to this section. An object
sitting on top of a thin film of fluid on a table can be surprisingly tricky to pick up. We
have already seen how a thin fluid layer can lead to normal forces which are much larger
than tangential ones.
Consider two parallel plane walls, with one of the walls moving away from the other as
shown:
y
h moving wall
We assume that
∂ ∂2
≪ ν 2,
∂t ∂y
so the flow is steady.
117
We need to solve the Reynolds lubrication equation,
1 ∂ 3 ∂p 1 dh
h =− U + V,
12µ ∂x ∂x 2 dx
together with the boundary conditions
u=v=0 on y = 0; u = 0, v = ḣ on y = h(t).
We want to calculate the force, F , resisting the motion. As for the slider bearing, we assume
the forces at the ends are negligible so the force exerted by the fluid on the stationary wall is
equal and opposite to that exerted on the moving wall. We can thus compute the resistance
by calculating minus the force exerted on the lower wall. We know that
Z a Z a
F = − σ22 |y=0 dx, (since vy |y=0 = 0 from cty) = p dx
−a −a
a
6µḣ 8µa3 ḣ
Z
= (x2 − a2 ) dx =− .
h3 −a h3
Note that force is negative so the fluid is trying to pull the moving wall back downwards.
The most interesting thing about this result is that F ∼ O(h−3 ) as h → 0, which implies a
huge resistance force from a very thin film!
118
z
z=h(x,y,t) W
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111 Fluid
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
V
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
U
0000000000000000000000000000000
1111111111111111111111111111111
h
y
∂p ∂2u
∼ µ 2,
∂x ∂y
1 n 2 o Uz
u= px z − zh + ,
2µ h
119
1 n 2 o Vz
v= py z − zh + .
2µ h
Again, we use the continuity equation to find the vertical velocity component:
Z z
−w = ux + vy dz.
0
Thus, at z = h,
Z h
−W = ux + vy dz
0
Z h Z h
∂ ∂
so −W = u dz + v dz − U hx − V hy .
∂x 0 ∂y 0
Using the above results for u and v, and integrating, we obtain the three-dimensional
version of the Reynolds lubrication equation:
∂ h 3 ∂p i ∂ h 3 ∂p i n o
h + h = 6µ − U hx − V hy + 2W (7.13)
∂x ∂x ∂y ∂y
Sphere Fluid
h
x
y
z
We assume that the sphere falls slowly enough to justify a steady lubrication analysis. In
other words, we assume that ∂/∂t ≪ ν∂ 2 /∂z 2 .
On the surface of the sphere u = v = 0 and w = ḣ, so U = V = 0 and W = ḣ. In this case,
Reynolds lubrication equation states that
∂ h h3 ∂p i ∂ h h3 ∂p i
+ = µW.
∂x 12 ∂x ∂y 12 ∂y
120
Appealing to the symmetry of the problem, we expect that p = p(r), where r 2 = x2 + y 2 .
It is therefore more sensible to discuss the problem in cylindrical polar coordinates.
r
θ
x
y
z
1 d h 3 dp i
rh = 12µW,
r dr dr
(integrating) rh3 pr = 6µW r 2 + c. (7.15)
ϕ a
r
h(r)
ho
Now, within the lubrication limit, we found above that h = h0 + r 2 /2a. Thus, when r is
large, but not large enough to destroy the validity of the lubrication approximation, we have
121
h ∼ O(r 2 ). Using the result above for the pressure, it follows that
r 1
pr ∼ O ∼ O for large r.
h3 r5
h i∞
We can thus set πr 2 p(r) = 0 and deduce that
0
∞ ∞
2 dp r3
Z Z
F = dr = −6µπW πr dr
0 dr 0 [h0 + r 2 /2a]3
12µπa2 W ∞ ξ
Z
(setting r 2 = 2ah0 ξ) = − dξ
ho 0 (1 + ξ)3
6µπa2 W
Z ∞
ξ 1
So F = − since dξ = .
ho 0 (1 + ξ)3 2
An interesting consequence of this is that the sphere never reaches the wall!
ho
mg
ḣ0
mḧ0 = −6µπa2 − mg.
h0
This can be written as
d
ḧ0 + κ
(log h0 ) = −g
dt
for positive constant κ. Integrating once, we find
h
0
ḣ0 + κ log = −gt + Ḣ0 ,
H0
122
Flow of a thin film down a sloped wall
A thin film of viscous fluid spreads along a plane wall which is inclined at an angle α to the
horizontal, as shown in the diagram. The nose of the film is at xN (t).
w
z
111111111111111111
000000000000000000
u
000000000000000000
111111111111111111
000000000000000000
111111111111111111
000000000000000000
111111111111111111 g
000000000000000000
111111111111111111
000
111
000000000000000000
111111111111111111
000
111
000000000000000000
111111111111111111
000
111
000000000000000000
111111111111111111
000
111
000000000000000000
111111111111111111
x111
000
N
α
x
Starting with the two-dimensional lubrication equations (7.3, 7.4, 7.5), we need to add body
force terms to account for gravity. We therefore write
where z is now playing the role of the vertical coordinate and w the vertical fluid speed.
Integrating the second equation,
p = ρg[h − z] cos α + pa ,
ii) the tangential stress must equal zero. Since we are implicitly assuming a rather slender
film, the normal vector points almost vertically. Consistent with the approximations already
made in deriving the lubrication equations, we take the normal to point exactly vertically6
and require
∂u ∂w
2µexz = µ + =0
∂z ∂x
as the zero tangential stress condition. Note that the tangential direction is almost in the
x direction, and is taken exactly in the x direction to within a consistent approximation.
6
This is clearly not true near the nose of the film where the surface is highly curved. As a first estimate,
we ignore this fact and proceed regardless. The results turn out to be surprisingly accurate, a not uncom-
mon feature of lubrication-type analyses. More accurate calculations can be made taking into account the
curvature at the nose, but these are beyond the scope of the current course.
123
But in the lubrication approximation,
∂u ∂w
≫ ,
∂z ∂x
and so we need ∂u/∂z = 0 on z = h.
Substituting for the pressure, (7.16) is now,
Again, appealing to the fact that the film is thin, we can drop hx as it is expected to be
small. This leaves
νuzz = −g sin α,
which is easily integrated. Applying no slip at z = 0 and the zero tangential stress condition
at z = h we therefore find
g sin α 1
u= (hz − z 2 ).
ν 2
As usual, we can find the vertical velocity component by integrating the continuity equation:
g sin α
wz = − hx z,
ν
g sin α
(using w = 0 at z = 0) =⇒ w = − hx z 2 .
2ν
Lastly we need a kinematic condition at the free surface z = h (see the section on boundary
conditions earlier in the course). In fact, we need:
D
(z − h) = 0,
Dt
∂h ∂h
i.e. w = +u on z = h(x, t).
∂t ∂x
124
which is only expected to be valid at large times.
Since the total volume of the film,
Z xN (t)
h(x, t) dx = A, say,
0
Note: Despite neglecting effects associated with the curvature around the nose of the
spreading film, the prediction for the rate of spread agrees well with experiments (see
Huppert (1986), Journal of Fluid Mechanics, vol. 173, pp. 557-594. This journal is available
in the Library).
125
wall
liquid
air
The thick-lined arrows indicate the direction of the driving Marangoni force γθ = −2 sin 2θ.
Notice that four viscous eddies are established circling the interior of the film.
The presence of the eddies might seem surprising at first. How can the fluid turn a corner
and not contravence the restrictions of the lubrication analysis? This stipulates that while
the streamwise velocity (here the azimuthal component) u is O(1), the vertical component
(here the radial component) should be O(δ). Yet as an eddy turns a corner, u is zero and
v is comparatively large.
Careful thought shows that the restrictions are not in fact contravened. All that is happen-
ing is that, as θ increases, u passes through zero. Provided u still varies slowly with respect
to θ, the requirements of lubrication theory are still met.
8. Boundary layers
The Reynolds number, defined as R = U L/ν, for typical velocity and lengthscales U , L,
126
appears in the non-dimensional form of the Navier-Stokes equations:
1 2
∇ · u = 0, u · ∇u = −∇p + ∇ u. (8.1)
R
When R is very small we have Stokes flow, which we have already discussed.
What happens when R is very large? Letting the Reynolds number in equations (8.1) tend
to infinity, it would appear that we can just drop the viscous term ∇2 u/R on the grounds
that it is very small. This would leave us with Euler’s equations governing an inviscid fluid.
In dimensional form these are
∇ · u = 0, ρu · ∇u = −∇p.
Suppose we ask the question: what is the drag on an object placed in a high Reynolds number
flow?
Ui
From section ??, D’Alembert’s Paradox tells us that the drag on the object is precisely zero!
But this is clearly at odds with experience. Experiments of high Reynolds number flows
past a circular cylinder, for example, predict a drag on the cylinder which is most certainly
not zero.
What is going on?
εf ′′ + f ′ = a, f = f (y), (8.2)
(1 − a)
f = ay + (1 − e−y/ε ).
(1 − e−1/ε )
127
For y = 0(1), i.e. not small,
f ∼ ay + 1 − a = f0 (y), say.
This can be thought of as an outer, ‘inviscid’ solution, valid in a region away from y = 0.
It can be seen to be the solution of equation (8.2) if we naı̈vely set ε to zero on the grounds
that it is small. Crucially, however, it satisfies the condition at y = 1 but not that at y = 0.
So it is not the complete solution.
Suppose we now focus on what happens when y is small. To this end, we introduce the
change of variables y = εY , so that if Y is of order one, y will be small. The equation
becomes
fY Y + fY = εa.
Since ε ≪ 1, to a first approximation we may set the right hand side to zero (as we did to
obtain f0 (y) above), leaving
fY Y + fY = 0.
Integrating this equation yields
It is noticeable that this function satisfies the condition at y = Y = 0 but not that at y = 1.
To summarize, we have an inner solution fi , which satisfies the boundary condition at y = 0,
and an outer solution, fo , which satisfies the boundary condition at y = 1. We can tie these
solutions together by using the following matching procedure
lim fi = lim f0 .
Y →∞ y→0
fi = (1 − a)(1 − e−Y ).
128
The following graphs illustrate these solutions when ε = 0.1 and ε = 0.01.
0.8
0.6
f
0.4
outer
0.2 inner
exact
0
0 0.2 0.4 0.6 0.8 1
y
0.8
0.6
f
0.4
0.2 outer
inner
exact
0
0 0.2 0.4 0.6 0.8 1
y
Clearly the inner solution only performs well near y = 0, while the outer works well away
from y = 0. The agreement between the approximate inner and outer solutions and the
exact solution improves as ε → 0.
The key feature of the above toy problem is the clear development of a region (here near
y = 0) in which there is a rapid change of the solution. The outer solution has only a slight
slope. However, near y = 0, the inner solution takes over and exhibits a dramatic change
over a very short distance. In this very thin region, it is not appropriate to drop the term
εf ′′
129
over most of the flow field. However, very close to the surface of the body there is a boundary
layer in which this term is no longer negligible.
The flow is taken to be steady. The external stream consists of inviscid flow moving with
velocity u = U i, with U constant.
Let L be a characteristic horizontal lengthscale. If the plate is finite, L can be chosen to be
its length. The Reynolds number is
UL
R= ≫ 1,
ν
and so is assumed to be very large.
Following the above discussion, we anticipate that the viscous term R−1 ∇2 u will only be
important in a very thin boundary layer of height δ(x) close to the plate.
We aim to derive a set of approximate equations valid in the boundary layer. We will first
do this using an argument based on the estimation of the order of magnitude of all the
terms. Later we will go through a more formal derivation.
The continuity equation states
∂u ∂v
+ = 0.
∂x ∂y
Since u varies from zero at the plate (no-slip) to U in the external stream, we can say that
u ∼ U, typically.
Changes in the boundary layer are expected to occur on a lengthscale of typical size L
parallel to the plate and δ perpendicular to it (i.e. there is much more rapid change in the
y direction than the x direction). We can therefore estimate
∂ 1 ∂ 1
∼ , ∼ .
∂x L ∂y δ
If V is a typical size for the vertical velocity v, a balance of terms in the continuity equation
suggests that
U V δ
∼ =⇒ V ∼ U.
L δ L
130
Now consider the x-momentum Navier-Stokes equation
1
uux + vuy = − px + ν[uxx + uyy ].
ρ
Prandtl’s Hypothesis
This states that the inertial term, uux , should balance the viscous one, νuyy , inside the
boundary layer, i.e.
U2 U
uux ∼ νuyy =⇒ ∼ ν 2.
L δ
Simplifying, we obtain the following estimate of the boundary layer thickness,
δ ν 1
2
∼ ,
L UL
which can be re-stated as
δ ∼ LR−1/2 .
So viscous effects are only expected to be important in a layer of thickness R−1/2 . Since R
is very large, this layer is very thin.
Continuing, we can estimate
V U δ/L U2
uvy ∼ U =U = ∼ uux ,
δ δ L
and so these two terms balance.
If we choose a scale for the pressure so that
1
px ∼ uux =⇒ p ∼ ρU 2 ,
ρ
131
then, in summary, the boxed terms in the x-momentum equation are comparable for large
R,
uux + vuy = − 1ρ px + ν(uxx + uyy ).
This states that the pressure does not vary across the boundary layer.
Collecting everything together, we have the Prandtl boundary layer equations
1
uux + vuy = − px + νuyy , (8.3)
ρ
py = 0,
ux + vy = 0.
u = v = 0 on y = 0; u→U as y → ∞.
132
No condition on v can be imposed as y → ∞ since the term vyy has been neglected.
The above all follows through if U is a function of x, that is if the flow at infinity (the
external stream) is non-uniform.
The special feature of the boundary layer flow, namely that py = 0, allows us to find the
pressure gradient term px by relating it to the external stream velocity. If we let y → ∞ in
equation (8.3), since u → U (x) in general, we obtain
1
U Ux + 0 = − px + 0.
ρ
The first zero follows since uy → ∂U/∂y = 0 and the second for a similar reason.
Hence,
px = −ρU Ux ,
and we can write the x-momentum boundary layer equation (8.3) as
x = Lx∗ , y = Ly ∗ , u = U u∗ , v = U v∗ , p = ρU 2 p∗ .
Here a ∗ denotes a dimensionless variable. Note that we have chosen ρU 2 as a suitable scale
for the pressure. Check for yourself that ρU 2 has dimensions of pressure.
We now enter these definitions into the Navier-Stokes equations. Let’s just consider the
x-momentum equation:
∂u ∂u 1 ∂p
u +v =− + ν∇2 u
∂x ∂y ρ ∂x
So, we obtain,
U 2 ∗ ∂u∗ ∗ ∂u
∗ U 2 ∂p∗ νU ∂ 2 u∗ ∂ 2 u∗
u + v = − + + ∗2 .
L ∂x∗ ∂y ∗ L ∂x∗ L2 ∂x∗2 ∂y
∂u∗ ∗ ∂u
∗ ∂p∗ ν ∂ 2 u∗ ∂ 2 u∗
u∗ + v = − + + ∗2 .
∂x∗ ∂y ∗ ∂x∗ U L ∂x∗2 ∂y
133
To save writing, we now drop the asterisks and write
∂u ∂u ∂p 1 ∂2u ∂2u
u +v =− + + 2 ,
∂x ∂y ∂x R ∂x2 ∂y
where
UL
R=
ν
is the Reynolds number.
A similar procedure follows for all of the momentum equations, and for the continuity
equation. Finally, we have the two-dimensional Navier-Stokes equations in dimensionless
form:
1
uux + vuy = −px + ∇2 u,
R
1
uvx + vvy = −py + ∇2 v,
R
ux + vy = 0,
with R = U L/ν and where u, v have been scaled on U , p on ρU 2 , and x, y on L.
In the boundary layer, we expect
∂
∼ R−1/2 for y ≪ 1.
∂y
134
These are the boundary layer equations in dimensionless form. They can be viewed as
the leading order asympotic approximation to the Navier-Stokes equations when R → ∞,
y ∼ R−1/2 .
Restoring the dimensions, we have
∂u ∂u 1 ∂p ∂2u
u +v = − +ν 2,
∂x ∂y ρ ∂x ∂y
∂p
0 = ,
∂y
∂u ∂v
+ = 0.
∂x ∂y
Since the flow is two-dimensional, it may be convenient to write the boundary layer equa-
tions in terms of a stream function ψ, where
∂ψ ∂ψ
u= , v=−
∂y ∂x
We find:
1
ψy ψxy − ψx ψyy = − px + νψyyy ,
ρ
py = 0.
Boundary conditions
i) No slip requires
u=v=0 on y = 0,
or, equivalently,
ψ = ψy = 0 on y = 0.
u → U (x) as y → ∞, so ψy → U as y → ∞,
which leads to
−px = U Ux ,
where the external stream velocity U is calculated from inviscid theory.
Notes
1. The equations apply on a curved surface, with x, y measuring distance along and normal
to the surface respectively.
135
y x
Body
2. The stream function version of the equations are 3rd order in y. Therefore we can only
apply 3 boundary conditions on ψ. The vertical velocity v = ∂ψ/∂x is fixed on solving the
equations.
3. The equations are parabolic in x. This means that if we specify ψ(= f (x), say) at x = x,
then ψ is determined for x > x. There can be no upstream influence. This means that the
flow at location x = x1 has no effect on that at x = x0 if x0 < x1 . So, if we introduce
a small disturbance (say an imperfection in the wall) at x = x0 , only the shaded region
downstream is affected.
x0
This opens up many interesting questions about boundary layer flows. How, for example,
does a boundary layer negotiate a small bump on the wall? Common experience suggests
that the streamlines will start to move upwards as they near the bump on the upstream
side.
1111111111
0000000000
0000000000
1111111111
But note 3 says that this is not allowed since the bump may not affect the flow upstream.
Famously, some experimentalists showed that upstream effects may occur when a shock
wave impinges on a boundary layer in a supersonic flow.
1111111111111111111111111111111111111111
0000000000000000000000000000000000000000
Shock wave
U
Upstream influence
Boundary layer
1111111111111111111111111111111111111111
0000000000000000000000000000000000000000
1111111111111111111111111111111111111111
0000000000000000000000000000000000000000
136
Blasius’ solution
The solution of the boundary layer equations for high speed flow past a flat plate was given
by Blasius in 1908. We take the case when the external stream, U , is constant, = U0 , say.
Since ∂U0 /∂x ≡ 0, there is no streamwise pressure gradient and the boundary layer problem
to be solved is as follows
∂u ∂u ∂2u
u +v = ν , (8.5)
∂x ∂y ∂y 2
∂u ∂v
+ = 0.
∂x ∂y
u=v=0 on y = 0; u → U0 as y → ∞.
We seek a similarity solution to this problem. To motivate this, consider briefly the flow
when the flat plate is of finite length, L.
From the above arguments, at any x station in the boundary layer the flow cannot be
affected by what is happening further downstream. In other words, the oncoming flow is
unaware that the plate is of length L. Any solution we try to construct should therefore
reflect this property.
We used a scaling argument to show that the typical boundary layer thickness is
νL 1/2 UL
δ ∼ LR−1/2 = , since R = .
U ν
At a general x station, we expect δ to depend on ν, U and x (but not L). Since
L2 L
[ν] = , [U ] = , [x] = L
T T
then, since [δ] = L, it must appear as the combination
νx 1/2
δ∼
U
which is independent of L.
Now, inside the boundary layer, δ provides a typical scale for y so
νx 1/2
y∼δ∼
U
i.e. y ∼ x1/2 .
137
This motivates us to try looking for a solution which is a function of the order one variable
U 1/2
η= y.
νx
∂ψ ∂ψ U0 1/2 ∂ψ
u= = = U0 .
∂y ∂η 2νx ∂η
∂ψ h∂ ∂ ih i
v = − =− + (2νU0 x)1/2 ψ
∂x ∂x ∂η
h ∂ψ ψ η ∂ψ i
= −(2νU0 x)1/2 + −
∂x 2x 2x ∂η
Furthermore,
∂u ∂u ∂η h U0 i1/2 ∂ 2 ψ
= = U0 2 .
∂y ∂η ∂y 2νx ∂η
2
∂ u h U0 i 3
∂ ψ
2
= U0 3 .
∂y 2νx ∂η
∂2ψ ∂2ψ h ∂ψ ∂ 2 ψ ∂ 2 ψ ∂ψ i
+ ψ = 2x − .
∂η 3 ∂η 2 ∂η ∂η∂x ∂η 2 ∂x
138
The boundary conditions are
∂ψ
ψ= = 0 on η=0 [u = v = 0 on y = 0],
∂η
and u → U0 as y → ∞ which becomes
∂ψ
→1 as η → ∞.
∂η
The boundary conditions suggest that we look for a solution which is purely a function of
η. Thus, letting ψ = f (η), we find that f satisfies
f ′′′ + f f ′′ = 0,
with conditions
f = f ′ = 0 on η = 0, f′ → 1 as η → ∞.
Notes
1. f is called the Blasius function. The existence and uniqueness of f was proved by Weyl
in 1942.
2. f must be calculated numerically.
3. In terms of f ,
u = U0 f ′ (η).
f′ → 1 as η → ∞.
We can do the adjusting of the value of f ′′ (0) using Newton-Raphson iteration, for example.
139
A neat trick
In fact it turns out that we can get away with only having to do one integration over η. To
see this, suppose we integrate the system
F ′′′ + F F ′′ = 0,
with conditions
F = F ′ = 0 and F ′′ (0) = 1 on η = 0.
F ′ → κ as η → ∞,
for some constant κ, which will not in general (unless we’re really lucky) equal one.
Now let f (η) = aF (aη) for constant a. Then f satisfies
f ′′′ + f f ′′ = 0,
and
f (0) = f ′ (0) = 0.
Furthermore,
f ′ (∞) = a2 F ′ (∞) = a2 κ.
So all we need to do is to solve the system for F , and we get the solution, f , to the Blasius
problem straightaway by taking
u/U
1
140
Shear stress on the plate
Suppose we wish to calculate the drag in the x direction on a portion of the plate of length
L. To do this, we will need to know
∂u U 1/2
0
σ12 |y=0 = µ =µ f ′′ (0).
∂y y=0 2νx
Displacement thickness
The displacement thickness provides a measure of the thickness of the boundary layer in
terms of the amount by which the oncoming streamlines have been displaced.
It is defined to be ∞
u
Z
δ1 = 1− dy,
0 U0
where u = U0 f ′ .
Physical interpretation
As stated above, the displacement thickness is a measure of the upwards displacement of
the streamlines of the flow by the boundary layer.
y
u(y)
u
U0
141
y
δ 1 1111111111111
0000000000000
1111111111111
0000000000000
0000000000000
1111111111111
0000000000000
1111111111111
0000000000000
1111111111111
0000000000000
1111111111111 u
U0
where δ1 is the corresponding ‘thickness’ of the uniform stream needed for the same flux
level.
The shaded areas in these two figures are the same.
Rearranging, we have
∞
u
Z
δ1 = 1− dy,
0 U0
for the displacement thickness.
In addition, we may define the ‘momentum thickness’
Z ∞
u u
θ1 = 1− dy.
0 U0 U0
Asymptotic form of f
We know that f ′ → 1 as η → ∞, but can we say anything more about the manner in which
f ′ approaches this unit value for large η?
Since f ′ → 1 as η → ∞, we have f → η − η as η → ∞ for constant η. Suppose we now
write
f → η − η + g(η) as η → ∞,
where g ≪ 1. We aim to determine the function g.
Substituting into the equation f ′′′ + f f ′′ = 0 we have
142
But, since g is small, we can neglect the nonlinear term and just take
g′′′ + (η − η)g′′ = 0.
φ′ + (η − η)φ = 0.
Hence,
2 /2
φ = Ae−(η−η) , constant A.
So, as η → ∞, we have
2 /2
f ∼ η − η + Ae−(η−η) .
In fact, had we not neglected the nonlinear terms above, we would find that, more precisely
the constant A should be replaced by
0.331
− .
(η − η)2
Falkner-Skan solutions
The Blasius solution is a special case of a broader class of similarity solutions to the boundary
layer equations known as Falkner-Skan solutions.
To motivate these, consider a general external stream velocity U (x), which is assumed to
be given. We seek solutions of the boundary layer equations of the form
y
ψ = U (x)g(x)f (η), where η = .
g(x)
U(x)
g(x)
This form is a more general version of that presented for the Blasius solution. By analogy
with that solution, g(x) here represents the ‘thickness’ of the layer.
The boundary layer equations to be solved are
143
with
ψ = ψy = 0 on y = 0; ψy → U (x) as y → ∞.
Now,
u = ψy = g(x)U (x)ηy fη = U f ′ ,
ηU ggx f ′
−v = ψx = (U g)x f − , since ηx = −ηgx /g,
g
U ′′ U ′′′
ψyy = f , ψyyy = f ,
g g2
1h i
ψxy = (U g)x f ′ − U gx (ηf ′′ + f ′ ) .
g
f ′′′ + αf f ′′ + β(1 − f ′2 ) = 0,
where
να = g(U g)x ,
νβ = g2 Ux .
f will satisfy an ordinary differential equation provided that α and β are constants. Oth-
erwise the equation depends on x and η. So henceforth, we assume α and β are constant.
Consider the expression
g2 U = ν(2α − β)(x − x0 ),
g2 Ux β n β
2
= = =⇒ n = ,
g U (2α − β)x x (2α − β)
144
Therefore, x n
U = U0 ,
L
for constant U0 , L. Also,
h ν(2α − β)x i1/2 h ν(2α − β)L i1/2 x 1−n
2
g= = .
U (x) U0 L
Taking α = 1, we have
2n
β=
n+1
and the Falkner-Skan solutions are given by
x n n + 1 1/2 U 1/2
U = U0 , η= y,
L 2 νx
2 1/2
ψ= (U νx)1/2 f (η),
n+1
where f satisfies
f ′′′ + f f ′′ + β(1 − f ′2 ) = 0,
with
f (0) = f ′ (0) = 0, f ′ (∞) = 1.
Notes
1. External flows with U ∝ xn arise for inviscid flow over a wedge of half-angle φ =
nπ/(n + 1).
x
n
x
U
2φ
145
x
stagnation point
In this case g(x) is a constant, and so the boundary layer at the wall has constant thickness.
This solution was discussed before, when we noted that it formed an exact solution of the
Navier-Stokes equations.
4. For some values of n there is more than one solution.
5. If U Ux < 0 so that px > 0 and the pressure gradient is adverse (so that the pressure
increases as you move downstream), eventually the boundary layer will tend to separate
from the wall.
xs
The boundary layer equations break down at the separation point xs . Here the flow en-
counters the Goldstein singularity and the boundary layer solution cannot be continued for
x > xs .
The jet shoots out of a narrow slot in a wall into another fluid which is at rest. The pressure
at infinity in the stationary fluid is asummed to be constant, equal to p0 .
We assume that the jet is very thin, so that the variation of velocity across the jet is
146
extremely rapid. Thus
∂u ∂u
≫ .
∂y ∂x
with
Thus,
u = ψy = xα+β f ′ ; −v = ψx = xα−1 {αf + βηf ′ },
ux = xα+β−1 {(α + β)f ′ + βηf ′′ }; uy = xα+2β f ′′ ; uyy = xα+3β f ′′′ .
xα+β f ′ .xα+β−1 [(α + β)f ′ + βηf ′′ ] − xα−1 [αf + βηf ′ ]xα+2β f ′′ = νxα+3β f ′′′ .
Simplifying,
So we need
α = β + 1. (8.7)
To obtain another relation between α and β, we integrate the momentum equation across
the jet as follows. Z ∞ Z ∞ h i∞
uux dy + vuy dy = ν uy .
−∞ −∞ −∞
147
Applying the boundary condition that u → 0 as y → ±∞, the right hand side vanishes.
Now, integrating by parts,
Z ∞ Z ∞
∞
vuy dy = [vu]−∞ − uvy dy
−∞ −∞
Z ∞
(since u(±∞) = 0) = − uvy dy
Z ∞−∞
(using cty) = uux dy.
−∞
and thus, Z ∞
u2 dy = M,
−∞
for constant M .
Inserting the form of the similarity solution and noting that, since the integration is per-
formed at a fixed x value,
∂η
dη = dy,
∂y
we find
Z ∞
x2α+β f ′2 dη = M, (8.8)
−∞
α = 1/3, β = −2/3.
3νf ′′′ + f f ′′ + f ′2 = 0.
148
Integrating again,
1 B2
3νf ′ + f 2 = ,
2 2
for constant B.
Rearranging,
df
Z Z
6ν = dη,
B − f2
2
We know that
1 1 + x
tanh−1 x = ln .
2 1−x
Hence, we have
6ν f
tanh−1 = η + C.
B B
for constant C. But f (0) = 0 =⇒ C = 0, and so
Bη
f (η) = B tanh .
6ν
Therefore
B2 Bη
f ′ (η) = sech2 ,
6ν 6ν
and so
∞
2 B3
Z
f ′2 dη = .
−∞ 9 ν
149
Horizontal velocity profile
0.2
0.15
0.1
u 0.05
0
-40 -20 0 20 40
y
The ‘width’ of the jet corresponds to where the argument of sech2 is O(1); that is, where
the velocity components are significant in magnitude. Thus, where
6ν
y∼O x2/3 .
B
9. Hydrodynamic stability
150
number. On this point, it is instructive to draw an analogy with the solution of some simple
algebraic equations.
Model problems
(i) Consider the following simple quadratic equation to be solved for u:
a − (u − u0 )2 = 0, a = R − Rc .
The solution is
u = u0 ± (R − Rc )1/2 .
There are three possible cases:
• R = Rc : One solution, u = u0
U0
R
Rc
No solns 2 solns
au − u3 = 0, a = R − Rc .
151
The solution is
u = 0, u = ±(R − Rc )1/2 .
There are two cases to consider:
• R ≤ Rc : The unique solution is u = 0
• R > Rc : There are three solutions: u = 0 and u = ±(R − Rc )1/2 .
If we plot u versus R, we have:
o Rc R
1 soln 3 solns
In the event of multiple solutions to the Navier-Stokes equations, the question arises:
Which solution would you observe in a laboratory experiment?
To answer this question, we need to assess a solution’s stability.
Stability theory
The basic idea is to take a solution of the equations and perturb it a little bit. Physically, we
can think of adding a small disturbance into a flow and observing whether the flow settles
back down to its original state, or if it changes dramatically into a different kind of flow
altogether.
We will only consider linear stability. This means that we will only consider infinitesimal
disturbances.
We call the flow whose stability is to be assessed the basic flow or basic state.
As usual, we take the flow to be governed by the Navier-Stokes equations, written here in
dimensionless form:
1
∇ · u = 0, ut + u · ∇u = −∇p + ∇2 u, (9.1)
R
152
for a suitably defined Reynolds number.
11111111111111111111111111111111
00000000000000000000000000000000 x
00000000000000000000000000000000
11111111111111111111111111111111
To this flow we add a small time-dependent disturbance, writing the new, perturbed flow
as
and ε ≪ 1.
We perturb the pressure field in a similar way, writing
p = PB + εp̂.
Substituting this, together with (9.3), into the Navier-Stokes equations (9.1), we obtain
∇ · U + ε∇ · û = 0,
1 2 1
εut + U · ∇U + εû · ∇U + εU · ∇û + ε2 û · ∇û = −∇PB − ε∇p̂ + ∇ U + ε ∇2 û.
R R
But, from (9.2),
1 2
∇ · U = 0, 0 = −∇pB + ∇ U
R
So,
1
εut + U · ∇U + εû · ∇U + εU · ∇û + ε2 û · ∇û = −ε∇p̂ + ε ∇2 û.
R
153
Also, U · ∇U = 0 since U(y) represents a unidirectional flow (see previous work on unidi-
rectional flows). Neglecting the term of O(ε2 ) as being smaller than the others, we are left
with
1 2
∇ · û = 0, ût + û · ∇U + U · ∇û = −∇p̂ + ∇ û. (9.4)
R
It is usual at this stage to assume that the disturbances û, v̂ manifest themselves in the
form of travelling waves. We therefore write
So the forms (9.5) represent travelling waves moving with speed cr of amplitude ekci t . Always
k > 0, so it follows that if ci > 0, the amplitude of the wave will grow, the disturbance will
get bigger. Eventually, the disturbance will grow so large, it will completely disrupt the
basic flow.
To summarise,
154
We first consider equation (9.4) in the limit R → ∞, and assume that it is all right to drop
the final term,
∇ · û = 0, ût + û · ∇U + U · ∇û = −∇P̂ .
Substituting in the normal modes, we have, writing the equations in component form,
−ik(c − U )u + U ′ v = −ikP,
−ik(c − U )v = −P ′ ,
iku + v ′ = 0.
It is convenient to try to eliminate the pressure from these equations. Differentiating the
first with respect to y, multiplying the second by −ik and adding, we obtain
Then, using the third equation to note that iku + v ′ = 0 and to write u′ = −v ′′ /ik, we have
and so
n U ′′ o
v ′′ + − k2 v = 0.
(c − U )
Example 1
Consider now the specific problem of parallel flow in a channel 0 ≤ y ≤ 1:
y
1111111111111111111111111111111
0000000000000000000000000000000
0000000000000000000000000000000
1111111111111111111111111111111
1111111111111111111111111111111
0000000000000000000000000000000 x
This is an eigenvalue problem for c. That is, given a particular disturbance of frequency k,
we seek values of c such that this problem has a solution.
155
Physically, we might think of a flexible ribbon being placed within the channel. The ribbon
is vibrated so that it produces fluctuations of fixed wavenumber k. We are interested in the
response of the fluid to this disturbance.
As above, we write c = ci + ici and say that the flow is stable if ci < 0 and unstable if
ci > 0.
The system (9.6) can be solved numerically. However, some analytical progress is still
possible. If we multiply (9.6) by the complex conjugate of v, namely v, and then integrate
over the channel width, we get
Z 1 Z 1n
U ′′ o
′′
v v dy + − k2 |v|2 dy = 0.
0 0 (c − U )
Integrating by parts,
1 1n
i1 U ′′
h Z Z o
′ 2
vv ′
− |v | dy + − k |v|2 dy = 0.
2
0 0 0 (c − U )
Now
1 1 cr − U − ici
= = ,
c−U cr + ici − U |c − U |2
and so the imaginary part of (9.7) gives,
Z 1 ′′ 2
U |v|
ci 2
dy = 0.
0 |c − U |
So if ci > 0, the only way for this equation to hold is if U ′′ changes sign at some point in the
interval [0, 1]. So U must have an inflection point in that interval. This gives a necessary
but not a sufficient condition for instability of the flow.
It is known as Rayleigh’s inflection point theorem.
156
As before we use the method of normal modes, writing
Therefore we have
iku + v ′ = 0, (9.9)
1
−ikcu + ikU u + vU ′ = −ikp + [−k2 u + u′′ ] (9.10)
R
1
−ikcv + ikU v = −p + [−k2 v + v ′′ ].
′
(9.11)
R
1 h d2 2
i2 d2
2
− k v = (U − c) − k v − U ′′ v. (9.12)
ikR dy 2 dy 2
This is the Orr-Sommerfeld equation.
1111111111111111111111111111111
0000000000000000000000000000000 x
The flow is driven by a constant negative pressure gradient dp/dx = −G, with G > 0 and
the basic solution is U (y) = y(1 − y). We can now recognise this as Plane Poiseuille Flow
discussed above in the Exact Solutions section.
The disturbance flow is governed by the Orr-Sommerfeld equation (9.12) with boundary
conditions
u = v = 0 on y = 0, 1.
Using equation (9.9) these can be re-expressed as
v = v ′ = 0 on y = 0, 1.
157
For this example we know that U = y(1 − y). We now select a wavenumber k (so we restrict
attention to disturbances of this wavenumber) and a Reynolds number, R. Then we have
to solve the problem
1 ′′′′
[v − 2k2 v ′′ + k4 v] = (y − y 2 − c)(v ′′ − k2 v) + 2v. (9.13)
ikR
with
v = v ′ = 0 on y = 0, 1,
for c. As before, writing c = cr + ici ,
Solving for c is generally a numerical task for a computer. Computing c in this way for many
values of k and R we can plot the neutral curve, which delineates the boundary between
stable and unstable solutions. For the current example the neutral curve looks like this:
k
cI < 0 stable
kc
cI > 0
unstable
R
Rc
Every point inside the hoop corresponds to a value of c with positive imaginary part, i.e.
an unstable flow. Every point outside has ci < 0, i.e. a stable flow. On the curve itself,
ci = 0. Modes corresponding to values of c lying on the curve are known as neutral modes.
Note: Both of the tails of the hoop asymptote to the R axis as R → ∞.
As we can see from the picture, there is a critical Reynolds number beneath which there
are no unstable modes. It is computed to be
Rc ≈ 5774.
The corresponding wavenumber at the very nose of the neutral curve is found to be
kc = 102.
158
A surprise
We have just looked at the stability of flow in a channel. Suppose we instead look at that
in a circular pipe. In cylindrical polars, the basic flow is now,
U = 1 − r2,
û = ûr̂ + ŵk,
(assuming there is no flow component in the azimuthal (i.e. θ) direction), we can write down
the corresponding form of the Orr-Sommerfeld equation for cylindrical polars, substitute in
for U , pick a wavenumber k and Reynolds number R and solve for c. Somewhat surprisingly
we find that ci is always negative. In other words, the flow is stable at any value of the
Reynolds number. As was mentioned earlier, this contradicts experiments which show that
the flow becomes turbulent at a critical Reynolds number. The resolution of this problem is
still not complete, but it is generally believed that the turbulence results from perturbations
of finite rather than infinitesimal amplitude.
We now summarize known stability results for the classic flows discussed in this course.
Plane Poiseuille flow
y
1111111111111111111111111111111
0000000000000000000000000000000
0000000000000000000000000000000
1111111111111111111111111111111
1111111111111111111111111111111
0000000000000000000000000000000 x
This flow (see above) is linearly unstable at sufficiently large Reynolds numbers
159
Capillary instability of a liquid jet
Lord Rayleigh (1879) conducted a famous study of the stability of a jet of liquid shooting
into an ambient gas.
perturbed interface
gas r
liquid
x
Fix a frame of reference moving with the jet, so that in this frame the jet appears stationary.
We assume that the flow in the liquid is incompressible and irrotational. This means that
we can model it using Laplace’s equation:
∇2 φ = 0 (9.14)
P − p0 = 2 γ κ (9.15)
where p0 is the pressure in the gas (e.g. atmospheric pressure) and P is the pressure in the
liquid. κ is the mean curvature of the interface. According to geometry, if the surface is
given by r = f (x), then
1 1 + f ′2 − f f ′′
κ= .
2 f (1 + f ′2 )3/2
So, for the basic flow, f = a and κ = 1/(2a).
To assess the jet’s stability, we perturb the interface so that now
for some small parameter ε. We assume that we can expand the potential like this:
∂ 2 φ1 ∂ 2 φ1 1 ∂φ1 g′
+ + = g ′′
+ − k2 g = 0.
∂x2 ∂r 2 r ∂r r
160
So,
This is a modified form of Bessel’s equation. Recall that Bessel’s equation of zeroth order
looks like this:
r 2 g′′ + rg′ + r 2 g = 0.
Equation (9.16) is called the modified Bessel’s equation of zeroth order. It has linearly
independent solutions I0 (kr), K0 (kr). We can thus write the solution as
D
(r − f ) = 0 on r = f.
Dt
This states the the interface at r = f (x, t) must move along with the liquid. Expanding,
D Dr Df ∂f ∂f
(r − f ) = − = u− −w , (9.17)
Dt Dt Dt ∂t ∂x
on r = f . But u = ∇φ and so
∂φ ∂φ
u= = αkI0′ (kr) cos kx, w= = −αkI0 (kr) sin kx.
∂r ∂x
Therefore, on r = f , (9.17) becomes
h da1 i
ε αkI0′ (ka) − cos kx + O(ε2 ) = 0,
dt
giving
da1
= αkI0′ (ka).
dt
To find α, we need to consider the pressure jump at the surface, given by (9.15). To find
the pressure in the liquid we apply Bernoulli’s equation (see Hydrodynamics I, II) in its
unsteady version:
∂φ 1 p
+ |∇φ|2 + = c(t) (9.18)
∂t 2 ρ
161
Assuming we can write the pressure as
p = P + εp̂ cos kx
p − p0 = (P − p0 ) + εp̂(a, t) cos kx
h1 a1 i
2 2
= 2γ κ = γ −ε (1 − a k ) cos kx .
a a2
So P = p0 + γ/a and
a1
p̂(a, t) = −γ (1 − a2 k2 ).
a2
Combining this with the previous result,
kI ′ (ka)
0 d2 a1
− p̂(a, t) = ,
ρI0 (ka) dt2
d2 a1 γ k̂I0′ (k̂)
− a1 (1 − k̂2 ) = 0,
dt2 ρa3 I0 (k̂)
162
or
d2 a1 2
γ k̂I ′ (k̂)
2 − σ a1 = 0, with σ2 = 0
(1 − k̂2 )
dt ρa3 I0 (k̂)
This equation tells us how the amplitude a1 of the disturbance to the jet varies in time.
The solution is
a1 = c1 exp (σt) + c2 exp (−σt).
NOTE: If a suitable wavenumber destabilizes the jet, the amplitude of the disturbance
will keep on growing. Eventually it will grow so large that the above theory is invalidated.
Remember that we had to insist that ε is small. At this stage the disturbance interacts
nonlinearly with the basic flow and disrupts it drastically. Eventually, the jet pinches into
a sequence of satellite drops like this:
In section 2(b) we noted that it is possible to rotate from one set of axes to a second in
which the strain tensor is a diagonal matrix with non-zero entries along its diagonal and
zeros everywhere else. Under this action, the trace of the matrix is invariant. In this
appendix, we will see why this is the case.
First, we investigate rotations of a vector. Define a transformation matrix R such that its
action on a general vector x is to map it to another vector y of the same length. Accordingly,
yT y = xT x.
163
Furthermore,
yT y = (Rx)T (Rx) = xT RT Rx = xT x.
So,
xT (RT R − I)x = 0.
Since x is arbitrary, we deduce that
RT R = I.
In other words, R is an orthonal matrix whose transpose is equal to its inverse. Also, since
the length of the vector does not change, R represents a pure rotation of that vector through
some determined angle.
Now, if a matrix A has components a′ij written with respect to one Cartesian frame, and
A′ has components a′ij has written with respect to a second frame rotated with respect to
the first, it can be shown that
a′ij = Rki Rlj akl .
Or,
T
a′ij = Rik Rlj akl ,
and so,
A′ = RT AR.
164
Hence,
I1 = λ1 + λ2 + λ3
1
(λ1 + λ2 + λ3 )2 − (λ21 + λ22 + λ23 )
I2 = λ1 λ2 + λ2 λ3 + λ1 λ3 =
2
I3 = λ1 λ2 λ3 .
Or,
1
(trA)2 − tr(A2 ) ,
I1 = trA, I2 = I3 = detA.
2
So, since I1 = trA, the trace of a matrix A is invariant under rotation. Similarly, the
determinant of a matrix A is invariant under rotation, as is
1
(trA)2 − tr(A2 ) .
2
The momentum integral equation (3.30) is derived as follows. First, we state Newton’s
second law for a fixed volume of constant density fluid V (t),
∂u
Z Z Z Z
ρ dV + ρ u · ∇u dV = ρF dV + ∇ · σ dV .
V ∂t V V V
Note that8
u · ∇u = ∇ · (uu),
since ∇ · u = 0. Using the divergence theorem then, we have
∂u
Z Z Z Z
ρ dV + ρ u(u · n) dS = ρF dV + σ · n dS.
V ∂t S V S
Thus,
∂
Z Z Z Z
(ρ u) dV = − ρ u(u · n) dS + ρF dV + σ · n dS,
V ∂t S V S
as in (3.30).
8
Note that uu is an example of a dyadic product. In index notation, it means the matrix ui uj .
165