2.1. Adaptive mesh refinement

An important extension of Nek5000 was implemented and developed by our group, namely the AMR technique (Offermans et al. Reference Offermans, Peplinski, Marin and Schlatter2020, Reference Offermans, Massaro, Peplinski and Schlatter2023; Massaro, Peplinski & Schlatter Reference Massaro, Peplinski and Schlatter2023d). The AMR technique is considered one of the relevant milestones in the future (and current) numerical simulations (Slotnick et al. Reference Slotnick, Khodadoust, Alonso, Darmofal, Gropp, Lurie and Mavriplis2014), providing significant meshing flexibility and dynamical mesh modifications during the simulation based on the instantaneous or time-averaged estimated computational error. It decreases significantly the count of computational degrees of freedom, and allows the study of flows without prior knowledge of their dynamics. The process involves the assessment of the error, with a suitable refinement criterion, and an appropriate refinement/coarsening technique. Offermans (Reference Offermans2019) introduced in Nek5000 an $h$-type refinement, where the number of degrees of freedom is modified by splitting the entire element. In this case, an oct-tree refinement is performed, replacing a single parent element with eight (in three dimensions) children (Burstedde, Wilcox & Ghattas Reference Burstedde, Wilcox and Ghattas2011). Following Kruse (Reference Kruse1997), the conforming-space/non-conforming-mesh approach is used, where the hanging nodes are not considered as real degrees of freedom, and non-conforming interfaces are treated with an interpolation operator. Massaro et al. (Reference Massaro, Peplinski and Schlatter2023d) showed that non-conforming interfaces neither affect the solution field nor introduce any further instabilities. Appendix A provides supplemental details about the error measurement and the mesh quality assessment. In particular, it is shown that the present configuration corresponds to fully resolved DNS throughout the domain.

2.2. Problem formulation

The numerical configuration for the flow around a three-dimensional stepped cylinder is presented here. First, the set of governing equations with the corresponding initial and boundary conditions is specified. Then the geometry of the problem is described, specifying the domain sizes. The current study performs DNS of the incompressible Navier–Stokes equations, which in non-dimensional form read

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

where $\boldsymbol {u}$ is the non-dimensional velocity vector, $P$ is the non-dimensional pressure, and $Re_D=UD/\nu$ is the Reynolds number based on the cylinder diameter $D$, the uniform inflow streamwise velocity component $U$, and the kinematic viscosity $\nu$. The set of equations is incomplete without proper initial and boundary conditions. In our set-up, the initial condition is a uniform distribution of the streamwise velocity with unitary value $\boldsymbol {u}(\boldsymbol {x},t_0)=(1,0,0)$, which, however, adapts quickly to a more physical solution. At the inflow, the same is set as a Dirichlet boundary condition $\boldsymbol {u}(\boldsymbol {x},t_0)=\boldsymbol {u}(\boldsymbol {x}_0,t)=(1,0,0)$. The outflow consists of natural boundary condition $(-p \boldsymbol {I} + \nu \,\boldsymbol {\nabla } \boldsymbol {u})\boldsymbol {\cdot } \boldsymbol {n} = 0$, where $\boldsymbol {n}$ is the normal vector. We focus on the junction and wake dynamics, rather than considering the end-plate effects. Thus our top/bottom boundary conditions assume an infinitely long cylinder, prescribing symmetry boundary conditions $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n} = 0$ with $(\boldsymbol {\nabla } \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {t}) \boldsymbol {\cdot } \boldsymbol {n}= 0$, where $\boldsymbol {t}$ is the tangent vector. For the front and back boundaries, mixed conditions are used, similar to the open boundary at the outflow, but prescribing zero velocity increment in non-normal directions. The stepped cylinder surface has a no-slip and impermeable wall.

The stepped cylinder presents two cylinders of different diameters joined at one extremity, with a ratio $D/d=2$ ($D=1$). An infinitely long cylinder is considered, using the boundary conditions discussed above. We have devised a configuration akin to that in Morton & Yarusevych (Reference Morton and Yarusevych2014b), in order to corroborate our findings through comparison with prior experiments. The stepped cylinder is located at the origin of the Cartesian frame oriented with the $z$-axis on the cylinder. The large and small cylinders span vertically for 10 and 12 diameters $D$, respectively. The computational box is a parallelepiped with extension $(-10D,30D)$ in $x$, $(-10D,10D)$ in $y$, and $(-10D,12D)$ in $z$, where $x,y,z$ are the streamwise, lateral and vertical (or spanwise) directions, respectively (see figure 1). To prevent any undesired side and boundary effects, the domain has been designed carefully. In the planning phase, tests were conducted to monitor closely both local and global quantities while varying the domain sizes. Furthermore, the implementation of the AMR technique allows for the usage of a substantially larger domain while ensuring DNS-like resolution. For more details, refer to Appendix A.

2.3. Consistency with experimental measures

Not many reference data of high quality exist for the present flow case. In order to verify the accuracy of our numerical configuration, local and integral measures were compared with experiments from Morton & Yarusevych (Reference Morton and Yarusevych2014b) and Roshko (Reference Roshko1961).

First, the spectral analysis of velocity signals is used to determine the vortex shedding frequency and investigate how it varies with the spanwise location. Five virtual velocity probes are placed at the same locations as in Morton & Yarusevych (Reference Morton and Yarusevych2014b). The time signals are acquired with the sampling frequency $f D/U = 100$. The Nyquist theorem requirements are definitely met since the highest frequency of the flow is expected to be the small cylinder vortex shedding at $f_s D/U \approx 0.4$. The time interval is long enough to capture the slowest dynamics exploring over approximately $200$ S-cell vortex shedding periods. The virtual probes have identical streamwise and lateral locations ($x=5D$ and $y=0.75D$), but different spanwise locations: $z_1=4D$, $z_2=-1D$, $z_3=-3D$, $z_4=-5D$ and $z_5=-9D$. At these positions, the power spectral density (PSD) of the streamwise velocity component is computed. The results are in excellent agreement with the measurements by Morton & Yarusevych (Reference Morton and Yarusevych2014b); see table 1. For each cell, a single peak is evident with a high energy content ($z_1, z_3, z_5$). At the upper (S–N) and lower (N–L) interfaces, a double peak appears. Here, the PSD highlights the interaction between two different shedding frequencies, and the twofold hills have a lower energy content; see figure 2(b).

Table 1. Strouhal numbers for different probe locations. The results from the current study and the experiments by Morton & Yarusevych (Reference Morton and Yarusevych2014b) are reported in the first and second rows, respectively.

Figure 1. Sketch of the three-dimensional stepped cylinder with planar cross-sections at $z=0$ and $x=0$. The frame of reference and the dimensions of the domain are reported. The black arrows indicate the uniform inflow with streamwise velocity $U$. The figure is adapted from Massaro et al. Reference Massaro, Peplinski and Schlatter2023c.

Figure 2. Logarithmic plot of power spectral density (PSD) of the streamwise velocity component at locations (a) $z_1=4D$, (b) $z_2=-1D$, (c) $z_3=-3D$, and (d) $z_5=-9D$. The green circles indicate the double peak at the N–L interface. The time signals are acquired with the sampling frequency $f D/U = 100$, for at least $200$ S-cell vortex shedding periods.

The frequency peak is used to calculate the Strouhal number, defined as $St=fD/U$. The frequency $f$ at the specific location and the unitary inflow velocity $U$ are considered. Once again, the agreement with the experimental measurements is excellent. Each wake region is characterised by a distinct shedding frequency and exhibits a single peak: $St_S=0.408$, $St_N=0.188$ and $St_L=0.201$ for the S-, N- and L-cells, respectively. At the S–N boundary ($z_2$), the two dominant peaks are related to the Strouhal numbers of each of the adjacent cells. Analogously, a double peak at the N–L boundary is observed. Morton & Yarusevych (Reference Morton and Yarusevych2014b) use the streamwise velocity spectra to analyse how the energy content varies in the different cells and the extension of the cell. Analogous conclusions can be drawn in this scenario: at the interface in figure 2(b), the energy content is diminished substantially compared to that within the single cells; see figures 2(a,c,d). Nevertheless, a more detailed discussion about the energy content for each cell is carried out in the POD section.

When studying the flow around a bluff body, the pressure and drag coefficients are useful integral indicators. Thus the estimated mean value and variance of the drag for each cylinder are computed. The dimensionless drag coefficient $c_D$ is defined as the component of the aerodynamic force aligned with the inflow flow (here horizontal) and normalised with $(1/2)\rho U^2 D L_z$, with $L_z$ being the vertical cylinder length. Then it is averaged in time to yield $\bar {c}_D$. The aerodynamic force components along the $y$ and $z$ directions are computed as well. For symmetry reasons, the lateral $y$ component of the coefficient $\bar {c}_y$ is expected to approach zero. In the current study, the results are $\bar {c}_{y_S}=0.0017$ and $\bar {c}_{y_L}=0.00053$, for the small and large cylinders, respectively. This assesses the good temporal convergence of our statistics for any integral quantities. When it comes to the $z$-averaged mean drag coefficients, the results are $\bar {c}_{D_S}=1.1084$ and $\bar {c}_{D_L}=0.9083$. The standard deviation of both signals is small, with value $\sigma _{c_D} \approx 0.016$. As predicted by Roshko (Reference Roshko1961), in the Reynolds number range $100< Re_D<10^5$, the drag coefficient is $c_D\approx 1$. Now the viscous and pressure contributions to the overall drag are computed. For both cylinders, the pressure term is dominant (approximately $90\,\%$) with respect to the viscous one, but the viscous component is slightly larger for the small cylinder, where $\bar {c}_{Dv} \sim 12.7\,\%$ instead of $\bar {c}_{Dv} \sim 9.6\,\%$. Because previous studies did not provide $c_D$ estimations at this specific Reynolds number and did not analyse individually drag contributions from the small and large cylinders, we undertake a comparison with the outcomes of the uniform cylinder by Roshko (Reference Roshko1961). The results are in good agreement, with integral difference less than 5 %. As expected, the junction does not affect substantially the integral drag values. However, as we will explore in the following sections, this is not true for three-dimensional time-averaged statistics, particularly in the near-wake area.

3.1. Characterisation of the junction region

Looking at the instantaneous vortical structures, we first focus on the junction region, i.e. the flat surface with sharp edges that joins the two cylinders. Aware that an unequivocal definition for a fundamental structure as the vortex does not yet exist (Cucitore, Quadrio & Baron Reference Cucitore, Quadrio and Baron1999), we rely on the $\lambda _2$ criterion by Jeong & Hussain (Reference Jeong and Hussain1995) to identify organised and coherent structures in the flow. The $\lambda _2$ definition corresponds to connected planar regions of pressure minimum, discarding the contributions of unsteady irrotational straining and viscous terms in the Navier–Stokes equations. An adequate normalisation is taken into account, and from now on each $\lambda _2$ value is normalised by $U^2/D^2$. The structures’ classification follows the convention by Massaro et al. (Reference Massaro, Peplinski and Schlatter2023c).

At $Re_D=1000$ (figure 4a), on the front of the junction, the impingement of the flow, due to the presence of the small cylinder, together with the leading edge separation, is mainly responsible for the horseshoe vortex H1 formation. By definition, the horseshoe vortex consists of a trio of interconnected vortices (Dargahi Reference Dargahi1989). These comprise a bound vortex encircling the small cylinder, along with two trailing vortices that stretch from the bound vortex into the wake area. The consequent recirculating bubble has the time-averaged pressure minimum at $(x,y,z)=(-0.354D,0D,0.046D)$. In addition, the blockage and the H1 rotation induce the incoming flow to divert laterally, where it spills over the edges of the step. Here, it rolls up into two streamwise edge vortices E, similar to the tip vortices on the free-end surface of a finite-length cylinder (Zdravkovich et al. Reference Zdravkovich, Brand, Mathew and Weston1989). These observations agree with the measurements first made by Dunn & Tavoularis (Reference Dunn and Tavoularis2006) at $Re_D=150$ and with the hydrogen bubble visualisations by Morton & Yarusevych (Reference Morton and Yarusevych2014b) at $Re_D=1000$. Recently, Tian et al. (Reference Tian, Jiang, Pettersen and Andersson2021) has pointed out some additional structures, confirmed by Massaro et al. (Reference Massaro, Peplinski and Schlatter2023c). In the laminar front part of the junction, they observe two horseshoe vortices, H3 and H4, with H4 appearing upstream of H3, and rotating in the same direction as H1. At $Re_D=1000$, the configuration is different. Figure 4(a) shows clearly that the horseshoe vortices H3 and H4 do not appear, and the leading edge separation generates only H1. This leads to a crucial difference in the E formation mechanism. As long as the cylinder shear layer remains stable, the process that generates the edge vortices agrees with the Dunn & Tavoularis (Reference Dunn and Tavoularis2006) observations at lower $Re_D=150$. The leading edge separation with the small cylinder impingement drives the recirculating vortex H1 formation that laterally deviates the incoming flow shaping the E. Furthermore, no hairpin vortices (HP) appear at $Re_D=1000$. On the other hand, a train of hairpins is clearly visible in figure 4(b), and their cyclic dynamics is described by Massaro et al. (Reference Massaro, Peplinski and Schlatter2023c). The second kind of structure reported by Tian et al. (Reference Tian, Jiang, Pettersen and Andersson2021) is named H2. Interestingly, this is present also at the current $Re_D=1000$ (figure 4). Less dynamically relevant, since it evolves at the base of the small cylinder and does not interact with other structures, this horseshoe vortex has a cross-flow width that continues to increase downstream, as for H1. The H2 structure is induced by the downward flow separation along the small cylinder wall, as confirmed by the time-averaged vertical velocity in figure 5. It is understandable that Morton & Yarusevych (Reference Morton and Yarusevych2014b) do not report these structures due to the challenges of handling hydrogen bubbles on the junction. More surprisingly, this vortex has not been classified before at lower Reynolds numbers (Vallès, Andersson & Jenssen Reference Vallès, Andersson and Jenssen2002; Dunn & Tavoularis Reference Dunn and Tavoularis2006; Tian et al. Reference Tian, Jiang, Pettersen and Andersson2020). Nevertheless, the downward flow separation along the small cylinder is not associated with the turbulent character of the flow, and H2 structures are expected to be observed at lower $Re$ as well. Simulating all three different flow regimes, and guaranteeing a well-resolved mesh around the junction with the AMR, the presence of H2 is confirmed. Based on the overall view, we assert that this is a common feature for all three flow regimes.

Figure 4. Top view of instantaneous $\lambda _2$ structures around the junction: isocontours of $\lambda _2$ coloured with the streamwise vorticity for (a) $Re_D=1000$ ($\lambda _2 =-4 U^2/D^2$) and (b) $Re_D=5000$ ($\lambda _2 =-40 U^2/D^2$). The horseshoe and hairpin vortices are indicated (H1, H2, H3, H4 and HP), along with the edge (E) and ring (R) vortical structures. In the circle, the E vortex rolling up around the edge and splitting is highlighted.

Figure 5. Time-averaged vertical velocity component in the $xz$ ($y=0$) and $xy$ ($z=0$) planes at(a) $Re_D=1000$ and (b) $Re_D=5000$.

All these vortices are stable and with only small fluctuations in time, as they develop in the laminar/transitional part of the flow. On the contrary, the back structures are weaker and more difficult to identify since they are placed in the turbulent wake of the small cylinder. This makes them almost indistinguishable from the small turbulent eddies. However, the heads of the two streamwise base vortices that begin to develop on the back part of the junction are identified; see also Massaro et al. (Reference Massaro, Peplinski and Schlatter2022). Furthermore, the $\lambda _2$ structure named ring R vortex points out the recirculation occurring at the trailing edge and is reported here for the first time. This structure is weaker than H1 and with a different formation mechanism: the E vortices roll up around the edges and split before diverging laterally; see figure 4(a). The low-pressure region behind the cylinder stabilises R on the back, whereas the strong downstream advection interrupts the ring connection laterally. At $Re_D=5000$, the ring vortex is not observed (independently on the $\lambda _2$ threshold) as the hairpins interaction modifies the flow topology substantially.

The analysis of time-averaged statistical quantities supports the description of the main junction mechanisms. In particular, figure 5 shows the vertical velocity component $\bar {W}$. On the forepart of the $xz$ plane, the leading edge separation at $x=-0.5D$ is in red, and the downward flow separation along the small cylinder wall is in blue. In the $xy$ plane, a symmetric behaviour is observed for positive and negative $y$: an upward flow induced by H1 surrounded by the downward flow separation along the small cylinder and the edge H2 induced flow. A larger number of structures at $Re_D=5000$ leads to a more complicated statistical footprint. This is always the case when an obstacle interacts with a wall, namely the small cylinder with the junction. Nevertheless, it is remarkable that these interactions are observed on such a narrow surface with only $0.25D$ radial widths; see figure 5(b). The $xy$ plane confirms the structures’ evolution detected above: the H3 appearance only at higher Reynolds numbers, and the E formation mechanism being $Re_D$-related.

3.2. Characterisation of the wake

As the Reynolds number increases, the wake becomes more intricate, and numerous turbulent eddies emerge in a range that varies from the cylinder diameter to the Kolmogorov scale. The increasing complexity is clearly visible in figure 6. At $Re_D=1000$, despite the turbulent wake behaviour, the vertical vortices with $z$-oriented axis are still visible. Superposed, smaller streamwise spaghetti-like structures are observed; see figure 6(a). In laminar vortex shedding, the $z$ vortices are distinguishable, making the vortex dislocations simple to track (Tian et al. Reference Tian, Jiang, Pettersen and Andersson2020). This is definitely not the case for a turbulent wake. Nevertheless, at the S–N boundary, the occurrence of the dislocation is still visible since the boundary is characterised by two largely different frequencies and is thus stable. Figure 6(a) shows the breaking of the initial in-phase vortex shed by the cylinder, and the consequent dislocation. The oblique shedding occurs in the modulation region with the dislocation at its boundary (Massaro et al. Reference Massaro, Peplinski and Schlatter2022). Contrary to the upper interface, the N–L layer is unstable, and dislocation detection is not possible. At the highest Reynolds number, the situation worsens even further; see figure 6(b). Massaro et al. (Reference Massaro, Peplinski and Schlatter2023c) characterise the wake of the latter case: secondary spanwise vortices appearing behind the small cylinder (caused by a Kelvin–Helmholtz instability), spanwise vortices stack organisation with the appearance of small- and large-scale braids (see figure 5 in Massaro et al. Reference Massaro, Peplinski and Schlatter2023c). Previous works developed different techniques to study the wake dynamics, e.g. for obtaining the phase and the shed position of vortices (Tian et al. Reference Tian, Jiang, Pettersen and Andersson2020). Due to the inherent complexity of turbulence, these methods, valid at lower Reynolds numbers, no longer hold. In the next section, a different approach is proposed, using the modal analysis to identify the energy content and the extension of the cells, and to study the large-scale vortex connections.

Figure 6. Instantaneous $\lambda _2$ structures coloured with the vertical vorticity component in the wake at (a) $Re_D=1000$ ($\lambda _2 =-4 U^2/D^2$) and (b) $Re_D=5000$ ($\lambda _2 =-40 U^2/D^2$), where (b) is adapted from Massaro et al. (Reference Massaro, Peplinski and Schlatter2023c).

Due to the turbulent nature of the wake, analysing its three-dimensional time-averaged statistics can guide us in distinguishing between the subcritical ($Re_D=1000$) and fully turbulent ($Re_D=5000$) regimes. As expected, the wakes of both the small and large cylinders exhibit similarities to the uniform cylinder wakes far away from the junction. Nevertheless, significant differences are noticeable in the region below the junction, which extends to a considerable distance of 4–5 times the diameter of the large cylinder. The time-averaged streamwise velocity component in the $xz$ planes ($y=0$) points out distinctly different extensions of the recirculation; see figure 7(a). At $Re_D=1000$, the negative $\bar {U}$ region extends for $0.63D$ and $2.43D$ behind the small and large cylinders, respectively. In the case of the larger cylinder, the recirculation gets stronger inside the N-cell. Qualitatively, this is clear from the more intense blue colour in figure 7(a). This is a meaningful characteristic of the modulation area only, and the connection with the downwash mechanism is discussed below. Quantitatively, albeit the trend of the streamwise velocity $\bar {U}$ does not vary for the two Reynolds numbers, the backflow phenomenon is more intense at $Re_D=5000$. The peak reaches 40 % of the inflow velocity, with the minimum occurring at a lower vertical location; see figure 7(b). Likewise, figure 7(c) shows that at higher Reynolds numbers, the mean vertical velocity $\bar {W}$ exhibits significantly more negative values. Despite the differences in the wake characteristics, the same trend is observed at distance approximately $5D$ from the cylinder (figure 7). Specifically, the wake profile begins to resemble that of a uniform cylinder at this point, regardless of the Reynolds number. This suggests that the extension of the N-cell structure, which is responsible for the oblique periodic shedding of vortices, is not affected significantly by changes in the Reynolds number and the transition of the cylinder shear layer. This is also confirmed by the POD results in the next section. Eventually, the variation in the vortex formation length is worth noting and deserves further discussion as it plays a crucial role in the dynamics of the wake.

Figure 7. (a) Time-averaged streamwise velocity component in the $xz$ plane at $y=0$ for $Re_D=1000$. The black dashed lines indicate the locations where the data are extracted for (b,c), corresponding to $z=0D$ (black), $z=-2D$ (blue), $z=-4D$ (green) and $z=-8D$ (red). The time-averaged (b) streamwise and (c) vertical velocity components are plotted at different spanwise locations. Both velocity components are normalised by the inflow velocity $U$. The solid and dashed lines refer to $Re_D=1000$ and $Re_D=5000$, respectively.

The spanwise variation of the vortex formation length refers to the change in the distance between the consecutive vortex shedding points in the wake, as we move along the spanwise direction. This variation can lead to the development of different modes of vortex shedding, which can affect the wake structure and its associated aerodynamic properties. Therefore, understanding the spanwise variation of the vortex formation length is essential to gain insight into the wake dynamics and its impact on the flow field. In the current study, the location of the maximum streamwise velocity fluctuation $\overline {uu}$ is used to define the vortex formation length (Norberg Reference Norberg1998; Govardhan & Williamson Reference Govardhan and Williamson2001). The vortex formation lengths behind the small cylinder for two different $Re_D$ values are observed to match well, with a noticeable reduction inside the N-cell region; see figure 8. Specifically, decreases of approximately 20 % and 15 % are measured at $Re_D=1000$ and $Re_D=5000$, respectively. This contraction has a significant connection with the downwash phenomenon, as discussed below. Indeed, the reduction in vortex formation length leads to a decrease in the distance between the vortices shed by the cylinder. This, in turn, creates a larger wake deficit, which has a substantial impact on the aerodynamic forces experienced by the cylinder.

Figure 8. Spanwise variation of the vortex formation length for Reynolds numbers $Re_D=1000$ and $5000$.

In summary, the two turbulent regimes at $Re_D=1000$ and $Re_D=5000$ present similar characteristics far away from the junction, where the wake resembles the uniform cylinder shape. The reduction of the vortex formation length inside the N-cell is also consistent. However, the ratio between the frequencies of N- and L-cells tends to diverge with the Reynolds number (Massaro et al. Reference Massaro, Peplinski and Schlatter2023c). The localised and stronger vertical velocity flux inside the modulation region at $Re_D=5000$ (with more intense $\overline {ww}$; see figure 9 in Massaro et al. Reference Massaro, Peplinski and Schlatter2023c) is the consequence of a less energetic downwash phenomenon behind the junction. Therefore, the instability of the cylinder shear layer results in an increase of the energy content for the modulation cell, in agreement with the trend of the POD spectrum (see figure 10 in Massaro et al. Reference Massaro, Peplinski and Schlatter2023c). As discussed below, the POD identifies the three cells as the most energetic structure, in agreement with the statistical analysis for both $Re_D$ values. Nevertheless, a remarkable difference in the downwash dynamics is observed between the stable and unstable shear layer cases.

4.1. Theoretical formulation

Given a set of $m$ snapshots $\boldsymbol {u}_j$, with $j=0,\ldots,m-1$, these are assembled in the snapshot matrix $\boldsymbol{\mathsf{U}}_m$:

(4.1)\begin{equation} \boldsymbol{\mathsf{U}}_m = [\boldsymbol{u}_0, \boldsymbol{u}_1,\ldots,\boldsymbol{u}_{m-1}] \in \mathbb{R}^{n \times m}, \end{equation}

where $n$ is the total number of grid points multiplied by the number of velocity components. The modal analysis allows us to split the system into space- and time-dependent parts. In this regard, the POD is employed (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993). The POD performs the modal decomposition of the nonlinear flow minimising the residual energy between the snapshots and their reduced linear representation. They are constructed to be the orthogonal basis that represents the optimal projection, in energy terms, of the most important dynamical features. It was first introduced in the fluid mechanics context by Lumley (Reference Lumley1970), and is usually named classical POD. Here, the snapshot POD is used, developed by Sirovich (Reference Sirovich1987) to handle the analysis of three-dimensional flows with a large set of grid points. The notation adopted in the following derivation follows Sarmast et al. (Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanell, Sørensen and Henningson2014).

When assembling the snapshots matrix, the collected snapshots are sampled according to the Nyquist criterion ($\Delta t =0.25 D/U$) and collected for an interval that is sufficiently long to embed the slowest frequency. As indicated in (4.1), the set of $m$ snapshots $\boldsymbol {u}_j$ is rearranged in the snapshot matrix $\boldsymbol{\mathsf{U}}_m$. The orthogonal basis function $\boldsymbol {X}$ and the corresponding time coefficients $\boldsymbol {T}$ are computed via the singular value decomposition of the snapshot matrix as

(4.2)\begin{equation} \boldsymbol{\mathsf{U}}_m = \boldsymbol{X} \boldsymbol{T} = \boldsymbol{X} \boldsymbol{\varSigma} \boldsymbol{\mathsf{W}}, \end{equation}

where $\boldsymbol {T} = \boldsymbol {\varSigma } \boldsymbol{\mathsf{W}}$, with $\boldsymbol{\mathsf{W}}$ being the matrix of the right singular vectors, and $\boldsymbol {\varSigma } \in \mathbb {R}^{n\times m}$ the matrix of singular values of $\boldsymbol{\mathsf{U}}_m$ (i.e. the energies). To guarantee the unit energy of each mode, the proper normalisation has to be taken into account: $\| \boldsymbol {X} \|_{{M}} = \boldsymbol{\mathsf{I}}$ and $\| \boldsymbol{\mathsf{W}} \|_{{N}} = \boldsymbol{\mathsf{I}}$, where $\boldsymbol{\mathsf{M}} \in \mathbb {R}^{n\times n}$ and $\boldsymbol{\mathsf{N}}\in \mathbb {R}^{n\times n}$ are the mass and the temporal weight matrix, respectively, and $\boldsymbol{\mathsf{I}}$ is the identity matrix. If the snapshots are collected with equidistant time intervals, then $\boldsymbol{\mathsf{N}}=\boldsymbol{\mathsf{I}}$. The above normalisation is obtained by considering $\boldsymbol{\mathsf{M}}^{1/2} \boldsymbol{\mathsf{U}}_m \boldsymbol{\mathsf{N}}^{1/2}$ and computing

(4.3) \begin{equation} \boldsymbol{\mathsf{M}}^{1/2} \boldsymbol{\mathsf{U}}_m \boldsymbol{\mathsf{N}}^{1/2} = \hat{\boldsymbol{\mathsf{X}}} \boldsymbol{\varSigma} \hat{\boldsymbol{\mathsf{W}}}, \end{equation}

where $\hat {\boldsymbol{\mathsf{X}}}$ and $\hat {\boldsymbol{\mathsf{W}}}$ are unitary matrices: $\hat {\boldsymbol{\mathsf{X}}}^{\rm T} \hat {\boldsymbol{\mathsf{X}}} =\boldsymbol{\mathsf{I}}$ and $\hat {\boldsymbol{\mathsf{W}}}^{\rm T} \hat {\boldsymbol{\mathsf{W}}} =\boldsymbol{\mathsf{I}}$.

Eventually, the modes are reconstructed as $\boldsymbol {X} = \boldsymbol{\mathsf{M}}^{-1/2} \hat {\boldsymbol{\mathsf{X}}}$, with unit energy and orthogonal to each other, while the time coefficients are $\boldsymbol {T} = \boldsymbol {\varSigma } \boldsymbol{\mathsf{W}}^{\rm T}$, where the energies are given by the diagonal matrix $\boldsymbol {\varSigma }$.

4.2. Convergence criterion

Before discussing the modal analysis results, the POD convergence is assessed and the set of snapshots is guaranteed to be adequately large. To this end, first the energy spectrum for progressively larger $m$ is computed, and then the correlation coefficient is estimated for two different snapshot sets. To improve the convergence of the decomposition, the $y$-symmetry of the stepped cylinder is exploited, resulting in a statistical symmetry for the flow (Berkooz et al. Reference Berkooz, Holmes and Lumley1993). In figure 9, the spectrum is computed for progressively larger data sets, from $m=1200$ to $m=2400$ snapshots. The last three cases (with $m=1800,2000,2400$) show excellent agreement for the first four modes. Starting with modes 5 and 6, the blue line ($m=1200$) deviates from the others. Contrary to these, the azure ($m=2000$) and black ($m=2400$) curves nearly overlap, especially for $j \leqslant 200$. As the POD orders the modes based on their fluctuating energy levels, we observe that the spectra of the first $200$ energetically relevant modes do not show significant deviations. In addition, good agreement between mode 0 (which corresponds to the expected mean value according to our notation) and the correspondent mean flow (estimated via statistical analysis in Massaro et al. 2023, personal communication) is appreciated.

Figure 9. The POD spectrum with the percentage energy contribution: (a) first $500$ and (b) first $20$ modes are shown. Different numbers of snapshots $m$ are considered. The orange dashed ellipses and solid circle indicate the three most energetic travelling modes and the downwash mode (single point), respectively.

To corroborate previous qualitative observations, a quantitative estimation of the POD convergence is performed by dividing the data sets into two blocks ($A$ and $B$), with a minimum of $10$ % overlapping. On each subset, the POD is carried out, and the correlation between the sets is computed. The criterion is adapted from previous works (Lesshafft et al. Reference Lesshafft, Semeraro, Jaunet, Cavalieri and Jordan2019; Abreu et al. Reference Abreu, Tanarro, Cavalieri, Schlatter, Vinuesa, Hanifi and Henningson2021) where the SPOD is used. Thus some adjustments are required. Let us consider the pair of modes $i$-$j$: they represent a unique travelling mode where the mode $j$ is the same mode $i$ with a 90-degree phase shift. For the two sets of snapshots $A$ and $B$, the travelling modes are $\varPhi _{ij}$ and $\hat {\varPhi }_{ij}$, respectively:

(4.4)\begin{gather} \varPhi_{ij} = {\rm Re} \left\{ \phi_i + \mathrm{i} \phi_j \right\}, \end{gather}

(4.5)\begin{gather}\hat{\varPhi}_{ij} = {\rm Re} \left\{ (\phi_i + \mathrm{i} \phi_j)\,{\rm e}^{\mathrm{i} \theta} \right\}, \end{gather}

where $\mathrm {i}$ is the imaginary unit, and the indices $i$ and $j$ refer to the POD mode number. The angle $\theta$ is needed to take into account the possible phase variation between $\varPhi$ and $\hat {\varPhi }$, as explained below. At this point, the correlation coefficient is introduced for a travelling mode $p$, made of two modes $i$ and $j$:

(4.6)\begin{equation} \rho_{p} = \frac{ \left \langle \varPhi_{ij} \hat{\varPhi}_{ij} \right\rangle }{\| \varPhi_{ij} \| \| \hat{\varPhi}_{ij} \|}, \end{equation}

e.g. the first travelling mode has $p=1$ with $i=1$ and $j=2$. Since the travelling modes $\varPhi _{ij}$ and $\hat {\varPhi }_{ij}$ arise from different snapshots, they have many different orientations in the complex plane, and the possible phase shift has to be taken into account. Thus $\hat {\varPhi }$ is multiplied by ${\rm e}^{\mathrm {i} \theta }$, with $\theta$ varying between 0 and $2{\rm \pi}$. When they are aligned, their correlation is going to be maximum. The phase difference varies among the travelling modes as reported in table 2.

Table 2. Correlation coefficient $\rho _p$ for all the main pairs $p$ of POD modes $i$-$j$. In the last column, the phase shift $\theta$ for the travelling mode $p$ is computed from two sets of snapshots.

When the correlation coefficient is close to unity, the POD convergence is achieved. For the first two pairs, the correlation is higher than $95\,\%$. Interestingly, the second pair has better convergence than the first. This can be explained by the fact the travelling mode 2-3 is related to the small cylinder, which has a vortex shedding frequency two times larger than the other frequencies. Consequently, for the same period $T$, double the number of vortices are shed. Considering the modulation region modes (5-6), even they show good convergence: $\rho _{p}>0.87$. This criterion does not apply to single mode 7, which is discussed below. Even the next pair, $8$-$9$, is highly correlated, indicating that the first ten modes are converged adequately ($\rho _{p}>0.87$) given the collected number of snapshots.

4.3. The most energetic structures

We begin by presenting the most relevant structures extracted via POD of the flow. They correspond to three main travelling modes in the streamwise direction; see figure 10. Each of them consists of a pair of modes with a $90$-degree phase shift, which characterises the L-, S- and N-cells, respectively. They are all made symmetric with respect to the $y$-axis, reflecting the inherent $y$-symmetry of the flow.

Figure 10. Isosurfaces of the velocity magnitude $0.1U$ for the (a) first, (b) third and (c) fifth POD modes at arbitrary phase. The second, fourth and sixth modes correspond to the phase-shifted versions of the shown modes. Dashed lines indicate the corresponding cells, and the solid ellipse highlights the S vortical structures deflection. Dotted lines point out the N-cell contraction in the streamwise direction.

The first pair of modes (1 and 2) corresponds to the large cylinder structures in the L-cell. Their energy contribution is estimated by computing the fraction $\lambda _j/\sum _{j>0}\lambda _j$, where $j$ is the mode index, and $\lambda$ is the singular value for mode $j$. They account for approximately ${\sim }7.5\,\%$ of the total (fluctuating) energy each. This means that the coherent structures behind the large cylinder (and below the modulation region) give an overall energy contribution of ${\sim }15\,\%$. The S-cell of the small cylinder is characterised by the second travelling mode 3-4, overall as counting for ${\sim }11\,\%$. It is noticeable that the small cylinder dynamics (with $Re_d = Re_D/2$) still has a significant energy contribution. The modulation region pair is formed by modes $5$ and $6$, whose contribution is half of the previous couple, i.e. ${\sim }2.8\,\%$, in total ${\sim }5.6\,\%$.

Morton & Yarusevych (Reference Morton and Yarusevych2014b) measure the extension of the cells using hydrogen bubble visualisations. The qualitative description of the flow carried out in § 3 through $\lambda _2$ structures corroborates their observation. In addition, the extracted POD modes allow us to define the spatial extension of each cell in a quantitative way for the turbulent regime. Indeed, each of the most energetic POD modes corresponds to a specific cell as it holds no value beyond the regions outlined below. The S region spreads behind the small cylinder and slightly below the $xy$ centre plane, at $z\approx - 0.5D$; see figure 10. On the other side, the large cylinder modes (1-2) extend for $z<-5D$. Regarding the N-cell, Dunn & Tavoularis (Reference Dunn and Tavoularis2006) observed a non-constant vertical extension in the streamwise direction. In particular, at the laminar Reynolds numbers $Re_D=150$ and $152$, this is also reported by Morton & Yarusevych (Reference Morton and Yarusevych2010). It is worth noting that such a contraction is persistent in the turbulent regime and depicted by the fifth mode in figure 10(c). The spatially coherent structure of the N-cell has a width that varies with streamwise location. In particular, the upper and lower envelope of the cell contract proportionally to $y_{up} \sim -0.1x$ and $y_{low} \sim 0.15x$, respectively. Moreover, the modulation region develops only behind the large cylinder $-0.5D< z<-4.5D$. At first glance, this may appear counterintuitive since one would expect that the intermediate region is straddling the junction, rather than being only behind one of the cylinders. A possible explanation and its connection with the downwash mechanism are discussed below. At the same time, between cells, the vortex dislocation mechanism by Williamson (Reference Williamson1989) occurs: the vortices in the cell with lower frequency are induced downstream by the high-frequency cell; see N–L interaction. That is not the case for the N–L boundary since their frequencies are close (table 1). Differently from the S–N stable interface where the POD modes have a net separation, the unstable layer N–L generates an overlap at $z\approx -5D$ between the modes 1-2 and 5-6. The spatial extension of the three cells is also confirmed by inspecting the vertical variation of the vortex shedding frequency from the nonlinear field.

In previous studies, the S-, N- and L-cells and their locations are distinguished by frequencies and convective streamwise velocities. Here, they are categorised by the energy content too. Many studies (Perrin et al. Reference Perrin, Braza, Cid, Cazin, Barthet, Sevrain, Mockett and Thiele2007; Siegel et al. Reference Siegel, Seidel, Fagley, Luchtenburg, Cohen and Mclaughlin2008; Kourentis & Konstantinidis Reference Kourentis and Konstantinidis2012) consider the distribution of energy for the uniform cylinder, and observe the energy of the first pair being approximately ${\sim }40\,\%$. In the stepped cylinder case, a similar amount is (not equally) distributed among the three different cells. Morton & Yarusevych (Reference Morton and Yarusevych2014a) use the particle image velocimetry technique to study the flow around a dual-stepped cylinder focusing on the near-wake region. They perform a POD of planar velocity data for the cases $0.2< L_j/D < 3$, $D/d=2$ at $Re_D=2100$ and $Re_D=1050$, where $L_j$ is the extension of the junction region. The evident geometrical differences and two-dimensional (instead of three-dimensional) modal analysis do not allow an adequate comparison (Hufnagel et al. Reference Hufnagel, Canton, Örlü, Marin, Merzari and Schlatter2018). Our studies confirm how the energy distribution varies with respect to the uniform case: a considerable energy reduction of the first POD travelling mode. The discontinuity creates three separated pairs of modes that correspond to the first pair of a uniform cylinder: two energetic regions (S and L), interrupted by a low-energy cell (the modulation N). Even when the Reynolds numbers are higher, three distinct POD modes are shown, corresponding to the L-, S- and N-cells at $Re_D=5000$. However, the energetic contributions of S- and N-cells become comparable (Massaro et al. Reference Massaro, Peplinski and Schlatter2023c).

The analysis of time coefficients $\boldsymbol {T}$ identifies the shedding frequency of each travelling mode, i.e. of each cell. Moreover, the spectral analysis performed via local velocity measures, and presented in § 3, is compared with the POD results. The PSD analysis of the time coefficients of the first three modes, computed as Welch's PSD, is presented in figure 11. Despite the POD modes being orthogonal only in space, and each mode containing multiple frequencies, the agreement between Welch's peak in figure 11 and table 1 is remarkable. Indeed, the first three pairs of modes, i.e. the travelling modes of the L-, S- and N-wake regions, have Strouhal numbers $St_{p=1}=0.2011$, $St_{p=2}=0.4062$ and $St_{p=3}=0.1877$, respectively. As each flow region is dominated by a single frequency, compared to alternative modal decomposition, e.g. dynamic mode decomposition (DMD), POD is equally able to locate these regions and reveal their dynamics (Rowley et al. Reference Rowley, Mezic, Bagheri, Schlatter and Henningson2009; Malm, Bagheri & Schlatter Reference Malm, Bagheri and Schlatter2011). We would not expect any better mode identification using either DMD or spectral POD in this particular case.

Figure 11. Welch's PSD estimate for the time coefficients $t_i$ of the most energetic three-dimensional POD modes: mode 1, mode 3 and mode 5. The same estimates are obtained for the corresponding paired modes 2, 4 and 6.

4.4. The downwash mode

The POD of the flow not only portrays the energy content in each cell but also improves our N-cell understanding. We comment here about two isolated modes visible in the POD spectrum. In figure 9(b), we notice that the seventh mode to be a single point in the spectrum, as opposed to all previous ones that pair. The energy content of mode 7 is significantly lower than the top ones, but still not negligible, ${\sim }1.4\,\%$. For larger $j$, a further pair of modes (8-9) is identified, and eventually another single point (mode 10). The modes for $j>10$ are energetically irrelevant and do not add any more to our analysis. In our effort to explain these single modes, the temporal dynamics of modes 7 and 10 are analysed, and the PSD via Welch's method is estimated. The PSD peaks are located at the same frequency: $f_{7}=f_{10}=0.0135$. All this suggests that we deal with a unique travelling mode, as discussed in the next subsection. Even more interestingly, the frequency $f_{7}=f_{10}$ is the exact difference between the L- and N-cell shedding frequencies computed via POD: $f_b = f_{p=1}-f_{p=3}=0.2012-0.1877=0.0135$, i.e. the beating frequency of the modulation region. Previously, Morton & Yarusevych (Reference Morton and Yarusevych2014b) measured a cyclic behaviour of the N-cell with a beating frequency $\tilde {f}_b=0.009$, which is fairly close to $f_b$. A little deviation is acceptable since they computed the PSD through a single point probe (a hot wire) and for a slightly different Reynolds number ($Re_D=1050$). In addition, they have a broad peak estimation in the spectrum. With regard to the downwash, they used flow visualisations to identify its location, and they measured its frequency with a velocity probe at $x=1D$, $y=0$, $z=-1.5D$. The result is noticeable: they observe the same value as the N–L beating $\tilde {f}=0.009$. The modal analysis is thus used to identify the dynamical mode related to the N-cell beating (mode 7), rather than relying only on local measurements. By visualising mode 7 (and 10), they exhibit two downward-deflecting velocity structures. Whereas $x$- and $y$-velocity components are negligible, the $z$-velocity has a spatial distribution that portrays the downward flow motion taking place at the back of the junction, into the N-cell (figure 12). Its location behind the step discontinuity (in the large cylinder wake) extends for $z \in (-2.3D,-0.46D)$, $y \in (-0.60D,0.60D)$ and $x\in (1.09D,2.58D)$, in agreement with the statistical measures discussed below.

Figure 12. The downwash mechanism behind the junction of a stepped cylinder. (a) Isosurface $W/U=-0.04$ of the vertical velocity component $W$ of the downwash mode 7 in blue. The edge vortex E at $y<0$ is reported (instantaneous $\lambda _2=-2$ vortical structure) and coloured with the streamwise vorticity (red is positive; see figure 4). At $y=0.45D$, in the background, the time-averaged vertical component of the vorticity (downwards in red) is shown. The edge vortices induce the downward flow motion, as indicated by the black arrows. This flow feature is depicted by the POD modes 7 and 10. The $W$ contour plots in (b) the $xz$ plane at $y=-0.45D$ and (c) the $yz$ plane at $x=2D$, of the POD mode 7.

Thus the POD technique allows us to depict the dynamical mode associated with the downwash mechanism. The downwash is the phenomenon occurring behind the free end of a cylinder, where the low pressure leads to a downward motion of the flow. From the mass conservation law, it follows that the flow displaced by the cylinder needs to fill the low-pressure region created behind it. However, even the flow above the free end of a finite-height circular cylinder (a geometry much simpler than the stepped cylinder) is far from being completely understood. Despite several models, experimental campaigns and numerical simulations, we lack complete knowledge of the flow topology and dynamics (Sumner Reference Sumner2013). Most of the studies are in agreement on some basic elements of the average flow field on the free-end surface, including, for example, the emergence of a persistent recirculation zone with distinct mean cross-stream structures such as an ‘arch’ or ‘mushroom’ vortex. This structure, together with the near-wake cross-stream vortex, has been identified as the main thing responsible for the downwash phenomenon. However, the arch vortex does not appear on the flat surface of the stepped cylinder, due to the presence of the small cylinder installed in the middle of the junction; see figure 4. On the other hand, we observe the edge vortices developing at the sides of the junction, analogously to the side-tip or tip/trailing streamwise vortices of the uniform cylinder. In the past, it has been debated about the (weak) influence of these structures on the downwash itself (Park & Lee Reference Park and Lee2000). Our findings support previous theories and results (Johnston & Wilson Reference Johnston and Wilson1996; Palau-Salvador et al. Reference Palau-Salvador, Stoesser, Frohlich, Kappler and Rodi2010). Indeed, at $Re_D=1000$ (subcritical regime), the two stable edge vortices generate a remarkable downward motion, depicted by POD modes 7-10 (a few percentages of the total fluctuating energy). As sketched in figure 12, the negative vertical velocity isosurface clearly highlights the downward flow in correspondence with the two edge vortices. As previously stated, the latter structures are stable, as related to the laminar shear layer. The time-averaged vertical vorticity confirms the statistical tendency of the flow to deflect behind the E vortices (figure 12a). Increasing the Reynolds number, the shear layer gets unstable and these structures get weaker, being surrounded by a train of hairpins (Massaro et al. Reference Massaro, Peplinski and Schlatter2023c). Consequently, the modal analysis at $Re_D=5000$ does not identify such a flow feature among the most energetic modes.

Furthermore, there exists a temporal relationship between the POD modes and the modulation region's frequency, where $f_{7}$ and $f_{10}$ coincide with the N beating frequency. The ‘downwash mode’, is thus defined by POD modes 7-10 and has an overall contribution to the fluctuating energy of approximately 2.2 %. Remarkably, by means of the modal analysis, we manage to isolate its dynamics and show the spatial and temporal connection with the modulation cell. The relation between the seventh and tenth modes is explained further in the next subsection, where the analytical superposition of two travelling waves is used to support our results.

4.5. Superposition of travelling modes

During the spectral analysis, we identify the presence of a distinct mode, mode 7, and investigate its potential connection with mode 10. Our current objective is to establish a robust link between them. Furthermore, we also consider the influence of the travelling mode 8-9.

To provide a better understanding, we will use the example of a one-dimensional travelling wave $T_i=\sin (k_ix-\omega _i t)$, with constant wave number $k_i$ and angular frequency $\omega _i$. Each travelling wave is characterised by a frequency $f$ ($\omega _i = 2 {\rm \pi}f_i$), a wavelength $\lambda$ ($k_i = 2 {\rm \pi}/\lambda _i$), and a speed $v=\lambda f$. Defining $\alpha =k_1x-\omega _1 t$ and $\beta =k_2x-\omega _2 t$, the superposition of two travelling waves is

(4.7)\begin{align} T_1 + T_2 & = \sin(\alpha)+\sin(\beta) \nonumber\\ & = 2\sin \left( \frac{\alpha + \beta}{2} \right) \cos \left( \frac{\alpha - \beta}{2} \right) \nonumber\\ & = 2\sin \left( \frac{k_1+k_2}{2}\,x - \frac{\omega_1+\omega_2}{2}\,t \right) \cos \left( \frac{k_1-k_2}{2}\,x - \frac{\omega_1-\omega_2}{2}\,t \right) \nonumber\\ & = 2\sin \left( \bar{k} x - \bar{\omega} t \right) \cos \left( \hat{k} x - \hat{\omega} t \right), \end{align}

where $\bar {k} = (k_1+k_2)/2$, $\bar {\omega } = (\omega _1+\omega _2)/2$, $\hat {k} = (k_1-k_2)/2$ and $\hat {\omega } = (\omega _1-\omega _2)/2$. The result is a product of two travelling waves: the sine part oscillates with the average frequency $\bar {f}$, and the cosine wave oscillates with the difference in frequencies $\hat {f}$. The latter controls the amplitude envelope of the wave, causing the beating phenomenon with the beat frequency $f_b$, which is twice $\hat {f}$.

Keeping this model as the basis, we now turn our attention to the modes obtained through modal analysis for the stepped cylinder, focusing specifically on the L- and N-cells. The temporal frequencies are obtained directly from the POD spectrum, $f_L=0.201$ and $f_N=0.188$, i.e. angular velocities $\omega _L=1.263$ and $\omega _N=1.179$. The approach differs for the spatial wavelengths. These modes travel only in the streamwise direction with a given streamwise wavenumber. The wavelengths $\lambda$ can be estimated by inspecting the extracted POD modes at fixed $y$ and $z$ locations to obtain a simplified one-dimensional travelling wave. The $xy$ planes at $z=-8D$ and $-3D$ are considered, i.e. in the cores of L- and N-cells, respectively, selecting the $y=0$ location. Finally, a constant wavelength in $x$ is assumed. In fact, a slight variation in the distance between consecutive peaks is observed in the near-wake region for $x<3D$. For this reason, $x>3D$ is considered to estimate the constant wavenumbers: $\lambda _L=2.02$, $k_L=3.109$, and $\lambda _N=2.30$, $k_N=2.730$. Despite our approximations, the simple model (4.7) improves our understanding of the superposition between the N and L travelling modes. In relation to $T_L+T_N$, the formula (4.7) gives the product of two waves with $\bar {\lambda } = 2.15$, $\bar {f} = 0.19$, $\hat {\lambda } = 33.14$ and $\hat {f} =0.0069$. The beating frequency is $f_b=2\hat {f}=0.0136$, which corresponds to the frequency peak observed in the POD spectrum of modes 7-10 ($\,f_{p=1}-f_{p=3}=0.2012-0.1877=0.0135$). The frequency of the downwash mode aligns precisely with the expected beating frequency of the N-cell, thereby establishing a clear connection between the downwash observed behind the junction and the beating oscillation of the N-cell. This finding offers compelling proof that the downwash phenomenon constitutes an essential aspect of the flow dynamics at the junction, connected intricately to the behaviour of the N-cell. In addition, the beating nodes are not fixed but rather evolve in time with wave speed $v=\lambda f$.

Using the formula (4.7), as explained above, the estimated $x$ wavelength of the beating mode is at least one order of magnitude bigger than all the others, $\hat {\lambda }=33D$. Eventually, the $T_L+T_N$ model predicts also the travelling mode 8-9 with high accuracy: the sinusoidal wave modulated by the beating has the averaged temporal frequency and streamwise wavenumber $\bar {\lambda } = 2.15$ and $\bar {f} = 0.19$. These results are in excellent agreement with the POD modes 8-9 ($\lambda =2.14$ and $f=0.1875$). Differently from the L–N interaction, the superposition of S–N travelling modes does not turn out to be relevant energetically. Indeed, all POD modes larger than 10 contribute to less than $0.5\,\%$. As observed by Tian et al. (Reference Tian, Jiang, Pettersen and Andersson2020), the phase difference accumulation (led by the slightly different frequencies) along with the varying convective velocities generates an unstable interface at the N–L boundary, differently from the stable S–N cell boundary.

4.6. Statistical footprint of the downwash

In § 3, the wake dynamics is also characterised using three-dimensional statistics. Now, concentrating on the subcritical regime, we search for any statistical traces of the downwash in the region outlined by the modal analysis. No previous work analyses statistical trends in this cell for the stepped cylinder; just a few studies describe the cell behind the free end of the three-dimensional uniform cylinder in terms of pressure and velocity variations. In the limit of a significant discontinuity $r=D/d \gg 1$ (e.g. $r=2$), the flow behind the junction is expected to resemble the free-end surface of a uniform cylinder, in terms of statistical behaviour.

Approaching the free end of a uniform cylinder, Park & Lee (Reference Park and Lee2000) observe a decrease in the turbulence intensity and the vortex formation length. A similar trend is observed in the modulation cell. The peak of the first diagonal component of the Reynolds stress tensor $\overline {uu}$ shows a reduced vortex formation length in the N-cell. Figure 8 indicates a decrease of approximately 16 %, as the peak of $\overline {uu}$ goes from $(x=2.70D,y=0.47D)$ at $z=-1.5D$ to $(x=2.28D,y=0.44D)$ at $z=-8.5D$. Similarly, the pressure is a further important downwash indicator, due to the interaction between the separated downwash flow and the vortices shed from the two sides of the cylinder (Park & Lee Reference Park and Lee2000). Zdravkovich et al. (Reference Zdravkovich, Brand, Mathew and Weston1989) and Williamson (Reference Williamson1989) confirm the appearance of a low-pressure cell. For the stepped cylinder, this characteristic cell is also present, as shown by the contraction of the three-dimensional pressure isosurface just behind the junction and in the upper part of the N-cell; see figure 13. In this area, the fluctuating field contributes to the mean in a manner distinct from any other parts of the flow; see the $\overline {pp}$ spatial distribution in figures 13(b,c). The fluctuating pressure field $\overline {pp}$ is located in the primary vortex cores everywhere except at $z \in (-0.5D,1.5D)$, where a third peak appears along the centreline at $x\approx 2.7D$. Related to this, the N-cell also shows large vertical velocity fluctuations $\overline {ww}$. The isosurface for $\overline {ww}= 71\,\%$ (of the maximum) is located in the modulation area. These fluctuations are expected as the downwash involves a vertical mass flow. More interestingly, they characterise the entire N-cell, not only the upper part where the POD mode 7 is prominent; see figure 14.

Figure 13. (a) Three-dimensional pressure isosurface $\bar {P}/\rho U^2=-0.18$ with a contraction behind the junction. Two-dimensional contour plots of the pressure fluctuations $\overline {pp}$ at (b) $z=-1.5D$ and (c) $z=-8.5D$.

Figure 14. (a) Time-averaged vertical velocity fluctuations in the $xz$ plane at $y=0$ with the isosurface $\overline {ww}/U^2=0.102$ (in red). Two-dimensional contour plots at different vertical locations: (b) $z=8.5D$,(c) $z=-1.5D$, and (d) $z=-8.5D$.

The modulation region turns out to have statistical characteristics similar to the free end of a uniform cylinder. On the whole, the N-cell is a region with a reduction of the base pressure, intense $\overline {pp}$ with large $\overline {ww}$ fluctuations, and a net reduction of the vortex formation length. Albeit the downwash mechanism is different with respect to the uniform cylinder, its statistical footprint is similar. First, by using the modal analysis and subsequently computing the statistics, it becomes possible to characterise this cell and illustrate the downward flow motion.

4.7. Reduced-order reconstruction

In the previous subsections, three distinct modes are clearly identified in the cylinder wake. Now we aim to describe their related dynamics. However, in a turbulent environment, this becomes unfeasible, and relying only on instantaneous flow visualisations is insufficient. Morton & Yarusevych (Reference Morton and Yarusevych2014b) suggest an aerodynamic model for the wake vortices interaction based on the idea that the circulation is inversely proportional to the frequency of vortex shedding and proportional to the square of the free stream velocity (Berger & Wille Reference Berger and Wille1972). Their experimental data seem to agree well with the model. Differently, in this case, the extracted POD modes are used to perform a reduced-order reconstruction that serves as a model that captures effectively the dynamics of the wake. The identification of such critical POD modes could facilitate the creation of a reduced-order model (ROM) by the Galerkin projection of the system dynamics on the first few, and most relevant, basis functions. The ROM can enable us to make accurate predictions for future scenarios of the subcritical regime.

The first six POD modes are essential to capture the L, S and N behaviour, i.e. the most energetic spatially coherent structures. Similarly, the modes from 7 to 10 are needed to encapsulate the interaction between the N-cell and the downwash:

(4.8)\begin{equation} \boldsymbol{u}(\boldsymbol{x}_i,t_j) = \sum_{k=0}^{\tilde{m}-1} \phi_k (\boldsymbol{x}_i)\, \boldsymbol{a}_k(t_j) , \end{equation}

where $\tilde {m}=11$, $\phi _k (\boldsymbol {x}_i)$ is a column of $\boldsymbol {X} = \boldsymbol{\mathsf{M}}^{-1/2} \hat {\boldsymbol{\mathsf{X}}}$, and $\boldsymbol {a}_k(t_j)$ is a row of $\boldsymbol {T} = \boldsymbol {\varSigma } \boldsymbol{\mathsf{W}}^{\rm T}$. In this way, the reconstructed field accounts for approximately $50\,\%$ of the total fluctuating energy. The comparison between the original and the reconstructed field at the same time confirms that the largest and more energetic structures are detected and visible. Filtering out the uncorrelated (small) scales enables us to discern the interactions occurring at the boundaries of the cell. Figure 15 shows the vortical structures for $\lambda _2=-2$ coloured by the vertical vorticity component. The red and blue colours refer to a $z$-vortex whose axis points upwards and downwards, respectively. The connection at the S–N and N–L boundaries is highlighted with dashed yellow lines. At the upper S–N interface, the large frequency difference between the S- and N-cells makes the vortex identification simpler.

Figure 15. Time sequence of the POD reduced-order reconstruction of flow around a three-dimensional stepped cylinder at $Re_D=1000$. The vortical structures ($\lambda _2 =-2\ U^2/D^2$) in the cylinder wake coloured with the vertical vorticity component at different times are shown. The time interval is $\Delta t = 0.25 D/U$, and the images are shown at the instants $t_0$, $t_0+5\Delta t$, $t_0 + 10\Delta t$, $t_0+120\Delta t$, $t_0+125\Delta t$ and $t_0+160\Delta t$. All times are normalised with the inflow velocity and the large cylinder diameter. In (d), the three different cells are clearly visible and highlighted by yellow dashed lines.

At the initial time $t_0$ (figure 15a), the N–L cells are in phase, whereas the S-cell starts to get out of phase. Just a little after, the dislocation occurs and the N-cell vortex splits into two sub-vortices that connect with two distinct S structures (at $t=t_0+5\Delta t$ with $\Delta t= 0.25$). The Y connection is clearly visible in figure 15(b). This unstable Y linking does not last long; see figure 15(c) at $t=t_0+10\Delta t$. When the half-loop connections between S-vortices happen, the N structure can form a half-loop with an adjacent N structure, or connect to another S vortex; see figures 15(c) and 15(d), respectively. The S, N and L vortical structures are displayed clearly by the dashed lines in figure 15(d). In addition, a larger penetration of S vortices into the modulation area is observed as the dislocation among the cells increases (Williamson Reference Williamson1989). The study of the N–L interface is more difficult since the interface is stable and the cells have almost identical frequencies ($St_N=0.188, 0.201$). However, the reduced-order reconstruction, based on a large number of snapshots, enables the characterisation of the N–L interaction. The N vortices form a direct connection with L for most of the cycles (figures 15a,c). The phase lag accumulation gets more and more significant, and eventually, the dislocation occurs; see figure 15(d) at $t=t_0+120\Delta t$. The N structures tend to split, and connections between consecutive vortices are visible in the L region (figure 15d). Completely out of phase in 15(e) at $t=t_0+125\Delta t$, the modulation region reaches its maximum extension. Eventually, the L- and N-cells are again in phase, and they behave like a single mode (figure 15f).

The POD technique applied to the flow field yields a reduced-order reconstruction that supports the aerodynamic model proposed by Morton & Yarusevych (Reference Morton and Yarusevych2014b). However, our results reveal some significant differences, particularly in the modulation region cycle. While Morton & Yarusevych (Reference Morton and Yarusevych2014b) reported a cyclic variation of the N-cell with its dimension changing over time, our observations suggest the occurrence of a direct connection between the L- and S-cells, when the phase lag is zero. This is supported also by the analytical model used for the superposition of two travelling modes; see § 4.5. The L and S alignment leads to the emergence of just two cells in the stepped cylinder wake for a brief period, representing a novel feature captured by our analysis through the inclusion of the downwash mode. Integrating this mode into the reconstruction enables us to grasp the behaviour that was previously overlooked, emphasising the significance of this phenomenon in the flow dynamics at the junction.