1Linear stability analysis

III Hydrodynamic Stability



1.5 Stratified flows
In the Kelvin–Helmholtz scenario, we had a varying density. Let’s now try to
model the situation in complete generality.
Recall that the set of (inviscid) equations for Boussinesq stratified fluids is
ρ
u
t
+ u · u
= −∇p
0
gρ
0
ˆ
z,
with incompressibility equations
· u = 0,
Dρ
Dt
= 0.
We consider a mean base flow and a 2D perturbation as a normal mode:
u =
¯
U(z)
ˆ
x + u
0
(x, z, t)
p = ¯p(z) + p
0
(x, z, t)
ρ = ¯ρ(z) + ρ
0
(x, z, t),
with
(u
0
, p
0
, ρ
0
) = [
ˆ
u(z), ˆp(z), ˆρ(z)]e
i(kxωt)
= [
ˆ
u(z), ˆp(z), ˆρ(z)]e
ik(xct)
.
We wish to obtain an equation that involves only
ˆw
. We first linearize our
equations to obtain
¯ρ
u
0
t
+
¯
U
u
0
x
+ w
0
d
¯
U
dz
=
p
0
x
¯ρ
w
0
t
+
¯
U
w
0
x
=
p
0
z
gρ
0
ρ
0
t
+
¯
U
ρ
0
x
+ w
0
d¯ρ
dz
= 0
u
0
x
+
w
0
z
= 0.
Plugging in our normal mode solution into the four equations, we obtain
ik¯ρ(
¯
U c)ˆu + ¯ρ ˆw
d
dz
¯
U = ikˆp
ik¯ρ(
¯
U c) ˆw =
dˆp
dz
g ˆρ. ()
ik(
¯
U c)ˆρ + w
0
d¯ρ
dz
= 0 ()
ikˆu +
w
0
z
= 0.
The last (incompressibility) equation helps us eliminate
ˆu
from the first equation
to obtain
¯ρ(
¯
U c)
d ˆw
dz
+ ¯ρ ˆw
d
¯
U
dz
= ikˆp.
To eliminate
ˆp
as well, we differentiate with respect to
z
and apply (
), and then
further apply () to get rid of the ˆρ term, and end up with
¯ρ(
¯
U c)
d
2
ˆw
dz
2
+ ¯ρ ˆw
d
2
¯
U
dz
2
+ k
2
¯ρ(
¯
U c) ˆw =
g ˆw
¯
U c
d¯ρ
dz
.
This allows us to write our equation as the Taylor–Goldstein equation
d
2
dz
2
k
2
ˆw
ˆw
(
¯
U c)
d
2
¯
U
dz
2
+
N
2
ˆw
(
¯
U c)
2
= 0,
where
N
2
=
g
¯ρ
d¯ρ
dz
.
This
N
is known as the buoyancy frequency, and has dimensions
T
2
. This is
the frequency at which a slab of stratified fluid would oscillate vertically.
Indeed, if we have a slab of volume
V
, and we displace it by an infinitesimal
amount ζ in the vertical direction, then the buoyancy force is
F = V gζ
ρ
z
.
Since the mass of the fluid is
ρ
0
V
, by Newton’s second law, the acceleration of
the parcel satisfies
d
2
ζ
dt
2
+
g
ρ
0
ρ
z
ζ = 0.
This is just simple harmonic oscillation of frequency N .
Crucially, what we see from this computation is that a density stratification
of the fluid can lead to internal waves. We will later see that the presence of
these internal waves can destabilize our system.
Miles–Howard theorem
In 1961, John Miles published a series of papers establishing a sufficient condition
for infinitesimal stability in the case of a stratified fluid. When he submitted
this to review, Louis Howard read the paper, and replied with a 3-page proof of
a more general result. Howard’s proof was published as a follow-up paper titled
“Note on a paper of John W. Miles”.
Recall that we just derived the Taylor–Goldstein equation
d
2
dz
2
k
2
ˆw
ˆw
(
¯
U c)
d
2
¯
U
dz
2
+
N
2
ˆw
(
¯
U c)
2
= 0,
The magical insight of Howard was to introduce the new variable
H =
ˆ
W
(
¯
U c)
1/2
.
We can regretfully compute the derivatives
ˆw = (
¯
U c)
1/2
H
d
dz
ˆw =
1
2
(
¯
U c)
1/2
H
d
¯
U
dz
+ (
¯
U c)
1/2
dH
dz
d
2
dz
2
ˆw =
1
4
(
¯
U c)
3/2
H
d
¯
U
dz
2
+
1
2
(
¯
U c)
1/2
H
d
2
¯
U
dz
2
+ (
¯
U c)
1/2
dH
dz
d
¯
U
dz
+ (
¯
U c)
1/2
d
2
H
dz
2
.
We can substitute this into the Taylor–Goldstein equation, and after some
algebraic mess, we obtain the decent-looking equation
d
dz
(
¯
U c)
dH
dz
H
k
2
(
¯
U c) +
1
2
d
2
¯
U
dz
2
+
1
4
d
¯
U
dz
2
N
2
¯
U c
= 0.
This is now self-adjoint. Imposing boundary conditions
ˆw,
d ˆw
dz
0 in the far
field, we multiply by H
and integrate over all space. The first term gives
Z
H
d
dz
(
¯
U c)
dH
dz
dz =
Z
(
¯
U c)
dH
dz
2
dz,
while the second term is given by
Z
k
2
|H|
2
(
¯
U c)
1
2
d
2
¯
U
dz
2
|H|
2
|H|
2
1
4
d
¯
U
dz
2
N
2
(
¯
U c
)
|
¯
U c|
2
dz.
Both the real and imaginary parts of the sum of these two terms must be zero.
In particular, the imaginary part reads
c
i
Z
dH
dz
2
+ k
2
|H|
2
+ |H|
2
N
2
1
4
d
¯
U
dz
2
|
¯
U c|
2
dz = 0.
So a necessary condition for instability is that
N
2
1
4
d
¯
U
dz
2
< 0
somewhere. Defining the Richardson number to be
Ri(z) =
N
2
(d
¯
U/dz)
2
,
the necessary condition is
Ri(z) <
1
4
.
Equivalently, a sufficient condition for stability is that Ri(z)
1
4
everywhere.
How can we think about this? When we move a parcel, there is the buoyancy
force that attempts to move it back to its original position. However, if we move
the parcel to the high velocity region, we can gain kinetic energy from the mean
flow. Thus, if the velocity gradient is sufficiently high, it becomes “advantageous”
for parcels to move around.
Piecewise-linear profiles
If we want to make our lives simpler, we can restrict our profile with piecewise
linear velocity and layered density. In this case, to solve the Taylor–Goldstein
equation, we observe that away from the interfaces, we have
N
=
U
00
= 0. So
the Taylor–Goldstein equation reduces to
d
2
dz
2
k
2
ˆw = 0.
This has trivial exponential solutions in each of the individual regions, and to
find the correct solutions, we need to impose some matching conditions at the
interface.
So assume that the pressure is continuous across all interfaces. Then using
the equation
¯ρ(
¯
U c)
d ˆw
dz
+ ρ ˆw
d
¯
U
dz
= ikˆp,
we see that the left hand seems like the derivative of a product, except the sign
is wrong. So we divide by
1
(
¯
Uc)
2
, and then we can write this as
ik
¯ρ
ˆp
(
¯
U c)
2
=
d
dz
ˆw
¯
U c
.
For a general
c
, we integrate over an infinitesimal distance at the interface. We
assume
ˆp
is continuous, and that
¯ρ
and (
¯
U c
) have bounded discontinuity.
Then the integral of the LHS vanishes in the limit, and so the integral of the
right-hand side must be zero. This gives the matching condition
ˆw
¯
U c
+
= 0.
To obtain a second matching condition, we rewrite the Taylor–Goldstein equa-
tion as a derivative, which allows for the determination of the other matching
condition:
d
dz
(
¯
U c)
d ˆw
dz
ˆw
d
¯
U
dz
g ¯ρ
ρ
0
ˆw
¯
U c

= k
2
(
¯
U c) ˆw
g ¯ρ
ρ
0
d
dz
ˆw
¯
U c
.
Again integrating over an infinitesimal distance over at the interface, we see
that the LHS must be continuous along the interface. So we obtain the second
matching condition.
(
¯
U c)
d ˆw
dz
ˆw
d
¯
U
dz
g ¯ρ
ρ
0
ˆw
¯
U c

+
= 0.
We begin by applying this to a relatively simple profile with constant density,
and whose velocity profile looks like
h
U
For convenience, we scale distances by
h
2
and speeds by
U
2
, and define
˜c
and α by
c =
U
2
˜c, α =
kh
2
.
These quantities
α
and
˜c
(which we will shortly start writing as
c
instead) are
nice dimensionless quantities to work with. After the scaling, the profile can be
described by
z = 1
z = 1
¯
U = 1
¯
U = z
¯
U = 1
III
II
I
ˆw = Ae
α(z1)
ˆw = Be
αz
+ Ce
αz
ˆw = De
α(z+1)
We have also our exponential solutions in each region between the interfaces,
noting that we require our solution to vanish at the far field.
We now apply the matching conditions. Since
¯
U c
is continuous, the first
matching condition just says ˆw has to be continuous. So we get
A = Be
α
+ Ce
α
D = Be
α
+ Ce
α
.
The other matching condition is slightly messier. It says
(
¯
U c)
d ˆw
dz
ˆw
d
¯
U
dz
g ¯ρ
ρ
0
ˆw
¯
U c

+
= 0,
which gives us two equations
(Be
α
+ Ce
α
)(1 c)(α) = (Be
α
Ce
α
)(1 c)(α) (Be
α
+ Ce
α
)
(Be
α
+ Ce
α
)(1 c)(α) = (Be
α
Ce
α
)(1 c)(α) (Be
α
+ Ce
α
).
Simplifying these gives us
(2α(1 c) 1)Be
α
= Ce
α
(2α(1 + c) 1)Ce
α
= Be
α
.
Thus, we find that
(2α 1)
2
4α
2
c
2
= e
4α
,
and hence we have the dispersion relation
c
2
=
(2α 1)
2
e
4α
4α
2
.
We see that this has the possibility of instability. Indeed, we can expand the
numerator to get
c
2
=
(1 4α + 4α
2
) (1 4α + 8α
2
+ O(α
3
))
4α
2
= 1 + O(α).
So for small
α
, we have instability. This is known as Rayleigh instability. On
the other hand, as
α
grows very large, this is stable. We can plot these out in a
graph
k
ω
ω
i
ω
r
0.64
We see that the critical value of
α
is 0
.
64, and the maximally instability point
is α 0.4.
Let’s try to understand where this picture came from. For large
k
, the
wavelength of the oscillations is small, and so it seems reasonable that we can
approximate the model by one where the two interfaces don’t interact. So
consider the case where we only have one interface.
z = 1
¯
U = z
¯
U = 1
II
I
Ae
α(z1)
Be
α(z1)
We can perform the same procedure as before, solving
ˆw
¯
U c
+
= 0.
This gives the condition A = B. The other condition then tell us
(1 c)(α)A = (1 c)αA Ae
α(z1)
.
We can cancel the A, and rearrange this to say
c = 1
1
2α
.
So we see that this interface supports a wave at speed
c
= 1
1
2α
at a speed
lower than
¯
U
. Similarly, if we treated the other interface at isolation, it would
support a wave at c
= 1 +
1
2α
.
Crucially these two waves are individually stable, and when
k
is large, so that
the two interfaces are effectively in isolation, we simply expect two independent
and stable waves. Indeed, in the dispersion relation above, we can drop the
e
4α
term when α is large, and approximate
c = ±
2α 1
2α
= ±
1
1
2α
.
When
k
is small, then we expect the two modes to interact with each other,
which manifests itself as the
e
4α
term. Moreover, the resonance should be the
greatest when the two modes have equal speed, i.e. at
α
1
2
, which is quite
close to actual maximally unstable mode.
Density stratification
We next consider the case where we have density stratification. After scaling,
we can draw our region and solutions as
z = 1
z = 0
z = 1
¯
U = 1
¯
U = z
¯
U = 1
III
IV
II
I
Ae
α(z1)
Be
α(z1)
+ Ce
α(z1)
De
α(z+1)
+ Ee
α(z+1)
F e
α(z+1)
¯ρ = 1
¯ρ = +1
Note that the
¯ρ
is the relative density, since it is the density difference that
matters (fluids cannot have negative density).
Let’s first try to understand how this would behave heuristically. As before,
at the I-II and III-IV interfaces, we have waves with
c
=
±
1
1
2α
. At the
density interface, we previously saw that we can have internal gravity waves
with c
igw
= ±
q
Ri
0
α
, where
Ri
0
=
gh
ρ
0
is the bulk Richardson number (before scaling, it is Ri
0
=
gρh
ρ
0
U
2
).
We expect instability to occur when the frequency of this internal gravity
wave aligns with the frequency of the waves at the velocity interfaces, and this
is given by
1
1
2α
2
=
Ri
0
α
.
It is easy to solve this to get the condition
Ri
0
' α 1.
This is in fact what we find if we solve it exactly.
So let’s try to understand the situation more responsibly now. The procedure
is roughly the same. The requirement that ˆw is continuous gives
A = B + C
Be
α
+ Ce
α
= De
α
+ Ee
α
D + E = F.
If we work things out, then the other matching condition
(
¯
U c)
d ˆw
dz
ˆw
d
¯
U
dz
Ri
0
¯ρ
ˆw
¯
U c

+
= 0
gives
(1 c)(α)A = (1 c)(α)(B C) (B + C)
(1 c)(α)F = (1 c)(α)(D E) (D + E)
(c)(α)(Be
α
Ce
α
) +
Ri
0
(c)
(Be
α
+ Ce
α
)
= (c)(α)(De
α
Ee
α
)
Ri
0
(c)
(De
α
+ Ee
α
).
This defines a 6
×
6 matrix problem, 4th order in
c
. Writing it as
SX
= 0 with
X
= (
A, B, C, D, E, F
)
T
, the existence of non-trivial solutions is given by the
requirement that det S = 0, which one can check (!) is given by
c
4
+ c
2
e
4α
(2α 1)
2
4α
2
Ri
0
α
+
Ri
0
α
e
2α
+ (2α 1)
2α
2
= 0.
This is a biquadratic equation, which we can solve.
We can inspect the case when
Ri
0
is very small. In this case, the dispersion
relation becomes
c
4
+ c
2
e
4α
(2α 1)
2
4α
2
.
Two of these solutions are simply given by
c
= 0, and the other two are just those
we saw from Rayleigh instability. This is saying that when the density difference
is very small, then the internal gravitational waves are roughly non-existent, and
we can happily ignore them.
In general, we have instability for all
Ri
0
! This is Holmboe instability. The
scenario is particularly complicated if we look at small
α
, where we expect the
effects of Rayleigh instability to kick in as well. So let’s fix some 0
< α <
0
.
64,
and see what happens when we increase the Richardson number.
As before, we use dashed lines to denote the imaginary parts and solid lines
to denote the real part. For each
Ri
0
, any combination of imaginary part and
real part gives a solution for
c
, and, except when there are double roots, there
should be four such possible combinations.
Ri
0
c
We can undersatnd this as follows. When
Ri
0
= 0, all we have is Rayleigh
instability, which are the top and bottom curves. There are then two
c
= 0
solutions, as we have seen. As we increase
Ri
0
to a small, non-sero amount, this
mode becomes non-zero, and turns into a genuine unstable mode. As the value
of
Ri
0
increases, the two imaginary curves meet and merge to form a new curve.
The solutions then start to have a non-zero real part, which gives a non-zero
phase speed to our pertrubation.
Note that when
Ri
0
becomes very large, then our modes become unstable
again, though this is not clear from our graph.
Taylor instability
While Holmboe instability is quite interesting, our system was unstable even
without the density stratification. Since we saw that instability is in general trig-
gered by the interaction between two interfaces, it should not be surprising that
even in a Rayleigh stable scenario, if we have two layers of density stratification,
then we can still get instability.
Consider the following flow profile:
z = 1
z = 1
¯
U = z
III
II
I
Ae
α(z1)
Be
αz
+ Ce
αz
De
α(z+1)
¯ρ = R 1
¯ρ = R + 1
As before, continuity in ˆw requires
A = Be
α
+ Ce
α
D = Be
α
+ Ce
α
.
Now there is no discontinuity in vorticity, but the density field has jumps. So
we need to solve the second matching condition. They are
(1 c)(α)(Be
α
+ Ce
α
) +
Ri
0
1 c
(Be
α
+ Ce
α
) = (1 c)(α)(Be
α
Ce
α
)
(1 c)(α)(Be
α
+ Ce
α
) +
Ri
0
1 + c
(Be
α
+ Ce
α
) = (1 c)(α)(Be
α
Ce
α
)
These give us, respectively,
(2α(1 c)
2
Ri
0
)B = Ri
0
Ce
2α
(2α(1 + c)
2
Ri
0
)C = Ri
0
Be
2α
.
So we get the biquadratic equation
c
4
c
2
2 +
Ri
0
α
+
1
Ri
0
2α
2
Ri
2
0
e
4α
4α
2
= 0.
We then apply the quadratic formula to say
c
2
= 1 +
Ri
0
2α
±
r
2Ri
0
α
+
Ri
2
0
e
4α
4α
2
.
So it is possible to have instability with no inflection! Indeed, instability occurs
when c
2
< 0, which is equivalent to requiring
2α
1 + e
2α
< Ri
0
<
2α
1 e
2α
.
Thus, for any
Ri
0
>
0, we can have instability. This is known as Taylor
instability.
Heuristically, the density jumps at the interface have speed
c
±
igw
= ±1
Ri
0
2α
1/2
.
We get equality when
c
= 0, so that
Ri
0
= 2
α
. This is in very good agreement
with what our above analysis gives us.