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 Reference Salwen, Cotton and Grosch1980; Meseguer & Trefethen Reference Meseguer and Trefethen2003).

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 Reference Darbyshire and Mullin1995; Duguet et al. Reference Duguet, Pringle and Kerswell2008a; Willis et al. Reference Willis, Cvitanović and Avila2013) or a constant pressure gradient (Wedin & Kerswell Reference Wedin and Kerswell2004; Shimizu & Kida Reference Shimizu and Kida2009). In what follows we will seek comparisons with ECSs enumerated in Willis et al. (Reference Willis, Cvitanović and Avila2013); 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. Reference Marensi, Ding, Willis and Kerswell2020). 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. Reference Marensi, Ding, Willis and Kerswell2020). 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 Reference Willis2017) 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 Reference Willis2017). 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 Reference Willis2017). 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. (Reference Willis, Cvitanović and Avila2013), 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 Reference Wedin and Kerswell2004; Pringle et al. Reference Pringle, Duguet and Kerswell2009; Willis et al. Reference Willis, Cvitanović and Avila2013). 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.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 (Reference Lucas and Yasuda2022). 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 (Reference Lucas and Yasuda2022). 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 Reference Meseguer and Trefethen2003), we attempt to stabilise known unstable travelling waves (Willis et al. Reference Willis, Cvitanović and Avila2013) 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 Reference Lucas and Yasuda2022, 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 Reference Viswanath2007, Reference Viswanath2009) packaged with the openpipeflow code (Willis Reference Willis2017). 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 Reference Nakajima1997).

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 Reference Lucas and Yasuda2022). 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. Reference Akervik, Brandt, Henningson, Hoepffner, Marxen and Schlatter2006; Shaabani-Ardali et al. Reference Shaabani-Ardali, Sipp and Lesshafft2017) 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 (Reference Lucas and Yasuda2022) 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. Reference Linkmann, Knierim, Zammert and Eckhardt2020, 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. Reference Akervik, Brandt, Henningson, Hoepffner, Marxen and Schlatter2006; Shaabani-Ardali et al. Reference Shaabani-Ardali, Sipp and Lesshafft2017), 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.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 (Reference Socolar, Sukow and Gauthier1994). 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 (Reference Ahlborn and Parlitz2004), 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. Reference Socolar, Sukow and Gauthier1994; Sukow et al. Reference Sukow, Bleich, Gauthier and Socolar1997).

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.1. Adaptive gain method

The speed-gradient method of Lehnert et al. (Reference Lehnert, Hövel, Flunkert, Guzenko, Fradkov and Schöll2011) 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. Reference Shaabani-Ardali, Sipp and Lesshafft2017) 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 Reference Nakajima1997). 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 (Reference Lucas and Yasuda2022), 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 (Reference Lucas and Yasuda2022). 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. Reference Pringle, Duguet and Kerswell2009), 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. Reference Willis, Cvitanović and Avila2013). 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.