1. Introduction
Flows past bluff bodies and their transitions with increasing Reynolds numbers from steady two-dimensional (2-D) wake flows, through unsteady and three-dimensional (3-D) flows, to fully turbulent wakes have attracted much attention over the years (Oertel Reference Oertel1990; Williamson Reference Williamson1996b ; Choi et al. Reference Choi, Jeon and Kim2008; Thompson et al. Reference Thompson, Leweke and Hourigan2021), as their relevance goes beyond the fundamental interest and encompasses several industrial applications, for example, in the field of vortex-induced vibrations (Williamson & Govardhan Reference Williamson and Govardhan2008).
1.1. Two-dimensional cylinders
Most studies on 2-D bluff bodies have focused on circular and square cylinders as prototypes to characterise the flow bifurcations. At the critical Reynolds number (based on the free-stream velocity
$U_\infty$
and cylinder diameter
$D$
)
$Re_c \approx 45{-}47$
the flow undergoes a Hopf bifurcation from a symmetric steady state towards a time-periodic asymmetric state (Noack & Eckelmann Reference Noack and Eckelmann1994) that gives origin to the von Kármán vortex shedding. The triggering mechanism is known to result from a global instability (Jackson Reference Jackson1987), which arises when the region of local absolute instability is large enough (Chomaz Reference Chomaz2005), and whose onset therefore depends on the size of the wake recirculation region and on the maximum reverse flow velocity (Chiarini et al. Reference Chiarini, Quadrio and Auteri2022c
). At larger Reynolds numbers the flow undergoes a secondary instability and becomes three dimensional. For the circular cylinder, mode A with spanwise wavelength
$\lambda \approx 3.9 D$
becomes unstable at
$Re \approx 190$
, while mode B with
$\lambda \approx 1.2 D$
becomes unstable at slightly larger
$Re$
(Barkley & Henderson Reference Barkley and Henderson1996; Williamson Reference Williamson1996a
,Reference Williamson
b
). A further quasi-periodic mode with
$\lambda \approx 2.5 D$
has been found to become unstable at larger Reynolds numbers (Blackburn & Lopez Reference Blackburn and Lopez2003; Blackburn et al. Reference Blackburn, Marques and Lopez2005; Blackburn & Sheard Reference Blackburn and Sheard2010).
For 2-D rectangular cylinders, the flow dynamics and the physical mechanism of the bifurcations vary with the length-to-height ratio
$L/H$
. For
$Re \geqslant 300$
, in fact, the Strouhal number of the flow increases in an almost stepwise manner with
$L/H$
, due to a pressure feedback loop that locks the shedding of vortices from the shear layers separating from the leading edge (LE) and trailing edge (TE) (Okajima Reference Okajima1982; Nakamura et al. Reference Nakamura, Ohya and Tsuruta1991; Mills et al. Reference Mills, Sheridan, Hourigan and Welsh1995; Hourigan et al. Reference Hourigan, Thompson and Tan2001). This stepwise increase is closely related to the number of large-scale vortices that fit along the sides of the rectangular cylinders. Two different regimes are possible depending on the relative phase of the LE and TE vortices (Chiarini et al. Reference Chiarini, Quadrio and Auteri2022d
). Using dynamic mode decomposition, Zhang et al. (Reference Zhang, Kareem, Xu and Jiang2023) observed that the feedback mechanism at
$Re=1000$
changes with the length-to-height ratio. For small
$L/H$
, the pressure feedback loop encompasses the whole separation region and the flow is controlled by the impinging shear-layer instability; for larger
$L/H$
, the feedback loop covers the entire chord of the cylinder and the flow is controlled by the LE vortex shedding instability. The first onset of three-dimensionality also depends on
$L/H$
: unlike short cylinders, where the flow becomes three dimensional via modes A, B and C of the wake (Robichaux et al. Reference Robichaux, Balachandar and Vanka1999; Blackburn & Lopez Reference Blackburn and Lopez2003), for long enough rectangular cylinders, Chiarini et al. (Reference Chiarini, Quadrio and Auteri2022b
) found that the first 3-D unstable mode is of quasi-subharmonic nature (mode QS) and is due to an inviscid mechanism triggered by the interaction of the vortices placed over the cylinder sides.
1.2. Axisymmetric 3-D bluff bodies
Less studies have considered flows past 3-D bluff bodies, although they are ubiquitous in human life and engineering applications, e.g. tall buildings, chimneys, pylons, cars, trains and other ground vehicles. These flows exhibit a richer physics and, already at low
$Re$
, the first bifurcations and their sequence depend on the body geometry.
Unlike 2-D cylinder wakes, the flow past a sphere becomes asymmetric prior to a transition to unsteady flow (Magarvey & Bishop Reference Magarvey and Bishop1961a
,Reference Magarvey and Bishop
b
; Magarvey & MacLatchy Reference Magarvey and MacLatchy1965). The wake remains steady and axisymmetric up to
$Re \approx 211$
(Johnson & Patel Reference Johnson and Patel1999), and then transitions to a steady asymmetric state through the regular bifurcation of an eigenmode of azimuthal wavenumber
$m=1$
(Tomboulides & Orszag Reference Tomboulides and Orszag2000). The resulting wake is characterised by a pair of steady streamwise vortices and exhibits a reflectional symmetry about a longitudinal plane of arbitrary azimuthal orientation. At higher
$Re$
, a Hopf bifurcation renders the flow oscillatory. Natarajan & Acrivos (Reference Natarajan and Acrivos1993) found with stability analysis that an unsteady
$m=1$
mode of the steady axisymmetric base flow becomes unstable at
$Re \approx 277$
, as confirmed by the experiments and simulations of Tomboulides et al. (Reference Tomboulides, Orszag and Karniadakis1993); Johnson & Patel (Reference Johnson and Patel1999); Tomboulides & Orszag (Reference Tomboulides and Orszag2000). This unsteady regime consists of hairpin vortices (HVs) shed downstream of the sphere in the same plane as that of the steady asymmetric regime. Citro et al. (Reference Citro, Siconolfi, Fabre, Giannetti and Luchini2017) performed a 3-D global stability analysis of the steady asymmetric flow to characterise the eigenmode responsible for this second bifurcation. They found a critical Reynolds number of
$Re = 272$
, and deduced from a structural sensitivity analysis (Giannetti et al. Reference Giannetti, Camarri and Luchini2010) that the instability is driven by the region outside the asymmetric wake. At larger Reynolds numbers,
$Re\gt 600$
, the wake loses its periodicity and the flow becomes more chaotic (Magarvey & Bishop Reference Magarvey and Bishop1961
Reference Magarvey and Bishop
b
).
A similar bifurcation scenario has been observed for the flow past a disk perpendicular to the incoming flow. At
$Re \approx 115$
, a regular bifurcation leads to a 3-D steady asymmetric state, with a reflectional symmetry Natarajan & Acrivos (Reference Natarajan and Acrivos1993); Fabre et al. (Reference Fabre, Auguste and Magnaudet2008); Meliga et al. (Reference Meliga, Chomaz and Sipp2009). At slightly larger
$Re$
, i.e.
$Re \approx 121$
, a Hopf bifurcation breaks the remaining reflectional symmetry and time invariance. At
$Re \approx 140$
, a third bifurcation occurs that preserves the flow unsteadiness but restores a reflectional symmetry orthogonal to that preserved by the first bifurcation.
Similarly, Bohorquez et al. (Reference Bohorquez, Sanmiguel-Rojas, Sevilla, Jiménez-González and Martínez-Bazán2011) investigated the stability of the flow past bullet-shaped bodies (slender cylindrical bodies with an elliptical LE and a blunt TE). They found the same sequence of bifurcations and observed that the critical
$Re$
increases with the aspect ratio. Compared with the sphere and disk wakes, the symmetry plane selected by the first regular bifurcation is preserved over a larger range of Reynolds numbers, up to
$Re \approx 500$
.
1.3. Non-axisymmetric 3-D bluff bodies
Many studies on non-axisymmetric 3-D bluff bodies have focused on two simple types of geometry: finite circular cylinders and rectangular prisms.
Inoue & Sakuragi (Reference Inoue and Sakuragi2008) investigated the vortex shedding past finite circular cylinders of span-to-diameter ratio
$0.5 \leqslant W/D \leqslant 100$
with direct numerical simulations (DNS) for
$Re \leqslant 300$
. They found that the flow changes drastically depending on
$W/D$
and
$Re$
, and identified five basic wake patterns: (i) periodic oblique vortex shedding, (ii) quasi-periodic oblique vortex shedding, (iii) periodic hairpin vortex falloff, (iv) two stable counter-rotating vortex pairs, and (v) alternating shedding of counter-rotating vortex pairs from the free ends. For the same geometry with
$W/D \in [0.5,2]$
, Yang et al. (Reference Yang, Feng and Zhang2022) obtained consistent results by DNS and linear stability analysis (LSA) for
$Re \leqslant 1000$
.
Rectangular prisms differ from finite circular cylinders, as the sharp corners set the location of flow separation. The bifurcation sequence and critical
$Re$
are expected to change with the length-to-height ratio
$L/H$
(as already observed for 2-D rectangular cylinders) and width-to-height ratio
$W/H$
. Marquet & Larsson (Reference Marquet and Larsson2015) investigated the linear stability of the steady flow past thin rectangular plates with
$L/H=1/6$
and various widths. For
$ 1 \leqslant W/H \leqslant 2$
, a pitchfork bifurcation breaks the top/bottom planar symmetry, similar to the flow past axisymmetric bodies like spheres and disks. For
$W/H \gt 2.5$
, a Hopf bifurcation breaks the top/bottom symmetry and leads to a time-periodic regime, similar to 2-D cylinders. Surprisingly, for intermediate widths
$ 2 \leqslant W/H \leqslant 2.5$
, a Hopf bifurcation breaks the left/right planar symmetry. Zampogna & Boujo (Reference Zampogna and Boujo2023) set the width to
$W/H=1.2$
and investigated the flow bifurcations for
$1/6 \leqslant L/H \leqslant 3$
. Their analysis yields a series of pitchfork and Hopf bifurcations. Two stationary modes that break respectively the top/bottom and left/right planar symmetries become first unstable. At larger
$Re$
, two oscillatory modes become unstable and again break either the top/bottom or left/right symmetry. The critical
$Re$
of the first bifurcation increases monotonically with
$L/H$
, similar to flows past 2-D rectangular cylinders (Chiarini et al. Reference Chiarini, Quadrio and Auteri2021) and axisymmetric bodies (disk, sphere and bullet-shaped bodies). For the flow past a cube, this bifurcation scenario is confirmed by the low-
$Re$
experiments of Klotz et al. (Reference Klotz, Goujon-Durand, Rokicki and Wesfreid2014) and DNS of Saha (Reference Saha2004) and Meng et al. (Reference Meng, An, Cheng and Kimiaei2021): at
$Re \approx 216$
the steady symmetric flow bifurcates towards a steady asymmetric regime with only one planar symmetry; at
$Re=265$
the flow undergoes a Hopf bifurcation and oscillates about the bifurcated asymmetric regime, with a shedding of HVs similar to that past a sphere and a disk.
1.4. Present study
In this work we take a step forward from the works of Marquet & Larsson (Reference Marquet and Larsson2015) and Zampogna & Boujo (Reference Zampogna and Boujo2023), and characterise the sequence of bifurcations for the flow past rectangular prisms in a 3-D space of parameters, varying simultaneously the length-to-height and width-to-height ratios and the Reynolds number. We use linear and weakly nonlinear (WNL) stability analyses to characterise the first bifurcations, and fully nonlinear DNS to describe the flow regimes at larger
$Re$
. We consider
$1 \leqslant L/H \leqslant 5$
,
$1.2 \leqslant W/H \leqslant 5$
and
$Re \lessapprox 700$
. The DNS focus on
$L/H=5$
, similar to the benchmark known as BARC (benchmark on the aerodynamics of a rectangular 5:1 cylinder; see https://www.aniv-iawe.org/barc-docs), which aims to characterise the flow past elongated bluff bodies and set standards for simulations and experiments in the turbulent regime.
The structure of the paper is as follows. Section 2 briefly describes the mathematical formulation and numerical methods. The low-
$Re$
steady flow is characterised in § 3.1, while §§ 3.2and 3.3 deal with the LSA and § 3.4 with the WNL stability analysis. Section 4 presents the DNS results and details the flow regimes after the first bifurcations for
$L/H=5$
. A concluding discussion is provided in § 5.
2. Methods
2.1. Flow configuration

Figure 1. Overview of the prism geometry for various lengths
$L$
and widths
$W$
.
The incompressible flow past rectangular prisms with aspect ratio
$1 \leqslant L/H \leqslant 5$
and
$1.2 \leqslant W/H \leqslant 5$
is considered; here
$L$
,
$W$
and
$H$
are the length, width and height of the prism, i.e. the sizes in the streamwise, spanwise and vertical directions; see figure 1. A Cartesian reference frame is placed at the LE of the prisms, with the
$x$
axis aligned with the flow direction and the
$y$
and
$z$
axes denoting the spanwise and vertical directions. The Reynolds number is based on
$H$
and on the free-stream velocity
$U_\infty$
. Unless otherwise stated, all quantities are made dimensionless with
$H$
and
$U_\infty$
. The flow is governed by the incompressible Navier--Stokes (NS) equations for the velocity
$\boldsymbol {u} = (u,v,w)$
and pressure
$p$
,

2.2. Mathematical formulation
The onset of instability is studied by LSA. The field
$\{\boldsymbol {u},p\}$
is divided into a time-independent base flow
$\{\boldsymbol {u}_0,p_0\}$
and an unsteady contribution
$\{\boldsymbol {u}_1,p_1\}$
of small amplitude
$\epsilon$
, i.e.

Using this decomposition in the NS equations (2.1), the steady base flow equations for
$\{\boldsymbol {u}_0,p_0\}$
are obtained at order
$\epsilon ^0$
, while the linearised NS equations governing the small perturbations
$\{\boldsymbol {u}_1,p_1\}$
are obtained at order
$\epsilon ^1$
. Using the normal-mode ansatz
$\{\boldsymbol {u}_1,p_1\}(\boldsymbol {x},t)=\{\hat {\boldsymbol {u}}_1,\hat {p}_1\}(\boldsymbol {x})e^{\lambda t} + c.c.$
(where
$c.c$
stands for complex conjugate) yields an eigenvalue problem for the complex eigenvalue
$\lambda = \sigma + i \omega$
and eigenvector
$\{\hat {\boldsymbol {u}}_1,\hat {p}_1\}$
:

The linear stability of the system is determined by the sign of the real part of
$\lambda$
, i.e
$\sigma$
. If all
$\sigma \lt 0$
, perturbations decay and the flow is stable. If
$\sigma \gt 0$
for at least one mode, perturbations grow exponentially. A mode is stationary or oscillatory if
$\omega =0$
or
$\omega \gt 0$
, respectively. For each geometry
$(L,W)$
, several modes may become unstable as
$Re$
increases; the critical Reynolds number
$Re_c$
of each mode is computed by cubic interpolation of
$\sigma (Re)$
with at least three values of
$Re$
that bracket
$\sigma = 0$
.
2.3. Numerical method
Two different numerical methods are used. The LSA is performed with the non-commercial, finite elements software FreeFem++ (Hecht Reference Hecht2012), while the DNS are performed with an in-house code based on finite differences.
For the LSA (§ 3.2), the steady base flow solution of (2.1) and the eigenvalue and eigenmodes solution of (2.3) are calculated with the same method as in Zampogna & Boujo (Reference Zampogna and Boujo2023). Given the two symmetry planes
$y=0$
and
$z=0$
of the prism, we consider only the quarter space
$y,z\geqslant 0$
and build in the numerical domain
$\{x,y,z \, | \, -10\leqslant x \leqslant 20; 0\leqslant y,z \leqslant 10\}$
an unstructured tetrahedral mesh with nodes strongly clustered near the prism. See Appendix A for convergence studies with respect to mesh and domain sizes. The weak form of the equations is discretised with Arnold--Brezzi--Fortin MINI-elements (P1 for pressure and P1b for each velocity component). For
$W/H \gt 2.5$
, dimensions in the
$y$
direction are normalised by
$W$
and
$y$
derivatives are scaled by a factor
$1/W$
. The nonlinear base flow equations are solved with a Newton method, and the eigenvalue problem with a Krylov--Schur method. Calculations are run in parallel with a domain-decomposition method, typically on 48 processes. Boundary conditions are as follows: free-stream velocity
${\boldsymbol {u}}_0=(1,0,0)^T$
and zero perturbation
${\boldsymbol {u}}={\boldsymbol {0}}$
at the inlet; then, for both the base flow and the perturbations, a no-slip condition
${\boldsymbol {u}}={\boldsymbol {0}}$
on the body, stress-free condition
$-p \boldsymbol {n} + Re^{-1} \boldsymbol {\nabla } \boldsymbol {u} \cdot \boldsymbol {n}=\boldsymbol {0}$
at the outlet and symmetric (
$S$
) conditions
$u_{n}=0, \partial _n \boldsymbol {u}_{t}=\boldsymbol {0}$
or antisymmetric (
$A$
) conditions
$\boldsymbol {u}_{t}=\boldsymbol {0}, \partial _n u_{n}=0$
on the symmetry planes
$y=0$
and
$z=0$
(where subscripts
$n$
and
$t$
denote normal and tangential directions, respectively), the specific choice depending on the considered field, i.e. doubly symmetric base flow
$S_y S_z$
or eigenmode belonging to one of the four possible families
$S_y S_z$
,
$S_y A_z$
,
$A_y S_z$
and
$A_y A_z$
. Adjoint modes, used to compute structural sensitivities, are obtained in a similar way (see Zampogna & Boujo (Reference Zampogna and Boujo2023) for more details).
The 3-D nonlinear simulations (§ 4) consider
$L/H=5$
and
$W/H=1.2,2.25$
and
$5$
. The NS equations are solved using a DNS code introduced by Luchini (Reference Luchini2016) and written in the computer programming language CPL (Luchini Reference Luchini2021), which employs second-order finite differences on a staggered grid in the three directions. This DNS code has been previously validated and used to compute the flow past rectangular cylinders in the laminar (Chiarini et al. Reference Chiarini, Quadrio and Auteri2022b
) and turbulent (Chiarini & Quadrio 2021, 2022; Chiarini et al. Reference Chiarini, Gatti, Cimarelli and Quadrio2022a
) regimes. The momentum equation is advanced in time by a fractional step method using a third-order Runge--Kutta scheme. The Poisson equation for the pressure is solved using an iterative SOR (successive over-relaxation) algorithm. The presence of the prism is dealt with by a second-order implicit immersed boundary method, implemented in staggered variables (Luchini Reference Luchini2013). The computational domain extends for
$-25H \leqslant x \leqslant 75H$
,
$-40 H \leqslant y \leqslant 40 H$
and
$-25H \leqslant z \leqslant 25H$
in the streamwise, spanwise and vertical directions, corresponding to domain sizes
$L_x=100H$
,
$L_y=80H$
and
$L_z=50H$
. The computational domain is discretised with
$N_x=1072$
and
$N_z=590$
points in the streamwise and vertical directions. The number of points in the spanwise direction is
$N_y=666$
,
$720$
and
$804$
for
$W/H=1.2, 2.25$
and
$5$
, respectively. The points are distributed with a geometrical progression in the three directions, to obtain a higher resolution close to the prisms, especially close to the corners and in the wake. Close to the corners the grid spacing is
$\Delta x = \Delta y = \Delta z \approx 0.008H$
. The adaptive time step is computed at every iteration to enforce a Courant--Friedrichs--Lewy (CFL) number below unity, yielding an average value of
$\Delta t \approx 0.0066 H/U_\infty$
. After the transient regime, all simulations are advanced for more than
$1000$
convective time units
$H/U_\infty$
to ensure convergence of all quantities of interest. See Appendix A for a convergence study with respect to the grid resolution and the time step.
3. Global LSA
3.1. Base flow

Figure 2. Base flow near the first bifurcation, visualised with isosurfaces of zero streamwise velocity (
$u_0 = 0$
) coloured by streamwise vorticity
$\omega _x$
.

Figure 3. Base flow properties: (a) recirculation length, (b) drag coefficient. Solid lines:
$l_r(Re)$
and
$C_x(Re)$
for different body lengths
$L=1/6, 0.5, 1 \ldots , 4.5, 5$
; dotted line:
$l_r$
and
$C_x$
at the onset of the first bifurcation according to the LSA (§ 3.2). Body width from left to right:
$W=1.2, 2.25, 5$
.
The steady symmetric base flow is shown in figure 2 by means of isosurfaces of zero streamwise velocity (
$u_0=0$
) coloured by streamwise vorticity
$\omega _x$
. The wake has a different shape depending on the prism geometry. Long and narrow bodies (e.g.
$L=5$
,
$W=1.2$
) tend to exhibit a convex backflow region centred about the vertical and horizontal symmetry planes
$y=0$
and
$z=0$
, somewhat reminiscent of axisymmetric wakes. For these bodies, the top/bottom and left/right recirculation regions originating from the LEs reattach on the body upstream of the TEs. By contrast, short and wide bodies (e.g.
$L=1$
,
$W=5$
) tend to exhibit a non-convex backflow region, with two lobes above and below the horizontal symmetry plane
$z = 0$
. Interestingly, moving downstream the backflow region enlarges in the vertical direction
$z$
and shrinks in the spanwise direction
$y$
, resulting in ‘peanut’-shaped
$y-z$
cross-sections. In addition, the upper/lower recirculation regions originating from the LEs extend all the way down to the TEs without reattaching on the body and are therefore connected with the main wake recirculation.
Some quantitative properties of the base flow are shown in figure 3. For all the body widths considered, the recirculation length
$l_r$
(length of the downstream recirculation region along the symmetry axis
$y=z=0$
measured from the TE, figure 3
$a$
) generally increases with
$Re$
, like most 2-D and 3-D bluff body wakes (Giannetti & Luchini Reference Giannetti and Luchini2007; Marquet & Larsson Reference Marquet and Larsson2015). Also,
$l_r$
decreases with
$L$
due to a decreasing LE separation angle, in agreement with the observations of Zampogna & Boujo (Reference Zampogna and Boujo2023) for prisms with
$W = 1.2$
and of Chiarini et al. (Reference Chiarini, Quadrio and Auteri2022c
) for 2-D bluff bodies of different shapes. The recirculation length is also seen to increase with
$W$
. Anticipating on the linear stability results (§ 3.2), we follow the evolution of the recirculation length at the first bifurcation (lowest critical
$Re$
; blue dashed line):
$l_r$
increases with
$W$
, while it remains fairly independent on
$L$
with an overall slightly decreasing trend. For
$W=2.25$
, we note that the value of
$l_r$
at the onset of the primary bifurcation is not continuous with
$L$
, because the mode that drives the primary instability changes with the length of the prism and the lowest
$Re_c$
jumps at
$L\simeq 1.5$
(see figure 5).
The base flow drag coefficient
$C_x =2F_x/(\rho U_\infty H W)$
, where
$F_x$
is the drag force, decreases or stays constant with the Reynolds number in almost the complete space of parameters; see figure 3
$(b)$
. For narrower bodies,
$W=1.2$
and
$2.25$
, the drag coefficient has a non-monotonic trend with
$L$
, first decreasing and then increasing. For longer bodies,
$W = 5$
, it monotonically increases with
$L$
for all
$Re$
. At the first bifurcation (blue dashed line),
$C_x$
varies with
$L$
following different trends depending on the body width: it decreases and then increases for
$W = 1.2$
, monotonically decreases for
$W = 2.25$
and remains essentially constant for
$W = 5$
.

Figure 4. Stability diagram for the first bifurcation. Filled symbols: pitchfork (stationary) bifurcation; open symbols: Hopf (oscillatory) bifurcation. Red: symmetry breaking in the horizontal
$y$
direction (
$A_yS_z$
); black: symmetry breaking in the vertical
$z$
direction (
$S_yA_z$
). Circles: actual calculations; squares: interpolations along
$L$
or
$W$
. By symmetry, for
$W=1$
, the
$S_y A_z$
and
$A_y S_z$
eigenmodes bifurcate simultaneously.
3.2. Linear stability
We now turn to the LSA. For each geometry, we consider the steady symmetric base flow and look at the eigenmodes that become linearly unstable as
$Re$
increases. In particular, we are interested in the symmetries of the modes and their stationary/oscillatory nature (pitchfork/Hopf bifurcation).
3.2.1. First bifurcation: stability diagram
We start considering the first bifurcation. Figure 4 summarises the effect of the geometry. Three main regions appear in the
$L$
-
$W$
plane. The flow past wide bodies (large
$W/L$
) first become unstable to oscillatory
$S_y A_z$
perturbations that break the temporal and top/bottom planar symmetries and preserve the left/right planar symmetry, leading to a periodic vortex shedding across the shorter body dimension. This is consistent with the 2-D limit
$W\rightarrow \infty$
, and resembles the results for the flow past finite-length circular cylinders (Yang et al. Reference Yang, Feng and Zhang2022). Increasing
$L/W$
, the flow first becomes unstable to stationary
$S_y A_z$
perturbations that break the top/bottom spatial symmetry and maintain the left/right one. This bifurcation leads to a static vertical deflection of the wake. Finally, further increasing
$L/W$
, the flow first becomes unstable to stationary
$A_y S_z$
perturbations that break the left/right symmetry while preserving the top/bottom one. This corresponds to a static horizontal deflection of the wake. By symmetry, the
$S_y A_z$
and
$A_y S_z$
modes become unstable simultaneously in the special case of bodies of square cross-section
$W/H=1$
, as observed by Marquet & Larsson (Reference Marquet and Larsson2015), Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) and Zampogna & Boujo (Reference Zampogna and Boujo2023). Interestingly, the oscillatory
$A_y S_z$
mode observed by Marquet & Larsson (Reference Marquet and Larsson2015) for thin plates
$L=1/6$
of intermediate width
$2 \leqslant W \leqslant 2.5$
seems restricted to a narrow region of the
$L$
-
$W$
plane (see also the neutral curves for
$W=2.25$
in figure 5 and for
$L=1/6$
and 1 in figure 7).
3.2.2. Eigenmodes and neutral curves
We now proceed to characterise in more detail the linear stability of the steady symmetric base flow for increasing values of
$Re$
and a few selected geometries. Figure 5 shows the critical Reynolds number and frequency of the first bifurcations as a function of the body length
$L$
for
$W=1.2$
,
$2.25$
and
$5$
. The eigenmodes of the first two bifurcations for these widths and the two body lengths
$L=1$
and
$5$
are shown in figures 6; see also the supplementary material. We focus on the first few bifurcations, as the eigenmodes that become unstable at significantly larger Reynolds numbers become less relevant.

Figure 5. (a) Neutral curves: critical Reynolds number
$Re_c$
as a function of body length
$L$
for different body widths
$W=1.2$
, 2.25 and 5. Filled symbols and solid lines indicate stationary (pitchfork) bifurcations; open symbols and dashed lines indicate oscillatory (Hopf) bifurcations. (b) Critical frequency as a function of body length
$L$
for the first few bifurcations: Strouhal number
$St_c=\omega (Re_c)/(2\pi )$
. Different lines correspond to different modes. In both
$(a)$
and
$(b)$
, thicker lines show the eigenmode that becomes unstable at the lowest Reynolds number.

Figure 6. First and second bifurcating eigenmodes (left and right, respectively), visualised with isosurfaces of streamwise velocity.
For narrow bodies,
$W=1.2$
, the first two bifurcations of the steady symmetric base flow occur in short succession (as
$W$
is close to 1), and correspond to stationary
$S_y A_z$
and
$A_y S_z$
modes that break either of the two planar symmetries (Zampogna & Boujo Reference Zampogna and Boujo2023). A pair of streamwise streaks of positive and negative streamwise perturbation velocity arise in the wake (figure 6), leading to a vertical or horizontal displacement of the wake. Other bifurcations occur at larger
$Re$
and correspond to oscillatory
$S_y A_z$
and
$A_y S_z$
modes. For small
$L$
, these modes are almost as unstable as their stationary counterparts, but they quickly become much more stable as
$L$
increases, i.e. the corresponding
$Re_c$
increases faster; accordingly, for long bodies, these oscillatory modes of the symmetric base flow are not directly relevant to the secondary stability of the deflected flow (§ 4.1.2). The associated angular frequency is of the order
$\omega \simeq 0.4-0.6$
(
$St = f H / U_\infty \simeq 0.06-0.09$
, where
$f = \omega /(2 \pi )$
) for
$L \leqslant 2$
. The increase of
$Re_c$
with
$L$
is consistent with the stationary bifurcation of the flow past axisymmetric bodies: a disk much thinner than its diameter
$L/D \ll 1$
, a sphere
$L/D=1$
, bullet-shaped bodies
$1 \leqslant L/D \leqslant 3.5$
(Natarajan & Acrivos Reference Natarajan and Acrivos1993; Fabre et al. Reference Fabre, Auguste and Magnaudet2008; Meliga et al. Reference Meliga, Chomaz and Sipp2009; Bohorquez et al. Reference Bohorquez, Sanmiguel-Rojas, Sevilla, Jiménez-González and Martínez-Bazán2011) and with the primary oscillatory bifurcation of the flow past 2-D bodies of different shape (Jackson Reference Jackson1987; Chiarini et al. Reference Chiarini, Quadrio and Auteri2022c
).
For wide bodies,
$W=5$
, the first bifurcation corresponds to an oscillatory
$S_y A_z$
mode for all
$L$
(as already observed in figure 4). The critical Reynolds number increases with
$L$
and is rather well separated from the following bifurcations. The angular frequency is
$\omega \simeq 0.5-0.6$
(
$St \simeq 0.08-0.10$
). The unstable eigenmode is stronger in the central region (small
$|y|$
); see figure 6. The second bifurcation corresponds to an oscillatory
$S_y S_z$
modes for all
$L$
. Like the leading
$S_y A_z$
mode, it is associated with vortex shedding in the upper and lower shear layers, with a similar frequency but different inter-layer phasing. Unlike the
$S_y A_z$
mode that breaks the top/bottom symmetry, the
$S_y S_z$
mode preserves both planar symmetries. In other words, in the vertical plane
$y=0$
the two leading modes can be categorised as ‘sinuous’ and ‘varicose’, respectively. To the best of our knowledge, unsteady doubly symmetric
$S_y S_z$
modes have not been reported in 3-D bluff body wakes in previous works except for thin plates (
$L=1/6$
,
$W\gt 3.1$
, Marquet & Larsson Reference Marquet and Larsson2015); see also § 4.2. At larger
$Re$
further bifurcations are observed, including some involving doubly antisymmetric modes (
$A_y A_z$
). All these bifurcations are oscillatory and the corresponding frequencies are almost independent of
$L$
, and differ by one order of magnitude depending on the bifurcation; see figure 5.
For bodies of intermediate width,
$W=2.25$
, the increasing trend of
$Re_c$
with
$L$
is similar, but the leading mode changes with
$L$
: stationary for
$L \gt 1.5$
(similar to narrower bodies) and oscillatory with
$\omega \simeq 0.3-0.6$
(
$St \simeq 0.05-0.10$
) for
$L \lt 1.5$
(similar to wider bodies). As mentioned previously (figure 4), for this specific width, the first bifurcation breaks the left/right symmetry for very short bodies (
$L\lt 0.5$
), which leads to the unusual vortex shedding across the larger body dimension observed by Marquet & Larsson (Reference Marquet and Larsson2015), and breaks the top/bottom symmetry otherwise, leading to the more usual vortex shedding across the smaller dimension (
$0.5 \lt L \lt 1.75$
) or to a static vertical wake deflection (
$L\gt 1.75$
). Interestingly, the stationary
$S_y A_z$
mode remains fully stable in the range of
$Re$
investigated when
$L\lt 1.5$
, while for
$L\gt 1.5$
, it is restabilised for large enough
$Re$
. Other bifurcations at larger
$Re$
draw a more complicated pattern than for
$W=1.2$
and
$5$
. Their order of appearance depends on
$L$
and some modes are unstable only in some interval of
$L$
. For almost all values of
$L$
, however, doubly symmetric
$S_y S_z$
and doubly antisymmetric
$A_y A_z$
modes only correspond to the third or higher bifurcations.
Figure 7 offers a complementary picture, with neutral curves shown as a function of the body width
$W$
for the three body lengths
$L=1/6$
,
$L=1$
and
$L=5$
. The first critical
$Re$
tends to decrease with
$W$
, albeit not monotonically: the first bifurcation lies on the envelope of different neutral curves corresponding to different modes, each mode being more unstable in a preferred interval of
$W$
.

Figure 7. Same as figure 5, as a function of body width
$W$
for different body lengths
$L=1/6$
(data from Marquet & Larsson Reference Marquet and Larsson2015),
$1$
and
$5$
.
3.3. Structural sensitivity
We now look at the structural sensitivity of the two leading modes. Bodies with
$L=1$
and
$5$
and
$W=1.2$
,
$2.25$
and
$5$
are considered. Figure 8 shows isosurfaces of
$||\hat {\boldsymbol {u}}_1||\,||\hat {\boldsymbol {u}}_1^{\dagger }||$
, the product of the Euclidean norms of the eigenmode (direct mode) and the associated adjoint mode. This quantity was introduced by Giannetti & Luchini (Reference Giannetti and Luchini2007) as an upper bound for the eigenvalue variation
$|\delta \lambda |$
induced by a specific perturbation of the linearised NS operator, namely a ‘force--velocity coupling’ representing feedback from a localised velocity sensor to a localised force actuator at the same location. In this sense, the structural sensitivity is an indicator of the eigenvalue sensitivity and identifies the wavemaker (Monkewitz et al. Reference Monkewitz, Huerre and Chomaz1993). The most sensitive regions, where the direct and adjoint modes overlap, are found in the near wake of the body for all modes, which is typical of bluff body wakes. As expected, these regions with large sensitivity are symmetric with respect to both
$y$
and
$z$
.
All the
$S_y$
modes, which preserve the left/right symmetry, are most sensitive in two regions located symmetrically with respect to the horizontal symmetry plane
$z=0$
. For oscillatory
$S_y$
modes, observed for smaller
$L/W$
, the structural sensitivity is maximum exactly on the
$u_0=0$
isosurface, in two flat and elongated regions separated by more than one body height; see also the cross-sections in the supplementary material. By contrast, for stationary
$S_y$
modes observed for larger
$L/W$
, the structural sensitivity is maximum inside the backflow region, in more compact regions located closer to the horizontal plane.
Stationary
$A_yS_z$
modes are most sensitive in two regions located symmetrically with respect to the vertical symmetry plane
$y=0$
, and the structural sensitivity is maximum inside the backflow region.
The peculiar oscillatory
$A_y S_z$
mode observed for short lengths and intermediate widths (e.g.
$L=1$
,
$W=2.25$
) has a more convoluted structural sensitivity, maximum outside and near the downstream end of the backflow region.
Overall, the structural sensitivity shows that the wavemaker of all the considered modes is located downstream, in relation with the wake recirculation region. For elongated bodies with
$L=5$
(§§ 3.4
--4), the structural sensitivity is null in the recirculating regions along the lateral sides of the prism (see the
$u_0=0$
isosurface in figure 8), meaning that these modes do not originate from the LE shear layer.

Figure 8. Structural sensitivity of the first (left) and second (right) bifurcating eigenmodes: representative isosurfaces of
$||\hat {\boldsymbol {u}}_1||\,||\hat {\boldsymbol {u}}_1^{\dagger }||$
(opaque orange). The isosurface
$u_0=0$
is reproduced from figure 2 (translucent grey).
3.4. Weakly nonlinear analysis
In § 3.2 we have investigated the linear stability of the
$S_yS_z$
steady base flow, characterising the different eigenmodes that become unstable. It is worth emphasising that, when
$Re$
increases past the first bifurcation, the linear stability of the base flow is not necessarily relevant. Rigorously speaking, to detect secondary bifurcations one should study the linear stability of the nonlinearly bifurcated state itself (e.g. static deflected wake or periodic vortex shedding) or use nonlinear simulations. For instance, the onset of the 3-D instability in the 2-D cylinder wake can be predicted using Floquet analysis of the 2-D limit cycle for
$Re\gt 47$
(Barkley & Henderson Reference Barkley and Henderson1996). In general, when the base flow is unstable to at least two eigenmodes, the nonlinear flow cannot be predicted only by comparing the growth rates or the critical
$Re$
of these modes. In some cases, WNL analysis is able to capture the dominant nonlinear interactions between competing modes, thus providing a very useful reduced-order model (a set of scalar amplitude equations) able to inform about the types of bifurcations (e.g. subcritical or supercritical), the different bifurcated flows and their linear stability, ranges of multistability and hysteresis, etc. For example, Zampogna & Boujo (Reference Zampogna and Boujo2023) found that for the flow past an Ahmed body (3-D rectangular prism with
$(L,W)=(3,1.2)$
and rounded LEs), where a steady
$S_y A_z$
mode becomes linearly unstable before a steady
$A_y S_z$
mode, the WNL analysis correctly predicts that the vertically and horizontally deflected wakes exchange their stability as
$Re$
increases.
In this section we perform a WNL analysis to study the flow behaviour in the vicinity of codimension-two points, where two eigenmodes become unstable simultaneously. Anticipating on the DNS (§ 4), we focus on
$(L,W)=(5,1.2)$
and
$(5,2.25)$
. For each of these two geometries, the top/bottom (
$S_y A_z$
) and left/right (
$A_y S_z$
) symmetry-breaking eigenmodes undergo a stationary bifurcation at close
$Re$
values. The method is the same as in Zampogna & Boujo (Reference Zampogna and Boujo2023) (see Appendix B for details about the derivation). In short, we use the technique of multiple scales. We introduce a small parameter
$\epsilon ^2$
quantifying the departure from criticality
$Re_c^{-1}-Re^{-1}$
, expand the flow field as a power series in
$\epsilon$
, introduce slow time scales
$\epsilon ^2 t$
,
$\epsilon ^4 t \ldots$
, and use a small-amplitude shift operator meant to bring together the two bifurcations at the same critical Reynolds number
$Re_c$
. By injecting all this in the NS equations and collecting like-order terms in
$\epsilon$
, we obtain a series of linear problems to be solved successively. A set of two coupled equations for the slowly varying amplitudes of the
$S_yA_z$
and
$A_yS_z$
modes, i.e.
$A$
and
$B$
, respectively, is then obtained by imposing non-resonance conditions (Fredholm alternative). Remarkably, although their coefficients are computed at
$Re=Re_c$
only, these ordinary differential equations can be solved very easily and yield stable and unstable solutions
$(A,B)$
as a function of
$Re$
(assumed sufficiently close to
$Re_c$
).
For
$(L,W)=(5,1.2)$
, we obtain the third-order amplitude equations


See Appendix B for the expressions of the
$\lambda , \chi$
and
$\eta$
coefficients. The above system has four possible equilibrium solutions (
$\mathrm {d}_t A=\mathrm {d}_t B=0$
): (i) symmetric base flow
$(A,B)=(0,0)$
described in § 3.1, (ii) pure vertical deflection
$(A,0)$
, (iii) pure horizontal deflection
$(0,B)$
, and (iv) mixed state
$(A,B)$
. Note that (3.1) and (3.2) is invariant under reflections
$A\rightarrow -A$
and
$B \rightarrow -B$
, so we discuss only
$A \geqslant 0$
and
$B\geqslant 0$
. The linear stability of each state is determined by computing the eigenvalues of the
$2\times 2$
Jacobian of (3.1) and (3.2) linearised about the state of interest. Here, we obtain the bifurcation diagram of figure 9
$(a)$
. The symmetric state
$(A,B)=(0,0)$
is stable up to
$Re \approx 352$
, when the
$S_yA_z$
eigenmode becomes unstable and the flow settles to a vertically deflected state
$(A,0)$
. Almost at the same
$Re$
, the
$A_y S_z$
eigenmode becomes unstable, but the horizontally deflected state
$(0,B)$
is initially unstable. Near
$Re \approx 353$
, there is an exchange of stability between states
$(A,0)$
and
$(0,B)$
, i.e. the wake leaves the vertically deflected state and enters a horizontally deflected state. In the very narrow range of bistability between these two states, a mixed state
$(A,B)$
exists, but is unstable. Qualitatively, the bifurcation scenario is similar to that found by Zampogna & Boujo (Reference Zampogna and Boujo2023) for Ahmed bodies of dimensions
$(L, W)=(3, 1.2)$
and a range of close-by geometries.

Figure 9. Bifurcation diagrams obtained from WNL analysis in the vicinity of codimension-two steady bifurcations. Here
$A$
and
$B$
are the amplitudes of the top/bottom (
$S_y A_z$
) and left/right (
$A_y S_z$
) symmetry-breaking eigenmodes, respectively. Solid and dashed lines denote stable and unstable branches, respectively. Results are shown for
$(a)$
$L=5, W=1.2$
and
$(b)$
$L=5, W=2.25$
.
For
$(L,W)=(5,2.25)$
, we find a subcritical pitchfork bifurcation of the
$A_y S_z$
mode. In this case, leading-order nonlinearities do not lead to saturation but instead make the system more unstable, so (3.1) and (3.2) does not yield a stable solution. This calls for the inclusion of higher-order nonlinearities to capture saturation, so we derive the fifth-order amplitude equations


The above system has the same type of equilibrium solutions as previously described. However, the pure deflected states may each have one or two different amplitudes
$(A_1,0)$
and
$(A_2,0)$
, and
$(0,B_1)$
and
$(0,B_2)$
; similarly, up to four mixed states may exist,
$(A_i,B_i)$
,
$i=1,2$
. Here, we obtain the bifurcation diagram of figure 9
$(b)$
. The symmetric state
$(A,B)=(0,0)$
is stable up to
$Re \approx 282$
, when the
$S_yA_z$
eigenmode becomes unstable and the flow settles to a vertically deflected state
$(A,0)$
. This first bifurcation is supercritical, so there exists only one such state. At
$Re \approx 286$
the
$A_yS_z$
mode becomes unstable, and this second bifurcation is subcritical: the branch
$(0,B)$
is initially unstable and moves towards lower Reynolds numbers, before becoming stable again and folding back to a larger Reynolds number. One mixed state
$(A,B)$
exists but is unstable. The state
$(A,0)$
becomes unstable at
$Re \approx 287$
, leaving
$(0,B)$
as the only stable state. Therefore, there is a narrow interval of bistability between the vertically and horizontally deflected states.
4. Three-dimensional nonlinear simulations
In this section we perform 3-D, unsteady, fully nonlinear simulations to explore the successive bifurcations observed in the flow past elongated prisms with
$L=5$
and
$W=1.2$
,
$2.25$
and
$5$
, up to
$Re=700$
. The specific length
$L=5$
is chosen as it defines the international BARC benchmark (Bruno et al. Reference Bruno, Salvetti and Ricciardelli2014). To distinguish the different regimes, we use the notation
$xY_yZ_z$
: the lower case
$x$
determines whether the regime is steady (
$s$
), periodic (
$p$
) or aperiodic (
$a$
); the upper case
$Y$
and
$Z$
determine the planar symmetries of the flow, i.e. whether the flow is symmetric (
$S$
) or antisymmetric (
$A$
) with respect to the
$y=0$
and
$z=0$
planes. For unsteady flows, the symmetry refers to the time-averaged flow. For example, the
$sS_yS_z$
regime refers to a steady regime that retains the left/right and top/bottom planar symmetries (like the symmetric base flow of § 3.1). For
$W=2.25$
and
$5$
, the periodic regimes
$pS_yS_zl$
and
$pS_yS_zt$
are characterised by vortex shedding from the LE (
$l$
) and TE (
$t$
) shear layers, respectively.
We use the aerodynamic forces as driving quantities to detect the different flow regimes. Indeed, the frequency spectra and the instantaneous/mean/root-mean-square values of the aerodynamics forces provide immediate information about the spatial and temporal symmetries of the flow, i.e. whether the flow is steady/periodic/aperiodic and whether the instantaneous and mean fields exhibit top/bottom or right/left planar symmetries. However, we have verified that the same information regarding the temporal symmetries of the flow is obtained by inspecting the velocity field and probing time signals at different locations of the flow. To further characterise the flow regimes, we also report visualisations of mean and instantaneous fields, phase-space diagrams and proper orthogonal decomposition (POD) modes.
4.1.
Narrow prism:
$W=1.2$

Figure 10. Dependence of the aerodynamic forces on the Reynolds number for
$L=5$
and
$W=1.2$
. (a) Streamwise force
$F_x$
, where the red line denotes
$8.47 \times Re^{-0.3799}$
. Here circles refer to the average value and bars to the root mean square of the fluctuations. (b) Lateral forces
$F_y$
and
$F_z$
. (c) Frequency spectra of
$F_x$
(left),
$F_y$
(centre) and
$F_z$
(right) in logarithmic scale.
For
$L=5$
and
$W=1.2$
, five different regimes are identified when the Reynolds number is increased up to
$Re=700$
. This is conveniently visualised in figure 10, where the aerodynamic forces are shown as a function of
$Re$
. For
$Re \lessapprox 352$
, the flow is steady and retains the
$S_yS_z$
symmetry, in agreement with the LSA (§ 3.2). At
$Re \approx 355$
the
$sS_yS_z$
regime is unstable and the flow experiences a pitchfork bifurcation towards the
$sS_yA_z$
regime. The flow loses the top/bottom planar symmetry but retains the left/right one. For
$Re \gtrapprox 370$
, the flow enters the
$sA_yS_z$
regime, i.e. recovers the top/bottom planar symmetry and loses the left/right one. This regime remains stable up to
$Re \approx 500$
, when the flow becomes unsteady. The flow first approaches the periodic
$pA_yS_z$
regime, characterised by an alternating shedding of HVs from the top/bottom LE shear layers. For
$Re \geqslant 515$
, the wake oscillates and the flow is in the aperiodic
$aA_yS_z$
regime.
4.1.1. Low
$Re$
: steady
$sS_yA_z$
and
$sA_yS_z$
regimes
The nonlinear simulations show that the primary bifurcation breaks the top/bottom planar symmetry leading to the steady
$sS_yA_z$
regime, similarly to shorter bodies (Marquet & Larsson Reference Marquet and Larsson2015; Zampogna & Boujo Reference Zampogna and Boujo2023). This bifurcation occurs at
$Re \approx 355$
, in agreement with the critical Reynolds number
$Re_{c}=353$
found in § 3.2. The simulations show that, for
$Re \gtrapprox 365$
, the steady
$A_yS_z$
mode dominates and the flow recovers the top/bottom planar symmetry but loses the left/right one, as predicted by the WNL analysis. We observe a small discrepancy between the linear and WNL stability analyses (
$Re_{c} \approx 355$
) and the fully nonlinear simulations (
$Re_{c} \approx 365$
).

Figure 11. Steady regimes for
$L=5$
and
$W=1.2$
. Isosurfaces of streamwise vorticity, with red/blue indicating
$\omega _x=\pm 0.075$
. Left:
$sS_yA_z$
regime at
$Re=355$
. Right:
$sA_yS_z$
regime at
$Re=370$
. Top: horizontal
$x-y$
plane. Centre: vertical
$x-z$
plane. Bottom:
$y-z$
plane.
To characterise the
$sA_yS_z$
and
$sS_yA_z$
regimes, figure 11 shows isosurfaces of positive (red) and negative (blue) values of the streamwise vorticity, and illustrates the loss of the planar symmetries. This is similar to the wakes past disks, spheres and bullet-shaped bodies at Reynolds numbers larger than
$Re \approx 115, 210$
and
$350$
(Johnson & Patel Reference Johnson and Patel1999; Tomboulides & Orszag Reference Tomboulides and Orszag2000; Fabre et al. Reference Fabre, Auguste and Magnaudet2008; Bohorquez et al. Reference Bohorquez, Sanmiguel-Rojas, Sevilla, Jiménez-González and Martínez-Bazán2011). In these steady regimes, the wake features a pair of counter-rotating streamwise vortices that are absent at lower
$Re$
. These vortices exhibit an eccentricity that increases with
$x$
. The eccentricity is in the
$z$
direction in the
$sS_yA_z$
regime and in the
$y$
direction in the
$sA_yS_z$
regime.
In both the
$sS_yA_z$
and
$sA_yS_z$
regimes we find that the flow asymmetry is confined in the wake, confirming that the
$A_yS_z$
and
$S_zA_y$
global modes are unstable modes of the wake. Indeed, the recirculating regions that arise along the sides of the prism retain the top/bottom and left/right planar symmetries in both regimes (not shown for brevity). These regions do not play a relevant role in the triggering mechanism of the primary instability, as confirmed by the results of the structural sensitivity (figure 8). A similar conclusion has been drawn for the primary, Hopf instability of the flow past elongated 2-D rectangular cylinders (Chiarini et al. Reference Chiarini, Quadrio and Auteri2021).

Figure 12. Dependence of the aerodynamic forces on
$Re$
in the
$sS_yS_z$
,
$sS_yA_z$
and
$sA_yS_z$
regimes for
$L=5$
and
$W=1.2$
. (a,b) Total forces. (c) Viscous and pressure components of the drag force. (d) Viscous and pressure components of
$F_\perp$
, the single non-zero force perpendicular to the incoming flow (
$F_z$
in the
$sS_yA_z$
regime and
$F_y$
in the
$sA_yS_z$
regime). The red line in (a) denotes
$8.47 \times Re^{-0.3799}$
.
Figure 12 shows that when the flow bifurcates to the
$sS_yA_z$
and
$sA_yS_z$
regimes, the drag decreases according to a power law,
$F_x \sim Re^{\alpha }$
with
$\alpha =-0.3799$
. This is similar to the flow past a sphere (Johnson & Patel Reference Johnson and Patel1999), a short circular cylinder (Yang et al. Reference Yang, Feng and Zhang2022) and a cube (Saha Reference Saha2004; Meng et al. Reference Meng, An, Cheng and Kimiaei2021); for the latter, Saha (Reference Saha2004) report
$\alpha =-0.372$
and Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) found a value between
$-0.4$
and
$-0.3$
. The decrease of
$F_x$
is entirely due to a decrease of the friction contribution
$F_{x,f}$
(see figure 12(c)), as the pressure contribution
$F_{x,p}$
does not change with
$Re$
. An increase of
$Re$
, indeed, results into a downstream shift of the flow reattachment point
$x_r$
over the sides of the prism, leading to a longer recirculating region and to a decrease of
$F_{x,f}$
(figure 13). We note that a power function fits well the relation between
$x_r$
and
$Re$
, being
$x_r \sim Re^{1.1651}$
for the left and right sides and
$x_r \sim Re^{0.9485}$
for the top and bottom sides. In contrast, the non-zero cross-stream component of the aerodynamic forces increases with
$Re$
. Unlike for
$F_x$
, this is due to an increase of both the friction and pressure contributions (see the bottom right panel of figure 12), and is associated with the increasing flow asymmetry.
The regular bifurcation towards an asymmetric steady flow is found for prisms with small
$W$
only, irrespective of
$L$
(§ 3.2). This is the case of many other 3-D bluff bodies: besides the largely studied axisymmetric sphere and disk, the existence of a regular bifurcation in non-axisymmetric flows has been reported by Marquet & Larsson (Reference Marquet and Larsson2015) for thin plates with
$W/D\lt 2$
, by Yang et al. (Reference Yang, Feng and Zhang2022) for short circular cylinders with
$W/D\lt 1.75$
and by Sheard et al. (Reference Sheard, Thompson and Hourigan2008) for short cylinders with hemispherical ends. The main difference is that the bifurcated wake retains a symmetry plane that is randomly oriented (selected by perturbations; Tomboulides & Orszag Reference Tomboulides and Orszag2000) for axisymmetric bodies, but dictated by the planar symmetries of the geometry for non-axisymmetric bodies.
4.1.2. Large
$Re$
:
$pA_yS_z$
and
$aA_yS_z$
regimes

Figure 13. Position of the reattachment point
$x_r$
on the top/bottom and right/left sides (i.e. size of the recirculating regions on the prism walls) as a function of
$Re$
for
$L=5$
and
$W=1.2$
. Red solid line:
$0.0023 \times Re^{1.1651}$
. Red dashed line:
$0.0095 \times Re^{0.9485}$
. In the unsteady regimes (
$Re\gt 500$
),
$x_r$
is the reattachment point of the time-averaged flow.
A further bifurcation occurs at
$Re \approx 510$
. It consists of a Hopf bifurcation of the
$sA_yS_z$
deflected wake; the flow becomes unsteady and starts oscillating about the deflected
$sA_yS_z$
regime. A similar regime, with the flow oscillating about an asymmetric steady state has been observed for the flow past a sphere (Johnson & Patel Reference Johnson and Patel1999), a cube (Saha Reference Saha2004), a disk (Meliga et al. Reference Meliga, Chomaz and Sipp2009) and a short cylinder (Pierson et al. Reference Pierson, Auguste, Hammouti and Wachs2019). Interestingly, when considering the time-average forces, the decreasing power law that fits
$F_x$
and
$Re$
in the
$sS_yS_z$
,
$sS_yA_z$
and
$sA_yS_z$
regimes holds also in the unsteady regime for intermediate Reynolds numbers up to
$Re=575$
(figure 10). For these
$Re$
, indeed, the time-average flow field -- and its dependence on
$Re$
-- resembles the steady
$sA_yS_z$
regime. The average size of the recirculating regions over the lateral sides increases with the Reynolds number, as shown in figure 13 by means of the reattaching point
$x_r$
. For larger
$Re$
, the structure of the time-average flow changes and
$F_x$
increases. For
$Re \geqslant 575$
, indeed, the average
$x_r(Re)$
deviates from the power law found at smaller
$Re$
: the lateral recirculating regions shrink as
$Re$
increases (figure 13).
The dynamics of the flow progressively changes with
$Re$
. At
$Re=510$
the flow enters the periodic
$pA_yS_z$
regime, which is characterised by an alternating shedding of HVs from the top and bottom LE shear layers (figure 14). At this
$Re$
, a single peak
$St \approx 0.277$
is found in the frequency spectrum
$S(F_z)$
, while
$F_y$
remains practically constant in time. The attractor draws a limit cycle in the phase space.
For
$Re\gtrapprox 515$
, a frequency
$St \approx 0.0776$
appears in
$S(F_y)$
and the flow becomes aperiodic. Here it is characterised by a superposition of two different modes, (i) the shedding of HVs from the top and bottom LE shear layers (with
$St \approx 0.28$
), and (ii) the oscillating motion of the wake in the
$y$
direction (with
$St \approx 0.08$
); see figure 15 and the related discussion. The two frequencies are not commensurate and, as shown in the force diagrams of figure 14, a torus replaces the limit cycle in the phase space. At this
$Re$
, the peak at
$St \approx 0.28$
is much larger than the other one, indicating that the shedding of HVs dominates.
For
$Re \geqslant 550$
, the wake oscillates also in the
$z$
direction and a further frequency
$St \approx 0.03$
appears in
$S(F_z)$
(figure 14). The different
$Re$
at which the two oscillating modes arise agree with the results of the LSA, which reveals that, for
$L \leqslant 2$
, the unsteady
$A_yS_z$
wake mode becomes unstable at smaller
$Re$
compared with the unsteady
$S_yA_z$
mode. As mentioned above, the presence of this additional mode leads to an increase of the average
$F_x$
(figure 10). When
$Re$
is further increased, the different modes nonlinearly interact and the flow becomes progressively more chaotic (figure 14). The relative height of the peaks at
$St \approx 0.28$
and
$St \approx 0.03$
in
$S(F_z)$
indicates that the flow dynamics is mainly driven by the LE vortex shedding at low
$Re$
, while the vertical wake oscillating motion takes over at large
$Re$
(not shown).

Figure 14. Unsteady regimes for
$L=5$
and
$W=1.2$
. (a) Frequency spectra of
$F_y$
(left) and
$F_z$
(right),
$510 \leqslant Re \leqslant 550$
. (b–c) Structure of the flow for
$Re=510$
, 535 and 575 (top to bottom). (b) Instantaneous isosurfaces
$\lambda_2 = -0.05$
coloured by streamwise vorticity in the range
$-1 \leqslant \omega_x \leqslant 1$
. (c) Force diagrams
$F_y-F_x$
(left) and
$F_z-F_x$
(right).

Figure 15. The POD analysis for
$L=5$
,
$W=1.2$
and
$Re=515$
. (a) Energy fraction
$E_j = \lambda _j/\sum \lambda _j$
of the first eight POD modes. Here
$\lambda _j$
denotes the
$j$
th eigenvalue of the snapshot covariance matrix, and is related with the
$j$
th singular value by
$\sigma _j = \lambda _j^2$
(see Appendix C). (b) Frequency associated with the first eight POD modes, ordered by decreasing
$\sigma _j$
(i.e. decreasing
$E_j$
). Dashed lines: main DNS frequencies, also shown in figure 10. (c) Visualisation of the POD modes 1 (left) and 3 (right). Isosurfaces of
$\lambda _2$
(arbitrary value) coloured by streamwise vorticity (symmetric blue-to-red colour map from negative to positive values).
To examine the spatial structure of the different modes and relate them with the flow frequencies, we use POD; see Appendix C for details. Figure 15 considers
$Re=515$
. As shown in Figure 15(a,b), the dominant POD modes are associated with the same frequencies found in the frequency spectra (figure 14) and with their multiples. The POD modes come in pairs (i.e. with the same singular values and frequencies), as typical for oscillating structures. Mode 1 is associated with
$St \approx 0.28$
, and a shedding of vortices from the top and bottom sides is clearly visible in its spatial structure. Mode 3, instead, is associated with
$St \approx 0.07$
. Its spatial structure is confined in the wake and breaks the left/right symmetry consistently with an oscillating mode of the wake in the
$y$
direction about a
$A_yS_z$
state. At this
$Re$
the energy fraction associated with the LE vortex shedding is approximately
$1.5$
times larger than that associated with the lateral wake oscillation (see Figure 15
a), consistently with the peak heights in the frequency spectra.
For elongated prisms with
$W=1.2$
and
$L=5$
, the characteristic frequency of the wake oscillation (
$St \approx 0.07$
) is smaller than for a sphere (
$St \approx 0.129$
at
$Re=270$
; Tomboulides & Orszag Reference Tomboulides and Orszag2000), a cube (
$St \approx 0.091-0.095$
at
$Re= 270-300$
; Saha Reference Saha2004), and a short cylinder (
$St \approx 0.125$
at
$Re=283$
; Yang et al. Reference Yang, Feng and Zhang2022). This is consistent with the LSA results (figure 5) that point out a decrease of
$St$
with
$L$
.
4.2.
Intermediate width:
$W=2.25$
For
$L=5$
and
$W=2.25$
, the scenario becomes more complicated, with many regimes identified up to
$Re = 700$
(figure 16). In agreement with the LSA (§ 3.2), the critical
$Re$
corresponding to the onset of the first bifurcation decreases with the width of the prism. For
$W=2.25$
, the primary instability consists of a pitchfork bifurcation towards the steady
$sS_yA_z$
regime, for
$Re$
between
$270$
and
$290$
. At
$Re \approx 300$
, another bifurcation restores the top/bottom planar symmetry and breaks the left/right one, thus leading to the
$sA_yS_z$
regime, similar to smaller
$W$
. The simulations show an hysteresis in the transition from
$sS_yA_z$
to
$sA_yS_z$
. Bistability is detected for some values of
$Re$
, in agreement with the subcritical nature of the bifurcation (see the WNL analysis in § 3.4). Interestingly, at larger
$Re$
the flow recovers both planar symmetries and settles again in a
$sS_yS_z$
regime, before becoming unsteady at
$Re \approx 375$
with a combination of two different modes: (i) oscillation of the wake and (ii) symmetric shedding of HVs from the top and bottom LE shear layers. As
$Re$
further increases, the wake oscillation mode stabilises, leading to the periodic
$pS_yS_zla$
and
$pS_yS_zlb$
regimes, and then destabilises again at larger
$Re$
, leading eventually to the aperiodic
$aS_yS_z$
regime.

Figure 16. As figure 10, for
$L=5$
and
$W=2.25$
. Red line in the top panel:
$4.33 \times Re^{-0.2905}$
.

Figure 17. Steady regimes for
$L=5$
and
$W=2.25$
. Back view of the
$y-z$
plane. The red and blue isosurfaces are for
$\omega _x= \pm 0.1$
. Results are shown for
$(a)$
$Re=300$
, regime
$sS_yA_z$
;
$(b)$
$Re=300$
, regime
$sA_yS_z$
;
$(c)$
$Re=330$
, regime
$sA_yS_z$
;
$(d)$
$Re=360$
, regime
$sS_yS_z$
. (e) Zoom of figure 16(b) in the range
$280 \leqslant Re \leqslant 370$
.
4.2.1. Low
$Re$
: steady regimes
$sS_yA_z$
and
$sA_yS_z$
The nonlinear simulations show that the flow is steady and retains the
$S_yS_z$
symmetry for
$Re \lessapprox 290$
, in agreement with the results of the LSA. At
$Re \approx 290$
the flow undergoes a regular bifurcation towards the steady
$sS_yA_z$
regime, qualitatively similar to
$W=1.2$
. Further increasing
$Re$
, the flow remains steady but recovers the top/bottom symmetry and loses the left/right one, i.e. enters the
$sA_yS_z$
regime. Unlike for smaller
$W$
, this bifurcation is hysteretic, and both the
$sS_yA_Z$
and
$sA_yS_z$
regimes are simultaneously stable for a small range of
$Re \approx 300$
; see figure 17(a--d) and the WNL analysis in § 3.4.
By increasing the Reynolds number in the
$320 \leqslant Re \leqslant 350$
range, the flow asymmetry in the
$y$
direction decreases, and at
$Re=350$
the flow recovers the left/right planar symmetry; for
$Re \approx 360$
, the flow is steady and
$F_y \approx F_z \approx 0$
. This symmetrisation is shown in figure 17. Like for
$W=1.2$
(figure 11), at
$Re=300$
a pair of counter-rotating streamwise vortices are found in the wake, with eccentricity in the
$y$
direction (figure 17
b). As
$Re$
increases, a new pair of counter rotating vortices arises, becomes stronger and progressively pushes the other pair away (figure 17
c). Eventually, at
$Re \approx 360$
the two pairs of vortices are placed symmetrically with respect to the
$z$
axis, and have the same intensity (figure 17
d). At this
$Re$
the flow exhibits again both top/bottom and left/right planar symmetries (
$sS_yS_z$
regime). This does not happen for
$W=1.2$
, as two pairs of streamwise vortices of characteristic size
$H=1$
cannot be simultaneously accommodated at the TE when
$W\lt 2$
.
In passing, we note that in the steady
$sS_yA_z$
,
$sA_yS_z$
and
$sS_yS_z$
regimes a decreasing power law fits well the relation between
$F_x$
and
$Re$
, where
$F_x \sim Re^{\alpha }$
with
$\alpha =-0.2905$
. However, the agreement does not extend to the unsteady regime, unlike
$W=1.2$
. Also, like for smaller
$W$
, the decrease of
$F_x$
with
$Re$
is mainly due to the enlargement of the side recirculating regions (not shown).
4.2.2. Intermediate
$Re$
: aperiodic
$aS_yA_z$
and periodic
$pS_yA_z$
regimes

Figure 18. Forces and frequency content for
$L=5$
,
$W=2.25$
and
$375 \leqslant Re \leqslant 450$
. (a) Frequency spectra for
$F_x$
(left) and
$F_z$
(right) for
$Re=375, 390, 420$
and
$450$
. For
$Re=420$
,
$S(F_z)$
is not visible as fluctuations in
$F_z$
are almost null. (b) Main frequencies for
$F_x$
(left) and
$F_z$
(right). (c) Zoom of figure 16(a,b) in the range
$370 \leqslant Re \leqslant 460$
. For
$440 \leqslant Re \leqslant 460$
, the bars are not visible as the oscillations of
$F_x$
are small.

Figure 19. Structure of the flow for
$L=5$
,
$W=2.25$
and
$Re=375$
,
$390$
,
$420$
and
$450$
(top to bottom). (a) Side view of instantaneous isosurfaces of
$\lambda _2=-0.25$
coloured with
$-1 \leqslant \omega _x \leqslant 1$
. (b) Force diagrams
$F_y-F_x$
and
$F_z-F_x$
. The lack of top/bottom symmetry in the top two views of panel (a) is highlighted with arrows.
At
$Re \approx 375$
the flow becomes unsteady and enters different regimes in short succession. This is conveniently visualised in figures 18 and 19, where the frequency content and the structure of the flow are shown for
$375 \leqslant Re \leqslant 470$
.
For
$375 \leqslant Re \leqslant 380$
, the flow approaches the aperiodic
$aS_yA_z$
regime: unlike for smaller
$W$
, the flow oscillates about a
$S_yA_z$
state that retains the left/right symmetry and breaks the top/bottom one. With two non-commensurate frequencies,
$St \approx 0.05$
and
$0.23$
, the attractor draws a torus in the phase space. As shown in the top panel of figure 19(a), at this
$Re$
the flow behaviour is the result of the superposition of two different modes: (i) asymmetric oscillating mode of the wake in the
$z$
direction (
$St \approx 0.05$
) and (ii) symmetric in-phase shedding of HVs from the top and bottom LE shear layers (
$St \approx 0.23$
). For
$Re \approx 385$
, instead, the wake oscillation mode stabilises and the flow approaches the periodic
$pS_yA_z$
regime. The wake remains deflected in the
$z$
direction and the unsteadiness is only due to the in-phase shedding of HVs, with a frequency
$St \approx 0.24$
that slightly increases with
$Re$
(see figure 18
c). A limit cycle replaces the torus in the phase space (figure 19).

Figure 20. The POD analysis for
$L=5$
,
$W=2.25$
and
$Re=380$
. (a) Energy fraction
$E_j = \lambda _j/\sum \lambda _j$
of the first eight POD modes. (b) Frequency associated with these modes. Dashed lines: main DNS frequencies as in figure 18. (c) The POD modes 1 (left) and 5 (right). Isosurfaces of
$\lambda _2$
(arbitrary value) coloured by streamwise vorticity.
We use POD to examine the structure of the flow at
$Re=380$
. In figure 20 we separate the modes responsible for the two non-commensurate frequencies detected in the frequency spectra. Here POD mode
$1$
is associated with
$St \approx 0.23$
and its spatial structure shows a
$S_yS_z$
shedding of HVs from the top and bottom sides of the prism. Mode
$5$
, instead, matches the
$St \approx 0.05$
frequency, with a
$S_yA_z$
spatial structure that agrees with an asymmetric oscillating mode of the wake in the
$z$
direction. This suggests that the aperiodic
$aS_yA_z$
regime found at
$Re \approx 380$
is the result of a superposition of an oscillatory
$S_yA_z$
mode of the wake and an unsteady
$S_yS_z$
mode of the LE shear layers. Interestingly, this interpretation is supported by the LSA, which indeed predicts an unsteady
$S_yS_z$
mode to become unstable for
$Re \gtrapprox 350$
for this geometry, with a frequency that is compatible with
$St \approx 0.2$
(figure 5). Also, the neutral curves (figure 5) show that, for
$L\leqslant 4$
, an unsteady
$S_yA_z$
mode of the wake with
$St \approx 0.06-0.07$
becomes unstable, but for a small range of
$Re$
only. For
$L \approx 4$
, for example, the unstable range is between
$330 \lessapprox Re \lessapprox 350$
. Recalling that these neutral curves refer to the low-
$Re$
$sS_yS_z$
steady base flow, it is possible that due to nonlinear effects, an increase of
$Re$
extends this curve to
$L=5$
, explaining this sudden destabilisation/stabilisation of the wake mode.
Figure 20 shows that at
$Re=380$
the energy fraction associated with LE vortex shedding is approximately
$25$
times larger than that associated with the asymmetric oscillating mode of the wake; at this
$Re$
the flow dynamics is mainly driven by the HVs shed by the LE shear layers. The effect of the Reynolds number on the spatial structure of the POD mode associated with the shedding of HVs from the LE shear layers is further discussed in Appendix C.2 (figure 32).
4.2.3. Intermediate
$Re$
: periodic
$pS_yS_zl$
regime
As
$Re$
is further increased, the flow recovers once again the top/bottom planar symmetry and, for
$390 \lt Re \leqslant 470$
, it is in the periodic and symmetric
$pS_yS_zl$
regime (figure 19). In this regime, two distinct behaviours can be further distinguished, based on the fluctuations of
$F_z$
, characterised by a synchronisation/anti-synchronisation of the vortex shedding from the top and bottom LE shear layers.
For
$390 \lt Re \leqslant 420$
, the flow shows an in-phase shedding of HVs from the top and bottom LE shear layers, and the flow instantaneously retains the top/bottom and left/right planar symmetries (third row of figure 19); this is regime
$pS_yS_zla$
. Accordingly, the unbalance of the viscous and pressure forces between the top/bottom and left/right sides is null at all times, leading to
$F_z \approx F_y \approx 0$
. For this regime, we detect the leading flow frequency by inspecting the
$S(F_x)$
spectrum. To the best of our knowledge, such a synchronous vortex shedding that preserves all spatial symmetries at all times and generates neither lift force nor side force has not been reported.
For
$420 \lt Re \leqslant 470$
, instead, the flow remains periodic with
$St \approx 0.24$
, but the shedding of vortices from the top and bottom LE shear layers is in phase opposition (last row of figure 19); this is regime
$pS_yS_zlb$
. In this regime the flow retains the top/bottom planar symmetry in an average sense only, and the oscillations of
$F_z$
are not null. Accordingly, unlike for
$380 \leqslant Re \leqslant 420$
, at these
$Re$
the POD mode associated with the LE vortex shedding is top/down antisymmetric (see Appendix C.2).
We now consider regime
$pS_yS_zla$
and detail the dynamics of the HVs shed from the top and bottom LE shear layers at
$Re \approx 400$
.

Figure 21. Flow structures of the doubly symmetric
$pS_yS_zla$
regime for
$L=5$
,
$W=2.25$
and
$Re=400$
: isosurfaces of
$\lambda _2 = -0.05$
coloured by
$-1 \leqslant \omega _x \leqslant 1$
. The four snapshots are separated in time by
$T/4$
, where
$T$
is the period of the shedding of HVs.
Figure 21 shows the flow structures at four different instants equispaced in one shedding period
$T=1/St$
. A pair of HVs of opposite vorticity is shed in phase from the top and bottom LE shear layers, with a spanwise wavelength
$\lambda _z$
dictated by the width of the prism
$\lambda _z \approx W$
. The top/bottom LE shear layers, indeed, periodically release tubes of negative/positive spanwise vorticity that, due to the finite width of the prism, are not exactly aligned with the spanwise direction. As such, they are modulated by the vertical and spanwise velocity gradients and form HVs: the central part of the tubes, farther from the wall, is convected downstream faster and forms the HV heads; the lateral parts, closer to the wall, are convected more slowly and form the HV legs. Accordingly, the HVs induce low-speed velocity (negative velocity fluctuations) between their legs and high-speed velocity (positive velocity fluctuations) in the outer region. Once generated, the two HVs are shed downstream and continue moving as a pair in the wake, where they progressively move apart in the
$z$
direction due to their mutual induction, until they are dissipated by viscosity. The shedding of HVs is accompanied by a periodic enlargement and shrinking of the top/bottom recirculation regions (not shown). The recirculating regions enlarge while the LE shear layers accumulate vorticity, and suddenly shrink when the HVs are shed.
The shedding of HVs from the LE shear layers is reminiscent of the flow past 2-D rectangular cylinders: spanwise vortices are periodically shed from the LE shear layers due to the so-called impinging-LE-vortex (ILEV) instability (Naudascher & Rockwell Reference Naudascher and Rockwell1994). The ILEV instability is a resonant oscillation of the fluid. Vortices are shed periodically from the LE shear layer and, when a LE vortex passes over the TE, it triggers the shedding of a new LE vortex (Chiarini et al. Reference Chiarini, Quadrio and Auteri2022d ). The flow frequency thus changes with the length of the body and the dynamics of the LE and TE vortex shedding are synchronised. In the 2-D case, therefore, the shedding of vortices from the LE is not the result of an absolute instability of the LE shear layer, but is interconnected with the dynamics of the TE vortex shedding (Hourigan et al. Reference Hourigan, Thompson and Tan2001; Chiarini et al. Reference Chiarini, Quadrio and Auteri2022d ).

Figure 22. Lateral view of the flow structures for rectangular prisms with (a)
$(L,W)=(\infty,2.25)$
and a sharp LE and (b)
$(L,W)=(5,2.25)$
and a rounded LE with
$R = 1/64$
at
$Re=400$
. Isosurfaces of
$\lambda _2 = -0.05$
coloured by
$-1 \leqslant \omega _x \leqslant 1$
. The flow approaches the regimes
$pS_yS_zla$
and
$pS_yA_z$
, respectively.
To investigate the role of the LE and TE shear layers in the mechanism that sustains the
$pS_yS_zla$
regime, two types of additional simulations have been carried out at
$Re=400$
for modified prism geometries. First, to isolate the LE shear layer from the interaction with the TE, we have considered a rectangular prism of infinite length and
$W=2.25$
. As shown in figure 22(a), the flow enters the
$pS_yS_zla$
regime in this case too, with pairs of HVs shed in phase from the top and bottom LE shear layers. Interestingly, the shedding frequency closely matches that found for
$L=5$
,
$St \approx 0.24$
. This indicates that, for 3-D prisms, the triggering mechanism of the LE vortex shedding does not require the presence of a TE, as instead observed for 2-D blunt bodies (Chaurasia & Thompson Reference Chaurasia and Thompson2011; Thompson Reference Thompson2012). A possible feedback mechanism that triggers the formation of the HVs may therefore be embedded within the recirculating regions that arise over the top and bottom sides of the prism.
Second, to support this hypothesis, we have considered rectangular prisms with
$L=5$
and
$W=2.25$
, but with a rounded LE. Five curvature radii
$R$
have been considered, ranging from
$R=1/2$
(semicircular LE) to
$R=1/64$
. For
$R=1/2$
, the flow does not separate from the LE, and recirculating regions do not form over the top and bottom sides of the prism. In this case the flow approaches a steady and asymmetric
$sS_yA_z$
regime and shedding of HVs from the LE is not detected. When a smaller
$R$
is considered, the flow separates from the LE and progressively larger recirculating regions form over the top and bottom sides of the prism. However, the shedding of HVs is only detected for the smallest radius
$R=1/64$
. In this case the flow enters a periodic regime around the asymmetric
$S_yA_z$
state, characterised by the in-phase shedding of HVs from the top and bottom LE shear layers (figure 22
b). Again, the shedding frequency matches well that found for the sharp configuration,
$St \approx 0.238$
.
In summary, these additional simulations support the hypothesis that in the
$pS_yS_zl$
regime the shedding of HVs from the LE shear layers is the result of a feedback mechanism that takes place in the recirculating regions over the top/bottom sides of the prism and does not require the presence of a TE.
4.2.4. Large
$Re$
: aperiodic regime
$aS_yS_z$

Figure 23. Unsteady regimes at larger
$Re$
for
$L=5$
and
$W=2.25$
. (a) Frequency spectra of
$F_y$
(left) and
$F_z$
(right) for
$470 \leqslant Re \leqslant 625$
. (b–c) Structure of the flow for
$Re=470$
,
$485$
and
$625$
(top to bottom). (b) Instantaneous isosurfaces of
$\lambda _2=-0.05$
coloured by
$-1 \leqslant \omega _x \leqslant 1$
. (c) Force diagrams
$F_y-F_x$
(left) and
$F_z-F_x$
(right).
For larger Reynolds numbers
$Re \geqslant 460$
, the lateral and vertical oscillating modes of the wake become unstable, and the flow progressively enters a chaotic regime (figure 23). Unlike for smaller
$W$
, for
$500 \leqslant Re \leqslant 700$
, the flow oscillates around a
$S_yS_z$
state, and the time-average value of the cross-stream forces is
$F_y=F_z=0$
. This recalls the results of Zdravkovich et al. (Reference Zdravkovich, Brand, Mathew and Weston1989) and Zdravkovich et al. (Reference Zdravkovich, Flaherty, Pahle and Skelhorne1998), who investigated the flow past a finite circular cylinder at
$Re \approx 10^5$
, varying the width-to-height ratio between
$0.25 \leqslant W \leqslant 10$
. In agreement with our results, they found an asymmetric flow pattern for
$1 \leqslant W \leqslant 2$
only.
According to our simulations, the flow loses the instantaneous left/right symmetry at
$Re \approx 465$
. For
$465 \leqslant Re \leqslant 470$
,
$S(F_y)$
(figure 23) shows a single peak at
$St \approx 0.134$
, which is half the frequency of the LE vortex shedding (see
$S(F_z)$
in the same figure). The flow enters a
$pS_yS_z$
regime. In this range of
$Re$
the flow remains periodic and the attractor is a limit cycle that draws a closed line in the phase space (see the force diagrams in figure 23). As
$Re$
increases, the bifurcated limit cycle becomes unstable. At
$Re=485$
, a new incommensurate frequency
$St \approx 0.079$
arises in
$S(Fy)$
; the flow is therefore aperiodic and the attractor draws a torus in the phase space.
For
$Re \geqslant 500$
, additional frequencies appear and the flow becomes chaotic. In the
$aS_yS_z$
regime the flow dynamics is dominated by the nonlinear interaction of three different modes: (i) shedding of HVs from top and bottom LE shear layers; (ii,iii) oscillating modes of the wake in the
$z$
and
$y$
directions. The frequency of the first two modes are detected with the peaks in
$S(F_z)$
, i.e.
$St \approx 0.25{-}0.30$
and
$St \approx 0.02{-}0.03$
. The frequency
$St \approx 0.04{-}0.05$
of the wake oscillating mode in the
$y$
direction, instead, is visible in
$S(F_y)$
. As
$Re$
increases, the frequency of the LE vortex shedding increases and the corresponding peak in the frequency spectrum becomes less evident and more broad-banded (not shown). Like for smaller
$W$
, indeed, at larger
$Re$
the LE vortex shedding weakens and the flow unsteadiness is mainly driven by the wake oscillating modes. Note that
$St \approx 0.03$
is close to what was found for
$W=1.2$
, indicating that the frequency of the wake oscillating mode in the
$z$
direction is only marginally influenced by the width of the prism for
$1.2 \leqslant W \leqslant 2.25$
. In contrast, the frequency of the wake oscillation mode in the
$y$
direction clearly decreases with
$W$
. Figure 23 shows that, for
$Re \gt 500$
, the shedding of HVs arises also over the lateral sides of the prism.
4.3.
Wide prism:
$W=5$
For
$L=5$
and
$W=5$
, three main regimes have been identified for
$Re \lt 700$
(figure 24). In agreement with the LSA, at
$Re \approx 250$
the steady
$S_yS_z$
flow experiences a Hopf bifurcation towards the periodic
$pS_yS_zt$
regime, which is characterised by the unsteady
$S_yA_z$
mode of the wake. For
$Re \gtrapprox 305$
, the flow approaches the aperiodic
$aS_yS_z$
regime. Here the wake oscillates in both the
$y$
and
$z$
directions, and vortices are shed from the top/bottom LE shear layers. The interaction between these modes largely changes with the Reynolds number. For some
$Re$
, wake oscillations and vortex shedding synchronise and the flow recovers periodicity.
4.3.1. Low
$Re$
: periodic regime
$pS_yS_zt$

Figure 25. Characterisation of the
$pS_yS_zt$
regime for
$L=5$
and
$W=5$
. (a) Isosurfaces of spanwise vorticity (red and blue for
$\omega _y= \pm 0.25$
) at
$Re=275$
. (b) Friction drag for
$225 \leqslant Re \leqslant 300$
. (c) Pressure drag for
$225 \leqslant Re \leqslant 300$
. (d) Flow frequency for
$250 \leqslant Re \leqslant 300$
.
The nonlinear simulations show that the unsteady
$S_yA_z$
mode of the wake becomes unstable at
$Re \approx 250$
and that, for
$250 \lessapprox Re \lessapprox 300$
, the flow experiences a periodic oscillation of the wake in the
$z$
direction around the steady
$S_yS_z$
state (figure 25
a). Unlike for smaller
$W$
, here the flow unsteadiness is driven by the wake oscillation and HVs are not shed by the LE shear layers. The hairpin-like structures observed in the wake are thus the result of an interaction between the top and bottom TE shear layers; this resembles the so-called hairpin shedding regime found for shorter bluff bodies of various shapes (Tomboulides & Orszag Reference Tomboulides and Orszag2000; Saha Reference Saha2004; Yang et al. Reference Yang, Feng and Zhang2022). The flow is periodic with a single frequency close to that found with LSA (
$St=0.09$
; § 3.2). As expected, the intensity of the wake oscillation increases with
$Re$
, as conveniently visualised in figure 24 by means of the root mean square of
$F'_z = F_z - \! \left \langle {F_z} \right \rangle \!_t$
(
$ \! \left \langle {\cdot } \right \rangle \!_t$
indicates average in time):
$F'_{z,rms} \approx 0.004$
at
$Re=250$
and
$F'_{z,rms} \approx 0.036$
at
$Re=300$
. In this regime, an increase of
$Re$
is accompanied by an increase of
$F_x$
, which is entirely ascribed to a monotonic increase of the pressure contribution
$F_{x,p}$
. In fact, similarly to
$W=1.2$
(figure 13), an increase of
$Re$
leads to larger side (time-average) recirculating regions and, thus, to a decrease of the friction contribution to
$F_x$
(figure 25
b).
4.3.2. Intermediate
$Re$
: aperiodic and frequency-locking regimes

Figure 26. First aperiodic regimes for
$L=5$
and
$W=5$
. (a) Frequency spectra of
$F_y$
(left) and
$F_z$
(right) for
$305 \leqslant Re \leqslant 345$
. (b-c) Structure of the flow for
$Re=315$
(top) and
$Re=335$
(bottom). (b) Instantaneous isosurfaces
$\lambda _2=-0.25$
coloured by
$-1 \leqslant \omega _x \leqslant 1$
. (c) Force diagrams
$F_y-F_x$
(left) and
$F_z-F_x$
(right).
As the Reynolds number increases, the limit cycle loses its stability via a Neimark--Sacker bifurcation, and new frequencies appear in the flow. At
$Re \approx 305$
a torus replaces the limit cycle. The unsteady
$A_yS_z$
wake mode predicted by the LSA of the steady
$S_yS_z$
base flow (§ 3.2) becomes unstable, and the flow experiences an oscillating motion in both
$y$
and
$z$
directions. The new frequency is
$St \approx 0.03$
, in agreement with the value predicted by the LSA. As
$Re$
increases, the two modes interact nonlinearly and different frequencies arise in the spectra (figure 26), resulting in a progressive loss of coherency. However, the frequency spectra (and the POD, not shown) indicate that the horizontal oscillations at
$St=0.03$
are less intense than the vertical oscillations at
$St=0.09$
. As a result, up to
$Re \approx 325$
the dynamics is mainly driven by the vertical oscillation of the wake and the resulting hairpin-like structures.
For
$Re \gtrapprox 325$
, a new mode arises with a frequency in the range
$0.16 \leqslant St \leqslant 0.2$
associated with the shedding of spanwise vorticity tubes from the top/bottom LE shear layers that roll up and generate HVs in the wake, similarly to what was found for smaller
$W$
. The interaction of these vortices with the wake unsteadiness is analysed later in this subsection. The LE vortex shedding frequency for
$W=5$
is smaller than that for
$W=1.2$
and
$2.25$
. This LE vortex shedding progressively breaks the coherency of the flow (see the force diagrams and the emergence of several peaks in the frequency spectra in figure 26). As
$Re$
increases, the importance of the LE vortex shedding on the flow dynamics progressively increases (see the peaks at
$St \approx 0.09$
and
$0.19$
in
$S(F_z)$
and the following POD analysis) and, for
$Re \gtrapprox 330$
, its contribution is comparable to that of the wake oscillation.

Figure 27. Aperiodic and frequency-locking regimes for
$L=5$
and
$W=5$
at larger Reynolds numbers
$345 \leqslant Re \leqslant 385$
. Same conventions as figure 26, but with the main flow frequencies given in the bottom plot of panel (a). (b,c) Results are shown for
$Re=353$
(top) and
$Re=375$
(bottom). Note that vortices are shed from the LE shear layer in phase opposition for
$Re=353$
(lock-in region I) and in phase for
$Re=375$
(lock-in region II).
For two narrow ranges of Reynolds numbers, the flow oscillations influence each other in a way that produces synchronisation into a periodic regime. This phenomenon, known as frequency locking (Iooss & Joseph Reference Iooss and Joseph1990; Kuznetsov Reference Kuznetsov2004), has already been observed by Chiarini et al. (Reference Chiarini, Quadrio and Auteri2022d
) for LE and TE vortex shedding in the flow past 2-D rectangular cylinders. The two sheddings, of periods
$T_{LE}$
and
$T_{TE}$
, synchronise, and different stable cycles (periodic orbits) arise in the torus, with a long-time period
$T_{lp}=p T_{TE}=q T_{LE}$
, where
$p,q \in \mathbb {N}$
. In the present case, the frequency locking is visualised in figure 27. The characteristic frequencies of the vertical wake oscillation (
$St \approx 0.09-0.10$
; squares) and LE vortex shedding (
$St \approx 0.18-0.20$
; circles) increase weakly with
$Re$
. The frequency of the LE vortex shedding is approximately twice that of the wake oscillation, meaning that two LE vortices are shed from the top/bottom LE shear layers during one wake oscillation period. The flow is periodic for
$353 \leqslant Re \leqslant 357$
(lock-in region I) and
$ 373 \leqslant Re \leqslant 375$
(lock-in region II): the LE vortex shedding synchronises with the vertical wake oscillation, and a limit cycle arises with
$T_{lp}=T_{TE} = 2 T_{LE}$
(i.e.
$p=1$
and
$q=2$
). Notably, the relative phase between the top and bottom vortex shedding is different in the two synchronised regimes. The LE vortices are indeed shed from the top and bottom LE shear layers in phase opposition in region I and in phase in region II (see figure 27 and the following POD analysis). This explains why the
$St = 1/T_{LE}$
frequency is not visible in
$S(F_z)$
for the lock-in region II: since the top and bottom LE vortices are shed in phase, they do not create any instantaneous unbalance of vertical force. By contrast, in region I the closed trajectory in the
$F_z-F_x$
force diagram is not symmetric with respect to
$F_z=0$
, meaning that the time-averaged
$F_z$
is non-zero. More generally, the time-averaged top/bottom symmetry is lost for
$345 \leqslant Re \leqslant 360$
(figure 24). Regarding the lateral aerodynamic force, the different phase locking in the two synchronised regimes results in a different behaviour:
$F_y \approx 0$
at all times in region I, while
$F_y$
oscillates with
$St = 1/T_{TE} = 1/(2 T_{LE})$
in region II.

Figure 28. Frequency-locking regime for
$L=5$
,
$W=5$
and
$Re=353$
(lock-in region I). Lateral view of isosurfaces
$\lambda _2=-0.05$
coloured by streamwise vorticity (blue-to-red colour map for
$-1 \leqslant \omega _x \leqslant 1$
). The snapshots are separated in time by
$T/8$
, where
$T$
is the period of the wake oscillation. Black/red labels refer to vortices shed in the considered/previous period. Here LEV1, LEV2 and LEV3 refer to LE vortices shed from the top and bottom LE shear layers, respectively; HV1, HV2 and HV3 indicate HVs that arise in the wake once LEV1, LEV2 and LEV3 cross the TE.
We now illustrate the interaction between the LE and TE vortex shedding in region I. Figure 28 shows eight snapshots of the flow at
$Re=353$
, equispaced in one period
$T_{TE}$
. Black and red labels refer to vortices shed from the prism in the current and previous periods, respectively. Vortices shed from the top and bottom LE shear layers are in phase opposition. At
$t = t_1$
a spanwise-aligned vortex, LEV1, is generated from the top LE shear layer and convected downstream (see
$t=t_2$
). At
$t=t_3$
, LEV1 crosses the TE corner in phase with the upwards motion of the wake and rolls up (
$t_4$
to
$t_6$
), generating a large head-up hairpin vortex (HV1) that is then shed downstream (
$t_7$
and
$t_8$
). The same dynamics is observed over the bottom side. A spanwise-aligned vortex LEV3 is shed from the LE shear layer at
$t=t_3$
. It crosses the bottom TE at
$t=t_5$
and rolls up, generating a head-down hairpin vortex HV3 at
$t=t_6$
. Moving again to the top side, at
$t=t_5$
the second spanwise-aligned vortex LEV2 is shed from the top LE shear layer. It crosses the TE corner at
$t=t_5$
in phase with the downwards motion of the wake and generates a hairpin vortex HV2 that is then connected with the HV3 hairpin vortex.
In region II, LE vortices interact with the vertical wake oscillation in a similar way to generate the large HVs that characterise the wake dynamics. However, the top/bottom LE vortices are shed in phase, and thus, cross the top/bottom TE in phase (figure 27 b, bottom row). When a top/bottom LE vortex crosses the top/bottom TE in phase with an upward/downward motion of the wake, a head-up/head-down hairpin-like vortex arises in the wake.
In the aperiodic regimes (
$Re \in [345,350]$
,
$[360,365]$
and
$[380,385]$
), a similar interaction between the LE vortex shedding and the wake unsteadiness as in the lock-in region I is observed, although not fully periodic.

Figure 29. The POD analysis for
$L=5$
and
$W=5$
in the range
$305 \leqslant Re \leqslant 385$
. (a) Circles and squares refer to the unsteady wake mode (
$St \approx 0.09$
) and LE vortex shedding mode (
$St \approx 0.19$
), respectively. When the two phenomena are distributed over several POD modes, the sum of their energy content is reported in the right panel. (b) Lateral view of the POD mode associated with LE vortex shedding. Isosurfaces of
$\lambda _2$
coloured by streamwise vorticity. From left to right: aperiodic regime at
$Re=345$
, periodic regimes at
$Re=357$
and
$Re=375$
.
We use POD to further examine the dependence of the flow structure on
$Re$
(figure 29). In agreement with the above discussion, the dominant POD modes are associated with the wake unsteadiness and with the LE vortex shedding, and exhibit a single frequency of
$St \approx 0.09$
and
$St \approx 0.19$
, respectively. The POD modes associated with the wake unsteadiness mainly evolve in the wake and feature large-scale head-up and head-down HVs downstream of the TE (not shown for brevity). The POD modes associated with the LE vortex shedding, instead, originate over the top and bottom sides where they feature spanwise-aligned structures (figure 29
b). The spatial structure of the LE vortex shedding POD modes changes with
$Re$
, as visible along the sides of the prism. In the aperiodic regime (
$Re=345, 365, 385$
) and in lock-in region I (
$Re=357$
), the modes are antisymmetric with respect to an inversion of the
$z$
axis, with
$\omega _x(x,y,z) = \omega _x(x,y,-z)$
. This agrees with the out-of-phase shedding of vortices from the top and bottom LE shear layers. In lock-in region II (
$Re = 375$
), instead, the modes exhibit a top/down symmetry with
$\omega _x(x,y,z) = -\omega _x(x,y,-z)$
, which is consistent with in-phase LE vortex shedding. Also, the POD confirms that in the aperiodic
$aS_yS_z$
regime the relative intensity of the LE vortex shedding increases with
$Re$
(see right plot of figure 29
a). However, the wake unsteadiness dominates in the periodic regimes.
4.3.3. Large
$Re$
: aperiodic regime

Figure 30. Unsteady regimes for
$L=5$
and
$W=5$
at larger Reynolds numbers
$400 \leqslant Re \leqslant 600$
. (a) Frequency spectra of the aerodynamic forces
$F_y$
(left) and
$F_z$
(right). (b–c) Structure of the flow at
$Re=400$
(top) and
$Re=500$
(bottom). (b) Isosurfaces of
$\lambda_2 = -0.05$
coloured by streamwise vorticity. (c) Force diagrams
$F_y-F_x$
(left) and
$F_z-F_x$
(right).
For
$Re \geqslant 400$
, the synchronisation between the LE and TE vortex shedding is lost, and the flow becomes progressively more chaotic. Like for lower
$Re$
, the frequency of the wake oscillation remains approximately constant (
$St \approx 0.095{-}0.1$
for
$Re\in [400,600]$
), while that of the LE vortex shedding keeps increasing (
$St \approx 0.21{-}0.32$
). The increase in
$Re$
is accompanied by a slight increase in the oscillation amplitude in the
$y$
direction, as visualised in figure 24. For
$Re \geqslant 500$
, a third peak appears in
$S(F_y)$
at
$St \approx 0.125$
, and a shedding of HVs occurs also along the lateral sides of the prism (figure 30).
5. Conclusion

Figure 31. Regimes observed in the flow past rectangular prisms with
$L=5$
for Reynolds numbers up to
$Re=700$
. Blue, green and red shaded areas refer to steady, periodic and aperiodic states, respectively. Different tones refer to different regimes, accordingly to figures 10, 16 and 24. The solid lines separating the regimes are a guide to the eye. Vertical axis not to scale. Steady flows
$(a)$
-
$(c)$
are represented with
$\omega _x = \pm 0.1$
isosurfaces and unsteady flows
$(d)$
-
$(g)$
with snapshots of
$\lambda _2$
isosurfaces coloured by streamwise vorticity.
In this study we have investigated the sequence of bifurcations of the laminar flow past 3-D rectangular prisms for length-to-height and width-to-height ratios
$1 \leqslant L \leqslant 5$
and
$1.2 \leqslant W \leqslant 5$
, and Reynolds number up to
$Re \approx 700$
.
Linear stability analysis shows that the nature of the primary bifurcation changes with the geometry. The flow past wide prisms with large
$W/L$
experiences first a Hopf bifurcation and becomes unstable to oscillatory perturbations that break the top/bottom planar symmetry, leading to a periodic vortex shedding across the smaller dimension of the body. For smaller
$W/L$
, the primary instability consists of a pitchfork bifurcation: for intermediate (small)
$W/L$
, the flow becomes unstable to stationary perturbations that break the top/bottom (left/right) planar symmetry and lead to a static vertical (horizontal) deflection of the wake. In all cases, the critical
$Re$
of the primary instability increases with
$L$
and depends non-monotonically on
$W$
.
A WNL analysis has been performed for
$(L,W)=(5,1.2)$
and
$(L,W)=(5,2.25)$
in the vicinity of the critical
$Re$
of the static modes. Two coupled amplitude equations are obtained, revealing the sequence of bifurcations close to the bifurcation points. For
$W=1.2$
, third-order amplitude equations yield a supercritical bifurcation scenario similar to that of Ahmed bodies (Zampogna & Boujo Reference Zampogna and Boujo2023), with an exchange of stability between the vertically and horizontally deflected states at
$Re \approx 353$
. For
$W=2.25$
, fifth-order amplitude equations yield a subcritical bifurcation for the horizontally deflected state, with a narrow interval of bistability where the vertically and horizontally deflected states coexist. Fully nonlinear DNS confirm the bifurcation sequence predicted by the LSA and the WNL analysis, including the bistability-induced hysteresis for
$W=2.25$
and the bifurcated states.
At larger
$Re$
, nonlinear DNS revealed a rich sequence of bifurcations (figure 31) that has been investigated by means of frequency spectra, force diagrams and POD. Overall, the flow dynamics is driven by five different modes: (i,ii) static deflections of the wake in the vertical and horizontal directions, (iii) vortex shedding of HVs from the LE shear layers, and (iv,v) unsteady flapping of the wake in the vertical and horizontal directions. Their nonlinear interaction changes with
$W$
and
$Re$
giving rise to several regimes, ranging from steady and periodic regimes at small
$Re$
to aperiodic and chaotic regimes at larger
$Re$
. For intermediate
$W$
, we observed, for the first time, a periodic LE shedding of HVs that instantaneously preserves all spatial symmetries and generates neither lift nor side forces. Interestingly, in some portions of the parameter space the different modes synchronise and give origin to periodic regimes also at relatively large
$Re$
.
The mechanism sustaining the shedding of the HVs from the LE shear layers has been investigated. Unlike in the flow past 2-D rectangular cylinders, this LE vortex shedding does not require the presence of a sharp TE, and does not necessarily lock with the TE vortex shedding. Indeed, for
$W=2.25$
, the same vortex shedding has been found for a 3-D rectangular flat plate without TE (
$L = \infty$
). By rounding the LE corners, we have also shown that the LE vortex shedding is the result of a feedback mechanism embedded within the recirculating regions that form over the lateral sides of the prisms.
Having characterised the flow stability and dynamics over a wide range of body widths and lengths, the present study serves as a stepping stone to further investigation in more complex settings. For instance, the laminar wakes of spheres and 2-D cylinders translating close to a solid wall have been studied extensively (Thompson et al. Reference Thompson, Leweke and Hourigan2021). Ground proximity is by definition an essential feature for ground vehicles such as cars and trains, and although many numerical simulations have been carried out, the stability of these flows remains largely unexplored in both the laminar and turbulent regimes, except for a few studies (e.g. simplified train geometry with
$(L,W)=(6.9,0.83)H$
and ground clearance 0.15
$H$
in Li et al. Reference Li, Demange, Chen, Wang, Liang, Schmidt and Oberleithner2024). Another topic of interest is that of body attitude (yaw and pitch): in the laminar regime a small misalignment may significantly alter the onset and sequence of bifurcations; in the turbulent regime ground vehicles may be subject to crosswind or varying front/rear load distribution. Finally, the effect of stochastic conditions (e.g. disturbances in the incoming flow) on the laminar bifurcations and wake dynamics could be studied, especially near codimension-two bifurcations where two stationary modes become unstable simultaneously (e.g.
$(L,W)=(3,1.2)$
at
$Re \approx 330$
,
$(5,1.2)$
at
$Re \approx 350$
and
$(5,2.25)$
at
$Re \approx 300$
), which may result in multistability, with the wake switching randomly between two or more static deflected states. A treatment in the spirit of Ducimetière et al. (Reference Ducimetière, Boujo and Gallaire2024) would yield a rigorously derived reduced-order model in the form of coupled stochastic amplitude equations, one for each symmetry-breaking mode, which could then be used to predict flow statistics.
Supplementary material.
Supplementary material is available at https://doi.org/10.1017/jfm.2025.165.
Acknowledgements.
A.C. acknowledges the computer time provided by the Scientific Computing & Data Analysis section of the Core Facilities at OIST. E.B. thanks François Gallaire for useful discussions.
Funding.
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Convergence study
In this appendix we report the sensitivity of the results on the grid resolution for the LSA and 3-D nonlinear DNS.
For the LSA, tables 1
--2 show how the critical Reynolds numbers of the first two bifurcations vary with the mesh size and domain size. The numerical domain is
$\{ x,y,z \, | \, x_{min} \leqslant x \leqslant x_{max}; \, 0 \leqslant y,z \leqslant y_{max}=z_{max} \}$
. The mesh density is
$n_1$
on the prism surface,
$n_2$
on the boundaries of the subdomain
$\{ x,y,z \, | \, -5 \leqslant x \leqslant 15; \, 0 \leqslant y,z \leqslant 2 \}$
and
$n_3=1$
on the outermost boundaries. Given the weak influence of the mesh size and domain size on
$Re_{c,1}$
and
$Re_{c,2}$
, we choose mesh M3 and domain D1 throughout the linear and WNL analyses, similar to Zampogna & Boujo (Reference Zampogna and Boujo2023).
Table 1. Influence of the mesh size on the first two critical Reynolds numbers for four different prism geometries. Domain D1:
$(x_{min},x_{max})=(-10,20)$
,
$y_{max}=z_{max}=10$
. Here
$N_{elmts}$
is in millions.

Table 2. Influence of the domain size on the first two critical Reynolds numbers for four different prism geometries. Mesh M3:
$(n_1,n_2)=(60,10)$
. Here
$N_{elmts}$
is in millions.

For the 3-D nonlinear DNS, sensitivity to the grid resolution and to the time step has been investigated with six additional simulations, two for each considered geometry. The number of points in the streamwise and vertical directions has been decreased from
$(N_x,N_z)=(1072,590)$
on the standard grids to
$(N_x,N_z)=(860,400)$
for all three geometries; in the spanwise direction, the number of points has been decreased from
$N_y=666,720$
and
$804$
to
$N_y=460,500$
and
$560$
for
$W=1.2,2.25$
and
$5$
respectively. On the coarser grid, two different time steps have been used,
$\Delta t \approx 0.0066$
and
$0.0033$
, leading to an average CFL number of approximately
$1$
and
$0.5$
. For each geometry, we considered a single Reynolds number:
$Re = 535$
for
$W=1.2$
(
$aA_yS_z$
regime),
$Re=450$
for
$W=2.25$
(
$pS_yS_zlb$
) and
$Re=385$
for
$W=5$
(
$aS_yS_z$
). In all cases, we find an excellent qualitative and quantitative agreement on the mean and fluctuating forces between the standard and coarse grids, and between the two
$\Delta t$
as well (table 3). We observe an excellent agreement in terms of frequencies too, with differences smaller than 0.7 % between the coarse and standard grids (not shown).
Table 3. Convergence study for the 3-D nonlinear simulations. Influence of the grid resolution on the aerodynamic forces for
$(L,W)=(5,1.2)$
and
$Re=535$
(
$aA_yS_z$
regime),
$(L,W)=(5,2.25)$
and
$Re=450$
(
$pS_yS_zlb$
regime), and
$(L,W)=(5,5)$
and
$Re=385$
(
$aS_yS_z$
regime).

Appendix B. Weakly nonlinear analysis
We detail here the derivation of the systems of third- and fifth-order amplitude equations (3.1), (3.2) and (3.3), (3.4) in the vicinity of a codimension-two bifurcation for two stationary modes
$\rm A$
and
$\rm B$
. The procedure is similar to that in Zampogna & Boujo (Reference Zampogna and Boujo2023).
B.1. Derivation of the amplitude equations
When the two modes of interest do not bifurcate exactly at the same Reynolds number, we introduce a reference critical Reynolds number
$Re_c$
, typically chosen between the critical Reynolds numbers of the two modes. Departure from criticality is measured as

which defines the small parameter
$0\lt \epsilon \ll 1$
and the order-one parameter
$\tilde \alpha =O(1)$
. We also define a shift operator
$\mathcal {S}$
such that, unlike the original linearised NS operator
$\mathcal {A}$
, the shifted linearised NS operator
$\widetilde {\mathcal {A}} = \mathcal {A} + \epsilon ^2 \mathcal {S}$
is singular exactly at
$Re=Re_c$
for both modes. In other words, denoting
$\mathcal {B} \partial _t \boldsymbol {q} + {\mathcal {A}} \boldsymbol {q} = \boldsymbol {0}$
, the linearised NS equations and
$\hat {\boldsymbol {q}} = \{\hat {\boldsymbol {u}},\hat {p}\}$
, one has
$\widetilde {\mathcal {A}} \hat {\boldsymbol {q}}_1^A = \widetilde {\mathcal {A}} \hat {\boldsymbol {q}}_1^B = \boldsymbol {0}$
, while
$\mathcal {A} \hat {\boldsymbol {q}}_1^A = -\sigma _A \mathcal {B} \hat {\boldsymbol {q}}_1^A$
and
$\mathcal {A} \hat {\boldsymbol {q}}_1^B = -\sigma _B \mathcal {B} \hat {\boldsymbol {q}}_1^B$
.
We use the method of multiple scales. We introduce slow time scales
$T_1 = \epsilon ^2 t$
,
$\quad T_2 = \epsilon ^4 t, \ldots$
and inject the expansion

in the NS equations at
$Re=Re_c$
, where the time derivative now reads
$\partial _t + \epsilon ^2 \partial _{T_1} + \epsilon ^4 \partial _{T_2} + \ldots$
. Collecting like-order terms yields the following series of problems.
B.1.1. Zeroth and first orders
At order
$\epsilon ^0$
we find the nonlinear NS equations and the zeroth-order field
$\boldsymbol {q}_0(\boldsymbol {x})$
is the steady
$S_yS_z$
base flow at
$Re=Re_c$
. At order
$\epsilon ^1$
we obtain the linearised and shifted NS equations

Since
$\widetilde {\mathcal {A}}$
is singular at
$Re=Re_c$
, the first-order field is a superposition of the
$S_yA_z$
and
$A_yS_z$
eigenmodes:

Here
$A(T_1,T_2,\ldots )$
and
$B(T_1,T_2,\ldots )$
are real-valued slowly varying amplitudes to be determined.
B.1.2. Second order
At order
$\epsilon ^2$
the field
$\boldsymbol {q}_2$
is a solution of the linearised and shifted NS equations

where

and introducing the notation
$\mathcal {C}( \boldsymbol a , \boldsymbol b) = (\boldsymbol a \cdot \nabla ) \boldsymbol b + (\boldsymbol b \cdot \nabla ) \boldsymbol a.$
Symmetry considerations show that (B5) can be inverted: all the terms in
$\boldsymbol F_2$
are either
$S_y S_z$
or
$A_y A_z$
and cannot resonate since
$\widetilde {\mathcal {A}}$
is singular to
$S_y A_z$
and
$A_y S_z$
modes. Therefore,

where the
$\boldsymbol q_2^*$
fields are solutions of

B.1.3. Third order
At order
$\epsilon ^3$
the field
$\boldsymbol {q}_3$
is a solution of the linearised and shifted NS equations

where

and where
$\tilde \sigma _A=\sigma _A/\epsilon ^2$
and
$\tilde \sigma _B=\sigma _B/\epsilon ^2$
. All the terms in
$\boldsymbol F_3$
are resonant with either
$\hat {\boldsymbol u}_1^A$
or
$\hat {\boldsymbol u}_1^B$
. To avoid secular terms, we impose a compatibility condition, i.e. we require the respective inner products with either
$\hat {\boldsymbol u}_1^{A\dagger }$
or
$\hat {\boldsymbol u}_1^{B\dagger }$
be zero. With the normalisation
$\langle \hat {\boldsymbol u}_1^{A\dagger }, \hat {\boldsymbol u}_1^A \rangle = \langle \hat {\boldsymbol u}_1^{B\dagger }, \hat {\boldsymbol u}_1^B \rangle = 1$
, we obtain


with the coefficients

and similar expressions for
$\tilde \lambda _B$
,
$\tilde \chi _B$
and
$\tilde \eta _B$
upon exchange of indices
$A$
and
$B$
.
Equations (B11) and (B12) constitute the leading-order system of amplitude equations. If the derivation is stopped at this order, one can reintroduce the fast time
$t=\epsilon ^2 T_1$
and obtain the system (3.1), (3.2), with the coefficients
$\lambda _A = \epsilon ^2 \tilde \lambda _A$
,
$\chi _A = \epsilon ^2 \tilde \chi _A$
, etc.:


To derive higher-order amplitude equations, the third-order field

is needed, where the
$\boldsymbol q_3^*$
fields are solutions of

with similar expressions for
$\boldsymbol q_3^{B}$
,
$\boldsymbol q_3^{B^3}$
and
$\boldsymbol q_3^{A B^2}$
upon exchange of
$A$
and
$B$
.
Because the kernel of
$\widetilde {\mathcal {A}}$
is not empty, the
$\boldsymbol q_3^*$
fields are defined up to an arbitrary constant along
$\hat {\boldsymbol u}_1^A$
or
$\hat {\boldsymbol u}_1^B$
. We remove this component such that
$\langle \hat {\boldsymbol u}_1^{A\dagger }, \boldsymbol u_3^A \rangle =\langle \hat {\boldsymbol u}_1^{A\dagger }, \boldsymbol u_3^{A^3} \rangle =\langle \hat {\boldsymbol u}_1^{A\dagger }, \boldsymbol u_3^{AB^2} \rangle = \langle \hat {\boldsymbol u}_1^{B\dagger }, \boldsymbol u_3^B \rangle =\langle \hat {\boldsymbol u}_1^{B\dagger }, \boldsymbol u_3^{B^3} \rangle =\langle \hat {\boldsymbol u}_1^{B\dagger }, \boldsymbol u_3^{A^2B} \rangle =0$
.
B.1.4. Fourth order
At order
$\epsilon ^4$
the field
$\boldsymbol {q}_4$
is a solution of the linearised and shifted NS equations

where symmetry considerations show that
$\boldsymbol F_4$
is not resonant. Therefore, one can invert (B18) and compute

where the
$\boldsymbol q_4^*$
fields are solutions of

with similar expressions for
$\boldsymbol q_4^{(B^2)'}$
,
$\boldsymbol q_4^{B^2}$
,
$\boldsymbol q_4^{A B^3}$
,
$\boldsymbol q_4^{B^4}$
upon exchange of
$A$
and
$B$
.
B.1.5. Fifth order
At order
$\epsilon ^5$
the field
$\boldsymbol {q}_5$
is a solution of the linearised and shifted NS equations

where all the terms in
$\boldsymbol F_5$
resonate with either
$\hat {\boldsymbol u}_1^A$
or
$\hat {\boldsymbol u}_1^B$
. Imposing a compatibility condition yields

where

A similar equation is obtained for
$\partial _{T_2}B$
upon exchange of
$A$
and
$B$
.
Terms proportional to
$\tilde \lambda _A$
,
$\tilde \chi _A$
, etc., come from the interaction of the eigenmodes with the fourth-order fields
$ {\boldsymbol q_4^{(A^2)'}}$
,
$\boldsymbol q_4^{(AB)'}$
and
$\boldsymbol q_4^{(B^2)'}$
, and are obtained after developing
$ \partial _{T_1}(A^2)$
,
$\partial _{T_1}(AB)$
and
$\partial _{T_1}(B^2)$
and using (B11) and (B12). By contrast, all the forcing terms of the form
$ \partial _{T_1}(A^m B^n) {\boldsymbol u_3^{*}}$
with
$m+n=1$
or 3 vanish when imposing the compatibility condition because of our earlier choice
$\langle \hat {\boldsymbol u}_1^{\dagger }, \boldsymbol u_3^* \rangle =0$
.
Defining coefficients
$\tilde *'$
, we obtain


Finally, the total fifth-order amplitude equations are obtained by reconstructing the total derivatives
$\mathrm {d}_t(A,B) = (\epsilon ^2 \partial _{T_1}+ \epsilon ^4 \partial _{T_2})(A,B)$
, i.e. combining (B11), (B12) and (B24), (B25):


Here now
$\lambda _A = \epsilon ^2\tilde \lambda _A + \epsilon ^4\tilde \lambda _A'$
,
$\chi _A = \epsilon ^2\tilde \chi _A + \epsilon ^4\tilde \chi _A'$
, etc. Therefore, compared with the third-order system, the fifth-order system has not only higher-order terms but also corrected coefficients.
B.2. Equilibrium solutions
B.2.1. Third-order system
The third-order system (3.1), (3.2) has the following four possible equilibrium solutions.
-
• Symmetric state:
$ A=B=0$ .
-
• Pure vertical symmetry breaking:
$A^2 = \lambda _A/\chi _A$ ,
$B=0$ .
-
• Pure horizontal symmetry breaking:
$A=0$ ,
$B^2 = \lambda _B/\chi _B$ .
-
• Mixed state (double symmetry breaking):
$A^2 = (\chi _B \lambda _A - \eta _A \lambda _B)/(\chi _A \chi _B -\eta _A\eta _B)$ ,
$B^2 = (\chi _A \lambda _B - \eta _B \lambda _A)/(\chi _A \chi _B -\eta _A\eta _B).$
B.2.2. Fifth-order system
The fifth-order system (3.3), (3.4) has four types of equilibrium states, with the following multiple possible amplitudes.
-
• Symmetric state:
$A=B=0$ .
-
• Pure vertical symmetry breaking:
$A^2 = \pm \sqrt { (\chi _A \pm \sqrt {\Delta })/(2\gamma _A) }$ ,
$B=0,$ where
$\Delta = \chi _A^2 - 4\lambda _A\gamma _A$ . Real-valued solutions exist when
$\Delta \gt 0$ and
$(\chi _A \pm \sqrt {\Delta })/\gamma _A\gt 0$ . Focusing on positive values, this yields up to two different amplitudes
$A_1$ ,
$A_2$ .
-
• Pure horizontal symmetry breaking:
$A=0$ ,
$B^2 = \pm \sqrt { (\chi _B \pm \sqrt {\Delta })/(2\gamma _B) },$ where
$\Delta = \chi _B^2 - 4\lambda _B\gamma _B$ . Real-valued solutions exist when
$\Delta \gt 0$ and
$(\chi _B \pm \sqrt {\Delta })/\gamma _B\gt 0$ . Again, there are up to two different positive amplitudes
$B_1$ ,
$B_2$ .
-
• Mixed states (double symmetry breaking)
$A\neq 0, B\neq 0$ . These states correspond to the intersections of conic sections, since
$u=A^2$ and
$v=B^2$ satisfy the bivariate quadratic equations
(B28)\begin{align} 0 &= \lambda _A - \chi _A u - \eta _A v + \kappa _A v^2 + \beta _A u v + \gamma _A u^2, \end{align}
(B29)\begin{align} 0 &= \lambda _B -\chi _B v -\eta _B u + \kappa _B u^2 + \beta _B u v + \gamma _B v^2. \end{align}
An analytic treatment of these equations can be found, for example, in chapter 11 of Richter-Gebert (Reference Richter-Gebert2011). There are in general zero, two or four real-valued intersections. Here, we compute these intersections numerically with the algorithm of Nakatsukasa et al. (Reference Nakatsukasa, Noferini and Townsend2015).

Figure 32. The POD analysis for
$L=5$
and
$W=2.25$
at
$Re=420$
(left) and
$Re=450$
(right). In both cases the flow is periodic and has a single dominant mode, of frequency
$St \approx 0.255$
and
$St \approx 0.265$
, respectively (see figure 19). (a) Energy fractions of the first eight modes and corresponding frequencies. (b) Structure of POD mode 1, to be compared with figure 20. Isosurfaces of
$\lambda _2$
coloured by streamwise vorticity
$\omega _x$
(blue-to-red colour map ranges from negative to positive values). For
$Re=420$
, the structure of the mode is the same as for
$Re=380$
, and shows the same symmetries. For
$Re=450$
, instead, the symmetries of the mode change.
Appendix C. Proper orthogonal decomposition
C.1. Numerical implementation
To perform the POD analysis of § 4, we used the method of snapshots (Sirovich Reference Sirovich1987) on a portion of the domain (
$-3.7 \leqslant x \leqslant 17.7$
,
$|y|,|z| \leqslant 3.2$
for
$W=1.2,2.25$
and
$-3.7 \leqslant x \leqslant 17.7$
,
$|y| \leqslant 4.2$
,
$|z| \leqslant 6.2$
for
$W=5$
) and, to avoid the need of weight matrices for scalar products, interpolated on a uniform cubic grid (with
$N_x=410$
and
$N_y=N_z=120$
for
$W=1.2,2.25$
and
$N_x=318$
,
$N_y=100$
and
$N_z=150$
for
$W=5$
). This yields approximately
$5.9 \times 10^6$
(for
$W=1.2,2.25$
) and
$4.8 \times 10^6$
(for
$W=5$
) grids points and, with three velocity components at each point,
$N_{dof} \approx 17.7 \times 10^6$
and
$N_{dof} \approx 14.4 \times 10^6$
degrees of freedom. The number of snapshots and their time separation has been chosen carefully: the sampling time
$\delta t_s \simeq 1/(20 f_l)$
ensures that the largest frequency
$f_l$
of interest detected in the spectra is captured; similarly, the total time
$T \simeq 8/f_s$
ensures that the smallest frequency
$f_s$
is captured. For example, for
$W=2.25$
and
$Re=380$
, this yields
$\delta t_s = 0.195$
and
$T= 600 \delta t_s$
(i.e.
$m=600$
snapshots).
With the snapshots collected in the
$N_{dof}\times m$
matrix
$S$
, the
$m\times m$
covariance matrix
$S^T S$
is formed and the eigenvalue problem
$S^T S Z = Z \Lambda$
is solved for the
$r$
largest eigenvalues
$\lambda _j$
(components of the diagonal matrix
$\Lambda$
) and corresponding eigenvectors
$\boldsymbol {z}_j$
(columns of
$Z$
). The positive square roots of the eigenvalues of
$S^TS$
are the singular values of
$S$
and the columns of the matrix
$\Phi = SZ$
are the POD modes. Frequencies associated with a given POD mode are detected by projecting the snapshots on that mode and performing a frequency analysis of the so-obtained time-dependent coefficient.
C.2. Effect of
$Re$
on the LE vortex shedding POD mode for
$W=2.25$
Figure 32 considers
$L=5$
and
$W=2.25$
and shows the effect of
$Re$
on the spatial structure of the POD mode associated with the LE vortex shedding.
The structure of that mode does not change when
$Re$
is increased from
$380$
to
$420$
. The flow is periodic and the unsteadiness is driven by the in-phase shedding of HVs from the top/bottom LE shear layers (figure 19). When
$Re$
is further increased, however, the symmetries of the mode change, in agreement with a shedding of HVs in phase opposition from the top /bottom LE shear layers: for
$380 \leqslant Re \leqslant 420$
, the POD mode is symmetric with respect to
$z=0$
, with
$\omega _x(x,y,z) = - \omega _x(x,y,-z)$
. whereas for
$440 \leqslant Re \leqslant 460$
, it is antisymmetric, with
$\omega _x(x,y,z) = \omega _x(x,y,-z)$
.