Optimising energy growth as a tool for finding exact coherent structures

TL;DR

This paper introduces an energy growth optimisation method to identify exact coherent structures in fluid flows, demonstrated through plane Couette flow and stratification effects, aiding understanding of transition to turbulence.

Contribution

It presents a novel approach using energy growth optimisation to find and analyze invariant solutions in fluid dynamics, including localised solutions and effects of stratification.

Findings

Successfully rediscovered spanwise-localised 'snake' solutions.

Probed phase space at low Reynolds numbers, approaching the laminar attractor.

Explored how stratification affects flow structures and their connectivity.

Abstract

We discuss how searching for finite amplitude disturbances of a given energy which maximise their subsequent energy growth after a certain later time $T$ can be used to probe 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 optimisation procedure will naturally select the one which lies nearest 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: a) rediscovering the spanwise-localised `snake' solutions of Schneider et al. (2010b); b) probing phase space at very low Reynolds numbers ($< 127.7$) where the constant linear-shear solution is believed to be the global attractor; and finally c) examining how the edge between laminar and turbulent flow evolves when stable stratification kills the turbulent attractor. We also show that the steady snake solution smoothly delocalises as unstable stratification is gradually turned on until it connects (via an intermediary global 3D solution) to 2D Rayleigh-Benard roll solutions.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p-{Ri}_b\rho \hat{\mathbf{y}}+\frac{1}{Re}{\nabla }^2\mathbf{u},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \rho }{\partial t}+\mathbf{u}\cdot \nabla \rho =\frac{1}{RePr}{\nabla }^2\rho$$

$$R{i}_{\text{b}}:=\frac{{\Delta }_{\rho }gh}{{\rho }_0{U}^2},Re:=\frac{Uh}{\nu },Pr:=\frac{\nu }{\kappa }.$$

$$u(x,\pm 1,z,t))=\pm 1\&\rho (x,\pm 1,z,t)=\mp 1.$$

$$\mathbf{u}=y\hat{\mathbf{x}}\&\rho =-y.$$

$$\tilde{\mathbf{u}}(x,y,z,t)=\mathbf{u}-y\hat{\mathbf{x}},\tilde{\rho }(x,y,z,t)=\rho +y,$$

$$E:=⟨\frac{1}{2}{\tilde{\mathbf{u}}}^2+\frac{1}{2}R{i}_b{\tilde{\rho }}^2⟩=⟨\frac{1}{2}(\mathbf{u}-y\hat{\mathbf{x}}{)}^2+\frac{1}{2}R{i}_b(\rho +y{)}^2⟩.$$

$$\mathcal{L}:=$$

$$\displaystyle+\int^{T}_{0}\langle\,\mbox{\boldmath$\nu$}(\mathbf{x},t)\cdot\biggl{[}\frac{\partial\mathbf{\tilde{u}}}{\partial t}+y\frac{\partial\mathbf{\tilde{u}}}{\partial x}+\tilde{v}\mathbf{\hat{x}}+\mathbf{\tilde{u}}\cdot\mbox{\boldmath$\nabla$}\mathbf{\tilde{u}}+\mbox{\boldmath$\nabla$}\tilde{p}+Ri_{b}\,\tilde{\rho}\mathbf{\hat{y}}-\frac{1}{Re}\nabla^{2}\mathbf{\tilde{u}}\biggr{]}\,\rangle\,dt$$

$$\displaystyle+\int^{T}_{0}\langle\,\tau(\mathbf{x},t)\cdot\biggl{[}\frac{\partial\tilde{\rho}}{\partial t}+y\frac{\partial\tilde{\rho}}{\partial x}-\tilde{v}+\mathbf{\tilde{u}}\cdot\mbox{\boldmath$\nabla$}\tilde{\rho}-\frac{1}{RePr}\nabla^{2}\tilde{\rho}\biggr{]}\,\rangle\,dt$$

$$\displaystyle\frac{\partial\mbox{\boldmath$\nu$}}{\partial t}+y\frac{\partial\mbox{\boldmath$\nu$}}{\partial x}-\nu_{1}\mathbf{\hat{y}}-\mbox{\boldmath$\nu$}\cdot(\mbox{\boldmath$\nabla$}\mathbf{\tilde{u}})^{T}+\mathbf{\tilde{u}}\cdot\mbox{\boldmath$\nabla$}\mbox{\boldmath$\nu$}+\mbox{\boldmath$\nabla$}\pi+\frac{1}{Re}\nabla^{2}\mbox{\boldmath$\nu$}$$

$$\displaystyle\frac{\partial\tau}{\partial t}+y\frac{\partial\tau}{\partial x}+\mathbf{\tilde{u}}\cdot\mbox{\boldmath$\nabla$}\tau+\frac{1}{RePr}\nabla^{2}\tau$$

$$\mathbf{\tilde{u}}(\mathbf{x},T)+\mbox{\boldmath$\nu$}(\mathbf{x},T)={\bf 0},\qquad Ri_{b}\,\tilde{\rho}(\mathbf{x},T)+\tau(\mathbf{x},T)=0,$$

$$\frac{\delta{\cal L}}{\delta\mathbf{\tilde{u}}(\mathbf{x},0)}:=\lambda\mathbf{\tilde{u}}(\mathbf{x},0)-\mbox{\boldmath$\nu$}(\mathbf{x},0)={\bf 0},\qquad\frac{\delta{\cal L}}{\delta\tilde{\rho}(\mathbf{x},0)}:=\lambda Ri_{b}\,\tilde{\rho}(\mathbf{x},0)-\tau(\mathbf{x},0)=0.$$

$$(\,\mathbf{\tilde{u}}^{(m+1)}(\mathbf{x},0),\tilde{\rho}^{(m+1)}(\mathbf{x},0)\,)\,=\,(\,\mathbf{\tilde{u}}^{(m)}(\mathbf{x},0),\tilde{\rho}^{(m)}(\mathbf{x},0)\,)+\frac{\epsilon}{\lambda}\biggl{(}\frac{\delta{\cal L}}{\delta\mathbf{\tilde{u}}(\mathbf{x},0)^{(m)}},\frac{1}{Ri_{b}}\frac{\delta{\cal L}}{\delta\tilde{\rho}(\mathbf{x},0)^{(m)}}\biggr{)}$$

$${\cal R}(m)\,:=\,\frac{\,\langle\delta{\cal L}/\delta\mathbf{\tilde{u}}(\mathbf{x},0)^{(m)}\,\rangle\,+\,\frac{1}{Ri_{b}}\langle\,\delta{\cal L}/\delta\tilde{\rho}(\mathbf{x},0)^{(m)}\,\rangle}{\langle\,\mbox{\boldmath$\nu$}(\mathbf{x},0)^{(m)}\rangle\,+\,\frac{1}{Ri_{b}}\langle\tau(\mathbf{x},0)^{(m)}\rangle}$$

$$\mathcal{G}(E_0,T):=\frac{E(T)}{E(0)}$$

$$(\tilde{\mathbf{u}},\tilde{\rho },\tilde{p})=(\tilde{\mathbf{u}}(y),\tilde{\rho }(y),\tilde{p}(y))e^{i\beta z+\sigma t}$$

$$\sigma \tilde{u}$$

$$\sigma \tilde{v}$$

$$\sigma \tilde{w}$$

$$0$$

$$\sigma \tilde{\rho }$$

$$\frac{\partial \mathbf{U}}{\partial \tau }+\mathbf{U}\cdot \nabla \mathbf{U}$$

$$\frac{\partial \Theta }{\partial \tau }+\mathbf{U}\cdot \nabla \Theta$$

$$\nabla \cdot \mathbf{U}$$

$$u(x,\pm 1,z,t))=\pm 1\&\theta (x,\pm 1,z,t)=\mp 1.$$

$$t=(RePr)\tau ,\mathbf{u}=\mathbf{U}/(RePr),\rho =\Theta ,p=P/(RePr{)}^2$$

Full text

I Introduction

Optimisation has proved a powerful tool to extract information from the Navier-Stokes equations. In the shear flow transition problem, optimising over all possible infinitesimal disturbances to find the one which maximises 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’ Reddy & Henningson (1993); Henningson & Reddy (1994), ‘nonmodal instability’ Schmid (2007) or ‘optimal perturbation theory’ Butler & Farrell (1992) (see the reviews Grossmann (2000); Schmid (2007) and book Schmid & Henningson (2001)), the approach reveals key aspects of the linearised dynamics around the reference state which has helped to interpret finite-time flow phenomena and pick apart what causes 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 (e.g. Trefethen et al. (1993); Trefethen & Embree (2005); Schmid (2007); Schmid & Brandt (2014)). 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 Waleffe (1995); Dauchot & Manneville (1997)

Conceptually, the remedy to this is simple: let competing disturbances seeking to maximise 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 Pringle & Kerswell (2010); Cherubini et al. (2010); Kerswell et al. (2014). This, however, doubles the number of parameters ($E_{0}$ joins $T$) over which the results must be interpreted and leads to a fully nonlinear, non-convex optimisation problem where much less is known about its possibly multiple solutions (local and 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 in the transition problem Pringle & Kerswell (2010); Cherubini et al. (2010, 2011); Monokrousos et al. (2011); Pringle et al. (2012); Rabin et al. (2012); Cherubini et al. (2012, 2013); Duguet et al. (2013); Pringle et al. (2015) (see the review Kerswell et al. (2014)). For example, one can ask what is the smallest (most ‘dangerous’) energy disturbance which can trigger transition by some time $T$, with the answer in the large $T$ limit labelled the minimal seed for transition Pringle et al. (2012); Rabin et al. (2012); Pringle et al. (2015). The minimal seeds which emerge from this procedure are fully-localised and are therefore realistic targets for experimental investigations (e.g. Pringle et al. (2015)).

The optimisation 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\rightarrow\infty$ (all other disturbances have to decay eventually). The new state to which these disturbances are drawn needn’t be a turbulent attractor and so the minimal disturbance to reach another simple stable state can also be calculated (see §6.2 in Rabin (2013)). What is not so clear is whether the approach can find unstable solutions although a similar line of reasoning seems to hold. The optimisation procedure would be expected to select disturbances from the energy hypersphere which 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 new solution. A sufficiently close visit should yield flow states which can then be converged to the new 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-localised ‘snake’ solutions of Schneider et al. Schneider et al. (2010b) building upon the prior exploratory work of Rabin (see §6.3 in Rabin (2013)) who identified the key role played by the choice of $E_{0}$. In a wide geometry, snake solutions coexist with repeated copies of Nagata’s well known solution Nagata (1990) in a narrow geometry and, not surprisingly, the stable manifolds of the (lower energy) snake solutions pass closer (in energy norm) to the simple shear solution in phase space than those of Nagata’s solutions. However, the latter offer the possibility of greater energy growth as they lead to a global flow state and hence are preferred by the optimisation algorithm if they pass close to the energy 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 (§6.3 Rabin (2013)). We complete this calculation here by recomputing these optimal disturbances and converging out both snake solutions.

Armed with this success, we then probe phase space of pCf at very low Reynolds numbers ($<127.7$ Waleffe (2003)) looking for new solutions where the (basic) constant shear solution is believed to be the global attractor (a proof only exists for $Re<20.7$ Joseph (1966)). We find evidence of solutions but these turn out only to be the ghosts of known solutions at higher Reynolds numbers. Finally, as another example of how the optimisation approach can be utilised, we examine how the edge between laminar and turbulent flow evolves when stable stratification suppresses the turbulence. With the turbulent attractor present, a ‘bursting’ phenomenon 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 the solution acting as the edge state ceases to exist. This then indicates that the bursting is directly related to the presence of the unstable manifold of the edge state directed away from the uniform shear solution in phase space.

The paper starts with the formulation of the stratified plane Couette flow problem in §II.1 which introduces the three non-dimensional 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 optimisation approach used and the iterative solution technique adopted are then described in §II.2. The results section §III is divided into 3 parts: a description of the wide domain computations to find the snake solutions is given first in §III.1; followed by a discussion of efforts to probe pCf at very low $Re$ in §III.2; and then the calculations examining the bursting phenomenon are presented in §III.3. A final discussion in §IV recaps the various results, provides some prospectives and then looks forward to future work.

II Formulation

II.1 Stratified plane Couette flow

The usual plane Couette flow set-up is considered in this paper of two (horizontal) parallel plates separated by a distance $2h$ with the top plate moving at $U\mathbf{\hat{x}}$ and the bottom plate moving at $-U\mathbf{\hat{x}}$. Stable stratification is added by imposing that the fluid density is $\rho_{0}-\Delta_{\rho}$ at the top plate and $\rho_{0}+\Delta_{\rho}$ at the bottom plate (gravity $g$ is normal to the plates and directed downwards from the top plate to the bottom plate). Using the Boussinesq approximation ($\Delta_{\rho}\ll\rho_{0}$), the governing equations can be non-dimensionalised using $U$, $h$ and $\Delta_{\rho}$ to give

[TABLE]

[TABLE]

[TABLE]

where the bulk Richardson number $Ri_{\textrm{b}}$, Reynolds number $Re$, and the Prandtl number $Pr$ (always set to 1 in this study) are respectively defined as:

[TABLE]

Here $\mathbf{u}=(u,v,w)$ is the velocity field, $\kappa$ the thermal diffusivity, the total dimensional density is $\rho_{0}+\rho\Delta_{\rho}$, $p$ is the pressure and $\nu$ is the kinematic viscosity. The boundary conditions are then

[TABLE]

which admit the steady 1D solution

[TABLE]

The (possibly large) disturbance fields away from this basic state,

[TABLE]

conveniently satisfy homogeneous boundary conditions at $y=\pm 1$. Periodic boundary conditions are used in both the ($x$) streamwise and ($z$) spanwise directions over wavelengths $L_{x}h$ and $L_{z}h$ so that the (non-dimensionalised) computational domain is $L_{x}\times 2\times L_{z}$. The total ‘energy’ of the disturbance is taken as

[TABLE]

where $\langle\,(\cdot)\,\rangle:=\frac{1}{V}\iiint\,(\cdot)\,dV$ is a volume average.

II.2 Methods

To find the largest energy growth that a (finite-amplitude) perturbation $(\mathbf{\tilde{u}},\tilde{\rho})$ can experience over a fixed time interval $[0,T]$ requires seeking the global maximum of the constrained Lagrangian

[TABLE]

where $\lambda$, $\pi$, $\nu$ and $\tau$ 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 which, beyond the aforementioned constraints, comprise of the ‘dual’ evolution equations for the fields, \mbox{\boldmath\nu}=\nu_{1}\mathbf{\hat{x}}+\nu_{2}\mathbf{\hat{y}}+\nu_{3}\mathbf{\hat{z}} and $\tau$,

[TABLE]

the temporal end conditions

[TABLE]

and the initial conditions

[TABLE]

To eliminate spatial boundary terms, $\nu$ and $\tau$ are taken to obey the same homogeneous boundary conditions as $\mathbf{\tilde{u}}$ and $\tilde{\rho}$, and $\nu$ is further assumed incompressible to automatically satisfy the Euler-Lagrange equation with respect to $\tilde{p}$. The solution strategy to find the global maximum of ${\cal L}$ is iterative, starting with a ‘guess’ for the initial perturbation $(\mathbf{\tilde{u}},\tilde{\rho})$ which is then time-stepped across the time interval $[0,T]$ via the Navier-Stokes equation. The final values of $\mathbf{\tilde{u}}$ and $\tilde{\rho}$ ‘initiate’ $\nu$ and $\tau$ (via conditions (14) ) for the time integration of the dual equations (12) and (13) backwards to $t=0$ where the fact that the equations (15) are generally not satisfied is used to update the form of the initial perturbation (subject to it staying of total energy $E_{0}$) in the direction of increasing ${\cal L}$. We use a simple steepest ascent method where

[TABLE]

where, for example, $\mathbf{\tilde{u}}(\mathbf{x},t)^{(m)}$ is the $m$th iterate, and $\lambda$ is subsequently chosen to ensure the new $(m+1)^{th}$ iterate has energy $E_{0}$ as discussed in Pringle & Kerswell (2010); Pringle et al. (2012); Rabin et al. (2012); Kerswell et al. (2014), but other approaches are possible (e.g. Monokrousos et al. (2011); Duguet et al. (2013); Cherubini et al. (2010, 2011, 2013)). This direct and adjoint looping method is now well used (e.g. see the reviews Luchini & Bottaro (2014); Kerswell et al. (2014)) but of course there is no guarantee that the global maximum always emerges for this fully nonlinear problem. The hoped output of the procedure is the optimal initial condition - the optimal disturbance - which experiences the largest growth over a time horizon $[0,T]$ of all initial conditions with the same initial total energy $E_{0}$ and so is a function of $T$ and $E_{0}$.

The status of the iterative procedure is monitored by computing the residual

[TABLE]

which should approach 0 for convergence. The time integrations of the Boussinesq equations forward in time and the dual equations backwards in time were carried out using an adapted version of the parallelized DNS code ‘Diablo’ (Taylor (2008) and http://www.damtp.cam.ac.uk/user/jrt51/files.html) which uses a third-order mixed Runge-Kutta-Wray/Crank-Nicolson timestepper. The horizontal directions are periodic and treated pseudospectrally, while a second-order finite-difference discretization is used in the cross-stream direction. The resolution used was typically 64 Fourier modes per $2\pi$ in $x$ and $z$ and 128 finite difference points in $y$. If needed, this resolution was doubled to ensure numerical accuracy. Diablo was coupled it to a Newton-Raphson-GMRES algorithm (Viswanath, 2007) which allowed ECS to be converged from a good guess and continued around in parameter space (most notably by varying $Ri_{b}$). In what follows, it is usually more convenient to plot the gain ${\cal G}$ of iterates which is defined as

[TABLE]

where $E_{0}=E(0)$ is the a priori-fixed initial perturbation energy: maximizing this, of course, is equivalent to maximizing the final total energy over all perturbations with given initial energy $E_{0}$.

III Results

III.1 Wide domain pCf: Rediscovering Snakes

Optimal energy growth calculations were performed at $Re=180$ with stratification turned off ($Ri_{b}=0$) in the two wide domains, $4\pi\times 2\times 8\pi$ and $4\pi\times 2\times 16\pi$ where the snake solutions are known to exist Schneider et al. (2010a, b). The calculations were initiated with random initial conditions (energy scattered in the lowest modes) normalised so that the total kinetic energy was $E_{0}$. If $E_{0}$ is too small, only the immediate neigbourhood of the constant shear solution is explored with a nonlinear version of the 2D linear optimal - a global set of streamwise rolls - emerging as the optimal. If $E_{0}$ is too large, the optimal perturbation leads to another global state resembling multiple copies of Nagata’s solution (e.g. see figure 6.13 in Rabin (2013)). Flows states from this optimal evolution could, presumably, be used to converge Nagata’s solution but this was not pursued (it would be more numerically efficient to treat a narrower domain which supports just one spanwise wavelength if that was an objective). At intermediary $E_{0}$, the optimal perturbation is more spanwise localised and stays spanwise localised as it evolves into a state suggestive of a snake solution (see figure 1).

Figure 1 (top left) shows the convergence features of the optimal growth calculation at this intermediary initial energy using $T=150$ in the $16\pi$ wide box. The iterative algorithm is clearly struggling to converge - the residual remains $O(1)$ - yet the gain has levelled off and most importantly, a plateau has emerged in the time evolution of the optimal iterates (see top right figure). This signals a close approach to a constant-energy saddle (either an equilibrium or travelling wave) and it is from here that we take a flow snapshot (specifically at $t=95$) which is spanwise localised: see figure 2(upper). This state converged in 26 Newton steps to a very similar looking steady solution - see figure 2(lower) - which is the equilibrium snake solution (hereafter referred to as EQ) of Schneider et al. (2010a).

The optimisation procedure proceeded much more slowly in the narrower domain $4\pi\times 2\times 8\pi$ for reasons which are unclear but again a plateau is eventually established in the optimal evolution: see figure 3(left). Two flow states were extracted from this - see figure 4(upper) - with one converging and one apparently not - see figure 4(lower) and the convergence behaviour in figure 3(right). The converged state this time was the travelling wave snake solution (hereafter referred to as TW) of Schneider et al. (2010a).

Once the snake solutions, EQ and TW, had been found, they could be traced around in $(Re,Ri_{b})$ parameter space ($Pr=1$ to keep things manageable) using the Newton-GMRES algorithm. Fixing $Ri_{b}=0$ and varying $Re$ in the $16\pi$ wide box reproduced the ‘snakes and ladders’ plot of Schneider et al. (2010b) (their figure 2) confirming the identity of the solutions: see figure 5(left). Interestingly, fixing $Re=180$ and varying $Ri_{b}$ also shows snaking in EQ. Further continuing this solution to negative $Ri_{b}$ (unstable stratification) - see figure 6 - reveals that EQ connects to a Nagata’s solution (of different spanwise wavenumber $\beta=1.5$ to that which the TW connects in figure 2 of Schneider et al. (2010b) where $\beta=1$ ) just before this bifurcates off a 2D convective roll solution familiar from the Rayleigh-Benard problem (the extra shear from the boundaries does not affect this solution except to determine its orientation). Salewski et al. Salewski et al. (2017) have also recently found this same bifurcation sequence in rotating plane Couette flow which is known to be closely related to the Rayleigh-Benard problem.

In connecting to the $\beta=1.5$ Nagata solution, the EQ snake solution has to delocalise and figure 7 shows this is a gradual process as $Ri_{b}$ decreases from 0 as opposed to that found for $Ri_{b}$ increasing from 0 when snaking occurs. The key observation for $Ri_{b}<0$ is that the (spatial) spanwise decay of the snake disappears once the threshold for convective instability at $Ri_{b}=-3.2943\times 10^{-3}$ is crossed (the critical Rayleigh number$\,:=-Ri_{b}Re^{2}Pr=1708/16$; Salewski et al. Salewski et al. (2017) see the same phenomenon in rotating plane Couette flow - see their figure 4). This can be understood by examining the linear operator about the linearly sheared state for the least-(spatially)-damped, temporally steady eigenfunction since the deviation away from this linearly-sheared state becomes vanishingly small in the spanwise tails of the snake Gibson & Brand (2014). The snake becomes streamwise-independent in its tail regions suggesting analysis of the linear eigenvalue problem for 2D disturbances independent of the streamwise direction i.e.

[TABLE]

so that

[TABLE]

Normally, $\beta$ is assumed real and the eigenvalue problem is scrutinised for complex $\sigma$ with $\Re e(\sigma)=0$ to find spatially-periodic, neutral eigenfunctions. Here, instead, the interest is in (real) $\sigma=0$ and complex $\beta$ to find spatially-decaying steady eigenfunctions (since EQ is steady). Of primary interest is the eigenfunction with the smallest amplitude of $\Im m(\beta)$ where $\sigma=0$ which indicates the probable rate of spatial evanescence in the spanwise direction of a steady ECS when the amplitude gets small (see §4.1 of Gibson & Brand (2014)). Figure 8 shows this neutral curve in the complex $\beta$ plane on the left and the usual neutral curve as viewed in the wavenumber$-$control parameter plane is shown on the right. For $Ri_{b}\leq-3.2943\times 10^{-3}$, neutral spatially-periodic eigenfunctions can exist but otherwise $\Im m(\beta)\neq 0$. Figure 9 shows that in the tail regions there is indeed good correspondence between the expected spatial decay and the numerically observed decay close to the linear instability threshold at $Ri_{b}=-3.2943\times 10^{-3}$.

So, in summary, by using optimal energy growth, we have managed to rediscover the snake solutions of Schneider et al. (2010a, b). Having shown that this approach works, we now turn our attention to a region of parameter space in pCf where no solutions are currently known beyond the linearly-sheared base state.

III.2 Very low $Re$ in pCf

In this subsection, we turn stratification off ($Ri_{b}=0$) and set $Re=100$ which is above the energy stability threshold of $20.7$ Joseph (1966) up to which the basic sheared state is provably unique and below $127.7$ which is the current best estimate of when other solutions start to exist Nagata (1990); Waleffe (2003). A geometry of $4\pi\times 2\times 2\pi$ and a shortish target time of $T=20$ were chosen and $E_{0}$ gradually increased until the (nonlinear versions of the) linear optimal perturbations ($LOP_{1}$ and $LOP_{2}$) shown in figure 10 were no longer found.

At $E_{0}=10^{-3}$, a new nonlinear optimal perturbation ($NLOP_{1}$) emerges which actually experiences larger growth than both the $LOP_{1}$ or $LOP_{2}$ at earlier times. This initial condition is 3D, localised towards one wall and evolves into 2 pairs of wavy fast-slow streaks: see figure 11. However, there is no discernable energy plateau in its evolution so $T$ was increased to 40 whereupon a different nonlinear optimal ($NLOP_{2}$) emerges at $E_{0}\approx 1.5\times 10^{-3}$: see figure 12 for its structure at $E_{0}=8\times 10^{-3}$. This optimal gives rise to an energy plateau in its subsequent evolution as the initial energy is increased: see figure 13 for the situation at $E_{0}=8\times 10^{-3}$. In this, a good candidate to initiate a convergence attempt is the flow state at $t=16$, however, this simply converges to the linearly sheared base state. Repeating the calculation at $Re=130$ (with $E_{0}=6\times 10^{-3}$), again taking the flow state at $t=16$, does converge but to Nagata’s solution. A number of other searches were done for $Re<128$ with all guesses converging to the base state and for $Re>128$ where all attempts converged smoothly to Nagata’s solution (the geometry is slightly sub-optimal in that the saddle node value for Nagata’s solution is $\approx 128$ rather than $127.7$). No evidence emerged of any other state beyond Nagata’s solution existing during these computations adding further weight to the view that the linearly-sheared base state is a global attractor up to $Re=127.7$.

Despite this apparent simplicity, the results of a systematic optimal energy growth analysis over the $(E_{0},T)$ plane are still quite rich. Figure 14 indicates the various global optimals found at $Re=100$ and $T=10$, $20$ and $40$ over the interval $E_{0}\in[0,0.01]$. At $T=20$, for example, 4 different optimals emerge with all being the global optimal at some $E_{0}$: see Figure 14(right). That the nonlinear energy growth problem is nontrivial even in the absence of any exact coherent structures, is presumably because phase space is already starting to structure itself to incorporate such states at slightly higher $Re$.

III.3 Stratified pCf for high $Ri_{b}$

As our third (and final) application of the optimal energy growth technique, we consider what happens to the edge, which in phase space separates states that become turbulent from those that relaminarise, when the addition of stable stratification kills the turbulent attractor. Direct numerical simulations Olvera & Kerswell (2017) indicate that large energy growth can still occur and so we seek an explanation why by examining how the optimal energy growth perturbation changes as $Ri_{b}$ increases from 0 in a geometry of $2\pi\times 2\times\pi$ where $Re=400$.

There are two important values of $Ri_{b}$: the value $Ri_{b}^{t}(Re)$ beyond which no turbulent attractor exists and the limiting value $Ri_{b}^{m}(Re)$ for the existence of the global equilibrium EQ7 Gibson et al. (2009) which is the edge state in this geometry (  $10^{-6}<Ri_{b}^{t}(400)\lesssim 10^{-2}$ and $Ri_{b}^{m}(400)\approx 0.057$). Previous work at $Re=1000$ by Rabin et al. Rabin et al. (2012) and Eaves & Caulfield Eaves & Caulfield (2015) has identified the minimal seed for transition at $Ri_{b}=0$ and $0<Ri_{b}<Ri_{b}^{t}$ respectively. Our focus here is $Ri_{b}^{t}\,<Ri_{b}$ but we start by looking at one value of $Ri_{b}=10^{-6}<Ri_{b}^{t}$ to make contact with this earlier work.

Figure 15 shows partially-converged optimal perturbations at $E_{0}=1.5\times 10^{-4}$ using two target times $T=50$ and $100$ which both indicate that turbulence is triggered. Taking $T=50$ and reducing $E_{0}$ until turbulence is no longer triggered by any iterate or the final converged optimal suggests that $E_{c}$, the critical energy threshold for transition is between $E_{0}=3\times 10^{-5}$ and $5\times 10^{-5}$: see figure 16. The form of the optimal at $3\times 10^{-5}$ is still the nonlinearly adjusted $LOP$ whereas the optimal at $5\times 10^{-5}$ resembles the minimal seed found by Rabin et al. (2012) at $Ri_{b}=0$ (see figure 5 of Rabin et al. (2012)). The optimal for $E_{0}=5\times 10^{-5}$ leads to an evolution from which a good enough starting guess can be extracted to converge EQ7 (not shown). This is another case where the approach has worked to identify an unstable state albeit only one which is minimally unstable because, as an edge state, it only has one unstable direction or eigenvalue.

Moving to stronger stratification, figure 17 shows converged optimals for $E_{0}>E_{c}$ with $Ri_{b}=0.01$ (which is just above $Ri_{b}^{t}$) and $Ri_{b}=0.04$. The former shows a transient turbulent episode while the latter only a ‘burst’ of energy growth which then decays away. Again if $E_{0}$ is near enough to $E_{c}$, the evolution of the optimal transiently visits the neighbourhood of the edge state sufficiently closely to be able to select a flow state which subsequently converges. For $Ri_{b}=0.04$, we find $2.102\times 10^{-4}\,<\,E_{c}\,<\,2.108\times 10^{-4}$ and a flow state taken at $t=95$ from the $E_{0}=2.108\times 10^{-4}$ optimal evolution converges in just 7 steps to EQ7 (not shown).

Increasing the stratification further, figure 18 compares the results of working at $Ri_{b}=0.055<Ri_{b}^{m}$ where a burst is still discernable (notice the upward curvature of the energy curve) and $Ri_{b}=0.06>Ri_{b}^{m}$ where it is not. Clearly, the ‘bursting’ is produced by the unstable manifold of the edge state (directed away in phase space from the linearly sheared state) and largely vanishes when the edge state disappears (i.e. the upward curvature disappears). In fact, phase space for $Ri_{b}\gtrsim Ri_{b}^{m}$ should still reflect the memory of the manifold but this ebbs away with increasing $Ri_{b}$. To confirm this simple explanation, we carried out a couple of checks. The first was to confirm that the state reached at the energy peak at pre- and post-$Ri_{b}^{m}$ stratifications is roughly the same - see figure 19. And the second was attempting to converge a flow state on an optimal energy plateau for the $Ri_{b}=0.06$ (taken from the optimal evolution at $E_{0}=4.2\times 10^{-4}$: see the red dot in figure 18(right)). This failed to converge at $Re=400$ but did converge to EQ7 at $Re=500$ as $Ri_{b}^{m}(400)\,<\,Ri_{b}=0.06\,<\,Ri_{b}^{m}(500)$: see the inset of figure 18(right).

The conclusion of this subsection is then that when the turbulent attractor disappears, one still can see ‘bursting’ which is caused by the flow trajectory being repelled out from the vicinity of the linearly-sheared base state by the unstable manifold of the edge state. This dominant (first) feature of the transition process therefore survives well after the suppression of turbulence by stable stratification and relies only on the edge state continuing to exist. To further confirm this picture, the minimal seed to trigger transition or latterly this bursting (as $Ri_{b}$ increases) remains essentially the same across all $Ri_{b}$ studied including $Ri_{b}=10^{-6}$ which is essentially the unstratified case of Rabin et al. (2012): see figure 20. However, what does change is the evolution of the minimal seed once it has experienced the burst in energy growth. Finally we collect together in figure 21 the data collected on $E_{c}$ - the minimal energy to reach ‘above’ the edge - as a function of $Ri_{b}$ (the curve stops at $Ri_{b}=0.06$ since it is difficult to identify an edge beyond this point).

IV Discussion

In this paper, we have demonstrated that an optimisation technique, in which the energy growth of a finite-amplitude disturbance to a known solution is maximised, can be used to generate flow fields subsequently convergeable via a Newton-GMRES algorithm to another ‘nearby’ solution of the Navier-Stokes equations. That this may be possible has been noticed before in the particular case of an edge state to which a flow initiated by the minimal seed will get infinitesimally close during its evolution Pringle et al. (2012); Rabin et al. (2012); Kerswell et al. (2014). Now this has been confirmed in section III.3. What was not clear before, however, was whether the technique could be used to find unstable solutions more generally located in phase space in the absence of an edge. The rediscovery of both steady and travelling wave snake solutions in section III.1 demonstrates that the technique can also work in this situation too.

The technique has then been used to probe very low $Re$ pCf with negative results adding further weight to the view that the linear-shear solution is indeed the global attractor in pCf up to $Re=127.7$. Finally, a ‘bursting’ phenomenon was investigated in stably-stratified pCf and found to be produced by the unstable manifold of the edge state: when the edge state vanished on increasing the stratification so did the bursting. The key here is that the optimisation technique allowed the bursting to be found if it existed as the stratification was varied.

Finally, the steady spanwise-localised snake solution has been been found to connect, via a global 3D state, to the 2D rolls of the Rayleigh-Benard problem. In doing so, a particularly simple delocalisation process has been found where the spanwise tails of the snake gradually reduce their spatial decay rate until this vanishes at the global linear instability threshold whereupon the state is then global. This would seem a very generic phenomenon where a localised state moves from a region of subcriticality to one of supercriticality or vice versa.

In terms of further work, the most obvious question is whether the optimisation technique can be used to find periodic orbits in which the energy varies in time. So far, only constant energy solutions (equilibria and travelling waves) have been sought and the identification of an energy plateau during the optimal’s evolution has been central to indicate a ‘close approach’. In principle, a mildly fluctuating plateau should also be recognisable and provided the period is not too long, convergence should still be feasible. Another issue is whether the approach can provide any insight into phase space around a linearly unstable solution. Here the answer surely depends on the time scale $\tau$ of the linear instability. Typically the growth rates of linear instabilities in shear flows are much smaller than the typical instantaneous growth rates of energy growth optimals (e.g. in plane Poiseuille flow, the transient energy growth is at least an order of magnitude more than that from the linear instability over O(100) advective times). This then suggests that the optimisation technique will simply ignore the linear instability for $T\ll\tau$ in preference to more potent mechanisms which give better, albeit transient, energy growth. One such could be a nearby stable manifold of another solution as explored here.

As a final comment, it’s worth emphasizing that the principles which make the optimisation approach work here in a fluid mechanical context hold true also for any dynamical system where solution multiplicity is suspected. What makes shear flows so interesting, of course, is this is now understood to be the generic situation where not only are there multiple solutions but these are typically unstable yet instrumental in determining the fluid dynamics of the flow (e.g. the bursting of section III.3). This optimisation approach then looks to be a valuable new addition to a theoretician’s toolbox.

Acknowledgements.

DO would like to thanks CONACYT for the award of a scholarship which has supported his PhD studies and the availability of free HPC time on the University of Bristol’s BlueCrystal supercomputer. Both DO and RRK would like to thank Tom Eaves for sharing his nonlinear optimal perturbation routine (here converted to one focussed on total energy rather than total dissipation), John Taylor for sharing his time-stepping code ‘Diablo’ and Dan Lucas for help parallelizing the GMRES algorithm used. We also are grateful for encouragement from the rest of the EPSRC-funded ‘MUST’ team at Cambridge received during the course of this work.

Appendix: Rayleigh Benard Convection

In the normal Rayleigh-Benard set-up (e.g. Drazin & Reid (1981) equations (8.6)-(8.8) with nonlinearities reinstated ) the fully nonlinear equations for the total flow field $\mathbf{U}$ and temperature field $\Theta$ are

[TABLE]

subject to boundary conditions

[TABLE]

Stratified plane Couette flow (as described by equations (1)-(3)  ) is retrieved under the following transformation

[TABLE]

so that $Ra=-Ri_{b}\,Re^{2}\,Pr$. The critical $Ra$ for linear instability is $Ra_{crit}:=1708/16$ which is unchanged by introducing a unidirectional shear (e.g. Kelly (1977)) (the factor of $1/16$ is because the half-channel width has been used to non-dimensionalize the system).