Stabilising nonlinear travelling waves in pipe flow using time-delayed feedback

Abstract

We demonstrate the first successful non-invasive stabilisation of nonlinear travelling waves in a straight cylindrical pipe using time-delayed feedback control working in various symmetric subspaces. By using an approximate linear stability analysis and by analysing the frequency-domain effect of the control using transfer functions, we find that solutions with well-separated unstable eigenfrequencies can have narrow windows of stabilising time delays. To mitigate this issue we employ a ‘multiple time-delayed feedback’ approach, where several control terms are included to attenuate a broad range of unstable eigenfrequencies. We implement a gradient descent method to dynamically adjust the gain functions in order to reduce the need for tuning a high-dimensional parameter space. This results in a novel control method where the properties of the target state are not needed in advance, and speculative guesses can result in robust stabilisation. This enables travelling waves to be stabilised from generic turbulent states and unknown travelling waves to be obtained in highly symmetric subspaces.

1. Introduction

In this study we consider the incompressible Navier–Stokes equations as a high-dimensional dynamical system, where simple invariant solutions, or exact coherent structures (ECSs), can be considered key building blocks of the spatio-temporal chaos (see, e.g. Kawahara, Uhlmann & van Veen 2012; Graham & Floryan 2021). Exact coherent structures take the form of equilibria, relative equilibria (or travelling waves), periodic orbits, relative periodic orbits and invariant tori. The study of ECS has shed light on the origin of turbulence statistics (Chandler & Kerswell 2013; Lucas & Kerswell 2015; Page et al. 2024b), the physical mechanisms that play a role in sustaining turbulence (van Veen, Kida & Kawahara 2006; Lucas & Kerswell 2017; Yasuda et al. 2019; Graham & Floryan 2021; McCormack, Cavalieri & Hwang 2024), subcritical transition to turbulence (Skufca, Yorke & Eckhardt 2006; Kreilos & Eckhardt 2012), mixing and layer formation in stratified flow (Lucas & Caulfield 2017; Lucas, Caulfield & Kerswell 2017), pattern formation in convection (Beaume, Bergeon & Knobloch 2011; Reetz & Schneider 2020) and drag reduction (Bengana et al. 2022) amongst other applications. The fundamental objective is to be able to describe a turbulent flow by using ECSs as a reduced-order description and as a simple way to predict its statistics. Exact coherent structures in various fluid flows have been successfully isolated using homotopy (Nagata 1990), bisection (Itano & Toh 2001; Duguet, Willis & Kerswell 2008b) and quasi-Newton shooting (Viswanath 2007, 2009; Chandler & Kerswell 2013). Newton–Krylov ‘shooting’ methods have proved the most successful and efficient methods for converging these unstable states from the turbulent attractor; however, as the Reynolds number or system size is increased, and the flow becomes increasingly disordered, well-conditioned guesses become harder to identify. Recent research has sought to address this issue, generally by identifying better starting guesses (e.g. using dynamic mode decomposition (Page & Kerswell 2020; Marensi et al. 2023), convolutional autoencoders (Page et al. 2024a) or improved recurrence conditions (Redfern, Lazer & Lucas 2024)) and/or by ‘preprocessing’ guesses by gradient descent (e.g. using automatic differentiation (Page et al. 2024b) to obtain derivatives). An alternative to shooting is to converge closed loops using variational methods to minimise cost functions (Azimi, Ashtari & Schneider 2022; Parker & Schneider 2022), which ensures the governing equations are satisfied. This approach has a larger radius of convergence but is much slower to converge than Newton shooting.

Here, our approach will be to dispense with an iterative or root-finding approach and instead control the underlying ECSs, thereby allowing one to simply time step a modified set of equations onto them. We will develop a time-delayed feedback control method to stabilise travelling waves from turbulence in pipe flow. This method serves to complement existing methods to find ECSs and provides new insight into flow control in general.

A relatively rich literature on ECSs in pipe flow exists, with many travelling-wave solutions (Faisst & Eckhardt 2003; Wedin & Kerswell 2004; Pringle, Duguet & Kerswell 2009; Viswanath 2009; Willis, Cvitanović & Avila 2013; Ozcakir et al. 2016; Ozcakir, Hall & Tanveer 2019) and relative periodic solutions (Duguet, Pringle & Kerswell 2008a; Budanur et al. 2017) reported. Moreover experimental observations suggest that these solutions to the governing equations, even with idealised boundary conditions (periodicity in the streamwise direction), do have relevance to real-world applications (Hof et al. 2006). Feedback control has been used to great effect by Willis et al. (2017) to obtain ‘edge states’ in pipe flow. However, when applied to channel flow, Linkmann et al. (2020) discovered that such control can induce new instabilities, even when the original unstable mode is stabilised. It remains to be seen if more generalised control methods can obtain other ECSs by stabilising more than one direction. For these reasons, pipe flow is an excellent candidate system in which to test and develop this control approach.

The method of time-delayed feedback (TDF) control, sometimes called Pyragas control (Pyragas 1992), is a well-known approach for stabilising invariant solutions in chaotic systems and has been used to great effect in a variety of dynamical systems (Lüthje, Wolff & Pfister 2001; Ushakov et al. 2004; Popovych, Hauptmann & Tass 2005; Yamasue & Hikihara 2006; Stich, Casal & Beta 2013). A finite-dimensional autonomous dynamical system with state vector $\boldsymbol {X}(t)$ can be expressed as

(1.1)\begin{equation} \frac{\mathrm{d}\boldsymbol{X}}{\mathrm{d}t} = \boldsymbol{f}(\boldsymbol{X}(t); \boldsymbol{p}) + \boldsymbol{F} , \end{equation}

where, in our case, the vector field $\boldsymbol {f}$ is given by the discretised Navier–Stokes equation with parameters $\boldsymbol {p}$; we denote $\boldsymbol {F}$ as the TDF ‘force’,

(1.2)\begin{equation} \boldsymbol{F}(t) = G(t) [\boldsymbol{X}(t-\tau) - \boldsymbol{X}(t)] , \end{equation}

where $\boldsymbol {X}(t)$ and $\boldsymbol {X}(t-\tau )$ are the current state vector and the time-delayed state vector, with delay period $\tau,$ and $G(t)$ is the control gain, which in general could be some matrix (thereby coupling feedback between degrees of freedom) but here is a function of time only. Note that when a time-periodic state with period $\tau,$ or a time-independent state, is stabilised successfully, the control force (1.2) will decay exponentially towards zero (assuming $G(t)$ remains bounded). This can only occur if the stabilised state is itself a solution of the uncontrolled system, so that TDF is considered to be a non-invasive control method.

One reason for the lack of application of TDF in fluid systems is likely due to the so-called ‘odd-number’ limitation. This claims that states with an odd number of unstable Floquet multipliers are unable to be stabilised by this method (Nakajima & Ueda 1998b; Just et al. 1999). Nakajima (1997) explains the issue from a bifurcation perspective. The non-invasive feature means that the number of solutions of period $\tau$ cannot vary with $G$; however, any stabilisation requires a change of stability and hence bifurcation. Such a bifurcation cannot, therefore, be of a pitchfork or saddle-node type (changing $\tau$-period solutions), and so must involve the crossing of a complex conjugate pair of exponents in a Hopf or period-doubling bifurcation. There have been numerous studies offering resolutions to this issue, including forcing oscillation of the unstable manifold through $G$ (Schuster & Stemmler 1997; Flunkert & Schöll 2011), an ‘act-and-wait’ approach (Pyragas & Pyragas 2019), using time-dependent gains more generically (Sieber 2016), using symmetries (Nakajima & Ueda 1998a; Lucas & Yasuda 2022) and a notable counter-example where a transcritical bifurcation can be stabilised under certain conditions (Fiedler et al. 2011).

Shaabani-Ardali, Sipp & Lesshafft (2017) report the application of Pyragas control to suppress vortex pairing in a periodically forced jet. This work approaches the control method as a frequency damping technique (Akervik et al. 2006), filtering out non-harmonic frequencies, leaving only $\tau$ behind. Here the odd-number issue can be viewed as a zero-mode limit where there is no incipient frequency to damp. Recently, Lucas & Yasuda (2022) applied this method to two-dimensional turbulent Kolmogorov flow, validating the stabilisation of the base flow via linear stability analysis and showing successful stabilisation of several equilibria and travelling waves. This was achieved by including the symmetries of the target solutions into the control force. In many examples, including the stabilisation of the laminar solution, this was an effective means to avoid the odd-number problem (not too dissimilar to, but different from, the ‘half-period’ approach in Nakajima & Ueda 1998a). An adaptive, gradient descent, approach is also used to obtain the relative translations of travelling waves so that the method can be successful without any foreknowledge of the ECS.

Our main objective in this paper is to present an improved method for using TDF to stabilise nonlinear travelling waves in pipe turbulence, where the flows are more physically relevant, more unstable and have more spatial complexity.

The vast majority of TDF research has been devoted to investigating the linear stability of a target solution, it being a necessary property for successful stabilisation (Nakajima 1997). However, this is not a sufficient condition as successful practical stabilisation can also depend on the initial conditions. In high-dimensional systems, it would be advantageous to design the control to also maximise the basin of attraction of the stabilised state. Furthermore it will be necessary for the TDF control to stabilise several unstable directions, each with differing eigenfrequencies, without knowing what these frequencies are before initiating the control. To achieve this, we develop an adaptive version of the multi-frequency damping TDF control method (Ahlborn & Parlitz 2004), i.e. multiple time-delayed feedback (MTDF), where several TDF terms are used. One drawback of the MTDF approach is the need to optimise a separate control gain for each delay period applied. In order to avoid a trial-and-error sampling of this high-dimensional parameter space, we apply a gradient descent method (Lehnert et al. 2011) to evolve each $G$ towards its stabilising value.

This paper is organised as follows. In § 2, we describe the governing equations for pipe flow with the control force, the numerical method and the continuous and discrete symmetries and define relevant flow measures. In § 3, after introducing single TDF for pipe flow and the adaptive translation method, we demonstrate some successful cases stabilising an unstable travelling wave at low Reynolds numbers and in certain symmetric subspaces. We predict the behaviour of TDF by using an approximate linear stability analysis and control theory in such a way as to identify the optimal control parameters for this solution. In § 4, we introduce MTDF control, demonstrate its effectiveness in stabilising more highly unstable states and analyse its behaviour from the frequency damping perspective. In § 5, we demonstrate successful stabilisation of travelling waves from relatively high Reynolds number turbulence, using MTDF alongside the optimisation methods for gains and translations. In the course of doing this, we demonstrate the stabilisation of two unknown solutions from highly symmetric subspaces.

2. Numerical formulation

2.1. Pipe flow with time-delayed feedback

In this study, we consider incompressible, viscous flow in straight, cylindrical, pressure-driven pipes. We treat the governing equations in cylindrical-polar coordinates $(r,\theta,z)$, where $r$ is the radius, $\theta$ is the azimuthal angle and $z$ is the streamwise (axial) position. We non-dimensionalise the equations with the Hagen–Poiseuille (HP) centreline speed, $U_{cl}$, and pipe radius, $R$. This yields the dimensionless incompressible Navier–Stokes equations

(2.1)$$\begin{gather} \frac{\partial}{\partial t} \boldsymbol{U} + (\boldsymbol{U} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{U} ={-} \boldsymbol{\nabla} P + \frac{1}{Re}\nabla^2 \boldsymbol{U} +\boldsymbol{F}, \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U} = 0 , \end{gather}$$

with the no-slip condition on the boundary

(2.3)\begin{equation} \boldsymbol{U}(1,\theta,z) = \boldsymbol{0} . \end{equation}

Here, $\boldsymbol {U}=(U_r,U_\theta,U_z)$ is the three-dimensional velocity vector, with $P$ being the pressure; $Re = R U_{cl}/\nu$ is the Reynolds number, $\nu$ is the kinematic viscosity of the fluid and time $t$ is defined in the unit of $R/U_{cl}$. Lastly, $\boldsymbol {F}$ is an external body force that, here, will include the TDF control terms. With the aforementioned non-dimensionalisation, the laminar HP flow is expressed by means of velocity and pressure as

(2.4)$$\begin{gather} \boldsymbol{U}_{HP}(r) = (1-r^2) \hat{\boldsymbol{z}} , \end{gather}$$

(2.5)$$\begin{gather}P_{HP}(z) ={-} 4 z/Re , \end{gather}$$

where $\hat {\boldsymbol {z}}$ is the unit vector in the streamwise direction. The laminar HP flow is found at low Reynolds numbers and is linearly stable even at very large (or possibly infinite) Reynolds numbers (Salwen, Cotton & Grosch 1980; Meseguer & Trefethen 2003).

Using (2.4) and (2.5), the velocity and pressure fields can be decomposed as

(2.6)$$\begin{gather} \boldsymbol{U}(r,\theta,z,t) = \boldsymbol{U}_{HP} (r) + \boldsymbol{u}(r,\theta,z,t), \end{gather}$$

(2.7)$$\begin{gather}P(r,\theta,z,t) = P_{HP} (z) + p(r,\theta,z,t) , \end{gather}$$

where $u=(u_r , u_\theta, u_z)$ and $p$ are the velocity and pressure deviations from the laminar fields, respectively. Fundamental studies of pipe flow typically consider two driving mechanisms: a constant mass flux (Darbyshire & Mullin 1995; Duguet et al. 2008a; Willis et al. 2013) or a constant pressure gradient (Wedin & Kerswell 2004; Shimizu & Kida 2009). In what follows we will seek comparisons with ECSs enumerated in Willis et al. (2013); therefore, we choose to use the same constant mass flux formulation. To consider this formulation, we further decompose the pressure deviation, $p$, as $p = \hat {p}(r,\theta,z) + \zeta (t) z$, leading to the expression $\boldsymbol {\nabla } p = \boldsymbol {\nabla } \hat {p} + \zeta (t) \hat {\boldsymbol {z}}$ (Marensi et al. 2020). Here, $\zeta (t) = - 4 \beta (t) / Re$ represents an additional pressure gradient necessary to maintain a constant mass flux at all times, whereas $\boldsymbol {\nabla }\hat {p}$ denotes a pressure gradient that has no mean streamwise component. The dimensionless variable, $\beta (t)$, is an additional pressure fraction that can be determined in experiments using the following equation:

(2.8)\begin{equation} 1 + \beta(t) = \frac{\langle \partial P / \partial z \rangle_V}{\langle \partial P_{HP} / \partial z \rangle_V}, \end{equation}

where

(2.9)\begin{equation} \langle ({\cdot}) \rangle_V := \int_0^L \int_0^{2{\rm \pi}} \int_0^1 ({\cdot}) r\,\mathrm{d}r \,\mathrm{d}\theta \,\mathrm{d} z . \end{equation}

In our simulations, we compute $\beta (t)$ through the spatially integrated force balance equation between pressure and viscous wall shear stress in the streamwise direction

(2.10)\begin{equation} \beta(t) ={-}\frac{1}{2}\!\left.\frac{ \langle \partial {u_z}\rangle_{\theta,z}}{\partial r}\right|_{r=1}, \end{equation}

where

(2.11)\begin{equation} \langle{({\cdot})}\rangle_{\theta,z} := \frac{1}{2 {\rm \pi}L}\int_0^L \int_0^{2{\rm \pi}} ({\cdot}) \,\mathrm{d} \theta \,\mathrm{d} z . \end{equation}

Note that there is no contribution from an external body force in (2.10) since $\boldsymbol {F}$ has no mean streamwise component in this study (cf. Marensi et al. 2020). Finally, the total pressure gradient, $\boldsymbol {\nabla }P$, can be transformed using (2.5) and $\beta (t)$ as

(2.12)\begin{align} \boldsymbol{\nabla} P &= \boldsymbol{\nabla} P_{HP} + \zeta(t)\hat{\boldsymbol{z}} + \boldsymbol{\nabla} \hat{p} \end{align}

(2.13)\begin{align} &={-}\frac{4}{Re}(1+\beta(t))\hat{\boldsymbol{z}} + \boldsymbol{\nabla} \hat{p}. \end{align}

2.2. Direct numerical simulations

We solve the system outlined in § 2.1 numerically using the open-source code openpipeflow (Willis 2017) which allows the relatively easy implementation of our control terms. Our computational domain is periodic in both azimuthal and streamwise directions, and the streamwise length of the periodic pipe is $L=2{\rm \pi} /\alpha$, e.g. with $\alpha =1.25$ corresponding to $L\approx 5 R$. Here, $\boldsymbol {u}$ and $\hat {p}$ are both expanded in discrete Fourier series in the streamwise and azimuthal directions. Spatial derivatives with respect to $r$ are evaluated based on a nine-point finite-difference stencil with a non-uniform mesh, for which first-/second-order derivatives are calculated to eighth/seventh order (Willis 2017). With these spatial discretisation schemes, $\boldsymbol {u}$ is expressed as

(2.14)\begin{equation} \boldsymbol{u}(r_n,\theta,z,t) = \sum_{|k|< K} \sum_{|m|< M} \tilde{\boldsymbol{u}}_{km}(r_n,t) \exp({{\rm i}(\alpha k z+ m \theta)}), \end{equation}

where $\tilde {\boldsymbol {u}}_{km}$ are the Fourier coefficients of $\boldsymbol {u},$ $r_n$ ($n=1 \ldots N$) denotes the radial grid points (non-uniformly distributed on [0, 1]) and $k$ and $m$ are the streamwise and azimuthal wavenumbers, respectively, with $K$ and $M$ being the de-aliasing cutoff (see Willis 2017). The resolution of a given calculation is described by a vector $(N, M, K)$ and is adjusted until the energy in the spectral coefficients drops by at least five and usually six decades. Time stepping is executed via a second-order predictor–corrector scheme, with the Euler predictor method for the nonlinear terms and the Crank–Nicolson method for the viscous diffusion.

2.3. Continuous and discrete symmetries

The equations of pipe flow are invariant under continuous translations in $z$

(2.15)\begin{equation} \mathcal{T}_z(s_z) [u_r,u_\theta,u_z,p](r,\theta,z) \rightarrow [u_r,u_\theta,u_z,p](r,\theta,z+s_z) , \end{equation}

continuous rotations in $\theta$

(2.16)\begin{equation} \mathcal{T}_\theta(s_\theta) [u_r,u_\theta,u_z,p](r,\theta,z) \rightarrow [u_r,u_\theta,u_z,p](r,\theta+s_\theta,z) , \end{equation}

and reflections about $\theta =0$

(2.17)\begin{equation} \sigma[u_r,u_\theta,u_z,p](r,\theta,z) \rightarrow [u_r,-u_\theta,u_z,p](r,-\theta,z) . \end{equation}

Following the approach by Willis et al. (2013), we will restrict our investigations to the dynamics restricted to the ‘shift-and-reflect’ symmetric subspace, $S=\sigma \mathcal {T}_z ({L}/{2} )$

(2.18)\begin{equation} S[u_r,u_\theta,u_z,p](r,\theta,z) \rightarrow [u_r,-u_\theta,u_z,p](r,-\theta,z-L_z/2), \end{equation}

and the discrete ‘rotate-and-reflect’ symmetry $Z_{m_p}=\sigma \mathcal {T}_\theta ({{\rm \pi} }/{m_p} )$

(2.19)\begin{equation} Z_{m_p}[u_r,u_\theta,u_z,p](r,\theta,z) \rightarrow [u_r,-u_\theta,u_z,p](r,{\rm \pi}/m_p-\theta,z) . \end{equation}

Note that $m_p=1$ denotes azimuthal periodicity (the full space), and $m_p$-fold rotational symmetry is enforced for $m_p\geqslant 2$ (Wedin & Kerswell 2004; Pringle et al. 2009; Willis et al. 2013). By imposing $S$, the flow is ‘pinned’ in $\theta$ and continuous rotations are prohibited. This means that our solutions are only permitted to travel in the streamwise direction due to $\mathcal {T}_z.$

2.4. Flow measures

In order to monitor the system behaviour, we consider the spatially integrated quantities from the energy budget equation

(2.20)\begin{equation} \frac{\mathrm{d} E}{\mathrm{d} t} = I - D + I_{TDF} , \end{equation}

where $E$ is the total kinetic energy, $I$ is the total energy input due to the imposed pressure gradient and $D$ is the total energy dissipation

(2.21a–c)\begin{equation} E := \frac{1}{2} \int |\boldsymbol{U}|^2 \,\mathrm{d}V ,\quad I := \int \boldsymbol{U}\boldsymbol{\cdot}(-\boldsymbol{\nabla} P)\,\mathrm{d}V ,\quad D := \frac{1}{Re} \int | \boldsymbol{\nabla} \times \boldsymbol{U} |^2 \,\mathrm{d}V , \end{equation}

and $I_{TDF}$ is the total energy input due to feedback force $\boldsymbol {F}$

(2.22)\begin{equation} I_{TDF} := \int \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{F} \,\mathrm{d}V . \end{equation}

For the laminar state, $I_{TDF}=0$ and the values of energy, input rate and dissipation rate, $E_{lam}$, $I_{lam}$ and $D_{lam}$ can be computed from (2.21a–c) with $\boldsymbol {U}=\boldsymbol {U}_{HP}$.

3. Time-delayed feedback

3.1. Formulation

In this section we outline the application of the TDF control (1.2) to the pipe flow. In order to effectively stabilise travelling-wave solutions, the delayed state must be translated by a streamwise shift $s_z$ in relation to the phase speed, $c_z$, such that $c_z={s_z}/{\tau }$ for a given time delay $\tau.$ Using the translation operator (2.15), the most basic form of TDF force in pipe flow may be formulated as

(3.1)\begin{equation} \boldsymbol{F}_{TDF}(r,\theta,z,t) = G(t) [ \mathcal{T}_z(s_z) \boldsymbol{u}(r,\theta,z,t-\tau) - \boldsymbol{u}(r,\theta,z,t)] , \end{equation}

where $\tau$ is the delay period, and $G(t)$ is the gain. Note that the gain here is a simple scalar function of time, but it may itself be an operator (matrix) or spatially inhomogeneous. In the full space a rotation $\mathcal {T}_{\theta }$ may also be applied. In order to avoid a discontinuity propagating through the solution when the control is initiated, we set $G(t)$ using the following sigmoid function:

(3.2)\begin{equation} G(t) = \begin{cases} \dfrac{G^{max}}{1+\exp[a(b+t_s-t)]} , & t>t_s\\ 0 & t\leqslant t_s. \end{cases} \end{equation}

Here, $t_s$ is the time TDF is initiated, $G^{max}$ is the maximum gain, $a$ determines the slope of $G(t)$ and the half-height time, $t_h=t_s+b,$ set by $b,$ is when $G(t_h)=G^{max}/2$.

The value of the translation, $s_z,$ for successful stabilisation will be, in general, unknown in advance, so we require a method to evolve $s_z$ towards the required value. We implement the adaptive method developed in Lucas & Yasuda (2022). This method allows $s_z$ to vary via gradient descent by solving an ordinary differential equation

(3.3)\begin{equation} \frac{\mathrm{d} s_z}{\mathrm{d} t} = \gamma_s \delta s_z, \end{equation}

where $\gamma _s$ is a parameter controlling the speed of the descent and $\delta s_z$ is, near a travelling-wave state, an estimate of the streamwise translation remaining between the current state and the delayed and translated state. Specifically, $\delta s_z$ is computed as

(3.4)\begin{equation} \delta s_z = s_z^{est} - s_z , \end{equation}

where $s_z^{est}$ is a streamwise translation that is dynamically estimated. Here, $s_z^{est}$ can be computed via the time series of $c_z(t)$ such that

(3.5)\begin{equation} s_z^{est}(t;\tau) = \int_{t-\tau}^{t} c_z(t^{\prime}) \,\mathrm{d} t^{\prime} . \end{equation}

We can compute $c_z(t)$ using complex phase rotations such that

(3.6)\begin{equation} c_z(t) = \frac{1}{N_{e}} \sum_{k={\pm} 1}\sum_{m={\pm} 1}\sum_{n=1}^{N-1} \frac{1}{{\rm i}\alpha k \Delta t} \arg \left[\frac{\tilde{u}_{k m}(r_n, t)}{\tilde{u}_{k m}(r_n, t-\Delta t)} \right], \end{equation}

where $N_{e}=4(N-1)$ is the number of the elements in the above summation and the wall points ($n=N$) are excluded. Note that $c_z(t)$ is considered an average phase speed over one time step, $\Delta t$, in the streamwise direction. We solve the ordinary differential equation (ODE) (3.3) while computing $s_z^{est}$ using (3.5) alongside the direct numerical simulation with a second-order Adams–Bashforth time-stepping scheme. As will be seen in the later sections, using long time delays can be essential in successful stabilisation of unstable travelling waves in a straight cylindrical pipe. Computing estimated translations using (3.5) is slightly different from the original version in Lucas & Yasuda (2022). The current version is found to perform better than the original when encountering long delays.

3.2. Validation – stabilising weakly unstable solutions

In the absence of any simpler solutions to use for validation (the laminar HP flow is linearly stable at least up to $Re=10^7$, see Meseguer & Trefethen 2003), we attempt to stabilise known unstable travelling waves (Willis et al. 2013) at low $Re$ first. For the parameter values of $\alpha =1.25$ ($L_z \approx 5R$) and $Re=2400$, we are able to stabilise the travelling waves ML, UB, S2U in their respective symmetric subspace (see table 1). These travelling waves have only complex unstable eigenvalues at $\alpha =1.25$ and $Re=2400$ and in their symmetric subspaces, meaning that the odd-number issue is not applicable, or rather is avoided by projection into those subspaces (see Lucas & Yasuda 2022, for a discussion). In order to investigate the effectiveness of TDF as the stability of a target solution changes, we concentrate on the UB solution (in $(S,Z_2)$) and vary $Re,$ returning to the other solutions in § 5.3. Figure 1 shows the bifurcation diagram with the continuation of the UB solution branch (blue) as a function of $Re$ (alongside ML and S2U). This is generated numerically using the Newton–generalised minimal residual method (GMRES)–hookstep method (Viswanath 2007, 2009) packaged with the openpipeflow code (Willis 2017). This confirms that UB originates in a saddle-node bifurcation at $Re\approx 1150$. On increasing $Re$, the stable upper branch (solid line in figure 1) of UB becomes unstable via a supercritical Hopf bifurcation at $Re\approx 1800$ (unstable solution in dashed blue).

Table 1. Table summarising the properties of all travelling waves studied in this paper. Here, $Re = R U_{cl}/\nu$ is the Reynolds number, $\alpha$ is the streamwise wavenumber ($L=2{\rm \pi} /\alpha$), $E$ is the total kinetic energy, $I$ is the total energy input rate, $Re_{\tau }=\sqrt {2 Re I/I_{lam}}$ is the friction Reynolds number, $c_z$ is the streamwise phase velocity, ‘sym.’ refers to the symmetric subspace, $\lambda _i = \mu _i \pm \mathrm {i} \omega _i$ are the eigenvalues of the solutions such that the final three columns are, respectively, the number and type of unstable eigenvalues, the most unstable growth rate and the smallest unstable eigenfrequency.

Figure 1. Bifurcation diagram showing solution branches for flows with $\alpha =1.25$.

In figure 2, we plot the real and imaginary parts of the complex growth rates of the UB solution, which demonstrates the expected increase in instability as the Reynolds number increases. A sequence of further bifurcations is observed, and new unstable directions are formed as more eigenvalues cross the imaginary axis. The solution has only complex eigenvalues at least up to $Re=3050$; therefore, it is free from the odd-number limitation (Nakajima 1997).

Figure 2. The $Re$-dependence of leading eigenvalues of the UB solution. (a) Real part and (b) imaginary part of unstable complex eigenvalues. No purely real eigenvalue exists; thus the ‘odd-number limitation’ is not encountered.

Using this stability information we attempt to stabilise UB at various $Re$ using TDF. The initial condition is UB such that $\boldsymbol {U}(r,\theta,z,0)=\boldsymbol {U}_{\mathsf {UB}}$. We set $t_s=10$ such that the trajectory is still close to the UB solution when TDF is activated, avoiding the need to consider if our initial condition falls inside the basin of attraction for now. In this sense this section verifies a necessary condition for stabilisation and we will investigate more generic initial conditions later. In all cases where UB is stabilised an initial $c_z(0)=0.65$ is used.

In order to quantify the size of the TDF force term, we define the relative residual, $R_{tot}$, between the current state and the translated and delayed state, as

(3.7)\begin{equation} R_{tot}= \frac{\|\mathcal{T}_z(s_z)\boldsymbol{U}(r,\theta,z,t-\tau)-\boldsymbol{U}(r,\theta,z,t)\|_2} {\|\boldsymbol{U}(r,\theta,z,t)\|_2}, \end{equation}

where $\| \boldsymbol {A} \|_2 = \sqrt {\frac {1}{2} \langle \boldsymbol {A}\boldsymbol {\cdot }\boldsymbol {A} \rangle _V}.$ Figure 3(a) shows time series in $R_{tot}$ for seven Reynolds numbers, $Re=1900$, 2000, 2100, 2200, 2300, 2400, 2500, with $\tau =2,$ $G^{max}=0.5,$ $\gamma _s=0.1,$ $t_s=10,$ $a=0.1,$ $b=100.$ At $Re=1900$, stabilisation of UB is achieved quickly such that $R_{tot}$ has decreased to $O(10^{-12})$ by $t=2000$. The rate of attraction of the travelling wave decreases as $Re$ increases with all other parameters held fixed, as one might expect.

Figure 3. Time series in the relative residual, $R_{tot}(t;\tau )$ for attempted TDF stabilisation of UB. (a) Shows $Re=1900,2000,2100,2200,2300,2400,2500$ for TDF parameters $\tau = 2$ and $G^{max}=0.5,$ stabilisation fails for $Re>2400$. (b) Shows the $Re=2500$ case with $G^{max}=0.1$ and 0.5, and $\tau =2$ (both of which fail to stabilise) and $\tau =5$ (which stabilises successfully), demonstrating consistency with the linear analysis. The simulations are initiated with UB at each Reynolds number. Note that recording $R_{tot}(t;\tau )$ initiates at $t=\tau$.

At $Re=2500$, stabilisation is unsuccessful; even by increasing $G^{max}$, we are unable to stabilise the UB state at this $Re,$ without adjusting the other TDF parameters. It turns out that the growth rate is not the key quantity preventing stabilisation in this case, and neither is it the appearance of an odd-numbered eigenvalue (in contrast to the findings in Lucas & Yasuda 2022). Here, UB at $Re=2500$ has three unstable eigenfrequencies, which are decreasing as a function of $Re;$ it is these values which fall outside the domain of stabilisation of our method. In the next subsection, we analyse the TDF control from the perspective of ‘frequency damping’ (Akervik et al. 2006; Shaabani-Ardali et al. 2017) and linear stability analysis.

3.3. Frequency damping and linear stability with TDF

Predicting the outcome of TDF in these scenarios is difficult due to the relatively high dimension of the unstable manifold (of the uncontrolled solution) and the highly nonlinear nature of the solution. Nevertheless we can approximate the effect of the TDF terms on the unstable part of the eigenvalue spectrum. Starting from (1.1)–(1.2), assuming $\boldsymbol {X} = \bar {\boldsymbol {X}} + \boldsymbol {v}{\rm e}^{\lambda t},$ with $\bar {\boldsymbol {X}}$ a steady-state solution to (1.1), assuming $G$ is now constant and linearising in the usual way results in the modified eigenvalue problem

(3.8)\begin{equation} \lambda \boldsymbol{v} = \mathcal{J}\boldsymbol{v} + G\boldsymbol{v}\left({\rm e}^{-\lambda \tau} - 1\right), \end{equation}

with eigenvalue $\lambda,$ eigenvector $\boldsymbol {v}$ and $\mathcal {J}$ being the Jacobian of $\boldsymbol {f}.$ The eigenvalue problem is now transcendental and so requires some numerical root searching to obtain the eigenvalues. Lucas & Yasuda (2022) tackled this difficulty by assuming $|\lambda \tau | \ll 1$ and expanding the exponential, however, we are unable to make this assumption here. Moreover, root searching for $\lambda$ in the full problem, i.e. some numerical linearisation of (2.1), would be impractical. Instead, we create a ‘synthetic’ $\mathcal {J},$ which we design to have the same unstable eigenvalues as the uncontrolled solution and some randomly chosen eigenvectors. For instance, at $Re=2500$, with three pairs of unstable, complex eigenvalues, this amounts to a six-dimensional real-valued system. We can then solve the transcendental characteristic equation numerically with a simple Newton search to obtain all of the roots and hence eigenvalues. This procedure is able to predict well the effect of TDF on the UB solution; across the parameters shown in figure 3(a), and using the uncontrolled eigenvalues shown in figure 2, the analysis finds that the largest eigenvalue crosses the imaginary axis at $Re \approx 2450.$ Moreover the method is sufficiently efficient that we can explore the $(G,\tau )$ parameter space. Figure 4 shows two plots of the largest real part of the eigenvalues of our approximated TDF system at $Re=2500$. We observe that for $\tau = 2$, no value of $G$ results in stability. However, increasing slightly to $\tau =5$ results in a region of stability for relatively small $G$; we see $G\approx 0.1$ is roughly optimal. This is slightly surprising as one would have assumed that as solutions become more unstable, the remedy would be to apply larger gains to move eigenvalues back across the imaginary axis. It should be noted that this analysis is only an approximation to the problem; in particular, it does not predict whether TDF may destabilise stable modes (while it is possible for TDF to destabilise stable modes or introduce new modes of instability, as found in Linkmann et al. 2020, this is never observed when the uncontrolled modes are also stabilised and only for large $G$ and $\tau$).

Figure 4. Dependence of the largest real part of the eigenvalue spectrum $\max _i \mu _i$ ($\lambda _i = \mu _i + \mathrm {i} \omega _i,$) with (a) $G$ and (b) $\tau,$ for the $Re=2500$ UB solution, according to the approximate linear theory.

We can verify these predictions by performing some additional numerical simulations with TDF in the full Navier–Stokes equations. Figure 3(b) demonstrates that stability is observed for $\tau =5$ and $G^{max}=0.1.$ Moreover $G^{max}=0.5$ is less strongly stable, and $\tau =2$, $G^{max}=0.1$ is more strongly unstable, all of which is consistent with the above analysis. While this stability analysis is useful in validating and optimising TDF, it does not offer much insight into, for instance, why $\tau \approx 5$ is optimal in this example. Following the ‘frequency damping’ interpretation of TDF (Akervik et al. 2006; Shaabani-Ardali et al. 2017), we define the transfer function

(3.9)\begin{equation} H_{TDF}(\mathrm{i}\omega,\tau) = \frac{1}{G} \frac{\mathcal{L}\{\boldsymbol{F}_{TDF}\}}{\mathcal{L}\{\boldsymbol{u}\}} = {\rm e}^{-\mathrm{i}\omega \tau}-1 , \end{equation}

with $\mathcal {L}\{.\}(\omega ) = \int _0^\infty. \,{\rm e}^{-\omega t}\,{\rm d}t$ the Laplace transform. The Laplace transform of the control term is $\mathcal {L}\{\boldsymbol {F}_{TDF}\}(\omega ) = -G(1-{\rm e}^{-\omega \tau }) \mathcal {L}\{ \boldsymbol {u} \}.$ Here, $H_{TDF}(\mathrm {i}\omega )$ allows us to ascertain the relative influence of the control term upon the temporal frequencies of the system in general (unlike a specific linear stability analysis). Figure 5(a) plots the magnitude of the transfer function against a normalised angular frequency. The first observation is that any zero mode is undamped by TDF control; this is an alternative explanation of the odd-number limitation where a purely real-valued eigenvalue cannot be stabilised in isolation. The next observation is that any harmonic of the feedback frequency ${2{\rm \pi} }/{\tau }$ is also unaffected by the control. The consequence of these facts means that any unstable eigenmode of the target ECS with a frequency that is either close to zero, or close to the feedback frequency (or harmonic thereof), will not be stabilised. For high $Re$ travelling waves, this means that careful tuning of the delay period $\tau$ may be necessary as the likelihood of an unstable eigenfrequency falling near a zero of the transfer function will increase.

Figure 5. Transfer functions for TDF ($\epsilon = 0$) and ETDF (various $\epsilon \neq 0$). (a) Shows $H_{ETDF}(\omega /\omega ^*)$ where $\omega ^* = 2{\rm \pi} /\tau$ demonstrating peaks at subharmonics and zeros at harmonics of $\tau.$ (b) Shows the product of $H_n = H_{ETDF}(\omega _n,\tau )$ for the three unstable eigenfrequencies of UB at $Re=2500.$ This indicates an optimal $\tau,$ for TDF, of around 5.5. This is consistent with the linear analysis shown in figure 4(b).

We may use this analysis to interrogate the behaviour of TDF across a range of delay periods for the UB solution at $Re=2500$ discussed above. Defining $H_n(\tau ) = |H_{TDF}(\mathrm {i}\omega _n,\tau )|$ with $\omega _n$ the unstable eigenfrequencies of the uncontrolled solution, as shown in figure 2(b), namely at $Re=2500,$ $\omega _1 = 0.75$, $\omega _2=0.48$, $\omega _3=0.32,$ we plot the product $H_1H_2H_3$ in figure 5(b). This demonstrates clear agreement with the linear analysis of figure 4 and therefore also the stabilisation of UB in figure 3(b). We see that the product of transfer functions is maximised at $\tau \approx 5$ and distinct zeros at harmonics of the eigenperiods, also coinciding with maxima of $\max (\mu )$ in figure 4(b). More specifically, the first zero of $H_1H_2H_3$ is at $\tau ={2{\rm \pi} }/{\omega _1} = 8.25$ in the same location as the first maximum of $\max (\mu )$ for these parameters, indicating, as expected, TDF will fail for poorly chosen time delays.

3.4. Diagnosing unstable frequencies

We have gained good insight into the behaviour of TDF in this example. However, this analysis relies on knowledge of the unstable spectrum of the target solution. One goal for a fully developed TDF method would be the ability to obtain new, unknown solutions, without any knowledge of their properties in advance. In the example above, we began our investigation under this assumption, i.e. the values of $\tau =2$ and $G^{max}=0.5,$ used in figure 3(a), are set speculatively after a little trial and error.

When $|\lambda \tau | \ll 1,$ we can substitute ${\rm e}^{-\lambda \tau } \approx 1-\lambda \tau$ in (3.8), and we see that the effect of TDF on eigenvalues of $\mathcal {J}$ is a rescaling by $1+G\tau$ (i.e. the eigenvalue equation becomes $(1+G\tau )\lambda \boldsymbol {v}=\mathcal {J}\boldsymbol {v}$). In the example above, since $\tau$ is too small to effectively damp the lowest-frequency eigenvalues, we observe an unstable oscillation with modified frequencies $\omega ' \approx {\omega }/({1+G\tau }).$ In the linear analysis of UB, instability occurs at $Re\approx 2450$ with the $\tau =2$ and $G^{max}=0.5$ TDF term active, with two complex conjugate pairs of eigenvalues crossing the imaginary axis with frequencies $\omega ' \approx 0.16 \approx \omega _2/2$ and $\omega '\approx 0.246 \approx \omega _2/2.$ This means, if we can capture an estimate of these frequencies in the numerical simulation where stabilisation has failed, it would be possible to infer the original, unperturbed frequencies $\omega$. This would enable one to perform the frequency-domain analysis described earlier and choose more carefully tailored $\tau.$

Taking the time series for the total kinetic energy $E(t)$ for the $Re=2500$ case with $\tau =2$ and $G^{max}=0.5,$ as shown in figure 3, and performing a Fourier transform, we observe, in figure 7(b) a low-frequency/zero mode associated with a transient growth/decay of energy, but also two dominant modes at $0.246$ and $0.16,$ as predicted. Similarly we show the case with $\tau =2$ and $G^{max}=0.1,$ which has peaks in the power spectrum at $\omega _2/1.2\approx 0.4$ and $\omega _3/1.2\approx 0.266,$ as we have predicted, but also $\omega _2/1.2-\omega _3/1.2,$ indicating some nonlinear interactions taking place at slightly larger amplitude from the bifurcation. This approach allows one to obtain more effective time delays when TDF fails to stabilise a solution and results in invasive, low-dimensional behaviour, near to a stabilisable state.

4. Generalising time-delayed feedback control

4.1. Extended time-delayed feedback

In the analysis of the previous section we observe relatively narrow operating windows for which TDF will stabilise the target solution, which, as expected, decrease as $Re$ increases. A common way to address this problem is via ‘extended time-delayed feedback control’ (ETDF). This feeds an extended historical record into the control term by including times $t-k\tau$ for increasing integers $k$ as time evolves. For our implementation, this would read, neglecting, for now, translations in $z$

(4.1)\begin{align} F_{ETDF}(t) &= G \left[ (1-\epsilon)\sum_{k=1}^{\infty} \epsilon^{k-1} u(t-k\tau) - u(t) \right] , \nonumber\\ &= G [ u(t-\tau) - u(t) ] + \epsilon F_{ETDF}(t-\tau) , \end{align}

where $0 \leqslant \epsilon \leqslant 1$, see Socolar, Sukow & Gauthier (1994). When $\epsilon =0$ the control reduces back to the standard Pyragas control. This has the effect of broadening and flattening the transfer function

(4.2)\begin{equation} H_{ETDF}(\mathrm{i}\omega,\tau) = \frac{1}{G} \frac{\mathcal{L}\{F_{ETDF}\}}{\mathcal{L}\{u\}} = \frac{{\rm e}^{-\mathrm{i}\omega \tau}-1}{1-\epsilon {\rm e}^{-\mathrm{i}\omega \tau}}; \end{equation}

see figure 5. This method is attractive as it is relatively simple to implement and in practice only requires storing one additional history array for $F_{ETDF}(t-\tau )$ (over a period). However, for our case, in particular for the analysis performed above with UB, we observe that the flattening of the transfer function results in much less distinct maxima, in particular, the peak at $\tau =5.5$ is much smaller (figure 5b) and is unlikely to result in the same stabilisation observed with regular TDF. Moreover, we have seen that attempting to stabilise states with several unstable eigenfrequencies is challenging, and ETDF gives only a marginal advantage in this regard, particularly if the frequencies are well separated.

4.2. Multiple time-delayed feedback

In the previous example for UB the instability observed at $Re=2500$ using $\tau =2$ was remedied following a careful analysis of the effect of TDF terms on the unstable eigenvalues of the uncontrolled system. It was noted that the frequencies are key in choosing an appropriate time delay and that they may be observed in an unsuccessful, invasive TDF dynamics.

To attenuate oscillations across several frequencies, we introduce MTDF, $\boldsymbol {F}_{MTDF}$, such that

(4.3)\begin{equation} \boldsymbol{F}_{MTDF}(r,\theta,z,t) = \sum_i^N G_i(t) [\mathcal{T}_{z}(s_z) \boldsymbol{u}(r,\theta,z,t-\tau_i)-\boldsymbol{u}(r,\theta,z,t)] , \end{equation}

where the $i$th time delay is $\tau _i$ ($i=1,2,\ldots,N$) and we translate the $i$th delayed state, $\boldsymbol {u}(r,\theta,z,t-\tau _i)$, by the streamwise shift, $s_z(t;\tau _i)$. We adjust each streamwise shift for every delay term using the gradient descent method outlined in § 3.1, however, we have verified that adapting the $i=1$ term only and using $s_z(t;\tau _i) = c_z(t)\tau _i$, to match the phases speeds of the other terms is as effective. Delayed feedback with multiple delays was previously used in the low-dimensional chaotic system of Chua's circuit by Ahlborn & Parlitz (2004), to stabilise a nonlinear steady solution. They reported that MTDF helped to expand the basin of attraction of equilibrium solutions, being more efficient than ETDF (Socolar et al. 1994; Sukow et al. 1997).

We can repeat the analysis of § 3.3 with MTDF by introducing a second TDF term into (3.8). We will investigate the unstable $Re=2500$ case, i.e. setting $\tau _1=2,$ $G_1=0.5$ and varying $\tau _2$ and $G_2.$ Given that we have observed $\omega _2/2$ and $\omega _3/2$ in the TDF example, it would be sensible to choose a $\tau _2$ to stabilise these modes. Via the transfer function analysis of the previous section, now with $H_n^* = |H_{TDF}(\mathrm {i}\omega _n/(1+G_1\tau _1),\tau _2)|,$ we are again able to predict an optimal $\tau _2$. The maxima of the product $H_2^*H_3^*$ line up with the minima of $\max (\mu _i)$ from the MTDF linear theory, indicating an optimal $\tau _2\approx 16$, see figures 6(b) and 7(a). The interpretation here is that the second term should act to attenuate the frequencies modified by the first term (or vice versa). It should also be noted that there is now a far larger range of $\tau _2$ and $G_2$ which stabilises UB, meaning less speculative tuning of parameters, particularly of $G_2.$ This result is confirmed in the numerical simulation, figure 8 shows very fast stabilisation of UB at $Re=2500$ with $\tau _1=2$ and $\tau _2=16.$ Success is also observed with $\tau _2=32$ although with a slower rate of attraction, and $\tau _2=40$ shows instability, all of which is consistent with the linear and frequency analysis shown in figures 6 and 7.

Figure 6. Dependence of the largest real part of the eigenvalue spectrum $\max _i \mu _i$ with (a) $G_2$ and (b) $\tau _2$, for the $Re=2500$ UB solution with two-term MTDF and $\tau _1=2$ and $G_1=0.5$. Most effective stabilisation is observed with $\tau _2\approx 16$ and the rather modest $G_2\approx 0.05$.

Figure 7. (a) Shows the product of transfer functions $H_1^*H_2^*$ as defined in the text, for the modified eigenfrequencies when attempting to stabilise UB at $Re=2500$ using $\tau =2$ and $G=0.5.$ Note that the peak coincides with figure 6(b), i.e. $\tau _2\approx 16$ is around the optimal. (b) Shows the power spectrum of energy, $\hat {E}$, for the unsuccessful TDF cases at $Re=2500$ with $\tau =2$ and $G = 0.5$ and 0.1. Vertical lines show that the peaks in these spectra follow the expected scaling by $1/(1+\tau G).$

Figure 8. The MTDF stabilisation of the UB travelling wave at $Re=2500$, shown via time series of (a) $E/E_{lam}$ and (b) $R_{tot}$. Comparison shows various choices of $\tau _2$ with $\tau _1=2,$ $G_1=0.5$ and $G_2=0.1.$ We see rapid stabilisation for the optimal $\tau _2=16,$ moderate stabilisation for $\tau _2=32$ and instability for $\tau =40$.

Note one could define a transfer function $H_{MTDF},$ for MTDF, given by

(4.4)\begin{equation} H_{MTDF}(\mathrm{i}\omega) = \frac{1}{\sum_i^N G_i} \frac{\mathcal{L} \{F_{MTDF}\}}{\mathcal{L} \{u\}} = \frac{\sum_i^N G_i\exp(-{\rm i}\omega\tau_i)}{\sum_i^N G_i}-1. \end{equation}

However, applying this to the case under discussion effectively yields a repeat of the single TDF result, in that the optimal $\tau _2$ to stabilise the three eigenfrequencies, with $\tau _1=2,$ is around 5; $H_{MTDF}$ does not accurately predict the effect of combinations of time delays upon a given state.

5. Stabilisation of nonlinear travelling wave from turbulence

Having described the ability of MTDF to successfully damp instability of nonlinear travelling waves, in this section we will demonstrate that such stabilised states can have suitably large basins of attractions and that it is possible to stabilise them from turbulence. We proceed as though the properties of the target solutions are unknown and allow the multiscale fluctuations of turbulence to fully develop before attempting stabilisation using MTDF. Turbulent fluctuations may include much longer time scales than the eigenperiods of the target solution. We imagine that some well-chosen delays can suppress sufficient spatio-temporal fluctuations in turbulence, without relaminarising the flow completely, at the same time as achieving stabilisation of the target solution. In other words, there is some motivation for including MTDF terms to widen the basin of attraction of the target state, as well as to force eigenvalues across the imaginary axis. However, this introduces a challenge: the more delays we use, the more parameters we need to adjust to stabilise nonlinear travelling waves successfully.

Without prior knowledge of the solution's instability, we seek to exploit some automatic techniques for obtaining stabilising parameter values. We have seen how $\tau _i$ may be chosen by analysing the data; in the next section we discuss the adaptive gain method of Lehnert et al. (2011), which automates the selection of $G$ and seeks to avoid a laborious parameter search.

5.1. Adaptive gain method

The speed-gradient method of Lehnert et al. (2011) seeks to find an optimal TDF gain by dynamically adjusting $G(t)$ by a ‘speed-gradient’ descent method. In order to exploit this method, we need to extend it to handle multiple terms and the translation operator, which is applied to the control term(s). We define the cost function, $Q_i(t)$, for the $i$th feedback term as

(5.1)\begin{equation} Q_i(t) = \tfrac{1}{2} \int | \mathcal{T}_z(s_z) \boldsymbol{u}(\boldsymbol{x},t-\tau_i)- \boldsymbol{u}(\boldsymbol{x},t)|^2 \,\mathrm{d}V , \end{equation}

where $\boldsymbol {x}$ is the position vector and successful control yields $Q_i(t) \rightarrow 0$ as $t \rightarrow \infty$. The speed-gradient algorithm in the differential form is given by

(5.2)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t}{G}_i(t) ={-} \gamma_i^G \boldsymbol{\nabla}_{G_i} \frac{\mathrm{d}{Q}_i}{\mathrm{d} t} , \end{equation}

where $\gamma _i^G>0$ is a free parameter controlling the descent rate and $\boldsymbol {\nabla }_{G_i}$ denotes $\partial /\partial G_i$. By taking time derivative of (5.1), we have

(5.3)\begin{equation} \frac{\mathrm{d}{Q}_i}{\mathrm{d} t} = \int g_i(\boldsymbol{x},t) \,\mathrm{d}V , \end{equation}

where

(5.4)\begin{equation} g_i(\boldsymbol{x},t) =[\mathcal{T}_z(s_z) \boldsymbol{u}(\boldsymbol{x},t-\tau_i)- \boldsymbol{u}(\boldsymbol{x},t)]\boldsymbol{\cdot} \left[\mathcal{T}_z(s_z)\frac{\partial \boldsymbol{u}}{\partial t}(\boldsymbol{x},t-\tau_i) -\frac{\partial \boldsymbol{u}}{\partial t}(\boldsymbol{x},t)\right], \end{equation}

and we have assumed that $\gamma _s$ is chosen such that ${\mathrm {d}s_z}/{\mathrm {d} t} \ll 1,$ hence $s_z$ is approximately constant for the purposes of updating $G_i(t).$ Using the Navier–Stokes momentum equations with MTDF in the form

(5.5)\begin{equation} \frac{\partial}{\partial t} \boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{f}(\boldsymbol{x},t) + \sum_i^N G_i [ \mathcal{T}_z(s_z) \boldsymbol{u}(\boldsymbol{x},t-\tau_i) - \boldsymbol{u}(\boldsymbol{x},t)] , \end{equation}

where $\boldsymbol {f}$ includes all the terms from the right-hand side of the momentum equation, alongside (5.2) and (5.3), we obtain the following formula:

(5.6)\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t}{G}_i(t) ={-} \gamma_i^G \int h_i(\boldsymbol{x},t) \,\mathrm{d} V, \end{equation}

where

(5.7)\begin{align} h_i(\boldsymbol{x},t)&= [ \mathcal{T}_z(s_z) \boldsymbol{u}(\boldsymbol{x},t-\tau_i)- \boldsymbol{u}(\boldsymbol{x},t)] \nonumber\\ &\quad \boldsymbol{\cdot} [\boldsymbol{u}(\boldsymbol{x},t) - 2 \mathcal{T}_z(s_z) \boldsymbol{u}(\boldsymbol{x},t-\tau_i)+ \mathcal{T}_z(2s_z) \boldsymbol{u}(\boldsymbol{x},t-2\tau_i)] . \end{align}

Note that, on taking $\boldsymbol {\nabla }_{G_i}$, only the MTDF terms contribute to $h_i.$ Equation (5.7) indicates that, when using this adaptive gain method, the code should store instantaneous velocity field data over 2 delay periods, $2 \tau _i,$ hence doubling the storage requirements. However, we implement a temporal interpolation procedure using cubic splines (Shaabani-Ardali et al. 2017) to enable us to store longer historical records of $\boldsymbol {u},$ albeit at the cost of some approximation error in the interpolation method.

5.2. Stabilising UB with MTDF and adaptive gain

Here, we attempt to stabilise UB at $Re=3000$ from a turbulent state, in the ($S,Z_2$)-symmetric subspace using MTDF. In this subspace, UB at $Re=3000$ has 4 pairs of complex unstable eigenvalues, and therefore it is unaffected by the odd-number limitation (Nakajima 1997). We approach this problem as if the eigenvalues are unknown, and we select $\tau _i$ values and other MTDF parameters without any linear analysis to steer our choices.

We consider first an MTDF case with two delays, $(\tau _1,\tau _2)=(2,9)$. We set the following parameters for $\tau _1$; $(t_s,G^{max},a,b,\gamma ^s)=(1000,0.5,0.1,100,0.1)$, i.e. TDF becomes active at $t=1000,$ giving a period of turbulent activity, and an initial sigmoid profile of $G(t)$. The adaptive gain method for $G_1(t)$, with $\gamma _1^G=0.1$, is then started at $t=2000$, and $G_1(t)$ evolves following the ODE (5.2). The second feedback term with gain, $G_2(t)$, evolves from zero at $t=2000$ through the adaptive gain method ($\gamma _2^G=0.1$). Without knowing successful values of $G_i$ and $s_z(t;\tau _i)$ in advance, we successfully stabilise UB at $Re=3000$ using this double-delay MTDF (see figure 9). In order to verify the non-invasive nature of the stabilisation in the MTDF cases, the definition of $R_{tot}$ requires updating

(5.8)\begin{equation} R_{tot} = \frac{\|\sum_i^N\left[\mathcal{T}_z(s_z)\boldsymbol{U}(r,\theta,z,t-\tau_i)-\boldsymbol{U}(r,\theta,z,t)\right]\|_2} {\|\boldsymbol{U}(r,\theta,z,t)\|_2}. \end{equation}

Successful, non-invasive stabilisation is quantitatively confirmed by the very small values of both $R_{tot}$ (figure 9a) and $I_{TDF}/I_{lam}$, the latter being of the order of $10^{-8}$ at $t=10\,000$ (see figure 9d). The phase speed, $c_z$, is also evolved onto its exact value through the adaptive translation method and we observe $G_i$ being adjusted by the speed-gradient method onto stabilising values. A snapshot of the stabilised UB at $t=20\,000$ is shown in figure 10(b) with a turbulent field at $t=0$ for comparison in figure 10(a).

Figure 9. Successful stabilisation of UB at $Re=3000$ from a turbulent state using MTDF (4.3) with two terms ($\tau _1=2, \tau _2=9$) and four terms ($\tau _1=2, \tau _2=9, \tau _3=17,\tau _4=33$). Panels show (a) $R_{tot}$, (b) $c_z$, (c, top) $G_i$ with two terms, (c, bottom) $G_i$ with four terms, (d) $I/I_{lam}$ and $I_{TDF}/I_{lam}$. The gain, $G_1(t)$, is switched on at $t=1000$. Here, $G_1(t)$ is initially increased following the sigmoid gain function (3.2) where $(t_s,G_1^{max},a,b)=(1000,0.1,0.1,100)$. At $t=2000$, the $G_i$ parameters begin to evolve using the adaptive gain method as described by (5.6), with a constant value of $\gamma _i^G=0.1$ applied to all delay terms. The adaptive shift method is switched on at $t=5000$, where $\gamma ^s=0.1$ and $c_z(0)=0.65$.

Figure 10. Two snapshots from the simulation using MTDF with two delays (see figure 9). (a) Turbulent field at $t=0$ and (b) stabilised UB at $t=20\,000$. The cyan isosurface denotes $u_z=-0.1$. The red and blue isosurfaces denote $\omega _z=0.15$ and $-0.15$, respectively, where $\omega _z = (\boldsymbol {\nabla }\times \boldsymbol {u})_z$ is the streamwise fluctuating vorticity.

Figure 9 also shows an MTDF case with four delays where $\tau _1=2$, $\tau _2=9$, $\tau _3=17$, $\tau _4=33$. Motivated by our frequency-domain analysis in § 3.3, these delays are carefully chosen to be non-commensurate and provide broad coverage of the temporal spectrum. As in the two-delay case, only $G_1(t)$ is switched on at $t=1000$ following the sigmoid gain function (3.2), where $(t_s,G_1^{max},a,b,\gamma ^s)=(1000,0.1,0.1,100,0.1)$. The adaptive gain method is switched on at $t=2000$ with $\gamma ^G=0.1$ for all the terms. Once the time-dependent gains are turned on, $I_{TDF}/I_{lam}$ starts to fluctuate with an amplitude smaller than 0.02 and then exhibits a damped oscillatory behaviour, similar to the case with two delays. There are no significant qualitative differences between the two- and four-term MTDF cases; in the four-delay case, $G_3$ is the second largest (recall $\tau _3=17$) and $G_2$ takes a smaller value at stabilisation than the two-delay case. Gain $G_4$ is the smallest but still makes a significant contribution. This suggests that the terms are all influencing the stabilisation, while the two-delay case demonstrates that not all terms are necessary for stabilisation. Figure 9(a) shows that the rate of attraction is of the same order in both cases, and saturation is reached at roughly the same time. This may suggest that the principal benefit of additional terms in MTDF is to widen the parameter windows under which UB is stabilised. If stabilising $G_i$ values are not obtained, it may be observed that the state is brought near the target solution but then moves away from it along some unstable manifold, similar to the blue curves in figure 8. The adaptive method gives the gains freedom to find values which succeed in stabilisation, instead of this close approach. It is worth noting that the energy injected by the forcing is only around 2 %–3 % of the energy injected by the imposed pressure gradient (see figure 9d). This means that, even at early times when the MTDF terms are largest, the overall energetic influence of MTDF on the flow is quite small. This may provide motivation for the development of TDF or MTDF in real-world experimental control situations; TDF does not need to intervene strongly to stabilise these travelling waves.

5.3. Travelling waves S2U, ML, N3, N4U, N5 and N7

Having shown that UB is able to be stabilised from the turbulent state by our adaptive MTDF approach, we now demonstrate the generality of these results by presenting the stabilisation of the other travelling waves outlined in table 1; the stabilising parameters can be found in table 2. First, we tackle S2U at $Re=2400$ in the $S$ symmetric subspace where the solution has one pair of unstable eigenvalues $\lambda _{\pm }=0.14\pm 0.13\mathrm {i}$ (note in the full space this solution violates the odd-number condition, see table 1). In theory it should be possible to stabilise this solution with a single TDF term, provided $\tau$ and $G$ take suitable values. However, as in the previous case, we continue with MTDF as though we did not know the stability information and start from the turbulent attractor. We use four terms with the same choices for $\tau _i,$ i.e. 2, 9, 17 and 33, as before. In the course of calibrating our adaptive methods, it was noted that starting $G_i$ from 0 could lead to some slight instability; individual ‘speed gradients’ begin with a large value (not only when targeting this solution but in general). A simple way to avoid this is to begin with small non-zero gains before starting to adapt, but again initialised with a sigmoid function (3.2). In this case we keep all other parameters the same as the four-term UB case of figure 9, only with $G_i^{max}=0.01$ for $i=2,3,4$ and in the $S$ subspace, see table 2. The stabilisation of S2U appears quite similar to that of UB. Figure 11(a) shows rapid stabilisation with the residual becoming very small $O(10^{-11}).$ Figure 11(b) shows the normalised energy input rate, which demonstrates the large-amplitude turbulent fluctuations before MTDF is activated and the settling into the final constant value at late times. It is observed that the residual only falls to small values once $c_z$ is dynamically adjusted to its exact value, shown in figure 11(c). In this example, the gains, plotted in figure 11(d), do not undergo any significant dynamical adjustment. Figure 12(a) shows snapshots of this stabilisation in the $(r,\theta )$ plane. The travelling wave ML is quite weakly unstable, in the $(S,Z_2)$ subspace at $Re=3000$ and the solution has only one unstable direction with $\lambda _{\pm }=0.0087\pm 0.019\mathrm {i}$ (see table 1). This provides a useful test case to demonstrate that, in this subspace at this Reynolds number, multiple solutions can be stabilised by only adjusting the MTDF parameters without requiring special treatment of the initial condition. Note that the unstable eigenfrequency is very small for this solution, which necessitates using larger delay periods, and tests with similar parameters used to stabilise UB either restabilise UB or relaminarise. In this instance we choose $\tau _i\in {2,9,33,150}$ and initiate with $c_z(0) = 0.71,$ the only other difference from earlier cases is that $\gamma ^G=0.5$ for all terms (see table 2). As is shown in figure 11, ML, despite being much less unstable than the other cases we have studied, shows weaker stabilisation. In addition, we see that $G_4$ undergoes significant growth once the speed-gradient method is initiated at $t=2000,$ becoming the largest gain of the four. Only once $G_4$ has grown do we observe the solution stabilising, indicating that this term is dominant in ensuring stabilisation. This fact is consistent with the frequency-domain analysis. In retrospect, a larger starting $G^{max}_4$ is likely to improve the stabilisation of ML. Nevertheless the adaptive approach has been able to determine stabilising gains automatically. Figure 12(b) shows snapshots of this stabilisation.

Table 2. Table summarising the parameters used in the four-term MTDF stabilisation of the travelling waves discussed. In all cases the start time $t_s=1000$ is used, $a=0.1,\,b=100$ for the sigmoid initialisation of $G_i$ with adaptivity started at $t=2000$ for all terms in all cases. Adaptivity of the translation (or phase speed) is initiated at $t=5000$, with the initial value shown in the table as $c_z(0)$.

Figure 11. Plots showing successful stabilisation of S2U, ML, N4U, N5, N7 in their respective subspaces, from a turbulent state using four-term MTDF (4.3). Specific parameter values can be found in table 2. (a) Shows the residual $R_{tot},$ (b) the energy input $I/I_{lam},$ (c) the phase speeds $c_z$ dynamically converging to their exact values and (d) the time-dependent gains $G_i(t).$

Figure 12. Snapshots of the $(r,\theta )$ plane at $z=0$ in the frame of reference translating with $c_z$ of the travelling waves for the MTDF cases outlined in table 2. From top to bottom (a) S2U, (b) ML, (c) N4U, (d) N5 and (e) N7. Coloured contours represent $u_z$ (each with 7 contours evenly spaced in $[-0.45,0.45]$ except S2U, which uses $[-0.35,0.35]$) and arrows represent the in-plane velocity vector. The final snapshot at $t=10\,000$ shows the stabilised state. Supplementary movies for these examples are available at https://doi.org/10.1017/jfm.2024.1188.

Noting that ML and UB have been successfully stabilised at the same Reynolds number, in the same symmetric subspace and with quite similar MTDF parameters, we have verified that taking the parameters used to stabilise ML (row 2 of table 2) and changing only $c_z(0)=0.65$ results in the stabilisation of UB. In other words both solutions can be obtained varying only $c_z(0).$

In $(S,Z_3)$ at $Re=2500$ and $\alpha =2.5$, a travelling wave, N3, is stable, meaning that stabilisation is not necessary. However, we have confirmed that, by applying the symmetry operator $SZ_3$ to the TDF terms, e.g.

(5.9)\begin{equation} \boldsymbol{F}_{MTDF}(r,\theta,z,t) = \sum_i^N G_i(t)[SZ_3\mathcal{T}_{z}(s_z) \boldsymbol{u}(r,\theta,z,t-\tau_i)-\boldsymbol{u}(r,\theta,z,t)] , \end{equation}

stabilisation can be obtained in the full space with some arbitrary (small) $\tau$ and $G_i,$ from a turbulent initial history. As explained in Lucas & Yasuda (2022), TDF or MTDF will simply drive the dynamics onto the symmetric subspace where the travelling wave is an attractor. This result indicates that the MTDF results shown earlier could be repeated in less-restrictive subspaces, with symmetry operators embedded in the MTDF terms, similarly to Lucas & Yasuda (2022). The travelling wave N4U is stabilised in the $(S,Z_4)$ subspace at $Re=2500$, $\alpha =1.7$ with the same MTDF parameters used in the S2U case, only now with $c_z(0)=0.5,$ as outlined in table 2. There are few significant differences in this case, $G_4$ and $G_1$ both grow to values ${\approx }0.3$ with $G_3$ also becoming relatively large. As with earlier examples, an analysis of the unstable eigenfrequencies would indicate which terms are necessary in this case, but in the interests of brevity we will not repeat a similar calculation here, noting that, even without this analysis, significant parameter tuning was not necessary to stabilise this solution. As with the other travelling waves, the residual $R_{tot}$ only begins to fall to values for which MTDF may be considered non-invasive as the phase speed converges to its final value, i.e. $t>5000$. Figure 12(c) shows snapshots of this stabilisation.

In an effort to provide some evidence for the assertion that this method will enable new, unknown, solutions to be obtained, we have conducted investigations in the $(S,Z_5)$ and $(S,Z_7)$ subspaces. While solutions have been reported in the $(S,Z_5)$ subspace before (Pringle et al. 2009), insufficient stability or phase speed data are tabulated or available in databases to enable stabilising parameters to be predetermined (the previous states discussed in this paper are available at Openpipeflow.org and/or are tabulated in Willis et al. 2013). In the $(S,Z_7)$ subspace, we are unaware of any solutions having been reported. Therefore, for these cases, we do not, or cannot, target any particular solutions, and so this is a strong test of the ability of the method to obtain at least ‘unknown’ solutions, if not ones which are new to the field. The only observations that we use to inform our attempt are that, as the rotational symmetry order parameter, $m_p$, increases, the minimal axial wavenumber $\alpha$ should also increase for such travelling waves. In addition, we might expect a reduction in the phase speed. Therefore for $(S,Z_5)$ we use $\alpha =2$ and the initial condition will be $c_z(0)=0.47.$ We leave all other parameters the same as the successful N4U case, e.g. $Re=2500$ and the same $G_i$ and descent parameters. This results in a successful stabilisation of a travelling wave, which we denote as N5. In this first attempt, we observe a very slow rate of stabilisation, which appears to be a result of oscillation in the $c_z(t)$ descent. Reducing $\gamma ^s$ to 0.01 prevents this oscillation and results in a more acceptable stabilisation rate; figure 11(a) shows the residual falling to $10^{-9}$ by time 20 000. Note that we choose this reduction in $\gamma ^s$ only so that the accompanying plots in figure 11 are comparable to the other cases and are more easily interpreted (i.e. we avoid computing and displaying very late times ${\sim }50\,000$). Figure 12(d) shows snapshots of this stabilisation.

In $(S,Z_7)$, we study $Re=3500$ and $\alpha =3$, where sustained turbulence is observed prior to starting control. We use the same MTDF parameters as for N5 only our first attempt used $c_z(0)=0.45.$ This case was unsuccessful, and it was clear from the behaviour of $c_z(t)$ that this initial condition was too large. Choosing $c_z(0)=0.4$ resulted in the successful case displayed in figure 11, labelled N7. This was the extent of the trial-and-error required to stabilise this new travelling wave. Figure 11(a) shows quite fast stabilisation, and the energy input rate settles to its final value quite early (see figure 11b). This may be, in part, due to 0.4 being a very good guess for the phase speed, as indicated by the evolution of $c_z(t)$ in figure 11(c). Figure 12(e) shows snapshots of this stabilisation.

In the case of the N5 and N7 solutions, since they are unknown, we converge them using the Newton–GMRES–hookstep method and obtain their leading eigenvalues by Arnoldi iteration, as shown in table 1. Taking the final snapshot as the starting guess and a translation of $s_z=\tau _1 c_z(T)$ from the MTDF result gives a starting Newton residual around $10^{-6}$ or $10^{-7}$ for these cases, with Newton converging in one step.

6. Conclusion

In this paper, we present the first successful non-invasive stabilisation of highly unstable nonlinear travelling waves (having multiple unstable directions) in a straight cylindrical pipe through the use of a control method involving TDF. Our novel TDF protocol allows for the stabilisation of multiple nonlinear travelling waves at a range of Reynolds numbers in a variety of symmetric subspaces from a generic turbulent history and with speculative control parameters.

Furthermore, our development of the TDF method has led to a deeper understanding of the principles governing TDF through an ‘approximate’ linear stability analysis and frequency-domain analysis (see e.g. § 3.3). We have shown that the effect of TDF and MTDF on the unstable part of solution's eigenvalue spectrum can be approximated surprisingly well, enough to point a parameter study in the right direction. Moreover, it has provided clearer insight into the frequency-domain interpretation of the control method, which in turn gives a helpful means to choose delay periods in these cases. Finally, we have shown that, if stabilisation is nearly achieved, it can be possible to diagnose eigenfrequencies and hence pick more appropriate time delays, without the need for an a priori stability analysis of the target solution.

In order to enhance the performance of our TDF method, we implemented several optimisation methods that enable the feedback term(s) to vanish, as depicted in figure 9(b). Because the bulk flow in pipe flow is driven by an imposed axial pressure gradient, all invariant solutions take the form of relative solutions, such as travelling waves and relative periodic orbits. Therefore, applying a translation operator to the delay term(s) is essential. To achieve this, we utilised the adaptive shift method (see § 3.1), which dynamically adjusts the translations to match the phase speeds of the target solution through a simple ODE (3.3). We have demonstrated here, for the first time, that, by changing the initial condition for $s_z$ or $c_z,$ different travelling waves can be stabilised.

In the chaotic pipe flow system, we find that MTDF control is effective at improving the control's ability to stabilise a wide range of unstable eigenfrequencies, as shown in § 4.2. A helpful consequence of MTDF is that it will serve to damp very slow temporal oscillations, which are typical in pipe flow (Shih, Hsieh & Goldenfeld 2016). As demonstrated in the results of §§ 5.2 and 5.3, successful stabilisations are found when MTDF is initiated with a short time delay ($\tau _1=2$ in all our cases), which acts for some time period before the rest of the terms become active. The effect is to suppress slow oscillations, importantly without relaminarising, giving the remaining terms a better chance at controlling the target travelling wave. In other words, this initial delay is an effective means to widen the basin of attraction. In some cases $\tau _1$ contributes directly to altering the linear stability of the travelling waves; in others it only provides a nonlinear effect.

We have also sought to avoid expensive parametric studies. To achieve this, we have introduced the adaptive gain method (§ 5.1) into our TDF protocol. This method automates our search for an appropriate gain, thereby avoiding an exhaustive parametric search for $G_i$ values. We have observed this approach to be highly effective; for instance in § 5.3 when stabilising ML, the longest time delay gain, $G_4$, is observed to grow significantly, compared with the other terms, showing that $\tau =150$ was important in ensuring stabilisation. This is consistent with our frequency analysis when noticing that this travelling wave has a very large unstable eigenperiod of around 330.

During this work we have demonstrated that TDF, or more accurately MTDF, can stabilise multiple states at the same parameter values. In other words, multiple attractors can coexist; UB and ML are stabilised by varying only the initial $c_z.$ In this example, the two states are relatively well separated in phase space (upper and lower branch solutions, see figure 1) so it is perhaps surprising that ML is stabilised from a turbulent initial condition. This is a promising result as it demonstrates that the method does not require significant intervention to move from one solution to the next. However, it does open up a number of interesting questions. In particular, how might one design a systematic search, or data-driven approach, to explore basins of attraction of various potential solutions? For instance, UB at $Re=2400$ in the $S$ subspace is highly unstable and has 9 complex pairs of eigenvalues, in theory this could be stabilised by MTDF, however, it will ‘compete’ with S2U, which is much less unstable (1 unstable direction), is also upper branch and has a very similar phase speed (see table 1). If any of the 9 eigenvalues has a particularly small eigenfrequency, necessitating a long time delay, then close proximity to the solution is likely to be necessary for any stabilisation to be successful. We have also seen that trialling various subspaces and/or embedding symmetry operators into TDF terms is a useful way to obtain different solutions, avoid the odd-number issue and avoid dealing with multiple attractors. However, we have not conducted an exhaustive search of all possible subspaces, pipe lengths and $Re.$ We expect that travelling-wave solutions, or even relative periodic orbits, can be stabilised with TDF, at a wide range of Reynolds numbers, pipe lengths and subspaces – not to mention in other wall-bounded shear flows or systems with additional physics such as stratification or rotation where new dynamically important temporal frequencies arise.

The experimental applicability of these results is not immediately obvious. In order to address this, at least in a proof-of-concept manner, we have trialled a handful of cases where the control force $\boldsymbol {F}_{TDF}$ is premultiplied by the Kronecker delta function $\delta _{r_i=r_{N-1}}$. In other words, TDF is no longer acting as a full-state control method but is applied only to the cylinder of grid points adjacent to the wall. For instance, the stabilisation of UB shown in figure 3 ($\tau =5, G=0.1$) is successful with this spatially localised control without any further parameter tuning (figures not shown for brevity). The fact that TDF can be successful when only applied to a subset of the degrees of freedom in the problem gives some optimism that an experimental version of TDF may be successful. It is also worth noting that TDF has been successful in an experimental Taylor–Couette flow (Lüthje et al. 2001) by applying a boundary forcing.

In this paper, we have not tackled the odd-number limitation (Nakajima 1997), which is a contentious issue in the TDF literature. Further improvements to this method are necessary in order to stabilise odd-number solutions. One possibility is the half-period TDF (Nakajima & Ueda 1998a), which is similar in spirit to the use of symmetries here and in Lucas & Yasuda (2022). Another is ‘act-and-wait’ TDF (Pyragas & Pyragas 2018, 2019), where a time-dependent switching of the TDF gain is applied, meaning uncontrolled dynamics is always used in the delay period, or ‘unstable ETDF’ (Pyragas 2001), where an additional unstable degree of freedom is introduced into the problem to create an even pair of exponents. We have demonstrated that there is some potential for this method in controlling nonlinear states in spatio-temporal chaos, which will hopefully serve as motivation for further developments tackling both the odd-number issue and even higher-dimensional problems at large Reynolds numbers and in large domains. Potentially, the most promising avenue for TDF in fluid flows is as a control method in a real physical system where a non-trivial flow state, perhaps of a specific dissipation or mixing rate, is desired but full spatio-temporal chaos is not. We have shown that, with careful application of a time-delayed approach, various kinds of target solutions can be relatively easily obtained with minimal intervention.