3.1. Spectral proper orthogonal decomposition (SPOD) of the flow

To obtain more insight, SPOD (Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018) is applied to the fluctuating flow field. The method will identify the structures beating at the dominant frequencies identified earlier. A total of $K = 100\,000$ snapshots of velocity perturbations obtained every $\delta t=0.01$ are stacked column by column to form matrix $Y$ ,

(3.1) \begin{equation} Y = \left [\begin{array}{cccc} {u^{\prime }_1}^{(1)} & {u^{\prime }_1}^{(2)} & {\cdots} & {u^{\prime }_1}^{(K)} \\[3pt] {v^{\prime }_1}^{(1)} & {v^{\prime }_1}^{(2)} & {\cdots} & {v^{\prime }_1}^{(K)} \\[3pt] \vdots & \vdots & {\cdots} & \vdots \\[3pt] {u^{\prime }_N}^{(1)} & {u^{\prime }_N}^{(2)} & {\cdots} & {u^{\prime }_N}^{(K)} \\[3pt] {v^{\prime }_N}^{(1)} & {v^{\prime }_N}^{(2)} & {\cdots} & {v^{\prime }_N}^{(K)} \end{array}\right ]\!. \end{equation}

The superscripts denote the snapshot index $k = 1, 2, \ldots , K$ and the subscripts the cell number $n = 1, 2, \ldots , N$ , where $N=84\,253$ (cell count of mesh 2). Matrix $Y$ is divided into 11 blocks with a 50 % overlap. For each frequency $\textit{St}_k$ , a cross-spectral density (CSD) matrix $ \boldsymbol{S}(St_k)$ is formed and the SPOD eigenvalue problem is solved to find the SPOD modes $ \varPhi ^{st_k}$ and their energies $ \lambda (St_k)$ . Since the CSD matrix is Hermitian, the eigenvalues are real but the eigenvectors (SPOD modes) are complex, i.e. $ \varPhi ^{St_k} = \varPhi _r^{St_k} + i \varPhi _i^{St_k}$ ; for more details, see Towne et al. (Reference Towne, Schmidt and Colonius2018).

Figure 8.

Variation of SPOD mode energies (eigenvalues) across frequency.

Figure 8 presents the variation of SPOD eigenvalues across frequency. Mode 1 refers to the leading SPOD mode at each frequency. Each curve represents the energy of the corresponding SPOD mode at different frequencies. Mode 1 exhibits two clear peaks at $0.063$ and $0.160$ that match closely $\textit{St}_1$ and $\textit{St}_2$ , respectively, of the $C_{\!D}$ spectrum. For the lower frequency band, the energy of mode 1 is markedly higher than that of mode 2, indicating that the gap-jet flapping motion is well approximated by just one mode. In the higher frequency band, the energy separation between modes 1 and 2 is not so pronounced. Note also the presence of a band with $0.080$ , the subharmonic of the $0.160$ .

Figure 9.

Contour plots of streamwise (left column) and cross-stream (right column) velocities of the real part of SPOD mode 1 at (a) $ St = 0.063$ , $\varPhi _{1r}^{0.063}$ , (b) $\textit{St}=0.160$ , $\varPhi _{1r}^{0.160}$ and (c) $ St = 0.080$ , $\varPhi _{1r}^{0.080}$ .

Figure 9 shows contours of the streamwise ( $u$ ) and transverse ( $v$ ) velocities of the real part of SPOD mode 1 at $\textit{St}=0.063$ , $0.160$ and the subharmonic $0.080$ . At $\textit{St} = 0.063$ , the $v$ contours peak at the gap centreline in the near wake. Both $u$ and $v$ reach their maxima at $x \approx 4$ , precisely where the jet-flapping motion is most discernible. At $\textit{St} = 0.160$ (the vortex-shedding frequency), a train of small-scale structures appears in the wake of each prism in the $v$ contours. The peaks now appear behind the prisms, not in the gap centreline. Notice that the peaks are displaced further away from the centreline as $x$ increases, indicating the wake expansion. The $v$ field is antisymmetric with respect to the gap centreline, but the $u$ field is symmetric. The $v$ mode attains maximum values in the near-wake region ( $x\lt 4$ ), as for $\textit{St}=0.063$ . However, the $v$ contours of the subharmonic at $\textit{St}=0.080=0.160/2$ ) grow rapidly in the near wake (where vortex merging occurs), peaks in the mid-wake $x \in [12,20]$ and then decays.

Overall, the dominant SPOD mode captures the three key phenomena identified earlier: jet-flapping, vortex shedding behind each cylinder and vortex pairing. As can be seen in the SPOD energy spectrum, gap-jet-flapping carries the largest energy, followed by vortex pairing and then asynchronous shedding.

Having characterised the main features of the flow, the next step is to consider the flow as a dynamical system, and analyse the Lyapunov exponents (LEs) and covariant Lyapunov vectors (CLVs). A key objective is to explore the footprint of the above-mentioned features onto the CLVs.

4.1. Lyapunov exponents

After discretisation, the system of governing equations (2.1) can be put in the general form,

(4.1) \begin{equation} \frac {{\rm d} \boldsymbol{u}}{{\rm d}t} = \boldsymbol{f}\!\left (\boldsymbol{u}\right ), \quad \boldsymbol{u}(0) = \boldsymbol{u}_0, \end{equation}

where $\boldsymbol{u}(t)$ is the state vector $\boldsymbol{u}(t) = \{u_1(t), v_1(t), \ldots , u_{N}(t), v_{N}(t)\} \in \mathbb{R}^{2N}$ , and $N$ is the number of cells. To get (4.1), we have used the continuity equation to eliminate pressure. This is done by constructing a Poisson equation for pressure (Pope Reference Pope2000), which can be notionally solved, and thus, $\boldsymbol{\nabla }\!P$ can be expressed in terms of velocities. As the dynamical system evolves, the state vector traces a trajectory in phase space. We define the evolution operator $\mathcal{L}^{t}$ that maps $\boldsymbol{u}(0)$ to $\boldsymbol{u}(t)$ under the dynamics (4.1), $\mathcal{L}^{t}(\boldsymbol{u}(0)) =\boldsymbol{u}(t)$ .

A key concept in dynamical systems is ergodicity, which implies that the time-average state of the system is independent of the initial condition (Birkhoff Reference Birkhoff1931). For ergodic systems, the long-time trajectory converges to a bounded region in phase space known as attractor (Ruelle Reference Ruelle1980). Attractors can be classified into three types: fixed (or equilibrium) points, limit cycles and strange attractors. Fixed-point attractors are stable points that trajectories converge to over time, limit-cycle attractors correspond to stable periodic orbits and strange attractors represent chaotic systems that exhibit complex, irregular patterns (as the one considered in this study).

We define our (assumed ergodic) dynamical system in the $2N$ -dimensional Riemannian manifold $\mathcal{M}$ . The Oseledets (Reference Oseledets1968) theorem establishes that a set of Lyapunov exponents (LEs) exists, provided there is a tangent space $\mathcal{T}_{\boldsymbol{u}} \mathcal{M}$ that describes the local geometry of the attractor for every state $\boldsymbol{u}(t)$ . Here, $\mathcal{T}_{\boldsymbol{u}} \mathcal{M}$ is a $2N$ -dimensional vector space consisting of all the possible directions in which the state $\boldsymbol{u}(t)$ can be perturbed. To study the response of a dynamical system to an arbitrary perturbation $\boldsymbol{\delta \boldsymbol{u}} \in \mathcal{T}_{\boldsymbol{u}} \mathcal{M}$ , we define the local expansion rate as

(4.2) \begin{equation} \gamma \!\left ( \boldsymbol{\delta \boldsymbol{u}}, \boldsymbol{u},t \right ) = \frac {\left \|D \mathcal{L}^t(\boldsymbol{u}) \boldsymbol{\delta \boldsymbol{u}} \right \|}{\| \boldsymbol{\delta \boldsymbol{u}}\|} ,\end{equation}

where $\|\boldsymbol{\cdot }\|$ denotes the Euclidean norm and $D \mathcal{L}^t(\boldsymbol{u})$ is the linear evolution operator of the tangent space. More specifically, $D \mathcal{L}^t(\boldsymbol{u})$ evolves the linearised form of (4.1) around $\boldsymbol{u}$ ,

(4.3) \begin{equation} \frac {\partial (\boldsymbol{\delta \boldsymbol{u}})}{\partial t}= \frac {\partial \! \boldsymbol{f} }{\partial \boldsymbol{u} } \delta \boldsymbol{u}, \end{equation}

over a time window $t$ . In other words, if the perturbation at the start of the time window is $\boldsymbol{\delta \boldsymbol{u}}$ , the notation $D \mathcal{L}^t(\boldsymbol{u}) \boldsymbol{\delta \boldsymbol{u}}$ is the result of the action of this operator, i.e. $D \mathcal{L}^t(\boldsymbol{u}) \boldsymbol{\delta \boldsymbol{u}}=\boldsymbol{\delta \boldsymbol{u}}(t)$ . Thus, the operator performs the mapping $\mathcal{T}_{\boldsymbol{u}} \mathcal{M} \rightarrow \mathcal{T}_{\mathcal{L}^t(\boldsymbol{u})} \mathcal{M}$ . For the Navier–Stokes equations, the linearised form is

(4.4) \begin{equation} \begin{aligned} \frac {\partial (\boldsymbol{\delta \boldsymbol{u}})}{\partial t}+ \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }(\boldsymbol{\delta \boldsymbol{u}}) +\boldsymbol{\delta \boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} & =-\boldsymbol{\nabla }(\delta p)+\frac {1}{\textit{Re}} \Delta (\boldsymbol{\delta \boldsymbol{u}}), \\ \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\delta \boldsymbol{u}}) & =0, \end{aligned} \end{equation}

where $\boldsymbol{\delta \boldsymbol{u}}= (\delta u, \delta v)$ .

The Lyapunov exponent $\lambda$ associated with the perturbation vector $\boldsymbol{\delta \boldsymbol{u}}$ is the long-time average,

(4.5) \begin{equation} \lambda (\boldsymbol{\delta \boldsymbol{u}}) = \lim _{t \rightarrow \infty } \frac {1}{t} \ln \frac {\left \|D \mathcal{L}^t(\boldsymbol{u}) \boldsymbol{\delta \boldsymbol{u}}\right \|}{\| \boldsymbol{\delta \boldsymbol{u}}\|} .\end{equation}

This exponent represents the average exponential rate of expansion (or contraction) of the perturbation $\boldsymbol{\delta \boldsymbol{u}}$ , i.e. $\left \| \boldsymbol{\delta \boldsymbol{u}}(t)\right \| \sim \| \boldsymbol{\delta \boldsymbol{u}}(0)\| \mathrm{e}^{\lambda t}$ . The Oseledets (Reference Oseledets1968) theorem establishes that $\lambda$ assumes a finite number of distinct values $\lambda _1 \gt \lambda _2 \gt {\cdots} \gt \lambda _m$ , where $m \leqslant 2N$ , as the perturbation vector $\boldsymbol{\delta \boldsymbol{u}}$ varies in the tangent space $\mathcal{T}_{\boldsymbol{u}} \mathcal{M}$ .

The set of LEs, known as the ‘Lyapunov spectrum’, is a characteristic of the dynamical system. In a chaotic dynamical system, there is at least one positive LE, i.e. $\lambda _1 \gt 0$ . This means that the perturbed trajectory $ \boldsymbol{\boldsymbol{u}}(t)+\boldsymbol{\delta \boldsymbol{u}}(t)$ will diverge exponentially from the reference trajectory $ \boldsymbol{\boldsymbol{u}}(t)$ at an average rate $\lambda _1$ .

When a randomly chosen set of $m \leqslant 2N$ infinitesimal perturbations $(\boldsymbol{\delta \boldsymbol{u}}_1, \boldsymbol{\delta \boldsymbol{u}}_2,{} \ldots , \boldsymbol{\delta \boldsymbol{u}}_m)$ are propagated in the tangent space using (4.3) (or (4.4) in our case), the angles between $\boldsymbol{\delta \boldsymbol{u}}_i$ reduce with time as all perturbations align along the most unstable direction. Thus, one can compute only the leading LE: identifying smaller LEs becomes numerically ill-conditioned. Benettin et al. (Reference Benettin, Galgani, Giorgilli and Strelcyn1980a ,Reference Benettin, Galgani, Giorgilli and Strelcyn b ) devised an algorithm in which the set of $m$ perturbations is propagated for a time interval $t_{\textit{step}}$ ; the evolved ${\boldsymbol{\delta \boldsymbol{u}}}_i$ are then orthonormalised and replaced by a set of orthonormal vectors that span the same subspace. These vectors are used as initial conditions and are propagated for another time segment $t_{\textit{step}}$ . Benettin et al. (Reference Benettin, Galgani, Giorgilli and Strelcyn1980a ) proved that using Gram–Schmidt QR decomposition, a set of coordinate-independent LEs can be obtained. Coordinate independence means that the computed LEs are intrinsic to the system dynamics and therefore independent of the chosen coordinate system, and thus the norm used to compute them, see also Bewley (Reference Bewley1999). The process is presented in algorithm 1.

Algorithm 1:

Calculation of Lyapunov exponents (Benettin et al. Reference Benettin, Galgani, Giorgilli and Strelcyn1980a )

Figure 10.

Convergence of LEs and forward-backward evolution to calculate CLVs for $t_{\textit{step}} = 0.3$ .

To compute the LEs of the flow under consideration, we evolve six orthonormal perturbations $\{\boldsymbol{\delta \boldsymbol{u}} \}_{i=1 \ldots 6}$ in the tangent space. The general perturbations $\boldsymbol{\delta \boldsymbol{u}}$ can be expressed as the Fourier integral in the spanwise direction,

(4.7) \begin{equation} \boldsymbol{\delta \boldsymbol{u}}(x,y,z,t)=\int _{-\infty }^{+\infty } \boldsymbol{\widehat {\delta {\boldsymbol{u}}}}(x,y,t; \beta ) e^{\iota \beta z}\, {\rm d}\beta ,\end{equation}

and similarly for $\delta p$ . In the present work, we take $\beta =0$ , i.e. we consider only 2-D perturbations. The perturbations are orthonormalised periodically using Gram–Schmidt QR decomposition. Two segment lengths are chosen $t_{\textit{step}} = 0.3, 10$ . The convergence of the LEs with the integration time for $t_{\textit{step}} = 0.3$ is shown in figure 10. The unstable manifold consists of two positive LEs, which indicates that, on average, perturbations grow exponentially along two directions at any point on the attractor. The leading LE is $\lambda _1 \approx 0.1$ . One LE is equal to 0 (neutral manifold) which indicates the absence of growth or decay of perturbations along the trajectory in phase space. Furthermore, the system is non-degenerate as all LEs are distinct (Ginelli et al. Reference Ginelli, Chaté, Livi and Politi2013).

Figure 11.

Lyapunov spectrum of the flow.

The LEs were also computed for $t_{\textit{step}} = 10$ , which is approximately one Lyapunov time unit, $ {1}/{\lambda _1}$ . The LE spectra for the two values of $t_{\textit{step}}$ are shown in figure 11. The results are almost identical, confirming the robustness of the calculation.

4.2. Covariant Lyapunov vectors

It was shown that LEs converge to constant values when averaged over a long time window. The corresponding directions of growth/decay are known as covariant Lyapunov vectors (CLVs), see Ginelli et al. (Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007, Reference Ginelli, Chaté, Livi and Politi2013). CLVs provide a geometric characterisation of the tangent space. They are unit-norm vectors that span the Gram–Schmidt (GS) subspace and represent the expanding, contracting and neutral directions at each point of the attractor. These directions correspond to the positive, negative and zero LEs, respectively.

The concept of CLVs was first introduced by Oseledets (Reference Oseledets1968) and later formalised as tangent directions of invariant manifolds by Ruelle (Reference Ruelle1979). However, for decades, there were no effective algorithms for their computation. Ginelli et al. (Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007) proposed a stable algorithm to extract CLVs, leveraging the geometric understanding of the tangent space of Ruelle (Reference Ruelle1979). The algorithm is based on a forward-backward evolution of the tangent space. A brief overview is presented in Algorithm 2; the detailed methodology of Ginelli et al. (Reference Ginelli, Chaté, Livi and Politi2013) is presented in Appendix A.

Algorithm 2:

Calculation of covariant Lyapunov vectors (Ginelli et al. Reference Ginelli, Chaté, Livi and Politi2013)

The four phases of the algorithm are shown in figure 10. The algorithm begins with the forward transient phase (in our case, $t \in [0, 2100]$ ), in which the dynamical system and six (initially random) orthonormal perturbations are evolved until the tangent subspace and the LE spectrum converge. In the forward dynamics phase ( $t \in [2100, 4500]$ ), the orthonormal Gram–Schmidt vectors $\boldsymbol{G}_{\!n} = (\boldsymbol{g}_n^{(1)}\mid \boldsymbol{g}_n^{(2)}\mid \ldots \mid \boldsymbol{g}_n^{(6)} )$ and the $R_n$ matrices (obtained from the QR decomposition) are stored for every time segment, $n$ . In the backward transient phase ( $t \in [5400, 4500]$ ), the CLV expansion coefficient matrix $\boldsymbol{C}_{\!n}$ is evolved backwards. To determine the length of this phase, $k_2 t_{\textit{step}}$ was increased and the angles between the CLVs in the backward dynamics phase were monitored. At $k_2 t_{\textit{step}} = 150$ , the angles between the CLVs stabilised, indicating convergence. In the final backward dynamics phase, the recorded GS matrices $\boldsymbol{G}_{\!n}$ are multiplied with the expansion coefficient matrices $\boldsymbol{C}_{\!n}$ over the middle time horizon $t \in [2100, 4500]$ to generate 8000 samples of six CLVs (see (A8) in Appendix A).

Each CLV $\boldsymbol{V}_{\!i}$ has dimension $2N$ and can be expressed as $\boldsymbol{V}_{\!i} = [V_i^u, V_i^v ]$ , where superscripts $u$ and $v$ denote the $x$ and $y$ velocity components, respectively. The magnitude of the $i$ th CLV, $ |\boldsymbol{V}_{\!i}|$ , is defined as the Euclidean norm, i.e. $|\boldsymbol{V}_{\!i}|=\sqrt {(V_i^u)^2 + (V_i^v)^2}$ .

Figure 12.

Contour plots of the magnitude of the leading CLV $|\boldsymbol{V}_1|$ . (a) Instantaneous $|\boldsymbol{V}_1|$ and (b) time-averaged $|\boldsymbol{V}_1|$ . The blue isolines mark the recirculation regions behind the cylinders and the red dots indicate the locations of peak magnitude. For an animation of $|\boldsymbol{V}_1|$ contours, see supplementary movie 2.

Figure 12 shows contours of the instantaneous and time-averaged magnitude of the leading unstable CLV, $|\boldsymbol{V}_1|$ (note that the two plots have different streamwise and cross-stream extents). The instantaneous contours (panel a) reveal that $|\boldsymbol{V}_1|$ takes the form of coherent structures close to the cylinders, but the values are decaying for $x\gt 8$ . Occasional high values are detected much further downstream and away from the centreline (at the time instant shown, they are in $x\in [28,34]$ ). These high values are concentrated around vortices that have formed close to the cylinders and are ejected away from the centreline by the sweeping motion of the flapping jet. They are then captured by the free stream and transported for long distances. The time-averaged contours are symmetric (as expected) and peak above and below the gap centreline at $x =2.5$ , slightly downstream of the recirculation zone. The large footprint extends up to $x \approx 8$ and then the values decay. The high values that were detected far downstream in the instantaneous plot are occasional and do not affect the time-average.

Figure 13.

Contour plot of the second unstable CLV: (a) instantaneous $|\boldsymbol{V}_2|$ and (b) mean $|\boldsymbol{V}_2|$ . The blue isolines mark the recirculation regions behind the cylinders and the red dots indicate the locations of peak magnitude. For an animation of $|\boldsymbol{V}_2|$ contours, see supplementary movie 3.

Figure 13 shows contours of the magnitude of the second unstable CLV, $|\boldsymbol{V}_2|$ . The instantaneous magnitude (visualised at the same time instant as that of $\boldsymbol{V}_1$ ) appears slightly weaker in the near-wake region ( $x \leqslant 4$ ) achieving higher intensity further downstream. Again, a region of high values appears further downstream and away from the centreline. The time-average peaks at $x = 2.9$ , but maintains high values until approximately $x = 14{-}16$ . This observation suggests that $\boldsymbol{V}_2$ has presence further downstream compared with $\boldsymbol{V}_1$ . It is clear that there is a link between the structures of the unstable CLVs and the flow characteristics analysed earlier. This connection will be elucidated in detail in § 5.

Figure 14.

Contour plot of the neutral CLV: (a) instantaneous $|\boldsymbol{V}_3|$ and (b) mean $|\boldsymbol{V}_3|$ .

Figure 14 presents contours of the neutral CLV, $\boldsymbol{V}_3$ . The instantaneous contours of $\boldsymbol{V}_3$ (visualised at the same time instant as those of $\boldsymbol{V}_1$ and $\boldsymbol{V}_2$ ) exhibit spatially distributed activity across the wake, but its footprint is slightly weaker close to the cylinders compared with $\boldsymbol{V}_1$ and $\boldsymbol{V}_2$ . The neutral $\boldsymbol{V}_3$ is aligned with the gradient of the velocity field, ${\partial \boldsymbol{u}}/{\partial t}$ , i.e. corresponds to perturbations along the solution trajectory in phase space. Perturbations in this direction neither grow nor decay, thereby the corresponding LE is 0. We have checked snapshots of the magnitudes of $\boldsymbol{V}_3$ and $ {\partial \boldsymbol{u}}/{\partial t}$ (approximated as $|( {\boldsymbol{u}^{n+1}-\boldsymbol{u}^n})/{\Delta t}|$ and normalised to unit norm) at several time instants, and confirmed that the plots are almost identical (results not shown for brevity). Contours of the time-averaged $\boldsymbol{V}_3$ magnitude (shown in panel b) display an almost uniform intensity along the streamwise direction ( $x$ ).

Figure 15.

Contour plot of the stable CLV: (a) instantaneous $|\boldsymbol{V}_6|$ and (b) mean $|\boldsymbol{V}_6|$ .

Figure 15 depicts contours of the magnitude of the stable CLV, $|\boldsymbol{V}_6|$ . The instantaneous contour (again visualised at the same time instant as that of $\boldsymbol{V}_1$ and $\boldsymbol{V}_2$ ) reveals negligible footprint until approximately $x \approx 25$ , with the intensity peaking further downstream, at $x \approx 35$ . The time-averaged contour shows that the strongest activity of this CLV is concentrated near the outlet boundary, where dissipation overwhelmingly dominates production, consistent with the negative Lyapunov exponent associated with $\boldsymbol{V}_6$ . These trends align with the findings of Ni (Reference Ni2019) who investigated the flow around a cylinder and reported a systematic downstream shift of the spatial CLV footprints as the LE values decrease.

4.3. Assessment of hyperbolicity

In uniformly hyperbolic systems, the tangent space can be decomposed distinctly into unstable ( $\boldsymbol{E}_{\boldsymbol{u}}^u$ ), stable ( $\boldsymbol{E}_{\boldsymbol{u}}^s$ ) and neutral ( $\boldsymbol{E}_{\boldsymbol{u}}^0$ ) subspaces,

(4.8) \begin{equation} \mathcal{T}_{\boldsymbol{u}} \mathcal{M}=\boldsymbol{E}_{\boldsymbol{u}}^u \oplus \boldsymbol{E}_{\boldsymbol{u}}^s \oplus \boldsymbol{E}_{\boldsymbol{u}}^0, \end{equation}

at any point on the attractor. These subspaces are spanned by the CLVs corresponding to positive, negative and zero LEs, respectively. For a system to be classified as uniformly hyperbolic, the three subspaces must be bounded away from each other, i.e. they should not align at any point on the attractor. For the flow examined, $\boldsymbol{V}_1$ and $\boldsymbol{V}_2$ span $\boldsymbol{E}_{\boldsymbol{u}}^u$ , $\boldsymbol{V}_3$ spans $\boldsymbol{E}_{\boldsymbol{u}}^0$ , and $\boldsymbol{V}_4$ , $\boldsymbol{V}_5$ and $\boldsymbol{V}_6$ span $\boldsymbol{E}_{\boldsymbol{u}}^s$ . Note that we have determined only three stable CLVs; therefore, our understanding of the stable subspace is limited.

The angle $ \phi ^{i,j}$ between $ \boldsymbol{V}_{\!i}$ and $ \boldsymbol{V}_{\!j}$ , is calculated as

(4.9) \begin{equation} \phi ^{i,j} = \cos ^{-1} \left ( |\boldsymbol{V}_{\!i} \boldsymbol{\cdot }\boldsymbol{V}_{\!j}| \right )\!, \end{equation}

where $ \boldsymbol{V}_{\!i} \boldsymbol{\cdot }\boldsymbol{V}_{\!j}$ is the dot product (defined as $\boldsymbol{V}_{\!i}^\top \boldsymbol{V}_{\!j}$ , where $()^\top$ denotes the transpose of a vector) and $ |\boldsymbol{\cdot }|$ is the absolute value. Equation (4.9) provides an angle between $0^{\circ }$ and $90^{\circ }$ (Ginelli et al. Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007; Xu & Paul Reference Xu and Paul2016). If $\phi ^{i,j}=0$ , the two CLVs are parallel to each other, while if $\phi ^{i,j}=90$ , they are orthogonal.

Figure 16.

(a) Variation of the angle between $\boldsymbol{V}_3$ and $\boldsymbol{V}_4$ over time, (b) contour plots of instantaneous $|\boldsymbol{V}_3|$ and $|\boldsymbol{V}_4|$ at time $ = 3000$ .

The angles $\phi ^{i,j}$ are computed from 8000 samples taken every $t_{\textit{step}} = 0.3$ time units in the time window $[2100, 4500]$ , see figure 10. Figure 16 illustrates the angle between $\boldsymbol{V}_3$ and $\boldsymbol{V}_4$ in this time window. Notice that $\phi ^{3,4}$ approaches zero at $t = 3000$ . Contours of the magnitude of the two CLVs at this time instant are shown in panel (b). The plots are almost identical confirming that the CLVs are parallel, which indicates violation of hyperbolicity at this time instant.

Figure 17.

Probability density functions (p.d.f.s) of angles between CLVs spanning (a) the unstable and neutral subspaces, (b) the neutral and stable subspaces, and (c) the unstable and stable subspaces. The vertical dashed lines indicate the angle $\phi =5^{\circ }$ .

Figure 17 shows the probability density functions (p.d.f.s) of the angles between pairs of CLVs that belong to different subspaces. The domain $[0^{\circ },90^{\circ }]$ is segmented into 18 bins (bin width $5^{\circ }$ ). We consider that samples located in the first bin $[0^{\circ }, 5^{\circ }]$ signify instances of hyperbolicity violation. Notice that CLVs that belong to different subspaces, but nearby indices $i$ and $j$ , have a higher probability of being parallel (or nearly-parallel) to each other. For example the p.d.f.s of $\phi ^{2,3}$ , $\phi ^{3,4}$ and $\phi ^{3,5}$ approach (or have samples) within the first bin.

However, when the indices are far apart, the CLVs have higher probability of orthogonality, see for example the p.d.f.s of $\phi ^{1,3}$ , $\phi ^{3,6}$ , $\phi ^{1,5}$ and especially $\phi ^{1,6}$ which has a sharp spike at $90^{\circ }$ . The orthogonality can be explained by the non-overlapping footprints of the CLVs. For example, $\boldsymbol{V}_1$ and $\boldsymbol{V}_6$ peak in the near wake and far-wake regions, respectively, see figures 12 and 15, thus, the overlap is minimal and the CLVs are almost orthogonal.

It was found that hyperbolicity is violated in $2.4\,\%$ of the $\phi ^{3,4}$ samples, $0.67\,\%$ of the $\phi ^{3,5}$ samples and $0.69\,\%$ of the $\phi ^{3,6}$ samples. These angle combinations correspond to CLVs that belong to neutral and stable subspaces, and indicate that hyperbolicity is mainly violated due to the alignment of these two subspaces. These observations are consistent with the findings of Xu & Paul (Reference Xu and Paul2016) and Ni (Reference Ni2019) who studied hyperbolicity for the flow around a circular cylinder and for Rayleigh–Bénard convection, respectively.

5.1. Analysis of the most unstable CLV, $\boldsymbol{V_1}$

Figure 18.

Variation of SPOD mode energies (eigenvalues) of $\boldsymbol{V}_1$ across frequency.

The energy spectrum of the first six SPOD modes of $ \boldsymbol{V}_1$ is shown in figure 18. Two distinct peaks are observed that lie within the lower and upper frequency bands. The energy separation between the first and second SPOD modes is strongest in the lower frequency band, as in the flow field modes.

Figure 19.

Contour plots of (a) streamwise velocity and (b) cross-stream velocity of the real part of SPOD mode 1 of $\boldsymbol{V}_1$ at $ St = 0.052$ (top row) and $ St = 0.156$ (bottom row) i.e. $\varPhi _{1r}^{0.052}$ and $\varPhi _{1r}^{0.156}$ .

Figure 20.

Contour plots of (a) streamwise velocity and (b) cross-stream velocity of the imaginary part of SPOD mode 1 of $\boldsymbol{V}_1$ at $ St = 0.052$ and $ St = 0.156$ i.e. $\varPhi _{1i}^{0.052}$ and $\varPhi _{1i}^{0.156}$ .

The real and imaginary parts of the dominant SPOD mode at these peaks are shown in figures 19 and 20, respectively. In each figure, contour plots of streamwise velocity ( $u$ ) and cross-stream velocity ( $v$ ) are shown in panels (a) and (b), respectively. Note that the imaginary part is spatially shifted with respect to the real part, consistent with a propagating structure beating at the frequency shown. At the lower peak $\textit{St} = 0.052$ , which is close to the flapping frequency of the flow (see figure 9), the SPOD mode exhibits structures characteristic of the flapping jet. The $v$ contours of both the real and imaginary parts show alternating positive and negative peaks along the gap centreline that decay downstream. The alternating signs in the near wake are consistent with the jet-flapping up and down motions. These features demonstrate that the low-frequency band of $ \boldsymbol{V}_1$ captures coherent structures in the near wake that peak in the gap centreline. In the $v$ contour of the higher band at $\textit{St} = 0.156$ , a train of small-scale structures is observed to shed from the cylinders and propagate downstream while expanding in the cross-stream direction. These structures correspond to the vortex shedding behind the cylinders and also closely resemble those of the SPOD mode 1 of the flow in the same frequency band. These observations confirm that the dominant SPOD mode of $ |\boldsymbol{V}_1|$ accurately captures the main instability features of the flow in both the lower and upper frequency bands. However, it does not contain any signature of the subharmonic.

5.2. Analysis of the second unstable CLV, $\boldsymbol{V_2}$

Figure 21.

Variation of SPOD mode energies (eigenvalues) of $\boldsymbol{V}_2$ across frequency.

The energy spectrum of the first six SPOD modes of the second unstable CLV, $ \boldsymbol{V}_2$ , is shown in figure 21. While the energy of $\boldsymbol{V}_1$ is spread across two frequency bands, the dominant SPOD mode of $\boldsymbol{V}_2$ exhibits a single prominent peak at $\textit{St} = 0.078$ , followed by a rapid decay. At this frequency, the energy of SPOD mode 1 is approximately three times that of mode 2, underscoring its importance. This frequency lies in the subharmonic band and is precisely half of the primary vortex shedding frequency $\textit{St} = 0.156$ , indicating that $\boldsymbol{V}_2$ captures the subharmonic instabilities resulting from interactions between paired vortices.

Figure 22.

Contour plots of (a) streamwise velocity and (b) cross-stream velocity of the real part of SPOD mode 1 of $ \boldsymbol{V}_2$ at $ St = 0.078$ i.e. $\varPhi _{1r}^{0.078}$ .

Contour plots of the real part of the streamwise ( $ u$ ) and cross-stream ( $ v$ ) velocity components of SPOD mode 1 at $ St = 0.078$ are shown in figure 22. Large alternating red and blue bands appear in the $ v$ contours in the mid-to-far wake region. These SPOD structures reveal a pairing instability that develops downstream, distinct from the near-wake structures captured by $ \boldsymbol{V}_1$ . A direct comparison with the structures of $\boldsymbol{V}_1$ at $\textit{St}=0.156 \,(=2 \times 0.078)$ , further supports this interpretation. As the energy and spatial extent of $\boldsymbol{V}_1$ structures begin to diminish around $x \approx 10$ , the structures associated with $\boldsymbol{V}_2$ begin to grow in intensity. This spatial shift signals an energy transfer from $\boldsymbol{V}_1$ to $\boldsymbol{V}_2$ , corresponding to the nonlinear pairing of vortices and the emergence of the subharmonic. In contrast to $\boldsymbol{V}_1$ , which represents the instability of individual vortices in the near wake, $\boldsymbol{V}_2$ reflects a slower-growing instability that dominates the far wake and emerges from vortex pairing.

The two unstable CLVs together capture all the distinct but interconnected mechanisms, that is, the primary instabilities due to vortex shedding and jet-flapping as well as the secondary instability driven by vortex merging.