4.1. Effect of $Re$

We start by illustrating the effect of the Reynolds number for a selected geometry, namely $L=3$ and $R=0$. Figures 3 and 4 show cuts of the streamwise velocity field and contours of zero streamwise velocity for different values of $Re$. There are several recirculation regions found around the prism: a ‘wake recirculation’ originating from the trailing edges, and ‘side recirculations’ originating from the leading edges. Figure 5(a) shows the variation in length of these regions: the wake recirculation length $l_w$, measured along the symmetry axis $y=z=0$, and the side recirculation lengths $l_s$, measured in the symmetry planes on the left/right faces ($|y|=W/2$, $z=0$) and on the upper/lower faces ($y=0$, $|z|=H/2$). Since $W>H$, side recirculations on the upper/lower faces are slightly longer than those on the left/right faces. Their length is seen to increase monotonically with $Re$ in the investigated range of Reynolds number, until they extend all the way down to the trailing edges ($l_s = L = 3$) and connect with the wake recirculation at $Re \simeq 500$. By contrast, the length of the wake recirculation varies non-monotonically, reaching a maximum at $Re \simeq 350$. The maximum backflow, i.e. the minimum streamwise velocity along the symmetry axis, $U_b = \min u_0(x,0,0)$, decreases in the investigated range of $Re$ from $-$0.18 to $-$0.26 (inset of figure 5a). Its non-monotonic variation seems correlated with that of $l_w$. In addition, figure 4(a) shows that although the width of the backflow region in the wake increases monotonically, its height varies non-monotonically. As a result, vertical cross-sections of the backflow region change from circles to wider-than-tall ellipses. It is important to note that 2-D cuts in the symmetry planes $y=0$ and $z=0$ provide only partial information about the 3-D structure of the flow. The 3-D views in figures 4(b,c) show, for instance, that reattachment occurs farther downstream in the symmetry planes than away from them, and that the backflow regions on the sides of the prism become non-convex for large enough $Re$. All this suggests that even before they merge at $Re \simeq 500$, the recirculation regions on the sides and in the wake undergo some 3-D reorganisation, which affects $l_w$ and $U_b$.

Figure 3.

Streamwise velocity $u_0$ of the base flow past a rectangular prism, $W=1.2$, $L=3$, $R=0$, in the planes $z=0$ (top view), $y=0$ (side view) and $x=2.5$ (rear view), for (a) $Re=150$, (b) $Re=250$, (c) $Re=350$, and (d) $Re=450$.

Figure 4.

Contours of zero streamwise velocity $u_0=0$ for the base flow past a rectangular prism, $W=1.2$, $L=3$, $R=0$: (a) in the planes $z=0$ (top view $a1$), $y=0$ (side view $a2$) and $x=2.5$ (rear view $a3$), for $Re=150$, 250, 350 and 450; and 3-D views for (b) $Re=250$ and (c) $Re=450$. While the width of the wake recirculation increases monotonically with $Re$ (plot $a1$), its length and height vary non-monotonically (plot $a2$).

Figure 5.

Base flow past a rectangular prism, $W=1.2$, $L=3$, $R=0$. (a) Length of the wake recirculation and side recirculations. Dashed line is $l_s=L=3$, when the side recirculations reach the end of the body and connect with the wake recirculation. Inset: maximum backflow $U_b$. (b) Drag coefficient.

At larger $Re$, the wake recirculation becomes substantially longer, wider and taller (not shown) as it merges with the side recirculations, and the backflow becomes stronger. We do not, however, focus on this topological change of the base flow because bifurcations of interest occur at smaller $Re$, as will become clear in §§ 5 and 6.

Figure 5(b) shows the evolution of the drag coefficient,

(4.3)\begin{equation} C_x = \frac{F_x}{\tfrac{1}{2} \rho U_\infty^2 W H},\quad \mbox{where } F_x ={-} \oint_{\varGamma_{b}} \left( \boldsymbol{\sigma}_0 \boldsymbol{\cdot} \boldsymbol{n} \right) \boldsymbol{\cdot} \boldsymbol{e}_x \,\mathrm{d}\varGamma, \end{equation}

where $\varGamma _b$ is the surface of the full body, and the stress tensor $\boldsymbol {\sigma }_0 = -p_0 \boldsymbol {I} + Re^{-1}\,\boldsymbol {\nabla } \boldsymbol u_0$ includes pressure and viscous effects. The drag coefficient decreases with $Re$, as do both pressure and viscous contributions, which is typically of bluff bodies in the laminar regime (see e.g. Henderson (Reference Henderson1995) for the 2-D circular cylinder wake). As $Re$ increases, the viscous contribution becomes much smaller than its pressure counterpart, because the side recirculations extend over a larger body surface area while the wake recirculation is not modified substantially.

4.2. Effect of $L$

We now investigate the effect of the length $L$ on the steady base flow past bodies with sharp edges ($R=0$). For the sake of illustration, we choose $Re=250$ (which will prove to be close to the critical Reynolds number in § 5). Figure 6 shows the streamwise velocity field for $L$ between $1/6$ and 3. For large $L$, the flow separates at the leading edges and reattaches onto the body, as described previously; the flow then re-separates at the trailing edges with a small angle with respect to the streamwise direction, which leads to a rather narrow and short wake recirculation. As $L$ decreases, the side recirculations and the wake merge, and the wake recirculation becomes longer and the backflow stronger. The overall flow also becomes increasingly different in the two symmetry planes $y=0$ and $z=0$: separation occurs with a larger angle from the upper/lower leading edges, and the wake becomes significantly taller, but not wider.

Figure 6.

Streamwise velocity $u_0$ of the base flow past rectangular prisms, $W=1.2$, $R=0$, at $Re=250$, in the planes $z=0$ (top view), $y=0$ (side view) and $x=2.5$ (rear view), for (a) $L=1/6$, (b) $L=0.5$, (c) $L=1$, (d) $L=1.5$, (e) $L=2$, and (f) $L=2.5$.

These qualitative observations are confirmed in figure 7(a): the length of the wake recirculation decreases with $L$, while the side recirculations (disconnected from the wake only when $L > l_w$) keep a fairly constant length. The backflow $U_b$ (inset) is much stronger when $L\lesssim 1.5$, i.e. when the side and wake recirculations are connected.

Figure 7.

Base flow past rectangular prisms, $W=1.2$, $R=0$, at $Re=250$. (a) Length of the wake recirculation and side recirculations. Dashed line is $l_s=L$, when the side recirculations reach the end of the body and connect with the wake recirculation. Inset: maximum backflow $U_b$. (b) Drag coefficient.

Figure 7(b) shows that the drag coefficient varies non-monotonically with $L$, reaching a minimum for $L \simeq 1.5$. This is the result of a competition between pressure and viscous effects. Pressure drag decreases with $L$, as longer bodies have narrower wakes, associated with a weaker front/rear pressure difference. Conversely, viscous drag increases with $L$, as longer bodies are subjected to positive wall shear stress over a wider surface area.

4.3. Effect of $R$

We now investigate the effect of the fillet radius $R$ on the steady base flow past bodies of lengths $L=1$ and $L=3$. Figure 8 shows the streamwise velocity field for $L=3$ and several fillet radius values. Rounding the front edges clearly suppresses the side recirculations. It turns out that a very small fillet, $R \gtrsim 0.05$, is sufficient to keep the flow attached to the body. For the shorter body, $L=1$, this also makes the wake recirculation narrower and slightly shorter. For the longer body, $L=3$, this has no visible effect on the wake recirculation. Figure 9 confirms that the wake recirculation length is slightly reduced for $L=1$ and barely affected for $L=3$. The backflow (inset) follows different trends for $L=1$ and $L=3$.

Figure 8.

Streamwise velocity $u_0$ of the base flow past rectangular prisms, $W=1.2$, $L=3$, at $Re=250$, in the planes $z=0$ (top view), $y=0$ (side view) and $x=2.5$ (rear view), for (a) $R=0$, (b) $R=0.05$, (c) $R=0.1$, (d) $R=0.2$, (e) $R=0.3472$, and (f) $R=0.5$.

Figure 9.

Length of the wake recirculation in the flow past rectangular prisms, $W=1.2$, at $Re=250$, for (a) $L=1$, and (b) $L=3$. Insets: maximum backflow.

Finally, figure 10 shows that $R$ has a strong effect on drag. For both $L=1$ and $L=3$, pressure drag decreases with $R$, as a result of the reduced front area. By contrast, viscous drag increases because boundary layers stay attached to the body. Overall, pressure effects are stronger and $C_x$ decreases.

Figure 10.

Drag coefficient of rectangular prisms, $W=1.2$, at $Re=250$, for (a) $L=1$, and (b) $L=3$.

6.1. Derivation of the amplitude equations

For the reference Ahmed body, figure 16 shows the eigenspectrum at $Re=300$, as well as the evolution of the growth rates of modes $A$ and $B$ for Reynolds numbers in the vicinity of $Re=300$. Mode $A$ becomes unstable slightly before mode $B$: $Re_c^A = 293$ and $Re_c^B = 304$. These two modes are depicted in figures 17(a,b), together with the associated adjoint modes in figures 17(c,d). As already observed, mode $A$ breaks the top/down symmetry, while mode $B$ breaks the left/right symmetry. Velocity perturbations extend far downstream and, as expected for stationary modes, do not show the wavepacket structures typical of oscillatory modes. Adjoint perturbations have the same symmetries and are localised in the recirculation region. We note that although eigenmodes are generally complex-valued, stationary modes like modes $A$ and $B$ are associated with real eigenvalues ($\lambda =\sigma$) and therefore can always be defined as purely real-valued. Figures 17(e,f) show the structural sensitivities of modes $A$ and $B$, which are very similar to those of the two stationary modes in the wake of the shorter prism $L=1$ without fillet (figure 12). Therefore, for the reference Ahmed body too, modes $A$ and $B$ are mostly sensitive to localised perturbations in the recirculation region.

Figure 16.

(a) Eigenvalue spectrum for the reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$) at $Re=300$, just between the first and second pitchfork bifurcations ($Re_c^A = 293$, $Re_c^B = 304$). (b) Symbols indicate growth rates of modes $A$ and $B$, the two leading stationary modes. Dashed lines indicate shifted growth rates $\sigma _A-\sigma _A(Re_c)$ and $\sigma _B-\sigma _B(Re_c)$ for $Re_c=300$ (see the WNL analysis in § 6).

Figure 17.

Direct modes, adjoint modes and structural sensitivities for the reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$) at $Re=300$, in the planes $z=0$ (top view), $y=0$ (side view) and $x=2.5$ (rear view): (a) $\hat {\boldsymbol {u}}_1^A$, (b) $\hat {\boldsymbol {u}}_1^B$, (c) $\hat {\boldsymbol {u}}_1^{A{\dagger} }$, (d) $\hat {\boldsymbol {u}}_1^{B{\dagger} }$, (e) $\|\hat {\boldsymbol {u}}_1^{A{\dagger} }\| \, \|\hat {\boldsymbol {u}}_1^{A}\|$, and (f) $\|\hat {\boldsymbol {u}}_1^{B{\dagger} }\| \, \|\hat {\boldsymbol {u}}_1^{B}\|$.

We assume that the Reynolds number is close to the bifurcation thresholds of modes $A$ and $B$. We introduce a reference critical Reynolds $Re_c$ that we choose typically in $[Re_c^A, Re_c^B]$, for example $Re_c=300$ for the reference Ahmed body. In practice, the specific choice of $Re_c$ has little effect, as shown in Appendix B. We quantify the departure from criticality with

(6.1)\begin{equation} \frac{1}{Re_c} - \frac{1}{Re} = \alpha = \epsilon^2 \tilde\alpha, \end{equation}

where $0 <\epsilon \ll 1$, and $\tilde \alpha$ is a parameter of order 1 (negative for $Re< Re_c$ and positive for $Re>Re_c$). At $Re_c$, the growth rates $\sigma _A$ and $\sigma _B$ are small but non-zero: mode $A$ is slightly unstable, and mode $B$ is slightly stable (figure 16). Following Meliga et al. (Reference Meliga, Chomaz and Sipp2009a), we introduce rescaled order 1 growth rates

(6.2a,b)\begin{equation} \tilde\sigma_A = \frac{\sigma_A}{\epsilon^2},\quad \tilde\sigma_B = \frac{\sigma_B}{\epsilon^2}, \end{equation}

together with a shift operator $\mathcal {S}$ such that

(6.3a,b)\begin{equation} \mathcal{S} \hat{\boldsymbol q}_1^A = \tilde \sigma_A(Re_c)\,\mathcal{B} \hat{\boldsymbol q}_1^A,\quad \mathcal{S} \hat{\boldsymbol q}_1^B = \tilde \sigma_B(Re_c)\,\mathcal{B} \hat{\boldsymbol q}_1^B, \end{equation}

and $\mathcal {S} \hat {\boldsymbol q}_1 = \boldsymbol 0$ for all other modes. The interest of the shift operator is that although the linearised NS operator is not singular at $Re=Re_c$ (since $\mathcal {A} \hat {\boldsymbol q}_1^A = -\sigma _A \mathcal {B} \hat {\boldsymbol q}_1^A \neq \boldsymbol 0$ and $\mathcal {A} \hat {\boldsymbol q}_1^B = -\sigma _B \mathcal {B} \hat {\boldsymbol q}_1^B \neq \boldsymbol 0$), the shifted linearised NS operator defined by

(6.4)\begin{equation} \widetilde{\mathcal{A}} = \mathcal{A} + \epsilon^2 \mathcal{S} \end{equation}

is exactly singular at $Re=Re_c$ for both modes $A$ and $B$, i.e.

(6.5)\begin{equation} \widetilde{\mathcal{A}} \hat{\boldsymbol q}_1^A =\widetilde{\mathcal{A}} \hat{\boldsymbol q}_1^B = \boldsymbol 0. \end{equation}

In other words, modes $A$ and $B$ are exactly neutral for $\widetilde {\mathcal {A}}$ at $Re=Re_c$, as shown by the shifted growth rates $\sigma _A-\sigma _A(Re_c)$ and $\sigma _B-\sigma _B(Re_c)$ in figure 16.

We perform our WNL analysis with the method of multiple scales, therefore we introduce the slow time scale $T = \epsilon ^2 t$. The flow field expansion

(6.6)\begin{equation} \boldsymbol q(\boldsymbol x,t,T) = \boldsymbol q_0 + \epsilon \boldsymbol q_1+ \epsilon^2 \boldsymbol q_2 + \epsilon^3 \boldsymbol q_3 + \cdots \end{equation}

is injected in the NS equations (5.2) at $Re=Re_c$, where $\partial _t$ is now transformed into $\partial _t + \epsilon ^2 \partial _T$. Collecting like-order terms in $\epsilon$ and making use of (6.4) yields a series of problems detailed hereafter.

6.1.1. Orders $\epsilon ^0$ and $\epsilon ^1$

The equations at order $\epsilon ^0$ are the nonlinear NS equations (2.1a,b), and the zeroth-order field $\boldsymbol q_0(\boldsymbol x)$ is the steady base flow computed in § 4 at the reference Reynolds number $Re_c$.

The equations at order $\epsilon ^1$ are the linearised NS equations (5.3) at $Re=Re_c$, now with the shifted linearised NS operator:

(6.7)\begin{equation} \mathcal{B}\,\partial_t {\boldsymbol q}_1 + \widetilde{\mathcal{A}} {\boldsymbol q}_1 = \boldsymbol{0}. \end{equation}

We assume that the first-order field $\boldsymbol q_1$ is a superposition of the two stationary modes $A$ and $B$ computed in § 5 and shown in figure 17:

(6.8)\begin{equation} {\boldsymbol q}_1(\boldsymbol x,T) = A(T)\,\hat{\boldsymbol q}_1^A + B(T)\,\hat{\boldsymbol q}_1^B, \end{equation}

where $A$ and $B$ are two real-valued slowly varying amplitudes, yet to be determined.

6.1.2. Order $\epsilon ^2$

At order $\epsilon ^2$, the second-order field $\boldsymbol q_2$ is a solution of the linearised NS equations

(6.9)\begin{equation} \mathcal{B}\,\partial_t {\boldsymbol q}_2 + \widetilde{\mathcal{A}} {\boldsymbol q}_2 = (\boldsymbol F_2,0)^{\rm T}, \end{equation}

forced by a term that depends on lower-order fields only,

(6.10)\begin{equation} \boldsymbol F_2={-} \tilde\alpha\,\nabla^2 \boldsymbol{u}_0 -(\boldsymbol{u}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u}_1.\end{equation}

The first term in (6.10) is due to Reynolds number variations (proportional to $\tilde \alpha$) in the base flow $\boldsymbol {u}_0$, and the second term is due to the transport of $\boldsymbol {u}_1$ by itself. With the expression (6.8) of the first-order field, this forcing reads

(6.11)\begin{equation} \boldsymbol F_2 ={-} \tilde\alpha\,\nabla^2 \boldsymbol{u}_0 -A^2 (\hat{\boldsymbol u}_1^A\boldsymbol{\cdot}\boldsymbol{\nabla})\hat{\boldsymbol u}_1^A -B^2 (\hat{\boldsymbol u}_1^B \boldsymbol{\cdot}\boldsymbol{\nabla}) \hat{\boldsymbol u}_1^B -AB \mathcal{C} (\hat{\boldsymbol u}_1^A , \hat{\boldsymbol u}_1^B). \end{equation}

All terms in (6.11) are potentially resonant since they are forcing at zero frequency an operator that is singular precisely at zero frequency, as expressed by (6.5). Considering the spatial symmetries of these terms, however, shows that none of them is resonant, in a manner reminiscent of the cross-junction studied by Bongarzone et al. (Reference Bongarzone, Bertsch, Renaud and Gallaire2021). Indeed, $\nabla ^2 \boldsymbol {u}_0$, $(\hat {\boldsymbol u}_1^A\boldsymbol {\cdot }\boldsymbol {\nabla })\hat {\boldsymbol u}_1^A$ and $(\hat {\boldsymbol u}_1^B\boldsymbol {\cdot }\boldsymbol {\nabla })\hat {\boldsymbol u}_1^B$ are $S_y S_z$-symmetric, and $(\hat {\boldsymbol u}_1^A \boldsymbol {\cdot }\boldsymbol {\nabla }) \hat {\boldsymbol u}_1^B + (\hat {\boldsymbol u}_1^B \boldsymbol {\cdot }\boldsymbol {\nabla }) \hat {\boldsymbol u}_1^A$ is $A_y A_z$-symmetric, whereas $\widetilde {\mathcal {A}}$ is singular to the $S_y A_z$-symmetric mode $A$ and the $A_y S_z$-symmetric mode $B$. It follows that the forced equation (6.9) can be inverted. We look for a second-order field of the form

(6.12)\begin{equation} \boldsymbol{q}_2 (\boldsymbol x, T) = \tilde\alpha \hat{\boldsymbol{q}}_2^{\alpha} + A^2 \hat{\boldsymbol{q}}_2^{A^2} + B^2 \hat{\boldsymbol{q}}_2^{B^2} + AB \hat{\boldsymbol{q}}_2^{AB}, \end{equation}

where each term is the response to the individual forcing terms in (6.11):

(6.13)\begin{gather} \mathcal{B}\,\partial_t {\boldsymbol q}_2^\alpha + \widetilde{\mathcal{A}} {\boldsymbol q}_2^\alpha ={-} \tilde\alpha (\nabla^2 \boldsymbol{u}_0, 0)^{\rm T}, \end{gather}

(6.14)\begin{gather}\mathcal{B}\,\partial_t {\boldsymbol q}_2^{A^2} + \widetilde{\mathcal{A}} {\boldsymbol q}_2^{A^2} ={-}A^2 ( (\hat{\boldsymbol u}_1^A\boldsymbol{\cdot}\boldsymbol{\nabla})\hat{\boldsymbol u}_1^A, 0)^{\rm T}, \end{gather}

(6.15)\begin{gather}\mathcal{B}\,\partial_t {\boldsymbol q}_2^{B^2} + \widetilde{\mathcal{A}} {\boldsymbol q}_2^{B^2} ={-}B^2 ( (\hat{\boldsymbol u}_1^B \boldsymbol{\cdot}\boldsymbol{\nabla}) \hat{\boldsymbol u}_1^B, 0)^{\rm T}, \end{gather}

(6.16)\begin{gather}\mathcal{B}\,\partial_t {\boldsymbol q}_2^{AB} + \widetilde{\mathcal{A}} {\boldsymbol q}_2^{AB} ={-}AB (\mathcal{C} (\hat{\boldsymbol u}_1^A , \hat{\boldsymbol u}_1^B), 0)^{\rm T}. \end{gather}

For each of the above problems, boundary conditions are similar to those used for the eigenmodes. In particular, the symmetry of the forcing terms imposes the symmetry of the second-order fields, therefore we use on the symmetry planes the corresponding boundary conditions (5.7) or (5.10). Figure 18 shows the second-order fields. As expected from the symmetries of the forcing terms, $\boldsymbol u_2^{A^2}$, $\boldsymbol u_2^{B^2}$ and $\boldsymbol u_2^\alpha$ are all doubly symmetric ($S_y S_z$), while $\boldsymbol u_2^{AB}$ is doubly antisymmetric ($A_y A_z$).

Figure 18.

Second-order fields, shown in the planes $z=0$ (top view), $y=0$ (side view) and $x=2.5$ (rear view), for (a) $\boldsymbol u_2^{A^2}$, (b) $\boldsymbol u_2^{B^2}$, (c) $\boldsymbol u_2^{AB}$, and (d) $\boldsymbol u_2^\alpha$.

6.1.3. Order $\epsilon ^3$

At order $\epsilon ^3$, the third-order field $\boldsymbol q_3$ is a solution of the linearised NS equations

(6.17)\begin{equation} \mathcal{B}\,\partial_t {\boldsymbol q}_3 + \widetilde{\mathcal{A}} {\boldsymbol q}_3 = (\boldsymbol F_3,0)^{\rm T}, \end{equation}

forced by a term that depends on lower-order fields only,

(6.18)\begin{equation} \boldsymbol F_3 ={-}\partial_T \boldsymbol{u}_1 - \tilde\alpha\,\nabla^2 \boldsymbol{u}_1-\mathcal{C}(\boldsymbol{u}_1 , \boldsymbol{u}_2) + \mathcal{S} \boldsymbol{u}_1. \end{equation}

The first term in (6.18) is due to the slow variation of the first-order field $\boldsymbol {u}_1$, i.e. the yet unknown amplitudes $A(T)$ and $B(T)$. The second term is due to Reynolds number variations in the first-order field $\boldsymbol {u}_1$. The third term is due to the mutual transport of the first- and second-order fields $\boldsymbol {u}_1$ and $\boldsymbol {u}_2$. The last term is due to the non-zero values of the growth rates $\sigma _A$ and $\sigma _B$ at $Re=Re_c$, accounted for by the action of the shift operator. With the expressions (6.8) and (6.12) of the first- and second-order fields, the forcing becomes

(6.19)\begin{align} \boldsymbol F_3 &= \left( -\partial_T A - A \tilde\alpha\,\nabla^2 + A \sigma_A \right) \hat{\boldsymbol u}_1^A+ \left(-\partial_T B- B \tilde\alpha\,\nabla^2+ B \sigma_B \right) \hat{\boldsymbol u}_1^B \nonumber\\ &\quad -\mathcal{C}(A \hat{\boldsymbol u}_1^A + B \hat{\boldsymbol u}_1^B, \tilde\alpha \hat{\boldsymbol{u}}_2^{\alpha}+ A^2 \hat{\boldsymbol{u}}_2^{A^2} + B^2 \hat{\boldsymbol{u}}_2^{B^2}+ AB \hat{\boldsymbol{u}}_2^{AB}). \end{align}

It turns out that all terms in (6.19) are either $S_y A_z$-symmetric or $A_y S_z$-symmetric, and are therefore resonant with mode $A$ or mode $B$. To avoid secular terms, i.e. to ensure that (6.17) can be inverted, a compatibility condition must be enforced. The Fredholm alternative states that the forcing $\boldsymbol {F}_3$ must be orthogonal to the kernel of the adjoint linearised NS operator. This leads to the following equations for the amplitudes $A$ and $B$:

(6.20)\begin{gather} \partial_t A = \lambda_A A - \chi_A A^3 - \eta_A B^2 A, \end{gather}

(6.21)\begin{gather}\partial_t B = \lambda_B B - \chi_B B^3 - \eta_B A^2 B, \end{gather}

where the fast time $t=T/\epsilon ^2$ has been reintroduced. All the WNL coefficients in (6.20)–(6.21), whose values are reported in Appendix A, are computed as scalar products between the adjoint modes $\hat {\boldsymbol u}_1^{A{\dagger} }$, $\hat {\boldsymbol u}_1^{B{\dagger} }$ and resonant forcing terms:

(6.22)\begin{gather} \lambda_A = \epsilon^2 \tilde\lambda_A = \sigma_A -\alpha \langle \hat{\boldsymbol{u}}_1^{A{\dagger}}, \mathcal{C}( \hat{\boldsymbol{u}}_2^\alpha, \hat{\boldsymbol{u}}_1^A) + \nabla^2 \hat{\boldsymbol{u}}_1^A \rangle, \end{gather}

(6.23)\begin{gather}\lambda_B = \epsilon^2 \tilde\lambda_B = \sigma_B -\alpha \langle \hat{\boldsymbol{u}}_1^{B{\dagger}}, \mathcal{C}( \hat{\boldsymbol{u}}_2^\alpha, \hat{\boldsymbol{u}}_1^B) + \nabla^2 \hat{\boldsymbol{u}}_1^B \rangle, \end{gather}

(6.24)\begin{gather}\eta_A = \epsilon^2 \tilde\eta_A = \epsilon^2 \langle \hat{\boldsymbol{u}}_1^{A{\dagger}}, \mathcal{C}(\hat{\boldsymbol{u}}_2^{B^2}, \hat{\boldsymbol{u}}_1^A) + \mathcal{C}(\hat{\boldsymbol{u}}_2^{AB}, \hat{\boldsymbol{u}}_1^B) \rangle, \end{gather}

(6.25)\begin{gather}\eta_B = \epsilon^2 \tilde\eta_B = \epsilon^2 \langle \hat{\boldsymbol{u}}_1^{B{\dagger}}, \mathcal{C}(\hat{\boldsymbol{u}}_2^{A^2}, \hat{\boldsymbol{u}}_1^B) + \mathcal{C}(\hat{\boldsymbol{u}}_2^{AB}, \hat{\boldsymbol{u}}_1^A) \rangle, \end{gather}

(6.26)\begin{gather}\chi_A = \epsilon^2 \tilde\chi_A = \epsilon^2 \langle \hat{\boldsymbol{u}}_1^{A{\dagger}}, \mathcal{C}( \hat{\boldsymbol{u}}_2^{A^2} , \hat{\boldsymbol{u}}_1^A ) \rangle, \end{gather}

(6.27)\begin{gather}\chi_B = \epsilon^2 \tilde\chi_B = \epsilon^2 \langle \hat{\boldsymbol{u}}_1^{B{\dagger}}, \mathcal{C}(\hat{\boldsymbol{u}}_2^{B^2} , \hat{\boldsymbol{u}}_1^B) \rangle, \end{gather}

where we recall that the adjoint and direct modes are normalised such that $\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$.

For any given solution of the amplitude equations (6.20)–(6.21) (see § 6.2), the total flow field up to second order can be reconstructed as

(6.28)\begin{equation} \boldsymbol{q} = \boldsymbol q_0 + \epsilon \left(A \hat{\boldsymbol q}_1^A + B \hat{\boldsymbol q}_1^B\right) + \epsilon^2 \left(\tilde\alpha \hat{\boldsymbol{q}}_2^{\alpha} + A^2 \hat{\boldsymbol{q}}_2^{A^2} + B^2 \hat{\boldsymbol{q}}_2^{B^2} + AB \hat{\boldsymbol{q}}_2^{AB}\right)+ O(\epsilon^3). \end{equation}

We note that the values of the WNL coefficients $\lambda _A, \ldots, \chi _B$, and of the amplitudes $A,B$, depend on the choices of $\epsilon$ and of the direct/adjoint modes normalisation, but the final reconstructed field (6.28) itself is independent of these specific choices.

6.2. Bifurcation diagram

The amplitude equations (6.20)–(6.21) are coupled ordinary differential equations. They are partly similar to the classic Stuart–Landau equation describing the pitchfork/Hopf bifurcation of a single stationary/oscillating mode, but differ by the two nonlinear coupling terms $\eta _A A B^2$ and $\eta _B B A^2$. Their structure is well known in the literature; see, for instance, Kuznetsov (Reference Kuznetsov2004) for a general mathematical treatment, and Zhu & Gallaire (Reference Zhu and Gallaire2017) and Bongarzone et al. (Reference Bongarzone, Bertsch, Renaud and Gallaire2021) for fluid mechanical examples with two modes, one stationary and one oscillatory, that break two different spatial symmetries. There are in general four equilibrium solutions ($\partial _t A = \partial _t B = 0$), as follows.

1. (i) Symmetric state: $A=B=0$.

2. (ii) Pure state $A$ (top/down symmetry breaking):

   (6.29a,b)\begin{equation} A^2 = \frac{\lambda_A}{\chi_A},\quad B=0. \end{equation}

3. (iii) Pure state $B$ (left/right symmetry breaking):

   (6.30a,b)\begin{equation} A=0,\quad B^2 = \frac{\lambda_B}{\chi_B}. \end{equation}

4. (iv) Mixed state $(A,B)$ (double symmetry breaking):

   (6.31a,b)\begin{equation} A^2 = \frac{\chi_B \lambda_A - \eta_A \lambda_B}{\chi_A \chi_B -\eta_A\eta_B},\quad B^2 = \frac{\chi_A \lambda_B - \eta_B \lambda_A}{\chi_A \chi_B -\eta_A\eta_B}. \end{equation}

We assess the linear stability of each state in two ways: (i) by computing the eigenvalues of the Jacobian $\boldsymbol{\mathsf{J}}$ of the system linearised about the state of interest,

(6.32)\begin{equation} \boldsymbol{\mathsf{J}} = \left[ \begin{array}{cc} \lambda_A - 3 \chi_A A^2 - \eta_A B^2 & - 2 \eta_A B A \\ - 2 \eta_B A B & \lambda_B - 3 \chi_B B^2 - \eta_B A^2 \end{array}\right],\end{equation}

that state being stable if both eigenvalues have a negative real part; (ii) by solving in time the amplitude equations (6.20)–(6.21) from an initial condition close to the state of interest, that state being stable if it is an attractor of the final solution.

We obtain the bifurcation diagram shown in figure 19 for the reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$), using $Re_c=300$ and $\epsilon =0.1$. We use a standard representation, where solid lines indicate linearly stable states, and dashed lines indicate unstable states. Pitchfork bifurcations are invariant under reflections $A \rightarrow -A$ and $B \rightarrow -B$, so it is sufficient to report results for $|A|$ and $|B|$ (i.e. for half of the symmetric diagram). The bifurcation sequence is as follows.

1. (1) For $Re_c < 293$, only the symmetric base flow is stable.

2. (2) At $Re = Re_c^A=293$, mode $A$ bifurcates, leading to the stable pure state $(A,0)$ and breaking the top/down symmetry of the base flow.

3. (3) At $Re = Re_c^B=304$, mode $B$ bifurcates, but the pure state $(0,B)$ is unstable at first, and $(A,0)$ remains the only stable state.

4. (4) At $Re=314$, the pure state $(0,B)$ becomes stable, thus leading to bistability: either single symmetry-breaking state (top/down or left/right) can be expected to be observed.

5. (5) At $Re=330$, $(A,0)$ becomes unstable, leaving $(0,B)$ as the only stable state at larger Reynolds numbers.

Figure 19.

Bifurcation diagram associated with the amplitude equations (6.20)–(6.21) obtained from a WNL analysis, for the reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$), using $Re_c=300$ and $\epsilon =0.1$. Stable and unstable states are shown with solid and dashed lines, respectively. (a) The 3-D view in the $(A,B,Re)$ space. Insets show a rear view of the reconstructed field (see figure 20). (b) The 2-D views in the $(|A|,Re)$ and $(|B|,Re)$ spaces. The shaded region corresponds to bistability and hysteresis.

We note that the mixed state $(A,B)$ exists in the bistable region but is never stable, meaning that this double symmetry-breaking state is not expected to be observed in a stable manner. Figure 20 shows the four states at $Re=325$, reconstructed up to second order according to (6.28). The wake is deflected clearly in the vertical and horizontal directions in the (stable) pure states $(A,0)$ and $(0,B)$, respectively, and in both directions in the (unstable) mixed state $(A,B)$.

Figure 20.

Streamwise velocity $u$ of reconstructed flows (6.28) at $Re=325$, from the WNL analysis up to second order, shown in the planes $z=0$ (top view), $y=0$ (side view) and $x=2.5$ (rear view), for states (a) $(0,0)$ (unstable), (b) $(A,B)$ (unstable), (c) $(A,0)$ (stable), and (d) $(0,B)$ (stable).

Interestingly, although mode $A$ bifurcates before mode $B$, state $(A,0)$ quickly becomes unstable, and eventually $(0,B)$ remains the only stable state. In other words, the top/down symmetry-breaking state cannot be observed except in the range $293 \leq Re \leq 330$, while the left/right symmetry-breaking state can be observed for $Re \geq 316$. This is precisely the same symmetry breaking as that observed in the turbulent regime, where Ahmed body wakes deflect along the longer dimension of the body. Therefore, unlike wakes past 2-D cylinders and 3-D disks, spheres and bullet-shaped axisymmetric bluff bodies, the preferred turbulent symmetry breaking in the Ahmed body wake is reminiscent not of the first laminar bifurcation, but of the second one. Our WNL analysis rationalises this observation by unveiling the bifurcation sequence. In § 6.4, we will demonstrate the robustness of this sequence for other geometries typical of Ahmed bodies.

6.3. Comparison with DNS

In this subsection, we perform fully nonlinear DNS of the flow past our reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$) at various Reynolds numbers, with the aim of confirming the bifurcation sequence obtained in § 6.2. In particular, we investigate whether the single symmetry-breaking states $(A,0)$ and $(0,B)$ are indeed observed, and if they appear and disappear in the same order. We also look for evidence of bistability and hysteresis, increasing and decreasing the Reynolds number in two independent sequences: $Re=290 \rightarrow 310 \rightarrow 330$, etc., and $Re=410 \rightarrow 390 \rightarrow 370$, etc. In each case, simulations are run for a sufficiently long time to allow the flow to reach a truly stationary state. We determine the flow state by monitoring the horizontal velocity $v(x,0,0)$ and the vertical velocity $w(x,0,0)$ on the symmetry axis. Specifically, we should obtain:

1. (i) $v=w=0$ for the symmetric base flow;

2. (ii) $v=0$, $w\neq 0$ for the top/down symmetry-breaking state $(A,0)$;

3. (iii) $v\neq 0$, $w=0$ for the left/right symmetry-breaking state $(0,B)$;

4. (iv) $v\neq 0$, $w\neq 0$ for the double symmetry-breaking state $(A,B)$, although we do not expect to observe this state in the stationary regime.

Figure 21 shows the velocities obtained in the stationary regime at the location $(x^*,0,0) = (2.5,0,0)$. In all simulations, the final flow is steady, without any oscillations. When increasing $Re$, the vertical velocity $w$ first becomes non-zero when $Re \simeq 300$, while $v=0$. In other words, state $(A,0)$ is observed. This state loses stability when $Re \simeq 340$, and the flow shifts to state $(0,B)$, with a non-zero horizontal velocity $v$, while $w=0$, for this and larger Reynolds numbers. When decreasing $Re$, state $(0,B)$ persists until $Re \simeq 320$, at which point the flow shifts back to state $(A,0)$. Therefore, this suggests hysteresis in the approximate range of Reynolds numbers $[315 \pm 10, 340 \pm 10]$.

Figure 21.

Vertical and horizontal velocities in the wake in $(x^*,0,0) = (2.5,0,0)$, obtained from DNS of the reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$). Upward-pointing symbols ($\vartriangle$) correspond to increasing Reynolds numbers, and downward-pointing symbols ($\triangledown$) correspond to decreasing Reynolds numbers. The shaded region shows bistability and hysteresis.

It is worth mentioning that the flow transitions between $(A,0)$ and $(0,B)$ via a state reminiscent of the double symmetry-breaking state $(A,B)$ (rather than via the symmetric base flow), although this state is visited only transiently. We also note that transitions are rather slow (of the order of $10^3$ convective time units), consistent with the small growth rates of the Jacobian (6.32) calculated near the corresponding transitions, i.e. near the boundaries between stages 3 and 4, and between stages 4 and 5, in the bifurcation diagram of figure 19.

The qualitative agreement between DNS and WNL analysis is excellent for the transition Reynolds numbers between stages 1, 2 and 3. At larger $Re$, the agreement deteriorates for the transition Reynolds numbers between stages 3, 4 and 5, which may be ascribed to the increasing departure from criticality. Nonetheless, the whole bifurcation sequence is well predicted by the WNL analysis.

As a concluding remark, we note that since the two bifurcations of interest are stationary, it would also be possible to compute the fully nonlinear states with a continuation method or, alternatively, with the self-consistent model proposed by Camarri & Mengali (Reference Camarri and Mengali2019).

6.4. Effect of $W$ and $L$

We have presented the bifurcation diagram of our reference Ahmed body ($W=1.2$, $L=3$, $R=0.3472$) in § 6.2. One may wonder whether this sequence is robust to geometric modifications. Indeed, Ahmed bodies studied in the literature show some variations in dimension, especially in width and length (table 2). Therefore we investigate the effect of $W$ and $L$, keeping the fillet radius constant ($R=0.3472$). Specifically, we consider narrower and wider bodies ($W \in [1, 1.35]$, $L=3$), as well as shorter and longer bodies ($W=1.2$, $L \in [2.6, 3.8]$). We compute the base flow and eigenmodes for these bodies, and determine the critical Reynolds numbers $Re_c^A$ and $Re_c^B$. We then set $Re_c$ by rounding $(Re_c^A+Re_c^B)/2$ to the nearest multiple of 5, and finally we compute the coefficients of the amplitude equations (6.20)–(6.21).

Table 2.

Typical Ahmed body widths and lengths (normalised by the body height) found in the literature.

A convenient way to summarise the results is to consider the different possible bifurcation diagrams for this type of amplitude equation. As detailed in Kuznetsov (Reference Kuznetsov2004), there are five topologically different bifurcation diagrams, corresponding to different regions of the $\theta$–$\delta$ plane, where $\theta$ and $\delta$ are mixed ratios of the nonlinear coefficients of (6.20)–(6.21):

(6.33a,b)\begin{equation} \theta = \frac{\eta_A}{\chi_B},\quad \delta = \frac{\eta_B}{\chi_A}. \end{equation}

Figure 22(a) shows the five regions, labelled I to V, for the case $\theta \geq \delta$ relevant here. As will soon become clear, regions I ($\theta >0$, $\delta >0$, $\theta \delta >1$) and II ($\theta >0$, $\delta >0$, $\theta \delta <1$) are of particular interest to our study. In region I, the bifurcation diagram is as sketched in figure 23(a), and is topologically equivalent to that obtained in § 6.2. We recall briefly the bifurcation sequence for increasing $Re$:

1. (1) only the trivial state $(0,0)$ is stable;

2. (2) mode $A$ bifurcates and the stable state $(A,0)$ appears;

3. (3) mode $B$ bifurcates and the unstable state $(0,B)$ appears;

4. (4) state $(B,0)$ becomes stable and the unstable mixed state $(A,B)$ appears;

5. (5) state $(A,B)$ disappears and state $(A,0)$ becomes unstable – only $(0,B)$ remains stable.

Figure 22.

(a) Five regions of the $\theta$–$\delta$ plane, associated with the five possible bifurcation diagrams for the amplitude equations (6.20)–(6.21). (b) Variation of $\theta$ and $\delta$ with $W$ and $L$ for rectangular prisms of fixed fillet radius $R=0.3472$. Dashed lines are a guide to the eye. In the considered range of widths $1 \leq W \leq 1.35$, and lengths $2.6 \leq L \leq 3.8$, representative of most Ahmed bodies, $\theta$ and $\delta$ stay in region I, i.e. the bifurcation diagram remains as in figures 19 and 23(a).

Figure 23.

Schematic bifurcation diagrams in (a) region I and (b) region II of the $\theta$–$\delta$ plane of figure 22. Stable and unstable states are shown with solid and dashed lines, respectively. (c) Corresponding phase portraits in the $A$–$B$ amplitude plane. Stable and unstable states are shown with filled and open symbols, respectively. Dashed arrows in stage 4 show how the mixed state $(A,B)$ moves with $Re$.

In the neighbouring region II, the bifurcation diagram is as sketched in figure 23(b), and corresponds to another bifurcation sequence:

1. (1) only the trivial state $(0,0)$ is stable;

2. (2) mode $A$ bifurcates and the stable state $(A,0)$ appears;

3. (3) mode $B$ bifurcates and the unstable state $(0,B)$ appears;

4. (4) state $(A,0)$ becomes unstable and the stable mixed state $(A,B)$ appears;

5. (5) state $(A,B)$ disappears and state $(0,B)$ becomes stable.

Although the initial stages (1–3) and the final stage (5) are identical in regions I and II, what happens in between (stage 4, shaded region in figure 23a) is fundamentally different: in region I the pure states $(A,0)$ and $(B,0)$ are simultaneously stable (bistability) while the mixed state $(A,B)$ cannot be observed; by contrast, in region II the pure states $(A,0)$ and $(B,0)$ are unstable and only the mixed state $(A,B)$ can be observed. This is also apparent in the phase portraits of figure 23(c). From stages 3 to 5, the stability of the pure states $(A,0)$ and $(0,B)$ is exchanged. In region I, this happens as the unstable mixed state $(A,B)$ is born in $(0,B)$ (making it stable) and dies when colliding with $(A,0)$ (making it unstable). In region II, the stable mixed state $(A,B)$ is born in $(A,0)$ (making it unstable) and dies when colliding with $(0,B)$ (making it stable).

Figure 22(b) shows that all the Ahmed bodies considered in this study fall in region I of the $\theta$–$\delta$ plane. Our reference Ahmed body corresponds to $\theta = 2.68$, $\delta =0.554$ (filled symbol). The effect of $L$ is very limited, with $\theta$ and $\delta$ barely varying compared to their reference values (inset). The effect of $W$ is much more significant, with $\theta$ and $\delta$ varying substantially: the point $(\theta,\delta )$ moves away from the boundary between regions I and II as the body becomes wider, and vice versa. By symmetry, the limiting case $W=1$ (body as wide as tall, seldom used in practice for Ahmed bodies) lies exactly on the boundary. We can conclude that in the WNL framework, the bifurcation sequence obtained in § 6.2 is robust to width and length variations commonly encountered for Ahmed bodies.

One must keep in mind that the WNL analysis is rigorously valid in the vicinity of $Re_c$ only. Therefore, one cannot completely rule out a fully nonlinear bifurcation diagram that differs from the WNL one, e.g. an incursion in region II. More generally, it is likely that the deflected wake undergoes one or several secondary instabilities at sufficiently large Reynolds number. Nonetheless, the fact that stages 1–3 and 5 are identical in regions I and II suggests that the transition from state $(A,0)$ to state $(0,B)$ is, in any case, a robust feature of this flow. Furthermore, the fully nonlinear DNS bifurcation sequence (§ 6.3) is in very good agreement with the WNL prediction for the reference Ahmed body, which attests to the reliability of the WNL analysis.