5.1.1. Single bluff-body wake

For $g^*<0.25$, where the base flow is in SS, the results from Floquet analysis show that the onset of three-dimensionality develops at extremely low ${\textit {Re}}$, which ranges between 30 %–50 % of $\{{\textit {Re}}_{cr-1}\}_S \approx 188.84$ (figure 5). The small ${\textit {Re}}_{cr-1}$ can be attributed to the weaker gap flow at small $g^*$, causing the cylinder group to resemble a single bluff object with an effective length scale of approximately $2D+g^*{D}$ rather than $D$ of its single-cylinder counterpart. The single bluff object introduces a wider wake (Roshko Reference Roshko1954) and, hence, a smaller ${\textit {Re}}_{cr-1}$ related to the onset of 3-D transition. This also explains the large $\lambda _{cr-1}$ of SS base flow shown in figure 6.

The spatial distribution of the $\hat {\omega }_x$ field (figure 7a) for the unstable mode (with real and positive $\mu$) is similar to the mode A of a single cylinder, where high concentrations of $\hat {\omega }_x$ are located close to the vortex cores. With increasing $g^*$, the gap flow and the outer shear layers are strengthened and weakened, respectively, resulting in a weaker single bluff-body wake generated by the two cylinders. The weakening process is first supported by the quantification of the strength of shed vortices with varying $g^*$ at ${\textit {Re}} = 80$. Figure 7(c) shows the instantaneous strengths of a pair of shed vortices, represented by the $|\omega _z|$ values at vortex cores (see $\varOmega _1$ and $\varOmega _2$ in figure 7a), at the phase where $C_{L-1}$ reaches local maxima. As $g^*$ increases, the strength of both negative ($\varOmega _1$) and positive ($\varOmega _2$) shed vortices decrease, indicating a weakening of the outer shear layers. This weakening process is further evident in the location of the wake recirculation bubble, obtained from the time-averaged flow fields over 20 vortex shedding periods. In the presence of gap flow the wake recirculation bubble detaches from the cylinder surfaces and forms two wake stagnation points, namely $x_{s1}$ and $x_{s2}$, as shown in figure 7(d) at $(g^*, Re) = (0.2, 80)$. The bubble moves downstream with increasing $g^*$ until it disappears when the spatio-temporal symmetry breaks (e.g. at $g^* = 0.3$ for ${\textit {Re}} = 80$). Due to the increased gap flow and the weakened outer shear layer with increasing the $g^*$, the two wake stagnation points are formed further downstream as displayed in figure 7(d). Therefore, we attribute the increase of ${\textit {Re}}_{cr-1}$ with $g^*$ for the SS base flow to the weakening process of the single bluff-body wake.

Figure 7. (a) The spatial distribution of $\hat {\omega }_x$ contours for a mode A-SS at $(g^*, Re, \beta ) = (0.1, 80, 0.8)$ overlapped on base flow $\omega _z$ isolines at the instant where the lift coefficient of the top cylinder $C_{L-1}$ reaches a local maximum (the solid dots in (b)). (c) Variations of $|\omega _z| - g^*$ at $Re = 80$, where the blue and red curves are quantified as the cores of negative vortex $\varOmega _1$ and positive vortex $\varOmega _2$ at the instant shown as the solid dots in (b). (d) Variations of the wake stagnation points, $x_{s1}$ and $x_{s2}$, with respect to $g^*$ at $Re = 80$.

As $g^*$ is increased to 0.25–0.6 in the ASS state, the wake transition to 3-D flow shows a distinct characteristic, as demonstrated in figure 8 at $g^* = 0.3$. The $|\mu | - \beta$ and $\lambda - {\textit {Re}}$ variations (figure 8a,b) suggest that the most unstable 3-D mode (real and positive $\mu$) occurs at ${\textit {Re}}_{cr-1} \approx 87.12$ with $\lambda _{cr-1} \approx 8.64$. The spatial and temporal distributions of the reconstructed flow field from Floquet analysis in figure 8(c,d) indicate that the $\hat {\omega }_x$ is slightly deflected upwards due to the downwards biased gap flow. High concentrations of $\hat {\omega }_x$ are located close to both vortex cores and braid shear layers. The $\hat {\omega }_x$ field in figure 8(c) differs from that of the typical mode A (figure 3b-i) in a sense that the $\hat {\omega }_x$ values in neighbouring braid shear layers are of the same signs and weakly holds an even RT symmetry along $y_c/D \approx 0.5$. The ${\textit {Re}}_{cr-1}$ decreases from around 102 to 60 as $g^*$ increases from 0.25 to 0.6 (figure 5), forming an opposite trend to that of the SS flow. The $\lambda _{cr-1}$ of ASS flow amounts to more than two times the $\{\lambda _{cr-1}\}_S \sim 4D$ due to the single bluff-body wake of ASS flow. The $\lambda _{cr-1}$ increases with $g^*$ and shows a rapid increase at $g^* = 0.4\unicode{x2013} 0.5$. The above variations of ${\textit {Re}}_{cr-1}$ and $\lambda _{cr-1}$ with $g^*$ suggest that the first 3-D mode of the ASS flow may be associated with the unique characteristics that exist in the base flow.

Figure 8. (a) Evolutions of the modulus of the Floquet multipliers $|\mu |$ varying with the spanwise wavenumber $\beta$ for the ASS flow at $g^* = 0.3$, where $\mu$ values with complex conjugate, real and positive and real and negative are marked as grey, blue and red colours, respectively. (b) Corresponding neutral stability curves plotted in the ${\textit {Re}}$ and wavelength ($\lambda = 2{\rm \pi} /\beta$) space. (c) Spatial distribution of $\hat {\omega }_x$ field for a mode A$'$-ASS at $(g^*, Re, \beta )= (0.3, 100, 0.8)$. (d) The temporal evolution of $\hat {\omega }_x$ on the line of $x/D = 12$ over four vortex shedding periods ($T_v$).

The above distinct features for the ASS flow (e.g. decreasing trend of ${\textit {Re}}_{cr-1}$ with $g^*$) are attributed to the disturbances of biased gap flow to the cluster shear layers around the cylinders. It is well known that the wake transition to 3-D flow is subcritical for an isolated cylinder (e.g. Williamson Reference Williamson1996b; Thompson et al. Reference Thompson, Hourigan, Ryan and Sheard2006) and the critical Reynolds number for the subcritical transition is dependent on the level of disturbances in the flow. Our DNS results (presented in § 5.3.3) confirm that the transition of 2-D ASS flow to 3-D flow is also subcritical. The disturbance level induced by the biased gap flow is likely dependent on the level of asymmetry of the base flow and the flow rate through the gap. To confirm our inference, we quantify the variations of flow asymmetry and flow rate with $g^*$ over $g^* = 0.2\unicode{x2013} 0.5$ in figure 9. The level of asymmetry of the base flow is quantified by

(5.1)\begin{equation} \varGamma=\sum_{i=1}^N \left|\iint_{\varOmega}\omega_z(t_i)\,\mathrm{d}\varOmega+ \iint_{\varOmega}\omega_z(t_{i+N/2})\,\mathrm{d}\varOmega\right|, \end{equation}

where $N = 32$ is the equi-spaced snapshots over $1T_v$; $\varOmega$ is the integration region at $[-1\leq x/D \leq 10,\ -5 \leq y/D \leq 5]$. According to (4.1), the terms of $\iint _{\varOmega }\omega _z(t_i)\,\mathrm {d}\varOmega \ne -\iint _{\varOmega }\omega _z(t_{i+N/2})\,\mathrm {d}\varOmega$ when the spatio-temporal symmetry is broken. The magnitude of $\varGamma$ effectively measures the degree of the wake asymmetry during one vortex shedding cycle. The flow rate through the gap ($x/D = 0$) is defined as

(5.2)\begin{equation} \overline{U_{gap}}=\int_{{-}0.5g^*}^{0.5g^*} \overline{u(y)}\,\mathrm{d}y, \end{equation}

where the $\overline {u(y)}$ is obtained from an averaged field over 100 non-dimensional time units. The results presented in figure 9(a) show that both $\varGamma$ and $\overline {U_{gap}}$ increase with increasing $g^*$, confirming our inference. The strengthened gap flow rate and the large asymmetry of the base flow provide a sustained disturbance to the flow, leading to the wake transition to 3-D flow at low ${\textit {Re}}$.

Figure 9. The variations of the extent of wake asymmetry $\varGamma$ (left $Y$ axis) and gap flow rate $\overline {U_{gap}}$ (right $Y$ axis) with respect to (a) $g^*$ and (b) ${\textit {Re}}$ values. The solid and half-filled symbols indicate the occurrence of first and secondary 3-D instabilities, respectively.

The first 3-D mode for the $\textrm {IP}_B$ base flow at $g^* = 0.6\unicode{x2013} 0.9$ forms the similar mode A$'$ structure, where the value of Floquet multiplier $\mu$ is real and positive. The large inclination of the upward biased gap flow initiates the mode A$'$ structure in the braid shear layers of the wake behind the bottom cylinder, as presented in the Floquet mode and linear reconstructed 3-D field in figures 10(a) and 10(b), respectively. In figures 5 and 6 we did not give the marginal curves of this mode. This is because, all the $\textrm {IP}_B$ base flows in the purple shaded region in figure 5 are unstable to the mode A$'$ instability, as confirmed by both Floquet analysis and 3-D DNS using $L_z = 28$. As the $\textrm {IP}_B$ flow below the purple shaded region in figure 5 cannot be simply generated by decreasing ${\textit {Re}}$ from a fully developed $\textrm {IP}_B$ flow, the exact ${\textit {Re}}_{cr-1} - g^*$ and $\lambda _{cr-1} - g^*$ curves for the $\textrm {IP}_B$ base flow are not provided. We anticipate that the ${\textit {Re}}_{cr-1}$ of mode A$'$ in $\textrm {IP}_B$ is likely lower than its dominant ${\textit {Re}}$ range in 2-D simulation and presumably follows the ${\textit {Re}}_{cr-1} - g^*$ trend of mode A$'$ in the ASS.

Figure 10. Same as figure 3 but for the mode A$'$-$\textrm {IP}_B$ instability at $(g^*, Re, \beta )= (0.8, 72, 0.6)$. (a) Spatial distribution. (b) Reconstructed 3-D field.

5.1.2. Flip-flop wake

For flows in the FF regime, we could not perform Floquet analysis as the wake is aperiodic over each vortex shedding cycle. Therefore, we employed a series of 3-D DNS at $g^* = 0.6\unicode{x2013} 2.0$ and ${\textit {Re}} = 55\unicode{x2013} 100$ to investigate its 3-D transition features. The 3-D DNS is conducted by using its fully developed 2-D counterpart as an IC; and the dimensionless length of the cylinder is selected as $L_z = 24$. The flow without and with 3-D instabilities are presented as hollow and filled circles in figure 5 in the region filled with FF flow. The development of 3-D flow is measured during the simulation through the evolutions of $\varepsilon _{xy}/\varepsilon _{t}$ in the near wake region of $[0 \leq x/D \leq 10,\ -5 \leq y/D \leq 5]$. In the present study the 3-D instability of the FF flow (mode FW) is identifiable when the $\varepsilon _{xy}/\varepsilon _{t} > 10^{-3}$; and we use this criteria to assess whether mode FW is developed.

The wake transition to 3-D flow for FF flow is rather interesting in figure 5. The ${\textit {Re}}_{cr-1}$ experiences a step jump from ${\sim }65$ to ${\sim }75$ when the base flow changes from ASS to FF at $g^* \sim 0.6$. Further increasing $g^*$ over the FF range, the ${\textit {Re}}_{cr-1}$ for mode FW decreases till about $g^* = 1.5$, followed by a slight increase for $g^* \sim 1.5\unicode{x2013} 1.8$. The discontinuity at $g^* \sim 0.6$ suggests that the onset of 3-D instability for the FF flow may have been triggered by different flow characteristics compared with that of the ASS flow. We attribute it to the distinct irregular wake evolutions that exist in FF flow, which lead to the development of mode FW at such low ${\textit {Re}}_{cr-1}$. The decreasing trend of ${\textit {Re}}_{cr-1}$ with increasing $g^*$ for $g^* \leq 1.5$ is possibly related to the increasing level of wake irregularity. These conjectures are supported by the characteristics of the lift coefficients on the top and bottom cylinders, i.e. $C_{L,1}$ and $C_{L,2}$, respectively, and will be elaborated in the following paragraphs.

As an example, figure 11 shows the development of mode FW from a 2-D FF flow at $g^* = 1.0$ and ${\textit {Re}} = 64\unicode{x2013} 100$. The force time histories in the blue shaded area (figure 11a–d) are those obtained from 2-D DNS; and the instantaneous flow field at the end of 2-D simulation (termed as $t_0$) is used for the subsequent 3-D DNS. Comparing the $C_{L,1}$ and $C_{L,2}$ time histories in 2-D and 3-D simulations, the lift coefficient evolutions in FF flows become much more complex as ${\textit {Re}}$ is increased. The development of mode FW does not change the synchronisation state of two wakes in the 2-D FF flow. In addition, the results in figure 11 demonstrate a strong correlation between the complexity of lift coefficients in the 2-D DNS (blue shaded area in figures 11a–d) and the development of the 3-D structure (flow fields in figures 11i–l). Specifically, the time history of $C_{L,1}$ and $C_{L,2}$ forms a transition from quasi-periodic evolutions to irregular evolutions at ${\textit {Re}}_{cr-1} \approx 66$, coinciding with the wake transition from 2-D FF flow to mode FW. To quantify this correlation, we calculate the fluctuation of the total lift coefficient peaks, $\varPsi$, which is defined as

(5.3)\begin{equation} \varPsi= RMS(C_{L,T}^{pks}-\overline{C_{L,T}^{pks}}), \end{equation}

where $C_{L,T}^{pks}$ represents the local maxima of $C_{L,T} = C_{L,1}+C_{L,2}$; $\overline {C_{L,T}^{pks}}$ is the mean value of $C_{L,T}^{pks}$. The $\varPsi$ value is calculated through 2-D results because we believe that the irregularity triggering the 3-D transition largely originates from the 2-D flow.

Figure 11. Time histories of the lift coefficients (a–d) and $\varepsilon _{xy}/\varepsilon _{t}$ (e–h) for FF flows at $g^*=1.0$ and ${\textit {Re}} = 64 -100$, where $t_0$ indicates the time instant where 3-D DNS starts. Instantaneous 3-D flows in the global view (i–l) and $x$–$z$ plane view (m–p) of each case, where the isosurfaces are $\omega _x = \pm 0.4$ (blue and red) and $\omega _z = \pm 0.5$ (grey).

The evolutions of $\varPsi$ with $g^*$ and ${\textit {Re}}$ are shown in figure 12(a) and 12(b), respectively, where flows that transition to mode FW are marked as solid symbols. Figure 12(a) suggests that at ${\textit {Re}} = 62$ the $\varPsi$ of FF flow increases with enlarging $g^*$ for $g^*\leq 1.7$, which agrees with our anticipation. Mode FW is observed when the $\varPsi$ manifests a large magnitude at $g^* = 1.3\unicode{x2013} 1.8$. The values of $\varPsi \approx 0$ at $g^* = 0.7$ and 2.0 are due to the flow transits to an IP state with periodic wakes behind two cylinders. The evolutions of $\varPsi$ with ${\textit {Re}}$ at different $g^*$ values (figure 12b) are not monotonic due to the complex flow state transitions. The initiations of mode FW with respect to ${\textit {Re}}$ are highly associated with the large magnitudes of $\varPsi$.

Figure 12. Evolutions of wake irregularity parameter, $\varPsi$, with (a) $g^*$ at ${\textit {Re}} = 62$ and (b) ${\textit {Re}}$ at $g^* = 1.0, 1.5, 2.0$. The $\varPsi$ values are calculated based on 2-D results. The symbols filled by solid colours are cases where we observed mode FW in the 3-D DNS. Cases circled in red at $g^* = 2.0$ only transit to mode FW by using $L_z = 48$.

The relationship between $\varPsi$ and the 3-D instability, except for $g^* = 2.0$, consistently shows that the 3-D instability occurs when $\varPsi > \sim 0.0245$ (the dashed grey line in figure 12). However, an exception is found at $g^* = 2.0$ where no identifiable 3-D instability is observed with $\varPsi = 0.376$ at ${\textit {Re}} = 56$. We suspect the discrepancy is caused by the use of $L_z = 24$ at large $g^*$, which may not be large enough for the mode FW to develop as the wavelength of the 3-D structure is expected to be large at low ${\textit {Re}}$ values. To confirm this inference, we conducted additional 3-D simulations for $g^* = 2.0$ using $L_z = 48$. As expected, the flows at ${\textit {Re}} = 54$ and 56, where $\varPsi > 0.0245$ in figure 12(b), transit into mode FW with weak 3-D structures ($\varepsilon _{xy}/\varepsilon _{t} \approx 10^{-3}$). No further investigations were conducted to seek the exact ${\textit {Re}}_{cr-1} - g^*$ boundaries. We anticipate the lower bound ${\textit {Re}}_{cr}$ at large $g^*$ may be slightly lower than that given in figure 5 and conversely, the upper bound ${\textit {Re}}_{cr-1}$ may be slightly higher.

We also note that the quantity $\varPsi$ is an indicative correlation between 2-D wake irregularity and the development of mode FW. This argument is because, firstly, the $\varPsi$ in (5.3) is proposed under the assumption that the onset of mode FW is triggered by flow instabilities originating from the 2-D base flow. As the lift coefficients on cylinders are global parameters, the $\varPsi$ may not be the only parameter that determines the 3-D instability. In addition, the critical value of $\varPsi \sim 0.0245$ is ad hoc selected through figure 12 and may not be valid in other 3-D transition problems.

Another phenomenon for the 3-D transition of FF flow in figure 5 is that, over the range of $g^* = 1.2\unicode{x2013} 1.8$, the flow transition sequence of ‘2-D FF $\to$ mode FW $\to$ 2-D IP $\to$ mode A-IP’ with increasing ${\textit {Re}}$. This phenomenon differs from the single-cylinder counterpart where the wake remains in 3-D flow and becomes increasingly irregular with increasing ${\textit {Re}}$. This transition sequence is further affirmed to be physical by conducting 3-D DNS for $g^* = 1.2$ and ${\textit {Re}} = 58\unicode{x2013} 190$, where 2-D and 3-D flows are plotted as discrete hollow and solid symbols in figure 5. The spanwise lengths $L_z = 24$ and $4{\rm \pi}$ are used for cases at ${\textit {Re}} < 100$ and ${\textit {Re}} > 100$, respectively. The distinct transition that exists here is because flow characteristics in FF and IP flows have different wake interference on the development of 3-D instability. Precisely, when the wake irregularity is relatively weak, mode FW transits back to a 2-D flow when ${\textit {Re}}$ approaches the upper bound of the FF flow region (dashed grey line in figure 5). Further increasing ${\textit {Re}}$ in the IP flow till ${\textit {Re}} \sim 180$, the elliptical instability of primary vortex cores in the IP wake develops, constituting the development of mode A-IP.

Detailed characteristics of mode FW are discussed using cases shown in figure 11. Since the transition of 2-D FF to mode FW cannot be predicted through Floquet analysis, we apply the Hilbert transform (Kress, Maz'ya & Kozlov Reference Kress, Maz'ya and Kozlov1989) to quantify the spatial distribution of spanwise length scales in mode FW and their evolutions with ${\textit {Re}}$. For each $(x, y)$ location, the spanwise wavelength is obtained by analysing the Hilbert transform of $\omega _y$ along $z$. The localised spanwise length scale $\lambda _z$ is estimated as

(5.4)\begin{equation} \frac{2{\rm \pi}}{\lambda_z} = \frac{\mathrm{d} \phi_{\omega_y}(z,t)}{\mathrm{d}z}, \end{equation}

where $\phi _{\omega _y}(z,t) \equiv \textrm {arg}[\omega _y(z,t)]$ is the spanwise phase of $\omega _y(z,t)$ calculated from the Hilbert transform as

(5.5)\begin{equation} \omega_y(z,t) \equiv A_{\omega_y}(z,t) \exp({{\rm i}\phi_{\omega_y}(z,t)}), \end{equation}

where $A_{\omega _y}(z, t) \equiv |\omega _y(z, t)|$ is the local amplitude. Following the work by Sarwar & Mellibovsky (Reference Sarwar and Mellibovsky2020), (5.4) and (5.5) are calculated by adopting the approximation of $\omega _y \approx \partial u/\partial z$. The probability density function (PDF) of $\lambda _z$ along the $z$ direction, namely $\lambda _z^P$, is then used to characterise the distributions of $\lambda _z$. This method is superior for characterising the length scale of the flow with localised flow features and has also been used in the studies of wake flows behind a bluff body (Sarwar & Mellibovsky Reference Sarwar and Mellibovsky2020; Ong & Yin Reference Ong and Yin2022).

As an example, figure 13(a) displays the $\lambda _z^P - t^*$ contour in a near wake region of $(x, y) = (8, 0)$ at $(g^*, {\textit {Re}}) = (1.0, 80)$, where the dark contours indicate the most representative length scales. The energetic spanwise structure in figure 13(a) has a length scale of $\lambda \sim 5D$. Figure 13(b) shows the averaged value of the PDF over 100 non-dimensional time units; and the maximum $\lambda _z^P$, $\lambda _{z,m}$, is used as the most representative spanwise length scale at ($x, y$). The $\lambda _{z,m}$ of mode FW along the wake centreline of $x/D = 2\unicode{x2013} 10$ for cases in figure 11 are shown in figure 13(c). At a fixed ${\textit {Re}}$, e.g. ${\textit {Re}} = 66$, the $\lambda _{z,m}$ decreases monotonically with $x/D$ for $x/D> 2$ due to the break up of large-scale vortices into small-scale vortices. With ascending ${\textit {Re}}$ from 66 to 100, $\lambda _{z,m}$ is reduced significantly and nearly unchanged over $x/D = 5\unicode{x2013} 10$. We acknowledge that the approximation of $\omega _y \approx \partial u/\partial z$ is used in estimating the characteristic length scale of mode FW. As demonstrated in Appendix B, this approximation may lead to a smaller $\lambda _{z,m}$. We recognise that the magnitude of $\lambda _{z,m}$ in figure 13(c) may vary slightly if different vorticity or velocity components are used. However, the trends of $\lambda _{z,m}$ with respect to $x/D$ and ${\textit {Re}}$ for mode FW would remain similar.

Figure 13. (a) Time evolution of PDF of the instantaneous spanwise wavelength $\lambda _z^P$ distribution at a point in the near wake region of $(x, y) = (8, 0)$ for mode FW at $(g^*, {\textit {Re}})= (1.0, 80)$. (b) Time-averaged PDF distribution over 100$t^*$. The maximum $\lambda _z^P$ ($\lambda _{z,m}$, the blue dot) is used as the representative spanwise length scale at this point. (c) The most probable wavelength, $\lambda _{z,m}$, along the wake centreline ($y/D = 0$) for mode FW at $(g^*, {\textit {Re}}) = (1.0, 66\unicode{x2013} 100)$.

The evolutions of instantaneous 3-D flow fields show that the 3-D structures of mode FW flap over intermittently between two wakes and persist different strengths over time. The evolution of 2-D flow fields indicates that the gap flow of FF undergoes oscillations over each vortex shedding cycle. At small $g^*$, the gap flow inclination leads to the weakening of vortex shedding around the cylinder in the direction of its inclination, while the vortex shedding on the other cylinder is strengthened. This strengthening/weakening process is also present in mode FW and is responsible for the flapping-over character of 3-D structures behind both wakes. This conjecture is validated through a case near the onset transition boundary, e.g. $(g^*, {\textit {Re}}) = (1.0, 66)$. Figure 14(a) displays the transverse velocity $v$ (grey line) for a point at the end of gap $(x/D, y/D, z/D) = (0.5, 0, 12)$. The averaged value of $v(t)$ between successive peaks and troughs (the dashed orange line), which represents the mean gap flow inclination direction in each half of the vortex shedding cycle, undergoes a slow and random flapping motion over several vortex shedding periods. Good correlations between the gap flow inclination (indicated by averaged $v(t)$) and strengthening/weakening of 3-D instabilities are observed. For instance, as the mean gap flow flaps to $y/D > 0$ from its equilibrium position (figure 14b–d), the vortex shedding from the top cylinder is weakened, characterised by an increase of wake formation length and the weakening of vortices (N1–N3) shed from the top shear layer. The 3-D structure in the wake of the top cylinder is weakened as well, as indicated by the weakening of $\omega _x$ fields in the shear layers of the top cylinder. During the same period of $t^* \approx 3010\unicode{x2013} 3030$, the strengthening of vortex shedding and 3-D instability from the bottom cylinder are observed. The observed behaviours over the next period of $t^* \approx 3030\unicode{x2013} 3050$, when the mean gap flow flaps to $y/D < 0$ (figure 14e–g) from its equilibrium position, are opposite to those in figure 14(c,d).

Figure 14. Evolutions of (a) the transverse velocity $v$ at the end of the gap of $(x/D, y/D, z/D) = (0.5, 0, 12)$ and (b–g) the flow fields $x$–$y$ plane at $z/D = 12$, sampled at the phases corresponding to the crests of $v(t)$ (blue circles) in (a) for the FF flow at $(g^*, {\textit {Re}})= (1.0, 66)$. The positive and negative $\omega _x$ fields are represented by the red and blue contours, respectively. The overlapped solid and dashed black lines are the corresponding positive and negative $\omega _z$ contours.

The above discussion suggests that, in two interacting wakes, the wake irregularity plays a pivotal role in triggering the 3-D instability. The onset transition from FF flow to mode FW develops at ${\textit {Re}}$ significantly lower than the single-cylinder counterpart $\{ Re_{cr-1} \}_S = 188.84$.

5.1.3. Coupled wakes

For the flows in IP, bistable IP/AP and AP states, the onset 3-D wake transitions are similar to the results of the flow around a single cylinder, where an unstable region with real and positive $\mu$ is identified at $\beta \approx 1.6\unicode{x2013} 1.8$. Mode A bifurcates slightly earlier in terms of ${\textit {Re}}$ than the single-cylinder case. With increasing ${\textit {Re}}$, the flow transition sequence is in ‘2-D $\to$ mode A $\to$ mode B’ for each type of base flow. The most unstable Floquet modes are similar to mode A of a single-cylinder case, except that the $\hat {\omega }_x$ fields behind both cylinders are modulated by the base flow symmetry conditions (4.1), (4.2). Similar 3-D transitions have been observed for steady flow past two side-by-side square cylinders (Choi & Yang Reference Choi and Yang2013). The effects of IP and AP wake couplings on the 3-D transitions of two circular cylinders in a side-by-side arrangement are discussed below.

The AP wake coupling destabilises the 3-D wake transition to a slightly lower ${\textit {Re}}$ level, as illustrated in figure 15. This phenomenon is attributed to the confinement effect of the AP flow, which strengthens the vortex shedding process of each wake. This conjecture is firstly supported by the vortex shedding frequency ($f_{St}$) shown in Ren et al. (Reference Ren, Cheng, Tong, Xiong and Tang2023), where $f_{St}$ is slightly larger than the single-cylinder value at the corresponding ${\textit {Re}}$. As the interaction of vortex cores in each wake of AP is stronger than the unconfined wake, it constitutes the earlier initiation of mode A. As $g^*$ ascends, the AP confinement effect reduces gradually and the ${\textit {Re}}_{cr-1}$ asymptotes to the values of a single cylinder. The differences of ${\textit {Re}}_{cr-1}$ in AP from $\{{\textit {Re}}_{cr-1} \}_S$ is around 2.3 % at $g^* = 3.0$, suggesting that there are negligible differences (within 3 % tolerance) in the 3-D transitions between an isolated cylinder and the AP wake for $g^* > \sim 3$. In addition, a similar variation of ${\textit {Re}}_{cr-1}$ over $g^*$ has been found in other 3-D transition studies on a confined cylinder wake, i.e. the study of steady flow past an isolated cylinder near a moving wall (Jiang et al. Reference Jiang, Cheng, Draper and An2017a,Reference Jiang, Cheng, Draper and Anb), in which the variation was attributed to the reduced confinement of the moving wall. Given the cylinder wakes in both scenarios are confined by the presence of a moving wall (Rao et al. Reference Rao, Thompson, Leweke and Hourigan2015; Jiang et al. Reference Jiang, Cheng, Draper and An2017a) or the other wake (the present scenario), we believe the increase trend of ${\textit {Re}}_{cr-1}$ over $g^*$ in AP flow is indeed affected by the weakening of wake confinement. To further compare the wall confinement effect on two scenarios, the onset ${\textit {Re}}_{cr-1} - g^*$ curves for the AP base flow and moving wall case (Jiang et al. Reference Jiang, Cheng, Draper and An2017a) are shown in figure 15. We observe that, at $g^* > 1.5$, both the magnitude of ${\textit {Re}}_{cr-1}$ and its variation trend with respect to $g^*$ for AP base flow and the moving wall case at an equivalent gap distance are close, with the largest difference of 1.27 % occurring at $g^*=2.0$. This good agreement also suggests the confinement in the AP wake is comparable to the cylinder next to a moving wall.

Figure 15. The influence of wake confinement on the critical ${\textit {Re}}$ of 3-D wake transitions of a steady flow past two side-by-side cylinders in the AP flow (present) and an isolated cylinder near a moving wall (Rao et al. Reference Rao, Thompson, Leweke and Hourigan2015; Jiang et al. Reference Jiang, Cheng, Draper and An2017a,Reference Jiang, Cheng, Draper and Anb). The gap-to-diameter ratio $G^*$ in Jiang et al. (Reference Jiang, Cheng, Draper and An2017a,Reference Jiang, Cheng, Draper and Anb) and Rao et al. (Reference Rao, Thompson, Leweke and Hourigan2015) is translated to $g^* = 2G^*$.

For the first 3-D instability of the IP base flow, i.e. mode A, the ${\textit {Re}}_{cr-1}-g^*$ curve is very close to that of the single cylinder. This result is expected because the two IP wakes oscillate together with almost no confinement effect between wakes, compared with the strong confinement effect in the AP state. Even though the IP state has stronger wake coupling than the AP state, judging based on the sequence of their emergence in terms of enlarging $g^*$, the coupling is mainly related to the vortex merging in the downstream, which has a minor effect on the vortex shedding processes in the near wake. Consequently, the development of the first 3-D instability (mode A-IP) is less affected.

The above observations also suggest the wake confinement in AP is more favourable of interfering the 3-D instability than the vortex merging in the IP state.

5.2.1. Single bluff-body wake

The second 3-D instability for the ASS flow develops at large $\beta$ number, as shown in figure 8(a), whose Floquet multiplier $\mu$ in the peak region is real and negative. Spatial and temporal distributions of a typical secondary unstable mode in ASS at $(g^*, Re, \beta ) = (0.3, 100, 2.4)$ are shown in figure 16. A period-doubling feature of the perturbed Floquet mode is identified. The symmetry condition of perturbed $\hat {\omega }_x$ field can be expressed as

(5.6)\begin{equation} \mbox{Mode C}': \begin{cases} \hat{\omega}_x(x,y,t)={-}\hat{\omega}_x(x,y,t+T_v), \\ \hat{\omega}_x(x,y,t)=\hat{\omega}_x(x,y,t+2T_v). \end{cases} \end{equation}

Equation (5.6) suggests that this mode takes two periods of vortex shedding in the $x$–$y$ plane to complete one cycle of spanwise structure evolution in the wake. The 3-D instability that holds the same RT symmetry condition as (5.6) is commonly denoted as the mode C-type instability in literature (Zhang et al. Reference Zhang, Fey, Noack, König and Eckelmann1995; Carmo et al. Reference Carmo, Sherwin, Bearman and Willden2008; Sheard et al. Reference Sheard, Fitzgerald and Ryan2009). For this reason, we refer to this mode to as ‘mode C$'$’. Similar to mode B of an isolated cylinder, the $\hat {\omega }_x$ field of mode C$'$ is concentrated in the braids.

Figure 16. (a) Spatial and (b) temporal distributions for mode C$'$-ASS at $(g^*, Re, \beta ) = (0.3, 100, 2.4)$. (c) The $|\mu | - \beta$ variations for two side-by-side cylinders with a diameter ratio of 1/0.8 at $g^* = 0.3$. (d) Spatial distribution for the unstable mode (mode C$'$ type) at $(Re, \beta ) = (100, 2.4)$.

Bearing in mind that mode C is absent in flow around an isolated cylinder, we anticipate that the formation of an asymmetric wake induced by the biased gap flow is responsible for the occurrence of mode C$'$ in the ASS flow. One evidence for this conjecture is from figure 9(b), where large magnitudes of $\varGamma$ and $\overline {U_{gap}}$ are observed when the ASS base flow forms mode C$'$ (half-filled symbols). To further confirm that the present mode C$'$ has similar features to the typical mode C that is triggered by asymmetric geometry configurations, we artificially reduce the diameter of the bottom cylinder to $0.8D$, while keeping the gap distance as $0.3D$. Floquet analyses are conducted over the range of ${\textit {Re}}$ ($= 90\unicode{x2013} 105$) to assess if mode C$'$ would develop. The $|\mu | - \beta$ variations in figure 16(c) show that an unstable region is identified in the vicinity of $\beta \sim 2.4$. Corresponding spatial distribution of the Floquet mode at $({\textit {Re}}, \beta ) = (100, 2.4)$ is shown in figure 16(d). The similarities in the $\hat {\omega }_x$ field and the $\beta$ ranges between the unstable modes in figure 16(c,d) and mode C$'$-ASS in figures 8(a) and 16(a) are clearly observed. In light of the above, it is therefore reasonable to conclude that mode C$'$ would develop when the wake asymmetry, induced by either the flow or geometry configuration, exceeds a critical value.

5.2.2. Coupled wakes

For the IP and AP coupled wakes, the second 3-D instability with real and positive $\mu$ is of mode B type behind each cylinder, where $\hat {\omega }_x$ field intensifies in the near wake region and rapidly diminishes downstream of the cylinders. Due to the IP and AP coupling between two wakes, the mode B structure behind each cylinder is modified by the base flow symmetry conditions so that the mode B of AP and IP flows show a reflection symmetry and a spatio-temporal symmetry along $y_c/D = 0$, respectively.

The onset of mode B for the AP base flow is delayed to much higher ${\textit {Re}}$, e.g. ${\textit {Re}}_{cr-2} \approx 359$ for $g^* = 1.4$ in figure 15. With an increase of $g^*$, the ${\textit {Re}}_{cr-2}$ of AP decreases monotonically to the value that is comparable to ${\textit {Re}}_{cr-2}$ of a single cylinder. The mechanism responsible for the large ${\textit {Re}}_{cr-2}$ is investigated by comparing the base flow characteristics of a single wake and the AP coupled wakes at ${\textit {Re}} = 300$ (see figures 3 and 4h). It is evident that the AP flow is featured by a stronger vortex shedding process than the single cylinder judging based on the strong vortices shed from the shear layers in figure 4(h). The braid shear layer of the AP flow, which connects adjacent vortex cores, is somewhat weakened. Since mode B is originated from the hyperbolic instability of the braid shear layer region (Williamson Reference Williamson1996a), the initiation of mode B is therefore delayed to much higher ${\textit {Re}}$ in AP. We attribute the weakening of the braid shear layer in AP to the wake confinement effect. This conjecture can be supported by the previous study of steady flow past an isolated cylinder near a moving wall (Rao et al. Reference Rao, Thompson, Leweke and Hourigan2015), where ${\textit {Re}}_{cr-2}$ gradually approaches $\{{\textit {Re}}_{cr-2}\}_S$ when the wall confinement is reducing (solid purple line in figure 15). The variation trend and magnitudes of ${\textit {Re}}_{cr-2} - g^*$ for the present AP flow are close to those given in Rao et al. (Reference Rao, Thompson, Leweke and Hourigan2015), suggesting the significantly large ${\textit {Re}}_{cr-2}$ value is due to the weakening of wake confinement.

The ${\textit {Re}}_{cr-2}$ for IP flow is close to the $\{ {\textit {Re}}_{cr-2} \}_S$ in figure 5. This is possibly because the IP flow is less affected by the confinement mechanism than that of the AP. With increasing $g^*$, the ${\textit {Re}}_{cr-2}$ of IP flow decreases monotonically due to the weakening of coupling strength in IP flow.

5.3.1. Nonlinear transitions

Nonlinear transitions are examined first for a case at $(g^*, {\textit {Re}})= (0.3, 95)$ where only mode A$'$ is unstable based on Floquet analysis. The spanwise length is selected as $L_z = 28D$, which is around four times of the dominant spanwise wavelength quantified in Floquet analysis.

The transition feature is demonstrated through the evolution of spanwise velocity $w$ at a point $(x/D, y/D, z/D) = (3, 0, 14.875)$ (figure 17a) and the instantaneous flow fields in figure 17(b) and figure 18. The nonlinear evolution for the mode A$'$ flow evolves into the following three distinct stages.

1. (i) Stage I at $t^* \lessapprox 2400$: the mode A$'$ structure gradually develops from a 2-D flow field. An example of the flow at $t^* = 2200$ is shown in figure 17(a), where regular 3-D mode A$'$ structures are observed and agrees well with the Floquet analysis in many aspects. Firstly, the streamwise vortices in the braid shear layers remain the same sign at a constant $z/D$ (see the $z/D = 12$ slice at $t^* = 2200$ in figure 17b-iv). Secondly, a total of four pairs of the mode A$'$ structure develop along the cylinder span of $28D$. The resulting spanwise length of ${\sim }7D$ is close to the $\lambda _{cr-1} = 8.64$ predicted from the Floquet analysis (see figure 8b).

2. (ii) Stage II at $t^* \approx 2400\unicode{x2013} 3100$: transient stage where spanwise disorder mode A$'$ structures are developed (see figure 17b-ii, b-iii).

3. (iii) Stage III at $t^* \gtrapprox 3100$: the 3-D structure saturates into a new stable state, featured by the periodic fluctuations of velocity and force evolutions.

Figure 17. (a) The evolution of spanwise velocity ($w$) in a point $(x/D, y/D, z/D) = (3, 0, 14.875)$ of a flow at $(g^*, {\textit {Re}})= (0.3, 95)$. (b) Flow evolution in the $x$–$z$ plane (top) and a slice at $z/D = 12$ (bottom). Isosurfaces of $\omega _x = \pm 0.1$ (blue and red) and $\omega _z = \pm 0.5$ (grey).

Figure 18. Instantaneous mode ASS structure for the case $(g^*, {\textit {Re}})= (0.3, 95)$ at $t^* = 3500$. Flow structures in the (a) 3-D view and (b) $x$–$z$ plane. Isosurfaces of $\omega _x = \pm 0.1$ (blue and red) and $\omega _z = \pm 0.5$ (grey). (c) A global view of the gap shear layers, represented by $\omega _z = \pm 0.5$ isosurfaces and coloured by $v$ velocity contours. (d) Contours of $\omega _x$ fields at $z/D = 7, 14, 21$, where yellow arrows represent the gap flow inclination directions.

The characteristics of the 3-D structure in stage III is investigated through the instantaneous flow field at $t^* = 3500$ (figure 18). The 3-D field in the downstream ($x/D > \sim 15$ in figure 18a,b) holds a similar spatio-temporal symmetry to the conventional mode A-type flow, with four pairs of 3-D structures identified along the cylinder span ($\lambda _z \sim 7D$). A distinct flow feature in the saturated stage III is discovered in the gap flow. As illustrated by the $|\omega _z| = 0.5$ isosurfaces coloured by $v$ velocity shown in figure 18(c) and $\omega _z$ and $\omega _x$ fields at different $z$ locations (figure 18d), the gap flow along the span actually flaps alternatively between $+y$ and $-y$ directions, forming an ordered wave-like pattern with a wavelength of ${\sim }14D$. Since the waviness of the gap shear layers in the spanwise direction has not been reported in the side-by-side cylinders, to the best of the authors knowledge, we refer to this mode as ‘mode ASS’.

The nonlinear evolution of the mode C$'$ in ASS flow is then assessed through a case at $(g^*, {\textit {Re}})= (0.3, 100)$ (slightly above the ${\textit {Re}}_{cr-2}$ of 97.8 at $g^* = 0.3$). A large $L_z = 28D$ is used to make sure the nonlinear interaction between the first and second 3-D modes can be well captured by the simulation. The velocity evolutions at $(x/D, y/D, z/D) = (3, 0, 11.81)$ (figure 19a) suggest mode C$'$ undergoes a similar transition to those discussed in mode A$'$, except that the mode C$'$ with period-doubling characteristics is predominant at the onset of 3-D flow (stage I). The sign of $\omega _x$ in the braid shear layers changes between two successive vortex shedding cycles in mode C$'$, e.g. flow fields at $t^* = 1500$ and $1500+T_v$ in figure 19(b-i, b-ii). The period-doubling characteristics can also be revealed by the evolution of $u(t)$ in figure 19(a). Ten pairs of mode C$'$ structure are identified along the cylinder span of $28D$, resulting in a spanwise length scale of ${\sim }2.8D$. The above features of mode C$'$ agree well with those predicted from Floquet analysis. However, the 3-D DNS indicates that mode C$'$ cannot persist and is gradually replaced by mode A$'$ at $t^* \gtrapprox 1700$ (see stage II in figure 19c) followed by spanwise dislocations due to the nonlinear interaction (see stage III). The flow finally saturates to mode ASS in figure 19(d).

Figure 19. (a) The evolution of streamwise ($u$) and spanwise velocities ($w$) at $(x/D, y/D, z/D) = (3, 0, 11.81)$ of the ASS flow at $(g^*, {\textit {Re}})= (0.3, 100)$. (b–d) Instantaneous flow fields in the $x$–$z$ plane (top) and a slice at $z/D = 12.5$ (bottom). Isosurfaces of $\omega _x = \pm 0.1$ (blue and red) and $\omega _z = \pm 0.5$ (grey).

5.3.2. Spatio-temporal symmetry of gap shear layers

The spatio-temporal symmetry for mode A$'$, mode C$'$ and mode ASS are explored through examining the evolutions of $|\omega _z| = 0.5$ isosurfaces coloured by $v$ velocity in the gap region ($x/D <3$) in figure 20. The striking features are summarised below.

Figure 20. Evolution of gap shear layers over one vortex shedding cycle. (a) Mode A$'$ $(g^*, {\textit {Re}})= (0.3, 95)$, where $T_v = 11.2$; (b) mode C$'$ $(g^*, {\textit {Re}})= (0.3, 100)$, where $T_v = 11$; and (c) mode ASS $(g^*, {\textit {Re}})= (0.3, 95)$, where $T_v = 12$. The evolutions of gap shear layers are represented by $|\omega _z| = 0.5$ isosurfaces coloured by $v$ velocity contours at $x/D <3$.

1. (i) The flow fields of mode A$'$ sampled at an interval of $1/2T_v$ ($T_v = 11.2$) (figure 20a) indicate that its gap shear layers are permanently biased towards one of the cylinders along the cylinder span. The periodicity of mode A$'$ is identical to the vortex shedding period in 2-D base flow.

2. (ii) Mode C$'$ is distinct in its $2T_v$ periodicity of the spanwise structure. As shown by instantaneous fields over $2T_v$ (figure 20b), the $\omega _z$ isosurfaces at two instants that are $1T_v$ apart from each other form opposite waviness with respect to the $y$ axis (see the arrows with the same colours in figure 20b). A close visualisation of flow fields over $2T_v$ confirms that the three vorticity components share the same spatio-temporal symmetry conditions,

   (5.7)\begin{equation} \mbox{Mode C}': \quad \boldsymbol{\omega}(x,y,z,t)=\boldsymbol{\omega}(x,y,z\pm\lambda_z/2,t+T_v). \end{equation}
   Note that the $2T_v$ periodicity of mode C$'$ is revealed on time histories of velocities at a sampling point (the purple arrow in figure 19a) but not on the spanwise-averaged forces. This is because the mode C$'$ pair cancels with each other and the $2T_v$ periodicity can not be identified through spanwise-averaged quantities. Another striking feature is that, even though the spanwise waviness of the gap shear layers fluctuates with the vortex shedding process, the gap flow at different $z$ is permanently biased towards one of the cylinders.

3. (iii) As for mode ASS (figure 20c), the flow through the gap forms significant waviness along the cylinder span and towards ${\pm }y$ directions. The $v$-velocity distributions that differ by 1$T_v$ ($= 12$) are identical, suggesting that mode ASS is periodic in each vortex shedding period. The flow fields sampled at an interval of $1/8T_v$ over 1$T_v$ indicate that this wave-like pattern is time independent. The spatio-temporal symmetry along the spanwise direction can be described by

   (5.8)\begin{equation} \mbox{Mode ASS:}\quad \boldsymbol{\omega}(x,y,z,t)= \boldsymbol{\omega}(x,-y,z\pm\lambda_z/2,t), \end{equation}
   where $\lambda _z$ has a wavelength around twice of that identified in mode A$'$.

We would like to point out that the nonlinear transition to mode ASS is dependent on the $L_z$ value adopted in the 3-D DNS. A series of simulations at $(g^*, {\textit {Re}})= (0.3, 95)$ and (0.3, 100) with $L_z$ varying from 7 to 42 were conducted to investigate this dependence. We found that mode ASS is reproduced only when $L_z$ is larger than twice of the $\lambda _z$ for mode A$'$ ($L_z \gtrapprox 14D$). We suspect that the flow transition to mode ASS is due to the nonlinear interaction of two mode A$'$ pairs.

5.3.3. Hysteresis effect

The existence of hysteresis for the fully developed flow, mode ASS, is explored. The hysteresis effect is first assessed through the Landau model (Landau & Lifshitz Reference Landau and Lifshitz1976) using the case at $(g^*, {\textit {Re}})= (0.3, 95)$ where the flow is unable to achieve mode ASS. The Landau equation, which considers the evolution of perturbation amplitude $A(t)$, is written as

(5.9)\begin{equation} \mathrm{d}A/\mathrm{d}t = (\sigma+{\rm i}\omega)A-l(1+{\rm i}c)|A|^2 A+\cdots, \end{equation}

where $\sigma$ and $\omega$ represent the linear growth rate and angular oscillation frequency of the perturbation; $c$ is the Landau constant. The sign of $l$, which defines the bifurcation nature of a transition, can be determined by the plot of $\mathrm {d}\log |A|/\mathrm {d}t$ over $|A|^2$, where the slope of the curve at $|A|^2 \sim 0$ gives the $-l$. The transition is supercritical (non-hysteresis) and subcritical (hysteresis) if $l>0$ and $l<0$, respectively. More details of this method can be found in Carmo et al. (Reference Carmo, Sherwin, Bearman and Willden2008, Reference Carmo, Meneghini and Sherwin2010) and Jiang & Cheng (Reference Jiang and Cheng2020). In the present case the evolution of $|A|$ is measured using the envelope of spanwise velocity $w$ at a point $(x/D, y/D, z/D) = (3, 0, 14.875)$ in a stage where mode A$'$ is developing, stage I shown in figure 17(a). The curve of $\mathrm {d}\log |A|/\mathrm {d}t$ over $|A|^2$ in figure 21(b) shows a positive slope close to the $y$ axis ($l>0$), suggesting that mode ASS is subcritical in nature.

Figure 21. The examination of the hysteresis effect for the flow at $(g^*, {\textit {Re}})= (0.3, 95)$. (a) A zoomed-in view of $w(t)$ in figure 17(a) at $t^* = 1500\unicode{x2013} 2000$, where the envelope of $w(t)$, $A$, are shown as the solid red line. (b) The derivative of the amplitude logarithm over the square of the amplitude, suggesting a subcritical transition ($l<0$).

This subcritical transition nature is further confirmed through 3-D DNS at $g^* = 0.3$ by initialising the flow from the following three ICs.

1. (1) Initial condition $Re_H$ (IC-$Re_H$): the flow is initialised from a fully developed 3-D flow at $Re = 95$ (in mode ASS at the time instant shown in figure 17) and gradually increases $Re$ by an interval of $\Delta Re = 5$ until a 2-D wake is developed.

2. (2) Initial condition $Re_L$ (IC-$Re_L$): the flow is initialised from a 2-D flow at $Re = 85$ and gradually increases $Re$ by an interval of $\Delta Re = 5$ until $Re = 95$.

3. (3) Initial condition 2-D (IC-2-D): the flow is initialised from a 2-D state at the corresponding $Re$.

The distributions of 2-D and 3-D wakes for the above three ICs are marked as discrete symbols in figure 22. The arrows on the IC-$Re_L$ and IC-$Re_H$ branches indicate the directions of successive initialisation. The wake is deemed a 3-D wake first by assessing the $\varepsilon _{xy}/\varepsilon _t > 0.001$. Figure 22 shows that three IC branches have different lower-bound ${\textit {Re}}$ values for the onset of wake three-dimensionality, forming a hysteresis loop. Precisely, the IC-2-D and IC-$Re_L$ both transit at ${\textit {Re}} \approx 87.5$ that agree well with the marginal curve estimated by Floquet analysis (${\textit {Re}}_{cr-1} = 87.12$), while the flow remains a 3-D flow till ${\textit {Re}} = 80$ for cases with IC-$Re_H$ due to the large perturbation introduced by the mode ASS flow. The hysteresis loop observed here confirms the subcritical transition identified from the Landau equation.

Figure 22. Comparisons of the wake state distributions at $g^* = 0.3$ by initialising the flow from different ICs. The 2-D and 3-D wakes are represented by horizontal lines. The thick grey line indicates the presence of a wake state based on 3-D DNS using a 2-D state at the corresponding $Re$ (IC-2-D). The orange and blue solid lines denote the dominance of wake states based on 3-D DNS initialised by a fully developed flow field from adjacent higher ${\textit {Re}}$ (IC-${\textit {Re}}_{H}$) and lower ${\textit {Re}}$ (IC-${\textit {Re}}_{L}$) values, respectively. The stating points and directions for successive initialisation are marked by the solid star symbols and arrows, respectively.