Hostname: page-component-7b7567c8c-2n9c9 Total loading time: 0 Render date: 2026-11-15T15:04:57.759Z Has data issue: false hasContentIssue false
  • English
  • Français

Surfing the edge: using feedback control to find nonlinear solutions

Published online by Cambridge University Press:  13 October 2017

A. P. WillisOpen the ORCID record for A. P. Willis[Opens in a new window],
Y. Duguet,
O. Omel’chenko and
M. Wolfrum
A. P. Willis*
Affiliation:
School of Mathematics and Statistics, University of Sheffield S3 7RH, UK
Y. Duguet
Affiliation:
LIMSI-CNRS, UPR 3251, Université Paris-Saclay, F-91403, Orsay, France
O. Omel’chenko
Affiliation:
Weierstrass Institute, Mohrenstrasse 39, 10117 Berlin, Germany
M. Wolfrum
Affiliation:
Weierstrass Institute, Mohrenstrasse 39, 10117 Berlin, Germany
*
†Email address for correspondence: a.p.willis@sheffield.ac.uk
Rights & Permissions [Opens in a new window]

Abstract

Many transitional wall-bounded shear flows are characterised by the coexistence in state space of laminar and turbulent regimes. Probing the edge boundary between the two attractors has led in the last decade to the numerical discovery of new (unstable) solutions to the incompressible Navier–Stokes equations. However, the iterative bisection method used to compute edge states can become prohibitively costly for large systems. Here we suggest a simple feedback control strategy to stabilise edge states, hence accelerating their numerical identification by several orders of magnitude. The method is illustrated for several configurations of cylindrical pipe flow. Travelling waves solutions are identified as edge states, and can be isolated rapidly in only one short numerical run. A new branch of solutions is also identified. When the edge state is a periodic orbit or chaotic state, the feedback control does not converge precisely to solutions of the uncontrolled system, but nevertheless brings the dynamics very close to the original edge manifold in a single run. We discuss the opportunities offered by the speed and simplicity of this new method to probe the structure of both state space and parameter space.

JFM classification

Nonlinear Dynamical Systems: Nonlinear Dynamical SystemsFlow Control: Instability controlInstability: Transition to turbulence

Information

Type
Papers
Information
Journal of Fluid Mechanics , Volume 831 , 25 November 2017, pp. 579 - 591
Check for updates
Copyright
© 2017 Cambridge University Press 

1 Introduction

In recent years, there has been increasing evidence for the dynamical importance of exact unstable solutions, such as travelling waves and periodic orbits, for weakly turbulent flows (Kawahara, Uhlmann & Van Veen Reference Kawahara, Uhlmann and Van Veen2012; Chandler & Kerswell Reference Chandler and Kerswell2013; Cvitanović Reference Cvitanović2013; Willis, Short & Cvitanović Reference Willis, Short and Cvitanović2016). This is best understood so far for the case of subcritical wall-bounded shear flows, where the base flow is linearly stable at flow rates where turbulence is also sustained. These unstable solutions are the only non-trivial solutions, and are all disconnected from the trivial laminar solution (Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007).

Most solutions discovered so far appear in saddle-node bifurcations and hence arise in pairs. The solution on the lower-energy branch often belongs to the laminar–turbulent separatrix, ‘the edge’, which separates initial conditions that lead to relaminarisation from those that experience transition to turbulence (Itano & Toh Reference Itano and Toh2001; Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006). The attractor within the edge separatrix is labelled the ‘edge state’. When a solution located on the edge has only one real unstable eigenvalue, it acts as an attractor within the edge. There may be more than one edge state for a given system. Upper-branch solutions are thought to be either embedded in the turbulent attractor, or, together with the lower-branch solution, bracket the turbulent dynamics (Gibson, Halcrow & Cvitanović Reference Gibson, Halcrow and Cvitanović2008).

Knowledge of an invariant solution within the edge has become the first step of a common strategy to unfold the whole bifurcation diagram of the system (Kreilos & Eckhardt Reference Kreilos and Eckhardt2012; Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013; Ritter, Mellibovsky & Avila Reference Ritter, Mellibovsky and Avila2016), but identification of invariant solutions, otherwise known as exact coherent states (ECS), remains a demanding task. In the absence of a well-understood sequence of bifurcations from the base flow, i.e. when all classical continuation methods fail (Tuckerman & Barkley Reference Tuckerman and Barkley2000), several strategies have been employed to find such solutions. The size of the systems under study, usually with $10^{4}$ to $10^{9}$ degrees of freedom, requires specialised iterative strategies, which in particular avoid any operation involving the storage of full matrices. A perhaps more challenging issue is that a good initial guess for the identification is usually unavailable. During the 1990s and early 2000s, the main strategy used to bypass this issue was homotopy (see e.g. Kerswell Reference Kerswell2005). This relies on an efficient root finder, typically a Newton–Raphson solver coupled to an arc-length continuation algorithm. This method may track the solutions in parameter space, but does not help to identify unconnected branches. To apply this method, a companion problem featuring a linear instability of the base flow must first be identified. Once a bifurcated solution is found, it is continued nonlinearly to a solution of the original problem. Success requires very good intuition for suitable homotopies that might work, plus involves the complexity of solving for the multiple systems. It is also unclear how continued solutions participate in the dynamics.

Later, bisection methods, which make use of timesteppers, began to gain popularity: the amplitude of an arbitrary perturbation to the base flow is rescaled repeatedly, until a trajectory is found which for a sufficiently long time stays away from both the laminar and the turbulent state (Itano & Toh Reference Itano and Toh2001; Skufca et al. Reference Skufca, Yorke and Eckhardt2006; Schneider, Eckhardt & Yorke Reference Schneider, Eckhardt and Yorke2007). Such a trajectory is usually chaotic, but with some luck, if there exists a regular solution with only one unstable eigendirection, there is hope that the bisection algorithm will converge to it (Schneider et al. Reference Schneider, Gibson, Lagha, De Lillo and Eckhardt2008). Imposing discrete symmetries on the dynamics often reduces the number of unstable eigendirections of the ECS contained in the associated subspace (Duguet, Willis & Kerswell Reference Duguet, Willis and Kerswell2008b ). As a result the likelihood of identifying symmetry-invariant solutions is increased, but their relevance to the non-symmetric dynamics is uncertain. In the more general case where the edge trajectory is chaotic, recurrence analysis can sometimes identify an approach to a simple state, which can be converged using a good Newton solver, possibly enhanced by some globalisation method (Viswanath Reference Viswanath2007; Duguet, Pringle & Kerswell Reference Duguet, Pringle and Kerswell2008a ). The main drawbacks of bisection methods are the hazardous chances of success, the difficulty to converge interesting recurrent parts of the dynamics using Newton methods, and more importantly their cost. Indeed, the bisection method requires $n\sim -\text{log}_{2}(\unicode[STIX]{x1D700}_{M})$ individual runs until machine precision $\unicode[STIX]{x1D700}_{M}$ is reached, at which point the process needs to be restarted $M$ times until a simple edge state has been reached. The total cost is hence of $Mn\geqslant$ 100 runs at least, with the runs getting longer as the accuracy improves. For instance, in Khapko et al. (Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016), each bisection required a total of $Mn\approx$ 400 runs, which corresponds to $O(10^{6}{-}10^{7})$ CPU hours. This high cost makes parametric studies infeasible in practice. Other numerical methods have been suggested as alternatives to the bisection–rootfinder combination, e.g. iterative adjoint optimisation methods (Farazmand Reference Farazmand2016; Olvera & Kerswell Reference Olvera and Kerswell2017) although they involve significant mathematical and computational complexity. In summary, it is always desirable to find simpler and less expensive alternative methods.

In the present paper we demonstrate how unstable edge states can be found numerically, by introducing into the original system a ‘control’ term that counteracts the edge instability and is able to stabilise unstable states without altering them significantly. This enables the system to dynamically approach the edge state in a single simulation of the controlled system. A key property of a control for this purpose is that it should influence only the stability, not the structure of the a priori unknown target, i.e. the control should act as non-invasively as possible. At the same time, it should efficiently force a large set of initial conditions towards the states of interest. Using control as a tool to numerically find and analyse unstable objects in a complex dynamical system has already been successfully applied to a variety of problems. A classic example is the time-delayed feedback control (Pyragas Reference Pyragas1992) able to stabilise certain periodic orbits provided their period is known. Another example is selective frequency damping, which, by filtering all non-zero temporal frequencies, may stabilise steady state solutions (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006). In the latter case, however, the steady solutions must not possess unstable non-oscillatory eigenvalues (Vyazmina Reference Vyazmina2010), which is frequently the case for edge solutions in subcritical shear flows. Note that none of these methods is designed to target the edge manifold.

In this article we propose a remarkably simple linear feedback control able to constrain the dynamics to the edge manifold and to stabilise invariant solutions that are stable within the edge. This scheme has recently been applied to track unstable chimera states in systems of non-locally coupled phase oscillators (Sieber, Omel’chenko & Wolfrum Reference Sieber, Omel’chenko and Wolfrum2014; Wolfrum, Omel’chenko & Sieber Reference Wolfrum, Omel’chenko and Sieber2015). In an interesting analogy to shear flows, in systems of coupled oscillators ‘chimeras’ are metastable chaotic states that coexist with a fully synchronised (‘laminar’) state (Panaggio & Abrams Reference Panaggio and Abrams2015). The chimera also appear in a saddle-node bifurcation, and the corresponding unstable lower branch also acts as an edge state separating the synchronous from the chimera regime. In the current fluid context, we first need to restrict the dynamics to the low-energy levels characteristic of the edge manifold. So far, the bisection method has revealed that the dynamics on the edge manifold is relatively low-dimensional, especially compared to turbulent flow at equivalent parameters values (Duguet et al. Reference Duguet, Willis and Kerswell2008b ). We wish to take advantage of this property to extract invariant solutions from the controlled flow, in the same way that solutions have been extracted from edge bisections. Our aim here, however, is to achieve this in a single simulation, rather than via the more expensive bisection approach.

In the following section we describe the feedback control, then in § 3 we present the numerical set-up. In § 4 we apply the control scheme to several different pipe flow cases, first in a restricted domain and then in an extended domain allowing for axial localisation. Finally, possible applications for the methods are discussed in the concluding § 5.

2 Stabilisation of lower-branch equilibria by feedback control

In its most simple form, the control scheme makes a system parameter $\unicode[STIX]{x1D707}$ state dependent by imposing a linear relation between the parameter and an observable $A(t)$

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D707}(t)=\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D705}(A_{0}-A(t)),\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is the control gain. (For the shear flow problems we will identify $\unicode[STIX]{x1D707}$ with $Re$ and $A$ with a component of the perturbation energy.) How such a proportional control can be used to stabilise the unstable part of a folded branch can be explained by applying it to the normal form for a saddle-node bifurcation,

(2.2) $$\begin{eqnarray}{\dot{x}}(t)=\unicode[STIX]{x1D707}-x^{2}(t),\end{eqnarray}$$

where equilibria lie on the parabola $\unicode[STIX]{x1D707}=x^{2}$ . Inserting the constraint (2.1) in to (2.2) and choosing the simplest possible observable $A\equiv x$ , we arrive at the equation

(2.3) $$\begin{eqnarray}{\dot{x}}=f(x,\unicode[STIX]{x1D705},\unicode[STIX]{x1D707}_{0},x_{0}):=\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D705}(x_{0}-x)-x^{2}.\end{eqnarray}$$

Geometrically, the resulting dynamics can be understood by considering its representation in the $(\unicode[STIX]{x1D707},x)$ plane. For $\unicode[STIX]{x1D705}=0$ we have the uncontrolled system (2.2), where for a fixed choice of the parameter the dynamics is constrained to a vertical line $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{0}$ in the $(\unicode[STIX]{x1D707},x)$ plane. Equilibria are located at the intersections of the vertical line with the parabola, $x=\pm \sqrt{\unicode[STIX]{x1D707}}$ . According to (2.2), $x$ increases with time in the region enclosed by the parabola, for $x^{2}<\unicode[STIX]{x1D707}$ , while it decreases outside of this region. Hence, we conclude that for the uncontrolled system the equilibria on the upper branch, $x=+\sqrt{\unicode[STIX]{x1D707}}$ , of the parabola are dynamically stable, while those on the lower branch, $x=-\sqrt{\unicode[STIX]{x1D707}}$ , are unstable.

Restricting now the dynamics to a slanted straight line, as imposed by (2.1) with $\unicode[STIX]{x1D705}>0$ , we obtain the controlled system (2.3), which has equilibria

(2.4) $$\begin{eqnarray}x_{\pm }={\textstyle \frac{1}{2}}\left(-\unicode[STIX]{x1D705}\pm \sqrt{\unicode[STIX]{x1D705}^{2}+4(\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D705}x_{0})}\right),\end{eqnarray}$$

provided that

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D705}^{2}+4(\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D705}x_{0})=(\unicode[STIX]{x1D705}+2x_{0})^{2}+4(\unicode[STIX]{x1D707}_{0}-x_{0}^{2})>0.\end{eqnarray}$$

Note that the equilibria (2.4) coincide with the equilibria of the uncontrolled system (2.2) for the particular choice $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D705}(x_{0}-x_{\pm })$ . The controlled system has derivative

$$\begin{eqnarray}\unicode[STIX]{x2202}_{x}f(x_{\pm },\unicode[STIX]{x1D705},\unicode[STIX]{x1D707}_{0},x_{0})=-\unicode[STIX]{x1D705}-2x_{\pm }=\mp \sqrt{\unicode[STIX]{x1D705}^{2}+4(\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D705}x_{0})}.\end{eqnarray}$$

Thus the upper intersection $x_{+}$ of the controlled system is stable, while $x_{-}$ is unstable. In particular, whilst the lower branch $x=-\sqrt{\unicode[STIX]{x1D707}}$ is unstable for $\unicode[STIX]{x1D705}=0$ , it is possible to choose a slanted straight line such that both intersections $x_{\pm }$ meet the parabola on the lower branch. The upper intersection $x_{+}$ is now a stabilised equilibrium on the lower branch. Further, by varying the control gain $\unicode[STIX]{x1D705}$ , we can now sweep the straight line about the pivot point $(\unicode[STIX]{x1D707}_{0},x_{0})$ to track the branch of stable equilibria of the controlled system. This is possible until the straight line becomes tangential to the parabola, when $x_{\pm }$ disappear via a saddle-node bifurcation for $|\unicode[STIX]{x1D705}+2x_{0}|>2\sqrt{x_{0}^{2}-\unicode[STIX]{x1D707}_{0}}.$

The proportional control (2.1) can be applied in the same way to a system with more complicated dynamics, shown schematically in figure 1, provided that it contains a saddle-node bifurcation. Such a folded branch of invariant manifolds generically develops a single transversal direction of instability, similar to the example above. Using an observable that varies along the fold, we can apply the control. Whenever we obtain a controlled state with constant observable, it will exactly coincide with a dynamical state of the uncontrolled system with fixed parameter related to the constant observable by (2.1). Only the stability properties are affected by the control, i.e. the control is non-invasive for such states. Controlled trajectories with periodic or chaotic fluctuations in $A(t)$ may still show good quantitative and qualitative agreement with an uncontrolled trajectory for a parameter value close to the time average $\langle \unicode[STIX]{x1D707}\rangle _{t}$ . The control may be described as non-invasive on average (Sieber et al. Reference Sieber, Omel’chenko and Wolfrum2014). In this way, in the shear flow problems it will be possible to obtain periodic orbits to be used as initial states for a corresponding Newton solver. Also, by sweeping the control gain dynamically from a known stable dynamical regime of the uncontrolled system, the controlled system may yield access to qualitatively different dynamical regimes along the lower branch, i.e. the edge.

Figure 1.

Sketch in an $(A,Re)$ representation of the stability of the uncontrolled (a) and controlled system using the feedback control method (b) based on the constraint (2.1) corresponding to the dashed slanted line. The labels ‘UB’ and ‘LB’ refer to the stable upper and unstable lower branch in the uncontrolled system, respectively. In (b) the dynamics is restricted to a dashed slanted line. Rotating the line about a pivot point, fixed points on LB are stabilised and can be tracked dynamically.

3 Numerical methodology

3.1 Direct numerical simulation of pipe flow

We consider here the case of cylindrical pipe flow with a fixed mass flux. The control parameter $\unicode[STIX]{x1D707}$ in equation (2.1) is here the Reynolds number $Re=UD/\unicode[STIX]{x1D708}$ , where $U$ is the constant bulk velocity, $D$ the pipe diameter and $\unicode[STIX]{x1D708}$ the kinematic viscosity of the fluid. Using $U$ and $D$ as the scales, the non-dimensional laminar state is characterised by an axial velocity profile $\boldsymbol{u}_{lam}(r)=2(1-(2r)^{2})\boldsymbol{e}_{z}$ for $0\leqslant r\leqslant 1/2$ , driven by a homogeneous pressure gradient $(\unicode[STIX]{x2202}_{z}p_{lam})\boldsymbol{e}_{z}$ . The full velocity field $\boldsymbol{u}=(u,v,w)$ , which contains the radial, azimuthal and axial components of the velocity, satisfies, together with the pressure $p$ , the incompressible Navier–Stokes equations at all times:

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u} & = & \displaystyle 0,\end{eqnarray}$$
(3.2) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)\boldsymbol{u} & = & \displaystyle -\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}.\end{eqnarray}$$

The flow satisfies the no-slip velocity conditions at the wall. The code ensures that this condition is satisfied to machine precision via the influence matrix method (Kleiser & Schumann Reference Kleiser and Schumann1980). We advance (3.2) in time using the hybrid spectral finite-difference code openpipeflow.org (Willis Reference Willis2017). The code employs Fourier expansions $\text{e}^{\text{i}(2\unicode[STIX]{x1D6FC}kz+mm_{0}\unicode[STIX]{x1D703})}$ , where $k$ and $m$ are respectively the axial and the azimuthal wavenumbers. This imposes streamwise periodicity with a wavelength $L=\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}$ in units $D$ . The integer $m_{0}$ indicates the degree of rotational symmetry of the flow field in the azimuthal direction, $m_{0}=2$ refers to a twofold symmetry while $m_{0}=1$ corresponds to the absence of any discrete rotational symmetry. In some cases we also impose mirror symmetry $Z\boldsymbol{u}=\boldsymbol{u}$ about the plane $\unicode[STIX]{x1D703}=0$ , where

(3.3) $$\begin{eqnarray}\displaystyle Z:(u,v,w,p)(r,\unicode[STIX]{x1D703},z)\rightarrow (u,-v,w,p)(r,-\unicode[STIX]{x1D703},z), & & \displaystyle\end{eqnarray}$$

and the shift-and-reflect symmetry $S\boldsymbol{u}=\boldsymbol{u}$ , where

(3.4) $$\begin{eqnarray}\displaystyle S:(u,v,w,p)(r,\unicode[STIX]{x1D703},z)\rightarrow (u,-v,w,p)(r,-\unicode[STIX]{x1D703},z+L/2). & & \displaystyle\end{eqnarray}$$

In §§ 4.1 and 4.2 we consider simulations with $m_{0}=2$ in domains with respectively $\unicode[STIX]{x1D6FC}=1.25$ and $0.12566$ . With dealiasing, variables are evaluated on grids in $r\times \unicode[STIX]{x1D703}\times z$ of respectively $64\times 48\times 72$ and $64\times 48\times 576$ points. In § 4.3, for $m_{0}=1$ and $\unicode[STIX]{x1D6FC}=1.25$ , the resolution is $64\times 96\times 72$ .

3.2 Implementation of the feedback control

For the Navier–Stokes system (3.1)–(3.2), the feedback control (2.1) is applied by controlling the Reynolds number

(3.5) $$\begin{eqnarray}Re(t)=Re_{0}+\unicode[STIX]{x1D705}(A_{0}-A(t))\end{eqnarray}$$

using the scalar observable

(3.6) $$\begin{eqnarray}A(t)=\int _{r_{i}}^{r_{w}}(u^{2}+v^{2})\,r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\,\text{d}z,\end{eqnarray}$$

where $r_{w}=0.5$ corresponds to the wall and $r_{i}$ is set to 0.35. This observable $A$ , considered only over the near-wall annular region $r_{i}\leqslant r\leqslant r_{w}$ , is a robust signature of the presence of active coherent structures contributing to the turbulence production. Several other observables were examined, such as energies that include all components of the velocity, integrated over the whole domain, or the energy input related to the average wall drag. No significant difference between those observables was noted in the application of the method, which indicates the robustness of the control approach.

The factor $1/Re$ appears in (3.2) as the coefficient of the viscous term, which is treated implicitly in our timestepping scheme, and the factor is used in the evaluation of time stepping matrices. Now that $Re$ depends on time, rather than recalculating these matrices every time step, we consider the governing equation (3.2) in the form

(3.7) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}-\frac{1}{Re_{r}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u} & ~=\ & \displaystyle \left(\frac{1}{Re}-\frac{1}{Re_{r}}\right)\,\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}-\left(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)\boldsymbol{u}-\unicode[STIX]{x1D735}p,\end{eqnarray}$$

where $Re_{r}$ is a reference value kept within 1 % of $Re(t)$ . The diffusion term on the left-hand side of (3.7) is treated implicitly. The small correction term with coefficient $(1/Re-1/Re_{r})$ , along with the rest of the terms on the right-hand side, is treated explicitly. No numerical instability has been observed.

We now explain how to choose the constants $A_{0}$ and $Re_{0}$ . $A_{0}$ was chosen to be zero so that $Re(t)\rightarrow Re_{0}$ as $A(t)\rightarrow 0$ , i.e. as the observable approaches zero, $Re(t)$ approaches a maximum value. Since at low $Re$ turbulence frequently relaminarises, the value of $Re_{0}$ must be chosen to be sufficiently high to avoid such relaminarisations. It is determined indirectly from the value of $\unicode[STIX]{x1D705}$ necessary to stabilise turbulence at lower $Re$ and lower amplitudes $A$ : For a trial $Re_{0}$ and knowing $(Re,A)$ for a given (uncontrolled) velocity field, from (2.1) we can estimate a starting value for $\unicode[STIX]{x1D705}$ . Starting controlled simulations with $\unicode[STIX]{x1D705}$ at two or three times this value, in order to enforce a lower $A(t)$ , we then adjusted $Re_{0}$ to avoid relaminarisations. Thereafter, the value $Re_{0}=10^{4}$ was found to be adequate for all further simulations, giving the pivot point $(Re_{0},A_{0})=(10^{4},0)$ .

Using a large value of $Re_{0}$ at first sight suggests that spatial resolution issues may arise for large $Re(t)$ . In practice, no fully turbulent simulation is run at values of $Re$ as high as $Re_{0}$ . Large $Re(t)$ occurs only for small energies $A(t)$ , for which the spectral drop-off is more rapid. It has been verified that the drop-off in energy spectra is sufficient that all stabilised states considered in this paper are well resolved.

Finally, we introduce time dependence to the control gain $\unicode[STIX]{x1D705}(t)$ . Continuation along the lower-energy branch towards increasing $Re$ requires the slope of the line given by (2.2) in the $(Re,A)$ plane to decrease, i.e. $\unicode[STIX]{x1D705}(t)$ must increase with time. Once the energy $A(t)$ has dropped below its level for the uncontrolled system, it becomes more difficult to change $\unicode[STIX]{x1D705}$ in steps without each step causing an initial wild fluctuation in $Re(t)$ , which in turn is likely to cause immediate relaminarisation. This issue is avoided by allowing $\unicode[STIX]{x1D705}(t)$ to vary smoothly. A simple exponential form for $\unicode[STIX]{x1D705}(t)$ proved effective

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D705}(t)=\unicode[STIX]{x1D705}_{0}\,\text{e}^{t/T},\end{eqnarray}$$

where $T$ is a chosen constant. $T$ is the expected time scale for the rate of decrease in the resulting controlled $A(t)$ , and the approach is quite insensitive to the precise choice of $T$ (example in § 4.1). As far as convergence to invariant solutions is concerned, $T$ should be longer than the decay time of the least-damped mode for the stabilised solution of the controlled system. Since this is unknown, $T$ is simply chosen to be long relative to the $O(1{-}10)$ time scale of the variability observed in energy and dissipation time series.

4 Results

4.1 Short periodic pipe with $m_{0}=2$

We present first the results with the feedback control for a short pipe of length $L=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}$ with $\unicode[STIX]{x1D6FC}=1.25$ , $m_{0}=2$ and no additional symmetry. This system was considered in Kerswell & Tutty (Reference Kerswell and Tutty2007), and was shown to possess a lower-branch travelling wave (TW) with only one unstable real eigenvalue. It was later verified in Duguet et al. (Reference Duguet, Willis and Kerswell2008b ) that bisection in this system could identify a travelling wave as an edge state, but further, that bisection could converge to two travelling waves from different families, depending on the initial condition. We describe now how the feedback controller described in the previous section successfully identified and stabilised another travelling wave solution on the edge.

Figure 2.

Application of the feedback control for $\unicode[STIX]{x1D6FC}=1.25$ and $m_{0}=2$ . (a) $(Re,A)$ representation, control with $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}\exp (t/T)$ , $T=250$ (blue). Cases with $T$ doubled/halved converge to the same stabilised lower branch (grey/cyan). In the inset, the solution set is compared with results from arc-length continuation (using the Newton method, pink) and continuation in the opposite direction using the control with $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}\exp (-t/T)$ (green). (b) $A(t)$ uncontrolled (red), controlled with time-varying $\unicode[STIX]{x1D705}(t)$ (blue) and controlled with constant $\unicode[STIX]{x1D705}$ for the stabilisation of the TW (green) (zoom in inset).

Figure 3.

(a) Cross-section of the stabilised TW solution with $\unicode[STIX]{x1D6FC}=1.25$ and $m_{0}=2$ . Streamwise velocity perturbation in colour (from dark to white) and cross-stream components (vectors). (b) Eigenvalue spectrum for the TW solution with (blue circles) and without control (red crosses).

Figure 4.

Application of the feedback control with $L=25D$ , $m_{0}=2$ and mirror symmetry imposed. (a) $(Re,A)$ representation. The arrow shows the direction of rotation of the slanted straight line, displayed here are two different times. (b) Time series of $Re(t)$ with time-varying (blue/dark) and constant (green/light) $\unicode[STIX]{x1D705}(t)$ for the stabilisation of the RPO. (c) $(E_{3d},\unicode[STIX]{x1D6FD})$ representation of the RPO converged using the feedback control (red), the unstable RPO of the uncontrolled system (blue) at $Re=1450$ , close to each other relative to the laminar state (black circle). (d) Spectrum of Floquet exponents for the RPO solution with (blue circles) and without control (red crosses).

We start from a turbulent snapshot at $Re=2000$ , then vary $\unicode[STIX]{x1D705}(t)$ according to (3.8) with $T=250$ and $\unicode[STIX]{x1D705}_{0}=4\times 10^{5}$ . As shown in figure 2(a,b), the time series of $A(t)$ initially displays erratic fluctuations until $t\approx 300~D/U$ , where the dynamics appears to be smoother and fluctuations are damped. From $(A,Re)\approx (5\times 10^{-3},1800)$ onwards, the backward fold in figure 2(a) indicates that the dynamics stays away from the laminar state ( $A=0$ ) and that the whole lower branch is stabilised as far as the simulation is run, up to $Re\approx 3000$ (blue). Halving $T$ (grey), more of the lower branch is quickly reached. Convergence to the same branch is also observed for $T$ doubled (cyan). For the first case, $T=250$ , the state $t=375\,D/U$ is used as an initial condition for a controlled simulation with $\dot{\unicode[STIX]{x1D705}}$ set to zero, i.e.  $\unicode[STIX]{x1D705}$ is held fixed at $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{1}=4\times 10^{5}\times \text{e}^{1.5}$ (green line in figure 2 b). The new trajectory converges towards a fixed point, i.e. a travelling wave solution, with constant $A$ . The associated velocity field in a cross-section is displayed in figure 3(a). While the symmetry was not imposed, the time integration with the control converged to a solution invariant under $S$ , being an S2 state of (Kerswell & Tutty Reference Kerswell and Tutty2007; Pringle, Duguet & Kerswell Reference Pringle, Duguet and Kerswell2009). It was verified using a Newton solver (Willis, Cvitanović & Avila Reference Willis, Cvitanović and Avila2013) that the corresponding solution is also solution to the original uncontrolled Navier–Stokes equations. The linear stability of the states are compared in their eigenspectrum $\{\unicode[STIX]{x1D70E}\}$ , calculated using an Arnoldi algorithm, with and without the control ( $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{1}$ and $\unicode[STIX]{x1D705}=0$ , respectively), figure 3(b), where $\text{Re}(\unicode[STIX]{x1D70E})>0$ indicates instability. Whereas the TW is unstable in the uncontrolled system, it is stabilised by the feedback control. Applying the control with $\dot{\unicode[STIX]{x1D705}}<0$ (here $\unicode[STIX]{x1D705}(t)=\unicode[STIX]{x1D705}_{0}\,\text{e}^{-t/T}$ ) the simulation also remains on the stabilised lower branch, as demonstrated by the green line in figure 2(a). Sweeping along the lower branch in either direction is therefore a very efficient alternative to numerical continuation with a Newton scheme for rapid numerical continuation of the TW solution in this region.

4.2 Long periodic pipe with $m_{0}=2$

We consider next the case $L=25D$ with $m_{0}=2$ and mirror symmetry $Z$ . This system, considered by Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013) and later by Chantry et al. (Reference Chantry, Willis and Kerswell2014), is known to possess a relative periodic orbit as an edge state, in the form of a weakly modulated travelling wave whose velocity field is axially localised. As in § 3.1, a simulation starting from a turbulent state at $Re=1900$ is first used with time-varying $\unicode[STIX]{x1D705}(t)$ following (3.8) with $\unicode[STIX]{x1D705}_{0}=1.62\times 10^{5}$ and $T=400$ . Figure 4(a) shows the dynamics in an $(Re,A)$ representation, where for $A\leqslant 0.03$ a lower branch has been captured by the scheme. At time $t=300\,D/U$ , $\unicode[STIX]{x1D705}$ is held fixed at $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{1}=1.62\times 10^{5}\times \text{e}^{0.75}$ . Close-ups on the time series of $A(t)$ and $Re(t)$ show weak periodic modulations around steady values, indicating convergence towards a relative periodic orbit (RPO) (figure 4 b). The converged RPO has a time-averaged Reynolds number $\langle Re\rangle _{t}\approx 1450$ . A point on this orbit was passed to the Newton solver for the uncontrolled system at  $Re=1450$ , and converged in one Newton step. The controlled and uncontrolled RPOs do not exactly coincide, but are close. Figure 4(c) shows a projection in ( $E_{3d},\unicode[STIX]{x1D6FD}$ ) (observe the scale and distance from the origin), where $E_{3d}=\int |\boldsymbol{u}-\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D703},z}|^{2}\,r\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\,\text{d}z$ and $1+\unicode[STIX]{x1D6FD}=\langle \unicode[STIX]{x2202}_{z}p\rangle /\unicode[STIX]{x2202}_{z}p_{lam}$ . Their periods are respectively $10.67$ and $10.61$ units $D/U$ . Stabilisation is substantiated in the altered Floquet exponents, displayed in figure 4(d). Unlike for the previous short pipe, where the travelling wave solution was an edge state, the controller has achieved convergence to an RPO that is not precisely a solution to the uncontrolled equations. The control stays approximately non-invasive on average, however, and the isolated solution is close (only a single Newton step was required to go from one solution to the other with a relative error of less than $10^{-6}$ ). It manages this approach in only one short computation.

4.3 Short periodic pipe with $m_{0}=1$

Bisections in periodic pipes without any rotational symmetry imposed ( $m_{0}=1$ ) have been performed (Schneider et al. Reference Schneider, Eckhardt and Yorke2007; Duguet et al. Reference Duguet, Willis and Kerswell2008b ) but to date no edge state with a simple dynamics has been reported for this case – the dynamics on the edge remains chaotic. Here we apply mirror symmetry to an $m_{0}=1$ computation in a short pipe with $\unicode[STIX]{x1D6FC}=1.25$ . The control by itself did not stabilise any TW solution in this case, but the dynamics on the edge is substantially less chaotic than the turbulent state. A narrow window of unsteady yet quiescent dynamics of the controlled system was identified in figure 5(a) near the intersection with the indicated straight line. Taking the state at this point as an initial guess in a Newton search converged to a new TW solution at $Re=2680$ , whose cross-section and eigenvalue spectrum are displayed in figures 5(b) and 5(c), respectively. Note that this new TW does not satisfy the shift-and-reflect symmetry. It is characterised by one real unstable eigenvalue and one pair of complex unstable eigenvalues. Using the control the real eigenvalue is stabilised, and the resulting solution becomes only weakly unstable. Also shown is a cross-section of the solution after numerical continuation to $Re=1428$ , where the solution appears in a saddle-node bifurcation at this $\unicode[STIX]{x1D6FC}=1.25$ .

Figure 5.

Application of the feedback control with $\unicode[STIX]{x1D6FC}=1.25$ , $m_{0}=1$ . (a) $(Re,A)$ representation. The arrow shows the direction of rotation of the slanted straight line, which is displayed at the time where an approach to an unstable TW is identified. (b) Cross-section of the stabilised TW solution (same colour coding as in figure 3), $Re=2680$ . (c) Eigenvalue spectrum for the TW solution with (blue circles) and without control (red crosses). (d) Cross-section after arc-length continuation to $Re=1428$ , $\unicode[STIX]{x1D6FC}=1.25$ .

5 Conclusions

We have demonstrated that the present feedback control method is able to stabilise invariant solutions that are edge states of the uncontrolled pipe flow system. Stabilised travelling waves correspond also to solutions of the uncontrolled Navier–Stokes equations. The control strategy is thus non-invasive in the case where the edge state of the original system is a travelling wave (Sieber et al. Reference Sieber, Omel’chenko and Wolfrum2014). Relative periodic orbits or chaotic regimes stabilised using the method, however, are not precisely invariant solutions of the uncontrolled Navier–Stokes equations, but are expected to be sufficiently close for rapid continuation towards the original system. In both cases, the feedback control method proves very efficient at bringing the system close to the original edge state in only one short computation, whereas the same task using the bisection method would require $O(100)$ times as many simulations.

A similar feedback control has recently been considered for turbulence in stratified plane Couette flow (Taylor et al. Reference Taylor, Deusebio, Caulfield and Kerswell2016) to attain a target energy. The target energy needs to match that of the invariant solutions, which is not known a priori in such large domains.

Coupled with a Newton–Krylov solver, this new method, thanks to its low cost, can be used to probe the bifurcation diagram of the original system and explore the edge manifold without prior bisection. Adapting this method to other parallel shear flows, or to any other spatio-temporal system with hysteresis, is straightforward. Since one run is sufficient to uncover the whole lower branch parametrised by $Re$ (or potentially another parameter such as e.g.  $\unicode[STIX]{x1D6FC}$ or an additional force), it can either be used to perform wider parametric studies, or to identify the interesting regimes within the edge when multiple bifurcations occur (Khapko et al. Reference Khapko, Duguet, Kreilos, Schlatter, Eckhardt and Henningson2014). We recommend implementation of the present control scheme as a simple precursor to the more expensive bisection and complex Newton–Krylov implementation.

Acknowledgements

Discussions and sections of this work were completed at the Kavli Institute for Theoretical Physics, supported in part under grant no. NSF PHY11-25915. A.W. acknowledges support by EPSRC EP/P000959/1. M.W. and O.O. acknowledge the support by DFG in the framework of SFB 910, ‘Control of self-organizing nonlinear systems: theoretical methods and concepts of application’.

References

Åkervik, E., Brandt, L., Henningson, D. S., Hœpffner, J., Marxen, O. & Schlatter, P. 2006 Steady solutions of the Navier–Stokes equations by selective frequency damping. Phys. Fluids 18 (6), 068102.CrossRefGoogle Scholar
Avila, M., Mellibovsky, F., Roland, N. & Hof, B. 2013 Streamwise-localized solutions at the onset of turbulence in pipe flow. Phys. Rev. Lett. 110 (22), 224502.CrossRefGoogle ScholarPubMed
Chandler, G. J. & Kerswell, R. R. 2013 Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow. J. Fluid Mech. 722, 554–595.CrossRefGoogle Scholar
Chantry, M., Willis, A. P. & Kerswell, R. R. 2014 Genesis of streamwise-localized solutions from globally periodic traveling waves in pipe flow. Phys. Rev. Lett. 112 (16), 164501.CrossRefGoogle ScholarPubMed
Cvitanović, P. 2013 Recurrent flows: the clockwork behind turbulence. J. Fluid Mech. 726, 1–4.CrossRefGoogle Scholar
Duguet, Y., Pringle, C. C. & Kerswell, R. R. 2008a Relative periodic orbits in transitional pipe flow. Phys. Fluids 20 (11), 114102.CrossRefGoogle Scholar
Duguet, Y., Willis, A. P. & Kerswell, R. R. 2008b Transition in pipe flow: the saddle structure on the boundary of turbulence. J. Fluid Mech. 613, 255–274.CrossRefGoogle Scholar
Eckhardt, B., Schneider, T. M., Hof, B. & Westerweel, J. 2007 Turbulence transition in pipe flow. Annu. Rev. Fluid Mech. 39, 447–468.CrossRefGoogle Scholar
Farazmand, M. 2016 An adjoint-based approach for finding invariant solutions of Navier–Stokes equations. J. Fluid Mech. 795, 278–312.CrossRefGoogle Scholar
Gibson, J. F., Halcrow, J. & Cvitanović, P. 2008 Visualizing the geometry of state space in plane Couette flow. J. Fluid Mech. 611, 107–130.CrossRefGoogle Scholar
Itano, T. & Toh, S. 2001 The dynamics of bursting process in wall turbulence. J. Phys. Soc. Japan 70, 703–716.CrossRefGoogle Scholar
Kawahara, G., Uhlmann, M. & Van Veen, L. 2012 The significance of simple invariant solutions in turbulent flows. Annu. Rev. Fluid Mech. 44, 203–225.CrossRefGoogle Scholar
Kerswell, R. R. 2005 Recent progress in understanding the transition to turbulence in a pipe. Nonlinearity 18 (6), R17.CrossRefGoogle Scholar
Kerswell, R. R. & Tutty, O. R. 2007 Recurrence of travelling waves in transitional pipe flow. J. Fluid Mech. 584, 69–102.CrossRefGoogle Scholar
Khapko, T., Duguet, Y., Kreilos, T., Schlatter, P., Eckhardt, B. & Henningson, D. S. 2014 Complexity of localised coherent structures in a boundary-layer flow. Eur. Phys. J. E 37, 32.Google Scholar
Khapko, T., Kreilos, T., Schlatter, P., Duguet, Y., Eckhardt, B. & Henningson, D. S. 2016 Edge states as mediators of bypass transition in boundary-layer flows. J. Fluid Mech. 801, R2.CrossRefGoogle Scholar
Kleiser, L. & Schumann, U. 1980 Treatment of incompressibility and boundary conditions in 3-d numerical spectral simulations of plane channel flows. In Proceedings of Third GAMMA Conference on Numerical Methods in Fluid Mechanics, pp. 165–173. Springer.CrossRefGoogle Scholar
Kreilos, T. & Eckhardt, B. 2012 Periodic orbits near onset of chaos in plane Couette flow. J. Nonlinear Sci. 22 (4), 047505.Google ScholarPubMed
Olvera, D. & Kerswell, R. R. 2017 Optimising energy growth as a tool for finding exact coherent structures. Phys. Rev. Fluids 2, 083902.CrossRefGoogle Scholar
Panaggio, M. J. & Abrams, D. M. 2015 Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators. Nonlinearity 28 (3), R67.CrossRefGoogle Scholar
Pringle, C. C., Duguet, Y. & Kerswell, R. R. 2009 Highly symmetric travelling waves in pipe flow. Phil. Trans. R. Soc. Lond. A 367 (1888), 457–472.Google ScholarPubMed
Pyragas, K. 1992 Continuous control of chaos by self-controlling feedback. Phys. Lett. A 170 (6), 421–428.CrossRefGoogle Scholar
Ritter, P., Mellibovsky, F. & Avila, M. 2016 Emergence of spatio-temporal dynamics from exact coherent solutions in pipe flow. New J. Phys. 18 (8), 083031.CrossRefGoogle Scholar
Schneider, T. M., Eckhardt, B. & Yorke, J. A. 2007 Turbulence transition and the edge of chaos in pipe flow. Phys. Rev. Lett. 99, 034502.CrossRefGoogle ScholarPubMed
Schneider, T. M., Gibson, J. F., Lagha, M., De Lillo, F. & Eckhardt, B. 2008 Laminar-turbulent boundary in plane Couette flow. Phys. Rev. E 78 (3), 037301.Google ScholarPubMed
Sieber, J., Omel’chenko, E. & Wolfrum, M. 2014 Controlling unstable chaos: stabilizing chimera states by feedback. Phys. Rev. Lett. 112 (5), 054102.CrossRefGoogle ScholarPubMed
Skufca, J. D., Yorke, J. A. & Eckhardt, B. 2006 Edge of chaos in a parallel shear flow. Phys. Rev. Lett. 96, 174101.CrossRefGoogle Scholar
Taylor, J. R., Deusebio, E., Caulfield, C. P. & Kerswell, R. R. 2016 A new method for isolating turbulent states in transitional stratified plane Couette flow. J. Fluid Mech. 808, R1.CrossRefGoogle Scholar
Tuckerman, L. S. & Barkley, D. 2000 Bifurcation analysis for timesteppers. In Numerical Methods for Bifurcation Problems and Large-scale Dynamical Systems, pp. 453–466. Springer.CrossRefGoogle Scholar
Viswanath, D. 2007 Recurrent motions within plane Couette turbulence. J. Fluid Mech. 580, 339–358.CrossRefGoogle Scholar
Vyazmina, E.2010 Bifurcations in a swirling flow. PhD thesis, École Polytechnique, France.Google Scholar
Willis, A. P. 2017 The Openpipeflow Navier–Stokes solver. SoftwareX 6, 124–127.CrossRefGoogle Scholar
Willis, A. P., Cvitanović, P. & Avila, M. 2013 Revealing the state space of turbulent pipe flow by symmetry reduction. J. Fluid Mech. 721, 514–540.CrossRefGoogle Scholar
Willis, A. P., Short, K. Y. & Cvitanović, P. 2016 Symmetry reduction in high dimensions, illustrated in a turbulent pipe. Phys. Rev. E 93 (2), 022204.CrossRefGoogle Scholar
Wolfrum, M., Omel’chenko, O. E. & Sieber, J. 2015 Regular and irregular patterns of self-localized excitation in arrays of coupled phase oscillators. J. Nonlinear Sci. 25 (5), 053113.Google ScholarPubMed
Figure 0

Figure 1.Sketch in an $(A,Re)$ representation of the stability of the uncontrolled (a) and controlled system using the feedback control method (b) based on the constraint (2.1) corresponding to the dashed slanted line. The labels ‘UB’ and ‘LB’ refer to the stable upper and unstable lower branch in the uncontrolled system, respectively. In (b) the dynamics is restricted to a dashed slanted line. Rotating the line about a pivot point, fixed points on LB are stabilised and can be tracked dynamically.

Figure 1

Figure 2.Application of the feedback control for $\unicode[STIX]{x1D6FC}=1.25$ and $m_{0}=2$. (a) $(Re,A)$ representation, control with $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}\exp (t/T)$, $T=250$ (blue). Cases with $T$ doubled/halved converge to the same stabilised lower branch (grey/cyan). In the inset, the solution set is compared with results from arc-length continuation (using the Newton method, pink) and continuation in the opposite direction using the control with $\unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}\exp (-t/T)$ (green). (b) $A(t)$ uncontrolled (red), controlled with time-varying $\unicode[STIX]{x1D705}(t)$ (blue) and controlled with constant $\unicode[STIX]{x1D705}$ for the stabilisation of the TW (green) (zoom in inset).

Figure 2

Figure 3.(a) Cross-section of the stabilised TW solution with $\unicode[STIX]{x1D6FC}=1.25$ and $m_{0}=2$. Streamwise velocity perturbation in colour (from dark to white) and cross-stream components (vectors). (b) Eigenvalue spectrum for the TW solution with (blue circles) and without control (red crosses).

Figure 3

Figure 4.Application of the feedback control with $L=25D$, $m_{0}=2$ and mirror symmetry imposed. (a) $(Re,A)$ representation. The arrow shows the direction of rotation of the slanted straight line, displayed here are two different times. (b) Time series of $Re(t)$ with time-varying (blue/dark) and constant (green/light) $\unicode[STIX]{x1D705}(t)$ for the stabilisation of the RPO. (c) $(E_{3d},\unicode[STIX]{x1D6FD})$ representation of the RPO converged using the feedback control (red), the unstable RPO of the uncontrolled system (blue) at $Re=1450$, close to each other relative to the laminar state (black circle). (d) Spectrum of Floquet exponents for the RPO solution with (blue circles) and without control (red crosses).

Figure 4

Figure 5.Application of the feedback control with $\unicode[STIX]{x1D6FC}=1.25$, $m_{0}=1$. (a) $(Re,A)$ representation. The arrow shows the direction of rotation of the slanted straight line, which is displayed at the time where an approach to an unstable TW is identified. (b) Cross-section of the stabilised TW solution (same colour coding as in figure 3), $Re=2680$. (c) Eigenvalue spectrum for the TW solution with (blue circles) and without control (red crosses). (d) Cross-section after arc-length continuation to $Re=1428$, $\unicode[STIX]{x1D6FC}=1.25$.