2.1. Computational set-up

In Wu et al. (Reference Wu, Axtmann and Rist2021), two sets of boundary-layer flows with rotating cylindrical roughness elements, i.e. co-rotating and counter-rotating cylinder pairs, are investigated by linear stability theory. It is demonstrated that the velocity streaks generated by ‘positive’ counter-rotating pairs, which accelerate the flow between two immediate neighbours, are able to effectively stabilise TS instabilities, since they create a high-speed streak in the centre. In the present work, the corresponding boundary-layer transition induced by positive counter-rotating cylindrical roughness elements is investigated. The numerical set-up is illustrated in figure 1. Hereafter, all physical quantities are non-dimensionalised by the dimensional roughness height $\bar {k}=0.01({\rm m})$ and the incoming free-stream velocity $u_\infty =0.937({\rm m\ s}^{-1})$. The roughness elements with height $k$ and diameter $D$ are placed at the location $x_k= 99.2$ from the flat-plate leading edge, under zero-pressure gradient in the streamwise direction. At this position, the non-dimensional displacement thickness is $\delta ^*/k = 0.6883$. The roughness-height-based Reynolds number is $Re_{kk} = u(k) k/\nu = 465.8$, which is intentionally chosen to be below the critical transitional Reynolds number as reviewed by Von Doenhoff & Braslow (Reference Von Doenhoff and Braslow1961). Here, $u(k)$ is the undisturbed Blasius velocity at the height of the roughness element. The free-stream-velocity-based Reynolds number is $Re_k= u_\infty k /\nu = 620$. The aspect ratio of the cylinder stub is defined as $\eta = D/k= 1$. The counter-rotating cylinder stubs are placed at spanwise locations $z=\pm 2$ such that they are spaced by $\lambda =4$, and the spanwise spacing between roughness pairs is $\varLambda = 12$. The origin of the local $x$,$y$,$z$-coordinate system is placed at the centre of the roughness pair. Depending on the rotational direction of the roughness pair, different types of downstream streaks are created. For the co-rotating roughness case, the generated downstream vortices rotate in the same direction, which will counteract the momentum exchange effect of neighbouring vortices. In contrast, this momentum exchange effect is intensified in the case of counter-rotating roughness elements. The strength of the generated downstream vortices depends on the rotation speed of the roughness elements. This is measured by the ratio of the tangential velocity at the top of the cylindrical roughness element to the incoming local velocity, i.e.

(2.1)\begin{equation} \varOmega_u= \frac{\varOmega D}{2u(k)}, \end{equation}

where $\varOmega$ is the angular velocity of the cylinder. In this work, the studied rotation speeds are in the range of $\varOmega _u = 0 \sim 2$.

Figure 1. Set-up of counter-rotating roughness pair embedded into a flat-plate boundary layer. High- and low-momentum flows induced by streamwise vortices are coloured in red and blue, respectively.

2.2. Steady and unsteady flow computation

The characteristic flow structures induced by a rotating cylindrical roughness element are a dominating inner vortex (DIV) encircled by a secondary inner vortex (SIV), as described in Wu et al. (Reference Wu, Axtmann and Rist2021) and Wu & Rist (Reference Wu and Rist2022) and also shown in figure 3 further down. The behaviour of the downstream vortical structures in response to different rotation rates is investigated in laminar flow. For this purpose, steady base flows have been computed with the method of selective frequency damping (SFD, Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006). In the SFD approach, a filtered state $\tilde {\boldsymbol {u}}$ is introduced and its evolutionary equation is solved together with the Navier–Stokes equations. This leads to the following incompressible non-dimensional governing equations:

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

(2.2b)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p +\frac{1}{Re_k}\triangle \boldsymbol{u} - \chi (\boldsymbol{u} - \tilde{\boldsymbol{u}}), \end{gather}$$

(2.2c)$$\begin{gather}\frac{\partial \tilde{\boldsymbol{u}}}{\partial t} = \omega_c (\boldsymbol{u} - \tilde{\boldsymbol{u}}), \end{gather}$$

where $Re_k = \bar {u}_\infty \bar {k} / \bar {\nu }$ is the free-stream-velocity-based Reynolds number, $\omega _c$ is the filter cutoff circular frequency and $\chi$ is the feedback control coefficient. Here, an overbar denotes a dimensional value, so that $t=\bar {t} \bar {u}_\infty / \bar {k}$ and $p=\bar {p}/(\bar {\rho } \bar {u}^2)$. The additional forcing term on the right-hand side of the momentum equation acts as a temporal low-pass filter. From theoretical analysis it is known that the cutoff frequency $\omega _c$ should be lower than the lowest eigenfrequency of instabilities, while the feedback control coefficient $\chi$ should be higher than the growth rate of that instability (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006).

The above governing equations are implemented and solved with the open source OpenFOAM solver icoSfdFoam (Wu et al. Reference Wu, Axtmann and Rist2021), which solves the incompressible Navier–Stokes equations using the pressure-implicit with splitting of operators algorithm. The unsteady terms of the governing equations are discretised with an implicit backwards scheme. For investigating the spatio-temporal evolution of instabilities and transition in the presence of rotating cylindrical roughness elements, unsteady simulations are performed with the SFD mechanism switched off. A multi-directional cell-limited gradient scheme is used to approximate the gradient term. The divergence term is handled with a limited linear scheme. A second-order deferred correction scheme is used for the Laplacian term. The pressure equation is solved by the geometric algebraic multi-grid solver and the velocity equation by the preconditioned bi-conjugate gradient solver.

A Blasius boundary-layer velocity profile is prescribed at the inlet according to the distance from the leading edge of the flat plate. No-slip wall boundary conditions are applied at the bottom wall. The velocity at the cylinder wall is calculated with $\boldsymbol {u}_w=\boldsymbol {\varOmega } \times \boldsymbol {r}$, with $\boldsymbol {r}$ the vector from the rotational axis of the cylinder to the wall surfaces and $\boldsymbol {\varOmega }$ the angular velocity vector. The standard OpenFoam implementation of the zero-gradient boundary condition ($\partial p / \partial n = 0$) is applied to the pressure condition for all wall boundaries. Here, $n$ denotes the direction normal to the boundary. At the top of the integration domain, a Dirichlet boundary is imposed for pressure $p$, and the following directional mixed boundary conditions are imposed for velocity:

(2.3a)$$\begin{gather} \frac{u}{u_\infty} = 1, \end{gather}$$

(2.3b)$$\begin{gather}\frac{\partial v}{\partial n} = \frac{\partial w}{\partial n} = 0. \end{gather}$$

At the outlet, any possible reflection is eliminated by using the following advective boundary condition:

(2.4)\begin{equation} \frac{{\rm D} q}{{\rm D}t} = \frac{\partial q}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} q = 0, \end{equation}

where $q(\boldsymbol {x},t)= \{u,v,w,p \}(\boldsymbol {x},t)$ are the flow quantities. In addition, the following sponge zone utilising the SFD mechanism is applied at the outlet. The damping strength is gradually increased towards the outlet

(2.5a)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} = \mathcal{N}\boldsymbol{u} - \chi B(x^*) (\boldsymbol{u}- \bar{\boldsymbol{u}}), \end{gather}$$

(2.5b)$$\begin{gather}B(x^*) = \begin{cases} 1 - \cos\left(\pi \dfrac{x^*-x_s}{x_e-x_s}\right) , & x^* \in [x_s;x_e], \\ 1, & x^* \in [x_e,\infty], \end{cases} \end{gather}$$

where $\mathcal {N}$ is the Navier–Stokes operator, $\bar {\boldsymbol {u}}$ is the mean flow, $B(x^*)$ is the ramp function and $x_s$, $x_e$ are the streamwise starting and ending positions of the ramp, respectively. Also, $\chi$ is the feedback control coefficient from the SFD solver, i.e. (2.2b). Both velocity and pressure at transverse planes are set to periodic conditions.

Investigations of the convergence of the above procedures in an unsteady DNS are reported in the grid convergence study in Appendix A. Furthermore, a comparison of the numerical computation with experimental measurement is performed. Hot-film measurements have been conducted in the laminar water channel at the Institute of Aerodynamics and Gas Dynamics of the University of Stuttgart. The channel provides a reproducible measurement environment for flat-plate laminar boundary-layer studies. The turbulence intensity is 0.05 % between 0.1 and 10 Hz (Wiegand Reference Wiegand1996; Puckert, Wu & Rist Reference Puckert, Wu and Rist2020). A Dantec 55R15 hot-film probe is connected to a Dantec Streamware bridge, which works according to the constant temperature anemometer principle. The analogue output voltage of the bridge is converted with a 16-bit National Instruments USB-6216 A/D device into a digital signal and then converted into velocity $u$ through King's law. A digital filter between 0.1 and 10 Hz is applied. Each discrete measurement has a measurement time of 180 s with a sampling rate of 100 Hz. For the present investigations, rotating roughness elements for the laminar water channel have been developed and used.

The hot-film probe is placed at the height of cylinder $k$ behind the cylinder at spanwise position $z = 2$ for both near- and far-wake measurements. The power spectral density (PSD) of disturbance root-mean-square values for cases $\varOmega _u= 0.75$ and $\varOmega _u =1.5$ are shown in figure 2 for four streamwise positions. Very good agreement of the PSD distribution is obtained for case $\varOmega _u =1.5$ with higher rotation speed and higher PSD amplitudes, while discrepancies appear for $\varOmega _u =0.75$. This can be explained by the lower induced kinetic energy at lower rotation speed of the cylinder which then remains close to the background eigen-disturbances of the experimental facility. It can be said that the PSD levels between experimental and numerical data in case $\varOmega _u=1.5$ match perfectly, despite some additional low-frequency signals ($\omega =0.2\unicode{x2013}0.3$) which stem from the water-channel background disturbances in case $\varOmega _u=0.75$. Consequently, the numerical simulation provides reliable results.

Figure 2. The PSD of velocity fluctuation $u'$ downstream of rotating cylinder for cases (a) $\varOmega _u = 0.75$ and (b) $\varOmega _u = 1.5$. Amplitude of experimental data is scaled for comparison with numerical data.

2.3. Dynamic mode decomposition

The DMD is a data-driven modal decomposition method with each identified mode having a single characteristic frequency (Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). It was first introduced by Schmid (Reference Schmid2010) to extract the coherent features of fluid flow from a DNS or experimental data. The Navier–Stokes operator is by nature nonlinear, but since DMD is demonstrated to be closely related to the Koopman operator (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009), it is capable of providing information about the dynamics of linear as well as nonlinear systems from snapshots of data, i.e. $q_i$. Here, the subscript denotes the $i$th time snapshot. In this decomposition method, it is assumed that two consecutive snapshots can be approximated by a linear constant projector $\boldsymbol {A}$ :

(2.6)\begin{equation} \boldsymbol{q}_{i+1} = \boldsymbol{Aq}_{i}. \end{equation}

Then, the time series of snapshots forms the following Krylov sequence:

(2.7)\begin{equation} \boldsymbol{Q}^N_1=\{\boldsymbol{q}_1,\boldsymbol{q}_2,\ldots,\boldsymbol{q}_N\} = \{\boldsymbol{q}_1,\boldsymbol{Aq}_1,\ldots,\boldsymbol{A}^{N-1}\boldsymbol{q}_1\}, \end{equation}

where the sub- and super-scripts of $\boldsymbol {Q}$ indicate the starting and ending snapshots. The dynamics of the underlying systems is thus described by the linear projector $\boldsymbol {A}$. The procedure of DMD analysis is then to determine the eigenvalues and corresponding eigenfunction of this projector, which can be approximated by the eigenvalues and eigenmodes of the following companion matrix:

(2.8)\begin{equation} \tilde{\boldsymbol{A}}=\boldsymbol{U}^*\boldsymbol{Q}_2^N\boldsymbol{V}\boldsymbol{\varSigma}^{{-}1}. \end{equation}

Here, the right-hand side matrices/vectors are from singular value decomposition $\boldsymbol {Q}^N_1=\boldsymbol {U} \boldsymbol {\varSigma } \boldsymbol {V}^*$. The eigenvalues of the companion matrix $\tilde {\boldsymbol {A}}$, i.e. $\tilde {\boldsymbol {A}} \boldsymbol {y}=\lambda \boldsymbol {y}$, are then taken as the DMD eigenvalues. Based on a more efficient numerical procedure by Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2013), the DMD mode corresponding to $\lambda$ is recovered from

(2.9)\begin{equation} \boldsymbol{\varphi}=\frac{1}{\lambda} \boldsymbol{Q}_2^N \boldsymbol{V} \boldsymbol{\varSigma}^{{-}1} \boldsymbol{y}. \end{equation}

As a last step, the complex frequency is computed from $\omega =\log (\lambda )/\Delta t$. The real part of the complex frequency $\omega _r$ corresponds to the angular frequency of the DMD mode, while the imaginary part $\omega _i$ determines its temporal growth rates. Then, the time series of snapshots $Q^N_1$ can be decomposed into the following general matrix form:

(2.10)\begin{equation} \boldsymbol{Q}^N_1 =\underbrace{[\boldsymbol{\varphi}_1,\boldsymbol{\varphi}_2,\ldots,\boldsymbol{\varphi}_N]}_{\boldsymbol{\varPhi}} \underbrace{ \begin{bmatrix} a_1 & & & \\ & a_2 & & \\ & & \ddots & \\ & & & a_N \end{bmatrix} }_{\boldsymbol{D}} \underbrace{ \begin{bmatrix} 1 & \lambda_1 & \cdots & \lambda^{N-1}_1 \\ 1 & \lambda_2 & \cdots & \lambda^{N-1}_2 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & \lambda_N & \cdots & \lambda^{N-1}_N \end{bmatrix} }_{\boldsymbol{\varLambda}}, \end{equation}

where $\boldsymbol {\varPhi }$ is the eigenfunction matrix and $\boldsymbol {\varLambda }$ is a Vandermonde matrix composed of the eigenvalues $\lambda _i$. DMD extracts a reduced-order representation of the linear projector $\boldsymbol {A}$, which maps a snapshot from the time series onto its successive snapshot. The diagonal amplitude matrix $\boldsymbol {D}$ determines the relative contribution of each mode to the particular realisation of the system represented by the analysed snapshots, which is calculated by a bi-orthogonal projection of the snapshots onto the DMD modes. More details can be found in Rowley et al. (Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009), Schmid (Reference Schmid2010) and Belson, Tu & Rowley (Reference Belson, Tu and Rowley2014), for instance.

The above numerical approach is provided from the parallelised python library modred of Belson et al. (Reference Belson, Tu and Rowley2014), which is object oriented and parallelised with mpi4py (Dalcín, Paz & Storti Reference Dalcín, Paz and Storti2005). Validation has been performed to confirm that DMD can recover the TS instabilities accurately (including frequency, mode shape and growth rate), see Wu & Rist (Reference Wu and Rist2020).

2.4. PKE analysis

A physical understanding of the instability mechanisms can be obtained from a PKE analysis. The basic idea is to calculate the energy transfer between perturbations and mean flow. The rate of change of PKE per unit mass, i.e. $e_k = 0.5 (\overline {u'^2} + \overline {v'^2} + \overline {w'^2})$, is governed by the following non-dimensional perturbation kinetic energy transport equation:

(2.11)\begin{equation} \frac{D e_k}{D t}={-} \underbrace{\overline{u'_{i}u'_{j}} \frac{\partial {\overline{u_i}}}{\partial x_{j}} }_{ \mathcal{P} } + \underbrace{ \frac{1}{2Re_k} \frac{\partial \overline{u'_i u'_i }}{ \partial x_j^2} }_{\mathcal{D}_\nu } - \underbrace{ \frac{1}{2} \frac{\partial \overline{u'_ju'_ju'_i}}{ \partial x_i } }_{\mathcal{T}_t} - \underbrace{ \frac{1}{Re_k} \overline{ \frac{\partial u'_{i}}{\partial x_{j}}\frac{\partial u'_{i}}{\partial x_{j}} } }_{ \mathcal{\varepsilon}_k} - \underbrace{ \frac{\partial \overline{u'_i p' }}{ \partial x_i} }_{\varPi_k}. \end{equation}

Here, $D/ D t = \partial /\partial t + \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla }$ is the material derivative, and the right-hand terms are: production ($\mathcal {P}$), viscous diffusion ($\mathcal {D}_\nu$), turbulent transport ($\mathcal {T}_t$), viscous dissipation ($\mathcal {\varepsilon }_k$) and pressure diffusion ($\varPi _k$). Here, the symbols prime ($'$) and overbar ($\,\bar{}\,$) denote perturbed and mean states, respectively. It is found that the production terms $\mathcal {P}$ are of particular interest. They represent the work of the Reynolds stress tensor against the mean-flow shear components, which are

(2.12)\begin{equation} \left.\begin{array}{c@{}} I_1 ={-}u'^2\dfrac{\partial \bar{u}}{\partial x}, \quad I_2={-}u'v'\dfrac{\partial \bar{u}}{\partial y}, \quad I_3 ={-}u'w'\dfrac{\partial \bar{u}}{\partial z}, \\[6pt] I_4 ={-}v'u'\dfrac{\partial \bar{v}}{\partial x}, \quad I_5={-}v'^2\dfrac{\partial \bar{v}}{\partial y}, \quad I_6 ={-}v'w'\dfrac{\partial \bar{v}}{\partial z},\\[6pt] I_7 ={-}w'u'\dfrac{\partial \bar{w}}{\partial x}, \quad I_8={-}w'v'\dfrac{\partial \bar{w}}{\partial y}, \quad I_9 ={-}w'^2\dfrac{\partial \bar{w}}{\partial z}. \end{array}\right\} \end{equation}

The sign of $I_i$ indicates a local production of PKE, with a positive sign denoting a destabilising and a negative sign a stabilising effect. The same analyses have been performed by Schmidt & Rist (Reference Schmidt and Rist2014). Similar analyses with the Reynolds–Orr equation, where the nonlinear terms have been dropped from the PKE transport equation, have been performed in the framework of local linear stability analysis (Cossu & Brandt Reference Cossu and Brandt2004; Chu et al. Reference Chu, Wu, Rist and Weigand2020; Wu et al. Reference Wu, Axtmann and Rist2021) and global linear stability analysis (Loiseau et al. Reference Loiseau, Robinet, Cherubini and Leriche2014).

4.1. Bifurcation diagram

In order to investigate the dynamic behaviour as the rotation speed $\varOmega _u$ of the cylinder increases, the perturbation velocities ($u', v', w'$) at downstream locations are recorded. As the DIV rotates with increasing cylinder rotation rate $\varOmega _u$, the perturbation's maximum in the $y$–$z$ plane changes its spatial location accordingly. Therefore, in figure 8 the maximum streamwise velocity perturbation $u'$ at $x=10$ is presented. Three distinctive dynamical regions are identified: laminar ($0 < \varOmega _u \leq 0.65$), transitional ($0.65 \leq \varOmega _u \leq 1.31$) and chaotic ($\varOmega _u \geq 1.375$). The corresponding flow features at different stages of transition are visualised by instantaneous contours of streamwise velocity $u$, as shown in figure 9. It is found that, below $\varOmega _u = 0.65$, the flow is stable and remains laminar despite a certain level of initial perturbations. The laminar flow is featured by stable velocity streaks with a high-speed streak in the centre, see figure 9(a). Beginning from the critical rotation speed $\varOmega _u^c = 0.65$, a self-sustained global instability is observed. The grey line in figure 8 is the best least-squares fit using $\sqrt {\varOmega _u - \varOmega _u^c}$ as the fitting function. This means that an increase of the cylinder rotation speed results in a supercritical Hopf bifurcation. The equilibrium state in the low $\varOmega _u$ region is broken and turns into a periodic limit-cycle oscillation at $0.65 \leq \varOmega _u \leq 1.31$. The transitional flow is characterised by braid-like structures between the low- and high-speed steaks, resembling the K–H instability, see figure 9(b). Interestingly, the boundary layer relaminarises in the range of $1.06 \leq \varOmega _u \leq 1.25$. Although the streamwise velocity perturbation at $x=10$ for $\varOmega _u=1.0$ is still larger than that for $\varOmega _u=0.75$, the overall relaminarisation effect is already distinctive, leaving the braid-like structure confined to the near-wake region, see figure 9(c). On further increase the rotation speed to $\varOmega _u=1.25$ in figure 9(d), the perturbation structure near the cylinders recedes to an almost imperceptible level. However, the perturbation amplitude for case $\varOmega _u >1.25$ climbs rapidly to follow the least-squares fitted level. The near-wake flow structure of case $\varOmega _u= 1.31$ in figure 9(e) is similar to that of case $\varOmega _u= 1.0$. With a further increase of the rotation speed to $\varOmega _u= 1.375$, the flow becomes chaotic. In figure 9( f), the regular braid-like patterns in both the near wake and far wake are broken, indicating that laminar–turbulent transition has already happened.

Figure 8. Bifurcation diagram of the global instability for $Re_{kk} =465.8$. Perturbation data obtained at $(x,y,z) = (10, 1, 2)$. Solid grey line represents best least-squares fit. Dashed line represents unstable equilibrium. The dash-dotted and dotted lines for $\varOmega _u \geq 1.0$ are unrelated to the dynamics; their purpose is to connect data.

Figure 9. Instantaneous contours of streamwise velocity $u$ at $y = 0.99$: (a) $\varOmega _u=0.25$, (b) $\varOmega _u=0.75$, (c) $\varOmega _u=1.0$, (d) $\varOmega _u=1.25$, (e) $\varOmega _u=1.31$, ( f) $\varOmega _u=1.375$. Yellow curves in the near wake ($x\leq 3$) indicate reverse-flow regions.

Figure 10 shows phase portraits of the perturbation velocities at probe location $(x,y,z) = (20,1,2)$, and the corresponding PSD of the streamwise perturbation velocity $u'$ for several characteristic $\varOmega _u$, as already used in figure 9. The cases $\varOmega _u = 0.25$ and $\varOmega _u = 1.25$ are illustrated in figures 10(a,b) and 10(g,h), where a fixed point in the phase portrait indicates that the flow is in a steady state. Although the perturbations in both cases are not absolutely zero, the dynamical energies as exhibited in the PSD are negligible. For $\varOmega _u = 0.75$, an ‘8’-shaped limit cycle in the phase portrait and a fundamental frequency ($\omega _1 \approx 0.9$) accompanied by its higher harmonics in the PSD plot indicate that the system undergoes a periodic limit-cycle oscillation. The rapid rise of perturbation $u'$ at $1.26< \varOmega _u<1.285$ in figure 8 is the result of a period-halving bifurcation and an amplification of the higher-harmonic components, as displayed in figures 10(i,j) and 10(k,l). It will be discussed in § 4.2 that the primary frequency $\omega _1$ is associated with hairpin vortex shedding induced by the rotating cylinders and the higher harmonics are an indication of nonlinear effects. At rotation speed $\varOmega _u = 1.31$, a second frequency $\omega _2 = 0.47$ begins to emerge in figure 10(n), albeit of low PSD magnitude. The incommensurate frequencies ($\omega _2 / \omega _1$ is irrational) lead to a quasi-periodic motion of the perturbations, turning the limit cycle into a limit torus in figure 10(m). Further increasing the rotation speed to $\varOmega _u = 1.375$, a cascade of period-halving and period-doubling bifurcations led to a chaotic attractor of the phase portrait in figure 10(o). Its PSD spectrum is dominated by the second basic frequency $\omega _2 = 0.465$.

Figure 10. Three-dimensional phase portraits of perturbations ($u', v', w'$) and their corresponding PSD at probe location $(x,y,z)=(20, 1, 2)$; (a,b) $\varOmega _u=0.25$, (c,d) $\varOmega _u=0.75$, (e,f) $\varOmega _u=1.0$, (g,h) $\varOmega _u=1.25$, (i,j) $\varOmega _u=1.2673$, (k,l) $\varOmega _u=1.285$, (m,n) $\varOmega _u=1.31$, (o,p) $\varOmega _u=1.375$. Colours in phase portraits vary with sampling time.

The dynamic behaviour of the perturbations as they evolve spatially is presented in figure 11. The transitional case $\varOmega _u=1.31$ is compared with the chaotic case $\varOmega _u=1.375$. For $\varOmega _u=1.31$, its closed loop phase portrait in figure 11(a) and the single dominant frequency $\omega _1=1.81$ at $x=1$ in figure 11(b) indicate that the perturbations undergo a period-1 limit-cycle oscillation. As the flow evolves downstream, the local perturbations undergo a period-halving bifurcation at $x=10$, resulting in the amplification of higher harmonics $2\omega _1$ and $3\omega _1$. Further downstream at $x=40$, a quasi-periodic oscillation of the perturbation is identified by the appearance of a second and a third frequency $\omega _2 = 0.465$ and $\omega _3 = 2.28$, respectively, in figure 11(j) and the limit torus in the phase portrait of figure 11(i). At $x=90$, the dominating frequency of the perturbations transits to the secondary mode $\omega _2$ in figure 11(n). The incommensurate frequencies turn the phase portrait into a non-overlapping limit torus in figure 11(m), while the primary frequency $\omega _1$ declines to an almost unnoticeable level, indicating a changeover of the flow feature from the near-wake vortex shedding to a low-frequency secondary instability in the far wake. Whereas for case $\varOmega _u=1.375$, the strange attractor at $x=1$ in figure 11(c) implies that the system experiences a chaotic motion in the very near wake, nevertheless, its PSD in figure 11(d) is still dominated by the primary frequency $\omega _1= 1.81$, which is related to vortex shedding, along with low energy harmonic $2\omega _1$ and low frequency component modulations. As the flow evolves to $x=10$, a cascade of period-doubling and period-halving bifurcations of the local perturbations are observed, which bring about the energy amplification of higher harmonic $2\omega _1$, a second frequency $\omega _2$ and a third frequency $\omega _3$, as shown in figure 11(h). In addition to its fractal structure in the phase portrait which is still observable in figure 11(g), the orbit is blended by irregular loops, demonstrating the chaotic nature of the perturbations. At $x=40$, a tertiary instability identified by $\omega _3$ dominates the PSD spectrum in figure 11(l), whereas the primary frequency $\omega _1$ diminishes. Further downstream, the discrete frequency spectrum converges rapidly to a continuous one in figure 11(p), indicating an evolution of the aforementioned quasi-periodic deterministic perturbations into a turbulent state (Borodulin & Kachanov Reference Borodulin and Kachanov2013).

Figure 11. Three-dimensional phase portraits of perturbations ($u', v', w'$) and their corresponding PSD at different probe locations: 1st row, $(x,y,z)=(1, 1, 2)$; 2nd row, $(x,y,z)=(10, 1, 2)$; 3rd row, $(x,y,z)=(40, 1, 2)$; 4th row, $(x,y,z)=(90, 1, 2)$. Left columns (a,b,e,f,i,j,m,n) $\varOmega _u=1.31$; right columns (c,d,g,h,k,l,o,p) $\varOmega _u=1.375$. Colours in phase portraits vary with sampling time.

Typically, laminar–turbulent transition is accompanied by an increasing drag coefficient. Figure 12 shows the streamwise evolution of the skin-friction coefficient $C_f = 2 \tau _{wall} / (\rho u^2_\infty )$ for various rotation speeds $\varOmega _u$. The three dynamical regions as distinguished in figure 8 are marked by contrasting colours. For laminar cases ($\varOmega _u = 0, 0.25, 0.46$), $C_f$ rises behind the roughness elements and decays gradually after reaching a mild peak at around $x= 20 \sim 30$. The transitional cases ($\varOmega _u = 0.65 \sim 1.25$) are characterised by a two-peak wavy progression of $C_f$, with the first peak at around $x\approx 5$ and the second peak at around $x\approx 20$. Starting from the near wake, the skin-friction coefficients of chaotic cases ($\varOmega _u \geq 1.375$) rise monotonically towards the value of a turbulent boundary-layer correlation. Here, the $1/7$ power law with experimental calibration (Schlichting & Gersten Reference Schlichting and Gersten2003) is given for reference. The different qualitative features of the skin-friction coefficient $C_f$ evolution substantiate the aforementioned categorisation of the boundary-layer flow, and are a consequence of symmetry breakings in the evolution of the vortical systems induced by the rotating cylindrical roughness elements.

Figure 12. Streamwise evolution of skin-friction coefficient $C_f$ compared with $1/7$ power law with experimental calibration (Schlichting & Gersten Reference Schlichting and Gersten2003). Colour of the plots are kept in accordance with figure 8: blue (laminar), red (transitional), black (chaotic). Note that the $x$-axis is a combination of linear ($x \leq 10$) and logarithmic regions ($x>10$).

4.2. Instantaneous flow-field structures

Having discussed the dynamic behaviour of the perturbations, the instantaneous flow structures as visualised by isosurfaces of $\lambda _2$ (Jeong & Hussain Reference Jeong and Hussain1995) will be analysed in this section. Figure 13 shows instantaneous snapshots of vortices emanating from the cylinders at a dimensionless rotation speed of $\varOmega _u=0.75$. The depicted frames denote the positions where the near-wake regions are extracted, and are shown in figure 14. For illustration purposes, an additional half-unit of the flow domain is added at $z>6$. The vortices obtained with SFD, which remain laminar, are included in blue colour for reference. In the very near wake ($x\leq 5$), the instantaneous vortices almost overlap with the laminar base-flow vortices. The flow is characterised by a DIV surrounded by a SIV, as shown in the enlarged figure 14. Unsteady behaviour appears first on the SIV at around $x\approx 5$, it develops into a tiny hairpin vortex with its head rotated around the DIV, which decays quickly at around $x=23$. Wavy structures are observed on the DIV from $x\approx 5$ as well, and the DIV meanders in the direction normal to the streamwise direction. Such a meandering dynamics is well reported in the initial stage of hairpin vortex generation (Schoppa & Hussain Reference Schoppa and Hussain2002; Wang, Huang & Xu Reference Wang, Huang and Xu2015). Slightly downstream at $x \approx 23$, an arch-shaped structure forms, which is tilted towards the high-speed streak. Its stronger leg is rooted in the DIV and its weaker leg stretches towards the low-speed streak. The slight cross-flow $w$ induced by the DIV, as shown in figure 7(d), pushes the tilted hairpin vortex packets away from the centre high-speed region. Thus, they merge at the spanwise boundary with the neighbour packets at $x\approx 50$, see figure 13. Further downstream, higher harmonics are triggered and quasi-streamwise vortices are generated. Another vortex combines at the periodic boundary and then evolves into a new packet of hairpin vortices at $x\approx 72$. Since the spanwise velocity $w$ is cancelled out by the neighbouring vortices here, this new hairpin vortex is symmetric with legs of identical strength. The primary frequency $\omega _1$, as plotted in figure 10(d), therefore, corresponds to the periodic shedding of the tilted hairpin vortices. Moreover, it turns out that the frequency and wavelength of the near-wake tilted hairpin vortex packet are identical to those of the far-wake symmetric hairpin vortex. In addition, the ejection ($Q2$, $u'<0$, $v'>0$) and sweep ($Q4$, $u'>0$, $v'<0$) events of a hairpin vortex are shown as well (Adrian Reference Adrian2007). The strength of both $Q2$ and $Q4$ events grows in the streamwise direction, proving that the hairpin vortex generation for this case is a streak-instability-based scenario (Schoppa & Hussain Reference Schoppa and Hussain2002). As shown in figure 14(a), the interaction of DIV and SIV in the near wake ($x \leq 25$) resembles the dynamics of a vortex pair (Leweke, Le Dizes & Williamson Reference Leweke, Le Dizes and Williamson2016), where a short wave 3-D elliptic instability (Kerswell Reference Kerswell2002) is responsible for the breakdown into turbulent motion.

Figure 13. Instantaneous isosurface of flow pattern for case $\varOmega _u=0.75$ by means of $\lambda _2=-200$, coloured by streamwise velocity $u$. (a) Top view, ejection ($Q2$, coloured yellow) and sweep ($Q4$, coloured green) events are shown by $u'v' = -0.003$, (b) side view. Depicted frames are enlarged in figure 14.

Figure 14. Enlarged views of the frames in figure 13. (a) Perspective view, (b) top view and (c) side view. Red curve indicates the rotation direction.

The instantaneous flow organisation of case $\varOmega _u=1.31$, which is closer to a chaotic state, is shown in figure 15. For clarity, the near-wake region is also enlarged in figure 16. As discussed above, the hairpin vortices for case $\varOmega _u= 0.75$ shed from the quasi-steady DIV. Here, the hairpin vortex is generated directly on the decelerating side of the rotating cylindrical roughness elements and regenerates into a hairpin packet. The hairpin packet has an interlaced structure which clearly proves its parent–offspring origin (Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999), i.e. hairpin vortices are regenerated from the preceding parent hairpin vortex. The accompanying $Q2$ and $Q4$ events confirm their regeneration mechanism as well. The enlarged views in figure 16 clearly show the spatial organisation of the hairpin vortex packet. The mean reverse-flow zone, as shown by $\bar {u}_x=0$, has a structure protruding into the decelerated boundary layer at height $y=0.6 \sim 0.85$, see the enlarged side view in figure 16(c). This protrusion splits the unsteady vortex on the decelerated side into two parts, which then develop into streamwise oriented legs of the first hairpin vortex. Due to a K–H instability, the inflectional shear flow rolls up into a spanwise vortex which bridges the two legs and forms the hairpin head. The local ejection events near the head and legs of the first hairpin vortex then regenerate the next pair of legs, thereafter an offspring hairpin vortex is created. Further downstream, this hairpin packet gradually decays, whereas a low-frequency hairpin starting at $x\approx 68$ is formed with its two legs rising from neighbouring units, see figure 14(a). It is to be noted that these far-wake low-frequency hairpin vortices exhibit asymmetry, i.e. symmetry breaking, which is a precursor of chaos.

Figure 15. Instantaneous isosurface of flow pattern for case $\varOmega _u=1.31$ by means of $\lambda _2= -200$, coloured by streamwise velocity $u$. (a) Top view, ejection ($Q2$, coloured yellow) and sweep ($Q4$, coloured green) events are shown by $u'v' = -0.003$, (b) side view. The mean reverse-flow zone $\bar {u}_0=0$ is marked by cyan colour. Depicted frames are enlarged in figure 16.

Figure 16. Enlarged views of the frames in figure 15. (a) Perspective view, (b) top view and (c) side view. Red curve in (b) indicates the rotation direction.

The flow statistics of the chaotic cases ($\varOmega _u = 1.5, 2$) are presented together with mean-flow profiles in figure 17. The streamwise mean-flow velocity profiles, which are non-dimensionalised by $\bar {u}^+=\bar {u} / u_\tau$ with regards to wall coordinate $y^+= y u_\tau / \nu$, are given in figure 17(a,d). The theoretical curves for the viscous sublayer $\bar {u}^+ = y^+$ and the logarithmic layer $\bar {u}^+ = 1/\kappa \ln y^+ +C$ with $\kappa = 0.41$ and $C=5.2$ are plotted with dashed lines. For reference, the fully developed turbulent data of Schlatter & Örlü (Reference Schlatter and Örlü2010) are depicted with solid red lines as well. The mean-flow profiles of both cases evolve gradually towards a turbulent one, with case $\varOmega _u = 1.5$ in an oscillatory manner and case $\varOmega _u = 2$ in a gradual manner. The streamwise components of Reynolds normal stress $\tau ^+_{xx}$ as well as shear stress $\tau ^+_{xy}$ are given. At $x=1$, a very high peak of both Reynolds stresses is observed at height $y^+= 30$, which corresponds to the height of the cylinder. For the streamwise normal stress $\tau _{xx}$, this peak quickly decays and shifts to $y^+= 45$ at station $x= 10$. Further downstream at $x=20$, a second peak appears and grows in height $y^+\approx 13$. Whereas this location conforms to the peak of the fully developed turbulence, the magnitude grows to approximately two times larger at $x=100$. For the shear stress $\tau ^+_{xy}$, the peak at $y^+= 30$ gradually spreads out towards the fully developed turbulent profile.

Figure 17. Evolution of mean velocity profiles (a,d) and turbulence statistics (b,c,e,f) along streamwise direction. Upper row corresponds to case $\varOmega _u = 1.5$ and lower row to case $\varOmega _u = 2$. Red solid line is from DNS data of Schlatter & Örlü (Reference Schlatter and Örlü2010). For clarity, shear stresses of near-wake ($x \leq 20$) and far-wake ($x \geq 40$) locations are plotted with offset.

4.3. DMD analysis

In this section, the transitional flows controlled by the rotating cylindrical roughness elements are investigated by DMD (Schmid Reference Schmid2010; Belson et al. Reference Belson, Tu and Rowley2014). The spatio-temporal process of the 3-D boundary-layer flow can be decomposed into a set of coherent spatial modes which feature corresponding frequencies. At high rotation speed, nonlinear modal interactions can develop. In such a case, the resulting DMD mode should be interpreted as one particular representative phase-dependent dynamical mode determined by the nonlinear interactions of the snapshots employed. In low rotation speed cases, however, such nonlinear effects can be neglected. The obtained DMD modes are therefore identical to a global mode as obtained from a global stability analysis (Loiseau et al. Reference Loiseau, Robinet, Cherubini and Leriche2014). Since only the asymptotic behaviour of the perturbations are considered in the current work, DMD is performed with equilibrium state snapshots. Approximately 200 snapshots are used to perform DMD for the investigated cases, which results in a temporal sampling rate that is more than twice the highest dominant frequency of instability. Figure 18(a) shows the complex DMD eigenvalues for case $\varOmega _u= 1.31$. Here, almost all modes are located on the unit circle by virtue of the equilibrium flow state. The few low amplitude modes inside the unit circle are most likely nonlinearly generated higher harmonics and numerical artefacts. The non-dimensional physical frequency $\omega _r$ and its corresponding mode amplitude coefficients $a_i$ are plotted in figure 18(b). The primary frequency $\omega _1$ and its higher harmonics as well as the second frequency $\omega _2$ as shown in figure 11 are recovered by DMD.

Figure 18. The DMD spectrum of case $\varOmega _u = 1.31$. Colour and size in (a) are related to coefficient amplitude. (a) Complex eigenvalues $\lambda$ and (b) amplitudes of DMD modes.

Figure 19 shows the leading DMD modes of the selected cases from the bifurcation diagram in figure 8. The modes consist of positive and negative patches of velocity mostly localised around the central low-speed region, i.e. the DIV structure in the laminar cases. In case $\varOmega _u = 0.75$, the leading DMD mode is a low-frequency spatially growing mode, whose mode amplitude is shown in figure 25(b,c) in Appendix A to saturate at $x\approx 50$. For $\varOmega _u = 1.0$, the leading DMD mode is a high-frequency mode which fades away after $x=50$, whereas the DMD mode for the relaminarised case $\varOmega _u = 1.25$ only emerges in the far wake ($x\approx 50$). It is to be noted that the amplitude of the leading mode in $\varOmega _u = 1.25$ is several orders of magnitude lower than the DMD modes of other cases. Further increase of the rotation speed to $\varOmega _u = 1.31$ leads to the occurrence of two leading DMD modes as shown in figure 18, i.e. a high-frequency mode (coloured red/blue, $\omega _1=1.8$) in the near wake and a low-frequency mode (coloured white/black, $\omega _2=0.47$) in the far wake. The asymmetrical spatial structure of the low-frequency mode implies the onset of chaotic flow.

Figure 19. Real part of dominant DMD modes. Frequencies for each mode: (a) $\omega _r=0.92$, (b) $\omega _r=2.2$, (c) $\omega _r=0.66$, (d) primary mode $\omega _1=1.8$ with isosurface in red/blue and second mode $\omega _2=0.46$ with isosurface in white/black. Isosurfaces depict $\hat {u}= \pm 10\,\%$ (in colour red/blue or white/black) of maximum amplitude of mode.

The structures of the leading modes and their spatial locations relative to the mean flow are visualised in slices at $x=10$ in figure 20. Here, only the positive-$z$ half of the modes is shown due to the symmetry of the DMD modes with regard to the $z=0$ axis. The 3-D high shear region is identified by the $I_2$-criterion (Meyer Reference Meyer2003). For cases $\varOmega _u = 0.75, 1.0, 1.25$, the mean flow is characterised by a tilted high shear region surrounding the strong streamwise vortex, i.e. the DIV. Typical instability modes behind a static cylindrical roughness element can be categorised by their spanwise symmetries, i.e. varicose or sinuous, with respect to the vertical axis where the cylinder is located. As shown in figure 20, such spanwise symmetrical characteristic is no longer a feature of the induced modes. Nevertheless, the DMD modes presented here are exclusively located within the high shear stress regions, indicating their inflectional nature. Similarly tilted modes are observed with skewed (oblique) roughness elements in Groskopf & Kloker (Reference Groskopf and Kloker2016), with distributed surface roughness in Di Giovanni & Stemmer (Reference Di Giovanni and Stemmer2018) and more intensively with cross-flow instability in Malik, Li & Chang (Reference Malik, Li and Chang1994); Wassermann & Kloker (Reference Wassermann and Kloker2002) and Shi, Mader & Martins (Reference Shi, Mader and Martins2021). As for case $\varOmega _u = 1.31$, the mean flow has already been distorted by the interaction of its leading modes, which are identified in figure 18.

Figure 20. Real part of leading DMD modes with primary frequency $\omega _r$ in $x=10$ plane shown by red/blue contours for cases (a) $\varOmega _u= 0.75$, (b) $\varOmega _u= 1.0$, (c) $\varOmega _u= 1.25$, (d) $\varOmega _u= 1.31$. Red/blue contours show positive/negative values of normalised spatial eigenfunction. Thin dash-dotted lines are isolines of mean streamwise velocity $\bar {u} = 0.1{-}0.99$. Thick cyan solid lines visualise high shear stress regions by means of $I_2$-criterion (Meyer Reference Meyer2003).

4.4. Instability mechanisms

The PKE analysis is performed to get an insight into the instability mechanisms for the above identified DMD modes, particularly, how the modes obtain their energy from the work of Reynolds stress against the mean-flow shear. Figure 21 presents the streamwise evolution of perturbation productions $\int _{y,z} I_i \,\textrm {d}y\,\textrm {d}z$ scaled by the local dissipation $\int _{y,z} |D| \,\textrm {d}y\,\textrm {d}z$. The purpose is to magnify the otherwise imperceptible energy production at the initial stage of the instability. It is evident that the energy production is dominated by the $\int _{y,z} I_2 \,\textrm {d}y\,\textrm {d}z$ and $\int _{y,z} I_3 \,\textrm {d}y\,\textrm {d}z$ terms for almost all cases, which is reasonable and common for a streaky boundary layer, see Cossu & Brandt (Reference Cossu and Brandt2004), Loiseau et al. (Reference Loiseau, Robinet, Cherubini and Leriche2014) and Wu et al. (Reference Wu, Axtmann and Rist2021) as well. However, it should be noted that in the roughness region ($-0.5 \leq x \leq 0.5$) the term $\int _{y,z} I_1 \,\textrm {d}y\,\textrm {d}z$ dominates the energy production for cases with $\varOmega _u \geq 1.0$. As shown in (2.12), the term $I_1$ represents the production of $e_k$ from the normal Reynolds stress $u'^2$ against the streamwise gradient of the mean-flow streamwise component $\partial \bar {u} / \partial x$. Since the sign of the normal shear $u'^2$ is always positive and the sign for the gradient $\partial \bar {u} / \partial x$ on the reverse-flow side of the cylinder is negative due to the deceleration effect, the production term $I_1$ then turns out to be invariably positive. The mechanism is hereafter termed a deceleration mechanism. The crucial importance of the deceleration effect on intensifying the growth rate of the shear-layer instability has been addressed by Gad-El-Hak et al. (Reference Gad-El-Hak, Davis, Mcmurray and Orszag1984) and Shtern & Hussain (Reference Shtern and Hussain2003).

Figure 21. Streamwise evolution of production terms $\int _{y,z} I_i / |D| \,\textrm {d}y\,\textrm {d}z$ of leading DMD modes for (a) $\varOmega _u= 0.75$, (b) $\varOmega _u= 1.0$, (c) $\varOmega _u= 1.25$, (d) $\varOmega _u= 1.31$. Dissipation $D$ is plotted at $I_i/D = -1$ for reference. Note that $x$-axis contains linear ($x \leq 10$) and logarithmic ranges ($x>10$). Vertical shaded region in the insets indicate the location of cylinders.

The real parts of the leading modes at slice $x=0$ are shown in the first row of figure 22, the second and third rows represent the local total energy production at slices $x=0.3$ and $x=2$, respectively. For higher rotation speed cases $\varOmega _u = 1, 1.25, 1.31$, both the DMD modes and the local total energy productions $\sum I_i$ are localised around the protruding structures as marked by the reverse-flow region ($\overline {u}_x=0$), where the deceleration mechanism gives birth to the DMD modes shown in figure 19(b–d). Nevertheless, the total energy production for case $\varOmega _u = 0.75$ is counterbalanced by a positive region (contribution of $I_1$) and a negative region (contribution of $I_3$), which implies a different instability mechanism for the DMD mode as shown in figure 19(a).

Figure 22. Real part of leading DMD modes in $x=0$ plane shown by red and blue contours for case (a) $\varOmega _u= 0.75$, (b) $\varOmega _u= 1.0$, (c) $\varOmega _u= 1.25$, (d) $\varOmega _u= 1.31$. Thin dash-dotted lines are isolines of mean streamwise velocity $\bar {u} = 0.1-0.99$. Recirculating region is marked by black solid lines. Total local energy production $\sum I_i$ at $x=0.3$ (second row) and $x=2$ (third row) for (e,i) $\varOmega _u= 0.75$, (f,j) $\varOmega _u= 1.0$, (g,k) $\varOmega _u= 1.25$, (h,l) $\varOmega _u= 1.31$. Yellow dash-dot lines mark the phase speed $c_{ph}$ of the corresponding mode. Recirculating regions are marked by magenta thick solid lines. Thin magenta dash-dotted lines are isolines of $\bar {u} = 0.99$. Thick cyan solid lines visualise shear regions by means of $I_2$-criterion (Meyer Reference Meyer2003). Yellow line in (h) for comparison with figure 23(c).

Figure 23. Elliptic flow regions ($\beta <1$, coloured blue) at $x=2$ for cases (a) $\varOmega _u=0.75$ and (b) $\varOmega _u=1$. Centrifugally unstable regions ($\gamma <0$, coloured red) at $y=1.07$ for case (c) $\varOmega _u=1.31$. The dashed circle in (c) marks the location of the cylinder. Yellow line in (c) for comparison with figure 22(h).

The third row of figure 22 shows the total energy productions at slice $x=2$, which is a position across the downstream reverse-flow region. For case $\varOmega _u = 1.31$ in figure 22(l), the flow has curved streamlines around the cylinder and the velocity is decelerated at the region of total energy production, which are flow conditions prone to a centrifugal instability (Floryan Reference Floryan2002; Lanzerstorfer & Kuhlmann Reference Lanzerstorfer and Kuhlmann2012). The argument is supported by the sufficient condition of Sipp & Jacquin (Reference Sipp and Jacquin2000), that a 2-D inviscid flow is centrifugally unstable if the following parameter is negative along a closed streamline $\psi = {\rm const}$.

(4.1)\begin{equation} \gamma := \frac{|\bar{\boldsymbol{u}}|\omega }{R} < 0, \end{equation}

where $\omega$ is the vorticity of the mean flow, and $R$ is the local algebraic radius of curvature, defined as

(4.2)\begin{equation} R=\frac{|\bar{\boldsymbol{u}}|^3 }{ \boldsymbol{\nabla} \psi \boldsymbol{\cdot} (\bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{u}}) }. \end{equation}

Figure 23(c) compares the regions where $\gamma <0$ at height $y=1.07$ for case $\varOmega _u=1.31$ with the mean flow. It can be observed that negative $\gamma$ occurs in regions of flow divergence and convergence in the tilted wake, exactly at the edges of the vortex train in figure 15(b). The total energy production in figure 22(h) is also observed to overlap with this region. For better comparison, reference spatial locations are marked by a short yellow line in figures 22(h) and 23(c).

The cases with $\varOmega _u = 0.75$ and $1.0$ are slightly different because the perturbation kinetic energy is produced in two regions: the high shear region as identified by the $I_2$-criterion and the centre of the elliptic vortex of the mean flow, see figure 22(i,j). The latter is an indication of an elliptic instability (Pierrehumbert Reference Pierrehumbert1986; Sipp & Jacquin Reference Sipp and Jacquin1998). This interpretation is supported by the theory of elliptic instability proposed by Pierrehumbert (Reference Pierrehumbert1986) and Waleffe (Reference Waleffe1990), that a strained vortex is subject to elliptic instability if the eccentricity parameter $\beta$ of the flow is less than 1, i.e.

(4.3)\begin{equation} \beta := \frac{2\epsilon}{ |\omega| } < 1. \end{equation}

Here, $\epsilon$ is the magnitude of strain and $\omega$ is the vorticity. Figure 23(a,b) shows the elliptic flow regions at $x=2$ for cases $\varOmega _u=0.75$ and $\varOmega _u=1$, respectively. Further, the maximum of the perturbation, which is located in the crescent-shaped high shear region (see figure 20a,b for instance) and the maximum of perturbation energy production $\sum I_i$ approximately coincide with the location of the critical layers as marked by the phase velocities $c_{ph}$. A criterion that meets the Rayleigh–Fjortøft necessary condition for inviscid inflectional instability (Schmid & Henningson Reference Schmid and Henningson2001). Therefore, the cases $\varOmega _u = 0.75$ and $1.0$ are subject to a combination of inviscid instability and elliptic instability in the near wake. The case for $\varOmega _u=1.25$ in figure 22(k) exhibits the same inviscid instability feature. However, the case is unique in that it is dominated by a relaminarisation effect which is not the topic of the current investigation. In the far wake, both centrifugal and elliptic effects fade away.