1. Introduction
The study of bluff-body wakes, particularly those formed around finite-length cylinders, has garnered significant attention due to the critical influence of coherent structures and vortex topology on unsteady flow dynamics. Such bluff bodies appear in a wide range of applications, including buildings, gas storage tanks and launch vehicles. With the recent advances in advanced air mobility, there is a growing need to understand these flow phenomena in the context of ensuring the safe and low-noise operation of future multi-rotor electric vertical take-off and landing systems (Zhou, Wang & Wang Reference Zhou, Wang and Wang2024). In this regard, it is essential to investigate the wake topology around finite cylinders under realistic conditions, and to explore advanced techniques for post-processing the turbulent wake dynamics.
A considerable body of research has investigated the wake dynamics of wall-mounted square cylinders (WMSCs), with particular attention to the role of aspect ratio (AR) and Reynolds number (
$Re$
) in shaping vortex topology. Early studies by Sakamoto & Arie (Reference Sakamoto and Arie1983) and Okamoto & Sunabashiri (Reference Okamoto and Sunabashiri1992) distinguished between Kármán-type and arch-type vortex shedding, emphasising that their occurrence is strongly dependent on AR, with Kármán vortices prevailing at high AR and arch-type at low AR. These works also showed that Strouhal number increases with AR and remains approximately uniform along the span for a given AR. Sumner, Heseltine & Dansereau (Reference Sumner, Heseltine and Dansereau2004) further reported the presence of coherent tip and base vortex structures, noting that base vortices become prominent only beyond a critical AR threshold. At lower AR, such as 3, the absence of both Kármán and base vortices highlighted a transitional regime distinct from higher-AR configurations.
Subsequent studies extended these insights by identifying complex three-dimensional vortex interactions within the wake. Wang & Zhou (Reference Wang and Zhou2009) proposed an arch-type model in which spanwise vortices connect to tip and base vortices to form a unified structure. Bourgeois, Sattari & Doolan (Reference Bourgeois, Sattari and Doolan2011) and Hosseini, Bourgeois & Martinuzzi (Reference Hosseini, Bourgeois and Martinuzzi2013) visualised alternating streamwise vortex loops and demonstrated that the wake topology is sensitive to boundary-layer thickness, which governs whether the wake assumes a dipole or quadrupole configuration. The influence of AR and boundary-layer interaction was further supported by direct numerical simulation and large-eddy simulation (LES) investigations (Saha Reference Saha2013; Rastan, Sohankar & Alam Reference Rastan, Sohankar and Alam2017; Sumner et al. Reference Sumner, Rostamy, Bergstrom and Bugg2017; Zhang et al. Reference Zhang, Cheng, An and Zhao2017), which reported variations in shedding coherence, loop formation and vortex symmetry as
$Re$
and inflow profiles were varied. More recently, Goswami & Hemmati (Reference Goswami and Hemmati2025a
,Reference Goswami and Hemmati
b
) performed LES of wall-mounted cylinders with varying AR, depth ratios and Re to further investigate the influence of these parameters on the stability of the shear layer separated from the leading edge of the cylinder. Their results revealed interactions between the Kelvin–Helmholtz-type vortices and the large-scale wake structures, through which momentum and energy are transferred from the Kelvin–Helmholtz vortices to the larger scales. Moreover, increasing the Reynolds number shifts the onset of the Kelvin–Helmholtz-type structures upstream, bringing them closer to the leading edge of the cylinder.
These findings collectively support the existence of a coherent vortex system spanning tip, base and spanwise vortices, typically characterised by a single dominant frequency. However, this assumption has been challenged by aeroacoustic studies. Moreau & Doolan (Reference Moreau and Doolan2013) first observed evidence of spanwise variation in shedding frequency, termed cellular shedding, particularly at high AR. Porteous, Moreau & Doolan (Reference Porteous, Moreau and Doolan2017) reinforced this observation, showing that different sections along a cylinder height could oscillate independently, forming discrete shedding cells. The number of these cells was shown to increase with AR, marking a significant departure from the notion of globally coherent shedding. Recent high-fidelity simulations by Wang et al. (Reference Wang, Ma, Hu and Thompson2024) have confirmed the presence of cellular shedding in WMSCs at moderate Reynolds numbers (
$Re \approx 1.5 \times 10^4$
). Their analysis demonstrated that both the aerodynamic force spectrum and the radiated acoustic field exhibit multiple dominant frequencies, each associated with distinct vertical regions of the cylinder. These findings underscore the importance of considering spanwise inhomogeneity in both physical modelling and acoustic analysis of finite-cylinder wakes.
According to previous research on wake topology and vortex dynamics behind finite-length cylinders, a variety of coherent structures, including tip, base and spanwise vortices, can emerge, often exhibiting complex interactions and spanwise variations in shedding frequency. These findings highlight the need for analysis techniques capable of capturing both spatial and temporal coherence in turbulent flows. Modal decomposition techniques, particularly proper orthogonal decomposition (POD) and its spectral counterpart (SPOD), have proven effective in extracting such coherent structures. Wang et al. (Reference Wang, Lam, Zu and Cheng2019) and Wang & Lam (Reference Wang and Lam2019) applied POD to WMSCs and demonstrated that the first two modes captured the dominant wake features across a range of ARs (1–6). Wang et al. (Reference Wang, Lam, Zu and Cheng2019) also found coexisting antisymmetric (von Kármán) and symmetric (arch-type) structures in POD modes, with their spatial dominance varying along the cylinder height. Cao et al. (Reference Cao, Liu, Zhou and Ren2022) used a combination of POD, SPOD and conditional POD to study a rectangular cylinder at
$Re = 7.4 \times 10^4$
, showing that SPOD better isolates frequency-dependent structures, avoiding temporal mixing seen in standard POD. These findings underscore the advantage of SPOD in revealing the temporal evolution and spectral composition of coherent structures in turbulent wakes. Acknowledging the spatial and temporal coherence of turbulent structures, Yauwenas et al. (Reference Yauwenas, Porteous, Moreau and Doolan2019) employed SPOD to analyse the wake dynamics behind a square cylinder. Their results demonstrated that increasing the AR shifts the dominant vortex shedding frequency to higher values. Specifically, for AR = 4.3, dominant Strouhal numbers of 0.085 and 0.107 were observed, whereas for AR = 10, the corresponding frequencies increased to 0.107 and 0.128. Maleki, Mankbadi & Golubev (Reference Maleki, Mankbadi and Golubev2025) performed a multi-plane SPOD analysis of the wake behind a WMSC (AR = 7), using data from 26 planes in three orthogonal directions. This approach enabled the identification of distinct coherent structures at multiple frequencies, including a secondary shedding cell near the base and a high-frequency tip vortex. A comparison between two- and three-dimensional SPOD revealed that three-dimensional SPOD may underrepresent localised, low-energy structures unless the domain is carefully restricted. In contrast, two-dimensional SPOD effectively captured both dominant and high-frequency low-energy modes using single-plane evaluations, highlighting the complementary strengths of each method in resolving complex wake dynamics.
The dynamics and coherence of wake structures behind bluff bodies are highly sensitive to inflow turbulence, significantly altering vortex formation, shedding behaviour and shear-layer interactions. In particular, Lee (Reference Lee1975) and Castro & Robins (Reference Castro and Robins1977) by investigating wall-mounted cylinders reported that increased turbulence intensity (TI) promotes reattachment of the separated shear layer, reduces the pressure energy associated with vortex shedding and diminishes the extent of the recirculation zone behind the cylinder. Bearman & Morel (Reference Bearman and Morel1983) also investigated the influence of free-stream turbulence on the flow around cylinders. They identified three primary mechanisms through which the incoming turbulence interacts with the bluff-body flow: (i) accelerated transition to turbulence, (ii) enhanced mixing and entrainment and (iii) distortion of the free-stream turbulence by the mean flow. Building upon these findings, Nakamura & Ohya (Reference Nakamura and Ohya1984) further demonstrated that small-scale turbulence enhances entrainment in the near wake, thereby accelerating the growth of the separated shear layer, whereas large-scale turbulence reduces spanwise coherence and suppresses organised vortex shedding. More recently, Hearst, Gomit & Ganapathisubramani (Reference Hearst, Gomit and Ganapathisubramani2016) conducted experiments on the wake of a wall-mounted cube under turbulent boundary layers with four levels of TI, ranging from
$0.5\,\%$
to
$8.9\,\%$
. While the mean streamlines remained largely unaffected, likely due to the relatively low turbulence levels, flow reattachment on the cube’s surfaces was observed, along with a reduction in the energy associated with the dominant shedding frequency. Song et al. (Reference Song, Ping, Zhu, Zhou, Bao, Cao and Han2022) investigated the influence of incoming turbulence on the flow over a wall-mounted circular cylinder and found that, although the distribution of the mean pressure coefficient (
$\overline {C_p}$
) under turbulent inflow conditions is qualitatively similar to that observed in the uniform inflow case, the magnitude of the pressure fluctuations (
$C_{p,{\textit{rms}}}$
) increases substantially. Specifically, the mean drag coefficient (
$\overline {C_d}$
) rises by approximately 29.4 %, while the root mean square of the lift coefficient exhibits an increase of nearly 500 %. The influence of inflow turbulence and boundary-layer thickness has also been examined numerically by Chen et al. (Reference Chen, Li, Sun and Liang2022, Reference Chen, Li, Xue, He, Wang and Liang2023), who employed the synthetic eddy method (Jarrin et al. Reference Jarrin, Benhamadouche, Laurence and Prosser2006) and demonstrated that both parameters significantly affect the frequency content and spatial coherence of the SPOD modes. Specifically, thicker boundary layers led to the merging of initially distinct frequency peaks, while the presence of turbulence reduced the energy of high-frequency tip vortices, indicating a dampening effect on coherent shedding. Zhang, Ooka & Kikumoto (Reference Zhang, Ooka and Kikumoto2021) conducted a three-dimensional SPOD analysis of a square cylinder (AR = 2), using inlet velocity fields obtained from a separate LES simulating atmospheric boundary-layer flow. However, no direct comparison was made between laminar and turbulent inflow conditions. Also as reported by Maleki et al. (Reference Maleki, Mankbadi and Golubev2025), three-dimensional SPOD may underrepresent low-energy, high-frequency structures unless a localised analysis domain is carefully selected.
Despite prior investigations into the interaction between incoming turbulence and bluff-body wakes, most studies have considered only a narrow range of turbulence intensities. Direct comparisons between cases with and without inflow turbulence fluctuations are rarely reported, and the underlying mechanisms by which turbulence modifies coherent structures, particularly the changes in their spatial organisation and spectral content, remain largely unexplored. These gaps motivate the present work, which aims to investigate the influence of incoming turbulence on the wake dynamics of a WMSC with finite AR. To this end, a momentum-based source term originally developed by Golubev et al. (Reference Golubev, Dreyer, Visbal and Golubev2009, Reference Golubev, Nguyen, Brodnick and Visbal2011) is employed and coupled with a synthetic turbulence generation method. This approach enables the injection of divergence-free turbulent fluctuations directly within the computational domain, avoiding inflow boundary complications and ensuring robust numerical performance across stretched grids. Recent developments by Jiang (Reference Jiang2023) and Shao, Jiang & Chen (Reference Shao, Jiang and Chen2024) advanced the formulation to enable exact gust injection in idealised aerofoil settings, validated against analytical models such as the Sears function. In contrast, the present study leverages the momentum–turbulence coupling to generate realistic inflow disturbances, which are assessed against the von Kármán energy spectrum and validated through comparison with experimental measurements in the wake of a wall-mounted bluff body. To characterise the impact of inflow turbulence on the resulting unsteady flow field, a combination of post-processing techniques is employed, including time-averaged flow analysis, spectral analysis of velocity signals, vorticity-based diagnostics and SPOD. These analyses reveal how varying levels of inflow TI (10 %–30 %) differentially alter the coherent structures around the cylinder, including their influence on vortex-shedding dynamics, shear-layer behaviour and recirculation topology, providing insight into the fundamental mechanisms governing turbulence–structure interactions in bluff-body wakes.
2. Numerical simulation
2.1. Large-eddy simulation
Given the low velocities considered in this study, where the maximum Mach number is below 0.2, the air can be reasonably treated as incompressible. The governing equations for continuity and momentum, applicable to an incompressible fluid with constant properties, are expressed as follows:
where
$u_i$
is the filtered velocity,
$p$
is the filtered pressure and
$\rho$
is the fluid density. In this work, the wall-adapting local eddy-viscosity (WALE) model was employed as the subgrid-scale closure, leveraging the Boussinesq approximation. The WALE model, initially introduced by Nicoud & Ducros (Reference Nicoud and Ducros1999), has proven effective in capturing near-wall flow behaviours and has been validated in multiple studies (Okaze et al. Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021; Wornom et al. Reference Wornom, Ouvrard, Salvetti, Koobus and Dervieux2011; Foroutan et al. Reference Foroutan, Tang, Heist, Perry, Brouwer and Monbureau2018; Zhao et al. Reference Zhao, Wang, Gong and Xu2023). A significant advantage of this model is that its eddy viscosity approaches zero in the vicinity of walls, eliminating the need for additional damping functions. To compute subgrid-scale turbulent stresses we have
where
$\widetilde {S}_{\textit{ij}}$
is the resolved strain-rate tensor, defined as
The eddy viscosity
$\nu _{\textit{SGS}}$
in the WALE model is calculated as
\begin{equation} \nu _{\textit{SGS}} = (C_W \varDelta )^2 \frac {\left ( S_{\textit{ij}}^d S_{\textit{ij}}^d \right )^{3/2}}{\left ( \widetilde {S}_{\textit{ij}} \widetilde {S}_{\textit{ij}} \right )^{5/2} + \left ( S_{\textit{ij}}^d S_{\textit{ij}}^d \right )^{5/4}}. \end{equation}
In (2.5),
$C_w$
is a constant specific to the WALE model with a value of 0.325 and
$\varDelta$
is the grid size that defines the subgrid length scale and calculated based on
$\varDelta = V^{1/3}$
, where
$V$
is the volume of a cell in the computational domain. Additionally, the term
$S_{\textit{ij}}^d$
in (2.5) represents the traceless symmetric part of the squared velocity gradient tensor and is defined as
\begin{equation} S_{\textit{ij}}^d = \frac {1}{2} \left ( \left (\frac {\partial u_i}{\partial x_{\!j}}\right )^2 + \left (\frac {\partial u_{\!j}}{\partial x_i} \right )^2\right ) - \frac {1}{3} \delta _{\textit{ij}}\left ( \frac {\partial u_k}{\partial x_k} \right )^2. \end{equation}
2.2. Momentum source term formulation
As mentioned earlier, this study extends the momentum source framework introduced by Golubev et al. (Reference Golubev, Dreyer, Visbal and Golubev2009, Reference Golubev, Nguyen, Brodnick and Visbal2011), who embedded a three-dimensional turbulence field into the Navier–Stokes equations via a localised source term. The current study implements the approach to generate turbulence spectra consistent with the von Kármán energy distribution.
To this end, a three-dimensional, divergence-free and convected perturbation velocity field
$u_i(x_{\!j}, y)_{i,j=1,2,3}$
is generated by prescribing a momentum source term within an arbitrary region upstream of the body in the computational domain. This synthetic turbulence field is described in terms of a Fourier spectrum constructed from a superposition of harmonic functions:
\begin{equation} u_i(x_{\!j},t) = \sum ^{N} _{n=1} \big(a^n _i \cos \big[k^n_{\!j} x_{\!j} - \omega _n t\big] + b^n _i \sin \big[k^n _{\!j} x_{\!j} - \omega _n t\big]\big), \end{equation}
where, without loss of generality, the synthesised unsteady disturbance field is assumed to be convected by a uniform mean flow
$u_\infty$
, aligned with the
$x_1$
direction, ensuring that the modal perturbation frequency
$\omega _n$
is related to the wavenumber
$k^n_1$
by
$\omega _n = u_\infty k^n_1$
. Additionally, the disturbance is considered incompressible (divergence-free), leading to the following coupling between the modal amplitudes and wavenumbers:
Interactions between the perturbation modes and other waves are considered negligible. Thus, the momentum source terms
$S_i$
generating the disturbance velocity field in (2.7) should satisfy the following relationships:
In Appendix A, it is shown that the solution for the momentum source terms imposed in the region
$| x_1 - x_{1s} | \leq \pi / b$
of the computational domain is then obtained in the form
\begin{equation} S_i(x_{\!j}, t) = \begin{cases} \sum ^\infty _{n=1} \omega _n K_n g(x_1) \left ( b_i^n \sin \sigma _n - a_i^n \cos \sigma _n \right ) & \text{for } i = 1, \\[5pt] \sum ^\infty _{n=1} U_\infty K_n g'(x_1) \left ( b_i^n \cos \sigma _n + a_i^n \sin \sigma _n \right ) & \text{for } i = 2,3, \end{cases} \end{equation}
where
$K_n$
is a constant given by
\begin{equation} K_n = \frac {\left (k_1^n\right )^2 - b^2}{b^2 \sin \left (\dfrac {\omega _n \pi }{b U_\infty }\right )} \end{equation}
and
$g(x_1)$
is the limiter function defined in this work as
\begin{equation} g(x_1)=\left \{\begin{array}{cc} \dfrac {1}{2}\left (1+\cos \left [b(x_1-x_{1s})\right ]\right ) & \mbox{for}\;\;\; |x_1-x_{1s}|\leq \dfrac {\pi }{b}, \\[8pt] 0 & \mbox{for}\;\;\; |x_1-x_{1s}|\geq \dfrac {\pi }{b}. \end{array}\right \} \end{equation}
In the above equations,
$g^\prime (x_1) = {\rm d}g/{\rm d}x_1$
and
$\sigma = \omega _nt-k^n_1 x_{1s} - k^n_2 x_2 - k^n_3 x_3$
. The momentum source terms specified in (2.10) are centred in an upstream region at
$x_1 = x_{1s}$
. In the numerical procedure, all modes with different wavelengths are superimposed at the downstream boundary of the source region. The nine coefficients
$a^n_i$
,
$b^n_i$
and
$k^n_i$
for
$i = 1, 2, 3$
in (2.10) are derived from the method outlined by Shur et al. (Reference Shur, Spalart, Strelets and Travin2014).
2.3. Synthetic turbulence model
According to the work of Shur et al. (Reference Shur, Spalart, Strelets and Travin2014), anisotropic turbulence can be generated in the source region using the synthetic turbulence generation (STG) method. Specifically, the turbulent flow field is constructed through a three-dimensional Fourier spectrum that aligns with the specified energy distribution for the perturbation velocity. Similarly to other STG methods, the velocity vector at a point
$\boldsymbol{r} = \{x, y, z\}$
is given by
where
$\overline {U}_{i}$
is the mean velocity vector from the Reynolds-averaged Navier–Stokes solution and
${u_i}'({r}, t)$
represents the vector of velocity fluctuations. The fluctuations
${u_i}'({r}, t)$
are defined such that their second moment
$\langle u_i' u_{\!j}' \rangle$
at the LES inflow matches the Reynolds stress tensor
$\boldsymbol{R}$
, obtained from the Reynolds-averaged Navier–Stokes solution. This is done using a Cholesky decomposition:
where
\begin{align} \hat {A}={a_{\textit{ij}}} = \begin{pmatrix} \sqrt {R_{11}} & 0 & 0 \\ R_{21}/a_{11} & \sqrt {R_{22} - a_{21}^2} & 0 \\ R_{31}/a_{11} & (R_{32} - a_{21} a_{31})/a_{22} & \sqrt {R_{33} - a_{31}^2 - a_{32}^2} \end{pmatrix} \end{align}
and
$R_{\textit{ij}} = \langle u_i' u_{\!j}' \rangle$
are components of the Reynolds stress tensor. The velocity fluctuations
$\boldsymbol{u}'(\boldsymbol{r}, t)$
are expressed using the auxiliary fluctuation vector
$\boldsymbol{v}'(\boldsymbol{r}, t)$
:
where
$\boldsymbol{v}'(\boldsymbol{r}, t)$
satisfies the conditions
$\langle v_i' v_{\!j}' \rangle = \delta _{\textit{ij}}$
to ensure the correct Reynolds stress. The auxiliary vector
$\boldsymbol{v}'(\boldsymbol{r}, t)$
is computed by superimposing spatiotemporal Fourier modes:
\begin{equation} \boldsymbol{v}'(\boldsymbol{r}, t) = 2 \sqrt {\frac {3}{2}} \sum _{n=1}^{N} \sqrt {q^n} \cos (k^n \boldsymbol{d}^n \boldsymbol{\cdot }\boldsymbol{r}' + \phi ^n). \end{equation}
In this context,
$N$
represents the total number of modes, and
$q_n$
is the normalised amplitude of the mode
$n$
, based on the local energy spectrum. The wavenumber vector’s amplitude for mode
$n$
is given by
$k_n$
, where
$K_n = k_n \boldsymbol{d}_n$
, with
$\boldsymbol{d}_n$
being a random unit direction vector uniformly distributed over a sphere. The unit vector
$\boldsymbol{\sigma }_n$
is orthogonal to
$\boldsymbol{d}_n$
(
$\boldsymbol{\sigma }_n \boldsymbol{\cdot }\boldsymbol{d}_n = 0$
) (Kraichnan Reference Kraichnan1970) and its direction in the plane normal to
$\boldsymbol{d}_n$
is determined by a random angle uniformly distributed in the range
$[0, 2\pi )$
. Additionally,
$\phi _n$
is the random phase of mode
$n$
, also uniformly distributed in the range
$[0, 2\pi )$
. These random distributions are defined for each mode
$n$
. The normalised amplitudes of the modes
$q^n$
in (2.17) are given by
\begin{equation} q^n = \frac {E(k^n) \Delta k^n}{\sum\limits _{n=1}^{N} E(k^n) \Delta k^n}, \quad \sum _{n=1}^{N} q^n = 1 ,\end{equation}
where
$E(k^n)$
represents the spatial spectrum of turbulent kinetic energy, described by a modified von Kármán spectrum, and defined by the following expression:
\begin{equation} E(k) = \frac {(k/k_e)^4}{\left [ 1 + 2.4 \left ( k/k_e \right )^2 \right ]^{17/6}} f_\eta f_{\textit{cut}} ,\end{equation}
where
$k_e$
is the wavenumber at the spectral peak and
$f_\eta$
and
$f_{{\textit{cut}}}$
are empirical correction functions. Additional details about the derivation and equations can be found in Shur et al. (Reference Shur, Spalart, Strelets and Travin2014).
2.4. Computational domain and boundary conditions
In this study, a wall-mounted rectangular cylinder with an AR of 2 is employed. This configuration has been previously investigated in both experimental and numerical studies (Okaze et al. Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021; Zhang et al. Reference Zhang, Ooka and Kikumoto2021). The current simulations are performed using the open-source computational fluid dynamics (CFD) software OpenFOAM-v2412 with the pisoFoam solver for incompressible flows. Four cases are considered, one undisturbed case and three disturbed cases with turbulence intensities of
$10\,\%$
,
$20\,\%$
and
$30\,\%$
(here, TI is defined as
$\text{TI} = {u^{\prime}_{\textit{rms}}/(y = h)}{U_h} \times 100\,\%$
, where
$U_h$
is the mean inflow velocity at the cylinder height). The undisturbed case is solved using the incompressible Navier–Stokes equations within a LES framework, employing the WALE model for subgrid-scale turbulence closure. For the disturbed cases, a momentum-based turbulence source term, defined in (2.10), is added to the right-hand side of the Navier–Stokes equations. This modification is implemented by customising the pisoFoam solver in OpenFOAM-v2412. A non-uniform inlet velocity profile is prescribed with a reference velocity of
$U_h = 2.89\,\mathrm{m}\,\mathrm{s}^{-1}$
, where
$U_h$
represents the undisturbed inflow velocity at the cylinder height. The velocity profile is defined as
$U = U_h (y/h)^{0.27}$
, which matches the experimental profile of Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021). This profile was also adopted by Zhang et al. (Reference Zhang, Ooka and Kikumoto2021), who simulated a similar case. The Reynolds number, based on the reference velocity and cylinder height (
$h = 0.2$
m), is
$Re_h = 43{\,}000$
.
As illustrated in figure 1, the cylinder is positioned at a distance of
$4h$
from the inlet and
$10h$
from the outlet, where
$h$
is the cylinder height and
$d$
is the cylinder width. Following the recommendations of Sohankar, Norberg & Davidson (Reference Sohankar, Norberg and Davidson1998) and Saha (Reference Saha2013), the blockage ratio
$\beta = (d \times h)/(H \times W)$
, where
$H$
and
$W$
are the computational domain dimensions in the
$y$
and
$z$
directions, respectively, must remain below 0.05, and the domain height should satisfy
$H \gt h + 5d$
. In compliance with these criteria, the blockage ratio is set to 0.016, and the normalised domain height,
$(H - h)/d$
, is set to 8. A structured mesh consisting of 17 blocks surrounding the cylinder model was generated using Pointwise v2202-1, as illustrated in figure 2. The width of the cylinder was discretised into 70 uniformly spaced cells, each with a size of
$d/70$
. It is worth noting that Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021) and Zhang et al. (Reference Zhang, Ooka and Kikumoto2021) employed 40 and 32 cells, respectively, across the cylinder. The cylinder height was meshed with the same resolution as its width, and a constant grid spacing was maintained from the ground up to the top of the cylinder. Away from the cylinder, the grid spacing was gradually increased towards the domain boundaries using an expansion ratio of 1.07 in all directions. Based on these parameters, the average
$y^+$
value on the cylinder surface ranges from 1.4 to 14, which roughly classifies the present simulations within the high-resolution LES regime, according to the criterion proposed by Wagner, Hüttl & Sagaut (Reference Wagner, Hüttl and Sagaut2007). In total, approximately
$ 5.6 \times 10^6$
grid cells were used in each case of the present simulations. It should also be noted that separate simulations were performed to analyse the sensitivity of the results to the domain size and to assess grid independence for the case with
$\text{TI} = 10\,\%$
. The corresponding results are presented in Appendix B.
Illustration of the computational domain, showing the location of the cylinder relative to the domain boundaries and the allocated source region.

Symmetry boundary conditions were applied at the right, left and top boundaries of the computational domain. At the outflow boundary, pressureInletOutletVelocity and totalPressure conditions were imposed for the velocity and pressure fields, respectively. The numerical scheme employs the backward time integration method, which is a second-order-implicit scheme. For the spatial gradients, the second-order leastSquares scheme is used. To achieve an optimal balance between low numerical dissipation and numerical stability, the Gauss filteredLinear2V 1 0 scheme is employed for computing the divergence of the velocity field. This scheme is based on the second-order central difference method (Gauss linear), but in regions with high gradients, where central differencing may lead to numerical oscillations, it locally blends with the first-order upwind scheme to suppress artificial fluctuations. This approach has previously been shown to yield improved results compared with the LUST divergence scheme (Okaze et al. Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021; Maragkos et al. Reference Maragkos, Verma, Trouvé and Merci2022). The Gauss linear scheme is used for the remaining terms. Convergence criteria for each time step are specified as scaled residuals of
$10^{-5}$
for pressure and
$10^{-6}$
for velocity. The time step size,
$\Delta t = 1 \times 10^{-4}$
, is chosen to ensure numerical stability and maintain the Courant–Friedrichs–Levy number below
$0.7$
.
Computational grid around the WMSC.

2.5. Spectral proper orthogonal decomposition parameters
In this study, SPOD was employed to analyse the turbulent flow field around a wall-mounted cylinder with an AR of 2. The detailed formulation of the SPOD methodology can be found in Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018) and Schmidt & Colonius (Reference Schmidt and Colonius2020). In the present research, the SPOD-Python code developed by He et al. (Reference He, Fang, Rigas and Vahdati2021) was employed. To perform SPOD on the flow around the cylinder, 17 planes were defined in each of the
$x$
-,
$y$
- and
$z$
-normal directions, as summarised in table 1 and shown in figure 3. The time history of all three velocity components was recorded over the final 30 s of the simulation, with the total simulation time being 40 s. The first 10 s were used to allow the flow to reach a statistically stationary state and correspond to approximately 25 dominant vortex-shedding cycles. The subsequent 30 s of data sampling cover around 75 dominant vortex-shedding cycles. It should also be noted that this 30 s sampling interval corresponds to more than 2000 cycles of the highest-frequency phenomenon, which is associated with Kelvin–Helmholtz-type structures forming along the top and side surfaces of the cylinder.
Definition of planes used in the SPOD analysis along the
$x$
-,
$y$
- and
$z$
-normal directions.

Data were sampled at intervals of 0.001 s, corresponding to a sampling frequency of 1000 Hz. For the SPOD analysis, the resulting 30 000 snapshots at each plane were divided into 28 overlapping blocks, each containing 2048 snapshots with a 50 % overlap. Similar to the approach of previous studies (Maleki et al. Reference Maleki, Mankbadi and Golubev2024, Reference Maleki, Singh and Golubev2025), the data in each snapshot were interpolated using a cubic interpolation scheme. However, due to the finer grid resolution in the present study, a uniform interpolation spacing of
$d/40$
was adopted in all three spatial directions:
$x$
,
$y$
and
$z$
. The frequency resolution of the SPOD spectra, given by the relation
$\Delta f = {f_s}/{N_{\textit{FFT}}}$
, where
$f_s$
and
$N_{\textit{FFT}}$
are the sampling frequency and the number of snapshots in each block, respectively, is
$0.488$
Hz (Schmidt & Colonius Reference Schmidt and Colonius2020). The Welch method is employed to reduce energy leakage in the statistical data analysis. Furthermore, the turbulent kinetic energy norm, as described by Schmidt & Colonius (Reference Schmidt and Colonius2020), was used as the weight matrix
$ W$
.
Schematic showing the SPOD analysis planes around the cylinder.

3. Results
To analyse the effect of incoming turbulence on the flow around the cylinder, four simulation cases were defined. The first case represents an undisturbed inflow with a non-uniform velocity profile given by
$U = U_h \times (y/h)^{0.27}$
. This configuration serves as the baseline, against which both the time-averaged and instantaneous CFD results, as well as the SPOD outcomes, are compared. The remaining three cases involve turbulent inflow conditions introduced through the momentum source term described previously. These cases correspond to varying turbulence intensities ranging from
$10\,\%$
to
$30\,\%$
. It is demonstrated that the case with
$10\,\%$
TI is generally sufficient to generate the necessary velocity fluctuations at the inflow. The results from this case meet the requirements related to inflow turbulence characteristics and velocity levels within the shear layer of the cylinder. The additional cases with higher turbulence intensities (
$20\,\%$
and
$30\,\%$
) are included to investigate the impact of stronger inflow disturbances on the overall flow field.
3.1. Momentum turbulence source term
The STG approach developed by Shur et al. (Reference Shur, Spalart, Strelets and Travin2014) is capable of producing anisotropic turbulent inflow conditions. In the present work, however, isotropic turbulence is assumed to simplify the implementation and reduce the computational cost associated with the volumetric source term. Despite this simplification, comparisons with experimental data demonstrate that the isotropic assumption provides satisfactory agreement in the mean velocity field. A similar volumetric momentum source formulation has been previously employed by Maleki, Mankbadi & Golubev (Reference Maleki, Mankbadi and Golubev2025a ,Reference Maleki, Grace, Mankbadi and Golubev b ) and Voropayev et al. (Reference Voropayev, Marques, Maleki, Grace, Lyrintzis, Mankbadi and Golubev2025) to generate synthetic atmospheric turbulence and time-harmonic gusts for studying the interaction of turbulent inflow and single-frequency disturbances with cylinder and rotor configurations, respectively. It is also worth emphasising that, in contrast to conventional STG approaches, where velocity fluctuations are imposed either directly at the inlet boundary (Huang, Li & Wu Reference Huang, Li and Wu2010; Aboshosha et al. Reference Aboshosha, Elshaer, Bitsuamlak and El Damatty2015) or in its immediate vicinity (Kim, Castro & Xie Reference Kim, Castro and Xie2013), the present formulation introduces the fluctuating velocity components via a volumetric momentum source term distributed within the computational domain. This implementation allows flexible positioning of the source region away from the inlet and lateral boundaries, thereby improving the robustness of the boundary conditions. Furthermore, since the convective transport of the generated turbulence is decoupled from the local grid resolution between the source region and the solid surfaces, the method accommodates the use of stretched grids and consequently enhances computational efficiency.
As illustrated in figure 1, the turbulence source term is prescribed at a distance of
$x/h = 1.25$
upstream of the cylinder. The von Kármán energy spectrum, defined in (2.19), is discretised into 999 Fourier harmonics. The nine coefficients
$a^n_i$
,
$b^n_i$
and
$k^n_i$
in (2.10) are obtained from the trigonometric expansion of (2.17). In a fully anisotropic turbulence field, these coefficients must be computed across a two-dimensional plane using the method proposed by Shur et al. (Reference Shur, Spalart, Strelets and Travin2014). However, under the assumption of isotropy, calculating the coefficients in a single point is sufficient to generate the perturbation velocity field. To reduce computational expense, this study employs the isotropic formulation by evaluating the coefficients at a single spatial location rather than over an entire plane. It is also important to note that, consistent with the isotropic nature of the imposed turbulence, the off-diagonal components of the Reynolds stress tensor are assumed to be zero, while the diagonal components are considered equal. These diagonal terms vary across the three turbulence cases in accordance with the prescribed TI.
The target point chosen for the calculation of the coefficients was selected at
$(x/h, y/h, z/h)=(2.75, 1, 0)$
. This point lies at the same vertical height as the cylinder top surface and is located
$1.25h$
upstream of the front wall of the model. To evaluate the characteristics of the generated turbulence field, a probe point was defined between the volume source region and the cylinder. As shown in figure 1, this point is located at
$(x/h, y/h, z/h) = (3.5, 1, 0)$
. Figure 4 presents the power spectral density (PSD) of the streamwise, wall-normal and spanwise velocity components at the selected probe point for three inflow cases, with turbulence intensities ranging from
$10\,\%$
to
$30\,\%$
. The x-axis shows the non-dimensional frequency, calculated based on the reference velocity (
$U_h = 2.89$
m s−1) and the cylinder width (
$d = 0.1$
m). The solid black line indicates a
$-5/3$
slope, characteristic of the inertial subrange in the von Kármán turbulence spectrum. Its position is fixed to highlight the effect of increasing TI across the three cases. The dotted black line corresponds to the prescribed turbulence spectrum at the target point
$(x/h, y/h, z/h) = (2.75, 1, 0)$
, which was used to calculate the required coefficients. The spectrum was computed using the formulation by Shur et al. (Reference Shur, Spalart, Strelets and Travin2014) and implemented in a Python code, and the red line is associated with the CFD results at the probe point.
Comparison of the PSD of the streamwise, wall-normal and spanwise velocity components at the target point
$(x/h, y/h, z/h) = (2.75, 1, 0)$
and the probe point located at
$(x/h, y/h, z/h) = (3.5, 1, 0)$
.

According to figure 4, the red line follows the general trend and behaviour of the target spectrum. However, in the cases with 20 % and 30 % TI, the resulting spectra extend beyond the range originally prescribed by the target spectrum. This deviation is primarily attributed to the role of the
$f_{{\textit{cut}}}$
function embedded in the von Kármán spectrum formulation (2.19). This function is designed to truncate the spectral energy at higher frequencies in the target spectrum generated via the Python code, target spectra. In contrast, the turbulence imposed in the simulation is governed by the Navier–Stokes equations, which do not include such an artificial spectral cutoff. As a result, the generated turbulence retains energy at higher frequencies, depending on the grid resolution and its ability to resolve and sustain small-scale motions. Consequently, the actual spectrum reflects both the imposed energy distribution and the numerical characteristics of the grid.
To quantify the effective turbulence levels achieved in the simulations, the LES spectra were integrated to obtain the corresponding turbulence intensities. The resulting values are 10.2 %, 19.2 % and 25.6 % for the cases nominally prescribed as 10 %, 20 % and 30 %, respectively. The first two cases show close agreement with their target intensities, while the highest-intensity case exhibits a modest underprediction, which, despite confirmation of grid adequacy through a grid-independence study, is attributed to the limited resolution of small-scale motions at such high fluctuation levels. Nevertheless, the overall spectral shape and energy content remain consistent with the prescribed von Kármán distribution, confirming that the implemented momentum source framework reproduces physically realistic inflow turbulence fields.
To better understand the behaviour of the turbulence source term, instantaneous vorticity vectors in the
$z$
,
$y$
and
$x$
directions, non-dimensional with respect to the reference velocity
$U_h = 2.89$
m s−1 and cylinder height
$h = 0.2$
m, are presented in figures 5–8. For consistency and comparative analysis, the non-dimensional vorticity contours are plotted within a fixed range of
$-10$
to
$10$
across all cases and directional components. Figure 5 illustrates the vorticity distribution around the cylinder in the undisturbed case. As seen in the figure, flow separation occurs at the leading edges of both the top and side surfaces of the structure. The unsteady separation, oscillation and rolling-up of the shear layer over these surfaces give rise to Kelvin–Helmholtz instabilities (Wei & Smith Reference Wei and Smith1986; Dong et al. Reference Dong, Karniadakis, Ekmekci and Rockwell2006).
Three components of non-dimensional vorticity in the undisturbed case, visualised in three orthogonal planes: (a) spanwise vorticity (
$\omega _z$
) in the
$z$
-normal plane at
$z/h = 0$
; (b) vertical vorticity (
$\omega _y$
) in the
$y$
-normal plane at
$y/h = 0.75$
; and (c) streamwise vorticity (
$\omega _x$
) in the
$x$
-normal plane at
$x/h = 4.25$
.

Figures 6–8 show the non-dimensional vorticity components corresponding to turbulence intensities of
$10\,\%$
,
$20\,\%$
and
$30\,\%$
, respectively. As evident from these figures, along with figure 1 which shows that the momentum source term occupies the entire domain downstream of
$x/h = 2.75$
, the generated vorticity structures extend fully to the top and side boundaries of the computational domain. However, the influence of this region, the volume momentum source term, can be modulated using the function in (A9) which is defined in Appendix A and originally introduced by Golubev et al. (Reference Golubev, Dreyer, Visbal and Golubev2009). According to figure 6 all three components of vorticity are successfully generated in the entire computational domain. Furthermore, as the TI increases to
$20\,\%$
and
$30\,\%$
, shown in figures 7 and 8, respectively, the magnitude of vorticity correspondingly increases. This increase is visually evident in the figures through the more intense coloration, indicating stronger rotational motion within the flow field. A noteworthy observation from these figures is the stretching and dissipation of vorticity near the domain boundaries and at downstream locations relative to the top surface of the cylinder. These effects primarily originate from the grid stretching and are not directly associated with the behaviour of the volume source term. However, this highlights the importance of maintaining a low grid expansion ratio throughout the computational domain to ensure more accurate representation of the turbulence field and to preserve vorticity structures.
Three components of non-dimensional vorticity in the disturbed case with a TI of
$10\,\%$
, visualised in three orthogonal planes: (a) spanwise vorticity (
$\omega _z$
) in the
$z$
-normal plane at
$z/h = 0$
; (b) vertical vorticity (
$\omega _y$
) in the
$y$
-normal plane at
$y/h = 0.75$
; and (c) streamwise vorticity (
$\omega _x$
) in the
$x$
-normal plane at
$x/h = 4.25$
.

Three components of non-dimensional vorticity in the turbulent inlet case with
$\text{TI}=20\,\%$
, visualised in orthogonal planes: (a) spanwise vorticity (
$\omega _z$
) in the
$z$
-normal plane at
$z/h = 0$
; (b) vertical vorticity (
$\omega _y$
) in the
$y$
-normal plane at
$y/h = 0.75$
; and (c) streamwise vorticity (
$\omega _x$
) in the
$x$
-normal plane at
$x/h = 4.25$
.

Three components of non-dimensional vorticity in the turbulent inlet case with
$\text{TI} = 30\,\%$
, visualised in orthogonal planes: (a) spanwise vorticity (
$\omega _z$
) in the
$z$
-normal plane at
$z/h = 0$
; (b) vertical vorticity (
$\omega _y$
) in the
$y$
-normal plane at
$y/h = 0.75$
; and (c) streamwise vorticity (
$\omega _x$
) in the
$x$
-normal plane at
$x/h = 4.25$
.

In this regard, to investigate the effect of grid stretching, an additional simulation was performed using isotropic grid cells with a uniform resolution of
$d/12.5$
in the
$x$
,
$y$
and
$z$
directions. This case was simulated in the absence of the cylinder to isolate and examine the behaviour of the convected vortices, and included a momentum source term corresponding to the disturbed inflow case with a TI of 10 %. The same numerical schemes and domain dimensions were retained. The resulting non-dimensional
$y$
component of vorticity at plane
$y_4$
is shown in figure 9. The contours indicate that vortices are generated and convected without significant attenuation in amplitude. However, it should be noted that this simulation corresponds to a TI of 10 %, and higher turbulence intensities may require finer grid resolution to adequately capture and preserve the magnitude of the vortical structures. Furthermore, in contrast to the results shown in figures 5–8, where the vortical structures appeared distorted due to grid-cell stretching, no such behaviour was observed in figure 9. It should also be noted that the stretching of vortices near the boundaries does not affect the cylinder flow field. This is further demonstrated in the following section through an additional simulation designed to verify this effect.
Non-dimensional
$y$
component of vorticity in the
$y$
-normal plane
$y/h=0.8$
.

3.2. Time-averaged results
To enable comparison with the experimental data of Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021), and to analyse the effect of incoming turbulence on the flow around the cylinder, the final 70 s of a total 80 s simulation were time-averaged for all four cases. This 70 s interval corresponds to approximately 250 dominant vortex-shedding cycles in the current simulation set-up. Figure 10 presents a comparison of the three non-dimensional velocity components and their root-mean-square (RMS) values along seven vertical lines located at
$x/h = 3.875$
,
$4.25$
,
$4.5$
,
$4.625$
,
$4.875$
,
$5.25$
and
$5.875$
on the
$z$
-normal plane at
$z/h = 0$
. The data from the first line at
$x/h = 3.875$
capture the time-averaged and RMS velocity components induced by atmospheric conditions, as reported in the experimental study by Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021). The comparison between the experimental results of Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021) and the simulated cases in this study at the line
$x/h = 3.875$
indicates that the disturbed case with
$10\,\%$
TI, implemented via the volumetric momentum source term, effectively reproduces the required level of velocity fluctuations. This observation, together with the PSD results shown in figure 4, confirms the effectiveness of the momentum source term in generating synthetic turbulence.
Time-averaged comparison of non-dimensional velocity components and their RMS values along seven vertical lines located at
$x/h = 3.875$
,
$4.25$
,
$4.5$
,
$4.625$
,
$4.875$
,
$5.25$
and
$5.875$
on the
$z$
-normal plane at
$z/h = 0$
, compared with the experimental data of Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021).

Analysis of velocity data along two vertical lines at
$x/h = 4.25$
and
$4.5$
indicates that the presence of inflow turbulence exerts a more pronounced influence on the flow field than variations in the TI level. In other words, the differences among the disturbed cases are relatively minor compared to the difference between the undisturbed case and the case with
$\text{TI} = 10\,\%$
. This trend is particularly evident at
$x/h = 4.5$
, located on the trailing edge of the cylinder, where all three turbulent inflow cases yield notably similar results. One explanation is that the dominant characteristics of the separated shear layer outweigh the influence of the imposed inflow turbulence. Additionally, the reattachment of flow on the top surface of the cylinder may further contribute to this behaviour, as is discussed in the following sections.
Another key observation concerns the behaviour of the streamwise velocity component along the top surface of the cylinder. In the undisturbed case, the streamwise velocity at
$x/h = 4.25$
initially takes on negative values, indicating the presence of reverse flow. In contrast, this reversed flow is significantly weaker and more localised near the surface in the three disturbed cases, particularly those with higher turbulence intensities. These findings suggest that the introduction of turbulent inflow suppresses or eliminates the recirculation region above the top surface at
$x/h = 4.25$
and farther downstream. This effect is further illustrated in figures 13 and 14. In the last four lines within the wake region in figure 10, a noticeable difference in the streamwise velocity component is observed between the undisturbed and disturbed cases, which can be attributed to the wake-shrinking effect induced by inflow turbulence, as illustrated in figure 13. Furthermore, based on these profiles, increasing the TI at the inlet has a negligible impact on the time-averaged RMS of velocity components in the wake region. This observation may be attributed to the fact that the turbulence generated intrinsically by the cylinder model dominates over that introduced via the synthetic momentum source. Consequently, the turbulence characteristics in the wake are governed primarily by flow separation and geometry-induced vortex shedding, rather than by the externally imposed inlet disturbances.
Two components of non-dimensional vorticity in the disturbed case with a TI of
$10\,\%$
, visualised in two orthogonal planes: (a) spanwise vorticity (
$\omega _z$
) in the
$z$
-normal plane at
$z/h = 0$
and (b) vertical vorticity (
$\omega _y$
) in the
$y$
-normal plane at
$y/h = 0.8$
.

Time-averaged distributions of the non-dimensional velocity components and their RMS values along seven vertical lines located at
$x/h = 3.875$
,
$4.25$
,
$4.5$
,
$4.625$
,
$4.875$
,
$5.25$
and
$5.875$
on the
$z$
-normal plane at
$z/h = 0$
. The results correspond to the case with
$\text{TI} = 10\,\%$
and compare the baseline, where the source term is defined in the entire volume, and the limited source term region.

Time-averaged comparison of streamlines and non-dimensional streamwise velocity contours in a
$z$
-normal plane at
$z/h = 0$
: (a) undisturbed case; (b) disturbed case with
$\text{TI} = 10\,\%$
; (c) disturbed case with
$\text{TI} = 20\,\%$
; and (d) disturbed case with
$\text{TI} = 30\,\%$
.

As mentioned in the previous section, the stretching of vortices near the boundaries does not influence the near-field flow around the cylinder. To further confirm this, an additional simulation with a limited source term was performed. In this simulation, (A9) was applied in the
$y$
and
$z$
directions to define the source region within a confined volume. Specifically, the momentum source term was restricted to a region extending only
$2h$
above the cylinder and
$1.5h$
laterally from its sides, in contrast to the baseline case where the source term was applied across the entire domain boundaries. This simulation was conducted for the disturbed case with a TI of
$10\,\%$
. Figure 11 presents the instantaneous vorticity fields on two orthogonal planes,
$z_1 = 0$
and
$y_4 = 0.8h$
, non-dimensionalised by the cylinder height and reference velocity. Furthermore, figure 12 compares the time-averaged velocity components and their RMS values between the fully defined and restricted source term cases. The results demonstrate that vortices generated near the far-field boundaries and stretched by grid expansion have no discernible effect on the flow field in the vicinity of the cylinder.
Figure 13 presents a comparison of the time-averaged streamwise velocity component, non-dimensionalised by
$U_h$
, along with streamlines in the
$z/h = 0$
plane for all four cases. As shown in figure 13(b), corresponding to a TI of
$10\,\%$
, the recirculation region in the wake of the cylinder is noticeably reduced. A further reduction is observed between the
$10\,\%$
and
$20\,\%$
TI cases; however, the difference between the
$20\,\%$
and
$30\,\%$
cases is marginal, indicating a diminishing influence of additional turbulence beyond a certain threshold on the size of the recirculation region. The second observation is related to the region above the top surface of the cylinder, which is more clearly illustrated in figure 14. In figure 14(a), corresponding to the undisturbed inflow case, the separated shear layer detaches from the leading edge along the top surface of the cylinder and predominantly sheds into the wake region. Although the streamlines appear to connect to the trailing edge, it remains unclear whether the flow fully reattaches along the top surface. In contrast, the introduction of inflow turbulence leads to a clear reattachment along the top surface, as evidenced by the time-averaged streamlines. According to figure 14, two key trends are evident: first, the size of the recirculation region decreases with the introduction of inflow turbulence and continues to shrink as the TI increases; second, with increasing TI, the streamlines in figure 14(b–d) show progressively earlier reattachment along the top surface of the cylinder. These observations suggest that an inflow with higher TI leads to earlier reattachment and contributes to the stabilisation of the separated flow over the cylinder.
Time-averaged comparison of streamlines and non-dimensional streamwise velocity contours in a
$z$
-normal plane at
$z/h = 0$
: (a) undisturbed case; (b) disturbed case with
$\text{TI} = 10\,\%$
; (c) disturbed case with
$\text{TI} = 20\,\%$
; and (d) disturbed case with
$\text{TI} = 30\,\%$
.

To better understand the vortex topology and the influence of inflow turbulence on the flow dynamics around the cylinder, it is essential to examine the vorticity distribution around the cylinder. Figure 15 presents time-averaged, non-dimensionalised vorticity components on
$x$
-normal planes at
$x/h = 4.25$
,
$4.75$
and
$5.25$
, shown from left to right, respectively. The rotation associated with each vorticity component about the coordinate axes is also illustrated for the plane at
$x/h = 4.25$
, which lies at the mid-span of the cylinder in the streamwise direction. As indicated by the time-averaged distributions of
$\omega _x$
,
$\omega _y$
and
$\omega _z$
, the coherent structures exhibit rotation about the
$x$
,
$y$
and
$z$
axes, respectively, underscoring the inherently three-dimensional nature of their rotational topology. The time-averaged streamwise vorticity
$\omega _x$
reveals three distinct rotational regions, indicating the presence of coherent vortical structures. One of these rotational regions appears near the top surface of the cylinder and is associated with a Kelvin–Helmholtz-type structure generated by the separation of flow from the leading edge. This structure corresponds to the so-called secondary tip vortex previously identified by Rastan et al. (Reference Rastan, Shahbazi, Sohankar, Alam and Zhou2021). The second rotational region is associated with the tip vorticity, or primary tip vortex, as described by Rastan et al. (Reference Rastan, Shahbazi, Sohankar, Alam and Zhou2021), which forms above the sidewalls and is oriented inward towards the centreline of the cylinder. The third region corresponds to the vorticity concentrated near the base–wall junction, commonly referred to as the base vortex. Similar to the tip vortex, the base vortex also exhibits inward rotation towards the centre of the cylinder. The rotation of the tip and base vortices observed in the streamwise vorticity component
$\omega _x$
suggests that the downwash and upwash motions in the wake, commonly attributed to the tip and base vortices, are primarily governed by the
$x$
component of vorticity. Although
$\omega _z$
also exhibits large magnitudes, rotation about the
$z$
axis induces only localised, sectional upwash and downwash motions, and does not produce a net vertical (
$v$
) velocity across a bulk region of the flow. Moreover, analysis of the time-averaged results on the
$z$
-normal plane at
$z/h = 0.3$
, located
$0.05h$
from the cylinder sidewall, reveals that the upwash phenomenon persists continuously up to approximately
$x/h = 1.06$
in the wake region. However, the
$\omega _x$
and
$\omega _z$
components of vorticity associated with the base vortex retain their coherence only up to
$x/h = 1.5$
. This indicates that the disappearance of the upwash coincides with the loss of coherence in the base vortex. Moreover, in contrast to base vorticity and upwash phenomenon that cease to exist after a distance in the cylinder wake, the tip vorticity and downwash motion of flow persists up to the end of the computational domain. It also should be noted that the streamwise vorticity in figure 15 reveals that, between
$x/h = 4.25$
and
$x/h = 5.25$
, the relative contribution of the base vortex compared with the tip vortex is diminished.
Non-dimensional contours of vorticity components on the
$x$
-normal plane at
$x/h = 4.25$
,
$x/h = 4.75$
and
$x/h = 5.25$
from left to right, respectively.

To assess the impact of inflow turbulence on the coherent structures around the cylinder, figure 16 compares contours of
$\omega _x$
for four simulated cases. The comparison is made across the same three
$x$
-normal planes at
$x/h = 4.25$
,
$4.75$
and
$5.25$
. In the undisturbed case, the horseshoe vorticity is located near the cylinder corners at
$x/h = 4.25$
and extends beyond the plotted region by
$x/h = 5.25$
. The introduction of inflow turbulence shifts the vorticity inwards towards the cylinder’s shear layer. The figure also shows that the
$\omega _x$
component is weakened at
$x/h = 4.25$
in the presence of inflow turbulence, and this attenuation persists up to the case with
$\text{TI} = 30\,\%$
. However, based on the vorticity distributions at
$x/h = 4.75$
and
$x/h = 5.25$
, the reduction in
$\omega _x$
magnitude is not clearly observed across the disturbed cases, particularly in the two cases with higher turbulence intensities, which can be further explained by figure 17.
Non-dimensional
$x$
component of vorticity on the
$x$
-normal planes at
$x/h = 4.25$
,
$x/h = 4.75$
and
$x/h = 5.25$
for the undisturbed and three disturbed cases.

Illustration of reversed flow on the top and sidewalls of the cylinder across the undisturbed and three disturbed cases: (a) undisturbed case; (b) disturbed case with
$\text{TI} = 10\,\%$
; (c) disturbed case with
$\text{TI} = 20\,\%$
; and (d) disturbed case with
$\text{TI} = 30\,\%$
.

Figure 17 presents the volumetric reversed flow, defined as the volume of fluid with a streamwise velocity component less than zero, around the side and top walls of the cylinder across the four investigated cases. As shown in figure 17(a), corresponding to the undisturbed case, the flow over the top wall of the cylinder reattaches to the surface before entering the wake region. This reattachment was not clearly captured in figure 14, where time-averaged streamlines were plotted. It is noteworthy that the observed reattachment on the top surface of a cylinder with a side ratio of one contrasts with the findings of Rastan et al. (Reference Rastan, Shahbazi, Sohankar, Alam and Zhou2021), who reported no such reattachment. Given that the present study and that of Rastan et al. (Reference Rastan, Shahbazi, Sohankar, Alam and Zhou2021) differ in terms of Reynolds number and magnitude of streamwise velocity over the top surface, this suggests that the occurrence of reattachment in such configurations may depend on additional flow parameters beyond the side ratio, including Reynolds number and inflow velocity. Furthermore, although a localised reattachment occurs near the corner of the sidewall, there is no reattachment in the majority of the shear layer along the sidewall and it directly sheds into the wake region behind the cylinder.
Figure 17(b) corresponds to the disturbed case with
$\text{TI} = 10\,\%$
. The flow shows that reattachment on the top surface occurs earlier compared with the undisturbed case, and this reattachment shifts further upstream with increasing TI. This behaviour is consistent with the reduction in the size of the separation bubble, as shown in figure 14. In addition to the earlier reattachment of the shear layer on the top surface, the introduction of inflow turbulence also leads to flow attachment along the sidewalls. In the disturbed cases, the shear layer on the sidewalls no longer sheds into the wake region but instead remains attached to the surface. This reattachment is consistent with the findings of Cao et al. (Reference Cao, Liu, Zhou and Ren2022), who investigated the flow dynamics of a square cylinder with an aspect ratio of 3 under disturbed inflow conditions. This reattachment occurs slightly earlier in the disturbed case with
$\text{TI} = 20\,\%$
; however, the shear layers along the sidewalls in the cases with
$\text{TI} = 20\,\%$
and
$\text{TI} = 30\,\%$
appear qualitatively similar.
The fact that shear layers along the sidewalls in the two disturbed cases with
$\text{TI} = 20\,\%$
and
$\text{TI} = 30\,\%$
appear qualitatively similar suggests that increasing the TI from
$20\,\%$
to
$30\,\%$
has no significant additional impact on the side-wall shear layer dynamics. Consequently, the
$\omega _x$
contours at the downstream planes of
$x/h = 4.75$
and
$x/h = 5.25$
remain largely unchanged between these two disturbed cases. It should be noted that, according to figures 8 and 10, the TI in the case with
$\text{TI} = 30\,\%$
is indeed higher; however, this increase does not appear to exert a noticeable influence on the side shear layer. The underlying cause of this limited sensitivity remains to be investigated in future work.
3.3. Instantaneous results
In the previous section, the procedure for generating inflow turbulence using a volume momentum source term was explained. Additionally, the effect of inflow turbulence at various intensity levels was analysed. The primary objective of this section is to examine the influence of incoming turbulence and its intensity on the instantaneous flow field, with a particular focus on coherent turbulent structures. To achieve this, various post-processing techniques are employed, including the vortex identification criterion (Jeong & Hussain Reference Jeong and Hussain1995)
$\lambda _2$
, the instantaneous pressure coefficient and SPOD. These tools enable the identification and characterisation of flow structures and coherent turbulence features within the wake of the cylinder.
Figure 18 presents the instantaneous vortex structures identified using the
$\lambda _2$
criterion, as defined by Jeong & Hussain (Reference Jeong and Hussain1995), for the undisturbed case and the three disturbed cases. The vortices are coloured by the instantaneous streamwise velocity, non-dimensionalised with respect to the reference velocity
$U_h$
. Figure 18(a) corresponds to the undisturbed case. As shown, the flow separates almost immediately upon encountering the top and side surfaces of the cylinder. The vortices that form along the sidewalls exhibit a forward-bending behaviour, which is attributed to the non-uniform inflow velocity profile. On the top surface, similar vortices are clearly visible. These vortical structures on the side and top surfaces of the cylinder originate from the flapping motion of the separated shear layers and are indicative of instabilities caused by strong velocity gradients near the surface. As is shown in the SPOD analysis section, these two structures occur at the same frequency and are associated with Kelvin–Helmholtz instability on top and side surfaces. It should be noted that the Kelvin–Helmholtz instability on the top surface is referred to as secondary tip vortices (STVs), consistent with the observations and terminology introduced in the studies by Rastan et al. (Reference Rastan, Shahbazi, Sohankar, Alam and Zhou2021) and Maleki et al. (Reference Maleki, Mankbadi and Golubev2025).
Isosurface of the instantaneous
$\lambda _2d^2/U^2_{h}=-30$
criterion coloured by the streamwise velocity component, non-dimensionalised with respect to
$U_h$
: (a) undisturbed case; (b) disturbed case with
$\text{TI} = 10\,\%$
; (c) disturbed case with
$\text{TI} = 20\,\%$
; and (d) disturbed case with
$\text{TI} = 30\,\%$
.

Figure 18(b) shows the instantaneous
$\lambda _2$
field for the case with
$\text{TI}=10\,\%$
. In this case, the vortical structures on top and sidewall are still present and distinguishable, although their characteristic forward-bending behaviour is less apparent compared with the undisturbed case. Moreover, unlike the undisturbed case, these vortices are unable to maintain their coherent form downstream. This instantaneous visualisation of the
$\lambda _2$
field illustrates that the imposed inflow turbulence disrupts the organised shear layer, resulting in the earlier breakdown of the Kelvin–Helmholtz vortices. Figures 18(c) and 18(d) show the instantaneous
$\lambda _2$
fields for turbulence intensities of 20 % and 30 %, respectively. As evident from these illustrations, increasing the TI leads to a significantly more chaotic flow field. While the STV structures remain partially discernible in the case with
$\text{TI} = 20\,\%$
, they become increasingly disorganised. In the case with
$\text{TI} = 30\,\%$
, the high level of turbulent fluctuations obscures coherent structures, making it difficult to extract meaningful information about the organised vortex dynamics. It should be noted that, due to the inherently unsteady nature of flow structures, a single instantaneous snapshot cannot fully capture their overall behaviour. The evolution of coherent structures is better understood by observing their development over one or more characteristic periods. Alternatively, advanced post-processing techniques such as SPOD can be employed to extract dominant energetic modes and provide a more comprehensive representation of the underlying flow dynamics.
In this regard, a sequence of six instantaneous snapshots, shown in figure 19, qualitatively captures the evolution of the STV structures forming on the top surface of the cylinder in the undisturbed case. The snapshots are taken at equal time intervals with
$\Delta t=0.001$
s, allowing visualisation of the temporal development of the structures, even though the vortex-shedding period has not been explicitly calculated. The data are extracted on a plane located at
$z/h = 0$
, and the pressure coefficient is computed using the relation
$C_p = {p - p_{\infty }}/{0.5 \, \rho \, U_{h}^2}$
, and the contours of
$C_p$
range from
$-2$
to
$1$
. According to this figure, the structure originates at the leading edge of the cylinder and progressively evolves and grows in size throughout the sequence. The structure remains trackable up to the mid-span of the cylinder; however, the results in this figure as well as additional snapshots at different phases of its evolution indicate that the structure does not retain a coherent form throughout its entire life cycle. Typically, the vortex interacts with preceding vortices, leading to pairing or merging events near the top surface. In some instances, the structure may reappear near the trailing edge of the cylinder. This behaviour of STV structures over the top surface of cylindrical models with a free end, including their pairing or merging with increasing distance from the leading edge, has also been reported in previous studies (Lander et al. Reference Lander, Moore, Letchford and Amitay2018; Moore, Letchford & Amitay Reference Moore, Letchford and Amitay2019). In the undisturbed case, the STV structure exhibits fully periodic behaviour throughout the flow domain, whereas, as will be shown for the disturbed cases with inflow turbulence, the structure becomes irregular and less coherent. As a final note on the pressure coefficient distribution in the undisturbed case, the above statements on the life cycle and persistence of STV structures only up to the cylinder mid-span should be interpreted with caution, as they are based solely on the
$C_p$
distribution post-processing method. These observations are therefore not definitive, and further analyses such as spectral and SPOD methods have been employed to more reliably evaluate their evolution.
Instantaneous pressure coefficient in (a–f) six snapshots at the
$z$
-normal plane,
$z/h = 0$
, for the undisturbed case.

Figures 20–22 present the pressure coefficient distributions for the disturbed cases with
$\text{TI} = 10\,\%$
,
$\text{TI} = 20\,\%$
and
$\text{TI} = 30\,\%$
, respectively. The range of
$C_p$
values was kept consistent across all cases to enable direct comparison. Additionally, the time interval between successive snapshots was increased from
$10\Delta t$
to
$25\Delta t$
to capture a broader temporal evolution of the flow over the top surface of the cylinder. The first notable observation from these figures is the increase in the intensity of the pressure variations compared with the undisturbed case. With the introduction of turbulence at the inflow, the flow structures become more intensified in terms of the local pressure gradients they generate. This is evident through stronger and more localised variations in
$C_p$
.
Instantaneous pressure coefficient in (a–f) six snapshots at the
$z$
-normal plane,
$z/h = 0$
, in the disturbed case with
$\text{TI}=10\,\%$
.

The second observation concerns the increasingly chaotic behaviour of the flow over the top surface of the cylinder. Unlike the undisturbed case, where the STV is generated and convected downstream in a highly periodic manner, the flow in the three disturbed cases behaves more randomly and chaotically. In some snapshots, particularly in figure 20, corresponding to the disturbed case with
$\text{TI} = 10\,\%$
, it is still possible to observe remnants of periodic STV behaviour. However, in other snapshots from this case, the structures lose their coherence and deviate from the regular paths seen in the undisturbed case. As the TI increases, the initial formation of structures at the leading edge remains observable; however, their subsequent evolution becomes increasingly unpredictable. In some instances, the structures merge with neighbouring vortices, forming larger and more complex patterns, while in several instances, they dissipate and vanish before being convected downstream.
By combining the information from the time-averaged streamwise velocity shown in figure 14 with the instantaneous pressure coefficient distributions presented in figures 19–22, it can be concluded that the STVs continue to form from the leading edge even in the presence of inflow turbulence. However, they lose their periodic nature under disturbed conditions. Furthermore, as shown in figures 14 and 17, the separated shear layer becomes smaller in the disturbed cases, and its size continues to decrease with increasing TI. This trend suggests that, over longer simulation times, the STV structures are more prone to dissipate quickly in disturbed cases compared with the undisturbed case. Nevertheless, concluding that these structures cease to exist would be inaccurate, as they are clearly still observable in figures 20–22.
Instantaneous pressure coefficient in (a–f) six snapshots at the
$z$
-normal plane,
$z/h = 0$
, in the disturbed case with
$\text{TI}=20\,\%$
.

Instantaneous pressure coefficient in (a–f) six snapshots at the
$z$
-normal plane,
$z/h = 0$
, in the disturbed case with
$\text{TI}=30\,\%$
.

3.4. Spectral analysis of velocity vector
As mentioned earlier, SPOD post-processing tools were also employed in this study. As summarised in table 1, SPOD was applied to 17 planes to decompose the velocity field into coherent turbulent structures characterised by distinct frequencies, energy content and spatial distributions. Figure 23 illustrates the SPOD spectral energy distribution of the undisturbed case at the
$z_1$
plane within Domain 1, which extends
$8d$
downstream and
$1d$
above the cylinder. The distribution is presented in terms of the modal eigenvalues
$\lambda$
, normalised by the total energy of all 28 modes in the entire frequency bands, and plotted as a function of the Strouhal number. The Strouhal number is defined as the modal frequency non-dimensionalised with respect to the cylinder width (
$ 0.1$
m) and the streamwise velocity at the top of the cylinder (
$ 2.89$
m s−1). In the SPOD analysis, the modes are ranked according to the energy they contain, with the first mode being the most energetic (the dominant mode) and the second mode being the second most energetic. Moreover, the SPOD framework suggests that, at a given frequency, if the dominant mode or several modes have substantially larger energy than the others, it signifies the presence of one or more physically coherent structures associated with that frequency (Zhang et al. Reference Zhang, Ooka and Kikumoto2021). In this context, if only the dominant mode exhibits energy at a particular frequency, the structure demonstrates single-mode behaviour. Conversely, if multiple modes exhibit energy at the same frequency, the structure is classified as exhibiting multi-mode behaviour.
Modal energy spectra in the
$z_1$
plane from SPOD analysis of the undisturbed case in Domain 1.

According to figure 23, the dominant mode exhibits a spectral peak within the range
${\textit{St}} = 0.084{-}0.101$
. For comparison, Zhang et al. (Reference Zhang, Ooka and Kikumoto2021) reported a non-dimensional shedding frequency of
${\textit{St}} = 0.0845$
for a similar cylinder configuration. Furthermore, da Silva et al. (Reference da Silva, Hahn, Sumner and Bergstrom2022) reported a comparable value of
${\textit{St}} = 0.085$
in an experimental study of low-AR square cylinders at
$Re = 7.5\times 10^4$
, where the cylinder was immersed within a thick boundary layer. Building on that work, da Silva, Sumner & Bergstrom (Reference da Silva, Sumner and Bergstrom2024) performed LES for a low-AR cylinder at
$Re = 1\times 10^4$
and similarly identified a dominant shedding frequency of
${\textit{St}} = 0.086$
under thick-boundary-layer conditions. In the present results, although the mode at
${\textit{St}} = 0.084$
carries slightly higher energy than that at
${\textit{St}} = 0.101$
, neither aligns exactly with the dominant frequency; rather, the true spectral peak is situated between these two values. In this regard, two points should be addressed. First, the spatial extent and structure of the dominant mode (first mode) at both
${\textit{St}} = 0.084$
and
${\textit{St}} = 0.101$
appear similar across three velocity components (
$u$
,
$v$
and
$w$
), suggesting that two dominant frequencies are not independent structures and the energy of flow is just divided between these two frequencies.
Second, the ability to accurately resolve the dominant frequency is constrained by the spectral resolution of the SPOD analysis. According to the relation
$\Delta f = {f_s}/{N_{{\textit{FFT}}}}$
, described in § 2.5, the frequency resolution
$\Delta f$
in SPOD is governed by the sampling frequency
$f_s$
and the number of samples per block
$N_{{\textit{FFT}}}$
. Assuming a fixed sampling frequency, improving the resolution requires increasing the number of samples per block. This can be accomplished either by doubling the total sampling duration or by reducing the number of blocks, since
$N_{{\textit{FFT}}}$
must be a power of two. However, the first approach substantially increases the computational cost, particularly in this study, which involves simulations of four different cases. As for the second approach, preliminary results indicate that using only 14 blocks is insufficient for the disturbed cases and does not yield converged SPOD results. Therefore, a choice of 28 blocks was adopted to strike a balance between computational cost and sufficient frequency resolution for accurate modal identification.
Figure 24 illustrates the spatial distribution of the three velocity components on three orthogonal planes –
$z_1$
,
$y_4$
and
$x_5$
– at the non-dimensional frequency
${\textit{St}} = 0.084$
. This figure highlights how each velocity component is distributed around the cylinder. As observed in this figure, the
$u$
and
$v$
components exhibit antisymmetric distributions, primarily originating from the shear layers along the lateral surfaces of the cylinder; however, the
$w$
component displays a symmetric pattern. Moreover, as observed from the distribution of the
$w$
component in the
$z_1$
plane, the structure exhibits a forward-bending behaviour. A similar tendency is also evident in the
$u$
and
$v$
components in the
$z_3$
plane, as shown in figure 36. This behaviour is attributed to the non-uniform velocity profile in the vertical direction. As the flow velocity increases with height above the ground, the upper portion of the structure is convected downstream more rapidly than the lower portion. This velocity gradient results in a forward bending of the overall structure.
Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.084$
in the undisturbed case. From left to right,
$z_1$
,
$y_4$
and
$x_5$
planes.

Maleki et al. (Reference Maleki, Mankbadi and Golubev2025) also conducted a similar analysis around a high-AR (7) wall-mounted cylinder using multiple planes. At the dominant vortex shedding frequency, the spatial distribution of the
$u$
and
$v$
velocity components shows qualitative similarity between the two studies; however, two key differences are observed. First, unlike the high-AR case in Maleki et al., where cellular vortex shedding was identified at the dominant frequency, the lower-AR configuration investigated in the present study does not exhibit any such cellular shedding patterns, consistent with physical expectations (Moreau & Doolan Reference Moreau and Doolan2013; Porteous et al. Reference Porteous, Moreau and Doolan2017). Furthermore, as mentioned earlier, the dominant modes at both
${\textit{St}} = 0.084$
and
${\textit{St}} = 0.101$
display similar spatial characteristics, as illustrated for
${\textit{St}} = 0.101$
in figure 37, further supporting the conclusion that these high-energy modes are not associated with cellular vortex shedding. Second, while the present study reveals a forward-leaning tendency in the vortex structure under non-uniform velocity profile, the analysis of Maleki et al. showed a backward-leaning structure under uniform velocity profile. This contrast underscores the influence of AR and vertical velocity gradients on the orientation and convection of coherent flow structures. To investigate the behaviour of the main vortex shedding in the streamwise direction, figure 25 illustrates the spatial distribution of the three velocity components at the non-dimensional frequency of
${\textit{St}} = 0.084$
across three streamwise-oriented planes:
$x_2$
,
$x_3$
and
$x_5$
. According to this figure, the
$u$
and
$v$
components of velocity originate from the shear layers along the sides and top of the cylinder. As the flow progresses from plane
$x_2$
to
$x_5$
, their spatial distribution extends across the full height of the cylinder, although the
$v$
component remains primarily concentrated near the top region. In contrast, the
$w$
component is initially concentrated near the top of the cylinder at plane
$x_2$
, but by plane
$x_5$
, it becomes distributed across both the top and the junction region between the cylinder and the ground.
According to figure 23, another peak observed in the modal analysis of the
$z_1$
plane occurs at a non-dimensional frequency of
${\textit{St}} = 2.4$
, and figure 26 presents the spatial distribution of the velocity components across the
$z_1$
,
$y_7$
and
$x_2$
planes at this particular frequency. This figure reveals the development of a coherent structure originating from the top leading edge of the cylinder, which is consistent with the pressure coefficient results shown in figures 19–22. It is evident that this coherent structure, at the identified frequency, is associated with fluctuations in the separated shear layer at the top of the cylinder, giving rise to a Kelvin–Helmholtz-type instability, named STV. According to figure 26, the structure at this frequency extends up to approximately
$d/2$
downstream of the leading edge of the cylinder, beyond which it loses coherence.
Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.084$
in the undisturbed case. From left to right,
$x_2$
,
$x_3$
and
$x_5$
planes.

Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=2.4$
in the undisturbed case. From left to right,
$z_1$
,
$y_7$
and
$x_2$
planes.

To better analyse the behaviour and energy of STV structure and due to the influence that coherent turbulent structures may have on urban air mobility (Hua et al. Reference Hua, Afari, Golubev and Mankbadi2025), Domain 2 was designated for SPOD analysis to specifically capture the coherent structures present near the top of the cylinder. Figure 27 presents the modal energy spectra in the
$z_1$
plane for Domain 2, and the comparison of the first SPOD modes in the
$z_1$
plane between Domain 1 and Domain 2. Domain 2 extends 3
$d$
downstream of the cylinder and 1
$d$
above the cylinder tip. According to the results, since this SPOD domain does not encompass the primary region of the main vortex shedding, the energy of the first mode at
${\textit{St}} = 0.084$
is substantially diminished. Moreover, as is shown in figure 27(b), the energy of the first mode in the Strouhal number range of approximately 0.3–1.1 increases. A similar change in the SPOD results of Domain 1 and Domain 2 was also reported by Maleki et al. (Reference Maleki, Mankbadi and Golubev2025), though it was not discussed in detail. In this context, it is important to note that SPOD captures the energy contributions of all coherent structures across the entire plane on which it is conducted. As such, it may not be the most suitable approach for analysing spatial variations in dominant frequency within a localised region or along a specific direction. Therefore, to enable a more localised investigation of the flow near the top surface of the cylinder, particularly to analyse the behaviour of flow structures within the Strouhal number range of
${\textit{St}} \in [0.3,\,1.1]$
, and to track the evolution of the dominant frequency from the leading edge to the trailing edge in the streamwise direction, a fast Fourier transform should be employed to complement the SPOD results.
(a) Modal energy spectra in the
$z_1$
plane based on SPOD analysis of the undisturbed case, in Domain 2. (b) Comparison of the first SPOD modes of the undisturbed case in the
$z_1$
plane between Domain 1 and Domain 2.

To do so, four probe points were defined on the top surface of the cylinder, positioned outside the separated shear layer, as shown in figure 28. These locations were selected based on the recommendations of Lander et al. (Reference Lander, Moore, Letchford and Amitay2018) and Moore et al. (Reference Moore, Letchford and Amitay2019), who suggest placing probes outside the shear layer to more accurately capture the frequency of coherent structures. On the right-hand side of figure 28, a fast Fourier transform using the Welch method (Hamming Reference Hamming2013) was applied to a 5 s signal of the wall-normal velocity component. The spectral density amplitudes have been vertically offset to facilitate comparison of the energy distribution among the probe points. The first probe point, located at
$(x/h, y/h) = (4.05, 1.05)$
, was positioned near the leading edge of the cylinder’s top surface. The spectral analysis at this location reveals a dominant peak at a non-dimensional frequency of
${\textit{St}} = 2.4$
, consistent with the SPOD results presented earlier. Moreover, the agreement between the spatial distribution of the SPOD mode at
${\textit{St}} = 2.4$
which shows a concentration of coherent structures near the leading edge (figure 26) and the position of probe point (1), supports the conclusion that this frequency corresponds to fluctuations, or flapping motion of the separated shear layer at the top of the cylinder.
(a) Time-averaged streamwise velocity, non-dimensionalised by
$U_h$
, with the positions of four probe points located outside the separated shear layer. (b) The PSD of the wall-normal velocity component at the selected probe points. The spectral density amplitudes of probe points 2, 3 and 4 have been vertically offset by factors of 10, 100 and 1000, respectively, to facilitate comparison of the energy distribution among the probe points.

Probe points (2) to (4) were placed progressively farther downstream. As shown in figure 28, the spectral energy associated with
${\textit{St}} = 2.4$
decreases with increasing distance from the leading edge. At probe point (4), the dominant spectral energy becomes distributed over a broader frequency range, approximately
${\textit{St}} \sim [0.3,\,1.2]$
, which aligns with the SPOD results obtained in Domain 2 (figure 27). This spectral behaviour observed in figures 27 and 28 is attributed to vortex pairing occurring along the top surface of the cylinder. In summary, these STVs are initially generated by the flapping motion of the separated shear layer and, as they convect downstream, they undergo pairing and merge into larger-scale vortices with longer wavelengths, which correspond to the lower-frequency content seen in the spectra. This interpretation is consistent with the findings of Lander et al. (Reference Lander, Moore, Letchford and Amitay2018) and Moore et al. (Reference Moore, Letchford and Amitay2019).
It should be noted that the SPOD results also support the occurrence of the vortex pairing phenomenon. Figure 29 illustrates the spatial distribution of two velocity components in the
$z_1$
plane at two Strouhal numbers,
${\textit{St}} = 0.625$
and
${\textit{St}} = 1.098$
, within the frequency range of
${\textit{St}} = 0.3{-}1.1$
. According to this figure, coherent structures similar to the STV structures shown in figure 26 at
${\textit{St}} = 2.4$
are detected. However, in contrast to the limited spatial extent of the first mode at
${\textit{St}} = 2.4$
, where no coherence was observed beyond
$d/2$
downstream from the leading edge, the structures in the range of
${\textit{St}} = 0.3{-}1.1$
shown in figure 29 retain their coherence and continue to shed into the wake region. Furthermore, as expected, the size and wavelength of the STV structures within this frequency range are larger than those of the STV structure at
${\textit{St}} = 2.4$
, which is consistent with the vortex pairing phenomenon observed on the top surface of the cylinder.
Spatial distribution of mode 1 for two velocity components at
${\textit{St}} = 0.625$
and
${\textit{St}} = 1.098$
in the undisturbed case and in the
$z_1$
plane.

Before analysing the SPOD results for the disturbed cases, one final point regarding the undisturbed case should be addressed. In figure 18(a), specifically in the instantaneous results section, it was mentioned that the coherent structures observed along the side surfaces of the cylinder resemble Kelvin–Helmholtz-type structures, similar to those on the top surface. Figure 30 presents the SPOD analysis in the
$y_4$
plane for three configurations: both sides of the cylinder, the positive side only and the negative side only. As seen in figure 30(a), a peak appears at the non-dimensional frequency of
${\textit{St}} \approx 2.4$
, consistent with the SPOD results in the
$z_1$
plane shown in figures 23 and 27. However, the results in the
$y_4$
plane indicate that two modes contribute to the formation of this structure when both sides of the cylinder are included in the SPOD domain. But when SPOD is performed on either the positive or negative side alone, the energy is concentrated in a single mode. This behaviour is attributed to the presence of Kelvin–Helmholtz-type structures occurring symmetrically on both sides of the cylinder. When the full domain is considered, SPOD captures both structures, each exciting a separate mode. In contrast, applying SPOD to only one side captures a single Kelvin–Helmholtz structure, thereby exciting only one dominant mode. The spatial extent of these structures along the side surfaces is similar to those observed on the top surface; therefore, further presentation of the results is considered redundant.
Modal energy spectra in the
$y_4$
plane based on SPOD analysis for undisturbed cases: (a) domain extended on both sides of the cylinder; (b) domain extended on the positive side; and (c) domain extended on the negative side.

To analyse the effect of incoming turbulence on the flow around the cylinder, SPOD post-processing was performed for all three disturbed inflow cases, with data sampled at 17 planes in each case, as summarised in table 1. The resulting SPOD spectral energy distributions in the
$z_1$
plane for turbulence intensities of 10 %, 20 % and 30 % are presented in figure 31. According to this figure, introducing a turbulent inflow leads to a reduction in the energy associated with the dominant vortex shedding. Nevertheless, the primary shedding frequency remains detectable. A second observation pertains to the cases with 20 % and 30 % TI where, as the intensity increases, the dominant peak becomes broader and less sharply defined, appearing more as a bump in the energy spectrum. In contrast to the undisturbed case shown in figure 23, where energy was concentrated around one or two distinct frequencies, the higher-intensity disturbed cases exhibit a broader energy distribution, with energy spread across a wider band of frequencies.
Modal energy spectra in the
$z_1$
plane based on SPOD analysis of disturbed cases: (a)
$\text{TI}=10\,\%$
; (b)
$\text{TI}=20\,\%$
; (c)
$\text{TI}=30\,\%$
.

Figure 32 illustrates the spatial distribution of the first SPOD mode for the three velocity components at
${\textit{St}} = 0.084$
in the disturbed case with
$10\,\%$
TI. This figure shows that although the introduction of disturbed inflow affects the velocity field, the overall shape and spatial extent of the modal structures are largely preserved. However, under turbulent inflow conditions, a single specific frequency may no longer be sufficient to represent the dominant coherent structures across all velocity components. This is consistent with the observation that, at 20 % and 30 % TI, the dominant spectral peak appears more as a broad bump than a sharp peak. Overall, by comparing this result with figure 24, which illustrates the corresponding mode in the undisturbed case, it can be concluded that while the energy and coherence of the flow structures are influenced by the inflow turbulence, their general shape remains qualitatively consistent.
Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.084$
in the disturbed case with
$\text{TI}=10\,\%$
. From left to right,
$z_1$
,
$y_4$
and
$x_5$
planes.

The final observation drawn from figure 31 concerns the frequency of fluctuations in the separated shear layer, which gives rise to the STV structures at the top of the cylinder. As shown in the figure, the introduction of turbulent inflow reduces the energy captured in the SPOD spectra for the disturbed case with
$10\,\%$
TI at
${\textit{St}} = 2.4$
; however, SPOD is still able to detect coherent energy at this frequency. It should be recalled that the
$10\,\%$
TI case successfully reproduced the required velocity fluctuations to mimic the atmospheric conditions observed in the experimental data. The fact that SPOD continues to capture energy at the non-dimensional frequency of
${\textit{St}} = 2.4$
under turbulent inflow conditions suggests that the flapping motion of the separated shear layer retains its temporal and spatial coherence even in the presence of inflow turbulence. This finding stands in contrast to the recent results of Chen et al. (Reference Chen, Li, Xue, He, Wang and Liang2023), in which no SPOD energy was detected under similar
$10\,\%$
TI conditions. These results emphasise the critical importance of carefully selecting both SPOD parameters and simulation set-up to ensure the accurate detection of coherent structures in turbulent environments.
Regarding the cases with higher TI, according to figure 31 no significant energy peak is observed in the two higher-intensity cases. As previously discussed for the instantaneous pressure coefficient distributions (figures 21 and 22), the introduction of turbulence at the inflow disrupts the formation and behaviour of the STVs near the top of the cylinder. The added turbulence renders these structures less periodic, altering them into more gust-like, random phenomena rather than a predictable sequence of coherent vortices. This loss of periodicity and increased randomness reduces their spatial and temporal coherence. Since SPOD is inherently based on extracting energy from structures that exhibit spatial and temporal coherence, these modified, less-organised flow features are not strongly detected in the SPOD analysis. However, it should be noted that although there are still phenomena related to the STVs and fluctuations in the separated shear layer, these structures or phenomena no longer exhibit the characteristics of classical Kelvin–Helmholtz-type instabilities.
4. Summary and conclusion
This study presents a high-fidelity numerical investigation of the wake dynamics behind a WMSC with an AR of
$ 2$
, subjected to varying inflow turbulence intensities. LES was employed in conjunction with SPOD to analyse the unsteady wake behaviour. A non-uniform inlet velocity profile with a reference velocity of
$U_h = 2.89\,\mathrm{m}\mathrm{s}^{-1}$
was imposed, corresponding to a Reynolds number of
$Re_h = 43{\,}000$
. To emulate realistic atmospheric turbulence, a momentum-based turbulence source term was coupled with the STG method. The resulting inflow field was validated by comparing the PSD of the velocity components at the probe location with the targeted von-Kármán spectrum, showing good agreement. Four inflow conditions were studied: one undisturbed case and three disturbed cases with turbulence intensities of 10 %, 20 % and 30 %. Time-averaged velocity fields were compared with experimental measurements from Okaze et al. (Reference Okaze, Kikumoto, Ono, Imano, Ikegaya, Hasama, Nakao, Kishida, Tabata and Nakajima2021).
The comparison between the undisturbed and disturbed inflow cases revealed a significant influence of incoming turbulence on the vortex-shedding behaviour. Notably, the inflow turbulence reduced the size of the recirculation region in the cylinder wake, promoted earlier flow reattachment on the top surface and induced reattachment along the sidewalls, a phenomenon absent in the undisturbed configuration. In contrast to the findings of Rastan et al. (Reference Rastan, Shahbazi, Sohankar, Alam and Zhou2021), who reported no reattachment of the top shear layer for a cylinder with a side ratio of one under undisturbed inflow, the time-averaged reversed flow in the present undisturbed case indicates a clear reattachment. It also should be noted that with the introduction of turbulence at the inflow, the onset of reattachment occurred earlier, and increasing TI further facilitated this process. In addition, the inflow turbulence altered the development of the sidewall shear layers. In the undisturbed case, the time average of reversed flow showed no reattachment along the sidewalls. However, with a 10 % TI, reattachment occurred, and increasing the intensity to 20 % further accelerated this process. Nevertheless, increasing the turbulence to 30 % did not lead to a substantial additional effect on the sidewall reattachment. Vorticity field analysis demonstrated that inflow turbulence suppressed the formation of coherent tip and base vortices. In particular, the streamwise vorticity component (
$ \omega _x$
) exhibited weakened structures, indicating dissipation of vortices responsible for the characteristic upwash and downwash dynamics. Furthermore, the distribution of pressure coefficient in the undisturbed case revealed periodic Kelvin–Helmholtz structures along the top and side surfaces. The introduction of inflow turbulence disrupted the formation of these structures. While the
$10\,\%$
case still exhibited identifiable periodic features, the
$20\,\%$
and
$30\,\%$
turbulence-intensity cases showed a marked degradation in both spatial and temporal coherence, owing to the high turbulence level in the incoming inflow, and transitioned towards irregular, gust-like flow without clear periodicity.
A SPOD analysis, based on 30 000 velocity snapshots, was performed on 17 orthogonal planes aligned with the
$ x$
,
$ y$
and
$ z$
axes. The analysis revealed two dominant energy frequencies: a primary peak at the non-dimensional frequency of
$ St = 0.084$
, associated with dominant vortex shedding, and a secondary peak at
$ St = 2.4$
, linked to Kelvin–Helmholtz structures forming along the top and side surfaces of the cylinder. The modal energy distributions showed that the dominant shedding mode at
${\textit{St}}=0.084$
exhibits antisymmetric distribution in the streamwise and wall-normal velocity components, and symmetric distribution in the spanwise component. The SPOD was independently conducted on two computational domains, and comparison revealed consistent trends in modal energy content. However, in Domain 2, a frequency band within the range
$ 0.3 \leq St \leq 1.1$
exhibited elevated energy content relative to Domain 1. Fast Fourier transform analysis of probe points located outside the top shear layer confirmed that this frequency band is associated with vortex pairing of Kelvin–Helmholtz structures on the top surface. Moreover, the distribution of SPOD modes within this band resembled those at
$ St = 2.4$
, further supporting this association. In this research it was demonstrated that the frequency
$ St = 2.4$
, observed in both domains, corresponds to the flapping motion of the separated shear layer on the top surface, and that the enhanced energy within the
$ 0.3 \leq St \leq 1.1$
band originates from vortex-pairing phenomena in this region. The SPOD results also revealed that, among the disturbed cases, only the case with 10 % inflow TI retained periodic Kelvin–Helmholtz structures. Higher turbulence intensities resulted in a significant loss of spatial and temporal coherence, suppressing the classical Kelvin–Helmholtz signature. However, although introducing inflow turbulence reduces the energy associated with the dominant vortex shedding at
$ St = 0.084$
, the frequency content at this Strouhal number remains preserved.
Funding
This research is based upon work supported by the NASA Aeronautics Research Mission Directorate (ARMD) University Leadership Initiative (ULI) under cooperative agreement number 80NSSC21M0057. Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Aeronautics and Space Administration. All simulations in this work were conducted using ERAU’s Vega HPC system.
Declaration of interests
The authors report no conflict of interest.
Data and code availability
The modified OpenFOAM solver and the corresponding case files used in this study are publicly available at OpenFoam-Momentum-Source-Turbulence. The repository contains the full implementation of the momentum source formulation for generating synthetic atmospheric turbulence, as well as the OpenFOAM case set-up used to reproduce the wall-mounted cylinder simulations presented in this work.
Appendix A
Without loss of generality, the momentum source generating a convected incompressible flow disturbance represented by (2.7) is applied in a region of uniform flow aligned with the
$x_1$
direction. The interactions between the perturbation and any other waves are assumed to be negligible. We now prove that the solution for the momentum source terms in the equations
producing the disturbance downstream of the region
$|x_1 {-} x_{1s}| \leq {\pi }/{b}$
of the computational domain can be obtained in the form
\begin{equation} S_i(x_{\!j}, t) = \begin{cases} \sum \limits _{n=1}^{\infty } \omega _n K_n g(x_1) \left ( b_i^n \sin \sigma _n - a_i^n \cos \sigma _n \right ) & \text{for } i = 1, \\[10pt] \sum \limits _{n=1}^{\infty } U_\infty K_n g'(x_1) \left ( b_i^n \cos \sigma _n + a_i^n \sin \sigma _n \right ) & \text{for } i = 2,3, \end{cases} \end{equation}
where
$\sigma _n = \omega _n t - k_2^n x_2 - k_3^n x_3 - k_1^n x_{1s}$
,
$g'(x_1) = \mathrm{d}g/\mathrm{d}x_1$
and
$K_n$
is given by
\begin{equation} K_n = \frac {\left (k_1^n\right )^2 - b^2}{b^2 \sin \left (\dfrac {\omega _n \pi }{b U_\infty }\right )}. \end{equation}
The function
$g(x_1)$
is given by
\begin{equation} g(x_1) = \begin{cases} \dfrac {1}{2} \left ( 1 + \cos [b(x_1 - x_{1s})] \right ), & |x_1 - x_{1s}| \leq \dfrac {\pi }{b}, \\[8pt] 0, & |x_1 - x_{1s}| \gt \dfrac {\pi }{b}. \end{cases} \end{equation}
Since (A2) represents a linear superposition of modes, we restrict the analysis to a single mode
$n$
from (A2) that produces the corresponding linear response in (A1). The analysis is then conducted in the complex plane where the modes of the velocity disturbance take the form
where
$\hat {u}_i^n (x_1)$
denote modal complex amplitudes,
$\psi _n = \omega _n t - k_2^n x_2 - k_3^n x_3$
and
$i^2 = -1$
. Substituting
$S_i(x_{\!j}, t)$
from (A2), (A1) reduces to
with
\begin{equation} \tilde {S}_i^n (x_1) = \begin{cases} -k_1^n g(x_1) K_n (a_i^n + i b_i^n) \exp (-i k_1^n x_{1s}), & \text{for } i = 1, \\ g'(x_1) K_n (b_i^n - i a_i^n) \exp (-i k_1^n x_{1s}) , & \text{for } i = 2, 3. \end{cases} \end{equation}
Formal solution to (A6) is of the form
Integrating (A8) and substituting the result in (A5), the perturbation velocity modes in (2.7) are recovered. For a compact source region, the final result for
$S_i(x_{\!j}, t)$
in (A2) may include limiter functions
$\lambda (x_{\!j})$
selected to provide a smooth transition to the undisturbed region of uniform flow in
$x_1$
and
$x_{\!j}$
directions, such as
where
$x_{\textit{sj}}$
determine the extension of the source region in the
$x_{\!j}$
direction. Note that a shear layer is generated in the region where
$\lambda (x_{\!j})$
vary, but the resulting pressures have generally very moderate amplitudes. The generated solutions are well behaved as long as the parameters are selected to avoid very high values of
$K_n$
.
It should be noted that the appendix presents an idealised derivation under uniform convection to obtain (2.7) and (A1). In the numerical implementation, however, the resulting source term is embedded in the full Navier–Stokes equations and acts on the actual non-uniform mean flow. This defines a locally isotropic, analytically uniform forcing embedded in a spatially non-uniform mean flow, providing a clear separation between the idealised construction of the source and its physical interaction with the realistic inflow.
Once introduced, the synthesised eddies are convected by the non-uniform mean velocity profile, such that the initially isotropic fluctuations evolve under local shear into a height-dependent turbulence field consistent with the boundary-layer-type variation observed in the LES. The non-uniformity upstream is therefore handled by the solver, while the uniform-flow assumption serves only as a device to obtain a closed-form expression for the forcing. In future extensions, the same forcing could be formulated in a piecewise (layered) manner to better capture strongly sheared or wall-bounded inflows, but that was beyond the scope of the present study.
Appendix B
To assess the grid independent study and domain sensitivity studies, separate simulations were performed. Regarding the grid independent study, two supplementary simulations were performed for the disturbed inflow case with
$10\,\%$
TI, using 50 and 60 cells across the cylinder width (compared with 70 in the baseline case). Figure 33 compares the time-averaged velocity components and their RMS values across these three meshes. The results show no significant differences, confirming that the baseline grid is sufficiently resolved for the present framework. Regarding the domain size sensitivity analysis, for the case with
$\text{TI} = 10\,\%$
, the inlet–cylinder distance was increased from
$8d$
to
$10d$
, and the cylinder–outlet distance was extended from
$20d$
to
$25d$
, while all other numerical settings were kept unchanged. According to figure 34, comparison of the time-averaged velocity components and their RMS values between the baseline and extended domain cases revealed no significant differences.
Time-averaged distributions of the non-dimensional velocity components and their RMS values along seven vertical lines located at
$x/h = 3.875$
,
$4.25$
,
$4.5$
,
$4.625$
,
$4.875$
,
$5.25$
and
$5.875$
on the
$z$
-normal plane at
$z/h = 0$
. The results correspond to the case with
$\text{TI} = 10\,\%$
and compare the baseline case with 70 cells against two additional cases with 50 and 60 cells.

To analyse the sensitivity of the SPOD results to the number of blocks, three additional spectral analyses were performed with 22, 25 and 31 blocks for the undisturbed case. A comparison of the first SPOD mode from the four post-processing configurations is illustrated in figure 35. The results show that the dominant frequencies and modal energy distributions remain unchanged as the number of blocks varies, indicating that the SPOD results are statistically converged and insensitive to the chosen block number.
Time-averaged distributions of the non-dimensional velocity components and their RMS values along seven vertical lines located at
$x/h = 3.875$
,
$4.25$
,
$4.5$
,
$4.625$
,
$4.875$
,
$5.25$
and
$5.875$
on the
$z$
-normal plane at
$z/h = 0$
. The results correspond to the case with
$\text{TI} = 10\,\%$
and compare the baseline case with 70 cells and the original domain size against additional cases with an extended domain.

Comparison of first SPOD mode in four conducted SPOD post-processing at the
$z_1$
plane and Domain 1.

Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.084$
at the
$z_3$
plane.

Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.101$
. From left to right,
$z_1$
,
$y_4$
and
$x_5$
planes.

Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.084$
in the disturbed case with
$\text{TI}=20\,\%$
. From left to right,
$z_1$
,
$y_4$
and
$x_5$
planes.

Spatial distribution of mode 1 for three components of velocity at
${\textit{St}}=0.084$
in the disturbed case with
$\text{TI}=30\,\%$
. From left to right,
$z_1$
,
$y_4$
and
$x_5$
planes.




















































































































































































