4A variational point of view
III Hydrodynamic Stability
4 A variational point of view
In this chapter, we are going to learn about a more robust and powerful way
to approach and understand transient growth. Instead of trying to think about
G(t) as a function of t, let us fix some target time T , and just look at G(T ).
For any two times t
i
< t
f
, we can define a propagator function such that
u(t
f
) = Φ(t
f
, t
i
)u(t
i
)
for any solution
u
of our equations. In the linear approximation, this propagator
is a linear function between appropriate function spaces. Writing Φ = Φ(
T, T
0
),
the gain problem is then equivalent to maximizing
G(T, T
0
) =
E(T )
E(T
0
)
=
hu
p
(T ), u
p
(T )i
hu
p
(T
0
), u
p
(T
0
)i
=
hΦu
p
(T
0
), Φu
p
(T
0
)i
hu
p
(T
0
), u
p
(T
0
)i
=
hu
p
(T
0
), Φ
†
Φu
p
(T
0
)i
hu
p
(T
0
), u
p
(T
0
)i
.
Here the angled brackets denote natural inner product leading to the energy
norm. Note that the operator Φ
†
Φ is necessarily self-adjoint, and so this
G
is maximized when
u
p
(
T
0
) is chosen to be the eigenvector of Φ of maximum
eigenvalue.
There is a general method to find the maximum eigenvector of a (self-adjoint)
operator Φ. We start with a random vector
x
. Then we have Φ
n
x ∼ λ
n
1
v
1
as
n → ∞
, where
λ
1
is the maximum eigenvalue with associated eigenvector
v
1
.
Indeed, if we write
x
as the linear combination of eigenvectors, then as we apply
Φ many times, the sum is dominated by the term with largest eigenvalue.
So if we want to find the mode with maximal transient growth, we only need
to be able to compute Φ
†
Φ. The forward propagator Φ(
T, T
0
) is something we
know how to compute (at least numerically). We simply numerically integrate
the Navier–Stokes equation. So we need to understand Φ(T, T
0
)
†
.
Let u(t) be a solution to the linear equation
du
dt
(t) = L(t)u(t).
Here we may allow
L
to depend on time. Let
L
†
be the adjoint, and suppose
v(t) satisfies
dv
dt
(t) = −L(t)
†
v(t).
Then the chain rule tells us
d
dt
hv(t), u(t)i = h−L(t)
†
v(t), u(t)i + hv(t), L(t)
†
u(t)i = 0.
So we know that
hv(T
0
), u(T
0
)i = hv(T ), u(T )i = hv(T ), Φ(T, T
0
)u(T
0
)i.
Thus, given a
v
0
, to compute Φ(
T, T
0
)
†
v
0
, we have to integrate the adjoint
dv
dt
=
−L
(
t
)
†
v
backwards in time, using the “initial” condition
v
(
T
) =
v
0
, and
then we have
Φ(T, T
0
)
†
v
0
= v(T
0
).
What does the adjoint equation look like? For a time-dependent background
shear flow, the linearized forward/direct equation for a perturbation
u
p
is given
by
∂u
p
∂t
+ (U(t) · ∇)u
p
= −∇p
p
− (u
p
· ∇)u(t) + Re
−1
∇
2
u
p
.
The adjoint (linearized) Navier–Stokes equation is then
−
∂u
a
∂t
= Ω × u
a
− ∇ × (U × u
σ
) − ∇p
a
+ Re
−1
∇
2
u
a
,
where we again have
∇ · u
a
= 0, Ω = ∇ × U.
This PDE is ill-posed if we wanted to integrate it forwards in time, but that
does not concern us, because to compute Φ
†
, we have to integrate it backwards
in time.
Thus, to find the maximal transient mode, we have to run the direct-adjoint
loop.
u
p
(T
0
) u
p
(T )
Φ
†
Φ(u
p
(T
0
)) u
p
(T )
Φ
Φ
†
,
and keep running this until it converges.
Using these analysis, we can find some interesting results. For example, in
a shear layer flow, 3-dimensional modes with both streamwise and spanwise
perturbations are the most unstable. However, in the long run, the Kelvin-
Helmholtz modes dominate.
Variational formulation
We can also use variational calculus to find the maximally unstable mode. We
think of the calculation as a constrained optimization problem with the following
requirements:
(i) For all T
0
≤ t ≤ T , we have
∂q
∂t
= D
t
q = Lq.
(ii) The initial state is given by q(0) = q
0
.
We will need Lagrangian multipliers to impose these constraints, and so the
augmented Lagrangian is
G =
hq
T
, q
T
i
hq
0
, q
0
i
−
Z
T
0
h
˜
q, (D
t
− L)qi dt + h
˜
q
0
, q(0) − q
0
i.
Taking variations with respect to
˜
q
and
˜
q
0
recover the evolution equation and
the initial condition. The integral can then be written as
Z
T
0
h
˜
q, (D
t
− L)qi dt =
Z
T
0
hq, (D
t
+ L
†
)
˜
qi dt + h
˜
q
0
, q
0
i − h
˜
q
T
, q
T
i.
Now if we take a variation with respect to q, we see that
˜
q has to satisfy
(D
t
+ L
†
)
˜
q = 0.
So the Lagrange multiplier within such a variational problem evolves according
to the adjoint equation!
Taking appropriate adjoints, we can write
G =
hq
T
, q
T
i
hq
0
, q
0
i
+
Z
T
0
hq, (D
t
+ L
†
)
˜
qi dt + h
˜
q
0
, q
0
i−h
˜
q
T
, q
T
i+ boundary terms.
But if we take variations with respect to our initial conditions, then
δG
δq
0
= 0
gives
˜
q
0
=
2hq
T
, q
T
i
hq
0
, q
0
i
2
q
0
.
Similarly, setting
δG
δq
T
= 0, we get
˜
q
T
=
2
hq
0
, q
0
i
q
T
.
Applying h·, q
0
i to the first equation and h·, q
T
i to the second, we find that
h
˜
q
0
, q
0
i = h
˜
q
T
, q
T
i.
This is the equation we previously saw for the adjoint equation, and this provides
a consistency condition to see if we have found the optimal solution.
Previously our algorithm used power iteration, which requires us to integrate
forwards and backwards many times. However, we now have gradient information
for gain:
δG
∂q
0
=
˜
q
0
−
2hq
T
, q
T
i
hq
0
, q
0
i
2
q
0
.
This allows us to exploit optimization algorithms (steepest descent, conjugate
gradient etc.). This has an opportunity for faster convergence.
There are many ways we can modify this. One way is to modify the inner
product. Note that we actually had to use the inner product in three occasions:
(i) Inner product in the objective functional J
(ii) Inner product in initial normalization/constraint of the state vector.
(iii) Inner product in the definition of adjoint operator.
We used the same energy norm for all of these, but there is no reason we have
to use the same norm for all of these. However, there are strong arguments that
an energy norm is natural for norm 2. It is a wide open research question as to
whether there is an appropriate choice for inner product 3. On the other hand,
investigation of variation of inner product 1 has been widely explored.
As an example, consider p-norms
J =
1
V
Ω
Z
Ω
e(x, T )
p
dΩ
1/p
, e(x, T ) =
1
2
|u(x, T )|
2
.
This has the attraction that for large values of
p
, this would be dominated by
peak values of e, not average.
When we set
p
= 1, i.e. we use the usual energy norm, then we get a beautiful
example of the Orr mechanism, in perfect agreement with what SVD tells us.
For, say, p = 50, we get more exotic center/wall modes.
Non-linear adjoints
In the variational formulation, there is nothing in the world that stops us from
using a non-linear evolution equation! We shall see that this results in slightly
less pleasant formulas, but it can still be done.
Consider plane Couette flow with
Re
= 1000, and allow arbitrary amplitude
in perturbation
u
tot
= U + u, ∂
t
u + (u + U) · ∇(u + U) = −∇p + Re
−1
∇
2
u,
where U = yˆx. We can define a variational problem with Lagrangian
L =
E(T )
E
0
−[∂
t
u+N(u)+∇p, v]−[∇·u, q]−
1
2
hu
0
, u
0
i − E
0
c−hu(0)−u
0
, v
0
i,
where
N(u
i
) = U
j
∂
j
u
i
+ u
i
∂
i
U
j
+ u
j
∂
j
u
i
−
1
Re
∂
j
∂
j
u
i
, [v, u] =
1
T
Z
T
0
hv, ui dt.
Variations with respect to the direct variable
u
can once again define a non-linear
adjoint equation
δL
δu
= ∂
t
v + N
†
(v, u) + ∇q +
u
E
0
− v
t=T
+ (v −v
0
)|
t=0
= 0,
where
N
†
(v
i
, u) = ∂
j
(u
j
, v
i
) − v
j
∂
i
u
j
+ ∂
j
(U
j
v
i
) − v
j
∂
i
U
j
− v
j
∂
i
U
j
+
1
Re
∂
j
∂
j
v
i
.
We also have
δL
δp
= ∇ · v = 0,
δL
δu
0
= v
0
− cu
0
= 0.
Note that the equation for the adjoint variable
v
depends on the direct variable
u
, but is linear in
v
. Computationally, this means that we need to remember
our solution to
u
when we do our adjoint loop. If
T
is large, then it may not be
feasible to store all information about
u
across the whole period, as that is too
much data. Instead, the method of checkpointing is used:
(i) Pick “checkpoints” 0 = T
0
< ··· < T
K
= T .
(ii) When integrating u, we remember high resolution data for u(x, T
K
).
(iii)
When integrating
v
backwards in the interval (
T
K−1
, T
K
), we use the data
remembered at
T
K−1
to re-integrate to obtain detailed information about
u in the interval (T
k−1
, T
K
).
This powerful algorithmic approach allows us to identify minimal seed of turbu-
lence, i.e. the minimal, finite perturbation required to enter a turbulence mode.
Note that this is something that cannot be understood by linear approximations!
This is rather useful, because in real life, it allows us to figure out how to modify
our system to reduce the chance of turbulence arising.