Olvera Cabrera, D., & Kerswell, R. (2017). Optimising energy growth
as a tool for finding exact coherent structures.
Physical Review Fluids
,
2
(8), [083902]. https://doi.org/10.1103/PhysRevFluids.2.083902
Publisher's PDF, also known as Version of record
Link to published version (if available):
10.1103/PhysRevFluids.2.083902
Link to publication record in Explore Bristol Research
PDF-document
This is the final published version of the article (version of record). It first appeared online via APS Physics at
https://journals.aps.org/prfluids/abstract/10.1103/PhysRevFluids.2.083902. Please refer to any applicable terms
of use of the publisher.
University of Bristol - Explore Bristol Research
General rights
This document is made available in accordance with publisher policies. Please cite only the
published version using the reference above. Full terms of use are available:
http://www.bristol.ac.uk/red/research-policy/pure/user-guides/ebr-terms/
PHYSICAL REVIEW FLUIDS 2, 083902 (2017)
Optimizing energy growth as a tool for finding exact coherent structures
D. Olvera
*
and R. R. Kerswell
†
School of Mathematics, Bristol University, Bristol BS8 1TW, United Kingdom
(Received 30 January 2017; published 16 August 2017)
We discuss how searching for finite-amplitude disturbances of a given energy that
maximize their subsequent energy growth after a certain later time T can be used to probe
the phase space around a reference state and ultimately to find other nearby solutions.
The procedure relies on the fact that of all the initial disturbances on a constant-energy
hypersphere, the optimization procedure will naturally select the one that lies closest
to the stable manifold of a nearby solution in phase space if T is large enough. Then,
when in its subsequent evolution the optimal disturbance transiently approaches the new
solution, a flow state at this point can be used as an initial guess to converge the solution
to machine precision. We illustrate this approach in plane Couette flow by rediscovering
the spanwise-localized “snake” solutions of Schneider et al. [Phys. Rev. Lett. 104, 104501
(2010)], probing phase space at very low Reynolds numbers (less than 127.7) where
the constant-shear solution is believed to be the global attractor and examining how the
edge between laminar and turbulent fl ow evolves when stable stratification eliminates
the turbulence. We also show that the steady snake solution smoothly delocalizes as
unstable s tratification is gradually turned on until it connects (via an intermediary global
three-dimensional solution) to two-dimensional Rayleigh-Bénard roll solutions.
DOI: 10.1103/PhysRevFluids.2.083902
I. INTRODUCTION
Optimization has proved a powerful tool to extract information from the Navier-Stokes equations.
In the shear flow transition problem, optimizing over all possible infinitesimal disturbances to
find the one that maximizes the subsequent energy growth after some pre-selected time T has
proven invaluable in exposing the generic energy amplification mechanisms present. Called variously
transient growth [1,2], nonmodal instability [3], or optimal perturbation theory [4] (see the reviews
in [3,5]aswellas[6]), the approach reveals key aspects of the linearized dynamics around the
reference state that have helped to interpret finite-time flow phenomena and pick apart what causes
a transition. The approach owes its popularity to its linearity, which means that there are multiple
ways to extract the optimals and the mathematics in each case is well understood (see, e.g., [3,7–9]).
The downside of the approach is that it can say nothing about finite-amplitude disturbances or, in
other words, what can happen a finite distance away from the reference state in phase space [10,11].
Conceptually, the remedy to this is simple: Let competing disturbances seeking to maximize the
energy growth after time T all have the same initial finite energy E
0
and use the fully nonlinear
Navier-Stokes equations as a constraint [12–14]. This, however, doubles the number of parameters
(E
0
joins T ) over which the results must be interpreted and leads to a fully nonlinear, nonconvex
optimization problem where much less is known about its possibly multiple solutions (local as well
as global maxima) or how to find them. So far, the solution technique has necessarily been iterative
and this has revealed a number of interesting new insights into the transition problem [12,14–22]
(see the review in [13]). For example, one can ask what the smallest (most “dangerous”) energy
disturbance is that can trigger a transition by some time T , with the answer in the large-T limit
*
do12542@bris.ac.uk
†
R.R.Kerswell@bris.ac.uk
2469-990X/2017/2(8)/083902(22) 083902-1 ©2017 American Physical Society
D. OLVERA AND R. R. KERSWELL
−2 0 2 4 6 8 10 12
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
X
2
X
1
0 2 4 6 8 10
0
2
4
6
8
10
12
|X|
t
FIG. 1. Sample trajectories initiated close to the stable manifold (X
2
= 1) of the saddle point in the 2D
model of [13]: dX
1
/dt =−X
1
+ 10X
2
and dX
2
/dt = X
2
(10e
−X
2
1
/100
− X
2
)(X
2
− 1) (left plot). The closer
the initial conditions, the more a plateau in the norm |X| (here a proxy for the energy) where X := (X
1
,X
2
)
develops near the value for the saddle point (marked as a dashed black line at
√
101) during the subsequent
evolution (see the right plot). The initial conditions are (X
1
,X
2
) = (0,1 − ), where = 10
−2
(red line), 10
−5
(magenta line), 10
−8
(cyan line), and 10
−11
(blue line).
labeled the minimal seed for a transition [20–22]. The minimal seeds that emerge from this procedure
are fully localized and are therefore realistic targets for experimental investigations (see, e.g., [21]).
The optimization approach works by naturally selecting disturbances on the energy hypersphere
if they lie outside the basin of attraction of the reference state since then the energy remains finite for
T →∞(all other disturbances have to decay eventually). The new state to which these disturbances
are drawn need not be a turbulent attractor and so the minimal disturbance to reach another simple
stable state can also be calculated (see Sec. 6.2 in [23]). What is not so clear is whether the approach
can find unstable solutions, although a similar line of reasoning seems to hold. The optimization
procedure would be expected to select disturbances from the energy hypersphere that lie nearest or
on the stable manifold of a nearby solution in phase space if T is large enough as this is the best
way to avoid energy decay. The difference now, however, is since the nearby solution is unstable,
T cannot be too large; otherwise even these disturbances will have decayed away (realistically, it is
improbable to stay on the stable manifold to converge in to the unstable state). Once such an optimal
disturbance has been found, its temporal evolution will show evidence of a transient approach to the
alternative solution. Such an approach can be identified by noticing that the energy of the optimal
disturbance plateaus near the energy of the different solution for some period of time. Figure 1 uses
the simple two-dimensional (2D) model of [ 13] to illustrate this: The subsequent evolution of an
initial condition taken closer and closer to the stable manifold of a nearby saddle point (actually
the edge state in this model system) spends longer and longer close to the saddle before ultimately
being repelled to the base state. A sufficiently close visit should yield flow states that can then be
converged to the solution. The main purpose of this paper is to demonstrate that this approach can
work.
We illustrate this in the context of plane Couette flow (PCF) by rediscovering the spanwise-
localized “snake” solutions of Schneider et al. [24] building upon the prior exploratory work of
Rabin (see Sec. 6.3 in [23]), who identified the key role played by the choice of E
0
.Inawide
geometry, snake s olutions coexist with repeated copies of Nagata’s well known solution [25]ina
narrow geometry and, not surprisingly, the stable manifolds of the (lower-energy) snake solutions
pass closer (in energy norm) to the constant-shear solution in phase space than those of Nagata’s
solution. However, the latter offer the possibility of greater energy growth as they lead to a global
flow state and hence are preferred by the optimization algorithm if they pass close to the energy
083902-2
OPTIMIZING ENERGY GROWTH AS A TOOL FOR . . .
hypersphere. As a result, Rabin found a threshold initial energy below which the optimal disturbance
appears to approach a snake solution and above which Nagata’s solution is approached (Sec. 6.3
[23]). We complete this calculation here by recomputing these optimal disturbances and converging
both snake solutions.
Armed with this success, we then probe the phase space of PCF at very low Reynolds numbers
(less than 127.7[26]) looking for alternative solutions where the (basic) constant shear solution is
believed to be the global attractor (a proof only exists for Re < 20.7[27]). We find evidence of
solutions, but these turn out only to be t he ghosts of known solutions at higher Reynolds numbers.
Finally, as another example of how the optimization approach can be utilized, we examine how the
edge between laminar and turbulent flow evolves when stable stratification s uppresses the turbulence.
With turbulence present, a bursting phenomenon sometimes forms a distinctive initial feature of the
transition process for disturbances above the edge. This bursting is found to change little when
the turbulent attractor vanishes under increasing stratification but disappears when a certain ECS
ceases to exist. This then indicates that the bursting is directly related to the presence of the unstable
manifold of this state directed away from the constant-shear solution in phase space.
The paper starts with the formulation of the stratified plane Couette flow problem in Sec. II A,
which introduces the three nondimensional parameters that fully specify the problem once the
computational box is chosen: the Reynolds number Re, the bulk Richardson number Ri
b
, and the
Prandtl number Pr with Pr = 1 throughout. The optimization approach used and the iterative solution
technique adopted are then described in Sec. II B. Section III is divided into three parts: A description
of the wide domain computations to find the snake solutions is given first in Sec. III A,followedbya
discussion of efforts to probe PCF at very low Re in Sec. III B, and then the calculations examining
the bursting phenomenon are presented in Sec. III C. A final discussion in Sec. IV recaps the various
results, provides some perspective, and then looks forward to future work.
II. FORMULATION
A. Stratified plane Couette flow
The usual plane Couette flow setup is considered in this paper of two (horizontal) parallel plates
separated by a distance 2h with the top plate moving at U
ˆ
x and the bottom plate moving at −U
ˆ
x.
Stable stratification is added by imposing that the fluid density is ρ
0
−
ρ
at the top plate and
ρ
0
+
ρ
at the bottom plate (gravity g is normal to the plates and directed downward from the top
plate to the bottom plate). Using the Boussinesq approximation (
ρ
ρ
0
), the governing equations
can be nondimensionalized using U , h, and
ρ
to give
∂u
∂t
+ u · ∇u =−∇p −Ri
b
ρ
ˆ
y +
1
Re
∇
2
u, (1)
∇ · u = 0, (2)
∂ρ
∂t
+ u · ∇ρ =
1
Re Pr
∇
2
ρ, (3)
where the bulk Richardson number Ri
b
, Reynolds number Re, and the Prandtl number Pr are
respectively defined as
Ri
b
:=
ρ
gh
ρ
0
U
2
, Re :=
Uh
ν
, Pr :=
ν
κ
. (4)
Here u = (u,v,w) is the velocity field, κ is the thermal diffusivity, the total dimensional density is
ρ
0
+ ρ
ρ
, p is the pressure, and ν is the kinematic viscosity. The boundary conditions are then
u(x, ± 1,z,t) =±1,ρ(x, ± 1,z,t) =∓1, (5)
083902-3
D. OLVERA AND R. R. KERSWELL
which admit the steady 1D solution
u = y
ˆ
x,ρ=−y. (6)
The (possibly large) disturbance fields away from this basic state,
˜
u(x,y,z,t) = u − y
ˆ
x, ˜ρ(x,y,z,t) = ρ + y, (7)
conveniently satisfy homogeneous boundary conditions at y =±1. Periodic boundary conditions
are used in both the streamwise (x) and spanwise (z) directions over wavelengths L
x
h and L
z
h,so
the (nondimensionalized) computational domain is L
x
× 2 × L
z
. The total energy of the disturbance
is taken as
E :=
1
2
˜
u
2
+
1
2
Ri
b
˜ρ
2
=
1
2
(u − y
ˆ
x)
2
+
1
2
Ri
b
(ρ + y)
2
, (8)
where (·) :=
1
V
(·)dV is a volume average.
B. Methods
To find the largest energy growth that a (finite-amplitude) perturbation (
˜
u, ˜ρ) can experience over
afixedtimeinterval[0,T ] requires seeking the global maximum of the constrained Lagrangian
L :=
1
2
˜
u(x,T )
2
+
1
2
Ri
b
˜ρ(x,T )
2
+ λ
1
2
˜
u(x,0)
2
+
1
2
Ri
b
˜ρ(x,0)
2
− E
0
+
T
0
π(x,t)∇ ·
˜
udt
+
T
0
ν (x,t) ·
∂
˜
u
∂t
+ y
∂
˜
u
∂x
+ ˜v
ˆ
x +
˜
u · ∇
˜
u + ∇
˜
p + Ri
b
˜ρ
ˆ
y −
1
Re
∇
2
˜
u
dt
+
T
0
τ (x,t) ·
∂ ˜ρ
∂t
+ y
∂ ˜ρ
∂x
− ˜v +
˜
u · ∇ ˜ρ −
1
Re Pr
∇
2
˜ρ
dt, (9)
where λ, π , ν, and τ are the Lagrange multiplier fields imposing the constraints that the initial
perturbation energy is E
0
, the perturbation is incompressible, and both the perturbation Navier-Stokes
equation and the density equation are satisfied respectively. Taking variations with respect to all
the degrees of freedom leads to the Euler-Lagrange equations that, beyond the aforementioned
constraints, comprise the dual evolution equations for the fields ν = ν
1
ˆ
x + ν
2
ˆ
y + ν
3
z and τ ,
∂ν
∂t
+ y
∂ν
∂x
− ν
1
ˆ
y − ν · (∇
˜
u)
T
+
˜
u · ∇ν + ∇π +
1
Re
∇
2
ν = τ ∇ ˜ρ − τ
ˆ
y, (10)
∂τ
∂t
+ y
∂τ
∂x
+
˜
u · ∇τ +
1
Re Pr
∇
2
τ = Ri
b
ν
2
, (11)
the temporal end conditions
˜
u(x,T ) + ν (x,T ) = 0, Ri
b
˜ρ(x,T ) + τ (x,T ) = 0, (12)
and the initial conditions
δL
δ
˜
u(x,0)
:= λ
˜
u(x,0) − ν (x,0) = 0,
δL
δ ˜ρ(x,0)
:= λ Ri
b
˜ρ(x,0) − τ (x,0) = 0. (13)
To eliminate spatial boundary terms, ν and τ are taken to obey the same homogeneous boundary
conditions as
˜
u and ˜ρ, and ν is further assumed incompressible to automatically satisfy the Euler-
Lagrange equation with respect to
˜
p. The solution strategy to find the global maximum of L is
iterative, starting with a guess for the initial perturbation (
˜
u, ˜ρ), which is then time stepped across
the time interval [0,T ] via the Navier-Stokes equation. The final values of
˜
u and ˜ρ initiate ν and τ
[via conditions (12)] for the time integration of the dual equations (10) and (11) backward to t = 0
where the fact that Eqs. (13) are generally not satisfied is used to update the form of the initial
083902-4
![FIG. 2. The top left shows the gain (bold red line), residual (black line), and (green line) for variational computations with initial energyE0 = 5 × 10−4 at Re = 180, Rib = 0, and T = 150 in a 4π × 2 × 16π domain [the computation is started with random noise, so the gain is O(10−9) after one step of the algorithm]. The top right shows the time evolution of the initial perturbation for iterations m = 15, 50, and 200. The state at t = 100 from the initial condition at m = 200 is used as an initial guess for the Newton GMRES method. A horizontal dotted line shows the level of kinetic energy of the subsequently converged solution shown in Fig. 3. The red dots indicate the evolution of the perturbation shown below at times t = 0, 20, 60, and 100 (top to bottom). The left column shows the xz cross section at y = 0 using eight contour levels (between −0.58 and 0.059) and the right column isocontours showing ±60% of the maximum streamwise perturbation velocity.](/figures/fig-2-the-top-left-shows-the-gain-bold-red-line-residual-2j8ugcka.webp)
