1Linear stability analysis
III Hydrodynamic Stability
1.2 Rayleigh–B´enard convection
The next system we want to study is something that looks like this:
T
0
T
0
+ ∆TT
0
+ ∆T
d
There is a hot plate at the bottom, a cold plate on top, and some fluid in between.
We would naturally expect heat to transmit from the bottom to the top. There
are two ways this can happen:
–
Conduction: Heat simply diffuses from the bottom to the top, without any
fluid motion.
–
Convection: The fluid at the bottom heats up, expands, becomes lighter,
and moves to the top.
The first factor is controlled purely by thermal diffusivity
κ
, while the latter
is mostly controlled by the viscosity
ν
. It is natural to expect that when the
temperature gradient ∆
T
is small, most of the heat transfer is due to conduction,
as there isn’t enough force to overcome the viscosity to move fluid around. When
∆
T
is large, the heat transfer will be mostly due to conduction, and we would
expect such a system to be unstable.
To understand the system mathematically, we must honestly deal with the
case where we have density variations throughout the fluid. Again, recall that
the Navier–Stokes equation
ρ
∂u
∂t
+ u · ∇u
= −∇P − gρ
ˆ
z + µ∇
2
u, ∇ · u = 0.
The static solution to our system is the pure conduction solution, where
u
= 0,
and there is a uniform vertical temperature gradient across the fluid. Since the
density is a function of temperature, which in the static case is just a function
of
z
, we know
ρ
=
ρ
(
z
). When we perturb the system, we allow some horizontal
and time-dependent fluctuations, and assume we can decompose our density
field into
ρ = ρ
h
(z) + ρ
0
(x, t).
We can then define the hydrostatic pressure by the equation
dP
h
dz
= −gρ
h
.
This then allows us to decompose the pressure as
P = P
h
(z) + p
0
(x, t).
Then our equations of motion become
∂u
∂t
+ u · ∇u = −
1
ρ
∇p
0
−
gρ
0
ρ
ˆ
z + ν∇
2
u, ∇ · u = 0
We have effectively “integrated out” the hydrostatic part, and it is now the
deviation from the average that leads to buoyancy forces.
An important component of our analysis involves looking at the vorticity.
Indeed, vorticity necessarily arises if we have non-trivial density variations. For
concreteness, suppose we have a interface between two fluids with different
densities, ρ and ρ + ∆ρ. If the interface is horizontal, then nothing happens:
ρ
ρ + ∆ρ
However, if the interface is tilted, then there is a naturally induced torque:
ρ
ρ + ∆ρ
In other words, vorticity is induced. If we think about it, we see that the reason
this happens is that the direction of the density gradient does not align with the
direction of gravity. More precisely, this is because the density gradient does not
align with the pressure gradient.
Let’s try to see this from our equations. Recall that the vorticity is defined
by
ω = ∇ × u.
Taking the curl of the Navier–Stokes equation, and doing some vector calculus,
we obtain the equation.
∂ω
∂t
+ u · ∇ω = ω · ∇u +
1
ρ
2
∇ρ × ∇P + ν∇
2
ω
The term on the left hand side is just the material derivative of the vorticity. So
we should interpret the terms on the right-hand-side as the terms that contribute
to the change in vorticity.
The first and last terms on the right are familiar terms that have nothing to
do with the density gradient. The interesting term is the second one. This is
what we just described — whenever the pressure gradient and density gradient
do not align, we have a baroclinic torque.
In general, equations are hard to solve, and we want to make approximations.
A common approximation is the Boussinesq approximation. The idea is that
even though the density difference is often what drives the motion, from the
point of view of inertia, the variation in density is often small enough to be
ignored. To take some actual, physical examples, salt water is roughly 4% more
dense than fresh water, and every 10 degrees Celsius changes the density of air
by approximately 4%.
Thus, what we do is that we assume that the density is constant except in
the buoyancy force. The mathematically inclined people could think of this as
taking the limit g → ∞ but ρ
0
→ 0 with gρ
0
remaining finite.
Under this approximation, we can write our equations as
∂u
∂t
+ u · ∇u = −
1
ρ
0
∇p
0
− g
0
ˆ
z + ν∇
2
u
∂ω
∂t
+ u · ∇ω = ω · ∇u +
g
ρ
0
ˆ
z × ∇ρ + ν∇
2
ω,
where we define the reduced gravity to be
g
0
=
gρ
0
ρ
0
Recall that our density is to be given by a function of temperature
T
. We
must write down how we want the two to be related. In reality, this relation is
extremely complicated, and may be affected by other external factors such as
salinity (in the sea). However, we will use a “leading-order” approximation
ρ = ρ
0
(1 − α(T − T
0
)).
We will also need to know how temperature
T
evolves in time. This is simply
given by the diffusion equation
∂T
∂t
+ u · ∇T = κ∇
2
T.
Note that in thermodynamics, temperature and pressure are related quantities.
In our model, we will completely forget about this. The only purpose of pressure
will be to act as a non-local Lagrange multiplier that enforces incompressibility.
There are some subtleties with this approximation. Inverting the relation,
T =
1
α
1 −
ρ
ρ
0
+ T
0
.
Plugging this into the diffusion equation, we obtain
∂ρ
∂t
+ u · ∇ρ = κ∇
2
ρ.
The left-hand side is the material derivative of density. So this equation is saying
that density, or mass, can “diffuse” along the fluid, independent of fluid motion.
Mathematically, this follows from the fact that density is directly related to
temperature, and temperature diffuses.
This seems to contradict the conservation of mass. Recall that the conserva-
tion of mass says
∂ρ
∂t
+ ∇ · (uρ) = 0.
We can than expand this to say
∂ρ
∂t
+ u · ∇ρ = −ρ∇ · u.
We just saw that the left-hand side is given by
κ∇
2
ρ
, which can certainly be
non-zero, but the right-hand side should vanish by incompressibility! The issue
is that here the change in
ρ
is not due to compressing the fluid, but simply
thermal diffusion.
If we are not careful, we may run into inconsistencies if we require
∇ · u
= 0.
For our purposes, we will not worry about this too much, as we will not use the
conservation of mass.
We can now return to our original problem. We have a fluid of viscosity
ν
and thermal diffusivity
κ
. There are two plates a distance
d
apart, with the
top held at temperature
T
0
and bottom at
T
0
+ ∆
T
. We make the Boussinesq
approximation.
T
0
T
0
+ ∆TT
0
+ ∆T
d
We first solve for our steady state, which is simply
U = 0
T
h
= T
0
− ∆T
(z − d)
d
ρ
h
= ρ
0
1 + α∆T
z − d
d
P
h
= P
0
− gρ
0
z +
α∆T z
2d
(z − 2d)
,
In this state, the only mode of heat transfer is conduction. What we would
like to investigate, of course, is whether this state is stable. Consider small
perturbations
u = U + u
0
T = T
h
+ θ
P = P
h
+ p
0
.
We substitute these into the Navier–Stokes equation. Starting with
∂u
∂t
+ u · ∇u = −
1
ρ
0
∇p
0
+ gαθ
ˆ
z + ν∇
2
u,
we assume the term u · ∇u will be small, and end up with the equation
∂u
0
∂t
= −
1
ρ
0
∇p
0
+ gαθ
ˆ
z + ν∇
2
u,
together with the incompressibility condition ∇ · u
0
= 0.
Similarly, plugging these into the temperature diffusion equation and writing
u = (u, v, w) gives us
∂θ
∂t
− w
0
∆T
d
= κ∇
2
θ.
We can gain further understanding of these equations by expressing them in
terms of dimensionless quantities. We introduce new variables
˜
t =
κt
d
2
˜
x =
x
d
˜
θ =
θ
∆T
˜p =
d
2
p
0
ρ
0
κ
2
In terms of these new variables, our equations of motion become rather elegant:
∂
˜
θ
∂
˜
t
− ˜w =
˜
∇
2
˜
θ
∂
˜
u
∂
˜
t
= −
˜
∇˜p +
gα∆T d
3
νκ
ν
κ
˜
θ
ˆ
z +
ν
κ
˜
∇
2
˜
u
Ultimately, we see that the equations of motion depend on two dimensionless
constants: the Rayleigh number and Prandtl number
Ra =
gα∆T d
3
νκ
Pr =
ν
κ
,
These are the two parameters control the behaviour of the system. In particular,
the Prandtl number measures exactly the competition between viscous and
diffusion forces. Different fluids have different Prandtl numbers:
–
In a gas, then
ν
κ
∼
1, since both are affected by the mean free path of a
particle.
–
In a non-metallic liquid, the diffusion of heat is usually much slower than
that of momentum, so Pr ∼ 10.
–
In a liquid metal,
Pr
is very very low, since, as we know, metal transmits
heat quite well.
Finally, the incompressibility equation still reads
˜
∇ ·
˜
u = 0
Under what situations do we expect an unstable flow? Let’s start with some
rather more heuristic analysis. Suppose we try to overturn our system as shown
in the diagram:
Suppose we do this at a speed
U
. For each
d × d × d
block, the potential
energy released is
P E ∼ g · d · (d
3
∆ρ) ∼ gρ
0
αd
4
∆T,
and the time scale for this to happen is τ ∼ d/U .
What is the friction we have to overcome in order to achieve this? The
viscous stress is approximately
µ
∂U
∂z
∼ µ
U
d
.
The force is the integral over the area, which is
∼
µU
d
d
2
. The work done, being
R
F dz, then evaluates to
µU
d
· d
2
· d = µUd
2
.
We can figure out
U
by noting that the heat is diffused away on a time scale
of
τ ∼ d
2
/κ
. For convection to happen, it must occur on a time scale at least
as fast as that of diffusion. So we have
U ∼ κ/d
. All in all, for convection to
happen, we need the potential energy gain to be greater than the work done,
which amounts to
gρ
0
α∆T d
4
& ρ
0
νκd.
In other words, we have
Ra =
gα∆T d
3
νκ
=
g
0
d
3
νκ
& 1.
Now it is extremely unlikely that
Ra ≥
1 is indeed the correct condition, as
our analysis was not very precise. However, it is generally true that instability
occurs only when Ra is large.
To work this out, let’s take the curl of the Navier–Stokes equation we
previously obtained, and dropping the tildes, we get
∂ω
∂t
= RaPr∇θ ×
ˆ
z + Pr∇
2
ω
This gets rid of the pressure term, so that’s one less thing to worry about. But
this is quite messy, so we take the curl yet again to obtain
∂
∂t
∇
2
u = RaPr
∇
2
θ
ˆ
z − ∇
∂θ
∂z
+ Pr∇
4
u
Reading off the z component, we get
∂
∂t
∇
2
w = RaPr
∂
2
∂x
2
+
∂
2
∂y
2
θ + Pr∇
4
w.
We will combine this with the temperature equation
∂θ
∂t
− w = ∇
2
θ
to understand the stability problem.
We first try to figure out what boundary conditions to impose. We’ve got a
6th order partial differential equation (4 from the Navier–Stokes and 2 from the
temperature equation), and so we need 6 boundary conditions. It is reasonable
to impose that w = θ = 0 at z = 0, 1, as this gives us impermeable boundaries.
It is also convenient to assume the top and bottom surfaces are stress-free:
∂u
∂z
=
∂v
∂z
= 0, which implies
∂
∂z
∂u
∂x
+
∂v
∂y
= 0.
This then gives us
∂
2
w
∂z
2
= 0. This is does not make much physical sense, because
our fluid is viscous, and we should expect no-slip boundary conditions, not
no-stress. But they are mathematically convenient, and we shall assume them.
Let’s try to look for unstable solutions using separation of variables. We see
that the equations are symmetric in x and y, and so we write our solution as
w = W (z)X(x, y)e
σt
θ = Θ(z)X(x, y)e
σt
.
If we can find solutions of this form with σ > 0, then the system is unstable.
What follows is just a classic separation of variables. Plugging these into the
temperature equation yields
d
2
dz
2
− σ +
∂
2
∂x
2
+
∂
2
∂y
2
XΘ = −XW.
Or equivalently
d
2
dz
2
− σ
XΘ + XW = −
∂
2
∂x
2
+
∂
2
∂y
2
XΘ.
We see that the differential operator on the left only acts in Θ, while those on
the right only act on X. We divide both sides by XΘ to get
d
2
dy
2
− σ
Θ + W
Θ
= −
d
2
dx
2
+
d
2
dy
2
X
X
.
Since the LHS is purely a function of
z
, and the RHS is purely a function of
X
,
they must be constant, and by requiring that the solution is bounded, we see
that this constant must be positive. Our equations then become
∂
2
∂x
2
+
∂
2
∂y
2
X = −λ
2
X
d
2
dz
2
− λ
2
− σ
Θ = −W.
We now look at our other equation. Plugging in our expressions for
w
and
θ
,
and using what we just obtained, we have
σ
d
2
dz
2
− λ
2
W = −RaPrλ
2
Θ + Pr
d
2
dz
2
− λ
2
2
W.
On the boundary, we have Θ = 0 =
W
=
∂
2
W
∂z
2
= 0, by assumption. So it follows
that we must have
d
4
W
dz
4
z=0,1
= 0.
We can eliminate Θ by letting
d
2
dz
2
− σ − λ
2
act on both sides, which converts
the Θ into W . We are then left with a 6th order differential equation
d
2
dz
2
− σ − λ
2
Pr
d
2
dz
2
− λ
2
2
− σ
d
2
dz
2
− λ
2
!
W = −RaPrλ
2
W.
This is an eigenvalue problem, and we see that our operator is self-adjoint. We
can factorize and divide across by Pr to obtain
d
2
dz
2
− λ
2
d
2
dz
2
− σ − λ
2
d
2
dz
2
− λ
2
−
σ
Pr
W = −Raλ
2
W.
The boundary conditions are that
W =
d
2
W
dz
2
=
d
4
W
dx
4
= 0 at z = 0, 1.
We see that any eigenfunction of
d
2
dz
2
gives us an eigenfunction of our big scary
differential operator, and our solutions are quantized sine functions,
W
=
sin nπz
with n ∈ N. In this case, the eigenvalue equation is
(n
2
π
2
+ λ
2
)(n
2
π
2
+ σ + λ
2
)
n
2
π
2
+ λ
2
+
σ
Pr
= Raλ
2
.
When
σ >
0, then, noting that the LHS is increasing with
σ
on the range [0
, ∞
),
this tells us that we have
Ra(n) ≥
(n
2
π
2
+ λ
2
)
3
λ
2
.
We minimize the RHS with respect to
λ
and
n
. We clearly want to set
n
= 1,
and then solve for
0 =
d
d[λ
2
]
(n
2
π
2
+ λ
2
)
3
λ
2
=
(π
2
+ λ
2
)
2
λ
4
(2λ
2
− π
2
).
So the minimum value of Ra is obtained when
λ
c
=
π
√
2
.
So we find that we have a solution only when
Ra ≥ Ra
c
=
27π
4
4
≈ 657.5.
If we are slightly about
Ra
c
, then there is the
n
= 1 unstable mode. As
Ra
increases, the number of available modes increases. While the critical Rayleigh
number depends on boundary conditions, but the picture is generic, and the
width of the unstable range grows like
√
Ra − Ra
c
:
kd
Ra
stable
unstable
We might ask — how does convection enhance the heat flux? This is
quantified by the Nusselt number
Nu =
convective heat transfer
conductive heat transfer
=
hwT id
κ∆T
= h˜w
˜
θi.
Since this is a non-dimensional number, we know it is a function of Ra and Pr.
There are some physical arguments we can make about the value of
Nu
. The
claim is that the convective heat transfer is independent of layer depth. If there
is no convection at all, then the temperature gradient will look approximately
like
There is some simple arguments we can make about this. If we have a purely
conductive state, then we expect the temperature gradient to be linear:
This is pretty boring. If there is convection, then we claim that the temperature
gradient will look like
The reasoning is that the way convection works is that once a parcel at the
bottom is heated to a sufficiently high temperature, it gets kicked over to the
other side; if a parcel at the top is cooled to a sufficiently low temperature, it gets
kicked over to the other side. In between the two boundaries, the temperature is
roughly constant, taking the average temperature of the two boundaries.
In this model, it doesn’t matter how far apart the two boundaries are, since
once the hot or cold parcels are shot off, they can just travel on their own until
they reach the other side (or mix with the fluid in the middle).
If we believe this, then we must have
hwT i ∝ d
0
, and hence
Nu ∝ d
. Since
only Ra depends on d, we must have Nu ∝ Ra
1/3
.
Of course, the real situation is much more subtle.