2Absolute and convective instabilities

III Hydrodynamic Stability



2 Absolute and convective instabilities
So far, we have only been considering perturbations that are Fourier modes,
(u
0
, p
0
, ρ
0
) = [
ˆ
u(z), ˆp(z), ˆρ(z)]e
i(k·xωt
.
This gives rise to a dispersion relation
D
(
k, ω
) = 0. This is an eigenvalue problem
for
ω
(
k
), and the
k
th mode is unstable if the solution
ω
(
k
) has positive imaginary
part.
When we focus on Fourier modes, they are necessarily non-local. In reality,
perturbations tend to be local. We perturb the fluid at some point, and the
perturbation spreads out. In this case, we might be interested in how the
perturbations spread.
To understand this, we need the notion of group velocity. For simplicity,
suppose we have a sum of two Fourier modes of slightly different wavenumbers
k
0
± k
. The corresponding frequencies are then
ω
0
± ω
, and for small
k
, we
may approximate
ω
= k
ω
k
k=k
0
.
We can then look at how our wave propagates:
η = cos[(k
0
+ k
)x (ω
0
+ ω
)t] + cos[k
0
k
)x (ω
0
ω
)t]
= 2 cos(k
x ω
t) cos(k
0
x ω
0
t)
= 2 cos
k
w
ω
k
k=k
0
t
!!
cos(k
0
x ω
0
t)
Since
k
is small, we know the first term has a long wavelength, and determines
the “overall shape” of the wave.
η
As time evolves, these “packets” propagate with group velocity
c
g
=
ω
k
. This
is also the speed at which the energy in the wave packets propagate.
In general, there are 4 key characteristics of interest for these waves:
The energy in a wave packet propagates at a group velocity.
They can disperse (different wavelengths have different speeds)
They can be advected by a streamwise flow.
They can be unstable, i.e. their amplitude can grow in time and space.
In general, if we have an unstable system, then we expect the perturbation to
grow with time, but also “move away” from the original source of perturbation.
We can consider two possibilities we can have convective instability, where
the perturbation “goes away”; and absolute instability, where the perturbation
doesn’t go away.
We can imagine the evolution of a convective instability as looking like this:
t
whereas an absolute instability would look like this:
t
Note that even in the absolute case, the perturbation may still have non-zero
group velocity. It’s just that the perturbations grow more quickly than the group
velocity.
To make this more precise, we consider the response of the system to an
impulse. We can understand the dispersion relation as saying that in phase space,
a quantity χ (e.g. velocity, pressure, density) must satisfy
D(k, ω) ˜χ(k, ω) = 0,
where
˜χ(k, ω) =
Z
−∞
Z
−∞
χ(x, t)e
i(kxωt)
dx dt.
Indeed, this equation says we can have a non-zero (
k, ω
) mode iff they satisfy
D
(
k, ω
) = 0. The point of writing it this way is that this is now linear in
˜χ
, and
so it applies to any χ, not necessarily a Fourier mode.
Going back to position space, we can think of this as saying
D
i
x
, i
t
χ(x, t) = 0.
This equation allows us to understand how the system responds to some external
forcing F (x, t). We simply replace the above equation by
D
i
x
, i
t
χ(x, t) = F (x, t).
Usually, we want to solve this in Fourier space, so that
D(k, ω) ˜χ(k, ω) =
˜
F (k, ω).
In particular, the Green’s function (or impulse response) is given by the response
to the impulse
F
ξ,τ
(
x, t
) =
δ
(
x ξ
)
δ
(
t τ
). We may wlog assume
ξ
=
τ
= 0,
and just call this
F
. The solution is, by definition, the Green’s function
G
(
x, t
),
satisfying
D
i
x
, i
t
G(x, t) = δ(x)δ(t).
Given the Green’s function, the response to an arbitrary forcing
F
is “just” given
by
χ(x, t) =
Z
G(x ξ, t τ)F (ξ, τ) dξ dτ.
Thus, the Green’s function essentially controls the all behaviour of the system.
With the Green’s function, we may now make some definitions.
Definition (Linear stability). The base flow of a system is linearly stable if
lim
t→∞
G(x, t) = 0
along all rays
x
t
= C.
A flow is unstable if it is not stable.
Definition
(Linearly convectively unstable)
.
An unstable flow is linearly con-
vectively unstable if lim
t→∞
G(x, t) = 0 along the ray
x
t
= 0.
Definition
(Linearly absolutely unstable)
.
An unstable flow is linearly absolutely
unstable if lim
t→∞
G(x, t) 6= 0 along the ray
x
t
= 0.
The first case is what is known as an amplifier, where the instability grows
but travels away at the same time. In the second case, we have an oscillator .
Even for a general
F
, it is easy to solve for
˜χ
in Fourier space. Indeed, we
simply have
˜χ(k, ω) =
˜
F (k, ω)
D(k, ω)
.
To recover χ from this, we use the Fourier inversion formula
χ(x, t) =
1
(2π)
2
Z
L
ω
Z
F
k
˜χ(k, ω)e
i(kxωt)
dk dω.
Note that here we put some general contours for
ω
and
k
instead of integrating
along the real axis, so this is not the genuine Fourier inversion formula. However,
we notice that as we deform our contour, when we pass through a singularity,
we pick up a multiple of
e
i(kxωt)
. Moreover, since the singularities occur when
D
(
k, ω
) = 0, it follows that these extra terms we pick up are in fact solutions to
the homogeneous equation
D
(
i∂
x
, i∂
t
)
χ
= 0. Thus, no matter which contour
we pick, we do get a genuine solution to our problem.
So how should we pick a contour? The key concept is causality the
response must come after the impulse. Thus, if
F
(
x, t
) = 0 for
t <
0, then we
also need
χ
(
x, t
) = 0 in that region. To understand this, we perform only the
temporal part of the Fourier inversion, so that
˜χ(k, t) =
1
2π
Z
L
ω
˜
F (k, ω)
D(k, ω; R)
e
t
dω.
To perform this contour integral, we close our contour either upwards or down-
wards, compute residues, and then take the limit as our semi-circle tends to
infinity. For this to work, it must be the case that the contribution by the
circular part of the contour vanishes in the limit, and this determines whether
we should close upwards or downwards.
If we close upwards, then
ω
will have positive imaginary part. So for
e
t
not to blow up, we need
t <
0. Similarly, we close downwards when
t >
0. Thus,
if we want
χ
to vanish whenever
t <
0, we should pick our contour so that it
it lies above all the singularities, so that it picks up no residue when we close
upwards. This determines the choice of contour.
But the key question is which of these are causal. When
t <
0, we require
χ
(
k, t
)
<
0. By Jordan’s lemma, for
t <
0, when performing the
ω
integral, we
should close the contour upwards when when perform the integral; when
t >
0,
we close it downwards. Thus for
χ
(
k, t
) = 0 when
t <
0, we must pick
L
ω
to lie
above all singularities of ˜χ(k, ω).
Assume that
D
has a single simple zero for each
k
. Let
ω
(
k
) be the corre-
sponding value of ω. Then by complex analysis, we obtain
˜χ(k, t) = i
˜
F [k, ω(k)]
˜
D
ω
[k, ω(k)]
e
(k)t
.
We finally want to take the inverse Fourier transform with respect to x:
χ(x, t) =
1
2π
Z
F
k
˜
F [k, ω(k)]
˜
D
ω
[k, ω(k)]
e
(k)t
dk.
We are interested in the case
F
=
δ
(
x
)
δ
(
t
), i.e.
˜
F
(
k, ω
) = 1. So the central
question is the evaluation of the integral
G(x, t) =
i
2π
Z
F
k
exp(i(kx ω(k)t)
˜
D
ω
[k, ω(k)]
dk.
Recall that our objective is to determine the behaviour of
G
as
t
with
V
=
x
t
fixed. Since we are interested in the large
t
behaviour instead of obtaining
exact values at finite
t
, we may use what is known as the method of steepest
descent.
The method of steepest descent is a very general technique to approximate
integrals of the form
H(t) =
i
2π
Z
F
k
f(k) exp ((k)) dk,
in the limit t . In our case, we take
f(k) =
1
˜
D
ω
[k, ω(k)]
ρ (k) = i (kV ω(k)) .
In an integral of this form, there are different factors that may affect the
limiting behaviour when
t
. First is that as
t
gets very large, the fast
oscillations in
ω
causes the integral to cancel except at stationary phase
ρ
i
k
= 0.
On the other hand, since we’re taking the exponential of
ρ
, we’d expect the
largest contribution to the integral to come from the
k
where
ρ
r
(
k
) is the greatest.
The idea of the method of steepest descent is to deform the contour so that
we integrate along paths of stationary phase, i.e. paths with constant
ρ
i
, so that
we don’t have to worry about the oscillating phase.
To do so, observe that the Cauchy–Riemann equations tell us
ρ
r
· ρ
i
= 0,
where in
we are taking the ordinary derivative with respect to
k
r
and
k
i
,
viewing the real and imaginary parts as separate variables.
Since the gradient of a function is the normal to the contours, this tells us
the curves with constant
ρ
i
are those that point in the direction of
ρ
r
. In other
words, the stationary phase curves are exactly the curves of steepest descent of
ρ
r
.
Often, the function
ρ
has some stationary points, i.e. points
k
satisfying
ρ
0
(
k
) = 0. Generically, we expect this to be a second-order zero, and thus
ρ
looks like (k k
)
2
near k
. We can plot the contours of ρ
i
near k
as follows:
k
where the arrows denote the direction of steepest descent of
ρ
r
. Note that
since the real part satisfies Laplace’s equation, such a stationary point must
be a saddle, instead of a local maximum or minimum, even if the zero is not
second-order.
We now see that if our start and end points lie on opposite sides of the ridge,
i.e. one is below the horizontal line and the other is above, then the only way to
do so while staying on a path of stationary phase is to go through the stationary
point.
Along such a contour, we would expect the greatest contribution to the
integral to occur when ρ
r
is the greatest, i.e. at k
. We can expand ρ about k
as
ρ(k) ρ(k
) +
1
2
2
ρ
k
2
(k
)(k k
)
2
.
We can then approximate
H(t)
i
2π
Z
ε
f(k)e
(k)
dk,
where we are just integrating over a tiny portion of our contour near
k
. Putting
in our series expansion of ρ, we can write this as
H(t)
i
2π
f(k
)e
(k
)
Z
ε
exp
t
2
2
ρ
k
2
(k
)(k k
)
2
dk.
Recall that we picked our path to be the path of steepest descent on both sides
of the ridge. So we can paramretrize our path by K such that
(iK)
2
=
t
2
2
ρ
k
2
(k
)(k k
)
2
,
where K is purely real. So our approximation becomes
H(t)
f(k
)e
(k
)
p
2π
2
00
(k
)
Z
ε
ε
e
K
2
dK.
Since
e
K
2
falls of so quickly as
k
gets away from 0, we may approximate that
integral by an integral over the whole real line, which we know gives
π
. So our
final approximation is
H(t)
f(k
)e
(k
)
p
2π
00
(k
)
.
We then look at the maxima of
ρ
r
(
k
) along these paths and see which has the
greatest contribution.
Now let’s apply this to our situation. Our ρ was given by
ρ(k) = i(kV ω(k)).
So k
is given by solving
ω
k
(k
) = V. ()
Thus, what we have shown was that the greatest contribution to the Green’s
function along the
V
direction comes from the the modes whose group velocity
is V ! Note that in general, this k
is complex.
Thus, the conclusion is that given any
V
, we should find the
k
such that (
)
is satisfied. The temporal growth rate along
x
t
= V is then
σ(V ) = ω
i
(k
) k
i
V.
This is the growth rate we would experience if we moved at velocity
V
. However,
it is often more useful to consider the absolute growth rate. In this case, we
should try to maximize
ω
i
. Suppose this is achieved at
k
=
k
max
, possibly
complex. Then we have
ω
i
k
(k
max
) = 0.
But this means that
c
g
=
ω
k
is purely real. Thus, this maximally unstable mode
can be realized along some physically existent V = c
g
.
We can now say
If ω
i,max
< 0, then the flow is linearly stable.
If ω
i,max
> 0, then the flow is linearly unstable. In this case,
If ω
0,i
< 0, then the flow is convectively unstable.
If ω
0,i
> 0, then the flow is absolutely unstable.
Here
ω
0
is defined by first solving
ω
k
(
k
0
) = 0, thus determining the
k
0
that
leads to a zero group velocity, and then setting
ω
0
=
ω
(
k
0
). These quantities
are known as the absolute frequency and absolute growth rate respectively.
Example.
We can consider a “model dispersion relation” given by the linear
complex Ginzburg–Landau equation
t
+ U
x
χ µχ (1 + ic
d
)
2
x
2
χ = 0.
This has advection (
U
), dispersion (
c
d
) and instability (
µ
). Indeed, if we replace
t
and
x
ik, then we have
i(ω + Uk)χ µχ + (1 + ic
d
)k
2
χ = 0.
This gives
ω = Uk + c
d
k
2
+ i(µ k
2
).
We see that we have temporal instability for
|k| <
µ
, where we assume
µ >
0.
This has
c
r
=
ω
r
k
= U + c
d
k.
On the other hand, if we force ω R, then we have spatial instability. Solving
k
2
(c
d
i) + Uk + ( ω) = 0,
the quadratic equation gives us two branches
k
±
=
U ±
p
k
2
4( ω)(c
d
i)
2(c
d
i)
.
To understand whether this instability is convective or absolute, we can complete
the square to obtain
ω = ω
0
+
1
2
ω
kk
(k k
0
)
2
,
where
ω
kk
= 2(c
d
i)
k
0
=
U
2(i c
d
)
ω
0
=
c
d
U
2
4(1 + c
2
d
)
+ i
µ
U
2
4(1 + c
2
d
)
.
These k
0
are then the absolute wavenumber and absolute frequency.
Note that after completing the square, it is clera that
k
0
is a double root at
ω
=
ω
0
. Of course, there is not nothing special about this example, since
k
0
was
defined to solve
ω
k
(k
0
) = 0!
I really ought to say something about Bers’ method at this point.
In practice, it is not feasible to apply a
δ
function of perturbation and see
how it grows. Instead, we can try to understand convective versus absolute
instabilities using a periodic and switched on forcing
F (x, t) = δ(x)H(t)e
f
t
,
where H is the Heaviside step function. The response in spectral space is then
˜χ(k, ω) =
˜
F (k, ω)
D(k, ω)
=
i
D(k, ω)(ω ω
f
)
.
There is a new (simple) pole precisely at
ω
=
ω
f
on the real axis. We can invert
for t to obtain
˜χ(k, t) =
i
2π
Z
L
ω
e
t
D(k, ω)(ω ω
F
)
dω =
e
f
t
D(k, ω
f
)
+
e
(k)t
(ω(k) ω
f
)
˜
D
ω
[k, ω(k)]
.
We can then invert for x to obtain
χ(x, t) =
e
f
t
2π
Z
F
k
e
ikx
D(k, ω)
dk
| {z }
χ
F
(x,t)
+
1
2π
Z
F
k
e
i[kxω(k)t]
[ω(k) ω
f
]
˜
D
ω
[k, ω(k)]
| {z }
χ
T
(x,t)
.
The second term is associated with the switch-on transients, and is very similar
to what we get in the previous steepest descent.
If the flow is absolutely unstable, then the transients dominate and invade
the entire domain. But if the flow is convectively unstable, then this term is
swept away, and all that is left is
χ
F
(
x, t
). This gives us a way to distinguish
between
Note that in the forcing term, we get contributions at the singularities
k
±
(
ω
f
).
Which one contributes depends on causality. One then checks that the correct
result is
χ
F
(x, t) = iH(x)
e
i(k
+
(ω
f
)xω
f
t)
D
k
[k
+
(ω
f
), ω
f
]
iH(x)
e
i(k
(ω
f
)xω
f
t)
D
k
[k
(ω
f
), ω
f
]
.
Note that
k
±
(
ω
f
) may have imaginary parts! If there is some
ω
f
such that
k
+
i
(
ω
f
)
>
0, then we will see spatially growing waves in
x >
0. Similarly, if
there exists some
ω
f
such that
k
i
(
ω
f
)
<
0, then we see spatially growing
waves in x < 0.
Note that we will see this effect only when we are convectively unstable.
Perhaps it is wise to apply these ideas to an actual fluid dynamics problem.
We revisit the broken line shear layer profile, scaled with the velocity jump, but
allow a non-zero mean U :
z = 1
z = 1
¯
U = 1 + U
m
¯
U = z + U
m
¯
U = 1 + U
m
III
II
I
Ae
α(z1)
Be
αz
+ Ce
αz
De
α(z+1)
We do the same interface matching conditions, and after doing the same compu-
tations (or waving your hands with Galilean transforms), we get the dispersion
relation
4(ω U
m
α)
2
= (2α 1)
2
e
4α
.
It is now more natural to scale with
U
m
rather than
U/
2, and this involves
expressing everything in terms of the velocity ratio
R
=
U
2U
m
. Then we can write
the dispersion relation as
D(k, ω; R) = 4(ω k)
2
R
2
[(2k 1)
2
e
4k
] = 0.
Note that under this scaling, the velocity for
z <
1 is
u
= 1
R
. In particular,
if
R <
1, then all of the fluid is flowing in the same direction, and we might
expect the flow to “carry away” the perturbations, and this is indeed true.
The absolute/convective boundary is given by the frequency at wavenumber
for zero group velocity:
ω
k
(k
0
) = 0.
This gives us
ω
0
= k
0
R
2
2
[2k
0
1 + e
4k
0
].
Plugging this into the dispersion relations, we obtain
R
2
[2k
0
1 + e
4k
0
] [(2k
0
1)
2
e
4k
0
] = 0.
We solve for
k
0
, which is in general complex, and then solve for
ω
to see if it
leads to
ω
0,i
>
0. This has to be done numerically, and when we do this, we see
that the convective/absolute boundary occurs precisely at R = 1.
Gaster relation
Temporal and spatial instabilities are related close to a marginally stable state.
This is flow at a critical parameter
R
c
with critical (real) wavenumber and
frequency: D(k
c
, ω
c
; R
c
) = 0 with ω
c,i
= k
c,i
= 0.
We can Taylor expand the dispersion relation
ω = ω
c
+
ω
k
(k
c
; R
c
)[k k
c
].
We take the imaginary part
ω
i
=
ω
i
k
r
(k
c
, R
c
)(k
r
k
c
) +
ω
r
k
r
(k
c
, R
c
)k
i
.
For the temporal mode, we have k
i
= 0, and so
ω
(T )
i
=
ω
i
k
r
(k
c
, R
c
)(k
r
k
c
).
For the spatial mode, we have ω
i
= 0, and so
0 =
ω
i
k
r
(k
c
, R
c
)[k
r
k
c
] +
ω
r
k
r
(k
c
, R
c
)k
(S)
i
.
Remembering that c
g
=
ω
r
k
r
. So we find that
ω
(T )
i
= c
g
k
(S)
i
.
This gives us a relation between the growth rates of the temporal mode and
spatial mode when we are near the marginal stable state.
Bizarrely, this is often a good approximation when we are far from the
marginal state.
Global instabilities
So far, we have always been looking at flows that were parallel, i.e. the base
flow depends on
z
alone. However, in real life, flows tend to evolve downstream.
Thus, we want to consider base flows U = U(x, z).
Let
λ
be the characteristic wavelength of the perturbation, and
L
be the
characteristic length of the change in
U
. If we assume
ε
λ
L
1, then we may
want to perform some local analysis.
To leading order in
ε
, evolution is governed by frozen dispersion relation at
each
X
=
εx
. We can then extend notions of stability/convective/absolute to
local notions, e.g. a flow is locally convectively unstable if there is some
X
such
that ω
i,max
(X) > 0, but ω
0,i
(X) < 0 for all X.
However, we can certainly imagine some complicated interactions between
the different region. For example, a perturbation upstream may be swept
downstream by the flow, and then get “stuck” somewhere down there. In general,
we can define
Definition
(Global stability)
.
A flow is globally stable if
lim
t→∞
G
(
x, t
) = 0 for
all x.
A flow is globally unstable if there is some
x
such that
lim
t→∞
G
(
x, t
)
.
For a steady spatially developing flow, global modes are
χ(x, t) = φ(x)e
G
t
.
The complex global frequency
ω
G
is determined analogously to before using
complex integration.
It can be established that
ω
G,i
ω
0,i,max
. This then gives a necessary
condition for global instability: there must be a region of local absolute instability
within the flow.
Sometimes this is a good predictor, i.e.
R
G
c
' R
t
. For example, with mixing
layers, we have
R
t
= 1
.
315 while
R
G
c
= 1
.
34. On the other hand, it is sometimes
poor. For example, for bluff-body wakes, we have Re
t
= 25 while Re
Gc
= 48.5.