2.1. The Galerkin framework

The fluid flow satisfies the non-dimensionalized incompressible Navier–Stokes equations

(2.1)\begin{equation} \partial_t \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} \otimes \boldsymbol{u} = \nu \triangle \boldsymbol{u} - \boldsymbol{\nabla} p, \end{equation}

where $p$ and $\boldsymbol {u}$ are respectively the pressure and velocity flow fields, $\nu = 1/Re$, with the Reynolds number $Re$. Here, $\partial _t$, $\boldsymbol {\nabla }$, $\triangle$, $\otimes$ and $\boldsymbol {\cdot }$ respectively denote the partial derivative in time, the nabla and Laplace operators as well as the outer and inner tensor products. All the variables have been non-dimensionalized, with the cylinder diameter $D$, the oncoming velocity $U$, the time scale $D/U$ and the density $\rho$ of the fluid.

It is assumed that there exists at least one steady solution $(\boldsymbol {u}_s,p_s)$, satisfying the steady Navier–Stokes equations

(2.2)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_s \otimes \boldsymbol{u}_s = \nu \triangle \boldsymbol{u}_s - \boldsymbol{\nabla} p_s. \end{equation}

For the Galerkin framework, the space of the square-integrable vector fields $\mathcal {L}^{2} ( \varOmega )$ is introduced in the observation domain $\varOmega$. The associated inner product for two velocity fields $\boldsymbol {u} (\boldsymbol {x})$ and $\boldsymbol {v} (\boldsymbol {x})$ reads

(2.3)\begin{equation} \left( \boldsymbol{u} , \boldsymbol{v} \right)_{\varOmega} := \int\limits_{\varOmega} \,\textrm{d}\boldsymbol{x} \, \boldsymbol{u}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{v}(\boldsymbol{x}). \end{equation}

The velocity field is decomposed in a basic mode $\boldsymbol {u}_0$ and a fluctuating contribution. The basic mode may be the steady Navier–Stokes solution $\boldsymbol {u}_s$ or the time-averaged flow $\bar {\boldsymbol {u}}$. The fluctuation is represented by a Galerkin approximation of $N$ orthonormal space-dependent modes $\boldsymbol {u}_i ( \boldsymbol {x} )$, $i=1,\ldots ,N$, with time-dependent amplitudes $a_i(t)$

(2.4)\begin{equation} \boldsymbol{u} (\boldsymbol{x},t ) = \sum_{i=0}^{N} a_i(t) \boldsymbol{u}_i ( \boldsymbol{x} ), \end{equation}

where the basic mode $\boldsymbol {u}_0$ is associated with $a_0 \equiv 1$ following Rempfer & Fasel (Reference Rempfer and Fasel1994). The orthonormality condition reads $( \boldsymbol {u}_i , \boldsymbol {u}_j )_{\varOmega } = \delta _{ij}$, $i,j \in \{1,\ldots ,N\}$.

The Galerkin expansion (2.4) satisfies the incompressibility condition and the boundary conditions by construction. The evolution equation for the mode amplitudes $a_i$ is derived by a Galerkin projection of the Navier–Stokes equation (2.1) onto the modes $\boldsymbol {u}_i$

(2.5)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} a_i = \nu \sum_{j=0}^{N} l_{ij}^{\nu} a_j + \sum_{j,k=0}^{N} q_{ijk}^{c} a_j a_k + \sum_{j,k=0}^{N} q_{ijk}^{p} a_j a_k, \end{equation}

with the coefficients $l_{ij}^{\nu } = ( \boldsymbol {u}_i, \triangle \boldsymbol {u}_j )_{\varOmega }$, $q_{ijk}^{c} = ( \boldsymbol {u}_i, \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_j \otimes \boldsymbol {u}_k )_{\varOmega }$ and $q_{ijk}^{p} = ( \boldsymbol {u}_i, -\boldsymbol {\nabla } p_{jk} )_{\varOmega }$ for the viscous, convective and pressure terms in the Navier–Stokes equations (2.1), respectively. Details are provided by Noack, Papas & Monkewitz (Reference Noack, Papas and Monkewitz2005). Thus, a linear-quadratic Galerkin system (Fletcher Reference Fletcher1984) can be derived,

(2.6)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} a_i = \nu \sum_{j=0}^{N} l_{ij}^{\nu} a_j + \sum_{j,k=0}^{N} \left[ q_{ijk}^{c} + q_{ijk}^{p}\right] a_j a_k. \end{equation}

2.2. Drag and lift forces on a body

Let ${\varGamma }$ be the boundary of the body in the flow domain $\varOmega$ and $\boldsymbol {n}$ the unit normal pointing outward from the surface element $\textrm {d}S$. The $\alpha$-component $F_\alpha ^{\nu }$ ($\alpha =x, y,z$) of the viscous force vector $\boldsymbol {F}^{\nu }$ on the boundary is expressed by

(2.7)\begin{equation} F_{\alpha}^{\nu} = \boldsymbol{F}^{\nu}\boldsymbol{\cdot} \boldsymbol{e}_\alpha = 2 \nu \oint\limits_{\varGamma} \sum_{\beta=x,y,z} {S}_{\alpha,\beta}\,n_{\beta}\,\textrm{d}S , \end{equation}

where $\boldsymbol {e}_\alpha$ is the unit vector in the $\alpha$-direction and $S_{\alpha ,\beta } = ( \partial _{\alpha } u_{\beta } + \partial _{\beta } u_{\alpha } )/2$ the strain rate tensor with indices $\alpha , \beta = x, y,z$.

Similarly, the $\alpha$-component of the global pressure force, exerted on an immersed body, is defined as

(2.8)\begin{equation} F_{\alpha}^{p} = \boldsymbol{F}^{p}\boldsymbol{\cdot} \boldsymbol{e}_\alpha ={-} \oint\limits_{\varGamma} \,\textrm{d}S \, n_{\alpha}p . \end{equation}

Without external forces, the viscous and pressure forces in $\varOmega$ counter-balance the inertial terms provided by the left-hand side of (2.1). The drag force is defined as the projection on $\boldsymbol {e}_x$ of the pressure and viscous forces exerted on the body

(2.9)\begin{equation} F_D(t) = F_{x}^{p}(t) + F_{x}^{\nu}(t). \end{equation}

The lift force is similarly defined as the projection on $\boldsymbol {e}_y$ of the resulting pressure and viscous forces exerted on the body

(2.10)\begin{equation} F_L(t) = F_{y}^{p}(t) + F_{y}^{\nu}(t). \end{equation}

The drag and lift coefficients read

(2.11a,b)\begin{equation} C_D (t)= \frac{2F_D (t)}{\rho U^{2}},\quad C_L (t)= \frac{2F_L (t)}{\rho U^{2}}. \end{equation}

Employing the Galerkin approximation (2.4), the viscous force (2.7) can be re-written as

(2.12)\begin{equation} F_{\alpha}^{\nu} = \sum_{j=0}^{N} q^{\nu}_{\alpha; j} a_j, \end{equation}

where $q^{\nu }_{\alpha ; j}$ can easily be derived from (2.7) with the corresponding $S_{\alpha ,\beta }$ of the velocity mode $\boldsymbol {u}_j$, with the form

(2.13)\begin{equation} q^{\nu}_{\alpha; j} = 2 \nu \oint\limits_{\varGamma} \sum_{\beta=x,y,z} {S}_{\alpha,\beta}(\boldsymbol{u}_j) n_{\beta}\, \textrm{d}S . \end{equation}

Note that the contribution of the viscous force is linear with respect to the mode amplitudes $a_j$.

Similarly, from the pressure Poisson equation

(2.14)\begin{equation} \nabla^{2} p = \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} \otimes \boldsymbol{u} \right ) ={-} \sum_{\alpha=x,y,z}\sum_{\beta=x,y,z}\partial_{\alpha} u_{\beta} \partial_{\beta} u_{\alpha}, \end{equation}

the expression for the pressure field is derived as

(2.15)\begin{equation} p(\boldsymbol{x},t) = \sum_{j,k=0}^{N} p_{jk} (\boldsymbol{x}) a_j(t) a_k(t), \end{equation}

with

(2.16)\begin{equation} \nabla^{2} p_{jk}= \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_j \otimes \boldsymbol{u}_k \right )={-}\sum_{\alpha=x,y,z}\sum_{\beta=x,y,z}\partial_{\alpha} u_{\beta} (\boldsymbol{u}_j) \partial_{\beta} u_{\alpha}(\boldsymbol{u}_k). \end{equation}

The boundary conditions for partial pressures $p_{jk}$ are discussed by Noack et al. (Reference Noack, Papas and Monkewitz2005). Integrating (2.8) with (2.15) shows that the pressure force is a quadratic polynomial of the $a_j$ values

(2.17)\begin{equation} F_{\alpha}^{p} = \sum_{j,k=0}^{N} q^{p}_{\alpha; j k} a_j a_k, \quad \mbox{where} \ q^{p}_{\alpha; j k} ={-} \oint\limits_{\varGamma} \,\textrm{d}S \, n_{\alpha} p_{jk}. \end{equation}

Taking the steady solution as the basic mode $\boldsymbol {u}_0 = \boldsymbol {u}_s$ with $a_0 \equiv 1$ implies that $\boldsymbol {a}$ with $a_i = \delta _{0i}$ is a fixed point of (2.6) and the total force can be expressed as a constant-linear-quadratic expression in terms of the mode coefficients

(2.18)\begin{equation} F_{\alpha} = F_{\alpha}^{\nu} + F_{\alpha}^{p} = c_{\alpha} + \sum_{j=1}^{N} l_{\alpha; j} a_j + \sum_{j,k=1}^{N} q_{\alpha; j k} a_j a_k, \end{equation}

where

(2.19a–c)\begin{equation} c_{\alpha} = q^{\nu}_{\alpha; 0} + q^{p}_{\alpha; 0 0}, \quad l_{\alpha; j} = q^{\nu}_{\alpha; j} + q^{p}_{\alpha; j 0} + q^{p}_{\alpha; 0 j}, \quad q_{\alpha; j k} = q^{p}_{\alpha; j k}. \end{equation}

The force expression in (2.18) can be alternatively derived from the residual of the Navier–Stokes equations in the flow domain $\varOmega$, as demonstrated in appendix A.

With constant $\rho$ and $U$, the drag and lift coefficients in (2.11a,b) can be rewritten in the form

(2.20a)\begin{gather} C_D = c_{x} + \sum_{j=1}^{N} l_{x; j} a_j + \sum_{j,k=1}^{N} q_{x; j k} a_j a_k, \end{gather}

(2.20b)\begin{gather}C_L = c_{y} + \sum_{j=1}^{N} l_{y; j} a_j + \sum_{j,k=1}^{N} q_{y; j k} a_j a_k. \end{gather}

A crucial step relies on the choice of the $\boldsymbol {u}_i$ modes for the decomposition of (2.4). These could be the POD modes, as usually considered in fluid flows. However, a better choice could be to decompose the flow field on a basis of modes that are becoming active when the system is undergoing a bifurcation. This choice of the so-called bifurcation modes will be investigated in § 3.3.

2.3. The Navier–Stokes equations under the $Z_2$-symmetry

When the fluid flow configuration exhibits a mirror symmetry, the Navier–Stokes equations (2.1) possess at least one symmetric steady solution $(\boldsymbol {u}_s,p_s)$, satisfying (2.2). The $Z_2$-symmetry of the velocity and pressure fields, with respect to the ($x,z$)-plane defined by $y=0$, implies

(2.21a)\begin{gather} \left.\begin{gathered} u^{s}(x,-y,z) = u^{s}(x,y,z), \quad v^{s}(x,-y,z) ={-}v^{s}(x,y,z), \\ p^{s}(x,-y,z) = p^{s}(x,y,z), \end{gathered}\right\} \end{gather}

(2.21b)\begin{gather} \left.\begin{gathered} u^{a}(x,-y,z) ={-}u^{a}(x,y,z), \quad v^{a}(x,-y,z) = v^{a}(x,y,z), \\ p^{a}(x,-y,z) ={-}p^{a}(x,y,z), \end{gathered}\right\} \end{gather}

where the symmetric components $(u^{s}, v^{s}, p^{s}) \in \mathcal {U}^{s}$ and the antisymmetric components $(u^{a}, v^{a}, p^{a}) \in \mathcal {U}^{a}$, $\mathcal {U}^{s}$ and $\mathcal {U}^{a}$ being respectively the symmetric and antisymmetric subspaces of the system. Other steady solutions can exist, which break the symmetry of the system. We will consider the symmetric steady solution $(\boldsymbol {u}_s,p_s)$ as the reference point of (2.1) in the Reynolds decomposition of the flow field as (2.25).

The dynamics under consideration can include transient and post-transient regimes. Here, we introduce the $T$-averaged flow fields $\bar {\boldsymbol {u}}_T(\boldsymbol {x},t)$ as

(2.22)\begin{equation} \bar{\boldsymbol{u}}_T(\boldsymbol{x},t) = \frac{1}{T}\int_{t-T/2}^{t+T/2}\boldsymbol{u}(\boldsymbol{x},\tau)\,\textrm{d}\tau, \end{equation}

where $T$ is a time scale to be chosen. When the flow field is oscillating in time, an appropriate choice for $T$ is the period of the local oscillation. The mean flow field is further defined as

(2.23)\begin{equation} \bar{\boldsymbol{u}}(\boldsymbol{x}) = \lim _{T\rightarrow \infty}\bar{\boldsymbol{u}}_T(\boldsymbol{x},t) \end{equation}

and only focuses on the post-transient limit.

When two mirror-conjugated attractors co-exist, it is convenient to introduce the ensemble-averaged flow field $\bar {\boldsymbol {u}}_T^{\bullet } (\boldsymbol {x},t)$ as

(2.24)\begin{equation} \bar{\boldsymbol{u}}_T^{{\bullet}} (\boldsymbol{x},t) = \tfrac{1}{2}(\bar{\boldsymbol{u}}_T^{+}(\boldsymbol{x},t)+\bar{\boldsymbol{u}}_T^{-}(\boldsymbol{x},t)), \end{equation}

where $\bar {\boldsymbol {u}}_T^{\pm }(\boldsymbol {x},t)$ are the $T$-averaged flow field on the way to each individual attractor. This definition could be readily extended to more than two conjugated attractors. As an ensemble average on mirror-conjugated attracting sets, the ensemble-averaged flow field $\bar {\boldsymbol {u}}_T^{\bullet } (\boldsymbol {x},t)$ belongs to the symmetric subspace $\mathcal {U}^{s}$.

At this point, it is most convenient to introduce the Reynolds decomposition of the flow field, in the form

(2.25)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \bar{\boldsymbol{u}}_T^{{\bullet}} (\boldsymbol{x},t) + \boldsymbol{u}'(\boldsymbol{x},t) = \boldsymbol{u}_s(\boldsymbol{x}) + \boldsymbol{u}_{\varDelta}(\boldsymbol{x},t) + \boldsymbol{u}'(\boldsymbol{x},t), \end{equation}

where the mean-field deformation $\boldsymbol {u}_{\varDelta }(\boldsymbol {x},t)$ accounts for the distortion of the flow field from the symmetric steady solution $\boldsymbol {u}_s(\boldsymbol {x})$ to the ensemble-averaged flow field $\bar {\boldsymbol {u}}_T^{\bullet } (\boldsymbol {x},t)$ as

(2.26)\begin{equation} \boldsymbol{u}_{\varDelta}(\boldsymbol{x},t)= \bar{\boldsymbol{u}}_T^{{\bullet}} (\boldsymbol{x},t) - \boldsymbol{u}_s(\boldsymbol{x}). \end{equation}

The fluctuation flow field $\boldsymbol {u}'(\boldsymbol {x},t)$ has a vanishing time average, meaning that $\boldsymbol {u}(\boldsymbol {x},t)$ is centred on $\bar {\boldsymbol {u}}_T^{\bullet } (\boldsymbol {x},t)$. By construction, $\bar {\boldsymbol {u}}_T^{\bullet } (\boldsymbol {x},t), \boldsymbol {u}_{\varDelta }(\boldsymbol {x},t), \boldsymbol {u}_s(\boldsymbol {x})$ belongs to the symmetric subspace $\mathcal {U}^{s}$ and $\boldsymbol {u}'(\boldsymbol {x},t)$ to the antisymmetric subspace $\mathcal {U}^{a}$. Thus, a symmetry-based decomposition of (2.1) results in a symmetric and an antisymmetric part, yielding

(2.27a)\begin{gather} \partial_t \boldsymbol{u}_{\varDelta} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ \boldsymbol{u}_s \otimes \boldsymbol{u}_{\varDelta} + \boldsymbol{u}_{\varDelta} \otimes \boldsymbol{u}_s + \boldsymbol{u}_{\varDelta} \otimes \boldsymbol{u}_{\varDelta} + \boldsymbol{u}^{\prime} \otimes \boldsymbol{u}^{\prime} \right] = \nu \triangle \boldsymbol{u}_{\varDelta} - \boldsymbol{\nabla} p_{\varDelta}, \end{gather}

(2.27b)\begin{gather}\partial_t \boldsymbol{u}^{\prime} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ \bar{\boldsymbol{u}}_T^{{\bullet}} \otimes \boldsymbol{u}^{\prime} + \boldsymbol{u}^{\prime} \otimes \bar{\boldsymbol{u}}_T^{{\bullet}} \right] = \nu \triangle \boldsymbol{u}^{\prime} - \boldsymbol{\nabla} p^{\prime}. \end{gather}

Integrating (2.27a) on the spatial domain $\varOmega$, both the left- and right-hand sides yield a time-evolving force vector aligned on $\boldsymbol {e}_y$, while integrating (2.27b) yields a time-evolving force vector aligned on $\boldsymbol {e}_x$. The former is the resulting lift force applying to the boundaries of the fluid domain, while the latter is the drag force. Thus, the $Z_2$-symmetry applied to equations (2.20a) and (2.20b) yields

(2.28a)\begin{gather} C_D = C_D^{{\circ}} + \sum_{j=1}^{N} \underbrace{ [l_{x; j} a_j]}_{{\in} \mathcal{U}^{s}} + \sum_{j,k=1}^{N} \underbrace{ [q_{x; j k} a_j a_k]}_{{\in} \mathcal{U}^{s}}, \end{gather}

(2.28b)\begin{gather}C_L = \sum_{j=1}^{N} \underbrace{ [l_{y; j} a_j]}_{{\in} \mathcal{U}^{a}} + \sum_{j,k=1}^{N} \underbrace{ [q_{y; j k} a_j a_k]}_{{\in} \mathcal{U}^{a}}, \end{gather}

where $C_D^{\circ }$ is the drag coefficient of the symmetric steady solution.

The vanishing terms in (2.28) can be easily derived from the definitions of $q^{\nu }_{\alpha ; j}$ and $q^{p}_{\alpha ; j k}$ in § 2.2 as

(2.29a)\begin{gather} l_{x; j} = q_{x;j}^{\nu} + q_{x; 0j}^{p} + q_{x; j0}^{p} = 0, \quad \boldsymbol{u}_j \in \mathcal{U}^{a} , \end{gather}

(2.29b)\begin{gather}l_{y; j} = q_{y;j}^{\nu} + q_{y; 0j}^{p} + q_{y; j0}^{p} = 0, \quad \boldsymbol{u}_j \in \mathcal{U}^{s} , \end{gather}

(2.29c)\begin{gather}q_{x; j k} = q_{x; j k}^{p} = 0, \quad \boldsymbol{u}_j \otimes \boldsymbol{u}_k \in \mathcal{U}^{a} , \end{gather}

(2.29d)\begin{gather}q_{y; j k} = q_{y; j k}^{p} = 0, \quad \boldsymbol{u}_j \otimes \boldsymbol{u}_k \in \mathcal{U}^{s} . \end{gather}

As a result, the drag contribution must come from the symmetric subsets of the constant-linear-quadratic polynomial functions, and from the antisymmetric subsets for the lift contribution.

3.1. The fluidic pinball

The geometric configuration, shown in figure 1, consists of three fixed cylinders of unit diameter $D$ mounted on the vertices of an equilateral triangle of side length $3D/2$ in the $(x, y)$ plane. The flow domain is bounded with a $[-6,+20]\times [-6,+6]$ box. The upstream flow, of uniform velocity $U_\infty$ at the input of the domain, is transverse to the cylinder axis and aligned with the symmetry axis of the cylinder cluster. All quantities will be non-dimensionalized with the cylinder diameter $D$, the velocity $U_\infty$ and the unit fluid density $\rho$. Considering the symmetry of this configuration, a Cartesian coordinate system will be used in the following discussion, with its origin in the middle of the rightmost two cylinders. In this study, no external force will be applied to these three cylinders. A no-slip condition is applied on the cylinders and the velocity in the far wake is assumed to be $U_\infty$. Here, the Reynolds number is defined as $Re = U_\infty D/\nu$, where $\nu$ is the kinematic viscosity of the fluid. A no-stress condition is applied at the output of the domain.

Figure 1. Configuration of the fluidic pinball and dimensions of the simulated domain. A typical field of vorticity is represented in colourwith $[-1.5, 1.5]$. The upstream velocity is denoted $U_\infty$.

The resolution of the Navier–Stokes equations (2.1) is based on a second-order finite-element discretization method of the Taylor–Hood type (Taylor & Hood Reference Taylor and Hood1973), on an unstructured grid of 4225 triangles and 8633 vertices and an implicit integration of the third-order in time. The instantaneous flow field is calculated with a Newton–Raphson iteration until the residual reaches a tiny tolerance prescribed. This approach is also used to calculate the steady solution, which is derived from the steady Navier–Stokes equations (2.2). The direct Navier–Stokes solver used herein has been validated in Noack et al. (Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003) and Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020), with a detailed technical report (Noack & Morzyński Reference Noack and Morzyński2017). The grid used for the simulations was shown to provide a consistent flow dynamics, compared to a refined grid, see Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020).

3.2. Flow features and the corresponding force dynamics

Different from Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020), where the viscous contribution to the forces has been ignored, the lift $C_L$ and drag $C_D$ coefficients are here calculated from the resulting force $\boldsymbol {F}$ of pressure and viscous components exerted on the three cylinders.

The flow characteristics depend on the Reynolds number $Re$. Following the literature on clusters of cylinders (Chen et al. Reference Chen, Ji, Alam, Williams and Xu2020), the characteristic length scale is chosen to be the cylinder diameter $D$ and not the transverse width $5 D/2$ of the configuration. This width loses its dynamic significance for the large distances considered in other studies.

For Reynolds numbers $Re < Re_{H}\approx 18$, the symmetric steady solution $\boldsymbol {u}_s(\boldsymbol {x})$ was found to be stable and is the only attractor of the system. A supercritical Hopf bifurcation occurs at $Re = Re_{H}$, associated with the cyclic release of counter-rotating vortices in the wake of the three cylinders from the shear layers that delimit the configuration, forming a von Kármán street of vortices. The corresponding Reynolds number based on the transverse width of the fluidic pinball is $45$, i.e. is well aligned with typical onsets of vortex shedding behind bluff bodies. For the critical value $Re = Re_{PF} \approx 68$, the system undergoes a supercritical pitchfork bifurcation. As a result, two additional (unstable) steady solutions occur, namely $\boldsymbol {u}_s^{+}(\boldsymbol {x})$ and $\boldsymbol {u}_s^{-}(\boldsymbol {x})$, which break the reflectional symmetry of the configuration, as shown with the lift coefficients of the steady solutions in figure 2. The mean field inherits the asymmetry of the steady solutions, with the jet between the two downstream cylinders being deflected upward or downward. As reported in Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020), at $Re=Re_{PF}$, the statistically symmetric limit cycle, associated with the statistically symmetric vortex shedding, becomes unstable with respect to two mirror-conjugated statistically asymmetric limit cycles, associated with statistically asymmetric von Kármán streets of vortices.

Figure 2. Lift coefficients at different Reynolds numbers (a) of the symmetric steady solutions $\bar {\boldsymbol {u}}_s$ (black curve), the asymmetric steady solutions $\bar {\boldsymbol {u}}_s^{-}$ (blue curve), the asymmetric steady solutions $\bar {\boldsymbol {u}}_s^{+}$ (red curve), exemplified with the vorticity field of $\bar {\boldsymbol {u}}_s^{+}$, $\bar {\boldsymbol {u}}_s$, $\bar {\boldsymbol {u}}_s^{-}$ at $Re = 100$ from top to bottom (b).

Figure 3 shows the time evolution of the lift and drag coefficients at $Re=80$, when the initial condition is either the symmetric steady solution $\boldsymbol {u}_s$ (figure 3a) or the asymmetric steady solution $\boldsymbol {u}_s^{+}$ (figure 3b). In both cases, the asymptotic regime is the same. However, when starting from the symmetric steady solution in figure 3(a), a long-living plateau of the drag coefficient is reached around $t\approx 775$, which corresponds to the transient exploration of the unstable limit cycle, centred on the symmetric $T$-averaged flow field $\bar {\boldsymbol {u}}_{98}(\boldsymbol {x},775)$. Note that, during the transient dynamics from the steady solution to the unstable limit cycle, the drag coefficient is monotonically increasing, before reaching the transient plateau. The drag coefficient is further increased when leaving the unstable limit cycle towards the asymptotically stable limit cycle, the latter being centred on the asymmetric mean flow field $\bar {\boldsymbol {u}}^{+}$.

Figure 3. Time evolution of the drag (a,b) and lift (c,d) coefficients, starting (a,c) from the symmetric steady solution $\boldsymbol {u}_s$, (b,d) from the asymmetric steady solution $\boldsymbol {u}_s^{+}$, at $Re=80$.

Figure 4 shows another representation of the transient dynamics for $Re=30$, 80 and 100, starting from different initial conditions in the plane ($C_L,\Delta C_D$), where $\Delta C_D = C_D - C_D^{\circ }$, $C_D^{\circ }$ being the drag associated with the symmetric steady solution at the Reynolds number under consideration. The black cross ($\times$) stands for the symmetric steady solution $\boldsymbol {u}_s$ while the asymmetric $\boldsymbol {u}_s^{+}$ and $\boldsymbol {u}_s^{-}$ steady solutions are respectively represented by a red circle and a blue square, when they exist, at $Re=80$ and 100. As can be observed in this figure, different from what happens at $Re=80$, the transient dynamics from the symmetric steady solution at $Re=100$ first reaches one of the two asymmetric steady solutions, before evolving toward the stable attracting limit cycle.

Figure 4. Trajectories in the $(C_L,\Delta C_D)$ plane, for Reynolds numbers $Re=30$, 80 and 100, starting, for the black trajectories, close to the symmetric steady solution $\boldsymbol {u}_s$ ($\times$), for the red trajectories close to the asymmetric steady solution $\boldsymbol {u}_s^{+}$ ($\bullet$) and for the blue trajectories close to the asymmetric steady solution $\boldsymbol {u}_s^{-}$ ($\blacksquare$). Here, $\Delta C_D = C_D - C_D^{\circ }$, where $C_D^{\circ }$ is the drag coefficient of the symmetric steady solution at the corresponding Reynolds number.

3.3. The bifurcation modes of the fluidic pinball

In the case of two subsequent supercritical Hopf and pitchfork bifurcations, Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020) have shown that the reduced-order model must comprise 5 modes

(3.1)\begin{equation} \boldsymbol{u} (\boldsymbol{x},t ) = \boldsymbol{u}_s ( \boldsymbol{x} ) + \sum_{j=1}^{5} a_j(t) \boldsymbol{u}_j ( \boldsymbol{x} ). \end{equation}

Hence, in the decomposition of (2.4), the number of modes is restricted to $N=5$. For dynamic interpretability, the basic mode $\boldsymbol {u}_0(\boldsymbol {x})$ is chosen to be the symmetric steady solution $\boldsymbol {u}_s(\boldsymbol {x})$. The first three modes $\boldsymbol {u}_{1,2,3}(\boldsymbol {x})$ are associated with the Hopf bifurcation, the last two modes $\boldsymbol {u}_{4,5}(\boldsymbol {x})$ with the pitchfork bifurcation. We will refer to these modes as the irreducible bifurcation modes of the system. Modes $\boldsymbol {u}_3(\boldsymbol {x})$ and $\boldsymbol {u}_5(\boldsymbol {x})$ are symmetric. The instability-related modes $\boldsymbol {u}_{1,2}(\boldsymbol {x})$ and $\boldsymbol {u}_4(\boldsymbol {x})$ are antisymmetric. Modes $\boldsymbol {u}_{1,2}(\boldsymbol {x})$ span the subspace associated with the limit cycle of the Hopf bifurcation, while $\boldsymbol {u}_4(\boldsymbol {x})$ accounts for the symmetry breaking of the pitchfork bifurcation. In Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020), modes $\boldsymbol {u}_{1,2}(\boldsymbol {x})$ are provided by the first two dominant POD modes, while mode $\boldsymbol {u}_4(\boldsymbol {x})$ is defined as

(3.2)\begin{equation} \boldsymbol{u}_4(\boldsymbol{x}) \propto \boldsymbol{u}_s^{+}(\boldsymbol{x})-\boldsymbol{u}_s^{-}(\boldsymbol{x}), \end{equation}

where $\boldsymbol {u}_s^{\pm }(\boldsymbol {x})$ are the two additional (asymmetric) steady solutions arising from the supercritical pitchfork bifurcation. Mode $\boldsymbol {u}_3(\boldsymbol {x})$ is slaved to $\boldsymbol {u}_{1,2}(\boldsymbol {x})$ while $\boldsymbol {u}_5(\boldsymbol {x})$ is slaved to $\boldsymbol {u}_4(\boldsymbol {x})$. The mode $\boldsymbol {u}_3(\boldsymbol {x})$ is usually defined as the shift mode from $\boldsymbol {u}_s(\boldsymbol {x})$ to the asymptotic mean flow field, $\boldsymbol {u}_3(\boldsymbol {x})\propto \bar {\boldsymbol {u}}(\boldsymbol {x})-\boldsymbol {u}_s(\boldsymbol {x})$, before being ortho-normalized to $\boldsymbol {u}_1(\boldsymbol {x})$ and $\boldsymbol {u}_2(\boldsymbol {x})$. Here, $\bar {\boldsymbol {u}}(\boldsymbol {x})$ will be restricted to the symmetric mean flow field, associated with the statistically symmetric limit cycle, whether this limit cycle is stable or unstable. Similarly to $\boldsymbol {u}_4(\boldsymbol {x})$, mode $\boldsymbol {u}_5(\boldsymbol {x})$ is defined as

(3.3)\begin{equation} \boldsymbol{u}_5(\boldsymbol{x}) \propto (\boldsymbol{u}_s^{+}(\boldsymbol{x})+\boldsymbol{u}_s^{-}(\boldsymbol{x}))-2\boldsymbol{u}_s(\boldsymbol{x}). \end{equation}

These two modes, together with modes $\boldsymbol {u}_{1,2,3}$, are shown in figure 5 after ortho-normalization by a Gram–Schmidt procedure, and the corresponding time-dependent amplitudes $a_i(t)$, $i=1,\ldots ,5$, in the full-flow dynamics are shown in figure 6 when starting from either the symmetric steady solution (figure 6a) or the asymmetric steady solution (figure 6b).

Figure 5. Vortical structure (colour) of the modes $\boldsymbol {u}_1(\boldsymbol {x})$, $\boldsymbol {u}_2(\boldsymbol {x})$, $\boldsymbol {u}_3(\boldsymbol {x})$ (a–c), $\boldsymbol {u}_4(\boldsymbol {x})$, $\boldsymbol {u}_5(\boldsymbol {x})$ (d,e), of the velocity field associated with the five elementary degrees of freedom $\{a_1(t) - a_5(t)\}$, at $Re=80$.

Figure 6. Mode amplitudes $a_i(t)$, $i=1,\dots ,5$ in the full-flow dynamics starting (a) from the symmetric steady solution $\boldsymbol {u}_s$, (b) from the asymmetric steady solution $\boldsymbol {u}_s^{+}$, at $Re=80$.

4.1. Force model associated with the supercritical Hopf bifurcation

The symmetric steady solution $\boldsymbol {u}_s \in \mathcal {U}^{s}$ is stable at low Reynolds numbers. At $Re \geqslant Re_{HP}$, it undergoes a supercritical Hopf bifurcation. The resulting Galerkin expansion reads

(4.1)\begin{equation} \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}_{\varDelta}}, \end{equation}

and the corresponding mean-field Galerkin system

(4.2a)\begin{gather} \textrm{d}a_1/\textrm{d}t = \sigma a_1 - \omega a_2, \end{gather}

(4.2b)\begin{gather}\textrm{d}a_2/\textrm{d}t = \sigma a_2 + \omega a_1, \end{gather}

(4.2c)\begin{gather}\textrm{d}a_3/\textrm{d}t = \sigma_3 a_3 + \beta_3 \left( a_1^{2} + a_2^{2} \right), \end{gather}

with $\sigma = \sigma _1 - \beta a_3$ and $\omega = \omega _1 + \gamma a_3$, where $\sigma _1$ and $\omega _1$ are the initial growth rate and frequency depending on the Reynolds number. For a direct supercritical Hopf bifurcation, $\sigma _1, \omega _1, \beta > 0$, $\sigma _3 < 0$ and $\beta _3>0$. We refer to Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020) for details.

Introducing (4.1) in (2.7) and (2.8), the total force can be written as (2.18) with $N=3$ degrees of freedom. From symmetry considerations, as $\boldsymbol {u}_{1,2} \in \mathcal {U}^{a}$ and $\boldsymbol {u}_{0,3} \in \mathcal {U}^{s}$, the coefficients $l_{x; 1}$, $l_{x; 2}$, $q_{x; 13}$, $q_{x; 23}$, $l_{y; 0}$, $l_{y; 3}$, $q_{y; 11}$, $q_{y; 12}$, $q_{y; 22}$, $q_{y; 33}$ are vanishing. Finally, the drag formulae (2.28) simplify to

(4.3a)\begin{gather} C_D = C_D^{{\circ}} + l_{x; 3} a_3 + q_{x; 1 1} a_1^{2} + q_{x; 1 2} a_1 a_2 + q_{x; 2 2} a_2^{2} + q_{x; 3 3} a_3^{2}, \end{gather}

(4.3b)\begin{gather}C_L = l_{y; 1} a_1 + l_{y; 2} a_2 + q_{y; 1 3} a_1 a_3 + q_{y; 2 3} a_2 a_3. \end{gather}

Here again, $C_D^{\circ }$ is the drag coefficient associated with the symmetric steady solution. The unknown parameters in the force model can be identified by a least-squares approach, according to the known force dynamics and the relevant mode amplitudes. However, for the mean-field Galerkin system (4.2), the slaving relation between the degree of freedom $a_3$ to the oscillating degrees of freedom $a_{1}$, $a_{2}$ imposes an additional sparsity in the force model. We employ the SINDy (sparse identification of nonlinear dynamics) algorithm (Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016) to arrive at simpler and more interpretable models. A $L1$-regularization can be introduced in the LASSO (least absolute shrinkage and selection operator) regression process, which includes the L1-norm of the vector of coefficients in calculating the loss function. Another option in the SINDy algorithm is the sequential thresholded least-squares regression, which iteratively applies the least-squares regression and eliminates terms with weight smaller than a given threshold. Both regression algorithms benefit from simplicity, only requiring one sparsity parameter $\lambda$. The optimal $\lambda$ balances the accuracy and complexity of the identified model. To evaluate the performance of the identified model, the complexity is presented with the number of non-zero coefficients and the accuracy by the coefficient of determination,denoted as the $r^{2}\ score$ (Draper & Smith Reference Draper and Smith1998). A detailed review of this sparsity parameter can be found in Loiseau & Brunton (Reference Loiseau and Brunton2018). A recent extension of the SINDy algorithm with physical constraints of energy-preserving quadratic nonlinearities successfully identifies the sparse model, benefiting from the Galerkin projection of the Navier–Stokes equations (Loiseau, Noack & Brunton Reference Loiseau, Noack and Brunton2018).

The LASSO algorithm is applied to a scenario starting with the unstable symmetric steady solution at $Re=30$. The training data used for the sparse regression are provided by the force coefficients and the mode amplitudes from the direct numerical simulation (DNS) starting with the symmetric steady solution to the final asymptotic regime. The resulting transient dynamics and the asymptotically attracting limit cycle are shown in the three-dimensional space of the time-delayed coordinates of $C_L$ and $C_D$ in figure 7.

Figure 7. Transient dynamics from the unstable symmetric steady solution $\boldsymbol {u}_s$ ($\times$) to the asymptotic limit cycle (statistically symmetric vortex shedding), at $Re=30$, in the time-delayed embedding space of the lift $C_L$ and drag $C_D$ coefficients, with $\tau =2$.

The possible over-fitting terms, such as the slaving relation between $a_3$ and $a_1^{2}, a_2^{2}$, can be suppressed with a larger $L1$-penalty parameter for the LASSO algorithm. The choice of the $L1$-penalty parameter drives the sparsity of the identified model. A too small $L1$ will lead to a complex model with few eliminated terms; on the contrary, a too large $L1$ can jeopardize accuracy. Both cases weaken the robustness of the identified model, and the same is observed for the sequential thresholded least-squares regression. The influence of the sparsity parameter $\lambda$ and the comparison of these two regression methods are presented in appendix B.

Gradually increasing the $L1$-penalty from $0$ to nearly $1$, the terms $a_1 a_2$, $a_3$, $a_2^{2}$, $a_1^{2}$ are eliminated subsequently in the drag model, while $a_3^{2}$ is always retained. The sparsity parameter $\lambda$, here the $L1$-penalty, is chosen as the largest value without any known over-fitting term. Hence, according to the order of elimination, $a_3$ is the over-fitting term in the drag model due to the slaving relation between $a_3$ and $a_1^{2}, a_2^{2}$. The details of this choice can be found in appendix B. Finally, the identified force model reads

(4.4a)\begin{gather} C_D = 4.82440448 - 0.00037484 a_1^{2} - 0.00098337 a_2^{2} + 0.01777408 a_3^{2}, \end{gather}

(4.4b)\begin{gather}C_L = 0.00867623 a_1 + 0.01397362 a_2 + 0.0166239 a_1 a_3 - 0.01302317 a_2 a_3. \end{gather}

The force model is highly accurate, as corrobororated by the $r^{2}$ scores of $0.9991$ and $0.9942$ for the drag and lift formulae, respectively. As shown in figure 8, the dynamics of the force model compares well with the real force transient dynamics, starting from the symmetric steady solution at $Re=30$.

Figure 8. Performance of the force model with the three elementary modes of the Hopf bifurcation. Time evolution of the drag $C_D$ (a) and lift $C_L$ (b) coefficients, in the full-flow dynamics (solid black line) and for the force model (red dashed line), at $Re=30$. Initial condition: symmetric steady solution.

In the drag model (4.4a), the coefficient of $a_3$ is vanishing. Mode $\boldsymbol {u}_3$ actually contributes to the increase of the drag through $a_3^{2}$, as evidenced by the positive coefficient of the $a_3^{2}$ term. This is an interesting result, since the effect of the bifurcation mode $\boldsymbol {u}_3$ is to decrease the length of the recirculation bubble in the $T$-averaged flow field $\bar {\boldsymbol {u}}_T(\boldsymbol {x},t)\approx \boldsymbol {u}_s(\boldsymbol {x})+a_3(t)\boldsymbol {u}_3(\boldsymbol {x})$, resulting in an increase of the drag through the quadratic term $a_3^{2}$. This quadratic dependency is also reported in Loiseau et al. (Reference Loiseau, Noack and Brunton2018).

It is also worth noticing that $a_3^{2}$ contributes to the mean value of $C_D$ while $a_1^{2}, a_2^{2}$ accounts for the instantaneous oscillations of $C_D$, as $C_D$ oscillates at twice the vortex shedding frequency. For $C_L$, the oscillatory pair $(a_1 , a_2)$ fits well with the phase of the initial transient part, while the pair $(a_1 a_3 , a_2 a_3)$ resolves the phase dependency of the post-transient part of the dynamics.

4.2. Force model associated with the supercritical pitchfork bifurcation

Next, we consider the supercritical pitchfork bifurcation, which breaks the symmetry of the symmetric steady solution $\boldsymbol {u}_s$ at $Re \geqslant Re_{PF}$. In this case the antisymmetric mode $\boldsymbol {u}_{4}$ describes the antisymmetric instability, which corresponds to an unstable eigenmode with a real eigenvalue. The resulting Galerkin expansion reads

(4.5)\begin{equation} \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}_{\varDelta}}, \end{equation}

and the corresponding mean-field Galerkin system

(4.6a)\begin{gather} \textrm{d}a_4/\textrm{d}t = \sigma_4 a_4 - \beta_4 a_4 a_5, \end{gather}

(4.6b)\begin{gather}\textrm{d}a_5/\textrm{d}t = \sigma_5 a_5 + \beta_5 a_4^{2}, \end{gather}

where $\sigma _4$ is the positive initial growth rate, which depends on the Reynolds number. For a direct supercritical pitchfork bifurcation, $\sigma _4, \beta _4 > 0$, $\sigma _5 < 0$ and $\beta _5>0$, see Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020) for details.

Substituting (4.5) in (2.7) and (2.8), with $N=2$ in (2.18), and with $\boldsymbol {u}_{4} \in \mathcal {U}^{a}$ and $\boldsymbol {u}_s, \boldsymbol {u}_5 \in \mathcal {U}^{s}$, the force model becomes

(4.7a)\begin{gather} C_D = C_D^{{\circ}} + l_{x; 5} a_5 + q_{x; 4 4} a_4^{2} + q_{x; 5 5} a_5^{2}, \end{gather}

(4.7b)\begin{gather}C_L = l_{y; 4} a_4 + q_{y; 4 5} a_4 a_5. \end{gather}

Five parameters, namely $l_{x; 0}$, $l_{x; 5}$, $q_{x; 44}$, $q_{x; 55}$, $l_{y; 4}$, $q_{y; 45}$ need to be identified.

In the fluidic pinball, the pitchfork bifurcation occurs after the primary Hopf bifurcation as the Reynolds number is increased. However, the transient dynamics observed at $Re=100$, when starting close to the symmetric steady solution, first exhibits the static symmetry breaking, which is typical of the pitchfork bifurcation, before developing the cyclic release of vortices, which is characteristic of the Hopf bifurcation. The early stage of the transient dynamics, starting from the symmetric steady solution and evolving toward one of the asymmetric steady solutions, is shown in figure 9. The time evolutions of the lift $C_L(t)$ and drag $C_D(t)$ coefficients are shown in figure 10.

Figure 9. Transient trajectories (solid and dashed lines) starting from two initial conditions close to the symmetric steady solution, at $Re=100$. Asymmetric steady solution $\boldsymbol {u}_s^{+}$ ($\bullet$), asymmetric steady solution $\boldsymbol {u}_s^{-}$ ($\blacksquare$).

Figure 10. Performance of the force model with the two elementary modes of the pitchfork bifurcation. Time evolution of the drag $C_D$ (a) and lift $C_L$ (b) coefficients in the full-flow dynamics (solid black line) and for the force model (red dashed line), at $Re=100$. Initial condition: symmetric steady solution.

Only the degrees of freedom associated with the pitchfork bifurcation are active in this early stage of the transient dynamics, as also shown in figure 17(a). The degrees of freedom associated with the Hopf bifurcation will only become active further in time during the transient dynamics, which will be further discussed in § 5.4. Accordingly, a force model is derived for the transition after a simple pitchfork bifurcation. The training data are the lift $C_L(t)$ and drag $C_D(t)$ coefficients and the relevant mode amplitudes in (4.7) from the early to final stage of the transient dynamics. The observed slaving of $a_5$ in $a_4^{2}$ may reduce the robustness of the identified model. Gradually increasing the $L1$-penalty parameter in the LASSO regression, the optimized force model reads

(4.8a)\begin{gather} C_D = 3.58248992 + 0.04367604 a_5 - 0.08525184 a_5^{2}, \end{gather}

(4.8b)\begin{gather}C_L ={-}0.13611053 a_4 + 0.09194312 a_4 a_5, \end{gather}

with $r^{2} =0.9949$ for the drag model and $r^{2} =0.9992$ for the lift model. The over-fitting term $a_4^{2}$ has been eliminated in the sparse formula of the drag force. Note that the mode $\boldsymbol {u}_5$ contributes to the drag through $a_5$, while $a_5^{2}$ acts in decreasing the drag, as indicated by the sign of their associated coefficients in (4.8a).

Figure 10 compares the evolutions of the drag and lift coefficients in the full-flow dynamics (solid black line) to their prediction by the force model (4.8) (red dashed curve), during the early stage of the transient dynamics at $Re=100$. The derived force model is well aligned with the real force dynamics using only two active degrees of freedom of the pitchfork bifurcation in the dynamics of the system.

5.1. Force model at $Re=80$

At $Re=80$, the system has already undergone a supercritical Hopf bifurcation and a supercritical pitchfork bifurcation. The trajectories issued from $\boldsymbol {u}_s$ and $\boldsymbol {u}_s^{\pm }$ are shown in the time-delayed embedding state space $(C_L(t),C_L(t-\tau ),C_D(t))$ of figure 11. The force model will rely on five degrees of freedom at minimum, namely the three degrees of freedom associated with the Hopf bifurcation $a_i$, $i=1,2,3$ and the two degrees of freedom $a_i$, $i=4,5$, associated with the pitchfork bifurcation. As a generalization of (4.3) and (4.7), the force model reads

(5.1a)\begin{align} C_D & = C_D^{{\circ}} + l_{x; 3} a_3 + q_{x; 1 1} a_1^{2} + q_{x; 1 2} a_1 a_2 + q_{x; 2 2} a_2^{2} + q_{x; 3 3} a_3^{2} \nonumber\\ & \quad + l_{x; 5} a_5 + q_{x; 4 4} a_4^{2} + q_{x; 5 5} a_5^{2} \nonumber\\ & \quad + q_{x; 1 4} a_1 a_4 + q_{x; 2 4} a_2 a_4 + q_{x; 3 5} a_3 a_5 , \end{align}

(5.1b)\begin{align} C_L & = l_{y; 1} a_1 + l_{y; 2} a_2 + q_{y; 1 3} a_1 a_3 + q_{y; 2 3} a_2 a_3 \nonumber\\ &\quad + l_{y; 4} a_4 + q_{y; 4 5} a_4 a_5 \nonumber\\ & \quad + q_{y; 1 5} a_1 a_5 + q_{y; 2 5} a_2 a_5 + q_{y; 3 4} a_3 a_4. \end{align}

Due to symmetry reasons, only 2 linear terms ($a_{3},a_{5}$) and 9 quadratic terms ($a_1^{2}$, $a_1a_2$, $a_1a_4$, $a_2^{2}$, $a_2a_4$, $a_3^{2}$, $a_3a_5$, $a_4^{2}$, $a_5^{2}$) are left in (5.1a) for the drag coefficient. For the lift coefficient, only 3 linear terms, $a_{1}$, $a_{2}$, $a_{4}$, and 6 quadratic terms, $a_1a_3$, $a_1a_5$, $a_2a_3$, $a_2a_5$, $a_3a_4$, $a_4a_5$, are left in (5.1b). The training data are taken from the DNS starting from the three steady solutions, with the real force dynamics, see the black curves in figure 12, and the relevant mode amplitudes, see figure 6. The coefficients of the force models are identified by the sequential thresholded least-squares regression with the optimal sparsity parameter $\lambda$. We note that the LASSO regression can also be used here. See appendix B for the comparison of these two methods. The resulting force model reads

(5.2a)\begin{align} C_D &= 3.77331204 + 0.05888312 a_5 - 0.01115552 a_1^{2} - 0.01088109 a_2^{2}\nonumber\\ &\quad + 0.01323449 a_3^{2} + 0.02949701 a_3a_5 - 0.25910470 a_5^{2} , \end{align}

(5.2b)\begin{align} C_L &= 0.00953160 a_1 + 0.00720164 a_2 - 0.10179203 a_4 \nonumber\\ &\quad - 0.00303677 a_1a_3 - 0.00197075 a_2 a_3 - 0.00200840 a_3a_4 \nonumber\\ &\quad + 0.05914386 a_4a_5. \end{align}

The good accuracy of the identified drag model can be determined from the high $r^{2}$ score of $0.9816$. The drag model of (5.2a) preserves both the basic forms of the drag model for the Hopf and pitchfork bifurcations and the signs of the coefficients. This indicates that the identified model is robust. The only remaining cross-term $a_3a_5$ provides the coupling relation between the degrees of freedom associated with both bifurcations.

Figure 11. Trajectories in the time-delayed embedding space of the lift $C_L$ and drag $C_D$ coefficients, with $\tau =2$, at $Re=80$. Black trajectories starting close to the symmetric steady solution $\boldsymbol {u}_s$ ($\times$); red trajectory starting close to the asymmetric steady solution $\boldsymbol {u}_s^{+}$ ($\bullet$); blue trajectory starting close to the asymmetric steady solution $\boldsymbol {u}_s^{-}$ ($\blacksquare$).

Figure 12. Performance of the force model with the five elementary modes. Time evolution of the drag $C_D$ (a,b) and the lift $C_L$ (c,d) coefficients in the full-flow dynamics (solid black line) and for the force model (red dashed line), at $Re=80$. Initial condition: (a,c) symmetric steady solution, (b,d) asymmetric steady solution.

A robust sparse formula for the lift model is more difficult to derive, due to the oscillating dynamics of the lift and the fact that $a_4$ and $a_5$ also oscillate at the fundamental frequency. With respect to the basic lift model of two bifurcations, a balanced method is used here to solve the difficulty of the identification. Starting with a large $L1$-penalty, the derived under-fitted system can figure out the most elementary features of the dynamics, eliminating $a_1a_5$, $a_2a_5$, $a_4a_5$. This is reasonable as $a_3$ is approximately ten times larger than $a_5$, which means that most of the mean-field distortion comes from $\boldsymbol {u}_3$. However, if the $L1$-penalty is too large, the term $a_4a_5$ can disappear from the lift model, making the resulting model inconsistent with (5.2b). In order to balance sparsity and robustness, $a_4a_5$ needs to be reintroduced into the library. The sparse formula of the lift model in (5.2b) is determined by least-squares regression, constraining the parameters of $a_1a_5$, $a_2a_5$ to zero. The $r^{2}$ score of the identified lift model is $0.9673$.

The identified force dynamics in (5.2) (dashed red line) is compared to the real force dynamics (solid black line) at $Re=80$ in figure 12. The force model based on the least-order model can reproduce the main features of the real force dynamics. The drag model of (5.2a) shows how the degrees of freedom of the Hopf ($a_1^{2}, a_2^{2}, a_3^{2}$) and pitchfork ($a_5, a_5^{2}$) bifurcations contribute to the drag force, as well as the coupling between these degrees of freedom ($a_3a_5$). The lift model of (5.2b) shows that the lift oscillations occur through the coupling of the oscillating degrees of freedoms $a_{1}$, $a_{2}$ to $a_3$, while the coupling between the degrees of freedoms $a_4$ and $a_5$ contribute to the mean value of $C_L$. Hence, the mean lift coefficient can be simplified with fewer terms, as $\bar {C_L} = l_{y; 4} a_4 + q_{y; 4 5} a_4 a_5 + q_{y; 3 4} a_3 a_4$, which is consistent with the Krylov–Bogoliubov assumption (Jordan & Smith Reference Jordan and Smith1999).

5.2. Assessing the predictive power of the force model

The time evolution of the drag and lift coefficients in the fluidic pinball are shown in figure 12 as solid black lines. The evolutions of the drag and lift coefficients in the model (5.2) are shown with dashed red lines. The model reproduces correctly the time scales of the force dynamics as well as the transient and asymptotic amplitudes of the forces. However, it is observed that the fine details of the transient dynamics, at the early stage of the linear instability, are not satisfactorily reproduced in the identification process (figure 12(a,b) at $t\approx 590$ and 475 respectively). The ranges of time concerned, in both cases, are also associated with oscillations in $a_4$, as observed during the initial stage at $t\approx 590$ in figure 6(a) and $t\approx 475$ in figure 6(b). This strongly suggests that the oscillations of $a_4$ be triggered by the degrees of freedom associated with the Hopf bifurcation. This means that the degrees of freedom of the pitchfork bifurcation are affected by the degrees of freedom of the Hopf bifurcation, at least when the distance from the bifurcation point is large enough, which is the case at $Re=80$.

In addition, as recalled in § 3.3, at $Re\approx 68$, both the steady symmetric solution and the symmetric-based limit cycle undergo a supercritical pitchfork bifurcation. We emphasize that this coincidence of two local pitchfork bifurcations might not occur by chance, as mentioned in Deng et al. (Reference Deng, Noack, Morzyński and Pastur2020). As a result of these two simultaneous bifurcations, the degrees of freedom involved in the pitchfork bifurcation of the fixed point might not coincide with those involved in the pitchfork bifurcation of the limit cycle. For this reason, it is reasonable to introduce two distinct sets of degrees of freedom for each of them, namely $a_4$, $a_5$ at the fixed point and $a_6$, $a_7$ at the limit cycle. These two additional degrees of freedom will complete the mean-field model with more details and will take into account the mean-field distortion during the transition from the fixed point to the limit cycle. The new resulting mean-field Galerkin system, with seven degrees of freedom, is derived in appendix E, while the new resulting force model is discussed in the next subsection.

5.3. The need for additional modes

All our attempts to smooth out the kinks observed at the beginning of the exponential growth, in both $C_D$ and $C_L$ in the frame of the force model (5.2), failed, even when over-fitting the model without any sparsity. This strongly indicates that five degrees of freedom might not be sufficient to account for the force evolution on the full time range.

Digging into this idea, it becomes manifest that the way $\boldsymbol {u}_{3}(\boldsymbol {x})$ is built, namely as the difference between the statistically symmetric mean flow field, associated with the unstable limit cycle, and the symmetric steady solution $\boldsymbol {u}_s(\boldsymbol {x})$,

(5.3)\begin{equation} \boldsymbol{u}_3(\boldsymbol{x}) = \bar{\boldsymbol{u}}_{T}(\boldsymbol{x},775) - \boldsymbol{u}_s(\boldsymbol{x}), \end{equation}

see figure 12(a), does not allow us to satisfactorily account for the complete dynamics of the lift and drag forces. This also indicates that $\boldsymbol {u}_3(\boldsymbol {x})$ gets distorted when the system is evolving along the manifold, which connects the unstable limit cycle to one of the two conjugated stable limit cycles. In other words, the mean-field distortion on the attractors associated with the two asymmetric mean flow fields $\bar {\boldsymbol {u}}^{\pm }$, namely

(5.4)\begin{equation} \boldsymbol{u}_3^{{\pm}}(\boldsymbol{x}) = \bar{\boldsymbol{u}}^{{\pm}}(\boldsymbol{x}) - \boldsymbol{u}_s^{{\pm}}(\boldsymbol{x}) \end{equation}

do not coincide exactly with $\boldsymbol {u}_3(\boldsymbol {x})$. The asymmetric mean flow fields $\bar {\boldsymbol {u}}^{\pm }$ only focus on the post-transient dynamics, as shown in figure 12(c), which can be expressed with $\bar {\boldsymbol {u}}_{T}^{\pm }(\boldsymbol {x},700)$. The difference between $\boldsymbol {u}_3^{\pm }(\boldsymbol {x})$ and $\boldsymbol {u}_3(\boldsymbol {x})$ is asymmetric and can be decomposed into a symmetric and an antisymmetric part, respectively $\boldsymbol {u}_6(\boldsymbol {x})$ and $\boldsymbol {u}_7 (\boldsymbol {x})$:

(5.5)\begin{equation} \boldsymbol{u}_3^{{\pm}}(\boldsymbol{x})-\boldsymbol{u}_3(\boldsymbol{x})={\pm} \boldsymbol{u}_6(\boldsymbol{x}) + \boldsymbol{u}_7(\boldsymbol{x}). \end{equation}

As a result, the modes $\boldsymbol {u}_6(\boldsymbol {x})$ and $\boldsymbol {u}_7(\boldsymbol {x})$ can be defined as,

(5.6a)\begin{gather} \boldsymbol{u}_6 (\boldsymbol{x}) \propto \bar{\boldsymbol{u}}^{+} (\boldsymbol{x}) - \bar{\boldsymbol{u}}^{-}(\boldsymbol{x}) , \end{gather}

(5.6b)\begin{gather}\boldsymbol{u}_7 (\boldsymbol{x}) \propto (\bar{\boldsymbol{u}}^{+} (\boldsymbol{x}) + \bar{\boldsymbol{u}}^{-}(\boldsymbol{x})) - 2 \bar{\boldsymbol{u}}_T(\boldsymbol{x},775). \end{gather}

After orthogonal normalization by a Gram–Schmidt procedure, the resulting modes are shown in figure 13, with their mode amplitudes in figure 14. When comparing the definitions of $\boldsymbol {u}_6$ and $\boldsymbol {u}_7$ in (5.6) and of $\boldsymbol {u}_4$ and $\boldsymbol {u}_5$ in (3.2) and (3.3), it is not surprising that the spatial structure of $\boldsymbol {u}_6$, respectively $\boldsymbol {u}_7$ (see figure 13), be so similar to the spatial structure of $\boldsymbol {u}_4$, respectively $\boldsymbol {u}_5$ (see figure 5). We also mention that, $\boldsymbol {u}_6$, $\boldsymbol {u}_7$ as defined in (5.6), would be equivalent to the pitchfork modes $\boldsymbol {u}_4$, $\boldsymbol {u}_5$ built on the periodic solutions instead of being built on the steady solutions. However, after the Gram–Schmidt procedure, the $\boldsymbol {u}_6$, $\boldsymbol {u}_7$ modes of figure 13 have been transformed into corrective modes of $\boldsymbol {u}_4$, $\boldsymbol {u}_5$ when departing from the steady solutions and approaching the asymptotic limit cycles. The corrective modes $\boldsymbol {u}_6$, $\boldsymbol {u}_7$ should be slaved to $\boldsymbol {u}_4$, $\boldsymbol {u}_5$ along the mean-field distortion of $\boldsymbol {u}_{3}$. The corresponding slaving relation will not be discussed in this paper. Hence, the combination of $\boldsymbol {u}_i$, $i=4,\ldots ,7$, works as a flexible pitchfork mode expansion, which adapts the whole phase space where all the invariant sets (steady/periodic) locate.

Figure 13. Vortical structure (colour) of the modes $\boldsymbol {u}_6(\boldsymbol {x})$ (a), $\boldsymbol {u}_7(\boldsymbol {x})$ (b), at $Re=80$.

Figure 14. Mode amplitudes $a_{6,7}(t)$ in the full-flow dynamics starting (a) from the symmetric steady solution $\boldsymbol {u}_s$, (b) from the asymmetric steady solution $\boldsymbol {u}_s^{+}$, at $Re=80$.

In figure 14, the transient dynamics of $a_6$, $a_7$ is shown to be also similar to $a_4$, $a_5$ in figure 6. Not surprisingly, the opposite initial bump of $a_6$, $a_7$ helps to better fit the dynamics on the manifold. Besides, $a_6$, $a_7$ show no contribution close to the steady solutions, as their role is to adapt the modes $\boldsymbol {u}_4$, $\boldsymbol {u}_5$ when approaching the stable limit cycle.

The force model identification is more challenging with these two additional modes. High robustness is required for our force model without losing the identified terms in § 5.1. Compared to the force formula (5.1) with five modes, 8 new terms are introduced in the drag formula, namely $a_7,a_1a_6,a_2a_6,a_4a_6,a_6^{2},a_3a_7,a_5a_7,a_7^{2}$, and 7 additional terms are considered in the lift formula, namely $a_6,a_3a_6,a_5a_6,a_1a_7,a_2a_7,a_4a_7,a_6a_7$. Due to the similar transient dynamics of $a_4$, $a_5$ and $a_6$, $a_7$, the corrective degrees of freedom $a_6$, $a_7$ can easily replace $a_4$, $a_5$ in the identified model. Hence, the original structure of the force model with five modes could be lost. To avoid possible over-fitting, we need to free the active terms gradually and constrain the parameters of $a_4$, $a_5$ during the sparse regression to ensure the robustness of the result. In addition, the newly introduced terms should work as a corrective function to the original force model with five degrees of freedom. In other words, the new force model with seven degrees of freedom should inherit the original structure of (5.2).

Based on the structure of the drag model (5.2a), the terms $a_7$, $a_7a_7$, $a_3a_7$, $a_4a_6$ and $a_5a_7$ are introduced in the extended model. The terms $a_1a_6,a_2a_6,a_6^{2}$ are firstly set to zero because their corresponding terms $a_1a_4,a_2a_4,a_4^{2}$ in (5.2a) are vanishing. In order to improve the robustness of the regression results, the terms $a_5$ and $a_5^{2}$ are constrained with the values from (5.2a). Increasing the $L1$-penalty of the LASSO regression, $l_{x; 7}$, $q_{x; 4 6}$ and $q_{x; 7 7}$ vanish successively, and an obvious under-fitting starts when losing $q_{x; 3 5}$. The introduced terms $q_{x; 3 7}$, $q_{x; 5 7}$ are robust with little possibility of over-fitting. Eventually, the drag model reads

(5.7)\begin{align} C_D &= 3.77331204 + 0.05888312 a_5 - 0.00169970 a_1^{2} - 0.00156775 a_2^{2} \nonumber\\ &\quad + 0.00513885 a_3^{2} + 0.00786294 a_3a_5 + 0.00950204 a_3a_7 \nonumber\\ &\quad - 0.25910470 a_5^{2} - 0.06264888 a_5a_7 . \end{align}

Equation (5.7) preserves the original form of (5.2a), with tiny changes to the coefficients. This extended model fits well the dynamics of the drag coefficient, with the $r^{2}$ score increasing to $0.9981$, which also can be seen with the red dashed curve of figure 15(a,b).

Figure 15. Performance of the force model with two additional slaved corrective modes. Time evolution of the drag $C_D$ (a,b) and lift $C_L$ (c,d) coefficients in the full-flow dynamics (solid black line) and for the force model (red dashed line), at $Re=80$. Initial condition: (a,c) symmetric steady solution $\boldsymbol {u}_s$, (b,d) asymmetric steady solution $\boldsymbol {u}_s^{+}$.

As already mentioned, the drag monotonically increases with the development of the vortex shedding. This is obvious, for instance, from figure 15, when the lift starts to oscillate and the drag to increase. The positive signs of $q_{y; 3 3}$, $q_{y; 3 5}$ and $q_{y; 3 7}$, in the drag model of (5.7), are responsible for this monotonic increase of the drag. Compared to the drag model with only $a_5$ in § 5.1, the contribution to the drag of $a_5$ and $a_7$ is more subtle. They contribute to an increase of the drag through $a_5$, $a_5a_3$ and $a_7a_3$, while they promote a decrease of the drag through $a_5^{2}$ and $a_5a_7$. As a non-trivial result, the statistically asymmetric (stable) limit cycles have a larger drag than the statistically symmetric (unstable) limit cycle, while the asymmetric steady solutions have a lower drag than the symmetric steady solution. This is obvious in figure 11 when considering the relative positions of the three steady solutions and three limit cycles along the $C_D$ axis. Note that the parameters $q_{x; 1 1}$, $q_{x; 2 2}$ and $q_{x; 5 7}$ all own negative signs but are relatively small. The two parameters $q_{x; 1 1}$, $q_{x; 2 2}$ solely contribute to the oscillating dynamics, as discussed in § 4.1, while $q_{x; 5 7}$ optimizes the fitting result when evolving toward the attracting limit cycles.

Analogously, for the lift model, the values of $l_{y; 4}$ and $q_{y; 4 5}$ are taken from the identified lift model in (5.2b), while $q_{y; 1 7}$, $q_{y; 2 7}$ are set to zero for consistency with the structure of (5.2b), in which $q_{y; 1 5}$, $q_{y; 2 5}$ are absent. Based on the structure of model (5.2b), the terms $a_6$, $a_6a_7$, $a_3a_6$, $a_5a_6$ and $a_4a_7$ are introduced in the extended model. The final sparse form is identified by the LASSO regression on gradually increasing the $L1$-penalty. A sparse lift model, compatible with the structure of (5.2b), is derived as

(5.8)\begin{align} C_L &= 0.00762433 a_1 + 0.01102097 a_2 - 0.10179203 a_4 - 0.03129798 a_6 \nonumber\\ &\quad - 0.00141416 a_1a_3 - 0.00289952 a_2 a_3 + 0.00656293 a_3a_4 - 0.01082375 a_3a_6\nonumber\\ &\quad + 0.05914386 a_4a_5 + 0.02365784 a_4a_7 - 0.03348935 a_5a_6. \end{align}

In addition to the lift model of (5.2b), the lift model of (5.8) contains the terms $a_{6}$, $a_3a_6$ and $a_5a_6$, as well as the coupling between $a_4$ to $a_7$. The $r^{2}$ score has increased to $0.9952$. Both the oscillating dynamics in the early stage and the symmetry-breaking stage are better reproduced for the lift coefficient, as the red dashed curve of figure 15(c,d) proves.

With the two additional degrees of freedom $a_6$, $a_7$, the time evolution of the drag and lift coefficients are well reproduced, as shown in figure 15. Without notable changes of the original lift structure, the phase of the lift dynamics is now correctly caught along with the complete transient dynamics.

5.4. Force model at $Re=100$

In § 4.2, we derived a basic force formula for the primary stage of the transient evolution at $Re=100$, when only the degrees of freedom of the pitchfork bifurcation were involved. We now consider the complete force evolution at $Re=100$. Figure 16 shows trajectories issued from the three different steady solutions in the three-dimensional time-delayed embedding space of $C_L$ and $C_D$. The black trajectory, issued from the symmetric steady solution $\boldsymbol {u}_s$ (black cross $\times$ in figure 16) first approaches the asymmetric steady solution $\boldsymbol {u}_s^{+}$ (red point) before escaping out of it and eventually reaching the stable (statistically asymmetric) limit cycles around $\bar {\boldsymbol {u}}^{+}$.

Figure 16. Trajectories in the time-delayed embedding space of the lift $C_L$ and drag $C_D$ coefficients, with $\tau =2$, at $Re=100$. Black trajectories starting close to the symmetric steady solution $\boldsymbol {u}_s$ ($\times$); red trajectory starting close to the asymmetric steady solution $\boldsymbol {u}_s^{+}$ ($\bullet$); blue trajectory starting close to the asymmetric steady solution $\boldsymbol {u}_s^{-}$ ($\blacksquare$).

The same mode decomposition strategy is proposed, resulting in a reduced-order model with 7 modes. The mode amplitudes from two DNS, starting from either the symmetric steady solution $\boldsymbol {u}_s$ (a) or the asymmetric steady solution $\boldsymbol {u}_s^{+}$ (b), are shown in figure 17.

Figure 17. Mode amplitudes $a_{1,\dots ,7}(t)$ in the full-flow dynamics starting (a) from the symmetric steady solution $\boldsymbol {u}_s$, (b) from the asymmetric steady solution $\boldsymbol {u}_s^{+}$, at $Re=100$.

As already observed in figure 10, the drag coefficient (solid black line) in figure 18(a) exhibits a minimal value for a transient state around $t\approx 700$. This transient state is the asymmetric steady solution $\boldsymbol {u}_s^{+}$ (red circle of figure 16). In the frame of our modal decomposition (3.1), $\boldsymbol {u}_s^{+}$ is approximated as

(5.9)\begin{equation} \boldsymbol{u}_s^{+} \approx \boldsymbol{u}_s + a_4(700)\boldsymbol{u}_4 + a_5(700)\boldsymbol{u}_5, \end{equation}

with only $a_4$ and $a_5$ being active in the dynamics of the fluidic pinball, as can be seen in figure 17(a). From (4.8a), the drag coefficient only depends on $a_5$ and $a_5^{2}$, which actually contribute to the transitory increase and an overall decrease on the drag. This is fully consistent with the transition of the drag coefficient observed in figures 10(a) and 18(a) from $t=300$ to $700$; $a_5$ is found to contribute to the initial rising of $C_D$, around $t\approx 650$, while $a_5^{2}$ contributes to the subsequent decrease of the drag coefficient, around $t=700$. The degrees of freedom associated with the Hopf bifurcation become active later during the transient dynamics, when the state space orbit leaves the unstable asymmetric steady solution $\boldsymbol {u}_s^{+}$ toward the stable attracting limit cycle around $\bar {\boldsymbol {u}}_s^{+}$.

Figure 18. Performance of the force model with two additional slaved corrective modes. Time evolution of the drag $C_D$ (a,b) and the lift $C_L$ (c,d) coefficients in the full-flow dynamics (solid black line) and for the force model (red dashed line), at $Re=100$. Initial condition: (a,c) symmetric steady solution $\boldsymbol {u}_s$, (b,d) asymmetric steady solution $\boldsymbol {u}_s^{+}$.

The training data are the real force coefficients and the mode amplitudes taken from the DNS starting with the three different steady solutions to the final asymptotic regimes. Following the same calibration procedure as for $Re=80$, we first apply the LASSO regression for the force model with the five leading degrees of freedom, and then introduce the two additional degrees of freedom $a_6$, $a_7$ into the regression for optimization. Performing the sparse regression in this way can prevent the elimination of $a_4$, $a_5$ and ensure the corrective effect of $a_6$, $a_7$, thereby improving the robustness of the identification. The force model at $Re=100$ reads

(5.10a)\begin{align} C_D &= 3.58248992 + 0.04367604 a_5 - 0.00302817 a_1^{2} - 0.00354079 a_2^{2} \nonumber\\ &\quad + 0.00158873 a_3^{2} + 0.02169661 a_3a_5 + 0.02223079 a_3a_7 \nonumber\\ &\quad - 0.08525184 a_5^{2} - 0.04763643 a_5a_7 , \end{align}

(5.10b)\begin{align} C_L &= 0.00346208 a_1 + 0.00269236 a_2 - 0.13611053 a_4 + 0.05962648 a_6\nonumber\\ &\quad + 0.00029274 a_1a_3 - 0.00045784 a_2a_3 + 0.00389912 a_3a_4 - 0.02102284 a_3a_6 \nonumber\\ &\quad + 0.09194312 a_4a_5 + 0.02056288 a_4a_7 - 0.10980990 a_5a_6, \end{align}

with $r^{2} = 0.9984$ for the drag model of (5.10a), and $r^{2} = 0.9901$ for the lift model of (5.10b). As shown in figure 18, the force model fits well the time evolution of the drag and lift coefficients. Moreover, (5.7), (5.8) and (5.10) own the same active terms. Henceforth, the drag force model preserves the same structure with the same signs of the active terms as the Reynolds number is increased. In addition, although the transient dynamics at $Re=80$ and 100 is qualitatively very different, with the seven degrees of freedom differently activated during the transient, the force model of (5.7) and (5.8) is still consistent at $Re=100$, with the correctly identified mean-field model. For the lift model (5.10b), we notice the same structure with the sign changes for the terms $a_1a_3$ and $a_6$, compared to (5.8), which is acceptable for the oscillating dynamics. Compatible with the basic lift force model, the lift force model with seven degrees of freedom also correctly identifies the force transitions, as shown in figure 18(c,d).