1Linear stability analysis
III Hydrodynamic Stability
1.1 Rayleigh–Taylor instability
In this section, we would like to investigate the stability of a surface between
two fluids:
η
ρ
2
ρ
1
We assume the two fluids have densities
ρ
1
and
ρ
2
, and are separated by a
smooth interface parametrized by the deviation η.
Recall that the Navier–Stokes equations for an incompressible fluid are
ρ
∂u
∂t
+ u · ∇u
= −∇P − gρ
ˆ
z + µ∇
2
u, ∇ · u = 0.
Usually, when doing fluid mechanics, we assume density is constant, and divide
the whole equation across by
ρ
. We can then forget the existence of
ρ
and carry
on. We can also get rid of the gravity term. We know that the force of gravity
is balanced out by a hydrostatic pressure gradient. More precisely, we can write
P = P
h
(z) + p(x, t),
where P
h
satisfies
∂P
h
∂z
= −gρ.
We can then write our equation as
∂u
∂t
+ u · ∇u = −∇
p
ρ
+ ν∇
2
u.
To do all of these, we must assume that the density
ρ
is constant, but this is
clearly not the case when we have tow distinct fluids.
Since it is rather difficult to move forward without making any simplifying
assumptions, we shall focus on the case of inviscid fluids, so that we take
ν
= 0.
We will also assume that the fluid is irrotational. In this case, it is convenient to
use the vector calculus identity
u · ∇u = ∇
1
2
|u|
2
− u × (∇ × u),
where we can now drop the second term. Moreover, since
∇ × u
= 0, Stokes’
theorem allows us to write
u
=
∇φ
for some velocity potential
φ
. Finally, in
each separate region, the density ρ is constant. So we can now write
∇
ρ
∂φ
∂t
+
1
2
ρ|∇φ|
2
+ P + gρz
= 0.
This is Bernoulli’s theorem.
Applying these to our scenario, since we have two fluids, we have two separate
velocity potentials
φ
1
, φ
2
for the two two regions. Both of these independently
satisfy the incompressibility hypothesis
∇
2
φ
1,2
= 0.
Bernoulli’s theorem tells us the quantity
ρ
∂φ
∂t
+
1
2
ρ|∇φ|
2
+ P + gρz
should be constant across the fluid. Our investigation will not use the full
equation. All we will need is that this quantity matches at the interface, so that
we have
ρ
1
∂φ
1
∂t
+
1
2
ρ
1
|∇φ
1
|
2
+ P
1
+ gρ
1
η
z=η
= ρ
2
∂φ
2
∂t
+
2
2
ρ
2
|∇φ
2
|
2
+ P
2
+ gρ
2
η
z=η
.
To understand the interface, we must impose boundary conditions. First of all
the vertical velocities of the fluids must match with the interface, so we impose
the kinematic boundary condition
∂φ
1
∂z
z=η
=
∂φ
2
∂z
z=η
=
Dη
Dt
,
where
D =
∂
∂t
+ u · ∇
We also make another, perhaps dubious boundary condition, namely that the
pressure is continuous along the interface. This does not hold at all interfaces.
For example, if we have a balloon, then the pressure inside is greater than the
pressure outside, since that is what holds the balloon up. In this case, it is the
rubber that is exerting a force to maintain the pressure difference. In our case,
we only have two fluids meeting, and there is no reason to assume a discontinuity
in pressure, and thus we shall assume it is continuous. If you are not convinced,
this is good, and we shall later see this is indeed a dubious assumption.
But assuming it for the moment, this allows us to simplify our Bernoulli’s
theorem to give
ρ
1
∂φ
1
∂t
+
1
2
ρ
1
|∇φ
1
|
2
+ gρ
1
η = ρ
2
∂φ
2
∂t
+
2
2
ρ
2
|∇φ
2
|
2
+ gρ
2
η.
There is a final boundary condition that is specific to our model. What we are
going to do is that we will start with a flat, static solution with
u
= 0 and
η
= 0.
We then perturb
η
a bit, and see what happens to our fluid. In particular, we
want to see if the perturbation is stable.
Since we expect the interesting behaviour to occur only near the interface,
we make the assumption that there is no velocity in the far field, i.e.
φ
1
→
0
as
z → ∞
and
φ
2
→
0 as
z → −∞
(and similarly for the derivatives). For
simplicity, we also assume there is no y dependence.
We now have equations and boundary conditions, so we can solve them. But
these equations are pretty nasty and non-linear. At this point, one sensible
approach is to linearize the equations by assuming everything is small. In addition
to assuming that
η
is small, we also need to assume that various derivatives such
as
∇φ
are small, so that we can drop all second-order terms. Since
η
is small,
the value of, say,
∂φ
1
∂z
at
η
should be similar to that at
η
= 0. Since
∂φ
1
∂z
itself is
already assumed to be small, the difference would be second order in smallness.
So we replace all evaluations at
η
with evaluations at
z
= 0. We are then left
with the collection of 3 equations
∇
2
φ
1,2
= 0
∂φ
1,2
∂z
z=0
= η
t
ρ
1
∂φ
1
∂t
+ gρ
1
η
z=0
= ρ
2
∂φ
2
∂t
+ gρ
2
η
z=0
.
This is a nice linear problem, and we can analyze the Fourier modes of the
solutions. We plug in an ansatz
φ
1,2
(x, z, t) =
ˆ
φ
1,2
(z)e
i(kx−ωt)
η(x, t) = Be
i(kx−ωt)
.
Substituting into Laplace’s equation gives
ˆ
φ
1,2
− k
2
ˆ
φ
1,2
= 0.
Using the far field boundary conditions, we see that we have a family of solutions
ˆ
φ
1
= A
1
e
−kz
,
ˆ
φ
2
= A
2
e
kz
.
The kinematic boundary condition tells us we must have
ˆ
φ
0
1,2
(0) = −iωB.
We can solve this to get
B =
kA
1
iω
= −
kA
2
iω
.
In particular, we must have A ≡ A
1
= −A
2
. We can then write
η =
kA
iω
e
ik(kx−ωt)
.
Plugging these into the final equation gives us
ρ
1
(−iωA) + gρ
1
kA
iω
= ρ
2
iωA + gρ
2
kA
iω
.
Crucially, we can cancel the
A
throughout the equation, and gives us a result
independent of
A
. This is, after all, the point of linearization. Solving this gives
us a dispersion relation relating the frequency (or phase speed
c
p
=
ω/k
) to the
wavenumber:
ω
2
=
g(ρ
2
− ρ
1
)k
ρ
1
+ ρ
2
.
If
ρ
2
ρ
1
, then this reduces to
ω
2
≈ gk
, and this is the usual dispersion relation
for deep water waves. On the other hand, if
ρ
2
− ρ
1
>
0 is small, this can lead
to waves of rather low frequency, which is something we can observe in cocktails,
apparently.
But nothing in our analysis actually required
ρ
2
> ρ
1
. Suppose we had
ρ
1
> ρ
2
. This amounts to putting the heavier fluid on top of the lighter one.
Anyone who has tried to turn a cup of water upside down will know this is highly
unstable. Can we see this from our analysis?
If
ρ
1
< ρ
2
, then
ω
has to be imaginary. We shall write it as
ω
=
±iσ
, where
σ > 0. We can then compute
σ =
s
g(ρ
1
− ρ
2
)k
ρ
1
+ ρ
2
.
Writing out our expression for
φ
1,2
, we see that there are
e
±σt
terms, and the
e
σt
term dominates in the long run, causing
φ
1,2
to grow exponentially. This is
the Rayleigh–Taylor instability.
There is more to say about the instability. As the wavelength decreases,
k
increases, and we see that
σ
increases as well. Thus, we see that short-scale
perturbations blow up exponentially much more quickly, which means the system
is not very well-posed. This is the ultra-violet catastrophe. Of course, we should
not trust our model here. Recall that in our simplifying assumptions, we not
only assumed that the amplitudes of our
φ
and
η
were small, but also that their
derivatives were small. The large
k
behaviour gives us large
x
derivatives, and
so we have to take into account the higher order terms as well.
But we can also provide physical reasons for why small scales perturbations
should be suppressed by such a system. In our model, we assumed there is no
surface tension. Surface tension quantifies the tendency of interfaces between
fluids to minimize surface area. We know that small-scale variations will cause
a large surface area, and so the surface tension will suppress these variations.
Mathematically, what the surface allows for is a pressure jump across the
interface.
Surface tension is quantified by
γ
, the force per unit length. This has
dimensions [
γ
] =
MT
−1
. Empirically, we observe that the pressure difference
across the interface is
∆p = −γ∇ ·
ˆ
n = 2γH = γ
1
R
x
+
1
R
y
,
where
ˆ
n
is the unit normal and
H
is the mean curvature. This is an empirical
result.
Again, we assume no dependence in the
y
direction, so we have a cylindrical
interface with a single radius of curvature. Linearizing, we have a pressure
difference of
(P
2
− P
1
)|
z=η
=
γ
R
x
= −γ
∂
2
η
∂x
2
1 +
∂η
∂x
2
3/2
≈ −γ
∂
2
η
∂x
2
.
Therefore the (linearized) dynamic boundary condition becomes
ρ
1
∂φ
1
∂t
+ gρ
1
η
z=0
+ γ
∂
2
η
∂x
2
= ρ
2
∂φ
2
∂t
+ gρ
2
η
z=0
.
If we run the previous analysis again, we find a dispersion relation of
ω
2
=
g(ρ
2
− ρ
1
)k + γk
3
ρ
1
+ ρ
2
.
Since
γ
is always positive, even if we are in the situation when
ρ
1
> ρ
2
, for
k
sufficiently large, the system is stable. Of course, for small
k
, the system is still
unstable — we still expect water to fall out even with surface tension. In the
case ρ
1
< ρ
2
, what we get is known as internal gravity-capillary waves.
In the unstable case, we have
σ
2
= k
g(ρ
2
− ρ
1
)
ρ
1
+ ρ
2
(1 − l
2
c
k
2
),
where
l
2
c
=
γ
g(ρ
1
− ρ
2
)
is a characteristic length scale. For
k`
c
>
1, the oscillations are stable, and the
maximum value of k is attained at k =
l
c
√
3
.
In summary, we have a range of wavelength, and we have identified a “most
unstable wavelength”. Over time, if all frequencies modes are triggered, we
should expect this most unstable wavelength to dominate. But it is rather
hopeful thinking, because this is the result of linear analysis, and we can’t expect
it to carry on to the non-linear regime.