Low-order model for successive bifurcations of the fluidic pinball

TL;DR

This paper introduces a minimal 5-dimensional Galerkin model that captures two successive supercritical bifurcations in a complex wake flow, advancing understanding of fluid flow transitions.

Contribution

It presents the first low-order Galerkin model incorporating symmetry considerations to describe successive bifurcations in fluid dynamics.

Findings

Model accurately reproduces bifurcation sequence in fluidic pinball

Successfully describes transition to chaos through bifurcations

Generalized mean-field approach applicable to other flow scenarios

Abstract

We propose the first least-order Galerkin model of an incompressible flow undergoing two successive supercritical bifurcations of Hopf and pitchfork type. A key enabler is a mean-field consideration exploiting the symmetry of the mean flow and the asymmetry of the fluctuation. These symmetries generalize mean-field theory, e.g. no assumption of slow growth-rate is needed. The resulting 5-dimensional Galerkin model successfully describes the phenomenogram of the fluidic pinball, a two-dimensional wake flow around a cluster of three equidistantly spaced cylinders. The corresponding transition scenario is shown to undergo two successive supercritical bifurcations, namely a Hopf and a pitchfork bifurcations on the way to chaos. The generalized mean-field Galerkin methodology may be employed to describe other transition scenarios.

Figures

Tables

Table 1. Table 2: Coefficients of the reduced-order model (ROM) at R ​ e = 30 𝑅 𝑒 30 Re=30 . See text for details.

${\sigma }_1$	${\omega }_1$	${\sigma }_3$	$\beta$	$\gamma$	$\kappa$
$3.80\times {10}^{-2}$	$5.00\times {10}^{-1}$	$-10{\sigma }_1$	$1.40\times {10}^{-2}$	$1.70\times {10}^{-2}$	$2.10\times {10}^{-1}$

Table 2. Table 3: Coefficients of the least reduced-order model ( 32 ) at R ​ e = 80 𝑅 𝑒 80 Re=80 .

${\sigma }_1$	$5.22\times {10}^{-2}$	$\beta$	$1.31\times {10}^{-2}$
${\omega }_1$	$5.24\times {10}^{-1}$	$\gamma$	$2.95\times {10}^{-2}$
${\sigma }_3$	$-5.22\times {10}^{-1}$	${\beta }_3$	$1.53\times {10}^{-1}$
${\sigma }_4$	$2.72\times {10}^{-2}$	${\beta }_4$	$2.45\times {10}^{-1}$
${\sigma }_5$	$-2.72\times {10}^{-1}$	${\beta }_5$	$2.14\times {10}^{-1}$

Table 3. Table 4: Coefficients of the reduced-order model at R ​ e = 80 𝑅 𝑒 80 Re=80 . See text for details.

${\sigma }_1$	$5.22\times {10}^{-2}$	$\beta$	$1.31\times {10}^{-2}$	$l_{14}$	$2.93\times {10}^{-1}$	$l_{24}$	$-4.87\times {10}^{-1}$
${\omega }_1$	$5.24\times {10}^{-1}$	$\gamma$	$2.95\times {10}^{-2}$	$q_{134}$	$-5.87\times {10}^{-2}$	$q_{234}$	$1.18\times {10}^{-1}$
${\sigma }_3$	$-5.22\times {10}^{-1}$	${\beta }_3$	$1.53\times {10}^{-1}$	$l_{41}$	$3.14\times {10}^{-2}$	$l_{42}$	$-5.14\times {10}^{-2}$
${\sigma }_4$	$2.72\times {10}^{-2}$	${\beta }_4$	$5.78\times {10}^{-2}$	$q_{413}$	$-7.56\times {10}^{-3}$	$q_{423}$	$1.28\times {10}^{-2}$
${\sigma }_5$	$-2.72\times {10}^{-1}$	${\beta }_5$	$1.91\times {10}^{-1}$	$q_{415}$	$2.99\times {10}^{-2}$	$q_{425}$	$1.71\times {10}^{-1}$
		${\beta }_{15}$	$-2.42\times {10}^{-2}$	$l_{35}$	$4.28$	$l_{53}$	$2.89\times {10}^{-2}$
		${\gamma }_{15}$	$1.70\times {10}^{-2}$	$q_{335}$	$-1.11$	$q_{533}$	$-7.22\times {10}^{-3}$
				$q_{355}$	$-5.13\times {10}^{-1}$	$q_{535}$	$1.48\times {10}^{-2}$
				$q_{314}$	$1.57\times {10}^{-2}$	$q_{514}$	$-9.44\times {10}^{-3}$

Equations

$$U_r=0\; on\; the\; cylinders\; and{U}_{\infty }={\mathbf{e}}_{x}at\; infinity.$$

$${\partial }_t\mathbf{u}+\nabla \cdot \mathbf{u}\otimes \mathbf{u}=\nu △\mathbf{u}-\nabla p,$$

$$\nabla \cdot {\mathbf{u}}_s\otimes {\mathbf{u}}_s=\nu △{\mathbf{u}}_s-\nabla p_s.$$

$${\left(\mathbf{u},\mathbf{v}\right)}_{\Omega }:=\underset{\Omega }{\int }d\mathbf{x}\mathbf{u}(\mathbf{x})\cdot \mathbf{v}(\mathbf{x}).$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}_0(\mathbf{x})+\sum _{i=1}^N{a}_i(t){\mathbf{u}}_i(\mathbf{x}).$$

$${\left({\mathbf{u}}_i,{\mathbf{u}}_j\right)}_{\Omega }={\delta }_{ij},i,j\in \{1,\dots ,N\}.$$

$$\frac{d}{dt}{a}_i=\nu \sum _{j=0}^N{l}_{ij}^{\nu }{a}_j+\sum _{j,k=0}^N{q}_{ijk}^c{a}_j{a}_k.$$

$$\frac{d}{dt}{a}_i=\sum _{j=1}^N{l}_{ij}{a}_j+\sum _{j,k=1}^N{q}_{ijk}{a}_j{a}_k,$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}^D(\mathbf{x},t)+{\mathbf{u}}^{\prime }(\mathbf{x},t),{\mathbf{u}}^D(\mathbf{x},t)={\mathbf{u}}_s(\mathbf{x})+{\mathbf{u}}_{\Delta }(\mathbf{x},t)$$

$${\mathbf{u}}^D(\mathbf{x},t)=⟨\mathbf{u}(\mathbf{x},t)⟩:=\frac{1}{T}\underset{t-T/2}{\overset{t+T/2}{\int }}d\tau \mathbf{u}(\mathbf{x},\tau ).$$

$${\mathbf{u}}^D(\mathbf{x},t)=\frac{1}{2}(⟨{\mathbf{u}}^{+}(\mathbf{x},t)⟩+⟨{\mathbf{u}}^{-}(\mathbf{x},t)⟩).$$

$$\bar{\mathbf{u}}(\mathbf{x})=\underset{T\to \infty }{\lim }\frac{1}{T}\underset{0}{\overset{T}{\int }}\mathbf{u}(\mathbf{x},\tau )d\tau .$$

$$u^s(x,-y)=u^s(x,y),v^s(x,-y)=-v^s(x,y)$$

$$u^a(x,-y)=-u^a(x,y),v^a(x,-y)=v^a(x,y).$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}^s(\mathbf{x},t)+{\mathbf{u}}^a(\mathbf{x},t).$$

$${\mathbf{u}}^D(\mathbf{x},t)={\mathbf{u}}^s(\mathbf{x},t),{\mathbf{u}}^{\prime }(\mathbf{x},t)={\mathbf{u}}^a(\mathbf{x},t).$$

$$\nabla \cdot {\mathbf{u}}^s\otimes {\mathbf{u}}^s$$

$$\nabla \cdot {\mathbf{u}}^a\otimes {\mathbf{u}}^a$$

$$\nabla \cdot {\mathbf{u}}^s\otimes {\mathbf{u}}^a$$

$$\nabla \cdot {\mathbf{u}}^a\otimes {\mathbf{u}}^s$$

$${\partial }_t{\mathbf{u}}^{\prime }+\nabla \cdot \left[{\mathbf{u}}^D\otimes {\mathbf{u}}^{\prime }+{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^D\right]=\nu △{\mathbf{u}}^{\prime }-\nabla p^{\prime }.$$

$${\partial }_t{\mathbf{u}}_{\Delta }+\nabla \cdot \left[{\mathbf{u}}_s\otimes {\mathbf{u}}_{\Delta }+{\mathbf{u}}_{\Delta }\otimes {\mathbf{u}}_s+{\mathbf{u}}_{\Delta }\otimes {\mathbf{u}}_{\Delta }+{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }\right]=\nu △{\mathbf{u}}_{\Delta }-\nabla p_{\Delta }$$

$${\partial }_t{\mathbf{u}}_{\Delta }+\nabla \cdot \left[{\mathbf{u}}_s\otimes {\mathbf{u}}_{\Delta }+{\mathbf{u}}_{\Delta }\otimes {\mathbf{u}}_s+⟨{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }⟩\right]=\nu △{\mathbf{u}}_{\Delta }-\nabla p_{\Delta }.$$

$$\nabla \cdot \left[{\mathbf{u}}_s\otimes {\mathbf{u}}_{\Delta }+{\mathbf{u}}_{\Delta }\otimes {\mathbf{u}}_s+⟨{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }⟩\right]=\nu △{\mathbf{u}}_{\Delta }-\nabla p_{\Delta }.$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}_s(\mathbf{x})+\underset{{\mathbf{u}}^{\prime }}{\underbrace{a_1(t){\mathbf{u}}_1(\mathbf{x})+a_2(t){\mathbf{u}}_2(\mathbf{x})}}+\underset{{\mathbf{u}}_{\Delta }}{\underbrace{a_3(t){\mathbf{u}}_3(\mathbf{x})}}.$$

$$d{a}_1/dt$$

$$d{a}_2/dt$$

$$a_3=\kappa \left(a_1^2+a_2^2\right)$$

$$d{a}_3/dt={\sigma }_3{a}_3+{\beta }_3\left(a_1^2+a_2^2\right),$$

$$d{a}_3/dt={\sigma }_3\left[a_3-\kappa \left(a_1^2+a_2^2\right)\right].$$

$$dr/dt={\sigma }_{1}r-\beta \kappa r^3,d\theta /dt={\omega }_1+\gamma \kappa r^2.$$

$$\mathbf{u}(\mathbf{x},t)={\mathbf{u}}_s(\mathbf{x})+\underset{{\mathbf{u}}^{\prime }}{\underbrace{a_4(t){\mathbf{u}}_4(\mathbf{x})}}+\underset{{\mathbf{u}}_{\Delta }}{\underbrace{a_5(t){\mathbf{u}}_5(\mathbf{x})}}.$$

Full text

1 Introduction

This study advances mean-field modeling for successive symmetry breaking due to Hopf and pitchfork bifurcations. The theoretical framework is applied to the transition of the flow around a cluster of circular cylinders, termed the fluidic pinball for the possibility to control fluid particle by cylinder rotation (Noack & Morzyński, 2017; Ishar et al., 2019).

Mean-field theory was pioneered by Landau (1944) and Stuart (1958) and is a singular triumph of nonlinear reduced-order modeling in fluid mechanics. Already the most simple mean-field model for the supercritical Hopf bifurcation reveals deep insights into the coupling between the fluctuations and the mean flow, e.g. the damping mechanism of unstable modes by Reynolds stress. In addition, Malkus (1956) principle of marginal stability for time-averaged flows, the square root growth law of fluctuation level with increasing Reynolds number, the cubic damping term from a linear-quadratic dynamics, the energetic explanation of this amplitude dynamics, and the slaving principle leading to manifolds driven by ensemble-averaged Reynolds stress are easily derived. Also, the idea of center manifold theory and the surprising success of linear parameter-varying models are analytically illustrated. Historically, Landau was the first to derive the normal form of the dynamics with Krylov-Bogoliubov approximation (an averaging method for spiral phase paths, see e.g. Jordan & Smith, 1999) while Stuart could explain how the cubic damping term arises from the distorted mean flow.

Mean-field models for a supercritical Hopf bifurcation with an unstable oscillatory eigenmode have been applied and validated for numerous configurations. The onset of vortex shedding behind a cylinder wake has been thoroughly investigated (Strykowski & Sreenivasan, 1990; Schumm et al., 1994; Noack et al., 2003). Even high-Reynolds number turbulent wake flow can display a distinct mean-field manifold and modeled by a noise-driven mean-field model (Bourgeois et al., 2013).

A supercritical pitchfork bifurcation similarly arises by an unstable eigenmode with a real eigenvalue. The onset of convection rolls in the Rayleigh-Bénard problem is a famous example (Zaitsev & Shliomis, 1971; Swift & Hohenberg, 1977; Cross & Hohenberg, 1993). The features of a pitchfork bifurcation are observed for the sidewise symmetry breaking of the time-averaged Ahmed body wake (Grandemange et al., 2012, 2013; Cadot et al., 2015; Bonnavion & Cadot, 2018) and more generally in three-dimensional wake flows (Mittal, 1999; Gumowski et al., 2008; Szaltys et al., 2012; Grandemange et al., 2014; Rigas et al., 2014). In contrast, the drag crisis of circular cylinder is associated with a subcritical bifurcation into two asymmetric sheddings with opposite mean lift values (Schewe, 1983).

Not surprisingly, numerous generalizations of mean-field models have been proposed. Landau (1944) and Hopf (1948) have conjectured that high-dimensional fully developed turbulence may be explained by an increasingly rapid succession of Hopf bifurcations. This idea has been discarded as unlikely (see, for instance, Landau & Lifshitz, 1987). The seminal paper by Ruelle & Takens (1971) showed that turbulence does not arise as a successive superposition of oscillators, but irregular chaotic behavior can already appear after few bifurcations. A second direction is the explanation of nonlinear coupling between two incommensurable shedding frequencies (Luchtenburg et al., 2009), also referred to as frequency crosstalk in the following. This amplitude coupling over the mean flow has been termed quasi-laminar in Reynolds & Hussain (1972) pioneering theoretical foundation of the triple decomposition. The advancements also include subcritical bifurcations (Watson, 1960). More specifically, the case of a codimension two bifurcation, involving both a pitchfork and a Hopf bifurcation, was addressed in Meliga et al. (2009), who derived the amplitude equation based on the weakly nonlinear analysis of the wake of a disk. Fabre et al. (2008) derived the same equation solely based on symmetry arguments for the wake of axisymmetric bodies. A resolvent analysis follows mean-field considerations in decomposing the flow in a time-resolved linear dynamics and a feedback-term with the quadratic nonlinearity (Gomez et al., 2016; Rigas et al., 2017).

Our study develops a generalized mean-field Galerkin model for the first two bifurcations of the fluidic pinball with increasing Reynolds number. The primary supercritical bifurcation leads to the periodic vortex shedding which is statistically symmetric. At higher Reynolds numbers, the resulting limit cycle undergoes a pitchfork bifurcation into a stable, asymmetric, mirror-symmetric pair of periodic solutions. This local bifurcation has a transverse effect resulting from the decoupling of these two bifurcations (see appendix E), which simultaneously leads to an identical local pitchfork bifurcation of the steady solution, into an unstable, asymmetric, mirror-symmetric pair of steady solutions. The underlying dynamics is modeled with a small number of assumptions. The key simplification results from exploiting the symmetry of the mean flow and the antisymmetry of the fluctuation. The generalized mean-field Galerkin methodology can be expected to be useful for describing other transition scenarios.

The manuscript is organized as follows. In § 2, the numerical plant is introduced and the Reynolds-number dependent flow behavior described. This phenomenology drives the mean-field modeling of the first two bifurcations in § 3. The resulting models for the Hopf and subsequent pitchfork bifurcation are present in § 4 and § 5, respectively. § 6 summarizes the results and outline future directions of research.

2 Flow configuration

In this section we describe the numerical toolkit and the flow features as the Reynolds number is increased. The direct Navier-Stokes solver with MATLAB interfaces, used for the simulation, is described in § 2.1. The fluidic pinball configuration and the flow features and route to chaos are described in § 2.2 and § 2.3, respectively.

2.1 Direct Navier-Stokes solver

The unsteady Navier-Stokes solver is based on fully implicit time integration and Finite-Element Method discretization (Noack & Morzyński, 2017; Noack et al., 2003, 2016). The time integration is third-order accurate while FEM discretization employs a second–order Taylor-Hood finite elements (Taylor & Hood, 1973). The solution is obtained iteratively, with the Newton-Raphson type approach. The tangent matrix is updated on 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 1(a)). To test the grid dependency of the solution we increased the number of triangles by nearly a factor 4 (26849 elements and 54195 nodes, see figure 1(b)). 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 center 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 a 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).

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 frequencies interactions over the past few decades (Price & Paidoussis, 1984; Sayers, 1987; Lam & Cheung, 1988; Tatsuno et al., 1998; Bansal & Yarusevych, 2017). 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, 2017) as this problem is a challenging task for control methods comprising several frequency crosstalk mechanisms (Noack & Morzyński, 2017). The centers 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:

[TABLE]

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 $\bm{u}_{s}(\bm{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

The steady solution $\bm{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, 1999; Barkley, 2006). Beyond $Re_{1}$, the steady solution becomes unstable with respect to vortices periodically and alternately shed at the top and bottom of the two right-most cylinders, following a Hopf bifurcation (instability of the fixed point via a pair of complex-conjugated eigenvalues, see e.g. Strogatz et al. (1994)). 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 $\bm{u}_{s}$, the mean flow $\overline{\bm{u}}$ and the instantaneous flow field $\bm{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 or the limit cycle, via a real eigenvalue, see e.g. Strogatz et al. (1994). 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 $\bm{u}_{s}^{\pm}$ and the two associated mean flows $\overline{\bm{u}}^{\pm}$, shown in figure 3 for $Re=80$, where the base-bleeding jet appears deflected to either upward or downward with respect to the symmetry plane. Note, however, that a symmetry-preserving mean flow ($\overline{\bm{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{\bm{u}}^{\pm}$, are mirror-conjugated and break the symmetry, while the last one, $\overline{\bm{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}$: $\bm{u}_{s}^{\pm}$ are mirror-conjugated and break the symmetry, while the last one, $\bm{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 $\bm{u}_{s}$, the flow regime arrives after a long transient on a limit cycle whose mean flow $\overline{\bm{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 (centered on either $\overline{\bm{u}}^{\pm}$). When the initial condition already breaks the symmetry of the flow configuration, the unstable “symmetry-centered” 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}$ illustrate a transverse action on the original state space resulting from the new active symmetric 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, 1994a, b) and from 3D Navier-Stokes simulations (Zhang et al., 1994). 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.

As a result, when the symmetry of vortex shedding is broken, the mean value $\overline{C}_{L}$ (solid line) of the pressure lift coefficient $C_{L}=2F_{L}/\rho U^{2}$, where $F_{L}$ is the total lift force from pressure, 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 about 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. 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 5(b)) of the main oscillator and thickens the limit cycle associated with the main oscillator (figure 5(c)). 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\geq 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{\bm{u}}$ is symmetric, as shown in figure 4 for $Re>115$. The time series exhibits neither periodic nor quasi-periodic features anymore (see figure 6(b)) and the phase portrait exhibits a much more complex dynamics (see figure 6(c)). 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 et al., 1978).

3 Low-dimensional modeling

We derive a mean-field Galerkin model for the primary and secondary bifurcations of the fluidic pinball. First (§ 3.1), the Galerkin method is recapitulated as a very general approach to reduced-order models. In § 3.2, the constitutive equations of the mean-field model are derived from a minimal set of assumptions. Then, mean-field Galerkin models are derived for the Hopf bifurcation (§ 3.3), the pitchfork bifurcation (§ 3.4) and the succession of both bifurcations (§ 3.5).

3.1 Galerkin method

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

[TABLE]

where $\nu=1/Re$. The velocity field satisfies the no-slip condition $\bm{u}=0$ on the cylinders, the free-stream condition $\bm{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 $\bm{u}_{s}$ satisfies the steady Navier-Stokes equations:

[TABLE]

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

[TABLE]

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

[TABLE]

Orthonormality implies:

[TABLE]

The projection of (5) on (2) leads to the linear-quadratic Galerkin system (Fletcher, 1984),

[TABLE]

Following Rempfer & Fasel (1994), $a_{0}=1$ is introduced. The coefficients $l_{ij}^{\nu}=\left(\bm{u}_{i},\triangle\bm{u}_{j}\right)_{\Omega}$ and $q_{ijk}^{c}=\left(\bm{u}_{i},\nabla\cdot\bm{u}_{j}\otimes\bm{u}_{k}\right)_{\Omega}$ 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 $\bm{u}_{0}=\bm{u}_{s}$. This implies that $\bm{a}=0$ is a fixed point of (7) and the constant term $\nu l_{i0}^{\nu}+q_{i00}^{\nu}=0$ vanishes as the projection of $\eqref{Eqn:NSE:Steady}$ onto the $i$th mode $\bm{u}_{i}$. In this case, (7) can be re-written as a linear-quadratic system of ordinary differential equations :

[TABLE]

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

3.2 Mean-field modeling

Mean-field modeling allows a dramatic simplification of a general Galerkin system (8) 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 $\bm{u}^{D}$ and fluctuation $\bm{u}^{\prime}$ with first-order (relaxational) or second-order (oscillatory) dynamics:

[TABLE]

Here, the mean-field deformation $\bm{u}_{\Delta}$ 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,

[TABLE]

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

[TABLE]

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

[TABLE]

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\pi]$. The somewhat loaded term “distorted mean flow” is directly adopted from the original publications of mean-field theory (Stuart, 1958). 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 $\bm{u}^{s}=\left(u^{s},v^{s}\right)\in\mathcal{U}^{s}$ with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}(\Omega)$.

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

[TABLE]

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

[TABLE]

Analogously, the symmetric part describes the distorted mean flow dynamics:

[TABLE]

Note that this symmetric component of the Navier-Stokes equations has not yet been averaged and the sum of Eqs. (3), (19) and (18) leads to the Navier-Stokes equations (2). To this point, all equations are strict identities for the symmetric and antisymmetric part of the Navier-Stokes dynamics.

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

[TABLE]

The mean-field deformation $\bm{u}_{\Delta}$ characterized by the scale $\delta$ is seen to respond linearly to the Reynolds stress force $-\nabla\cdot\left\langle\bm{u}^{\prime}\otimes\bm{u}^{\prime}\right\rangle$ scaling with $\epsilon^{2}$. Hence, $\delta\sim\epsilon^{2}$.

Summarizing, Eqs. (18) and (20) 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 $\sigma=\alpha(Re-Re_{c})$ and can be assumed to be small. In this case, $\partial_{t}\bm{u}_{\Delta}\in O(\sigma\delta)$, i.e. the time derivative of (20) can be neglected with respect to the other terms $\in O(\delta)=O(\varepsilon^{2})$. This leads to the steady linearized Reynolds equation:

[TABLE]

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 $\bm{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 $\bm{u}_{s}\in\mathcal{U}^{s}$ is observed. Periodic vortex shedding sets in with the occurrence of an unstable oscillatory antisymmetric eigenmode at $Re\geq Re_{1}$. The Reynolds-number dependent initial growth rate and frequency are denoted by $\sigma_{1}$ and $\omega_{1}$, respectively. The real and imaginary parts of this eigenmode are $\bm{u}_{1}$ and $\bm{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 (20). The mean flow deformation is described by the symmetric shift mode $\bm{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:

[TABLE]

Moreover, polar coordinates are introduced $a_{1}(t)=r(t)\cos\theta(t)$, $a_{2}(t)=r(t)\sin\theta(t)$, $d\theta/dt=\omega(t)$ where $r$ and $\omega$ are assumed to be slowly varying functions of time.

Substituting (22) in (18), projecting on $\bm{u}_{i}$, $i=1,2$ and applying the Krylov-Bogoliubov (Jordan & Smith, 1999) averaging method yields:

[TABLE]

Here, $\sigma_{1},\omega_{1},\beta>0$ for a supercritical Hopf bifurcation. We refer to Noack et al. (2003) for details.

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 percent error.

Substituting (22) in (21) replaces (25) by the mean-field manifold:

[TABLE]

with derivable proportionality constant $\kappa$.

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

[TABLE]

where $\sigma_{3}<0$ and $\beta_{3}>0$ are necessary for a globally stable limit cycle. Note that (25) can be rewritten as:

[TABLE]

Now, the slaving process which leads to the mean-field manifold of (24) can be appreciated from the mean-field Galerkin system. If $|\sigma_{3}|\gg\sigma_{1}$, the timescale of slaving $a_{3}$ to the fluctuation level $a_{1}^{2}+a_{2}^{2}$ is much smaller than the timescale of the transient and $da_{3}/dt$ can be set to zero.

Eqs. (23) and (24) yield the famous Landau equations with the cubic damping term:

[TABLE]

The Landau oscillator leads to a stable limit cycle with $r^{\circ}=\sqrt{\sigma_{1}/\beta\kappa}$, frequency $\omega^{\circ}=\omega_{1}+\sigma_{1}\gamma/\beta$ and shift-mode amplitude $a_{3}^{\circ}=\sigma_{1}/\beta$. The three nonlinearity parameters $\beta$, $\gamma$ and $\kappa$ can be uniquely derived from the limit cycle parameters $r^{\circ}$, $\omega^{\circ}$, and $a_{3}^{\circ}$. The growth rate $\sigma_{3}$ needs to be chosen sufficiently large, e.g. $\sigma_{3}=-10\sigma_{1}$ to ensure slaving on the manifold.

Eqs. (23), (25) are the mean-field Galerkin system, while Eqs. (23), (24) characterize the original mean-field model, i.e. the slaved Galerkin system. Near the Hopf bifurcation, when $\sigma_{1}(Re)$ crosses the zero line at $Re_{1}$, the growth rate is approximated by $\sigma_{1}=\alpha\left(Re-Re_{1}\right)$ implying the square-root law $r^{\circ}=\sqrt{\alpha/\beta\kappa}\sqrt{Re-Re_{1}}$.

The Landau equation has been proposed by Landau (see, e.g. Landau & Lifshitz, 1987), derived from the Navier-Stokes equation by Stuart (1958), generalized for Galerkin systems by Noack et al. (2003), and validated in numerous simulations and experiments for cylinder wakes (Schumm et al., 1994) and other soft onsets of oscillatory flows. 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, likes closeness to the Hopf bifurcation or the need for frequency filtered Navier-Stokes equations.

3.4 Supercritical pitchfork bifurcation

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

[TABLE]

Substituting (28) in Eqs. (18) and (20), and exploiting the symmetry of the modes, yields:

[TABLE]

Note that a linear $a_{5}$ term and quadratic $a_{4}a_{4}$ and $a_{5}a_{5}$ terms in (29a) are ruled out by symmetry. Similarly, a linear $a_{4}$ term or a mixed quadratic term $a_{4}a_{5}$ in (29b) are prohibited by symmetry. The quadratic $a_{5}a_{5}$ term is not consistent with the linearized Reynolds equation (20). The pitchfork bifurcation can be considered as a Hopf bifurcation with $\omega=0$ and a single mode. Replacing $a_{1}$ by $a_{4}$, $a_{3}$ by $a_{5}$ and setting $a_{2}=0$ yields (29a) from (23) and (29b) from (25)—modulo names of the coefficients.

Eqs. (29a), (29b) are the mean-field Galerkin system. Substituting (28) in (21) yields the manifold:

[TABLE]

with $\kappa_{5}=-\beta_{5}/\sigma_{5}$. The asymmetric steady solutions read $a_{4}^{\pm}=\pm\sqrt{\sigma_{4}/\beta_{4}\kappa_{5}}$, $a_{5}=\sigma_{4}/\beta_{4}$. The two nonlinearity parameters $\kappa_{4}$, $\beta_{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 $\sigma_{5}=-10\sigma_{4}$ to ensure slaving on the manifold.

From (29a) and (30), the famous unstable dynamics with cubic damping term is obtained:

[TABLE]

where $\kappa_{5}\beta_{4}>0$ for a supercritical bifurcation. Eqs. (29a), (29b) are the mean-field Galerkin system.

Near the secondary pitchfork bifurcation, $\sigma_{4}=\alpha_{2}\left(Re-Re_{2}\right)$ 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 $\sigma_{5}=-10\sigma_{4}$ will ensure the slaving of (30). 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 modeling 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}$ is 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:

[TABLE]

describes the flow where $a_{1},a_{2},a_{3}\in O(1)$ and $a_{4}\in O(\varepsilon)$ and $a_{5}\in O(\delta)$, $\varepsilon$, $\delta$ being smallness parameters associated with the pitchfork bifurcation. We project Eqs. (18) and (20) on (31). The $O(1)$ terms encapsulate the original Hopf model while the low-pass filtered $O(\varepsilon,\delta)$ terms yield the original pitchfork system. This yields the following generalized mean-field system:

[TABLE]

The linear instability parameters $\sigma_{1}$, $\omega_{1}$, $\sigma_{4}$ are obtained from the corresponding global stability analysis. Slaving is ensured with $\sigma_{3}=-10\sigma_{1}$ and $\sigma_{5}=-10\sigma_{4}$. The nonlinearity parameters $\beta$, $\gamma$, $\beta_{3}$, $\beta_{4}$ and $\beta_{5}$ are determined from the limit cycle parameters $r^{\circ}$, $\omega^{\circ}$ and $a_{3}^{\circ}$ and pitchfork parameters $a_{4}^{\pm}$ and $a_{5}^{\pm}$ in the asymptotic regime.

As the amplitude of the pitchfork bifurcation grows, the smallness argument does not hold and we get cross-terms, like $\sigma=\sigma_{1}-\beta a_{3}-\beta_{15}a_{5}$. We shall not pause to elaborate on the possible generalizations now, but will return to the topic in the result section.

3.6 Sparse Galerkin model from mean-field considerations

The mean-field Galerkin system (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 $\sigma$ may also depend on the pitchfork-related shift mode amplitude $a_{5}$. The most general Galerkin system (8) 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 $\chi_{i}=0$ for symmetric mean flow modes $\bm{u}_{i}$, $i=3,5$ and $\chi_{i}=1$ for the antisymmetric fluctuation modes $\bm{u}_{i}$, $i=1,2,4$. The linear coefficients

[TABLE]

can be shown to vanish if $\mod(\chi_{i}+\chi_{j},2)=1$. In other words, the coefficients $l_{ij}$ vanish if the modes $\bm{u}_{i}$ and $\bm{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 $\chi_{i}\not=\mod(\chi_{j}+\chi_{k},2)$. In other words, the $q_{ijk}$ vanishes if the symmetry of $\bm{u}_{i}$ does not coincide with the symmetry of quadratic term $\bm{u}_{j}\otimes\bm{u}_{k}$. The quadratic term is symmetric if the modes $\bm{u}_{j}$ and $\bm{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 $\bm{u}_{i}$, $i=1,2$ have oscillatory behaviour with angular frequency $\omega$, 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 $da_{4}/dt$. 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 $da_{4}/dt$ but $q_{113}a_{1}a_{3}$ can contribute to the oscillatory behaviour of $da_{1}/dt$. Similarly, $a_{1}a_{2}$ generates a second harmonic, $a_{1}a_{2}=r\cos\omega t\times r\sin\omega t=(\nicefrac{{1}}{{2}})\>r^{2}\>\sin 2\omega t$ does not have a steady contribution, so $q_{312}$ can be set to zero.

From symmetry and Krylov-Bogoliubov considerations, only 9 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 8 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. $\sigma=\sigma_{1}-\beta a_{3}-\beta_{15}a_{5}$ and $\omega=\omega_{1}+\gamma a_{3}+\gamma_{15}a_{5}$, introducing $\beta_{15}$ and $\gamma_{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 Primary flow regime

The primary flow regime covers the range of Reynolds numbers $Re_{1}<Re<Re_{2}$. We consider the flow and reduction of the dynamics at $Re=30$, as a representative case. In § 4.1, a linear stability analysis is done on the steady solution and the three degrees of freedom of the flow dynamics are identified. In § 4.2, we propose a least-order model of the flow dynamics at $Re=30$ and compare its performance with respect to the full flow dynamics.

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 $\bm{u}^{D}=\bm{u}_{s}+\bm{u}_{\Delta}$ changes, as expected by Eq. (3) and (19). The shift mode $\bm{u}_{3}$, involved in $\bm{u}_{\Delta}$ at $Re=30$, is shown in figure 8(c). In the permanent (time-periodic) flow regime, the distorted mean flow $\bm{u}^{D}$ eventually matches the asymptotic mean flow field $\overline{\bm{u}}$, and vortex shedding is well established with frequency $8.7\times 10^{-2}$.

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

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

[TABLE]

The von Kármán modes $\bm{u}_{1,2}$ could be chosen as the real and imaginary part 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 $\bm{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 to 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 $\bm{u}_{1,2}$ are preferred to the two leading eigenmodes, because we demand an accurate representation of the asymptotic periodic dynamics. Following Eq. (23)-(24), the dynamical system resulting from the Galerkin projection of ansatz (33) on the Navier-Stokes equations, after Krylov-Bogoliubov simplifications, reads

[TABLE]

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

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 $\bm{u}_{1,2}$, the dynamics on the asymptotic (permanent) limit cycle is well described in amplitude $r$ and angular frequency $\omega$ by the ROM. Moreover, the ROM also captures the transient dynamics on the parabolic manifold $a_{3}\equiv\kappa(a_{1}^{2}+a_{2}^{2})$. Henceforth, although all coefficients but one are fixed, the Galerkin system (33) is able to reproduce the most salient dynamical features of the flow in both the transient and the permanent regimes.

5 Secondary flow regime

The secondary flow regime ranges over $Re_{2}<Re<Re_{3}$. For illustration, we focus on the flow at $Re=80$, i.e. at a finite distance from the secondary bifurcation. In § 5.1 a linear stability analysis of the resulting three steady solutions is performed. A least-order model is proposed and discussed in § 5.2.

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 $\bm{u}_{s}$, unstable to the periodic vortex shedding beyond $Re_{1}$, and two mirror-conjugated asymmetric steady solutions $\bm{u}_{s}^{\pm}$, also unstable to vortex shedding but only existing beyond $Re_{2}$, see figure 3.

The linear stability analysis on $\bm{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 eigenmode has become unstable and can only be visited transiently in time.

The linear stability analysis on $\bm{u}_{s}^{-}$ (resp. $\bm{u}_{s}^{+}$) reveals two pairs of complex-conjugated eigenmodes with positive growth rate, centered 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.

5.2 Reduced-order model in the secondary flow regime

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

[TABLE]

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

[TABLE]

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

[TABLE]

and further orthonormalized to $u_{1}$, $u_{2}$, $u_{3}$ by a Gram-Schmidt procedure. The resulting mode $\bm{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 $\bm{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 $\bm{u}_{4}$.

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

[TABLE]

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

Close to the pitchfork bifurcation, the resulting dynamical system is described by Eqs. (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 et al., 2018). We avoid this complication by model identification from simulation data. The coefficients of the mean-field system (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 timescale 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 value 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 (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:

1. Step 1:

Keep the five-dimensional linear-quadratic form of the dynamical system from the Galerkin projection with 25 linear and 75 quadratic terms. 2. 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. 3. 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 $\sigma_{1}$ and frequency $\omega_{1}$ characterize the initial growth and angular frequency of $a_{1}$, $a_{2}$ and $\sigma_{4}$ represents the linear growth rate for $a_{4}$. 4. Step 4:

The slaving of the Reynolds-stress-induced modes $\bm{u}_{3}$ and $\bm{u}_{5}$ to the fluctuation level is imposed by setting the damping rate 10 times larger than the growth rate of the corresponding fluctuation: $\sigma_{3}=-10\sigma_{1}$, $\sigma_{5}=-10\sigma_{4}$. This strong damping rate quickly forces the trajectory onto the mean-field manifold. 5. 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}$ are governed by $\sigma_{1}$, $\omega_{1}$, $\beta$, $\beta_{15}$, $\gamma$, $\gamma_{15}$, and $\beta_{3}$. 6. Step 6:

Impose the asymptotic dynamics of the unstable symmetric limit cycle by fixing $\beta$, $\gamma$, and $\beta_{3}$. 7. Step 7:

Apply the SINDy algorithm (Brunton et al., 2016) 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 (2015). 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. 8. 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.

The resulting ROM reads:

[TABLE]

where $r^{2}=a_{1}^{2}+a_{2}^{2}$. The coefficients are summarized in table 4. The dynamics of the system (41) (dashed red line) is compared to the full flow dynamics of the fluidic pinball (solid blue line) in figure 13. Now the initial stage of the dynamics is much better reproduced, as well as the asymptotic oscillations of $a_{3}$, $a_{4}$ and $a_{5}$. Interestingly, the faster growth in $a_{1}$ to $a_{3}$, on the time-range around 600, could not be completely corrected. It is worthwhile noticing that this range of time also corresponds to oscillations in $a_{4}$, which could not be reproduced by any cross-terms compatible with the symmetries of the system. Noticing that $a_{4}\neq 0$ on this range of time may also question our choice for $\sigma_{1}=5.22\times 10^{-2}$, which is the linear growth rate of the leading eigenmode around the symmetric steady solution $\bm{u}_{s}$. Indeed, although the initial condition is close to this point, the large amplitude oscillations of $a_{4}$ on the time range around 600 mean that the trajectory transiently escapes the symmetric subspace not only along $a_{3}$, but also along $a_{4}$ and $a_{5}$, before coming back close to the $a_{3}$ axis, on the time range from 700 to about 800. Therefore, it would be reasonable here to keep all the coefficients of the model unconstrained, but the number of free parameters is now too large for the identification process to be computationally tractable — for instance, it provides positive $\sigma_{3}$ when it must necessarily be strongly negative.

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, 2015). Second, the Galerkin method assumes fixed expansion modes. Yet, least-order models often have deformable modes changing with the fluctuation level (Tadmor et al., 2011). 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 $\sigma_{3}$ and $\sigma_{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, timescales, transient and asymptotic behaviour of the full dynamics.

6 Conclusions and outlooks

Reduced-order models (ROM) serve a number of purposes. For instance, ROMs facilitate a deeper understanding of the physical mechanisms at play in a flow configuration, by extracting the low-dimensional manifold on which evolve the dynamics (Manneville, 2010). In that respect, the linear stability analysis of the steady solution shows the nonlinear amplitude saturation mechanism through the distorted mean flow. The difference between steady solution and the mean flow is caused by the Reynolds stress, captured by the shift mode $\bm{u}_{\Delta}$, and affects the stability properties (Noack et al., 2003; Barkley, 2006; Sipp & Lebedev, 2007; Turton et al., 2015). Deeper investigations would certainly deserve to be carried out on the relation between the stability analysis of fixed points, Floquet analysis of limit cycles and Lyapunov exponents of chaotic flow regimes. In addition, because Hopf and pitchfork bifurcations are generic bifurcations in fluid flows, the nonlinear dynamics identified in this study are expected to be extended and generalizable to other flows exhibiting similar bifurcations. Last but not least, a ROM provides fast estimators for predicting the forward evolution of the system. Such estimators could be used for control purposes (Brunton & Noack, 2015; Rowley & Dawson, 2017). All these considerations motivated the present study, whose main results are summarized in § 6.1. Outlooks of this work are listed in § 6.2.

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, 1999; Gumowski et al., 2008; Szaltys et al., 2012; Grandemange et al., 2014) or bluff body wake flows (Grandemange et al., 2012, 2013; Cadot et al., 2015; Bonnavion & Cadot, 2018; Rigas et al., 2014). 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 modeling of Hopf bifurcations was already addressed in Noack et al. (2003) for the cylinder wake flow, while Meliga et al. (2009) 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. (2008) 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 bifurcation could be well-captured by a five-dimensional model whose degrees of freedom couple through quadratic non-linearities, 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.

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. (2003). 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 a 8-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 8-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 $\sigma=\sigma_{1}(Re)-\beta(Re)a_{3}$. The transients can be expected to be much more accurately resolved by the mean flow dependent modes, e.g. $\bm{u}_{i}(Re,a_{3},a_{5},\vec{x})$, $i=1,2$ for the resolution of vortex shedding (Loiseau et al., 2018). More generally, the flow lives on a low-dimensional manifold which includes mode deformations (Noack, 2016). Locally linear embedding (LLE) is a powerful technique for identifying the dimension and a parameterization in an automatic manner (Roweis & Saul, 2000). 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 (1994) 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. (2008); Luchtenburg et al. (2009). Such multi-frequency models may be extended to resolve broadband frequency dynamics taking, for instance, the most dominant DMD modes (Rowley et al., 2009; Schmid, 2010). Strengths and weaknesses of techniques currently used for model reduction are discussed in Taira et al. (2017). While mean-field consideration expressly ignores non-trivial triadic interactions, their quantitative effect on the frequency crosstalk 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 (2018) on turbulent flow modeling.

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., 2016). The choice of the state spaces might be facilitated by manifold learning from many solution snapshots (Gorban & Karlin, 2005). The authors actively pursue all the mentioned directions.

Acknowledgements

This work is supported by a public grant overseen by the French National Research Agency (ANR) as part of the “Investissement d’Avenir” program, through the “iCODE Institute project” funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02, by the ANR grants ‘ACTIV_ROAD’ and ‘FlowCon’ (ANR-17-ASTR-0022), and by Polish Ministry of Science and Higher Education (MNiSW) under the Grant No.: 05/54/DSPB/6492.

We appreciate valuable stimulating discussions with Steven Brunton, Alessandro Bucci, Nathan Kutz, Onofrio Semeraro, Yohann Duguet, Laurette Tuckerman and the French-German-Canadian-American pinball team: François Lusseyran, Guy Cornejo-Maceda, Jean-Christophe Loiseau, Robert Martinuzzi, Cedric Raibaudo, Richard Semaan, and Arthur Ehlert.

Appendix A Asymmetric steady solutions

For a Reynolds number larger than the critical value of the pitchfork bifurcation $Re_{2}$, we can obtain two additional asymmetric steady solutions, one associated with the base-bleeding jet deflected upward, the other downward. These two asymmetric steady solutions are obtained by the steady Navier-Stokes solver initialized with a flow field (snapshot) with the same state of the base-bleeding jet. From the two asymmetric vortex shedding, their corresponding asymmetric steady solutions can be obtained by the following steps:

1. Step 1:

Run the unsteady Navier-Stokes solver with a time scale larger than the vortex shedding period, initialized from a snapshot of the asymmetric vortex shedding. The vortex shedding will quickly vanish for this artificially large time scale and approach the corresponding approximation of the steady state. 2. Step 2:

Run the steady Navier-Stokes solver, restarted from this vortex shedding vanished solution to further refine the steady state.

Initialized with the three steady solutions with different states of the base-bleeding jet at $Re=80$, the steady solutions for other Reynolds numbers can be obtained by the steady Navier-Stokes solver. At $Re=68$, all three steady solutions converge to a unique solution, which indicates that the critical value of the pitchfork bifurcation of the steady solution is between $68$ to $69$, as shown in figure 15 based on the pressure lift coefficient of the steady solutions.

Appendix B Linear stability analysis

The linear stability problem for a base flow $\left(\bm{U}(\bm{x}),P(\bm{x})\right)$ with small perturbations $\left(\bm{u^{\prime}}(\bm{x},t),p^{\prime}(\bm{x},t)\right)$ is governed by the linearized Navier-Stokes equations, which read:

[TABLE]

The assumption of small perturbation allows to linearize the equations, and we can separate the time and space dependence as:

[TABLE]

By introducing the linear operator $\bm{L}(\bm{U})$, all the terms except the time derivative term and the continuity equation can be cast as $\bm{L}(\bm{U})\bm{u^{\prime}}$, we can rewrite (42) as

[TABLE]

Introducing (43) into (44), the equations can be written as:

[TABLE]

We use subspace iteration to solve this eigenvalue problem. A detailed review can be found in Morzynski et al. (1999).

The global stability analysis of the steady solutions at different Reynolds numbers has been performed on a Krylov subspace of dimension 9-20. This converges after 50-100 iterations. The linear stability analysis on the symmetric steady solution $\bm{u}_{s}$ reveals a pair of conjugated eigenvalues with positive real part first appearing as the Reynolds number is changing from 18 to 19, see figure 16(a). A real eigenvalue becomes positive between $Re=68$ and $69$, see figure 16(b). This confirms that a Hopf bifurcation occurs on the symmetric steady solution at $Re_{1}\approx 18$, and a pitchfork bifurcation at $Re_{2}\approx 68$.

Appendix C Floquet stability analysis

Similar to the linear stability framework, the Floquet stability problem works with a T-periodic base flow $\left(\bm{U}(\bm{x},t),P(\bm{x},t)\right)$. The linear operator now reads:

[TABLE]

The linear operator $\bm{L}(\bm{U}(\bm{x},t))$ is T-periodic because of the base flow $\bm{U}(\bm{x},t)$. The solutions to (46) are seeked as:

[TABLE]

with the T-periodic Floquet modes $\bm{\hat{u}}(\bm{x},t)$ and the corresponding Floquet exponents $\sigma+i\omega$. We define the Floquet operator as the time-integrated $\bm{L}(\bm{U}(\bm{x},t))$ with the pre-stored periodic solutions over one period (Barkley & Henderson, 1996; Schatz et al., 1995), which reads:

[TABLE]

The method used to solve the eigenproblem is the same as the linear stability analysis. The Floquet multipliers of $A_{F}$ can be written as $\lambda_{F}=e^{(\sigma+i\omega)T}$. In our case, we consider the multipliers $\lambda=e^{(\sigma+i\omega)}$.

We performed the iteration on a Krylov subspace of dimension 9, using a symmetry-constrained T-periodic base flow. Below the critical Reynolds number, the base-bleeding jet is approximately steady and symmetric. This symmetry is also enforced at higher Reynolds numbers to compute the unstable periodic solution. The constraint is imposed on the central line as:

[TABLE]

As (49) restricts the vertical velocity of the nodes, the symmetry-constrained periodic solution is very close to the symmetry-preserving periodic solution. Normally, 10-20 iterations are enough to get a converged leading eigenvalue, initialized with a random vector or a Ritz vector computed at nearby Reynolds number with the periodic solution of the same symmetry. We do not attempt to calculate a complete, converged spectrum of all the eigenvalues, as the leading eigenvalue is associated with the instability of interest. The multipliers from the Floquet analysis around the critical value of the pitchfork bifurcation $Re_{2}$ are shown in figure 17. When increasing $Re$, the leading real eigenvalue crosses the unit cycle at $(+1,0)$ as $Re$ changes from 69 to 70. The critical value of the pitchfork bifurcation is, therefore, $Re_{2}\approx 69$, identical to the critical value for the steady solution at the precision of the numerics. Both of them have the same eigenmode. At $Re=80$, the Floquet modes of both the unstable symmetric periodic solution and the stable asymmetric periodic solutions are shown in figure 18.

Overall, combined with the result of the linear stability analysis of the steady solutions, the bifurcation scenario at low Reynolds numbers can be shown in figure 19. The linear stability analysis on the steady solution and the periodic solution show a highly consistent result: the same kind of bifurcation with nearly the same critical Reynolds number, and the same eigenmodes. Besides, at $Re=80$, the growth rates of the real eigenmode are very close: 0.0272 from the symmetric steady solution, 0.0247 from the symmetric periodic solution. This similarity is understood as the result of a transverse effect of the symmetric subspace.

Appendix D Transient dynamics from different steady solutions

In this section, we show in figure 20 some typical transient dynamics starting with the unstable symmetric/asymmetric steady solutions at different Reynolds numbers, based on the pressure lift coefficient $C_{L}(t)$ from the resulting force on the three cylinders. Combining $C_{L}(t)$ with the pressure drag coefficient $C_{D}(t)$ and the time-delayed lift coefficient $C_{L}(t-\tau)$ in which $\tau$ is a quarter period, provides the phase portraits of figure 21. Three comparative numerical simulations are shown, starting with the symmetric and the two mirror-conjugated asymmetric steady solutions at the same $Re$ respectively. The mirror-conjugated initial conditions provide mirror-conjugated transient dynamics.

Figure 20(a) shows the transient dynamics at the critical Reynolds number $Re_{2}$, initialized with the three steady solutions at $Re=75$. The lift coefficient starts oscillating quickly, and eventually reaches a unique oscillating state with zero mean value. This is consistent with the Floquet analysis at $Re=68$, where only one stable symmetry-centered limit cycle exists and any other state will eventually converge to this stable state.

Figure 20(b) shows three different scenarios at $Re=75$ depending on the initial condition: from the symmetric steady solution $\bm{u}_{s}$ to the symmetry-centered limit cycle, from the symmetric centered limit cycle to the asymmetry-centered limit cycles, and from the asymmetric steady solutions $\bm{u}_{s}^{\pm}$ to the asymmetry-centered limit cycles. Starting with the symmetric steady solution, it will first reach the unstable symmetry-centered limit cycle, before asymptotically approaching one of two stable asymmetry-centered limit cycles. However, starting with the asymmetric steady solutions, it will directly reach the corresponding stable asymmetry-centered limit cycle. If the initial perturbation introduced to the symmetric steady solution has a certain bias of symmetry, a transition from the symmetric steady solution $\bm{u}_{s}$ to one of the two asymmetry-centered limit cycles will occur.

As we keep increasing the Reynolds number up to $Re=100$, there still exist six states, but the transient scenario from the symmetric steady solution $\bm{u}_{s}$ is different. It will first reach one of two unstable asymmetric steady solutions $\bm{u}_{s}^{\pm}$, before asymptotically approaching the corresponding asymmetry-centered limit cycle, as shown in figure 20(c).

All these transient dynamics mentioned above are highlighted with the red arrows in figure 19.

The phase portraits starting with different steady solutions, as shown in figure 21, reveal the above-mentioned transient dynamics which is affected by the initial condition and the Reynolds number. At the same time, it also reflects the global effect of the pitchfork bifurcation at $Re_{2}$, splitting the state space (see figure 21(a)) into a symmetric sub-space and two mirror-conjugated asymmetric sub-spaces (see figure 21(b),(c)).

Appendix E On the simultaneous instability of the fixed point and the limit cycle

In this section, we exemplify the transverse effect of the pitchfork bifurcation on a three-dimensional dynamical system equivalent to the system of Eq. (32).

Dynamical system

The dynamical system reads:

[TABLE]

with $\mu_{1}=1,\mu_{2}=2,\omega_{0}=1$. This system undergoes a supercritical Hopf bifurcation in the $(x,y)$-plane at $\mu=\mu_{1}$ and a supercritical pitchfork bifurcation along the $z$-axis at $\mu=\mu_{2}$. For $\mu>\mu_{1}$, the stable fixed point at ${(0,0,0)}$ becomes unstable, and the limit cycle around $(0,0,0)$ with radius $r=\sqrt{\mu-\mu_{1}}$ and angular frequency $\omega=\omega_{0}+\mu-\mu_{1}$ is stable in the (x,y)-plane. Increasing $\mu$ until $\mu>\mu_{2}$, the fixed point undergoes a secondary instability, as well as the limit cycle. Three unstable fixed points ${(0,0,0)}$, ${(0,0,\pm\sqrt{\mu-\mu_{2}})}$, and three limit cycles around these fixed points with radius $r=\sqrt{\mu-\mu_{1}}$ and angular frequency $\omega=\omega_{0}+\mu-\mu_{1}$ in the (x,y)-plane are found, as shown in figure 22.

Linear stability analysis

The ODE (50) can be writen as:

[TABLE]

and we note $\bm{q}_{s}$ is the steady state,that is, $\bm{F}(\bm{q}_{s})=0$. Consider a small perturbation $\bm{q^{\prime}}$ around the steady state $\bm{q}_{s}$ by

[TABLE]

We derived the linearized evolution equation:

[TABLE]

where $\bm{DF(\bm{q}_{s})}$ is the Jacobian matrix of the considered steady state $\bm{q}_{s}$.

The stability of this steady state is determined by the eigenvalues $\sigma+i\omega$ of the Jacobian matrix. The eigenspectrum of the fixed point $\bm{q}_{s}^{0}=(0,0,0)$ is shown in figure 23. The growth rate and angular frequency of the pair of conjugated eigenvalues are $\mu-\mu_{1}$ and $\omega_{0}$ respectively. The growth rate of the real eigenvalue is $\mu-\mu_{2}$.

The eigenspectrum of the steady solution $\bm{q}_{s}^{\pm}=(0,0,\sqrt{\mu-\mu_{2}})$ for $\mu>\mu_{2}$ is shown in figure 26(right). The growth rate is $\mu-\mu_{2}$ at the steady solution $\bm{q}_{s}^{0}$ and $-2(\mu-\mu_{2})$ at the steady solution $\bm{q}_{s}^{\pm}$.

Floquet stability analysis

Now, we consider the periodic solution $\bm{q}_{p}(t)$ of the system of (50), which can be writen as:

[TABLE]

Consider a small perturbation $\bm{q^{\prime}}$ around the periodic solution by

[TABLE]

The first variational form reads:

[TABLE]

where $\bm{DF}(\bm{q}_{p}(t))$ is the Jacobian matrix of the considered periodic solution $\bm{q}_{p}(t)$, but now, the linear equation has periodic coefficients.

The monodromy matrix can be written as:

[TABLE]

or equivalently:

[TABLE]

The stability of this periodic solution is determined by the multipliers $\lambda=\exp((\sigma_{F}+i\omega_{F})/T)$, being the eigenvalues of the monodromy matrix. As shown in figure 25, the periodic solution $\bm{q}_{p}^{0}$: $(\sqrt{\mu-\mu_{1}}\cos(\theta),\sqrt{\mu-\mu_{1}}\sin(\theta),0)$ for $\mu\geq\mu_{2}$ becomes unstable. A real multiplier $\lambda$ crosses the unit cycle at $+1$. The Floquet exponent of this real multiplier is $\mu-\mu_{2}$, equal to the growth rate at the fixed point $\bm{q}_{s}^{0}$.

The multipliers of the other two periodic solutions $\bm{q}_{p}^{\pm}$: $(\sqrt{\mu-\mu_{1}}\cos(\theta),\sqrt{\mu-\mu_{1}}\sin(\theta),$ $\pm\sqrt{\mu-\mu_{2}})$ is shown in figure 26(right). The Floquet exponent of the leading multiplier is $-2(\mu-\mu_{2})$, equal to the growth rate at the fixed point $\bm{q}_{s}^{\pm}$.

As a consequence, if the additional degree of freedom $z$, introduced by the pitchfork bifurcation, do not couple, at the onset of the bifurcation, to the primary degrees of freedom $x,y$ associated with the Hopf bifurcation (see § 3.5), then it is easy to understand that both the symmetric steady solution $(x_{s},y_{s},z_{s})\equiv(0,0,0)$ and the statistically symmetric periodic solution $(x_{c}(t+T),y_{c}(t+T),0)=(x_{c}(t),y_{c}(t),0)$, $T$ being the period of the limit cycle, will both undergo an instability with respect to the symmetry breaking provoked by the pitchfork bifurcation ($z\neq 0$). This simultaneous instability of both the symmetric fixed point and the statistically symmetric limit cycle looks like an instability of the subspace ($x,y$) with respect to the transverse direction $z$, associated with the active degree of freedom introduced by the pitchfork bifurcation.