2.1 Direct Navier–Stokes solver

The unsteady Navier–Stokes solver is based on fully implicit time integration and finite-element method (FEM) discretization (Noack & Morzyński Reference Noack and Morzyński2017; Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003, Reference Noack, Stankiewicz, Morzyński and Schmid2016). The time integration is third-order accurate while FEM discretization employs second–order Taylor–Hood finite elements (Taylor & Hood Reference Taylor and Hood1973). The solution is obtained iteratively, with the Newton–Raphson type approach. The tangent matrix is updated at each iteration and computations are carried out until the residual is under a prescribed tolerance. The steady solution is obtained in a similar Newton–Raphson iteration for the steady Navier–Stokes equations. The convergence to one of the three steady solutions with different states of the base-bleeding jet is triggered by appropriate ‘initial’ conditions in the iteration, see appendix A. The solver quickly converges to one of the steady states and a final, near-zero residual confirms that this is indeed the steady flow solution sought. The computational domain is discretized on an unstructured grid. Pinball configuration uses a grid with 4225 triangles and 8633 vertices (see figure 1a). To test the grid dependency of the solution we increased the number of triangles by nearly a factor 4 (26 849 elements and 54 195 nodes, see figure 1b). The flow patterns shown in figure 2 develop from a steady solution at $Re=100$ subjected to an instantaneous rotation of cylinders at $T=0.2$. The upper cylinder rotates counterclockwise, the lower one clockwise and the centre cylinder also in a clockwise direction – all with unit circumferential velocity, i.e. the velocity of oncoming flow $U_{\infty }$. This configuration and boundary conditions result in the vortex shedding shown in figure 2 for the time instance $t=200$. Both simulations prove grid independence and yield dynamically consistent results (see figure 2).

Figure 1.

Computational grid for the fluidic pinball, with 8633 (a) and 54 195 vertices (b).

Figure 2.

Direct numerical simulation (DNS) computations on a grid with 8633 nodes (a) and 54 195 nodes (b). Vorticity depicted with colour is $[-1.5,1.5]$, $t=200$. The initial kick is provided by a rotation of all three cylinders at $t=0.2$, see text for details.

2.2 Pinball configuration

We refer to the configuration shown in figure 1 as the fluidic pinball as the rotation speeds allow one to change the paths of the incoming fluid just as flippers manipulate the ball of a conventional pinball machine. The fluidic pinball is a set of three equal circular cylinders with radius $R$ placed parallel to each other in a viscous incompressible uniform flow at speed $U_{\infty }$. The flow over a cluster of three parallel cylinders has been experimentally studied involving heat transfer, fluid–structure interactions and multiple frequency interactions over the past few decades (Price & Paidoussis Reference Price and Paidoussis1984; Sayers Reference Sayers1987; Lam & Cheung Reference Lam and Cheung1988; Tatsuno, Amamoto & Ishi-i Reference Tatsuno, Amamoto and Ishi-i1998; Bansal & Yarusevych Reference Bansal and Yarusevych2017). For the fluidic pinball, the cylinders can rotate at different speeds, creating a kaleidoscope of vortical structures or variety of steady flow solutions. The configuration is used for evaluation of flow controllers (Cornejo Maceda Reference Cornejo Maceda2017) as this problem is a challenging task for control methods comprising several frequency cross-talk mechanisms (Noack & Morzyński Reference Noack and Morzyński2017). The centres of the cylinders form an equilateral triangle with side length $3R$, symmetrically positioned with respect to the flow. The leftmost triangle vertex points upstream, while the rightmost side is orthogonal to the oncoming flow. The origin of the Cartesian coordinate system is placed in the middle of the top and bottom cylinder. The fluidic pinball computational domain, shown in figure 1, is bounded by the rectangle $[-6,20]\times [-6,6]$.

Without forcing, the boundary conditions comprise a no-slip condition on the cylinders and a unit velocity in the far field,

(2.1)$$\begin{eqnarray}U_{r}=0\quad \text{on the cylinders and }U_{\infty }=\boldsymbol{e}_{x}\text{ at infinity}.\end{eqnarray}$$

The far field boundary conditions are exerted on the inflow, upper and lower boundaries while the outflow boundary is assumed to be a stress-free one, transparent for the outgoing fluid structures. A typical initial condition is the unstable steady Navier–Stokes solution $\boldsymbol{u}_{s}(\boldsymbol{x})$.

In this study, all three cylinders remain static as we are interested in the natural dynamics of the flow as the Reynolds number is increased.

2.3 Flow features

Figure 3.

Flow states at different values of the Reynolds number: the stable states are labelled with a solid-line box and the unstable states with a dashed-line box. Steady solutions $\boldsymbol{u}_{s}(\boldsymbol{x})$ (on the left side): three steady solutions $\boldsymbol{u}_{s}(\boldsymbol{x})$ (symmetric), $\boldsymbol{u}_{s}^{\pm }(\boldsymbol{x})$ (asymmetric), exist at $Re=80$. For $Re=30$ and $Re=80$, the steady solutions are unstable: the permanent regime $\boldsymbol{u}(\boldsymbol{x},t)=\overline{\boldsymbol{u}}(\boldsymbol{x})+\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$ is unsteady with mean flow field $\overline{\boldsymbol{u}}(\boldsymbol{x})$. At $Re=80$, the instantaneous flow field $\boldsymbol{u}(\boldsymbol{x},t)=\overline{\boldsymbol{u}}(\boldsymbol{x})+\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$, transiently explored when starting close to $\boldsymbol{u}_{s}(\boldsymbol{x})$, is unstable with respect to either $\boldsymbol{u}^{\pm }(\boldsymbol{x},t)=\overline{\boldsymbol{u}}^{\pm }(\boldsymbol{x})+\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$.

The steady solution $\boldsymbol{u}_{s}$, shown in figure 3 for different values of the Reynolds number $Re$, is stable up to the critical value $Re_{1}\approx 18$. This value corresponds to $5/2\times Re_{1}\approx 45$ with respect to the actual body height $5R$, which is consistent with the critical value of the Reynolds number found for a single cylinder (Ding & Kawahara Reference Ding and Kawahara1999; Barkley Reference Barkley2006). Beyond $Re_{1}$, the steady solution becomes unstable with respect to vortices periodically and alternately shed at the top and bottom of the two rightmost cylinders, following a Hopf bifurcation (instability of the fixed point via a pair of complex-conjugated eigenvalues, see e.g. Strogatz et al. (Reference Strogatz, Friedman, Mallinckrodt and McKay1994)). In addition to the resulting von Kármán street of vortices, the gap between the cylinders makes possible the formation of a jet at the base of the two outer cylinders. The steady solution $\boldsymbol{u}_{s}$, the mean flow $\overline{\boldsymbol{u}}$ and the instantaneous flow field $\boldsymbol{u}$, are shown in figure 3, for $Re=30$.

The flow passing between the two rearward cylinders, the base-bleeding flow, has a critical impact on the successive bifurcations undergone by the system on the route to chaos. Indeed, beyond a secondary critical value $Re_{2}\approx 68$ of the Reynolds number, the system undergoes a pitchfork bifurcation, which affects both the fixed point and the limit cycle, via a real eigenvalue, see e.g. Strogatz et al. (Reference Strogatz, Friedman, Mallinckrodt and McKay1994). As a result, the symmetry of both the steady solution and the mean flow is broken with respect to the symmetry plane defined by $y=0$. This is illustrated by the two mirror-conjugated steady solutions $\boldsymbol{u}_{s}^{\pm }$ and the two associated mean flows $\overline{\boldsymbol{u}}^{\pm }$, shown in figure 3 for $Re=80$, where the base-bleeding jet appears deflected either upward or downward with respect to the symmetry plane. Note, however, that a symmetry-preserving mean flow ($\overline{\boldsymbol{u}}$ in figure 3 for $Re=80$) still exists beyond the secondary bifurcation, so that three mean flows exist beyond $Re_{2}$: two of them, $\overline{\boldsymbol{u}}^{\pm }$, are mirror conjugated and break the symmetry, while the last one, $\overline{\boldsymbol{u}}$ preserves the symmetry. This bifurcation of the limit cycle is coincident with the bifurcation of the fixed point, as three steady solutions can be found beyond $Re_{2}$: $\boldsymbol{u}_{s}^{\pm }$ are mirror conjugated and break the symmetry, while the last one, $\boldsymbol{u}_{s}$ preserves the symmetry. Yet, all three of them are unstable with respect to the cyclic shedding of von Kármán vortices, in which symmetry properties of the steady solution are succeeded in the resulting mean flow. When the initial condition is close to the symmetric steady solution $\boldsymbol{u}_{s}$, the flow regime arrives after a long transient on a limit cycle whose mean flow $\overline{\boldsymbol{u}}$ is symmetric, as illustrated in figure 3 for $Re=80$. However, the dynamics of this limit cycle is only transient, indicating that it is not a stable state. After a new transient, depending on the details of the initial condition, the flow regime eventually reaches one of the two mirror-conjugated limit cycles (centred on either $\overline{\boldsymbol{u}}^{\pm }$). When the initial condition already breaks the symmetry of the flow configuration, the unstable ‘symmetry-centred’ limit cycle is not explored and the system reaches directly one of the two stable limit cycles. The transient dynamics between these six typical states beyond $Re_{2}$ illustrates a transverse action on the original state space resulting from the new active symmetry-breaking mode decoupling with the primary Hopf bifurcation, as detailed in appendix D. The new active degree of freedom introduced by the pitchfork instability is responsible for the two simultaneous local pitchfork bifurcations, of both the steady solution and the periodic solution, as shown in the appendices B and C. The simultaneous bifurcation of the steady and periodic solution has also been observed for the cylinder wake transition from stability analyses (Noack & Eckelmann Reference Noack and Eckelmann1994a,Reference Noack and Eckelmannb) and from three-dimensional Navier–Stokes simulations (Zhang, Noack & Eckelmann Reference Zhang, Noack and Eckelmann1994). A further discussion about this non-generic situation is recorded in appendix E. Besides, the linear stability analysis (see appendix B) and the Floquet stability analysis (see appendix C) around $Re_{2}$ have been performed to prove these two simultaneous bifurcations.

Figure 4.

Bifurcation diagram based on the absolute value of the mean lift coefficient $|\overline{C}_{L}|$ (circles $+$ solid line) and its standard deviation (squares $+$ dashed line). By symmetry, each non-vanishing mean lift value is associated with a positive and negative sign for the two attractors. The vertically distributed black dots that are visible for $Re>Re_{3}$, are median values $\overline{c}_{n}=(C_{n}+C_{n+1})/2$ between successive local optima $C_{n}=C_{L}(t_{n})$ of $C_{L}(t)$, at a given Reynolds number, where the $t_{n}$ are times at which ${\dot{C}}_{L}(t_{n})=0$. Transition to unsteadiness occurs at $Re_{1}\approx 18$ (Hopf bifurcation), the average symmetry is broken beyond $Re_{2}\approx 68$ (pitchfork bifurcation), a secondary (incommensurable) frequency rises in the power spectrum at $Re_{3}\approx 104$ (Neimark–Säcker bifurcation) and transition to chaos occurs at $Re_{4}\approx 115$. Note that the symmetry is statistically recovered in the chaotic regime ($\overline{C}_{L}\approx 0$).

Figure 5.

Quasi-periodic dynamics at $Re=105$ displayed by (a) the power spectral density on time series of length $T_{data}=400$ (red curve), $T_{data}=900$ (black curve), (b) the time series and (c) the phase portrait of the lift coefficient $C_{L}$.

Figure 6.

Chaotic dynamics at $Re=130$ displayed by (a) the power spectral density on time series of length $T_{data}=400$ (red curve), $T_{data}=900$ (black curve), (b) the time series and (c) the phase portrait of the lift coefficient $C_{L}$.

As a result, when the symmetry of vortex shedding is broken, the mean value $\overline{C}_{L}$ (solid line) of the lift coefficient $C_{L}=2F_{L}/\unicode[STIX]{x1D70C}U^{2}$, where $F_{L}$ is the total lift force, no longer vanishes, as shown in figure 4. At the precision of our investigation, both the Hopf and pitchfork bifurcations were found to be supercritical.

The fluctuation amplitude of the lift coefficient is minimum for $Re\approx 80>Re_{2}$, as shown in figure 4 (dashed curve). It starts to decrease around $Re=30$, when the jet starts to grow at the base of the two outer cylinders. Henceforth, the growth of the base-bleeding jet, as the Reynolds number is increased, seems to be fed with the energy of the fluctuations. Transfers of energy between the dynamically dominant degrees of freedom will be made clear in § 3.

When the Reynolds number is further increased up to a critical value $Re_{3}\approx 104$, a new frequency rises in the power spectrum of the lift coefficient $C_{L}(t)$. This frequency is approximately one order of magnitude smaller than the natural frequency of the vortex shedding, as illustrated in figure 5(a) for $Re=105$ and in the movie QP.MP4 of the additional materials available online at https://doi.org/10.1017/jfm.2019.959. The low and natural frequencies couple to generate combs of sharp peaks in the power spectrum, while the background level depends on the length of the time series. The new frequency is associated with modulations of the base-bleeding jet around its deflected position. A visual inspection of both the time series and the phase portrait at $Re=105$ indicates that the new frequency also modulates the amplitude (figure 5b) of the main oscillator and thickens the limit cycle associated with the main oscillator (figure 5c). All these features are typical of a quasi-periodic dynamics, indicating that the system has most likely undergone a Neimark–Säcker bifurcation, e.g. a secondary Hopf bifurcation, at $Re=Re_{3}$, after which two mirror-conjugated 2-tori exist in the state space of the system.

At even larger values of the Reynolds number, $Re\geqslant Re_{4}\approx 115$, the main peak in the power spectral density of the lift coefficient widens significantly, as shown in figure 6(a) for $Re=130$. In this new regime, the instantaneous flow field is characterized by random switches between an upward or downward base-bleeding jet, see the movie CHAOS.MP4 in the additional materials, and the mean flow $\overline{\boldsymbol{u}}$ is symmetric, as shown in figure 4 for $Re>115$. The time series exhibits neither periodic nor quasi-periodic features anymore (see figure 6b) and the phase portrait exhibits a much more complex dynamics (see figure 6c). The dynamical regime henceforth exhibits many features of a chaotic regime, indicating that the system has most likely followed the Ruelle–Takens–Newhouse route to chaos (Newhouse, Ruelle & Takens Reference Newhouse, Ruelle and Takens1978).

3.1 Galerkin method

The starting point is the non-dimensionalized incompressible Navier–Stokes equations

(3.1)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\otimes \boldsymbol{u}=\unicode[STIX]{x1D708}\triangle \boldsymbol{u}-\unicode[STIX]{x1D735}p,\end{eqnarray}$$

where $\unicode[STIX]{x1D708}=1/Re$. The velocity field satisfies the no-slip condition $\boldsymbol{u}=0$ on the cylinders, the free-stream condition $\boldsymbol{u}=(1,0)$ at the inflow, a no-slip condition at the top and bottom boundary and the no-stress condition at the outflow. The steady solution $\boldsymbol{u}_{s}$ satisfies the steady Navier–Stokes equations

(3.2)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}\otimes \boldsymbol{u}_{s}=\unicode[STIX]{x1D708}\triangle \boldsymbol{u}_{s}-\unicode[STIX]{x1D735}p_{s}.\end{eqnarray}$$

The Galerkin method is based on an inner product in the space of square-integrable vector fields ${\mathcal{L}}^{2}(\unicode[STIX]{x1D6FA})$ in the observation domain $\unicode[STIX]{x1D6FA}$. The standard inner product between $\boldsymbol{u}(\boldsymbol{x})$ and $\boldsymbol{v}(\boldsymbol{x})$ reads

(3.3)$$\begin{eqnarray}(\boldsymbol{u},\boldsymbol{v})_{\unicode[STIX]{x1D6FA}}:=\int _{\unicode[STIX]{x1D6FA}}\,\text{d}\boldsymbol{x}\,\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{v}(\boldsymbol{x}).\end{eqnarray}$$

A traditional Galerkin approximation with a basic mode $\boldsymbol{u}_{0}$, for instance, the steady solution and $N$ orthonormal expansion modes $\boldsymbol{u}_{i}(\boldsymbol{x})$, $i=1,\ldots ,N$ with time-dependent amplitudes $a_{i}(t)$, reads

(3.4)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}_{0}(\boldsymbol{x})+\mathop{\sum }_{i=1}^{N}a_{i}(t)\boldsymbol{u}_{i}(\boldsymbol{x}).\end{eqnarray}$$

Orthonormality implies

(3.5)$$\begin{eqnarray}(\boldsymbol{u}_{i},\boldsymbol{u}_{j})_{\unicode[STIX]{x1D6FA}}=\unicode[STIX]{x1D6FF}_{ij},\quad i,j\in \{1,\ldots ,N\}.\end{eqnarray}$$

The projection of (3.4) on (3.1) leads to the linear–quadratic Galerkin system (Fletcher Reference Fletcher1984)

(3.6)$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}a_{i}=\unicode[STIX]{x1D708}\mathop{\sum }_{j=0}^{N}l_{ij}^{\unicode[STIX]{x1D708}}a_{j}+\mathop{\sum }_{j,k=0}^{N}q_{ijk}^{c}a_{j}a_{k}.\end{eqnarray}$$

Following Rempfer & Fasel (Reference Rempfer and Fasel1994), $a_{0}=1$ is introduced. The coefficients $l_{ij}^{\unicode[STIX]{x1D708}}=(\boldsymbol{u}_{i},\triangle \boldsymbol{u}_{j})_{\unicode[STIX]{x1D6FA}}$ and $q_{ijk}^{c}=-(\boldsymbol{u}_{i},\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{j}\otimes \boldsymbol{u}_{k})_{\unicode[STIX]{x1D6FA}}$ parametrize the viscous and convective Navier–Stokes terms. The pressure term vanishes for sufficiently large domains and is neglected in the following.

In the following, the steady solution is taken as the basic mode $\boldsymbol{u}_{0}=\boldsymbol{u}_{s}$. This implies that $\boldsymbol{a}=0$ is a fixed point of (3.6) and the constant term $\unicode[STIX]{x1D708}l_{i0}^{\unicode[STIX]{x1D708}}+q_{i00}^{\unicode[STIX]{x1D708}}=0$ vanishes as the projection of (3.2) onto the $i$th mode $\boldsymbol{u}_{i}$. In this case, equation (3.6) can be rewritten as a linear–quadratic system of ordinary differential equations

(3.7)$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}a_{i}=\mathop{\sum }_{j=1}^{N}l_{ij}a_{j}+\mathop{\sum }_{j,k=1}^{N}q_{ijk}a_{j}a_{k},\end{eqnarray}$$

where $l_{ij}=\unicode[STIX]{x1D708}l_{ij}^{\unicode[STIX]{x1D708}}+q_{ij0}^{c}+q_{i0j}^{c}$ and $q_{ijk}=q_{ijk}^{c}$ for $i,j,k\in \{1,\ldots ,N\}$.

3.2 Mean-field modelling

Mean-field modelling allows a dramatic simplification of a general Galerkin system (3.7) close to bifurcations. In this section, we derive constitutive equations with a small number of more general assumptions.

In the spirit of the Reynolds decomposition, the velocity field is decomposed into a slowly varying distorted mean flow $\boldsymbol{u}^{D}$ and fluctuation $\boldsymbol{u}^{\prime }$ with first-order (relaxational) or second-order (oscillatory) dynamics

(3.8a,b)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}^{D}(\boldsymbol{x},t)+\boldsymbol{u}^{\prime }(\boldsymbol{x},t),\quad \boldsymbol{u}^{D}(\boldsymbol{x},t)=\boldsymbol{u}_{s}(\boldsymbol{x})+\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}(\boldsymbol{x},t).\end{eqnarray}$$

Here, the mean-field deformation $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$ is the difference between the distorted mean flow and the steady solution. For the oscillatory dynamics considered, the distorted mean flow can be defined as an average over one local fluctuation period $T$ denoted by $\langle \cdot \rangle$. Thus,

(3.9)$$\begin{eqnarray}\boldsymbol{u}^{D}(\boldsymbol{x},t)=\langle \boldsymbol{u}(\boldsymbol{x},t)\rangle :=\frac{1}{T}\int _{t-T/2}^{t+T/2}\,\text{d}\unicode[STIX]{x1D70F}\,\boldsymbol{u}(\boldsymbol{x},\unicode[STIX]{x1D70F}).\end{eqnarray}$$

After the pitchfork bifurcation into two mirror-conjugated flows $\boldsymbol{u}^{+}$, $\boldsymbol{u}^{-}$, a symmetric distorted mean flow is enforced via

(3.10)$$\begin{eqnarray}\boldsymbol{u}^{D}(\boldsymbol{x},t)={\textstyle \frac{1}{2}}(\langle \boldsymbol{u}^{+}(\boldsymbol{x},t)\rangle +\langle \boldsymbol{u}^{-}(\boldsymbol{x},t)\rangle ).\end{eqnarray}$$

We note that $\boldsymbol{u}^{D}(\boldsymbol{x},t)$ is not the mean flow, which is defined by the post-transient limit

(3.11)$$\begin{eqnarray}\bar{\boldsymbol{u}}(\boldsymbol{x})=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}\boldsymbol{u}(\boldsymbol{x},\unicode[STIX]{x1D70F})\,\text{d}\unicode[STIX]{x1D70F}.\end{eqnarray}$$

The distorted mean flow coincides with mean flow for the post-transient phase before the pitchfork bifurcation. Technically, a harmonic fluctuation is assumed and this one period average is computed as an average of all phases in $[0,2\unicode[STIX]{x03C0}]$. The somewhat loaded term ‘distorted mean flow’ is directly adopted from the original publications of mean-field theory (Stuart Reference Stuart1958). J. T. Stuart considers this flow as ‘distorted’ from the steady solution by the Reynolds stress associated with the instability mode(s).

For a nominally symmetric cylindrical obstacle, the distorted mean flow can be expected to be symmetric while the dominant fluctuation is antisymmetric. This leads to a symmetry-based decomposition of the flow into a symmetric contribution $\boldsymbol{u}^{s}=(u^{s},v^{s})\in {\mathcal{U}}^{s}$ with

(3.12a,b)$$\begin{eqnarray}u^{s}(x,-y)=u^{s}(x,y),\quad v^{s}(x,-y)=-v^{s}(x,y),\end{eqnarray}$$

and an antisymmetric component $\boldsymbol{u}^{a}=(u^{a},v^{a})\in {\mathcal{U}}^{a}$ satisfying

(3.13a,b)$$\begin{eqnarray}u^{a}(x,-y)=-u^{a}(x,y),\quad v^{a}(x,-y)=v^{a}(x,y).\end{eqnarray}$$

Here, ${\mathcal{U}}^{s}$ and ${\mathcal{U}}^{a}$ denote the set of symmetric and antisymmetric vector fields, respectively. The resulting decomposition reads

(3.14)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}^{s}(\boldsymbol{x},t)+\boldsymbol{u}^{a}(\boldsymbol{x},t).\end{eqnarray}$$

In the sequel, we will identify the distorted mean flow with the symmetric component and the fluctuation with the antisymmetric one

(3.15a,b)$$\begin{eqnarray}\boldsymbol{u}^{D}(\boldsymbol{x},t)=\boldsymbol{u}^{s}(\boldsymbol{x},t),\quad \boldsymbol{u}^{\prime }(\boldsymbol{x},t)=\boldsymbol{u}^{a}(\boldsymbol{x},t).\end{eqnarray}$$

This identification is justified for symmetry-breaking bifurcations with first- or second-order dynamics with neglected higher harmonics. For brevity, ${\mathcal{U}}^{s}$ and ${\mathcal{U}}^{a}$ are introduced as symmetric and antisymmetric subsets of ${\mathcal{L}}^{2}(\unicode[STIX]{x1D6FA})$.

The convective term is easily shown to have the following symmetry properties:

(3.16a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{s}\otimes \boldsymbol{u}^{s}\in {\mathcal{U}}^{s}, & \displaystyle\end{eqnarray}$$

(3.16b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{a}\otimes \boldsymbol{u}^{a}\in {\mathcal{U}}^{s}, & \displaystyle\end{eqnarray}$$

(3.16c)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{s}\otimes \boldsymbol{u}^{a}\in {\mathcal{U}}^{a}, & \displaystyle\end{eqnarray}$$

(3.16d)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{a}\otimes \boldsymbol{u}^{s}\in {\mathcal{U}}^{a}. & \displaystyle\end{eqnarray}$$

The antisymmetric component is derived starting with (3.1), subtracting the steady version of (3.2) and exploiting the symmetry of $\boldsymbol{u}^{D}$, the antisymmetry of $\boldsymbol{u}^{\prime }$ as well as the symmetry relations (3.16). The fluctuation dynamics reads:

(3.17)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{\prime }+\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\boldsymbol{u}^{D}\otimes \boldsymbol{u}^{\prime }+\boldsymbol{u}^{\prime }\otimes \boldsymbol{u}^{D}]=\unicode[STIX]{x1D708}\triangle \boldsymbol{u}^{\prime }-\unicode[STIX]{x1D735}p^{\prime }.\end{eqnarray}$$

Analogously, the symmetric part describes the distorted mean flow dynamics

(3.18)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\boldsymbol{u}_{s}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}+\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\otimes \boldsymbol{u}_{s}+\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}+\boldsymbol{u}^{\prime }\otimes \boldsymbol{u}^{\prime }]=\unicode[STIX]{x1D708}\triangle \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}-\unicode[STIX]{x1D735}p_{\unicode[STIX]{x1D6E5}}.\end{eqnarray}$$

Note that this symmetric component of the Navier–Stokes equations has not yet been averaged and the sum of equations (3.2), (3.18) and (3.17) leads to the Navier–Stokes equations (3.1). To this point, all equations are strict identities for the symmetric and antisymmetric parts of the Navier–Stokes dynamics.

Next, we follow mean-field arguments and consider $\boldsymbol{u}^{\prime }$ and $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$ as small perturbations around the fixed point $\boldsymbol{u}_{s}$. Let $\boldsymbol{u}^{\prime }\in O(\unicode[STIX]{x1D700})$ and $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\in O(\unicode[STIX]{x1D6FF})$ where $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D6FF}$ are smallness parameters. Hence, $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\in O(\unicode[STIX]{x1D6FF}^{2})$ can be neglected in comparison to the $O(\unicode[STIX]{x1D6FF})$ terms $\boldsymbol{u}_{s}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$, $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\otimes \boldsymbol{u}_{s}$. We follow Stuart’s original idea to separate the fluctuation $\boldsymbol{u}^{\prime }$ driven by the instability and the resulting mean-field deformation $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$ and arrive at the unsteady linearized Reynolds equation

(3.19)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\boldsymbol{u}_{s}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}+\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\otimes \boldsymbol{u}_{s}+\langle \boldsymbol{u}^{\prime }\otimes \boldsymbol{u}^{\prime }\rangle ]=\unicode[STIX]{x1D708}\triangle \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}-\unicode[STIX]{x1D735}p_{\unicode[STIX]{x1D6E5}}.\end{eqnarray}$$

The mean-field deformation $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$ characterized by the scale $\unicode[STIX]{x1D6FF}$ is seen to respond linearly to the Reynolds stress force $-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle \boldsymbol{u}^{\prime }\otimes \boldsymbol{u}^{\prime }\rangle$ scaling with $\unicode[STIX]{x1D716}^{2}$. Hence, $\unicode[STIX]{x1D6FF}\sim \unicode[STIX]{x1D716}^{2}$.

Summarizing, equations (3.17) and (3.19) are the constitutive equations of mean-field theory exploiting only symmetry and smallness of the mean-field deformation.

Close to the critical Reynolds number, $Re_{c}$, the temporal growth rate can be Taylor expanded to $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D6FC}(Re-Re_{c})$ and can be assumed to be small. In this case, $\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\in O(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FF})$, i.e. the time derivative of (3.19) can be neglected with respect to the other terms $\in O(\unicode[STIX]{x1D6FF})=O(\unicode[STIX]{x1D700}^{2})$. This leads to the steady linearized Reynolds equation

(3.20)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\boldsymbol{u}_{s}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}+\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}\otimes \boldsymbol{u}_{s}+\langle \boldsymbol{u}^{\prime }\otimes \boldsymbol{u}^{\prime }\rangle ]=\unicode[STIX]{x1D708}\triangle \boldsymbol{u}_{\unicode[STIX]{x1D6E5}}-\unicode[STIX]{x1D735}p_{\unicode[STIX]{x1D6E5}}.\end{eqnarray}$$

This equation is also true for the post-transient solution, e.g. the limit cycle of a Hopf bifurcation or the asymmetric state of a pitchfork bifurcation. Often, the distorted mean flow $\boldsymbol{u}^{D}$ quickly responds to the Reynolds stress even far away from the bifurcation.

3.3 Supercritical Hopf bifurcation

At low Reynolds numbers, a symmetric stable steady solution $\boldsymbol{u}_{s}\in {\mathcal{U}}^{s}$ is observed. Periodic vortex shedding sets in with the occurrence of an unstable oscillatory antisymmetric eigenmode at $Re\geqslant Re_{1}$. The Reynolds number dependent initial growth rate and frequency are denoted by $\unicode[STIX]{x1D70E}_{1}$ and $\unicode[STIX]{x1D714}_{1}$, respectively. The real and imaginary parts of this eigenmode are $\boldsymbol{u}_{1}$ and $\boldsymbol{u}_{2}$, respectively, both antisymmetric modes. In the following, these modes are assumed to be orthonormalized.

This oscillation generates a Reynolds stress, which changes the mean flow via (3.19). The mean flow deformation is described by the symmetric shift mode $\boldsymbol{u}_{3}$ with unit norm. By symmetry, the first two modes are orthogonal with respect to the shift mode. Thus, the modes form an orthonormal basis. The resulting Galerkin expansion reads as follows:

(3.21)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}_{s}(\boldsymbol{x})+\underbrace{a_{1}(t)\boldsymbol{u}_{1}(\boldsymbol{x})+a_{2}(t)\boldsymbol{u}_{2}(\boldsymbol{x})}_{\boldsymbol{u}^{\prime }}+\underbrace{a_{3}(t)\boldsymbol{u}_{3}(\boldsymbol{x})}_{\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}}.\end{eqnarray}$$

Moreover, polar coordinates are introduced $a_{1}(t)=r(t)\cos \unicode[STIX]{x1D703}(t)$, $a_{2}(t)=r(t)\sin \unicode[STIX]{x1D703}(t)$, $\text{d}\unicode[STIX]{x1D703}/\text{d}t=\unicode[STIX]{x1D714}(t)$ where $r$ and $\unicode[STIX]{x1D714}$ are assumed to be slowly varying functions of time.

Substituting (3.21) in (3.17), projecting on $\boldsymbol{u}_{i}$, $i=1,2$ and applying the Krylov–Bogoliubov (Jordan & Smith Reference Jordan and Smith1999) averaging method yields

(3.22a)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{1}/\text{d}t=\unicode[STIX]{x1D70E}a_{1}-\unicode[STIX]{x1D714}a_{2},\quad \unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}, & \displaystyle\end{eqnarray}$$

(3.22b)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{2}/\text{d}t=\unicode[STIX]{x1D70E}a_{2}+\unicode[STIX]{x1D714}a_{1},\quad \unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}a_{3}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D714}_{1},\unicode[STIX]{x1D6FD}>0$

Krylov–Bogoliubov averaging implies a harmonic balancing on the slowly varying amplitude and frequency of oscillatory $a_{1,2}$ and the slowly varying $a_{3}$ dynamics. The corresponding original theorem includes a convergence proof of this approximation for a second-order ordinary differential equation for oscillations in the limit of small nonlinearity. We cannot perform this limit but justify the operation on the a priori observation that quadratic Galerkin system terms $q_{ijk}$ are typically two orders of magnitude smaller than the dominant linear coefficients, i.e. describe a small nonlinearity. A posteriori the operation is justified by the results, i.e. by obtaining amplitudes and frequencies with up to a few per cent error.

Substituting (3.21) in (3.20) replaces (3.24) by the mean-field manifold

(3.23)$$\begin{eqnarray}a_{3}=\unicode[STIX]{x1D705}\,(a_{1}^{2}+a_{2}^{2}),\end{eqnarray}$$

with derivable proportionality constant $\unicode[STIX]{x1D705}$.

Alternatively, the mean-field manifold may be obtained from the Galerkin system. Substituting (3.21) in (3.19) and projecting on $\boldsymbol{u}_{3}$ yields

(3.24)$$\begin{eqnarray}\text{d}a_{3}/\text{d}t=\unicode[STIX]{x1D70E}_{3}a_{3}+\unicode[STIX]{x1D6FD}_{3}(a_{1}^{2}+a_{2}^{2}),\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{3}<0$ and $\unicode[STIX]{x1D6FD}_{3}>0$ are necessary for a globally stable limit cycle. Note that (3.24) can be rewritten as

(3.25)$$\begin{eqnarray}\text{d}a_{3}/\text{d}t=\unicode[STIX]{x1D70E}_{3}[a_{3}-\unicode[STIX]{x1D705}(a_{1}^{2}+a_{2}^{2})].\end{eqnarray}$$

Now, the slaving process which leads to the mean-field manifold of (3.23) can be appreciated from the mean-field Galerkin system. If $|\unicode[STIX]{x1D70E}_{3}|\gg \unicode[STIX]{x1D70E}_{1}$, the time scale of slaving $a_{3}$ to the fluctuation level $a_{1}^{2}+a_{2}^{2}$ is much smaller than the time scale of the transient and $\text{d}a_{3}/\text{d}t$ can be set to zero.

Equations (3.22) and (3.23) yield the famous Landau equations with the cubic damping term

(3.26a,b)$$\begin{eqnarray}\text{d}r/\text{d}t=\unicode[STIX]{x1D70E}_{1}r-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D705}r^{3},\quad \text{d}\unicode[STIX]{x1D703}/\text{d}t=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}r^{2}.\end{eqnarray}$$

The Landau oscillator leads to a stable limit cycle with $r^{\circ }=\sqrt{\unicode[STIX]{x1D70E}_{1}/\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D705}}$, frequency $\unicode[STIX]{x1D714}^{\circ }=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D70E}_{1}\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D6FD}$ and shift-mode amplitude $a_{3}^{\circ }=\unicode[STIX]{x1D70E}_{1}/\unicode[STIX]{x1D6FD}$. The three nonlinearity parameters $\unicode[STIX]{x1D6FD}$, $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D705}$ can be uniquely derived from the limit cycle parameters $r^{\circ }$, $\unicode[STIX]{x1D714}^{\circ }$ and $a_{3}^{\circ }$. The growth rate $\unicode[STIX]{x1D70E}_{3}$ needs to be chosen sufficiently large, e.g. $\unicode[STIX]{x1D70E}_{3}=-10\unicode[STIX]{x1D70E}_{1}$, to ensure slaving on the manifold.

Equations (3.22), (3.24) are the mean-field Galerkin system, as illustrated in table 1, while equations (3.22), (3.23) characterize the original mean-field model, i.e. the slaved Galerkin system. Near the Hopf bifurcation, when $\unicode[STIX]{x1D70E}_{1}(Re)$ crosses the zero line at $Re_{1}$, the growth rate is approximated by $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D6FC}(Re-Re_{1})$ implying the square-root law $r^{\circ }=\sqrt{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D705}}\sqrt{Re-Re_{1}}$.

The Landau equation has been proposed by Landau (see, e.g. Landau & Lifshitz Reference Landau and Lifshitz1987), derived from the Navier–Stokes equation by Stuart (Reference Stuart1958), generalized for Galerkin systems by Noack et al. (Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003) and validated in numerous simulations and experiments for cylinder wakes (Schumm et al. Reference Schumm, Berger and Monkewitz1994) and other soft onsets of oscillatory flow. We note that the proposed derivation from symmetry considerations constrains the model to symmetric obstacles but liberates the mean-field Galerkin model from typical assumptions, such as closeness to the Hopf bifurcation or the need for frequency filtered Navier–Stokes equations.

3.4 Supercritical pitchfork bifurcation

Table 1.

Symmetries and hierarchy of equations.

Next, the symmetry-breaking pitchfork bifurcation of a steady symmetric Navier–Stokes equation is considered. Now, mode $\boldsymbol{u}_{4}$ describes the antisymmetric instability with positive growth rate $\unicode[STIX]{x1D70E}_{4}$. The shift mode $\boldsymbol{u}_{5}$ prevents unbounded exponential growth. The corresponding Galerkin expansion reads

(3.27)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}_{s}(\boldsymbol{x})+\underbrace{a_{4}(t)\boldsymbol{u}_{4}(\boldsymbol{x})}_{\boldsymbol{u}^{\prime }}+\underbrace{a_{5}(t)\boldsymbol{u}_{5}(\boldsymbol{x})}_{\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}}.\end{eqnarray}$$

Substituting (3.27) in equations (3.17) and (3.19), and exploiting the symmetry of the modes, yields

(3.28a)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{4}/\text{d}t=\unicode[STIX]{x1D70E}_{4}a_{4}-\unicode[STIX]{x1D6FD}_{4}a_{4}a_{5}, & \displaystyle\end{eqnarray}$$

(3.28b)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{5}/\text{d}t=\unicode[STIX]{x1D70E}_{5}a_{5}+\unicode[STIX]{x1D6FD}_{5}a_{4}^{2}. & \displaystyle\end{eqnarray}$$

$a_{5}$

$a_{4}a_{4}$

$a_{5}a_{5}$

$a_{4}$

$a_{4}a_{5}$

$a_{5}a_{5}$

$\unicode[STIX]{x1D714}=0$

$a_{1}$

$a_{4}$

$a_{3}$

$a_{5}$

$a_{2}=0$

Equations (3.28a), (3.28b) are the mean-field Galerkin system, as shown in table 1. Substituting (3.27) in (3.20) yields the manifold

(3.29)$$\begin{eqnarray}a_{5}=\unicode[STIX]{x1D705}_{5}a_{4}^{2},\end{eqnarray}$$

with $\unicode[STIX]{x1D705}_{5}=-\unicode[STIX]{x1D6FD}_{5}/\unicode[STIX]{x1D70E}_{5}$. The asymmetric steady solutions read $a_{4}^{\pm }=\pm \sqrt{\unicode[STIX]{x1D70E}_{4}/\unicode[STIX]{x1D6FD}_{4}\unicode[STIX]{x1D705}_{5}}$, $a_{5}=\unicode[STIX]{x1D70E}_{4}/\unicode[STIX]{x1D6FD}_{4}$. The two nonlinearity parameters $\unicode[STIX]{x1D705}_{4}$, $\unicode[STIX]{x1D6FD}_{4}$ are readily determined from the two asymptotic values $a_{4}$ and $a_{5}$. The growth rate can be set in analogy to the previous model to $\unicode[STIX]{x1D70E}_{5}=-10\unicode[STIX]{x1D70E}_{4}$ to ensure slaving on the manifold.

From (3.28a) and (3.29), the famous unstable dynamics with cubic damping term is obtained

(3.30)$$\begin{eqnarray}\displaystyle \text{d}a_{4}/\text{d}t=\unicode[STIX]{x1D70E}_{4}a_{4}-\unicode[STIX]{x1D705}_{5}\unicode[STIX]{x1D6FD}_{4}a_{4}^{3}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{5}\unicode[STIX]{x1D6FD}_{4}>0$ for a supercritical bifurcation. Equations (3.28a), (3.28b) are the mean-field Galerkin system.

Near the secondary pitchfork bifurcation, $\unicode[STIX]{x1D70E}_{4}=\unicode[STIX]{x1D6FC}_{2}(Re-Re_{2})$ and $a_{4}\propto \sqrt{Re-Re_{2}}$. The parameters of the pitchfork Galerkin system can be derived from the eigenmode and the asymptotic state in complete analogy to § 3.3. The growth rate $\unicode[STIX]{x1D70E}_{5}=-10\unicode[STIX]{x1D70E}_{4}$ will ensure the slaving of (3.29). We emphasize that this pitchfork model is derived primarily from symmetry considerations and does not require closeness to the critical parameter.

3.5 Pitchfork bifurcation of periodic solution

In the final modelling effort, a low-dimensional model from a primary supercritical Hopf bifurcation at $Re=Re_{1}$ and a secondary supercritical pitchfork bifurcation at $Re=Re_{2}>Re_{1}$ are derived following the numerical observations of the fluidic pinball in § 2. For simplicity, closeness to the secondary bifurcation is assumed. For the same reason, the mean-field Galerkin system shall still describe the periodic solution. In this case, the generalized 5-mode mean-field expansion

(3.31)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{u}_{s}(\boldsymbol{x})+\mathop{\sum }_{i=1}^{5}a_{i}(t)\boldsymbol{u}_{i}(\boldsymbol{x})\end{eqnarray}$$

describes the flow where $a_{1},a_{2},a_{3}\in O(1)$ and $a_{4}\in O(\unicode[STIX]{x1D700})$ and $a_{5}\in O(\unicode[STIX]{x1D6FF})$, $\unicode[STIX]{x1D700}$, $\unicode[STIX]{x1D6FF}$ being smallness parameters associated with the pitchfork bifurcation. We project equations (3.17) and (3.19) on (3.31). The $O(1)$ terms encapsulate the original Hopf model while the low-pass filtered $O(\unicode[STIX]{x1D700},\unicode[STIX]{x1D6FF})$ terms yield the original pitchfork system. This yields the following generalized mean-field system:

(3.32a)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{1}/\text{d}t=\unicode[STIX]{x1D70E}a_{1}-\unicode[STIX]{x1D714}a_{2},\quad \unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}, & \displaystyle\end{eqnarray}$$

(3.32b)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{2}/\text{d}t=\unicode[STIX]{x1D70E}a_{2}+\unicode[STIX]{x1D714}a_{1},\quad \unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}a_{3}, & \displaystyle\end{eqnarray}$$

(3.32c)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{3}/\text{d}t=\unicode[STIX]{x1D70E}_{3}a_{3}+\unicode[STIX]{x1D6FD}_{3}(a_{1}^{2}+a_{2}^{2}), & \displaystyle\end{eqnarray}$$

(3.32d)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{4}/\text{d}t=\unicode[STIX]{x1D70E}_{4}a_{4}-\unicode[STIX]{x1D6FD}_{4}a_{4}a_{5}, & \displaystyle\end{eqnarray}$$

(3.32e)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{5}/\text{d}t=\unicode[STIX]{x1D70E}_{5}a_{5}+\unicode[STIX]{x1D6FD}_{5}a_{4}^{2}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D70E}_{1}$

$\unicode[STIX]{x1D714}_{1}$

$\unicode[STIX]{x1D70E}_{4}$

$\unicode[STIX]{x1D70E}_{3}=-10\unicode[STIX]{x1D70E}_{1}$

$\unicode[STIX]{x1D70E}_{5}=-10\unicode[STIX]{x1D70E}_{4}$

$\unicode[STIX]{x1D6FD}$

$\unicode[STIX]{x1D6FE}$

$\unicode[STIX]{x1D6FD}_{3}$

$\unicode[STIX]{x1D6FD}_{4}$

$\unicode[STIX]{x1D6FD}_{5}$

$r^{\circ }$

$\unicode[STIX]{x1D714}^{\circ }$

$a_{3}^{\circ }$

$a_{4}^{\pm }$

$a_{5}^{\pm }$

As the amplitude of the pitchfork bifurcation grows, the smallness argument does not hold and we get cross-terms, like $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}-\unicode[STIX]{x1D6FD}_{15}a_{5}$. We shall not pause to elaborate on the possible generalizations now, but will return to the topic in the results section.

3.6 Sparse Galerkin model from mean-field considerations

The mean-field Galerkin system (3.32) with decoupled Hopf and pitchfork dynamics can, by construction, only be expected to hold near the pitchfork bifurcation $Re\approx Re_{2}$. At higher Reynolds numbers $Re>Re_{2}$, cross-terms will appear, e.g. the growth rate $\unicode[STIX]{x1D70E}$ may also depend on the pitchfork-related shift-mode amplitude $a_{5}$. The most general Galerkin system (3.7) contains $5\times 5=25$ linear terms and $5\times 5\times 6/2=75$ quadratic terms.

The assumed symmetry of the modes excludes roughly half of these 100 coefficients. Let $\unicode[STIX]{x1D712}_{i}=0$ for symmetric mean flow modes $\boldsymbol{u}_{i}$, $i=3,5$ and $\unicode[STIX]{x1D712}_{i}=1$ for the antisymmetric fluctuation modes $\boldsymbol{u}_{i}$, $i=1,2,4$. The linear coefficients

(3.33)$$\begin{eqnarray}\displaystyle l_{ij}=\unicode[STIX]{x1D708}(\boldsymbol{u}_{i},\triangle \boldsymbol{u}_{j})_{\unicode[STIX]{x1D6FA}}-(\boldsymbol{u}_{i},\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{s}\otimes \boldsymbol{u}_{j})_{\unicode[STIX]{x1D6FA}}-(\boldsymbol{u}_{i},\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{j}\otimes \boldsymbol{u}_{s})_{\unicode[STIX]{x1D6FA}} & & \displaystyle\end{eqnarray}$$

can be shown to vanish if $\text{mod}\;(\unicode[STIX]{x1D712}_{i}+\unicode[STIX]{x1D712}_{j},2)=1$. In other words, the coefficients $l_{ij}$ vanish if the modes $\boldsymbol{u}_{i}$ and $\boldsymbol{u}_{j}$ have opposite symmetries. This excludes 12 of the 25 linear coefficients. Analogously, the quadratic coefficient $q_{ijk}$ can be shown to vanish if $\unicode[STIX]{x1D712}_{i}\not =\text{mod}\;(\unicode[STIX]{x1D712}_{j}+\unicode[STIX]{x1D712}_{k},2)$. In other words, the $q_{ijk}$ vanishes if the symmetry of $\boldsymbol{u}_{i}$ does not coincide with the symmetry of quadratic term $\boldsymbol{u}_{j}\otimes \boldsymbol{u}_{k}$. The quadratic term is symmetric if the modes $\boldsymbol{u}_{j}$ and $\boldsymbol{u}_{k}$ are both symmetric or both antisymmetric and is antisymmetric if both modes have opposite symmetries. In summary, $q_{ijk}$ vanishes if one or three modes are antisymmetric.

An additional sparsity of the coefficients arises from the temporal dynamics. Modes $\boldsymbol{u}_{i}$, $i=1,2$ have oscillatory behaviour with angular frequency $\unicode[STIX]{x1D714}$, while the other modes show first-order dynamics, i.e. relaxation to asymptotic values. We apply the Krylov–Bogoliubov approximation with oscillatory $a_{1},a_{2}$ and slow $a_{3},a_{4},a_{5}$ dynamics. Thus, for instance, $l_{41}a_{1}$ vanishes on a one-period average and should not contribute to $\text{d}a_{4}/\text{d}t$. The linear coefficient $l_{41}$ can hence be set to zero. Taking the quadratic terms for example, $a_{1}a_{3}$ generates a first harmonic. Hence, $q_{413}a_{1}a_{3}$ cannot contribute to $\text{d}a_{4}/\text{d}t$ but $q_{113}a_{1}a_{3}$ can contribute to the oscillatory behaviour of $\text{d}a_{1}/\text{d}t$. Similarly, $a_{1}a_{2}$ generates a second harmonic, $a_{1}a_{2}=r\cos \unicode[STIX]{x1D714}t\times r\sin \unicode[STIX]{x1D714}t=(1/2)r^{2}\sin 2\unicode[STIX]{x1D714}t$ does not have a steady contribution, so $q_{312}$ can be set to zero.

From symmetry and Krylov–Bogoliubov considerations, only nine coefficients contribute to the linear term: $l_{11}$, $l_{12}$, $l_{21}$, $l_{22}$, $l_{33}$, $l_{35}$, $l_{44}$, $l_{53}$, $l_{55}$. Note that the oscillator equations contain the $2\times 2$ block, while the shift-mode equations $i=3,5$ have cross-terms and the pitchfork amplitude dynamics $i=4$ has no cross-terms. Similarly, only 16 quadratic coefficients survive. The first eight coefficients $q_{113}$, $q_{115}$, $q_{123}$, $q_{125}$, $q_{213}$, $q_{215}$, $q_{223}$, $q_{225}$ are consistent with the Landau oscillator but with cross-terms to the pitchfork-related shift-mode amplitude, i.e. $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}-\unicode[STIX]{x1D6FD}_{15}a_{5}$ and $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}a_{3}+\unicode[STIX]{x1D6FE}_{15}a_{5}$, introducing $\unicode[STIX]{x1D6FD}_{15}$ and $\unicode[STIX]{x1D6FE}_{15}$ as new coefficients. The first and second shift-mode equations $i=3,5$ may contain six quadratic terms $q_{311}$, $q_{322}$, $q_{344}$, $q_{511}$, $q_{522}$, $q_{544}$ from the Reynolds stresses. The amplification of the pitchfork dynamics is affected by the shift modes via $q_{443}$, $q_{445}$.

4.1 Eigenspectra of the steady solution

The steady solution becomes unstable beyond $Re=Re_{1}$, as reported in § 2. A linear stability analysis indicates that one pair of complex-conjugated eigenmodes have a positive growth rate on the range $Re_{1}<Re<Re_{2}$, as shown in figure 7 for $Re=30$. These two leading eigenmodes are associated with vortical structures shed downstream in the wake, at the angular frequency $1/2$. As the instability grows, the distorted mean flow $\boldsymbol{u}^{D}=\boldsymbol{u}_{s}+\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$ changes, as expected by equations (3.2) and (3.18). The shift mode $\boldsymbol{u}_{3}$, involved in $\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}$ at $Re=30$, is shown in figure 8(c). In the permanent (time-periodic) flow regime, the distorted mean flow $\boldsymbol{u}^{D}$ eventually matches the asymptotic mean flow field $\overline{\boldsymbol{u}}$, and vortex shedding is well established with frequency $8.7\times 10^{-2}$.

Figure 7.

(a) Eigenspectrum resulting from the linear stability analysis of the steady solution $\boldsymbol{u}_{s}$, together with (b) the first three leading eigenmodes, at $Re=30$. Only the real part of the complex eigenmodes is shown. Red colour and solid contours are positive values of the vorticity, blue colour and dashed contours are negative values.

Figure 8.

First two leading POD modes $\boldsymbol{u}_{1,2}$ at $Re=30$ (a) & (b) and shift mode $\boldsymbol{u}_{3}$ (c). Red colour and solid contours are positive values of the vorticity, blue colour and dashed contours are negative values.

4.2 Reduced-order model (ROM) of the primary flow regime

As introduced in § 3, the Galerkin ansatz for the Hopf bifurcation reads

(4.1)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)\approx \boldsymbol{u}_{s}(\boldsymbol{x})+a_{1}(t)\boldsymbol{u}_{1}(\boldsymbol{x})+a_{2}(t)\boldsymbol{u}_{2}(\boldsymbol{x})+a_{3}(t)\boldsymbol{u}_{3}(\boldsymbol{x}).\end{eqnarray}$$

The von Kármán modes $\boldsymbol{u}_{1,2}$ could be chosen as the real and imaginary parts of the first eigenmode, respectively. This choice would make sense to describe the transient dynamics close to the steady solution. A better choice for describing the dynamics on the asymptotic limit cycle is to choose $\boldsymbol{u}_{1,2}$ as the first two modes of a proper orthogonal decomposition (POD) of the limit cycle data. The first two POD modes actually contribute almost 95 % of the total fluctuating kinetic energy at $Re=30$ and are clearly associated with the von Kármán street of shed vortices, as shown in figure 8(a,b). For the construction of the ROM, POD modes $\boldsymbol{u}_{1,2}$ are preferred to the two leading eigenmodes, because we demand an accurate representation of the asymptotic periodic dynamics. Following equations (3.22)–(3.23), the dynamical system resulting from the Galerkin projection of ansatz (4.1) on the Navier–Stokes equations, after Krylov–Bogoliubov simplifications, reads

(4.2)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{1}/\text{d}t=\unicode[STIX]{x1D70E}a_{1}-\unicode[STIX]{x1D714}a_{2}, & \displaystyle\end{eqnarray}$$

(4.3)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{2}/\text{d}t=\unicode[STIX]{x1D70E}a_{2}+\unicode[STIX]{x1D714}a_{1}, & \displaystyle\end{eqnarray}$$

(4.4)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{3}/\text{d}t=\unicode[STIX]{x1D70E}_{3}(a_{3}-\unicode[STIX]{x1D705}(a_{1}^{2}+a_{2}^{2})), & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}$ and $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}a_{3}$. The value of the coefficients at $Re=30$ for the resulting ROM are summarized in table 2. Note that all coefficients but $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D70E}_{3}$ are fixed by either the linear stability analysis or the asymptotic dynamics, see § 3. The coefficient $\unicode[STIX]{x1D70E}_{3}$ can be chosen arbitrarily large as $a_{3}$ is slaved to $a_{1},a_{2}$ (here we chose $\unicode[STIX]{x1D70E}_{3}=-10\unicode[STIX]{x1D70E}_{1}$), while $\unicode[STIX]{x1D6FE}$ had to be calibrated in order to better match the asymptotic angular frequency.

Table 2.

Coefficients of the reduced-order model at $Re=30$. See text for details.

The dynamics of both the fluidic pinball (solid blue curve) and the ROM (dashed red curve) are compared in the three-dimensional subspace spanned by $a_{1}$, $a_{2}$, $a_{3}$, see the top of figure 9. In figure 9 are also shown the individual time series of $a_{1}$ to $a_{3}$ for both the fluidic pinball and the ROM (same representation). As expected from the POD modes $\boldsymbol{u}_{1,2}$, the dynamics on the asymptotic (permanent) limit cycle is well described in amplitude $r$ and angular frequency $\unicode[STIX]{x1D714}$ by the ROM. Moreover, the ROM also captures the transient dynamics on the parabolic manifold $a_{3}\equiv \unicode[STIX]{x1D705}(a_{1}^{2}+a_{2}^{2})$. Henceforth, although all coefficients but one are fixed, the Galerkin system (4.1) is able to reproduce the most salient dynamical features of the flow in both the transient and the permanent regimes.

Figure 9.

(a) Three-dimensional state space spanned by $a_{1}$, $a_{2}$, $a_{3}$, at $Re=30$. From direct numerical simulations of the fluidic pinball (solid blue line) and from the mean-field reduced-order model (dashed red line). The initial condition, identical in both systems, starts close to the fixed point (steady solution) before evolving on the parabolic manifold $a_{3}=\unicode[STIX]{x1D705}(a_{1}^{2}+a_{2}^{2})$ toward the asymptotic limit cycle. (b) Corresponding time series for $a_{1}$ to $a_{3}$.

5.1 Eigenspectra of the steady solutions

As a result of the pitchfork bifurcation, there exist three steady solutions beyond $Re_{2}$: the symmetric steady solution $\boldsymbol{u}_{s}$, unstable to the periodic vortex shedding beyond $Re_{1}$, and two mirror-conjugated asymmetric steady solutions $\boldsymbol{u}_{s}^{\pm }$, also unstable to vortex shedding but only existing beyond $Re_{2}$, see figure 3.

The linear stability analysis of $\boldsymbol{u}_{s}$ reveals two pairs of complex-conjugated eigenmodes with positive growth rate, and one eigenmode of zero frequency, see figure 10(a). The steady eigenmode is antisymmetric and reflects the symmetry broken by the pitchfork bifurcation. It is clearly associated with the base-bleeding jet, with all its energy concentrated in the near field. The two pairs of complex-conjugated eigenmodes are each associated with von Kármán streets of shed vortices. Both pairs of complex eigenmodes are antisymmetric and have quite similar angular frequencies. A closer view of the second pair of complex eigenmodes indicates that its growth rate cancels when the real eigenmode crosses the zero axis. This indicates that the new oscillatory mode is intimately connected to the symmetry breaking occurring at $Re_{2}$. At $Re>Re_{2}$, this gives rise to the only stable limit cycle for the flow dynamics, while the limit cycle associated with the leading pair of complex eigenmodes has become unstable and can only be visited transiently in time.

The linear stability analysis of $\boldsymbol{u}_{s}^{-}$ (respectively $\boldsymbol{u}_{s}^{+}$) reveals two pairs of complex-conjugated eigenmodes with positive growth rate, centred on an asymmetric mean flow, see figure 10(b). All eigenmodes are asymmetric, a property inherited from the steady solution.

In the permanent regime, the mean flow field will inherit the symmetry of one of the three (unstable) steady solutions, depending on the details of the initial perturbation.

Figure 10.

Eigenspectrum (top) and real part of the eigenvectors (bottom) of the symmetry-preserving steady solution $\boldsymbol{u}_{s}$ (a), of the symmetry-breaking steady solution $\boldsymbol{u}_{s}^{-}$ (b), both at $Re=80$. The red colour and solid contours in the eigenvectors are positive values of the vorticity, blue colour and dashed contours are negative values.

5.2 Reduced-order model in the secondary flow regime

Figure 11.

Additional modes arising from the pitchfork bifurcation, at $Re=80$, (a) mode $\boldsymbol{u}_{4}$, (b) mode $\boldsymbol{u}_{5}$. Red colour and solid contours are positive values of the vorticity, blue colour and dashed contours are negative values. See text for details about the computation of these two modes.

As discussed in § 3, an ansatz of the flow state can now be written as

(5.1)$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)\approx \underbrace{\boldsymbol{u}_{s}(\boldsymbol{x})}_{}+\underbrace{a_{1}(t)\boldsymbol{u}_{1}(\boldsymbol{x})+a_{2}(t)\boldsymbol{u}_{2}(\boldsymbol{x})}_{\mathit{leading}\;\mathit{POD}\;\mathit{modes}\;\mathit{at}~Re=80}+\underbrace{a_{3}(t)\boldsymbol{u}_{3}(\boldsymbol{x})}_{\mathit{shift}\;\mathit{mode}}+\underbrace{a_{4}(t)\boldsymbol{u}_{4}(\boldsymbol{x})+a_{5}(t)\boldsymbol{u}_{5}(\boldsymbol{x})}_{\mathit{pitchfork}\;\mathit{degrees}\;\mathit{of}\;\mathit{freedom}}.\end{eqnarray}$$

It is worthwhile noticing that, in the frame of this ansatz, the two asymmetric steady solutions $\boldsymbol{u}_{s}^{\pm }$ are related to the symmetric steady solution $\boldsymbol{u}_{s}$ via the additional antisymmetric mode $\boldsymbol{u}_{4}$

(5.2)$$\begin{eqnarray}\boldsymbol{u}_{s}^{\pm }=\boldsymbol{u}_{s}\pm \overline{a}_{4}\boldsymbol{u}_{4}+\overline{a}_{5}\boldsymbol{u}_{5},\end{eqnarray}$$

where $\overline{a}_{4}$ and $\overline{a}_{5}$ are the time-averaged coefficients in the permanent regime. Consequently, $\boldsymbol{u}_{4}$ can be easily computed as

(5.3)$$\begin{eqnarray}\boldsymbol{u}_{4}\propto (\boldsymbol{u}_{s}^{+}-\boldsymbol{u}_{s}^{-}),\end{eqnarray}$$

and further orthonormalized to $u_{1}$, $u_{2}$, $u_{3}$ by a Gram–Schmidt procedure. The resulting mode $\boldsymbol{u}_{4}$ is shown in figure 11(a). A comparison with the eigenmode associated with the real eigenvalue, in figure 10(a), shows that the shift mode $\boldsymbol{u}_{4}$ is just the real eigenmode against which the symmetric steady solution is unstable at $Re=80$, as expected by the definition of mode $\boldsymbol{u}_{4}$.

In a similar way, the additional mode $\boldsymbol{u}_{5}$ can be constructed as

(5.4)$$\begin{eqnarray}\boldsymbol{u}_{5}\propto (\boldsymbol{u}_{s}^{+}+\boldsymbol{u}_{s}^{-})/2-\boldsymbol{u}_{s}.\end{eqnarray}$$

Mode $\boldsymbol{u}_{5}$ is shown in figure 11(b) after orthonormalization.

Table 3.

Coefficients of the least reduced-order model (3.32) at $Re=80$.

Figure 12.

Time evolution of coefficients $a_{1}$ to $a_{5}$ in the full flow dynamics (solid blue line) and the ROM (dashed red line) without cross-terms and coefficients fixed by the linear stability analysis and the asymptotic dynamics. The initial condition for both systems is the same.

Close to the pitchfork bifurcation, the resulting dynamical system is described by (3.32). At the threshold, the degrees of freedom $a_{4}$, $a_{5}$ associated with the pitchfork bifurcation are expected to be fully uncoupled to the degrees of freedom $a_{1}$, $a_{2}$, $a_{3}$ associated with the Hopf bifurcation, and reciprocally, see § 3.3. In this case, an accurate linear and nonlinear dynamics from a Galerkin projection relies on deformable modes from eigenmodes near the fixed point to POD modes near the limit cycle (Loiseau, Noack & Brunton Reference Loiseau, Noack and Brunton2018). We avoid this complication by model identification from simulation data. The coefficients of the mean-field system (3.32) reported in table 3 were identified from the linear stability analysis of the symmetric steady solution and the asymptotic dynamics on the unstable symmetric preserving limit cycle and the stable symmetric-breaking limit cycle (see §§ 3.3 to 3.5). The resulting ROM dynamics is compared to the flow dynamics in figure 12. Inspection of figure 12 shows that such a model, reduced to only five degrees of freedom, is able to reproduce many features of the original dynamics: both the early transient and asymptotic dynamics of $a_{1}$ to $a_{3}$ are well reproduced, as well as the large time scale evolution of $a_{4},a_{5}$. However, the growth of coefficients $a_{1}$ to $a_{3}$ appears to be faster for the ROM, and $a_{4}$, $a_{5}$ also reach their asymptotic values significantly sooner than their values for the full simulations of the fluidic pinball. In addition, the transient kick in $a_{4},a_{5}$ is absent from the mean-field system (3.32), as well as the transient and asymptotic oscillations of $a_{4}$ and $a_{5}$ visible in figure 12, which would require coupling to $a_{1}$ or $a_{2}$, or both. All these features indicate that at Reynolds number $Re=80$, the Krylov–Bogoliubov assumption of pure harmonic behaviour with slowly varying amplitude and frequency no longer holds.

As a consequence, cross-terms must be included in the ROM. The assumption of non-oscillatory dynamics of the shift-mode amplitude $a_{3}$ and of the two pitchfork modes $a_{4}$, $a_{5}$ is relaxed to reproduce the oscillatory behaviour evidenced in figure 12. Following § 3.6, the model identification process reads as follows.

Step 1:

keep the five-dimensional linear–quadratic form of the dynamical system from the Galerkin projection with 25 linear and 75 quadratic terms.

Step 2:

remove vanishing terms arising from the symmetry of the modes. Thus, only 13 linear and 36 quadratic terms are left to be determined.

Step 3:

enforce the linear dynamics of the unstable Hopf and pitchfork eigenmodes from stability analysis in the Galerkin system. This implies that the growth rate $\unicode[STIX]{x1D70E}_{1}$ and frequency $\unicode[STIX]{x1D714}_{1}$ characterize the initial growth and angular frequency of $a_{1}$, $a_{2}$ and $\unicode[STIX]{x1D70E}_{4}$ represents the linear growth rate for $a_{4}$.

Step 4:

the slaving of the Reynolds-stress-induced modes $\boldsymbol{u}_{3}$ and $\boldsymbol{u}_{5}$ to the fluctuation level is imposed by setting the damping rate 10 times larger than the growth rate of the corresponding fluctuation: $\unicode[STIX]{x1D70E}_{3}=-10\unicode[STIX]{x1D70E}_{1}$, $\unicode[STIX]{x1D70E}_{5}=-10\unicode[STIX]{x1D70E}_{4}$. This strong damping rate quickly forces the trajectory onto the mean-field manifold.

Step 5:

enforce phase invariance for $a_{1}$, $a_{2}$ in the first two equations. This is implied by the mean-field theory and is found to be a good approximation from numerical inspection. Thus, the oscillatory dynamics of $a_{1}$, $a_{2}$ is governed by $\unicode[STIX]{x1D70E}_{1}$, $\unicode[STIX]{x1D714}_{1}$, $\unicode[STIX]{x1D6FD}$, $\unicode[STIX]{x1D6FD}_{15}$, $\unicode[STIX]{x1D6FE}$, $\unicode[STIX]{x1D6FE}_{15}$ and $\unicode[STIX]{x1D6FD}_{3}$.

Step 6:

impose the asymptotic dynamics of the unstable symmetric limit cycle by fixing $\unicode[STIX]{x1D6FD}$, $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FD}_{3}$.

Step 7:

apply the SINDy algorithm (Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016) to the remaining unknown terms, i.e. 6 linear terms and 28 quadratic terms. This step alone typically fails to yield a physics-based globally stable Galerkin system. This is not surprising in view of the necessary and sufficient conditions for global boundedness of Galerkin systems by Schlegel & Noack (Reference Schlegel and Noack2015). A Galerkin system identification with $\ell 1$-norm penalization of the coefficients is found to have a performance similar to SINDy for the chosen constraints.

Step 8:

two additional simplifying physics-based constraints are found to make the dynamic system identified by SINDy robust for a large range of initial conditions. Enforcing $q_{145}=q_{245}=0$ avoids the initial oscillatory dynamics being influenced by the initial kick of $a_{5}$, and setting $q_{344}=q_{434}=0$ leads to the right asymptotic transition of the pitchfork bifurcation.

Figure 13.

Performance of the ROM with cross-terms. Time evolution of coefficients $a_{1}$ to $a_{5}$ in the full flow dynamics (solid blue line) and for the ROM (red dashed line). The initial condition is the same for the ROM and the full flow dynamics.

The resulting ROM reads

(5.5a)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{1}/\text{d}t=a_{1}(\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}-\unicode[STIX]{x1D6FD}_{15}a_{5})-a_{2}(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}a_{3}+\unicode[STIX]{x1D6FE}_{15}a_{5})+l_{14}a_{4}+q_{134}a_{3}a_{4}, & \displaystyle\end{eqnarray}$$

(5.5b)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{2}/\text{d}t=a_{2}(\unicode[STIX]{x1D70E}_{1}-\unicode[STIX]{x1D6FD}a_{3}-\unicode[STIX]{x1D6FD}_{15}a_{5})+a_{1}(\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D6FE}a_{3}+\unicode[STIX]{x1D6FE}_{15}a_{5})+l_{24}a_{4}+q_{234}a_{3}a_{4}, & \displaystyle\end{eqnarray}$$

(5.5c)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{3}/\text{d}t=\unicode[STIX]{x1D70E}_{3}a_{3}+\unicode[STIX]{x1D6FD}_{3}r^{2}+l_{35}a_{5}+q_{314}a_{1}a_{4}+q_{335}a_{3}a_{5}+q_{355}a_{5}^{2}, & \displaystyle\end{eqnarray}$$

(5.5d)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{4}/\text{d}t=\unicode[STIX]{x1D70E}_{4}a_{4}-\unicode[STIX]{x1D6FD}_{4}a_{4}a_{5}+a_{1}(l_{41}+q_{413}a_{3}+q_{415}a_{5})+a_{2}(l_{42}+q_{423}a_{3}+q_{425}a_{5}),\qquad & \displaystyle\end{eqnarray}$$

(5.5e)$$\begin{eqnarray}\displaystyle & \displaystyle \text{d}a_{5}/\text{d}t=\unicode[STIX]{x1D70E}_{5}a_{5}+\unicode[STIX]{x1D6FD}_{5}a_{4}^{2}+l_{53}a_{3}+q_{514}a_{1}a_{4}+q_{533}a_{3}^{2}+q_{535}a_{3}a_{5}, & \displaystyle\end{eqnarray}$$

$r^{2}=a_{1}^{2}+a_{2}^{2}$

$a_{3}$

$a_{4}$

$a_{5}$

$a_{1}$

$a_{3}$

$a_{4}$

$a_{4}\neq 0$

$\unicode[STIX]{x1D70E}_{1}=5.22\times 10^{-2}$

$\boldsymbol{u}_{s}$

$a_{4}$

$a_{3}$

$a_{4}$

$a_{5}$

$a_{3}$

$\unicode[STIX]{x1D70E}_{3}$

Model identification can be very challenging for a number of reasons. First, the conditions for global boundedness of the attractor for the linear–quadratic Galerkin system are rather restrictive (Schlegel & Noack Reference Schlegel and Noack2015). Second, the Galerkin method assumes fixed expansion modes. Yet, least-order models often have deformable modes changing with the fluctuation level (Tadmor et al. Reference Tadmor, Lehmann, Noack, Cordier, Delville, Bonnet and Morzyński2011). The von Kármán vortex shedding metamorphosis from stability modes to POD modes may serve as an example. These deformations may also affect the structure of the dynamical system. Third, the Navier–Stokes dynamics may live on a strongly attracting manifold. This restriction of the state space may make certain Galerkin system coefficients numerically unobservable, like $\unicode[STIX]{x1D70E}_{3}$ and $\unicode[STIX]{x1D70E}_{5}$ in our case. Despite these challenges, the model identified in table 4, as exemplified by figure 13, constitutes a faithful least-order model for the fluidic pinball at $Re=80$. Although the model is only five-dimensional, it can reproduce most of the key features, time scales, transient and asymptotic behaviours of the full dynamics.

Table 4.

Coefficients of the reduced-order model at $Re=80$. See text for details.

6.1 Concluding remarks and discussion

Flow configurations undergoing successive Hopf and pitchfork bifurcations are common in fluid mechanics. This is, for instance, the case of three-dimensional wake flows such as spheres (Mittal Reference Mittal1999; Gumowski et al. Reference Gumowski, Miedzik, Goujon-Durand, Jenffer and Wesfreid2008; Szaltys et al. Reference Szaltys, Chrust, Przadka, Goujon-Durand, Tuckerman and Wesfreid2012; Grandemange et al. Reference Grandemange, Gohlke and Cadot2014) or bluff body wake flows (Grandemange et al. Reference Grandemange, Cadot and Gohlke2012, Reference Grandemange, Gohlke and Cadot2013; Cadot et al. Reference Cadot, Evrard and Pastur2015; Bonnavion & Cadot Reference Bonnavion and Cadot2018; Rigas et al. Reference Rigas, Oxlade, Morgans and Morrison2014). The drag crisis and stalled flows are also characterized by the pitchfork bifurcation of a primarily Hopf-bifurcated flow, but the secondary transition is subcritical in this case.

In this study, we have considered the fluidic pinball on its way to chaos and have identified least-order models of the flow dynamics in the primary Hopf-bifurcated and secondary pitchfork-bifurcated flow regimes. Reduced-order modelling of Hopf bifurcations was already addressed in Noack et al. (Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003) for the cylinder wake flow, while Meliga et al. (Reference Meliga, Chomaz and Sipp2009) derived the amplitude equation for a codimension two bifurcation (pitchfork and Hopf) based on a weakly nonlinear analysis in the wake of a disk, and Fabre et al. (Reference Fabre, Auguste and Magnaudet2008) derived the same equation solely based on symmetry arguments, in the wake of axisymmetric bodies. In the present contribution, we could demonstrate that the dynamics resulting from the successive Hopf and pitchfork bifurcations could be well captured by a five-dimensional model whose degrees of freedom couple through quadratic nonlinearities, as imposed by the Navier–Stokes equations.

For the fluidic pinball, the route to chaos is characterized by a primary supercritical Hopf bifurcation at $Re\approx 18$, followed by a secondary supercritical pitchfork bifurcation at $Re\approx 68$. The Hopf bifurcation corresponds to the destabilization of the steady solution with respect to vortex shedding, while the pitchfork bifurcation occurs when the mean flow breaks the symmetry with respect to the mirror plane. The fluctuation amplitude of the von Kármán street is reduced, over a finite range of the Reynolds number around $Re_{2}$, when the base-bleeding jet is rising. This means that energy is withdrawn from the fluctuations to feed the mean flow transformation. Before the next transition occurs, the fluctuation amplitude starts to grow again, up to the largest value of the Reynolds number considered in this work.

There is strong evidence that the next transition is a Neimark–Säcker bifurcation. The resulting flow regime is most likely quasi-periodic over the range $[Re_{3},Re_{4}]$. In this regime, a new oscillatory phenomenon takes place, characterized by slow oscillations of the base-bleeding jet. Three additional degrees of freedom might be necessary to deal with the newly arising oscillator. The flow dynamics eventually bifurcates into a chaotic regime, characterized by the random switching of the base-bleeding jet between two symmetric deflected positions. The overall route to chaos is summarized in the phenomenogram of figure 14.

Figure 14.

Phenomenogram. The route to chaos develops along the $Re_{D}$-axis. On this route are highlighted the least-order models identified at $Re=30$ and $Re=80$, where the manifolds on which the dynamics takes place are schematically drawn. The branches of the steady solution and periodic solution are presented by the black and blue curves, with the solid/dashed curves for the stable/unstable states. The degrees of freedom $a_{1}$, $a_{2}$ span the limit cycle subspace, $a_{3}$ is the axis of the parabolically shaped manifold, while $a_{4}$ is transversally associated with the pitchfork bifurcation, together with $a_{5}$, which slightly bends the surface to which the steady solutions $\boldsymbol{u}_{s}^{\pm }$ and $\boldsymbol{u}_{s}$ belong. The grey shaded shadows of the quasi-periodic and chaotic regimes are represented for the sake of illustration. Grey shaded inserts provide the least reduced-order models at the threshold for both the Hopf and pitchfork bifurcations under the constraint of the Navier–Stokes equations. Also shown at the top are the figure snapshots of the steady solutions, mean flow fields and instantaneous flow field, at either $Re=30$ or $Re=80$. The dotted arrows connect these snapshots to the corresponding points in the phenomenogram. In addition, individual degrees of freedom $\boldsymbol{u}_{1}$ to $\boldsymbol{u}_{5}$ are shown at the left.

The reduced-order models derived by Galerkin projections of the Navier–Stokes equations, based on the symmetry of the individual degrees of freedom, under Krylov–Bogoliubov simplifications, faithfully extract the manifolds on which the flow dynamics sets in. The ROM for the primary flow regime is only three-dimensional: two degrees of freedom are associated with the asymptotic stable limit cycle resulting from the vortex shedding. The third degree of freedom is a mode slaved to the two dominant modes and is mandatory for the description of the transient flow dynamics from the unstable steady solution to the post-transient mean flow, as already demonstrated in Noack et al. (Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003). The least-order model in the secondary flow regime has only five degrees of freedom, three of which are associated with the Hopf bifurcation, the two remaining degrees of freedom being associated with the pitchfork bifurcation. In the phenomenogram of figure 14 are reported the structure of both reduced-order models close to the Hopf and the pitchfork bifurcations. When the two sets of degrees of freedom are fully uncoupled, some features of the flow dynamics are well reproduced (asymptotic mean behaviour, parabolic manifolds of the Hopf and pitchfork bifurcations), but many details are missing. To reproduce most of the transient and asymptotic flow features far from the bifurcation point, additional cross-terms have been included in the model, which relaxes the steadiness constraint usually assumed for the shift modes.

6.2 Outlook

The current generalized mean-field model captures the Hopf bifurcation and subsequent pitchfork bifurcation of the steady solution and limit cycles. The following onset of a quasi-periodic regime with slow oscillations of the deflected base-bleeding jet might presumably be incorporated by another Hopf bifurcation, leading to an eight-dimensional mean-field Galerkin model. The transition to chaos is accompanied by a return to a statistically symmetric flow, i.e. the base-bleeding jet oscillates around one asymmetric state before it stochastically switches to the other mirror-symmetric one. This behaviour is reminiscent of the transition to chaos of a harmonically forced Duffing oscillator. In the case of the fluidic pinball, the vortex shedding would constitute a forcing. Hence, one may speculate that the transition to chaos may already be resolved by the eight-dimensional Galerkin model in which the effect of vortex shedding on the jet oscillation becomes stronger with an increasing Reynolds number.

An alternative direction is to increase the accuracy of the mean-field Galerkin model. While the structure of the Galerkin system prevails for a large range of Reynolds numbers, the modes and all Galerkin system coefficients change in a non-trivial manner, e.g. the growth-rate formula should read $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{1}(Re)-\unicode[STIX]{x1D6FD}(Re)a_{3}$. The transients can be expected to be much more accurately resolved by the mean flow dependent modes, e.g. $\boldsymbol{u}_{i}(Re,a_{3},a_{5},\boldsymbol{x})$, $i=1,2$ for the resolution of vortex shedding (Loiseau et al. Reference Loiseau, Noack and Brunton2018). More generally, the flow lives on a low-dimensional manifold which includes mode deformations (Noack Reference Noack2016). Locally linear embedding (LLE) is a powerful technique for identifying the dimension and a parameterization in an automatic manner (Roweis & Saul Reference Roweis and Saul2000). The normal form of the bifurcations can be expected to coincide with the dynamics on the LLE feature coordinates.

A third direction follows an observation of Rempfer (Reference Rempfer1994) that Galerkin systems of many fluid flows can be considered as nonlinearly coupled oscillators. For two incommensurable shedding frequencies, this observation has been formalized in a generalized mean-field model by Noack et al. (Reference Noack, Schlegel, Ahlborn, Mutschke, Morzyński, Comte and Tadmor2008), Luchtenburg et al. (Reference Luchtenburg, Günter, Noack, King and Tadmor2009). Such multi-frequency models may be extended to resolve broadband frequency dynamics taking, for instance, the most dominant dynamic mode decomposition (DMD) modes (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Schmid Reference Schmid2010). Strengths and weaknesses of techniques currently used for model reduction are discussed in Taira et al. (Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). While mean-field consideration expressly ignores non-trivial triadic interactions, their quantitative effect on the frequency cross-talk may still be well approximated by the mean flow interaction terms. Such multi-frequency mean-field models may eventually describe the effect of open-loop forcing on turbulence, see for instance the recent thorough review by Jiménez (Reference Jiménez2018) on turbulent flow modelling.

A fourth direction aligned with the large success of machine learning/artificial intelligence is the automated learning of state spaces, modes and dynamical systems. For the latter, SINDy provides an established elegant framework (Brunton et al. Reference Brunton, Proctor and Kutz2016). The choice of the state spaces might be facilitated by manifold learning from many solution snapshots (Gorban & Karlin Reference Gorban and Karlin2005). The authors actively pursue all the mentioned directions.