On the onset of long-wavelength three-dimensional instability in the cylinder wake

Abstract

We study the onset of the three-dimensional mode A instability in the near wake behind a circular cylinder. We show that long-wavelength perturbations organise in a time-shifting pattern such that the in-plane velocity in each streamwise slice corresponds to the base flow solution at shifted times. This observation introduces an additional unifying characteristic for certain mode A type instabilities. We then analyse the mechanisms which control the growth or decay of these perturbations and highlight the crucial role played by the tilting mechanism which operates via non-local interactions in a manner similar to Biot–Savart induction. We characterise its domain of influence using a Green's function-based approach which allows us to rationalise the non-trivial dependence of the growth rate on the spanwise wavenumber. We connect this behaviour to the subtle balance between the local growth of the perturbations as they are swept along by the flow and the feedback on the perturbations that are generated during the next period of the time-periodic base flow. Finally, we discuss generalisations of our findings to other types of flows.

1. Introduction

The flow of an incompressible viscous fluid around an infinitely long circular cylinder is characterised by the Reynolds number, ${\textit {Re}}=U_\infty d/\nu$ (defined by the free stream velocity $U_\infty$, the cylinder diameter $d$ and the kinematic viscosity $\nu$). With an increase of ${\textit {Re}}$, the flow undergoes several stages of stability loss before it becomes turbulent (Williamson 1996c). A key feature of the flow is the von Kármán vortex street which appears soon after the primary instability of the flow at ${\textit {Re}}={\textit {Re}}_{0}$ when the two-dimensional steady flow becomes time-periodic via supercritical Hopf bifurcation. Critical Reynolds numbers observed in experiments and obtained using theoretical analysis agree, ${\textit {Re}}_{0}\approx 46\unicode{x2013}47$ (Mathis, Provansal & Boyer 1984; Jackson 1987; Dušek, Le Gal & Fraunié 1994). This instability does not immediately lead to the appearance of the von Kármán vortex street – the formation of the vortices happens at a slightly larger Reynolds number far in the wake (approximately 100 diameters downstream) (Heil et al. 2017). As the Reynolds number is increased further, the two-dimensional time-periodic flow becomes unstable to two distinct modes (A and B) of three-dimensional instability (Barkley & Henderson 1996), which are also observed experimentally (Williamson 1988). The modes have different spatio-temporal structure and length scale (approximately four and one diameter of the cylinder, respectively).

The mode A instability arises at a critical Reynolds number of ${\textit {Re}}_A=188.5\pm 1$ and a wavelength of $\lambda _A=3.96\pm 0.02$ (Barkley & Henderson 1996) via a subcritical bifurcation (Henderson & Barkley 1996; Behara & Mittal 2010; Akbar, Bouchet & Dušek 2011). These theoretical predictions agree with experiments, see discussions by Miller & Williamson (1994), Williamson (1996a), Akbar et al. (2011) and Jiang et al. (2016b). Barkley (2005) demonstrated that the instability originates in the vortex formation region by applying a Floquet stability analysis to the various flow subregions. This showed that at ${\textit {Re}}=190$, the confined flow in the vortex formation region ($0\le x\le 3$ and $|y|\le 1.5$) still exhibited a mode A instability (as manifested by the same dependence of the Floquet multiplier on the spanwise wavenumber as for the entire computational domain), whilst the developed wake (region $2.25\le x\le 25$ and $|y|\le 4$) turned out to be stable. This finding is supported by Giannetti, Camarri & Luchini (2010), who performed a sensitivity analysis of the dominant Floquet modes to localised structural perturbation and also provided time-resolved details of the most sensitive subregions of the flow.

A distinctive characteristic of mode A behind a circular cylinder is its degeneration into the neutral two-dimensional mode in the limit of infinite spanwise wavelength, as highlighted by Barkley & Henderson (1996). It is well known that, in general, periodic solutions $\boldsymbol {U}(\boldsymbol {x},t)$ of autonomous problems admit neutrally stable Floquet modes in the form of $\partial \boldsymbol {U}(\boldsymbol {x},t)/\partial t$ (Iooss & Joseph 1990, § VII.6.2). Therefore, given that mode A shares its symmetry with this two-dimensional neutral mode, it inherits the symmetry of the base flow $\boldsymbol {U}(\boldsymbol {x},t)$. Three-dimensional instabilities linked to such neutral modes also occur in other problems, and we show that this general mathematical fact has implications for the kinematics of long-wavelength three-dimensional instabilities, thus elucidating a perturbation pattern for a class of instabilities.

Mode A type instabilities are also observed in other flows, e.g. behind elongated bluff cylinders, oscillating cylinders, rotating cylinders, square and elliptic cylinders, airfoils, and behind cylinders moving near a wall (Ryan, Thompson & Hourigan 2005; Leontini, Thompson & Hourigan 2007; Luo, Tong & Khoo 2007; Sheard, Fitzgerald & Ryan 2009; Lo Jacono et al. 2010; Leontini, Lo Jacono & Thompson 2015; Rao et al. 2015; Agbaglah & Mavriplis 2017; He et al. 2017; Rao et al. 2017; Thompson, Leweke & Hourigan 2021). However, it is interesting to note that there is no universally accepted definition that allows one to classify a particular three-dimensional instability as being mode A. One possible way to do this is by constructing a continuous transformation between different problems and tracking the relevant solution branch; see, e.g. Leontini et al. (2015). A less rigorous but common approach is to compare what are thought to be ‘intrinsic’ attributes of the mode A pattern, such as its critical wavelength, the local distribution of the perturbations and the spatio-temporal symmetry of the perturbations. Yet, mode A type perturbations can emerge on the background of non-symmetric base flow, and their spanwise wavelength can be of the order of tens of diameters of the cylinder; see, e.g. the flow around an elliptic and rotating cylinder (Rao et al. 2015,  2017).

Over the years, many attempts have been made to explain the physical mechanism responsible for the onset of the mode A instability, e.g. by analysing simplified flows that have certain key features observed in the actual, usually much more complicated, flow with the aim of predicting the pattern and critical parameters of the instability. The best known attempt of this type exploits the similarity of the perturbed base flow vortices with the structures that appear in the course of an elliptic instability of a stationary two-dimensional flow with elliptic streamlines (Lagnado, Phan-Thien & Leal 1984; Landman & Saffman 1987; Waleffe 1990; Kerswell 2002). This similarity was first noted by Williamson (1996b), who hypothesised that the mode A instability arises via the elliptic instability of the developing vortices in the vortex formation region. The hypothesis was supported by Leweke & Williamson (1998b) and Thompson, Leweke & Williamson (2001). In the latter work, the hypothesis was extended to a cooperative elliptic instability of two counter-rotating forming vortices (shedding from both sides of the cylinder) based on the resemblance with data by Leweke & Williamson (1998a) on three-dimensional instability of a vortex pair. The analysis provides an estimate for the spanwise wavelength of the mode A instability (of approximately three diameters of the cylinder) which agrees well with experimental observations. Ryan et al. (2005) and Leontini et al. (2007) found other correlations with the elliptic instability hypothesis for flows around other bluff bodies.

However, the hypothesis does not take into account the self-excited nature of the instability, i.e. the fact that the three-dimensional perturbations created in the forming vortex during a certain phase of the time-periodic base flow not only undergo local growth (while being advected by the flow), but also provide positive or negative feedback on the development of the instability during the next period. It is this balance between local growth and feedback that is at the heart of the instability mechanism – within the framework of Floquet analysis, it is characterised by the value of the Floquet multiplier. Furthermore, the flow in the forming vortex core is non-stationary, non-uniform and interacts with perturbations in other parts of the flow, and is, therefore, significantly more complex than assumed in the simplified models. This means that the role of the intensive growth of perturbations outside the vortex core is still not clear. Indeed, it is known that the growth of perturbations has two distinct phases that occur when the perturbations grow predominantly in the forming vortex and in the braid shear layer (Williamson 1996b; Leweke & Williamson 1998b; Thompson et al. 2001; Aleksyuk & Shkadov 2018,  2019). The elliptic instability hypothesis assumes that the amplification of perturbations during the second phase only has a secondary effect on the instability. Some support for this interpretation is provided by Thompson et al. (2001) and Julien, Ortiz & Chomaz (2004).

An alternative view on the local mechanisms for the instability was proposed by Giannetti et al. (2010) and Giannetti (2015), which takes into account the self-excited nature of the instability. Giannetti (2015) performed a stability analysis, based on applying the Lifschitz–Hameiri theory (Lifschitz & Hameiri 1991) in the limits ${\textit {Re}}\rightarrow \infty$ and $\gamma \rightarrow \infty$, along the closed periodic orbits found in the vortex formation region. They demonstrated that the local evolution of perturbations along a specific orbit could reproduce the instability characteristics of modes A and B. However, the quantitative agreement of the predictions with experimental observations is poor, presumably because of the strong assumptions on ${\textit {Re}}$ and $\gamma$. Indeed, Jethani et al. (2018) carried out a similar analysis that included finite ${\textit {Re}}$ and $\gamma$ corrections and obtained better agreement with the critical parameters for mode B and suggested that the mode B instability could be a manifestation of the local instability on the closed orbits. To our knowledge, there is still no quantitative agreement on mode A.

For a discussion of other, earlier hypotheses regarding the development of the mode A instability, based on the Benjamin–Feir instability (Leweke & Provansal 1995) or the centrifugal instability (Brede, Eckelmann & Rockwell 1996), say, we refer to Leweke & Williamson (1998b) and Thompson et al. (2001).

The aim of this paper is to clarify the mechanisms for the onset of mode A instability, specifically, the paper addresses two questions.

1. (i) What is the explanation for the pattern of mode A at the early (linear) stage of its development (§ 5)?

2. (ii) What physical mechanisms define whether this pattern is unstable at a specific Reynolds number and spanwise wavelength (§ 6)?

The structure of the paper is as follows. In §§ 2–4, we describe the problem formulation, the two-dimensional time-periodic base flow and the three-dimensional linear stability analysis performed to obtain the dominant Floquet modes. In § 5, we answer question (i) by considering a simplified case of small spanwise wavenumbers. Then, in § 6, we address question (ii) by describing perturbations in terms of perturbations to the in-plane vorticity. The results are summarised in § 7. Appendix A provides details of the numerical simulations. In Appendices B and C, we discuss the action of the basic physical mechanisms for the change of the in-plane vorticity of a fluid particle and the derivation of the Green's function for the screened Poisson equation to describe non-local interactions of perturbations.

2. Problem formulation

The flow of an incompressible viscous fluid around an infinitely long circular cylinder is described in the Cartesian coordinate system $\boldsymbol {x}=(x, y, z)$ with the $z$-axis coinciding with the axis of the cylinder and the $x$-axis aligned with the incoming flow velocity. All quantities are considered in non-dimensional form based on the diameter of the cylinder $d$, the free stream velocity $U_\infty$ and the fluid density $\rho _\infty$:

(2.1a–d)\begin{equation} t=\frac{U_\infty \tilde{t}}{d},\quad \boldsymbol{x}=\frac{\tilde{\boldsymbol{x}}}{d},\quad p=\frac{\tilde{p}}{\rho_\infty U_\infty^2},\quad \boldsymbol{u}=\frac{\tilde{\boldsymbol{u}}}{U_\infty}. \end{equation}

Here, $t$, $p(\boldsymbol {x}, t)$ and $\boldsymbol {u}(\boldsymbol {x}, t)=(u, v, w)$ are time, pressure and the velocity vector; a tilde is used to distinguish dimensional variables from their non-dimensional equivalents.

The solution $p$, $\boldsymbol {u}$ depends on only one parameter – the Reynolds number ${\textit {Re}}=U_\infty d/\nu$ (where $\nu$ is the coefficient of kinematic viscosity), and satisfies the Navier–Stokes equations

(2.2a,b) \begin{cases}{} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \\ \frac{\partial{\boldsymbol{u}}}{\partial t}+\boldsymbol{N}(\boldsymbol{u},\boldsymbol{u})={-}\boldsymbol{\nabla} p+\frac{1}{{\textit{Re}}}\nabla^2\boldsymbol{u} \end{cases}

subject to no-slip boundary condition $\boldsymbol {u}=(0,0,0)$ at the surface of the cylinder and $\boldsymbol {u}\rightarrow (1,0,0)$ as $\boldsymbol {r}=(x,y)\rightarrow \infty$. Here, the nonlinear advection term is expressed using $\boldsymbol {N}(\boldsymbol {u},\boldsymbol {v})=[(\boldsymbol {u} \boldsymbol {\cdot }\boldsymbol {\nabla }) \boldsymbol {v}+(\boldsymbol {v} \boldsymbol {\cdot }\boldsymbol {\nabla }) \boldsymbol {u}]/2$. We use the arguments $\boldsymbol {r}=(x,y)$ and $\boldsymbol {x}=(x,y,z)$ to indicate a function's dependence on the in-plane and full three-dimensional coordinates, respectively.

3. Two-dimensional base flow

The base flow velocity vector $\boldsymbol {U}(\boldsymbol {r},t)=(U,V,0)$ and pressure $P(\boldsymbol {r},t)$ satisfy (2.2), which we solved numerically using a second-order stabilised finite element method on triangular meshes with a second-order discretisation in time (see Appendix A).

In the range of the Reynolds number we consider in this paper ($50\le {\textit {Re}}\le 220$), the base flow in the near wake is $T$-periodic in time, e.g. $\boldsymbol {U}(\boldsymbol {r},t+T)=\boldsymbol {U}(\boldsymbol {r},t)$, and possesses the following symmetry:

(3.1)\begin{equation} \begin{pmatrix} U\\ V\\ P \end{pmatrix}(x,y,t+T/2)=\begin{pmatrix} U\\ -V\\ P \end{pmatrix}(x,-y,t). \end{equation}

As an example, figure 1 illustrates the base flow solution at ${\textit {Re}}=220$. Figure 1(a) shows the contours of the vorticity, $\varOmega =\partial V/\partial x-\partial U/\partial y$, and highlights where vortices are created and where they reach their fully formed state; the contours in figure 1(b) show the positive eigenvalue $S$ of the strain rate tensor and the associated principal directions – the latter indicate the direction of maximum stretching in the flow; finally, figure 1(c) shows the ratio $\kappa = 2S/|\varOmega |$, where $\varOmega /2$ is the local rate of rotation. Thus, $\kappa$ is a measure of the relative importance of stretching and rotation, and the lines $\kappa =1$ define the boundaries between hyperbolic (stretching-dominated) and elliptic (rotation-dominated) regions.

Figure 1. Plots of the base flow at ${\textit {Re}}=220$ in terms of (a) the vorticity $\varOmega$; (b) the positive eigenvalue $S$ of the strain rate tensor and its principal direction $\varPhi$ (shown by red line segments); (c) the ratio $\kappa =2S/|\varOmega |$ on a logarithmic scale. Solid lines correspond to the boundaries between elliptic and hyperbolic regions, $\kappa =1$. The time $t=0$ corresponds to the maximum of the lift coefficient. Panel (a) also identifies the key flow regions: the forming vortex, the braid shear layer and the fully formed vortex.

4. Dominant Floquet modes of three-dimensional perturbations

To elucidate the mechanisms responsible for the onset of the three-dimensional instability, we consider the initial stages of its development when the deviation from the two-dimensional time-periodic base flow is small. The perturbation velocity vector $\boldsymbol {u}'(\boldsymbol {x},t)=(u',v',w')$ and pressure $p'(\boldsymbol {x},t)$ satisfy the linearised Navier–Stokes equations,

(4.1a,b) \begin{cases}{} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}'=0, \\ \frac{\partial{\boldsymbol{u}'}}{\partial t}+2\boldsymbol{N}(\boldsymbol{U},\boldsymbol{u}')={-}\boldsymbol{\nabla} p'+\frac{1}{{\textit{Re}}}\nabla^2\boldsymbol{u}', \end{cases}

and homogeneous boundary conditions $\boldsymbol {u}'=(0,0,0)$ at the surface of the cylinder and as $\boldsymbol {r}\rightarrow \infty$. We seek perturbations with spanwise wavenumber $\gamma$:

(4.2)\begin{equation} \begin{pmatrix} u'\\ v'\\ w'\\ p' \end{pmatrix}(\boldsymbol{x},t)=\begin{pmatrix} \hat{u}\\ \hat{v}\\ \mathrm{i}\hat{w}\\ \hat{p} \end{pmatrix}(\boldsymbol{r},t)\exp({\mathrm{i}\gamma z})+\begin{pmatrix} \hat{u}^*\\ \hat{v}^*\\ -\mathrm{i}\hat{w}^*\\ \hat{p}^* \end{pmatrix}(\boldsymbol{r},t)\exp({-\mathrm{i}\gamma z}), \end{equation}

which satisfy

(4.3a,b) \begin{cases}{} \hat{\boldsymbol{\nabla}}^*\boldsymbol{\cdot}\hat{\boldsymbol{u}} = 0,\\ \frac{\partial{\hat{\boldsymbol{u}}}}{\partial t}+2\boldsymbol{N}(\boldsymbol{U},\hat{\boldsymbol{u}}) ={-}\hat{\boldsymbol{\nabla}}\hat{p}+ \frac{1}{{\textit{Re}}}\left(\nabla^2\hat{\boldsymbol{u}}-\gamma^2\hat{\boldsymbol{u}}\right), \end{cases}

where $\hat {\boldsymbol {u}}(\boldsymbol {r},t)=(\hat {u}, \hat {v}, \hat {w})$, $\hat {\boldsymbol {\nabla }}=(\partial /\partial x,\partial /\partial y,\gamma )$ and $\hat {\boldsymbol {\nabla }}^*=(\partial /\partial x,\partial /\partial y,-\gamma )$; $\hat {u}^*$, $\hat {v}^*$, $\hat {w}^*$ and $\hat {p}^*$ are complex conjugates of $\hat {u}$, $\hat {v}$, $\hat {w}$ and $\hat {p}$.

The system of equations (4.3) is linear and has $T$-periodic coefficients. Therefore, we represent each function with a hat, say $\hat {u}(\boldsymbol {r},t)$, as $\exp (\sigma t)u_{p}(\boldsymbol {r},t)$, where $u_{p}(\boldsymbol {r},t)$ is a $T$-periodic function and $\sigma =\sigma _{r}+\mathrm {i}\sigma _{i}$ is a complex number. These modes are either real or come in conjugate pairs since the coefficients of the system are real. Perturbations at a given $\gamma$ correspond to the combination of waves travelling along the $z$-axis with speed $\pm \sigma _{i}/\gamma$, for instance,

(4.4)\begin{align} u'(\boldsymbol{x},t) &= a(\boldsymbol{r},t)\,{\rm e}^{\sigma_{r}t} \left[C_1 \exp\left({\mathrm{i}\left[\sigma_{i}t+\gamma z+\phi(\boldsymbol{r},t)\right]}\right) +C_1^* \exp\left({-\mathrm{i}\left[\sigma_{i}t+\gamma z+\phi(\boldsymbol{r},t)\right]}\right)\right. \nonumber\\ &\quad \left.+C_2^*\exp\left({\mathrm{i}\left[\sigma_{i}t-\gamma z+\phi(\boldsymbol{r},t)\right]}\right)+C_2 \exp\left({-\mathrm{i}\left[\sigma_{i}t-\gamma z+\phi(\boldsymbol{r},t)\right]}\right)\right], \end{align}

where the $T$-periodic part of the solution is expressed using the amplitude $a(\boldsymbol {r},t)$ and argument $\phi (\boldsymbol {r},t)$: $u_{p}(\boldsymbol {r},t)=a\exp (\mathrm {i}\phi )$; the constants $C_1$ and $C_2$ appear as coefficients in a linear combination of complex-conjugate solutions. When $\sigma$ is real, the solution degenerates into a standing wave.

If at least one Floquet multiplier $\mu =\exp (\sigma T)$ lies outside the unit circle ($|\mu |>1$), the flow is unstable. Given ${\textit {Re}}$ and $\gamma$, we seek only the dominant mode with the largest $|\mu |$ using the numerical method described in Appendix A. Figure 2 shows the dependence of the dominant Floquet multiplier on $\gamma$ at ${\textit {Re}}=220$. There are three intervals, which correspond to modes A ($0<\gamma \le 2.6$) and B ($\gamma \ge 8.5$) with real Floquet multipliers, and quasi-periodic modes with complex $\mu$ in the intermediate range of $\gamma$. The flow is unstable ($|\mu |>1$) to perturbations of mode A for $1.1<\gamma <2.1$ (hatched region).

Figure 2. Dominant Floquet multiplier at ${\textit {Re}}=220$ (obtained by two methods, see Appendix A.2) and comparison with the data of Barkley & Henderson (1996). The hatched yellow area highlights unstable perturbations.

Figure 3 illustrates the changes in the corresponding eigenfunctions with $\gamma$ by plotting the distribution of perturbation kinetic energy

(4.5)\begin{equation} e=\frac{1}{L}\int_{0}^{L}\frac{1}{2}\left(u'^2+v'^2+w'^2\right){\rm d} z=|\hat{\boldsymbol{u}}|^2, \end{equation}

and in-plane perturbation velocity vectors; here $L=2{\rm \pi} /\gamma$ is the wavelength. The pattern of the perturbations remains qualitatively similar over the range of $\gamma$ considered. An increase in $\gamma$ causes the perturbations inside the formed vortex ($x>3$) to shift outside of it, but there is little change to the perturbation pattern in the vortex formation region (highlighted by the yellow shaded region).

Figure 3. Pattern of mode A perturbations at ${\textit {Re}}=220$ and $0\le \gamma \le 2.2$: perturbation energy $e$ (greyscale colour contours) and in-plane ($z=0$) perturbation velocity (arrows). Solid lines are the base flow vorticity isolines $\varOmega =\pm 1$. All plots are snapshots at $t=0.5T$, corresponding to the minimum of the lift coefficient. Perturbations at $\gamma =0$ are obtained by time differentiation of the base flow solution, see (5.2). The yellow shaded regions show the vortex formation region. Note that the greyscale contour levels were adjusted manually to highlight the similarities and differences of the perturbation patterns.

5. The pattern of long-wavelength perturbations

Given that the overall features of the flow field, particularly in the vortex formation region, do not change qualitatively with variations in the wavenumber (see, e.g. figure 3 at ${\textit {Re}}=220$ and $0\le \gamma \le 2.2$), we analyse the pattern of the three-dimensional perturbations in the small-$\gamma$ (i.e. long-wavelength) regime.

We start with the case $\gamma =0$. Taking the time-derivative (denoted by an over-dot) of the Navier–Stokes equations and boundary conditions for the base flow $(\boldsymbol {U}, P)$ leads to

(5.1a,b) \begin{cases} \boldsymbol{\nabla}\boldsymbol{\cdot}\dot{\boldsymbol{U}}=0,\\ \displaystyle\frac{\partial{\dot{\boldsymbol{U}}}}{\partial t}+2\boldsymbol{N}(\boldsymbol{U},\dot{\boldsymbol{U}})={-}\boldsymbol{\nabla} \dot{P}+\frac{1}{{\textit{Re}}}\nabla^2\dot{\boldsymbol{U}}, \end{cases}

with homogeneous boundary conditions. Comparison with (4.1), which governs the small-amplitude perturbations to the base flow, shows that for $\gamma =0$,

(5.2)\begin{equation} \begin{pmatrix} u', v', w', p' \end{pmatrix}=\tau_0\begin{pmatrix} \dot{U}, \dot{V}, 0, \dot{P} \end{pmatrix} \end{equation}

is a valid two-dimensional perturbation to the base flow. (The amplitude $\tau _0\ll 1$ is introduced to ensure that the perturbations are sufficiently small to justify the linearisation that leads to (4.1).) Since $(\dot {\boldsymbol {U}},\dot {P})$ are time-periodic, the perturbations (5.2) are too, implying that they are neutrally stable, $\mu =1$, consistent with the numerical results shown in figure 2.

The perturbations (5.2) correspond to a small temporal shift in the flow field since

(5.3)\begin{equation} u(\boldsymbol{x},t) = U(\boldsymbol{r},t)+\tau_0 \dot{U}(\boldsymbol{r},t) = U(\boldsymbol{r},t+\tau_0)+O(\tau_0^2), \end{equation}

where we have used the Taylor expansion of $U(\boldsymbol {r},t+\tau _0)$. This reflects the fact that the two-dimensional time-periodic base flow is only determined up to an arbitrary temporal phase shift, here represented by $\tau _0$.

Given the explicit expression (5.2) for perturbations with zero wavenumber, $\gamma =0$, we now pose a perturbation expansion in the regime $0<\gamma \ll 1$. We start by noting that in this regime, the Floquet multiplier $\mu$ is real; therefore, the solution has the standing wave form (see (4.4))

(5.4)\begin{equation} (u' , v' , w', p')=\tau\left[\left(u_{p}, v_{p},0, p_{p}\right)\cos(\gamma z)- \left(0,0, w_{p}, 0\right)\sin(\gamma z)\right], \end{equation}

where $\tau (t)=\tau _0 \,{\rm e}^{\sigma t}$ and the subscript ‘$p$’ indicates that a function is $T$-periodic. We assume that $\tau _0$ is sufficiently small to ensure that $\tau (t) \ll 1$; this is consistent with the tacit assumption that the exponential growth of the instability has not increased its amplitude to a level that would invalidate the linearisation underlying the derivation of (4.1).

Substituting (5.4) into the linearised Navier–Stokes equations (4.1) yields

(5.5a,b and c)\begin{cases} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{p}-\gamma w_{p}= 0,\\ \frac{\partial{\boldsymbol{u}_{p}}}{\partial t}+2\boldsymbol{N}(\boldsymbol{U},\boldsymbol{u}_{p})={-}\boldsymbol{\nabla} p_{p}+\frac{1}{{\textit{Re}}} \left(\nabla^2\boldsymbol{u}_{p} -\gamma^2 \boldsymbol{u}_{p} \right)- \sigma \boldsymbol{u}_{p},\\ \frac{\mathcal{D} w_{p}}{\mathcal{D}t} ={-}\gamma p_{p}+\frac{1}{{\textit{Re}}}\left( \nabla^2w_{p} - \gamma^2 w_{p}\right) -\sigma w_{p},\end{cases}

where $\boldsymbol {u}_{p}(\boldsymbol {r},t)=(u_{p}, v_{p},0)$ and $\mathcal {D}/\mathcal {D}t$ is the linearised substantial derivative $\mathcal {D}/\mathcal {D}t=\partial /\partial t + (\boldsymbol {U}\boldsymbol {\cdot }\boldsymbol {\nabla })$.

Using the explicit solution for two-dimensional perturbations (5.2), we obtain that in the limit $\gamma \rightarrow 0$, $(u_{p}, v_{p}, p_{p})$ must tend to $(\dot {U}, \dot {V},\dot {P})$ while $\sigma$ and $w_{p}$ must both tend to 0. This initially suggests the following expansions for the $T$-periodic functions $\boldsymbol {u}_{p}(\boldsymbol {r},t)=(u_{p},v_{p},0)$, $w_{p}(\boldsymbol {r},t)$, $p_{p}(\boldsymbol {r},t)$, and the growth rate $\sigma$:

(5.6) \begin{equation} \begin{alignedat}{4} \boldsymbol{u}_{p}(\boldsymbol{r},t)=~&\dot{\boldsymbol{U}}(\boldsymbol{r},t)&~+~&\gamma \boldsymbol{u}_1(\boldsymbol{r},t)&~+~&\gamma^2 \boldsymbol{u}_2(\boldsymbol{r},t) &~+~& O(\gamma^3),\\ w_{p}(\boldsymbol{r},t)=~&~&~&\gamma w_1(\boldsymbol{r},t) &~+~& \gamma^2 w_2(\boldsymbol{r},t)&~+~& O(\gamma^3),\\ p_{p}(\boldsymbol{r},t)=~&\dot{P}(\boldsymbol{r},t)&~+~&\gamma p_1(\boldsymbol{r},t)&~+~&\gamma^2 p_2(\boldsymbol{r},t) & ~+~ & O(\gamma^3),\\ \sigma=~&~&~&\gamma \sigma_{1}&~+~&\gamma^2 \sigma_{2} & ~+~ & O(\gamma^3). \end{alignedat} \end{equation}

We now note that since for both signs of $\gamma$ the dominant standing-wave mode (5.4) is the same, $u_p$, $p_p$ and $\sigma$ must be even in $\gamma$ and $w_p$ must be odd. Hence,

(5.7)\begin{equation} \left(\begin{array}{c} u' \\ v' \\ w' \\ p' \end{array}\right)=\tau \left[\left(\begin{array}{c} \dot{U} +\gamma^2 u_2 + O(\gamma^4) \\ \dot{V} +\gamma^2 v_2 + O(\gamma^4)\\ 0 \\ \dot{P} +\gamma^2 p_2 + O(\gamma^4) \end{array}\right) \cos(\gamma z) -\left(\begin{array}{c} 0 \\ 0 \\ \gamma w_1 + O(\gamma^3) \\ 0 \end{array}\right) \sin(\gamma z)\right]. \end{equation}

To demonstrate the consistency of this expansion with the numerical results, we define the functions $\chi _1={\left \lVert w _{p}\right \rVert }/{\left \lVert u _{p}\right \rVert }$ and $\chi _2={\left \lVert v _{p}\right \rVert ^2}/{\left \lVert u _{p}\right \rVert ^2}$. The expansion (5.7) then implies that for $\gamma \ll 1$, we have

(5.8a,b)\begin{equation} \chi_1= \frac{\left\lVert w_1\right\rVert}{\left\lVert\dot{U}\right\rVert}\gamma+O(\gamma^3),\quad \chi_2=\frac{\left\lVert\dot{V}\right\rVert^2}{\left\lVert\dot{U}\right\rVert^2} \left[1+2\left(\frac{\langle v_2, \dot{V}\rangle}{\left\lVert\dot{V}\right\rVert^2}- \frac{\langle u_2, \dot{U}\rangle}{\left\lVert\dot{U}\right\rVert^2}\right)\gamma^2+ O(\gamma^4)\right], \end{equation}

where $\left \lVert {\cdot }\right \rVert$, $\langle {{\cdot },{\cdot }}\rangle$ are the $L^2$-norm and inner product, respectively (calculated for the yellow shaded region shown in figure 3).

The symbols in figure 4 show $\chi _1$ and $\chi _2$ computed from the numerical results; the continuous lines are the approximations $\chi _1^{[fit]}=k \gamma$ and $\chi _2^{[fit]}=c + a \gamma ^2$, where we fitted $a$, $c$ and $k$ using the numerical data for $\gamma =0, 0.05, 0.1$ and $0.15$. The numerical data can be seen to be well described by the predictions from (5.8a,b); the fitted constant $c$ differs by less than $1.2\,\%$ from the value ${\left \lVert \dot {V}\right \rVert ^2}/{\left \lVert \dot {U}\right \rVert ^2}$.

Figure 4. Plot of the ratios $\chi _1={\left \lVert w _{p}\right \rVert }/{\left \lVert u _{p}\right \rVert }$ and $\chi _2={\left \lVert v _{p}\right \rVert ^2}/{\left \lVert u _{p}\right \rVert ^2}$ at ${\textit {Re}}=220$. The symbols represent the values obtained from the numerical simulations; the solid lines are fits based on the functional form (5.8a,b).

Having established that the leading-order terms in the expansion (5.7) provide a good description of the three-dimensional perturbations, we note that the Taylor expansion employed to derive (5.3) now shows that

(5.9)\begin{equation} \left(\begin{array}{c} u(\boldsymbol{x},t) \\ v(\boldsymbol{x},t) \\ p(\boldsymbol{x},t) \end{array}\right)=\left(\begin{array}{c} U(\boldsymbol{r},t) + u'(\boldsymbol{x},t) \\ V(\boldsymbol{r},t) + v'(\boldsymbol{x},t) \\ P(\boldsymbol{r},t) + p'(\boldsymbol{x},t) \end{array}\right)=\left(\begin{array}{c} U(\boldsymbol{r},t+\tau\cos(\gamma z))+ O(\tau^2,\tau \gamma^2)\\ V(\boldsymbol{r},t+\tau\cos(\gamma z))+ O(\tau^2,\tau \gamma^2)\\ P(\boldsymbol{r},t+\tau\cos(\gamma z))+ O(\tau^2,\tau \gamma^2) \end{array}\right). \end{equation}

This implies that, to leading order, long-wavelength perturbations to the two-dimensional base flow self-organise so that the flow in each streamwise slice corresponds to the base flow at shifted times, where the amount of shift depends on the spanwise coordinate, $z$. This is illustrated in the conceptual sketch in figure 5.

Figure 5. Illustration of the time-shifting pattern for the three-dimensionally perturbed flow: the flow in each streamwise slice is given by the two-dimensional flow at a slightly different time; the time shift depends on the spanwise coordinate $z$. (a) Two-dimensional base flow within streamwise slices at slightly different times. (b) Resulting three-dimensional perturbed flow.

Furthermore, substituting (5.6) into (5.5) shows that the equation for $w_1$ is uncoupled from the other perturbations,

(5.10)\begin{equation} \frac{\mathcal{D} w_1}{\mathcal{D}t} ={-}\dot{P}+\frac{1}{{\textit{Re}}}\nabla^2w_1,\end{equation}

and hence the leading-order spanwise flow is driven exclusively by the pulsations of the base flow pressure, $\dot {P}(\boldsymbol {r},t)$.

The perturbation to the vorticity is given by

(5.11)\begin{equation} \left(\begin{array}{c} \omega_x' \\ \omega_y' \\ \omega_z' \end{array}\right)=\tau \left[\left(\begin{array}{c} 0 \\ 0 \\ \dot{\varOmega} + O(\gamma^2) \end{array}\right)\cos(\gamma z) + \left(\begin{array}{c} \gamma\left(\dot{V}-{\partial w_1}/{\partial y}\right) + O(\gamma^3)\\ \gamma\left({\partial w_1}/{\partial x}-\dot{U}\right) + O(\gamma^3)\\ 0 \end{array}\right) \sin(\gamma z)\right]. \end{equation}

This shows that for small wavenumbers, the perturbations to the vorticity are dominated by the spanwise component, $\omega '_z$, which is largest in regions where the time derivative of the base flow vorticity, $\dot {\varOmega }$, is large. This is consistent with the observation that, in the course of the mode A instability, the vortex cores in the base flow undergo considerable spanwise wavy deformations (here due to the $\cos (\gamma z)$ term); see, for example, Barkley & Henderson (1996) and Jiang et al. (2016a).

The comparison of the in-plane perturbation velocity for mode A at $\gamma =0$ (obtained by time differentiation of the base-flow solution) with cases at $\gamma \neq 0$ in figure 3 shows that the perturbation pattern in the vortex formation region is still qualitatively similar to $(\dot {U}, \dot {V})$ even when $\gamma$ is not small and ${\textit {Re}}>{\textit {Re}}_A$ (at least up to $\gamma =2.2$ and ${\textit {Re}}=220$ according to figure 3). Furthermore, figure 6 shows that the pattern of the perturbations remains unchanged even at lower ${\textit {Re}}$. Although, formally, the leading-order approximation is no longer valid when $\gamma$ is not small (e.g. in figure 2, it is evident that the leading order term does not describe the dependence $\sigma (\gamma )$ for large $\gamma$), the persistence of the time-shifting pattern in the vortex formation region indicates its significant role in the onset of mode A. This suggests that the spatial structure of the mode A instability can be explained by the mechanism for the formation of the time-shifting pattern discussed above.

Figure 6. Pattern of perturbations at ${\textit {Re}}=100,150$ and $\gamma =0, 0.8$, and $1.6$: perturbation energy $e$ (greyscale colour contours) and in-plane ($z=0$) perturbation velocity (arrows). Solid lines are the base flow vorticity isolines $\varOmega =\pm 1$. All plots are snapshots at $t=0.5T$, corresponding to the minimum of the lift coefficient. Perturbations at $\gamma =0$ are obtained by time differentiation of the base flow solution, see (5.2). One should not directly compare the magnitude of perturbations in the different cases; it is defined up to a constant factor which we adjusted manually to highlight the similarities and differences of the perturbations patterns.

The symmetry of mode A behind a circular cylinder is inherited from the two-dimensional base flow. Williamson (1996b) observed it experimentally and gave a physical explanation based on the suggested self-sustaining process. Barkley & Henderson (1996) extracted the symmetry relations by examining numerically obtained eigenfunctions from the linear stability analysis:

(5.12)\begin{equation} \begin{pmatrix} u_p\\ v_p\\ w_p\\ p_p \end{pmatrix}(x,y,t+T/2)=\begin{pmatrix} u_p\\ -v_p\\ w_p\\ p_p \end{pmatrix}(x,-y,t). \end{equation}

For a base flow with the symmetry (3.1), only two types of synchronous bifurcations to the three-dimensional flow are allowed (Marques, Lopez & Blackburn 2004; Blackburn, Marques & Lopez 2005): preserving (like mode A) and breaking (like mode B) the base-flow symmetry. In this context, the instabilities with the Floquet branch connected to the neutral mode at $\gamma =0$ belong to the former group: the neutral mode has the base-flow symmetry and, consequently, the time-shifting pattern inherits it as well (e.g. see relation (5.6)). However, it is not evident whether all the modes that exhibit the symmetry (5.12) are caused by a common physical mechanism.

We note that none of the above analysis relies on the geometry of the cylinder, implying that our results are equally applicable to flows past other bluff bodies for which mode A instabilities are observed. In particular, the time-shifting pattern is observed even when the base flow is non-symmetric, e.g. in the flow past an elliptic cylinder at an incidence angle (Rao et al. 2017) and in the flow past a rotating cylinder (Rao et al. 2015). Therefore, the time-shifting pattern can serve as an additional unifying characteristic of certain mode A type three-dimensional instabilities.

6. Physical mechanisms for flow instability

The previous section showed that in the small-wavenumber limit, small-amplitude three-dimensional perturbations to the two-dimensional time-periodic base flow are dominated by a simple time shifting of that base flow. Comparison against the numerical solution of the perturbation equations showed that this pattern persists up to wavenumbers at which the base flow becomes unstable to the mode A instability. The approach, therefore, successfully predicts the flow pattern at the onset of the three-dimensional instability, but it does not explain why these perturbations grow for a specific range of wavenumbers at fixed Reynolds number (e.g. at ${\textit {Re}}=220$, it corresponds to $1.1<\gamma <2.1$, as illustrated in figure 2).

To address this issue, we now analyse the various physical mechanisms that affect the growth or decay of such three-dimensional perturbations. For this purpose, we define the in-plane perturbation velocity $\boldsymbol {v}(\boldsymbol {r}, t)=(v_x,v_y, 0)$ and vorticity $\boldsymbol {\zeta }(\boldsymbol {r}, t)=(\zeta _x,\zeta _y, 0)$ by the relations

(6.1) \begin{equation} \left.\begin{gathered} \boldsymbol{u}'\left(\boldsymbol{x}, t\right)=(v_x,v_y,0)\cos(\gamma z)+\left(0, 0, v_z\right)\sin(\gamma z),\\ \boldsymbol{\omega}'\left(\boldsymbol{x}, t\right)=(\zeta_x,\zeta_y,0)\sin(\gamma z)+ \left(0,0, \zeta_z\right)\cos(\gamma z). \end{gathered}\right\} \end{equation}

Using the definition of the three-dimensional vorticity, $\boldsymbol {\omega }'=\boldsymbol {\nabla }\times \boldsymbol {u}'$, and the fact that the three-dimensional perturbation velocity $\boldsymbol {u}'$ is divergence-free shows that these two-dimensional fields are related via

(6.2)\begin{equation} \boldsymbol{\zeta}=\gamma (v_{y}, -v_{x},0)+ \frac{1}{\gamma}\left(-\frac{\partial{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}}}{\partial y}, \frac{\partial{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}}}{\partial x},0\right). \end{equation}

The rate of change of the two-dimensional perturbation vorticity $\boldsymbol {\zeta }$ is governed by the linearised vorticity transport equation

(6.3)\begin{equation} \frac{\mathcal{D} \boldsymbol{\zeta}}{\mathcal{D}t}= \underbrace{\boldsymbol{\mathsf{E}} \boldsymbol{\cdot} \boldsymbol{\zeta}}_{stretching} \underbrace{+\frac{1}{2}\boldsymbol{\varOmega}\times\boldsymbol{\zeta}}_{rigid\ rotation} \underbrace{+\frac{1}{{\textit{Re}}}\left(\underbrace{\nabla^{2} \boldsymbol{\zeta}}_{\textit{in-plane}}- \underbrace{\gamma^{2} \boldsymbol{\zeta}}_{spanwise}\right)}_{viscous\ diffusion} \underbrace{-\gamma \varOmega \boldsymbol{v}}_{tilting}, \end{equation}

where $\boldsymbol {\varOmega } = (0,0,\varOmega )$. Each term on the right-hand side of (6.3) has a clear physical interpretation and explains the material rate of change of the perturbation vorticity $\boldsymbol {\zeta }$ in terms of vortex stretching by the base flow rate-of-strain field $\boldsymbol{\mathsf{E}}$; the re-orientation of the vorticity vector by the rigid body rotation of fluid particles in the base flow; the in-plane and spanwise viscous diffusion of the perturbation vorticity; and the tilting of the base flow vortex due to spanwise shear. See Appendix B for more details.

To facilitate the subsequent analysis, we combine (6.2) and (6.3) by exploiting that the perturbation vorticity, $\boldsymbol {\omega }'$, is divergence-free and that for an incompressible fluid, $\nabla ^2\boldsymbol {u}'=-\boldsymbol {\nabla }\times \boldsymbol {\omega }'$. This implies that

(6.4)\begin{equation} \nabla^2\boldsymbol{v}-\gamma^2\boldsymbol{v}=\gamma\boldsymbol{\zeta}_{\bot}+\frac{1}{\gamma}\boldsymbol{\zeta}_{\Delta}, \end{equation}

where

(6.5a,b)\begin{equation} \boldsymbol{\zeta}_{\bot}(\boldsymbol{r}, t)=\left(\zeta_{y},-\zeta_{x},0 \right),\quad \boldsymbol{\zeta}_{\Delta}(\boldsymbol{r}, t) =\left(-\boldsymbol{\nabla}\boldsymbol{\cdot}\frac{\partial\boldsymbol{\zeta}}{\partial y}, \boldsymbol{\nabla}\boldsymbol{\cdot}\frac{\partial\boldsymbol{\zeta}}{\partial x},0\right). \end{equation}

The screened Poisson equation (6.4) determines the in-plane velocity perturbation $\boldsymbol {v}$ in terms of the in-plane perturbation to the vorticity, $\boldsymbol {\zeta }$. An explicit relation between the two fields can therefore be obtained by introducing the Green's function $G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$, which satisfies

(6.6) \begin{equation} \left.\begin{gathered} \nabla^2G_\gamma-\gamma^2G_\gamma=\delta(\boldsymbol{r}-\boldsymbol{r}'),\\ G_\gamma=0\quad \text{at } r=0.5,\\ G_\gamma\rightarrow0\quad \text{as } r\rightarrow\infty. \end{gathered}\right\} \end{equation}

We show in Appendix C that the solution to (6.6) is given by

(6.7)\begin{equation} G_\gamma(\boldsymbol{r}, \boldsymbol{r}')={-}\frac{1}{2{\rm \pi}}K_0(\gamma|\boldsymbol{r}- \boldsymbol{r}'|)+\frac{1}{2{\rm \pi}}\sum_{m={-}\infty}^{\infty} \frac{I_m(\gamma/2)K_m(\gamma r)K_m(\gamma r')}{K_m(\gamma/2)}\cos m(\varphi-\varphi'), \end{equation}

where $I_m(r)$ and $K_m(r)$ are the modified Bessel functions of the first and second kind; $\boldsymbol {r}=r(\cos \varphi, \sin \varphi )$ and $\boldsymbol {r}'=r'(\cos \varphi ', \sin \varphi ')$.

Using this expression, the in-plane perturbation velocity $\boldsymbol {v}$ is given by

(6.8)\begin{equation} \boldsymbol{v}(\boldsymbol{r},t)=\int_{D}\gamma G_\gamma(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{\zeta}_{\bot}\left(\boldsymbol{r}',t\right)+\frac{1}{\gamma}G_\gamma(\boldsymbol{r},\boldsymbol{r}') \boldsymbol{\zeta}_{\Delta}\left(\boldsymbol{r}',t\right){\rm d} \boldsymbol{r}', \end{equation}

where $D$ is the exterior of the cylinder. Substituting this into (6.3) yields

(6.9)\begin{align} \frac{\mathcal{D} \boldsymbol{\zeta}}{\mathcal{D}t} &= \underbrace{\boldsymbol{\mathsf{E}} \boldsymbol{\cdot} \boldsymbol{\zeta}}_{stretching} +\underbrace{\frac{1}{2}\boldsymbol{\varOmega}\times\boldsymbol{\zeta}}_{rigid\ rotation} +\underbrace{\frac{1}{{\textit{Re}}}\left(\underbrace{\nabla^{2} \boldsymbol{\zeta}}_{{in\text{-}plane}}- \underbrace{\gamma^{2} \boldsymbol{\zeta}}_{spanwise}\right)}_{viscous\ diffusion} \nonumber\\ &\quad \underbrace{-\,\varOmega \int_{D}\gamma^2 G_\gamma(\boldsymbol{r}, \boldsymbol{r}') \boldsymbol{\zeta}_{\bot}\left(\boldsymbol{r}',t\right)+G_\gamma(\boldsymbol{r}, \boldsymbol{r}')\boldsymbol{\zeta}_{\Delta}\left(\boldsymbol{r}',t\right){\rm d} \boldsymbol{r}'}_{tilting}, \end{align}

which describes the evolution of perturbations to the flow entirely in terms of the perturbations to the in-plane vorticity, $\boldsymbol {\zeta }$. The equation shows that the first three physical mechanisms are local in the sense that their contribution to the rate of change of $\boldsymbol {\zeta }$ depends only on $\boldsymbol {\zeta }$ or its spatial derivatives. Conversely, tilting is a global effect – the rate of change of $\boldsymbol {\zeta }$ due to the final term depends on $\boldsymbol {\zeta }$ and its derivatives throughout the domain. Furthermore, (6.9) shows how variations in the two parameters ${\textit {Re}}$ and $\gamma$ affect the various mechanisms. The wavenumber only affects the spanwise diffusion and the tilting mechanism. The effect of variations in the Reynolds number is more subtle: it has a direct effect on the strength of the viscous diffusion but also affects the base flow and, thus, the stretching, rigid rotation and tilting mechanisms (via $\boldsymbol{\mathsf{E}}$ and $\varOmega$). We will now analyse the importance of these mechanisms in detail.

6.1. Effect of the viscous diffusion and the base flow

The Reynolds number simultaneously affects the base flow and the intensity of the in-plane and spanwise viscous diffusion. To study the contribution of these three effects separately, we replace the Reynolds number in front of the diffusion terms in (6.9) by ${\textit {Re}}'$ and ${\textit {Re}}''$, and thus write the evolution equation for the in-plane perturbation to the vorticity, $\boldsymbol {\zeta }$, as

(6.10)\begin{align} \frac{\mathcal{D} \boldsymbol{\zeta}}{\mathcal{D}t}&= \underbrace{\boldsymbol{\mathsf{E}} \boldsymbol{\cdot} \boldsymbol{\zeta}}_{stretching} +\underbrace{\frac{1}{2}\boldsymbol{\varOmega}\times\boldsymbol{\zeta}}_{rigid\ rotation} +\underbrace{\underbrace{\frac{1}{{\textit{Re}}'}\nabla^{2} \boldsymbol{\zeta}}_{{in\text{-}plane}}- \underbrace{\frac{\gamma^{2}}{{\textit{Re}}''} \boldsymbol{\zeta}}_{spanwise}}_{viscous\ diffusion} \nonumber\\ &\quad \underbrace{-\,\varOmega \int_{D}\gamma^2 G_\gamma(\boldsymbol{r}, \boldsymbol{r}')\boldsymbol{\zeta}_{\bot}\left(\boldsymbol{r}',t\right)+ G_\gamma(\boldsymbol{r}, \boldsymbol{r}')\boldsymbol{\zeta}_{\Delta}\left(\boldsymbol{r}',t\right){\rm d} \boldsymbol{r}'}_{tilting}. \end{align}

Here the ${\textit {Re}}$-dependent base flow affects the base-flow rate-of-strain tensor, $\boldsymbol{\mathsf{E}}$, and the base-flow vorticity, $\varOmega$.

Figure 7 illustrates the contributions that the mechanisms discussed so far make to the destabilisation of the flow as the Reynolds number is increased from 180 to 200. The two solid lines show the Floquet multipliers $\mu$ for the actual flow (i.e. when ${\textit {Re}}={\textit {Re}}'={\textit {Re}}''$) at Reynolds numbers ${\textit {Re}}=180$ and $200$. The remaining broken lines show the destabilising effects of the modification only in the base flow (blue line, ${\textit {Re}}=200,{\textit {Re}}'={\textit {Re}}''=180$), in-plane viscous diffusion (red line, ${\textit {Re}}'=200,{\textit {Re}}={\textit {Re}}''=180$) and spanwise viscous diffusion (green line, ${\textit {Re}}''=200,{\textit {Re}}={\textit {Re}}'=180$). (We obtained the curves corresponding to distinct ${\textit {Re}}$ and ${\textit {Re}}'$ by modifying the input data, feeding our stability code with the pre-computed base flow at ${\textit {Re}}$, while using ${\textit {Re}}'$ in the stability equations; the impact of the spanwise viscous diffusion (${\textit {Re}}''$) was assessed explicitly, see below.) For the wavelengths over which the mode A instability arises, the modification to the base flow and the in-plane viscous diffusion can be seen to have a considerable (and comparable) effect on the destabilisation of the flow, whereas the spanwise viscous diffusion only plays a minor role in this process.

Figure 7. Influence of the Reynolds number on the dominant Floquet multiplier near the onset of instability (${\textit {Re}}_A\approx 190$). Two solid black lines correspond to the actual Floquet multiplier $\mu$ (${\textit {Re}}={\textit {Re}}'={\textit {Re}}''$), other lines represent the Floquet multiplier obtained as a result of independent variation of the base flow (${\textit {Re}}$; blue), in-plane (${\textit {Re}}'$; red) and spanwise (${\textit {Re}}''$; green) viscous diffusion.

The overall effect of the spanwise viscous diffusion can be taken into account explicitly using the change of variables

(6.11)\begin{equation} \boldsymbol{\tilde{\zeta}}(\boldsymbol{r}, t) = \boldsymbol{\zeta}(\boldsymbol{r}, t) \exp(\gamma^2t/{\textit{Re}}), \end{equation}

which transforms (6.9) into an identical equation for $\boldsymbol {\tilde {\zeta }}$, but with the spanwise diffusion term removed. Equation (6.11), therefore, implies that in the absence of spanwise viscous diffusion, the Floquet multiplier $\mu$ would change to

(6.12)\begin{equation} \tilde{\mu} = \mu \exp(\gamma^2 T/Re) >\mu, \end{equation}

meaning that spanwise viscous diffusion is always stabilising, and that it does not have an effect on the spatial pattern of the perturbations. An increase in Reynolds number or a decrease in wavenumber both reduce the stabilising effect of the spanwise viscous diffusion. (We note that the period of the vortex shedding, $T$, also depends on the Reynolds number; however, for the regime considered here, $T$ decreases with ${\textit {Re}}$ and, therefore, does not affect our statement.)

6.2. Effect of the tilting mechanism

Figure 2 shows that the dependence of the dominant Floquet multiplier, $\mu$, on the wavenumber $\gamma$ is non-monotonic: $\mu (\gamma =0)=1$, corresponding to the neutral stability of the base flow to the time-shifting pattern discussed in § 5. The dominant Floquet multiplier then decreases with increasing $\gamma$ before it rises again, ultimately leading to the onset of the mode A instability for $1.1<\gamma <2.1$ at ${\textit {Re}}=220$.

Only the tilting and spanwise diffusion depend on the wavenumber and, as discussed above, the latter effect is always stabilising; more so, in fact, as $\gamma$ is increased. The destabilisation of the base flow at sufficiently large $\gamma$ must, therefore, be due to the tilting of the base flow vorticity due to spanwise shear ($\boldsymbol {v}\cos (\gamma z)$). Equation (6.8) shows that this mechanism is a non-local effect: $\boldsymbol {v}$ is induced by the perturbation to the in-plane vorticity, $\boldsymbol {\zeta }$, everywhere in the flow, in a manner similar to Biot–Savart induction.

The strength of this non-local interaction is defined by two kernel functions $G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$ and $\gamma ^2 G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$ in (6.9). They are determined by the problem geometry and describe how much the perturbations at point $\boldsymbol {r}$ are affected by the in-plane vorticity and its second derivatives at point $\boldsymbol {r}'$. Both kernel functions are singular at $\boldsymbol {r} = \boldsymbol {r}'$ and decay with an increase in $|\boldsymbol {r} - \boldsymbol {r}'|$. We illustrate the spatial variation of $G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$ in figure 8(b) for the case where $\boldsymbol {r}$ (identified by the red star symbol) is located at the instantaneous local maximum of the perturbation vorticity in the flow shown in figure 8(a). An increase in $\gamma$ makes both kernel functions more localised, as illustrated by the orange isolines, along which $G_\gamma (\boldsymbol {r}, \boldsymbol {r}') = -10^{-1.6}$. We note that an increase in $\gamma$ causes $G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$ to decrease throughout the domain; conversely, the magnitude of $\gamma ^2 G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$ increases in the vicinity of $\boldsymbol {r}$, enhancing the influence that the vorticity in the proximity of a given point has on the growth of the perturbations at that point.

Figure 8. The contribution of perturbation distribution to their growth or decay at the local maximum of $\zeta$ (marked with the star symbol) through the tilting mechanism at ${\textit {Re}}=220$, various $\gamma$ and $t=0.44T$. (a) In-plane perturbation vorticity: $\log _{10}\zeta$. (b) Kernel function $\log _{10} |G_\gamma (\boldsymbol {r},\boldsymbol {r}')|$. (c) Contribution to the tilting mechanism at point $\boldsymbol {r}$ (the star symbol): $\mathcal {T}(\boldsymbol {r},\boldsymbol {r}',t)$. The orange lines highlight the isoline of kernel function $G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$, shown in panel (b). Solid lines in panels (a,c) are isolines $\kappa =1$ (the boundaries of the elliptic regions). Perturbation vorticity is normalised so that the local maximum of $\zeta$ (star symbol) equals 1.

To determine the effect of the tilting mechanism on the actual growth rate of perturbations in a given flow, we multiply (6.9) by $\boldsymbol {\zeta }$ to obtain

(6.13)\begin{equation} \frac{1}{2}\frac{\mathcal{D} \zeta^2}{\mathcal{D}t}= \underbrace{\boldsymbol{\zeta}\boldsymbol{\cdot}\boldsymbol{\mathsf{E}} \boldsymbol{\cdot} \boldsymbol{\zeta}}_{stretching} \underbrace{+\frac{1}{{\textit{Re}}}\left(\boldsymbol{\zeta}\boldsymbol{\cdot}\nabla^{2} \boldsymbol{\zeta}-\gamma^{2} \zeta^2\right)}_{viscous\ diffusion} \underbrace{+ \int_{D}\mathcal{T}(\boldsymbol{r},\boldsymbol{r}', t)\,{\rm d} \boldsymbol{r}'}_{tilting}, \end{equation}

where $\zeta = |\boldsymbol {\zeta }|$ and

(6.14)\begin{equation} \mathcal{T}(\boldsymbol{r},\boldsymbol{r}',t) ={-}G_\gamma(\boldsymbol{r}, \boldsymbol{r}') \varOmega(\boldsymbol{r}, t)\boldsymbol{\zeta}(\boldsymbol{r}, t)\boldsymbol{\cdot}\left[\gamma^2 \boldsymbol{\zeta}_{\bot}\left(\boldsymbol{r}',t\right)+\boldsymbol{\zeta}_{\Delta}\left(\boldsymbol{r}',t\right)\right]. \end{equation}

When $\mathcal {T}(\boldsymbol {r},\boldsymbol {r}',t)$ is positive/negative, the perturbations in the vicinity of point $\boldsymbol {r}'$ tend to increase/decrease the magnitude of perturbations $\zeta$ at point $\boldsymbol {r}$.

Let us consider examples of the distribution of $\mathcal {T}$ at various $\gamma$ and fixed ${\textit {Re}}=220$ at time $t=0.44T$, which corresponds to the early stage of perturbation development in the forming vortex (cf. the illustration of the entire cycle in figure 9a). The greyscale contours in figure 8(a) illustrate the magnitude of the perturbation to the in-plane vorticity (on a logarithmic scale), normalised so that its local maximum at the location indicated by the star (the local maximum of $\zeta$ in the forming vortex) is equal to 1 in all cases. An increase in $\gamma$ leads to a strong increase in the magnitude of $\zeta$ as the perturbation develops; compare, e.g. the magnitude of $\zeta$ in the braid region for $\gamma =0.4$ and $\gamma =2.2$. This increase of $\zeta$ (and the derived quantities, $\boldsymbol {\zeta }_{\bot }$ and $\boldsymbol {\zeta }_{\Delta }$) is counteracted by the increasing localisation of the kernel functions, as illustrated by the isolines $G_\gamma (\boldsymbol {r}, \boldsymbol {r}') = -10^{-1.6}$ from figure 8(b). (The contribution of $\gamma ^2 G_\gamma (\boldsymbol {r}, \boldsymbol {r}') \boldsymbol {\zeta }_{\bot }$ to the tilting integral in (6.13) turns out to be negligible, so the isoline of $G_\gamma (\boldsymbol {r}, \boldsymbol {r}')$ gives a good indication of the domain of influence for the tilting mechanism.)

Figure 9. Local growth of perturbations at ${\textit {Re}}=220$ and $\gamma =1.6$: (a) in-plane perturbation vorticity $\log _{10}\zeta$; (b) local maximum $\zeta _{max}(t)$ (red line), which corresponds to the star symbol in panel (a). In panel (b), we also show similar curves for stable cases at ${\textit {Re}}=220$ (dashed black lines) and ${\textit {Re}}=50$ (dotted blue line). Solid lines in panel (a) are isolines $\kappa =1$ (the boundaries of the elliptic regions). Perturbation vorticity is normalised so that at $t=0.44T$, the local maximum of $\zeta$ (star symbol) equals one.

The red/blue colours (representing positive/negative values of $\mathcal {T}$) in figure 8(c) highlight which regions in the flow destabilise/stabilise the flow at point $\boldsymbol {r}$ and thus identify the significant non-local interactions of perturbations in various flow subregions. A key feature in figure 8(c) is that previously formed perturbations outside the forming vortex (particularly in the hyperbolic region) have a noticeable effect on the development of newly created perturbations in the vortex core. This raises questions about the validity of simplified models that attempt to explain the instability based on isolated flow features.

Figure 8(c) shows that, as expected, when $\gamma$ increases, the non-local interactions become more localised: the contribution to the growth of $\zeta$ at point $\boldsymbol {r}$ weakens more rapidly with the distance, even though surrounding perturbations become more intense at higher $\gamma$ (cf. figure 8a). Incidentally, another confirmation of more localised interactions of perturbations can be found in figures 3 and 6, where, with an increase in $\gamma$, the induced in-plane perturbation velocity, $\boldsymbol {v}$, can be seen to become smaller far from the regions where the perturbations in $\zeta$ concentrate; cf. (6.8).

Thus, the range of wavenumbers, $\gamma$, for which three-dimensional perturbations are unstable is determined by the tilting mechanism (recall that the uniformly stabilising action of the spanwise diffusion is explicitly taken into account by (6.12)). The tilting mechanism operates via non-local interactions between perturbations in different parts of the domain. The strength of these interactions is controlled by the spanwise wavenumber $\gamma$ through the kernel functions $G_\gamma$ and $\gamma ^2 G_\gamma$, which become more localised with an increase in $\gamma$.

6.3. Local growth and feedback

We will now use these results to analyse the subtle balance between the local growth of perturbations and their feedback on the newly developing perturbations that are generated during the next period of the time-periodic base flow. As discussed in § 1, our analysis can be confined to the vortex formation region, recognised as the origin of three-dimensional instability (Barkley 2005; Giannetti et al. 2010).

For this purpose, figure 9(a) shows the time evolution of the three-dimensional perturbations, characterised by logarithmic contours of the perturbation to the in-plane vorticity, $\zeta$, over one period of the time-periodic base flow for ${\textit {Re}}=220$ and $\gamma =1.6$. We note that for these values, the flow is unstable with a Floquet multiplier of $\mu = 1.25$. The symbols in figure 9(a) indicate the location $\boldsymbol {r}_M(t)$ of the local maximum of $\zeta$, which we follow from its inception in the forming vortex (figure 9ai) as it is swept through the domain (figure 9aii–viii). We normalised the perturbation so that at $t=0.44T$, the maximum is equal to 1, and show the subsequent evolution of $\zeta _{max}(t) = \zeta (\boldsymbol {r}_M(t),t)$ using the red line in figure 9(b).

The maximum is initially located in the elliptic region where the flow is dominated by rotation (figure 9ai–iv). This is the region in which the elliptic instability is commonly assumed to operate, and the perturbations can be seen to grow rapidly. At $t \approx 0.97T$, $\boldsymbol {r}_M(t)$ crosses the boundary of the elliptic region (identified by the white solid lines in figure 9a) and enters the hyperbolic region, where the flow is dominated by strain. Interestingly, $\zeta _{max}(t)$ can be seen to undergo a second phase of strong growth when it enters this region, and over one period of the time-periodic base flow (starting at $t=0.44T$), $\zeta _{max}(t)$ increases by a factor of $6$.

In this particular case, the tilting-induced feedback to the region where the next maximum is formed is strong enough to result in an overall flow instability with a Floquet multiplier of $\mu = 1.25$. Hence, at $t=1.44T$, when the cycle is about to repeat itself, the amplitude of the next maximum of the perturbation in the vortex formation region is $\mu = 1.25$ times bigger than during the previous cycle. This process then repeats itself with every new cycle (at least until nonlinearities, which are not included in our analysis, become significant).

The other lines in figure 9(b) illustrate the temporal evolution of the maximum perturbation for other values of the parameters $\gamma$ and ${\textit {Re}}$. An increase in $\gamma$ generally leads to stronger local growth, but this does not necessarily strengthen the overall instability. For instance, when $\gamma = 2.2$, $\zeta _{max}(t)$ increases by a factor of 8 over one period. Yet, despite this more rapid local growth, the overall flow is actually restabilised ($\mu = 0.92$). This can be explained by the weakening of the feedback using the analysis presented in § 6.2. Indeed, the variation of $\gamma$ primarily affects the flow (and the amount of feedback in particular) through the tilting mechanism, i.e. through non-local interactions of perturbations. Figure 8 illustrates that these interactions become more localised as $\gamma$ increases. Consequently, although the perturbations become more intense, they are not sufficiently felt during the formation of the next local maximum, resulting in a weakened feedback.

Similarly, it is interesting to note that even when the Reynolds number is reduced to ${\textit {Re}} = 50 < {\textit {Re}}_A\approx 190$, a regime where the flow is strongly stable, with $\mu = 0.28$, $\zeta _{max}(t)$ still undergoes a noticeable local growth. However, it remains too weak to provide sufficient positive feedback to the formation of perturbations during the next cycle.

7. Conclusions

We studied the onset of three-dimensional mode A instability in the near wake behind a circular cylinder. Our analysis showed that long-wavelength perturbations organise in a time-shifting pattern so that the in-plane part $(u, v, p)(x,y,z,t)$ of the perturbed velocity field is connected to the two-dimensional base flow $(U, V, P)(x,y,t)$ by the relation

(7.1)\begin{equation} (u, v, p)(x,y,z,t)\approx(U, V, P)\left(x,y,t+\tau\cos(\gamma z)\right), \end{equation}

where $\tau (t)$ describes the exponential growth of the instability, here assumed to remain small enough to justify the use of linearised equations for the perturbations. The spanwise component of the velocity perturbation is driven by fluctuations in the base flow pressure, $\dot {P}$.

While these predictions are based on the assumption of a small wavenumber, $\gamma \ll 1$, comparisons against the results from a numerical Floquet analysis showed that they provide a good qualitative description of the three-dimensional flow over a range of wavenumbers and Reynolds numbers, including the regime where the flow changes from being stable to being unstable to three-dimensional perturbations.

We therefore analysed the mechanisms which control the growth or decay of three-dimensional perturbations and established their dependence on the two non-dimensional parameters (the Reynolds number ${\textit {Re}}$ and the wavenumber $\gamma$). It turned out that near the onset of the instability (${\textit {Re}}\sim {\textit {Re}}_A$), changes in the base flow and the intensity of in-plane viscous diffusion with ${\textit {Re}}$ are essential in flow destabilisation (having comparable contributions); in contrast, the relative effect of spanwise viscous diffusion is negligible. The analysis also highlighted the crucial role played by the tilting mechanism, which operates via non-local interactions, similar to Biot–Savart induction. We characterised its domain of influence using a Green's function-based approach, which allowed us to rationalise the non-trivial dependence of the growth rate on the wavenumber $\gamma$: for $\gamma =0$, the base flow is neutrally stable, corresponding to a Floquet multiplier of $\mu =1$; for small positive values of $\gamma$, the growth rate of three-dimensional perturbations becomes negative, $\mu < 1$; it then reaches a minimum before increasing, finally reaching the onset of the mode A instability when $\mu > 1$. We attributed this behaviour to two competing effects: an increase in $\gamma$ leads to a more rapid local growth of the perturbations as they are swept along by the flow. However, the action of the tilting mechanism becomes more localised, weakening the feedback from existing perturbations on the perturbations that are being generated during the next period of the time-periodic base flow.

While our analysis was performed for flows past circular cylinders, none of the technical details rely on the specific geometry and, therefore, can be useful in studying other instabilities that transform into the time-shifting mode as the spanwise wavenumber tends to zero. For example, additional simulations (not shown here) indicate that the long-wavelength modes $\hat {\textrm {A}}$ and G behind an elliptic and rotating cylinder exhibit the same time-shifting pattern. This is consistent with the Floquet analyses of Rao et al. (2013) and Leontini et al. (2015).