3.1. Nonlinear base flow

The evolution of infinitesimal perturbations to a base state constitutes the scope of linear stability analysis. Thus, the base flow about which the NS equations are linearised needs to be extracted. As no analytical solution exists, the steady base flow ($\boldsymbol {U}, P$) at $Re_D = 50$ is extracted numerically from nonlinear NS equations. The selective frequency damping (SFD) by Åkervik et al. (Reference Åkervik, Brandt, Henningson, Hæpffner, Marxen and Schlatter2006) is adopted. The SFD technique consists of an additional forcing term $\boldsymbol {f}$, which damps the oscillations of the solution using a temporal low-pass filter. The forcing term is defined as $\boldsymbol {f} = -\chi (\boldsymbol {u}-\boldsymbol {\omega })$, where $\boldsymbol {u}$ is the flow solution and $\boldsymbol {\omega }$ is the temporally low-pass-filtered velocity obtained by a differential exponential filter $\boldsymbol {\omega }_t=(\boldsymbol {u}-\boldsymbol {\omega })/\varDelta$, with $\varDelta$ determining the filter width. The base flow is extracted when the tolerance $\varepsilon = \| \boldsymbol {u}-\boldsymbol {\omega } \|_{L^2(\varOmega )}$ falls below a value of $\varepsilon < 10^{-7}$. The robustness of the results is assessed by comparing the results with additional base flows obtained with a lower tolerance $\varepsilon < 10^{-9}$. From the nonlinear NS equations

(3.1) \begin{equation} \left. \begin{array}{c@{}} \displaystyle \dfrac{\partial \boldsymbol{u} }{\partial t} + ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{u} ={-} \boldsymbol{\nabla} p + \dfrac{1}{Re_D} \nabla^2 \boldsymbol{u} + \boldsymbol{f}, \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{array} \right\}\end{equation}

the base flow $\boldsymbol {u}$ is calculated and, hereinafter, labelled $\boldsymbol {U}$. For each case, the base flow is symmetric about the $y=0$ plane, but the streamwise velocity $u/U$ shows significant differences, see figure 1. In particular, the recirculation region ($u/U<0$) becomes progressively thinner behind the small cylinder as $r$ increases. As expected, the extension behind the large cylinder is almost unchanged. It is also worth observing the vertical velocity $W/U$ variation in the front part and in the near wake of the cylinder. The spanwise variation at $(x/D=-0.6,y/D=0)$ indicates an upward flow in front of the junction, where the impingement of the large cylinder gets the flow wrapped around the junction. This becomes gradually stronger as $r$ increases, i.e. as the junction surface gets larger. Similarly, in the near wake $(x/D=3,y/D=0)$, the downwash phenomenon is visible, as highlighted by the negative velocity component $W/U$ in figure 2.

Figure 1. Streamwise velocity component of the base flow ($y=0$ plane) at $Re_D=50$ for four different ratios $r=1.1$, $1.2$, $2$ and $4$ in (a), (b), (c) and (d), respectively. The colour maps are skew-symmetric to have $u/U=0$ in white.

Figure 2. Spanwise variation of the vertical velocity component of the base flow at $Re_D=50$ for the ratios $r=1.1$, $1.2$, $2$ and $4$. Two different locations are considered: $(x/D=-0.6,y=0)$ and $(x/D=3,y=0)$ with solid and dashed lines, respectively.

Interestingly, the spanwise extension of the downward flow is slightly affected by $r$ and corresponds to the N cell observed at higher Reynolds numbers (Massaro et al. Reference Massaro, Peplinski and Schlatter2023c). Considering the finding by Massaro et al. (Reference Massaro, Peplinski and Schlatter2023d), who documented the connection between the downwash phenomenon and the N cell, the downward motion could entail the formation of an N cell for any $r$. This would be in contrast with the classification by Lewis & Gharib (Reference Lewis and Gharib1992) in the laminar vortex shedding regime, who described a direct modes interaction (no N cell) for $r<1.25$. Further evidence is reported in the following sections. Note that the intensity of the downwash is larger when $r$ increases, with the negative peak of $W$ being around 1.8 %, 3 %, 11 % and 17 % of the inflow velocity for $r=1.1$, $1.2$, $2$ and $4$, respectively. The effect of the spanwise extensions of the small and large cylinders has been carefully assessed in the Appendix. It is also worth noting how the shape of the vertical velocity profile changes (see figure 2). For $r=1.1$ and $r=1.2$, the velocity profile is symmetric about the $z=0$ plane, in contrast to $r=2$ and $r=4$, where it is deflected behind the large cylinder ($z<0$). Given the connection between the downwash and the modulation cell, the symmetric downwash for $r=1.1,1.2$ suggests the appearance of a second modulation cell, as discussed later.

3.2. Global stability analysis

The governing equations for the perturbations ($\boldsymbol {u}',p'$) are obtained by linearising the NS equations about the extracted base flow $\boldsymbol {U}$

(3.2) \begin{equation} \left. \begin{array}{c@{}} \dfrac{\partial \boldsymbol{u}'}{\partial t} + (\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{U} + (\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}'={-}\boldsymbol{\nabla} p' + \dfrac{1}{Re_D} \nabla^2 \boldsymbol{u}', \\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}'=0, \end{array} \right\}\end{equation}

completed by the proper initial and boundary conditions. The dual solution ($\boldsymbol {u}^{{\dagger} },p^{{\dagger} }$) for the linearised adjoint set of equations is also computed

(3.3) \begin{equation} \left. \begin{array}{c@{}} -\dfrac{\partial \boldsymbol{u}^{{{\dagger}}}}{\partial t} - (\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla})\ \boldsymbol{u}^{{{\dagger}}} + (\boldsymbol{\nabla} \boldsymbol{U})^T \boldsymbol{u}^{{{\dagger}}} = \boldsymbol{\nabla} p^{{{\dagger}}} + \dfrac{1}{Re_D} \nabla^2 \boldsymbol{u}^{{{\dagger}}}, \\ -\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^{{{\dagger}}}=0, \end{array} \right\}\end{equation}

and their global stability is studied by evaluating the eigenvalues, and corresponding eigenvectors, for the linearised direct and adjoint NS operators. Under the normal-mode hypothesis, we express the velocity and pressure disturbance of the linear problem (3.2) as

(3.4) \begin{equation} \left.\begin{array}{c@{}} \boldsymbol{u}'(\boldsymbol{x},t) = \hat{\boldsymbol{u}}(\boldsymbol{x}) \,{\rm e}^{\lambda t}\quad \lambda \in \mathbb{C} \\ p'(\boldsymbol{x},t) = \hat{p}(\boldsymbol{x}) \,{\rm e}^{\lambda t } \quad \lambda \in \mathbb{C} \end{array}\right\},\end{equation}

with $\lambda =\sigma +\textrm {i} \omega$. The same approach is applied to the problem (3.3). Afterwards, by substituting (3.4) in (3.2), a generalised eigenvalue problem is formulated as

(3.5)\begin{equation} \lambda \mathcal{R} \hat{\boldsymbol{r}} = \mathcal{J} \hat{\boldsymbol{r}},\end{equation}

where

(3.6a–c)\begin{equation} \hat{\boldsymbol{r}} = \begin{pmatrix} \hat{\boldsymbol{u}} \\ \hat{p} \end{pmatrix}, \quad \mathcal{R} = \begin{pmatrix} \mathcal{I} & 0 \\ 0 & 0 \end{pmatrix}, \quad \mathcal{J} = \begin{pmatrix} -\boldsymbol{\nabla} \boldsymbol{U} -\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} + \dfrac{1}{Re_D} \nabla^2 & -\boldsymbol{\nabla} \\ \boldsymbol{\nabla} \boldsymbol{\cdot} & 0 \end{pmatrix}. \end{equation}

However, assembling these matrices is computationally prohibitive. Thus, the problem (3.5) is recast as an initial value problem for the velocity only by exploiting the incompressibility constraint. The solution reads

(3.7)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = {\rm e}^{\mathcal{L} t} \boldsymbol{u}_0(\boldsymbol{x}),\end{equation}

where $\mathcal {L}$ is the projection of $\mathcal {J}$ on a divergence-free space and $k$ are the eigenvalues of the matrix exponential $\textrm {e}^{\mathcal {L} t}$, related to those of $\mathcal {J}$ by the expressions

(3.8a,b)\begin{equation} \sigma = \dfrac{\ln (|k|)}{\Delta t}, \quad \omega=\dfrac{\text{arg}(k)}{\Delta t},\end{equation}

where $\Delta t$ is the time interval between the Krylov vectors generated in the time-stepper approach (Tuckerman & Barkley Reference Tuckerman and Barkley2000). To compute the eigenpairs, the implicitly restarted Arnoldi method (IRAM), proposed by Sorensen (Reference Sorensen1992), is used. It is implemented in the software package ARPACK (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998), and integrated in the KTH framework for Nek5000 (Peplinski, Schlatter & Henningson Reference Peplinski, Schlatter and Henningson2015; Massaro et al. Reference Massaro, Peplinski, Stanly, Mirzareza, Lupi, Mukha and Schlatter2024). The Arnoldi method approximates eigenpairs by searching for solutions within the Krylov subspace. In this work, the Krylov subspace has a dimension of $m=100$. We compute the initial $20$ eigenpairs for both the direct and adjoint problems, ensuring a residual tolerance of $10^{-6}$ for the eigenvalue calculation. Note that larger Krylov subspaces were also considered to assess the convergence of results.

For all the $r$ ratios, the spectra are observed to have only one unstable pair of complex conjugate eigenvalues ($\sigma > 0$) at $Re_D=50$, see figure 3. Thus, the flow undergoes a supercritical Hopf bifurcation between $Re_D= 40$ and $Re_D= 50$, always with a growth rate $\sigma \approx 0.01$ for any $r$. The growth rate slightly decreases when $r$ increases, likely due to the more prominent effect of the junction surface, as discussed later. The angular frequencies are $\omega = 0.746, 0.742, 0.728$ and $0.721$, corresponding to a Strouhal $St \approx 0.11$, i.e. $T \approx 8$ convective time units, based on $D$ and $U$. The results are in agreement with the circular cylinder by Giannetti & Luchini (Reference Giannetti and Luchini2007) and local velocity probes in the nonlinear simulation confirm the Strouhal obtained by the linear analysis as the flow undergoes from a steady state to a periodic regime. The critical Reynolds number is linearly interpolated and the neutral curve (figure 5) shows a critical $Re_D \approx 47$ for any geometries. The adjoint spectrum nearly overlaps the direct one, supporting the convergence of our numerical results (Robinson Reference Robinson2020). All the results are fairly close to the circular cylinder, indicating that the global stability mechanism is still two-dimensional and independent of the diameter ratios.

Figure 3. Portion of the spectrum for the linearised direct and dual NS operators for $r=1.1$, $1.2$, $2$ and $4$, in (a), (b), (c) and (d), respectively. Note that for $r=2$ and $r=4$, the mode corresponding to the small cylinder instability is not present at $Re_D=50$ since the local Reynolds number is significantly lower than the critical value.

3.3. Structural sensitivity

The character of the NS linearised operator is highly non-normal, leading to the substantial spatial separation between the direct and adjoint eigenmodes. Given that, to identify the origin of the instability, Giannetti & Luchini (Reference Giannetti and Luchini2007) introduced the idea of a wavemaker for global modes for the flow around a circular cylinder. The wavemaker pinpoints regions in the flow where the instability mechanism acts to give rise to self-sustained oscillations, i.e. where generic structural modifications of the stability problem result in the most significant shift in the leading eigenvalue. This shift implies that the structural perturbation directly affects the fundamental aspects of the instability mechanism. By computing the structural sensitivity function $\eta$, we can specify the locations where the feedback is most pronounced, i.e. where the instability mechanism is active. Giannetti & Luchini (Reference Giannetti and Luchini2007) define $\eta$ as

(3.9)\begin{equation} \eta(\boldsymbol{x}) = \dfrac{\| \hat{\boldsymbol{u}}^{{{\dagger}}}(\boldsymbol{x}) \| \| \hat{\boldsymbol{u}}(\boldsymbol{x}) \|}{\displaystyle \int_{\varOmega} \hat{\boldsymbol{u}}^{{{\dagger}}}(\boldsymbol{x}) \hat{\boldsymbol{u}}'(\boldsymbol{x}) \,\text{d}\varOmega},\end{equation}

where $\varOmega$ is the computational domain and $\| {\cdot } \|$ is the magnitude.

Figure 4 shows the wavemaker, together with the unstable direct and dual eigenmode for $r=2$. Similar results are obtained for $r=1.1$, $2$ and $4$. The wavemaker region is completely enclosed in the L cell and characterised by two lobes symmetrically placed across the separation bubble (see the enlarged view in figure 4), matching the area enclosed by the two streamlines that separate from the surface of the cylinder. These findings are in excellent agreement with circular cylinder (Giannetti & Luchini Reference Giannetti and Luchini2007), confirming the two-dimensional nature of the global instability mechanism. Note that the wavemaker identifies the S and N cells when computed for the second and third least unstable eigenmode for various $r$.

Figure 4. Isosurfaces of the negative and positive streamwise velocity (50 % of the maximum) of the unstable eigenmode of the direct (purple/orange) and dual (blue/red) operators. In green, the isosurface of the wavemaker $\eta =0.8$ ($0<\eta <1$). In the enlarged view, the adjoint solution at $z/D=-8$ with a black isoline for $\eta =0.8$ can be seen. The black arrow indicates the direction of the homogeneous inflow $U$.

3.4. Neutral curve

The neutral curve in the ($Re_D$, $r$) space is presented in figure 5. The first least-stable eigenmode has a critical Reynolds number of approximately $47$ for any geometries, with negligible deviation from the two-dimensional prediction (dashed line in figure 5). For the second least-stable eigenmode, the critical Reynolds number ($Re_{D,cr2}$) significantly varies with $r$. An increase of approximately 30 % is observed from $r=1.1$ to $r=4$. Specifically, $Re_{D,cr2}=50.97$, $57.13$, $60.97$ and $66.88$ for $r=1.1$, $1.2$, $2$ and $4$, respectively. This trend differs from the expected values ($\widehat {Re}_{d,cr2} = Re_{D,cr} \cdot r$) as $r$ increases. For $r=1.1$ and $r=1.2$, the critical Reynolds number is similar to the estimation based on the local diameter ($\widehat {Re}_{d,cr2} \approx 51.7$ and $56.4$). The first two red squares align with the dot-dashed line in figure 5. In contrast, for $r=2$ and $r=4$, the estimated critical Reynolds number is substantially lower than the expected critical Reynolds number based on the local diameter $d$ ($\widehat {Re}_{d,cr2} \approx 94$ and $188$, for $r=2$ and $r=4$, respectively). Indeed, the dot-dashed line diverges from the results of global stability analysis in these cases.

Figure 5. The critical Reynolds number as a function of the diameter ratios for the first (black) and second (red) least stable modes. The dashed and dot-dashed lines indicate the two-dimensional critical Reynolds number based on $D$ ($Re_D = 47$) and $d$ ($\widehat {Re}_{d,cr2} = Re_{D,cr} \cdot r$), respectively. Observe that for $r=2$ and $r=4$, the values are drastically higher, $Re_d =94$ and $Re_d =188$. The light blue area below the neutral curve is the stable region in the ($Re_{D}$, $r$) space.

Therefore, although the discontinuity does not impact the onset of transition, as the region behind the large cylinder undergoes transition with a Reynolds number close to the two-dimensional cylinder, the transition of the small cylinder is notably influenced by the presence of the junction. The junction introduces three-dimensional effects that trigger an earlier transition compared with an equivalent uniform cylinder with a diameter $d$. The unstable wake behind the large cylinder globally destabilises the flow, with the transition of the flow behind the small cylinder occurring at a local Reynolds number much below the critical threshold.

3.5. Unstable eigenmodes and wake cells

As shown in figure 4, the first least-stable eigenmode develops behind the large cylinder, resembling the vortex shedding in the L cell observed in both the laminar (Dunn & Tavoularis Reference Dunn and Tavoularis2006; Morton & Yarusevych Reference Morton and Yarusevych2010; Tian et al. Reference Tian, Jiang, Pettersen and Andersson2017) and turbulent (Massaro et al. Reference Massaro, Peplinski and Schlatter2023d; Tian et al. Reference Tian, Zhu, Holmedal, Andersson, Jiang and Pettersen2023) regime. It is interesting to observe that the modal linear mechanism leading toward chaos remains discernible at a higher Reynolds number. By conducting global stability analysis at higher $Re_D$, we aim to establish if the remaining cells exhibit similar mechanisms and assess their dependence on $r$.

For $r=1.1$ and $r=1.2$, the second least-stable eigenmode localises in the wake of the small cylinder (S cell), whereas the third least-stable eigenmode, at higher Reynolds numbers, predominantly lives in the modulation region (N cell). Conversely, for $r=2$ and $4$, the roles reverse; the second least-stable eigenmode resembles the N cell and the third is in the S cell. Regardless of the order, which is likely related to the approaching uniformity for ($r\rightarrow 1$), we find that for all diameter ratios, there exists an eigenmode leading to the N cell formation through a bifurcation. This finding contrasts with the classification proposed by Lewis & Gharib (Reference Lewis and Gharib1992), who identified only two dominant modes when $r<1.25$. In their experiments on laminar flow at low Reynolds numbers across various diameter ratios, the formation of the N cell was not detected for $r<1.25$. However, the global stability analysis reveals three supercritical Hopf bifurcations occurring in three distinct cells. This result is crucial as it suggests the presence of the N cell across all $r$ regimes, whereas the dominance of the S and L cells increases as the stepped cylinder approaches uniformity.

Another noteworthy observation pertains to the geometries with $r=1.1$ and $1.2$. As discussed in the base flow characterisation, although the downwash is limited in these cases, its effect is not negligible. In addition, the vertical velocity profile exhibits symmetry with respect to the $z=0$ plane. Conversely, for $r=2$ and $4$, the downwash occurs only behind the larger cylinder, as depicted in figure 2. The downwash phenomenon was previously examined by Massaro et al. (Reference Massaro, Peplinski and Schlatter2023d) in the turbulent regime. Although the current Reynolds number is considerably lower, some parallels can be drawn. Indeed, if the downwash initiated the formation of the N cell, a symmetric downwash in the base flow would likely result in a dual modulation region. Specifically, an analogous mode to the green/blue mode in figure 6(b), mirrored across the plane $z=0$, is expected. Our global stability analysis supports this conjecture, revealing a fourth least-stable eigenmode for both $r=1.1$ and $1.2$, which largely resembles a new modulation cell predominantly behind the smaller cylinder; see figure 6(a), hereinafter named N2 cell. A cell resembling such an eigenmode has not been observed yet. In this regard, further studies at various Reynolds numbers are needed, especially for $r=1.1$ and $1.2$ which have been barely studied in the past. However, this is out of the scope of the current work.

Figure 6. Isosurfaces of the negative and positive streamwise velocity (40 % of the magnitude, normalised by its maximum) of various unstable eigenmodes resembling the wake cells. For (a) $r=1.1$, the fourth least-stable eigenmode of the direct operator resembles a modulation region, but in the wake of the small cylinder (purple/green isosurfaces). In (b–d) the second least-stable eigenmode (purple/orange isosurfaces) of the direct operator for $r=1.2$, $r=2$ and $r=4$, respectively. In addition, in (b) the green/blue isosurfaces indicate the negative and positive streamwise velocity of the third least-stable eigenmode for $r=1.2$. The black arrows indicate the flow direction.