3.1.1. Flow decomposition with mean-field consideration

The starting point of the HiCROM is the triple decomposition of the flow field similar to Reynolds & Hussain (Reference Reynolds and Hussain1972)

(3.1)\begin{equation} \boldsymbol{u} (\boldsymbol{x}, t) = \underbrace{ \langle \boldsymbol{u} (\boldsymbol{x},t)\rangle _T }_{\omega \ll \omega_c} + \underbrace{ \tilde{\boldsymbol{u}} ( \boldsymbol{x}, t) }_{\omega \sim \omega_c} + \underbrace{ \boldsymbol{u}^{\prime} ( \boldsymbol{x}, t) }_{\omega \gg \omega_c}, \end{equation}

where the dominant angular frequency $\omega _c$ is defined as the dominant peak in the Fourier spectrum of the velocity field. Here, the velocity field is decomposed into a slowly varying mean-flow field $\langle \boldsymbol {u}\rangle _T$, a coherent component on time scales of order $2{\rm \pi} /\omega _c$, involving coherent structures $\tilde {\boldsymbol {u}}$, and the remaining non-coherent small scale fluctuations $\boldsymbol {u}^{\prime }$. This kind of decomposition can also be found in the low-order Galerkin models of Tadmor et al. (Reference Tadmor, Lehmann, Noack, Cordier, Delville, Bonnet and Morzyński2011) and the weakly nonlinear modelling of Rigas, Morgans & Morrison (Reference Rigas, Morgans and Morrison2017a).

The slowly varying mean-flow field $\langle \boldsymbol {u}\rangle _T$ can be defined as the average of the velocity field $\boldsymbol {u}$ over one local period $T\approx 2{\rm \pi} /\omega _c$ of the coherent structures,

(3.2)\begin{equation} \langle \boldsymbol{u} ( \boldsymbol{x}, t ) \rangle _T := \frac{1}{T} \int_{t-T/2}^{t+T/2} \, {\rm d} \tau \boldsymbol{u} (\boldsymbol{x}, \tau ), \end{equation}

which eliminates both the coherent contribution from $\tilde {\boldsymbol {u}}$ and the non-coherent contribution from $\boldsymbol {u}^{\prime }$. Unlike the mean-flow field defined by the post-transient limit,

(3.3)\begin{equation} \bar{\boldsymbol{u}} (\boldsymbol{x}) = \lim_{T\rightarrow \infty }\frac{1}{T}\int_{0}^{T} \boldsymbol{u} (\boldsymbol{x}, \tau) \, {\rm d} \tau, \end{equation}

the finite-time-averaged flow field considered in this study owns a slowly varying dynamics. From the mean-field theory of Stuart (Reference Stuart1958), the slowly varying mean-flow field evolves out of the steady solution under the action of the Reynolds stress associated with the most unstable eigenmode(s). The mean-flow field deformation $\boldsymbol {u}_{\varDelta }$ is used to describe the difference between the slowly varying mean-flow field and the invariant steady solution $\boldsymbol {u}_s (\boldsymbol {x})$, which reads

(3.4)\begin{equation} \langle \boldsymbol{u} (\boldsymbol{x},t)\rangle _T = \boldsymbol{u}_s (\boldsymbol{x}) + \boldsymbol{u}_{\varDelta} (\boldsymbol{x},t). \end{equation}

3.1.2. Clustering algorithm

We consider the state vectors, for instance, the velocity fields $\boldsymbol {u} (\boldsymbol {x},t)$ in the computational domain $\varOmega$, which is sampled at times $t^m = m\Delta t$ with a time step $\Delta t$, where the superscript $m=1,\ldots,M$ is the snapshot index. The clustering process aims at partitioning the $M$ time-discrete states (snapshots) $\boldsymbol {u}^m = \boldsymbol {u} (\boldsymbol {x},t^m)$ into $K$ clusters $\mathcal {C}_k, k = 1,\ldots,K$. Snapshots of a given cluster share similar attributes featured by its cluster centroid $\boldsymbol {c}_k$. The distance between the snapshot $\boldsymbol {u}^m$ and the centroid $\boldsymbol {c}_k$ is defined as

(3.5)\begin{equation} D^{m}_{k}:= \|\boldsymbol{u}^m - \boldsymbol{c}_k\|_{\varOmega}. \end{equation}

Each snapshot is partitioned to the cluster of the closest centroid by $\hbox {argmin}_k D^{m}_{k}$, and the characteristic function is defined as

(3.6)\begin{equation} \chi_{k}^{m} := \begin{cases} 1, & \mbox{if }\boldsymbol{u}^m\in \mathcal{C}_k,\\ 0, & \mbox{otherwise}. \end{cases}\end{equation}

A cluster index $k^m, m=1,\ldots,M$, indicates the cluster assignment of the corresponding snapshot with $\boldsymbol {u}^m \in \mathcal {C}_k$, and records the visited clusters consecutively. The number of snapshots $n_k$ in cluster $k$ is given by

(3.7)\begin{equation} n_k = \sum_{m=1}^M \chi_{k}^{m}.\end{equation}

The cluster centroids $\boldsymbol {c}_k$ are defined as the average of the snapshots belonging to the cluster $\mathcal {C}_k$:

(3.8)\begin{equation} \boldsymbol{c}_k= \frac{1}{n_k} \sum_{m=1}^M \chi_{k}^{m} \boldsymbol{u}^m. \end{equation}

The performance of clustering is judged by the within-cluster variances

(3.9)\begin{equation} J (\boldsymbol{c}_1,\ldots,\boldsymbol{c}_K ) = \sum_{k=1}^K \sum_{m=1}^M \chi_{k}^{m}\Vert\boldsymbol{u}^m - \boldsymbol{c}_k \Vert_{\varOmega}^2. \end{equation}

The clustering algorithm minimizes $J$ and determines the optimal centroid positions,

(3.10)\begin{equation} \boldsymbol{c}_1^{{opt}}, \ldots,\boldsymbol{c}_K^{{opt}} = \hbox{argmin}_{\boldsymbol{c}_1,\ldots,\boldsymbol{c}_K}\, J (\boldsymbol{c}_1,\ldots,\boldsymbol{c}_K ), \end{equation}

by iteratively updating the characteristic function and the centroid positions.

To solve the optimization problem (3.10), we use the $k$-means$++$ algorithm (Arthur & Vassilvitskii Reference Arthur and Vassilvitskii2006). Comparing with the traditional $k$-means algorithm, the $k$-means$++$ algorithm selects the initial centroids as far away as possible to avoid any bias from the initial conditions. The remaining steps of the two algorithms are the same. At each iteration, the snapshots are divided into clusters of the nearest newly determined centroids. The optimal centroids are obtained by iterating until either convergence or when the maximum number of iterations is reached.

3.1.3. Cluster-based network model

Based on the clustering result, Kaiser et al. (Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014) derived a CMM, which provides a probabilistic representation of the system using a Markov process, with the assumption that the fluid system is memoryless. Nair et al. (Reference Nair, Yeh, Kaiser, Noack, Brunton and Taira2019) removed the transitions residing in the same cluster and emphasized the non-trivial transitions between two different clusters. In these two works, the transitions are only characterized by probabilities. The CNM proposed in Fernex et al. (Reference Fernex, Noack and Semaan2021) and Li et al. (Reference Li, Fernex, Semaan, Tan, Morzyński and Noack2021), inherited the idea of focusing on the non-trivial transitions, and further introduced time scale characteristics by recording the transition times. We here briefly review some concepts of the CNM, as they will be used in our benchmark of HiCNM.

The $M$ consecutive snapshots define $M-1$ transitions, containing trivial transitions staying in the same cluster and non-trivial transitions between two different clusters. The number of transitions from $\mathcal {C}_j$ to $\mathcal {C}_i$ reads

(3.11)\begin{equation} n_{ij} = \sum_{m=1}^{M-1} \chi_{j}^{m} \chi_{i}^{m+1}. \end{equation}

Considering the non-trivial transitions, ${n_{j}}$ is the total number of departing snapshots from $\mathcal {C}_j$, with $n_j = \sum _{i=1}^K (1-\delta _{ij}) n_{ij}$. The direct transition probability $P_{ij}$ reads

(3.12)\begin{equation} P_{ij}= \frac{(1-\delta_{ij}) n_{ij}}{n_{j}}, \quad i, j = 1, \ldots, K, \end{equation}

where the non-migrating transition ${n_{jj}}$ is eliminated. All the non-trivial transitions are identified with the direct transition matrix ${\boldsymbol{\mathsf{P}}}$.

The residence time matrix ${\boldsymbol{\mathsf{T}}}$ relies on the time information of the snapshots. After clustering, each snapshot $\boldsymbol {u}^m= \boldsymbol {u}(t^m)$ is associated with the closest centroid with a cluster index $k^m$. Assuming that $N$ $(N < M)$ non-trivial transitions occur along the trajectory, the moments of transition $t_n, n=1,\ldots,N$ – including the initial time $t_0=t^1$ – are defined as the time entering into a new cluster,

(3.13)\begin{equation} t_n = t^m \mbox{ if }\boldsymbol{u}^{m-1}\in \mathcal{C}_k \ \& \ \boldsymbol{u}^{m} \notin \mathcal{C}_k, \end{equation}

with ascending order $t_0 < t_1 < \cdots < t_N$. The cluster index $k^m$ remains unchanged in a time range $[t_n, t_{n+1} )$. The sequence of visited clusters over time can be simplified with the first entering snapshot for each non-trivial transition with $k^m$ taken from the moments of transition $t_n$.

For a simple transition from $\mathcal {C}_j$ to $\mathcal {C}_i$ as illustrated in figure 5, where the trajectory first enters in $\mathcal {C}_j$ at time $t_{n}$, and leaves $\mathcal {C}_j$ for $\mathcal {C}_i$ at time $t_{n+1}$, the residence time in $\mathcal {C}_j$ is defined as

(3.14)\begin{equation} T_{ij}:= t_{n+1} - t_{n}. \end{equation}

In the case of multiple trajectories of transition from $\mathcal {C}_j$ to $\mathcal {C}_i$, the residence time will be averaged according to the number of trajectories.

Figure 5. An illustration of the residence time in the cluster-based network model; $\bullet$ remarks the entering time into new clusters.

At this point, the time-resolved snapshots $\boldsymbol {u}^m$ can be represented by cluster centroids $\boldsymbol {c}_{k^m}$ with the time evolution of the cluster index $k^m$. The transient dynamics is described by both the direct transition matrix of non-migrating transitions ${\boldsymbol{\mathsf{P}}}$ and the residence time matrix ${\boldsymbol{\mathsf{T}}}$.

3.2.1. Hierarchical clustering inspired by the triple decomposition

To better understand the global and local properties of the data, we present a novel hierarchical clustering algorithm inspired by the triple decomposition introduced in § 3.1.1. The principle of hierarchical clustering is to divide the snapshots into layers of clusters. Snapshots belonging to clusters of the parent layer are further partitioned into clusters of the child layer. Hierarchical clustering algorithms are generally divided into two categories:

1. (a) The agglomerating (‘bottom-up’) hierarchical clustering begins with the smallest clusters at the bottom, each snapshot being an elementary cluster. The two closest clusters are merged to generate a new cluster according to certain criteria, introducing an additional layer from the bottom. This merging is repeated until all snapshots belong to one cluster at the top of the hierarchy.

2. (b) The divisive (‘top-down’) hierarchical clustering starts with only one cluster, which owns all the snapshots. In our case, the centroid of the top cluster would be the mean-flow field from ensemble averaging. From the top to the bottom, the snapshots in each cluster of the parent layer are divided into multiple clusters in the child layer, according to certain criteria. The bottom layers will be associated with the small scale fluctuations of the flow field. The division can be continued until each snapshot is a cluster.

We employ a divisive hierarchical clustering to distil the different features in a hierarchy, which is consistent with the triple decomposition of (3.1). For instance, fluid flows characterized by multiple frequencies require only a finite number of layers to describe the different components bounded by frequency.

Transient and post-transient dynamics are statistically non-homogeneous due to the existence of multiple invariant sets. If so, a scale sub-division of the flow-field decomposition like in (3.1) is used during the clustering process. Accounting for the Reynolds stress contribution, the slowly varying mean-flow field $\langle \boldsymbol {u}\rangle _T$ is enough to describe the global trend. Next, the local dynamics around $\langle \boldsymbol {u}\rangle _T$ can be zoomed in, considering the coherent structures involved in $\tilde {\boldsymbol {u}}$. This scale sub-division can still be extended to a hierarchical structure with more layers, which involve secondary frequencies in the case of quasi-periodic dynamics or turbulence from $\boldsymbol {u}^{\prime }$, as illustrated in figure 6. In order to clearly describe the clusters in the hierarchy, we systematically name the clusters from top to bottom. The sole cluster on the top $\mathcal {L}_0$ contains the ensemble of input data, and we define it symbolically as $\mathcal {C}_0$. The sub-division of this cluster leads to $K_1$ subclusters in the first layer $\mathcal {L}_1$, named as $\mathcal {C}_{0, k_1}, k_1=1,\ldots, K_1$. The first subscript $k_0 = 0$ can be ignored because there is only one cluster in $\mathcal {L}_0$, and the second subscript $k_1$ indicates the index of the subcluster in $\mathcal {L}_1$. The second sub-division works on each cluster $\mathcal {C}_{k_1}$ separately, and generates refined $K_2$ subclusters for each of them. The cluster index in the current layer $\mathcal {L}_2$ is presented by an additional subscript $k_2 = 1, \ldots, K_2$, which is written as $\mathcal {C}_{k_1, k_2}$. For a higher layer number $\mathcal {L}_{L \in \mathbb {N}}, L \geqslant 3$, more subscripts $k_L, l = 1, \ldots, L$, are needed to record the cluster index in each layer $\mathcal {L}_l$ from the top to the bottom, written as $\mathcal {C}_{k_1, \ldots, k_L}$. This naming method can clearly trace out all clusters in the hierarchy, and also works for other properties of clusters, e.g. the centroids $\boldsymbol {c}_{k_1, \ldots, k_L}$ and the characteristic function $\chi _{k_1, \ldots, k_L}^{m}$. In this work, two or three layers ($L \leqslant 3$) will be enough to extract the transient dynamics out of multiple invariant sets.

Figure 6. An illustration of the hierarchical structure with different scales in the triple flow decomposition in § 3.1.1. Layer 0: the top layer is characterized by the invariant mean flow $\bar {\boldsymbol {u}}$. Layer 1: the global trend is described by the slowly varying mean-flow field $\langle \boldsymbol {u}\rangle _T$. Layer 2: the coherent part $\tilde {\boldsymbol {u}}$ is added for the local dynamics around the varying mean-flow field. Layer 3: the non-coherent part $\boldsymbol {u}^{\prime }$ is considered in the case of turbulent flow.

3.2.2. Hierarchical network modelling

The starting point is the hierarchical Markov model of Fine et al. (Reference Fine, Singer and Tishby1998), which introduces the hierarchical structure to describe the stochastic processes, comparing with the standard Markov model. Each state of a Markov model in the parent layer is considered separately, and a new Markov model of the sub-states of a state is built in the child layer. As the layer increases, the state is continuously sub-divided. Therefore, the hierarchical Markov model records a sequence of states in different layers. In our case, each cluster is seen as a state. When a cluster in the parent layer is activated, its subclusters in the child layer turns activated recursively. Meanwhile, the refined dynamics between the subclusters can be described by a Markov model. Hence, the hierarchical Markov model can more effectively solve the problem of subsets.

In this work, we derive the HiCNM by replacing the Markov model by the network model. The hierarchical structure is identical, and the only change is the way to describe the transient dynamics between clusters.

A typical structure between the parent and child layers is shown in figure 7. We start with cluster $\mathcal {C}_{k_1 , \ldots, k_{l-1} }$ in the parent layer $\mathcal {L}_{l-1}$, where the leading subscript $k_1 , \ldots, k_{l-2}$ refers to the cluster number in each upper layer. For convenience, when the context is unambiguous, the cluster will be only referenced by its number in the current layer, e.g. $\mathcal {C}_{k_{l-1}}$. As indicated with the sequence of cluster numbers in its complete name, this cluster comes from the sub-division of the cluster $\mathcal {C}_{k_1 , \ldots, k_{l-2} }$ in the parent layer $\mathcal {L}_{l-2}$. We suppose that $M$ snapshots $\boldsymbol {u}^{m}, m=1,\ldots,M$, exist in this cluster and are divided into $K_{l-1}$ subclusters by a sub-division clustering algorithm. The subcluster $\mathcal {C}_{k_{l-1}}$ contains $n_{k_{l-1} }$ snapshots, calculated from (3.7) with the characteristic function $\chi _{k_{l-1} }^{m}$. A standard network model for the cluster $\mathcal {C}_{k_1 , \ldots, k_{l-2} }$ can be derived with the direct transition matrix ${\boldsymbol{\mathsf{P}}}_{k_1 , \ldots, k_{l-2} }$ and the residence time matrix ${\boldsymbol{\mathsf{T}}}_{ k_1 , \ldots, k_{l-2} }$, as recorded in § 3.1.3, which describe the dynamics between the subclusters $\mathcal {C}_{k_{l-1}}$.

Figure 7. An illustration of the transitions between the parent and child layers in the hierarchical network model. The trajectories pass through the cluster in the parent layer: the entering and exiting snapshots are marked with red dot and the blue dot. After clustering, a classic network model is built between $N$ subclusters, with transition probability $P$. The vertical transitions indicate the ports of entry and exit of the subclusters with probabilities $Q_{o,j}$ and $Q_{e,j}$.

In the following, we focus on the trajectories passing through cluster $\mathcal {C}_{k_{l-1}}$. The snapshots entering and leaving from $\mathcal {C}_{k_{l-1}}$ are marked out for each trajectory, with the following characteristic function:

(3.15) \begin{equation} \left.\begin{array}{c@{}} \chi_{o, k_{l-1}}^{m} := \begin{cases} 1, & \mbox{if }\boldsymbol{u}^{m-1} \notin \mathcal{C}_{k_{l-1}} \ \& \ \boldsymbol{u}^{m} \in \mathcal{C}_{k_{l-1}}, \\ 0, & \mbox{otherwise}. \end{cases} \\ \chi_{e, k_{l-1}}^{m} := \begin{cases} 1, & \mbox{if }\boldsymbol{u}^{m}\in \mathcal{C}_{k_{l-1}} \ \& \ \boldsymbol{u}^{m+1} \notin \mathcal{C}_{k_{l-1}}, \\ 0, & \mbox{otherwise}. \end{cases} \end{array}\right\} \end{equation}

The entering snapshots are denoted by the subscript ‘$o$’, and the exiting snapshots by the subscript ‘$e$’. The number of entering snapshots $n_o$ and of exiting snapshots $n_e$ read

(3.16a,b)\begin{equation} n_o = \sum_{m=1}^M \chi_{o,k_{l-1}}^{m}, \quad n_e = \sum_{m=1}^M \chi_{e,k_{l-1}}^{m}. \end{equation}

In the child layer $\mathcal {L}_{l}$, the $n_k$ snapshots in the cluster $\mathcal {C}_{k_{l-1}}$ have been divided into the subclusters $\mathcal {C}_{k_1 , \ldots, k_{l} }$ , $k_{l} =1,\ldots,k_L$. Without loss of generality, a standard network model for the cluster $\mathcal {C}_{k_{l-1}}$ can be built with its subclusters with the direct transition matrix ${\boldsymbol{\mathsf{P}}}_{k_1 , \ldots, k_{l-1} }$ and the residence time matrix ${\boldsymbol{\mathsf{T}}}_{ k_1 , \ldots, k_{l-1} }$.

The snapshots $\boldsymbol {u}^m$ are approximated by the time evolution of the cluster centroids $\boldsymbol {c}_{k_1^m, \ldots, k_{l}^m }$ in $\mathcal {L}_{l}$,

(3.17)\begin{equation} \hat{\boldsymbol{u}}^m_{\mathcal{L}_{l}} = \boldsymbol{c}_{k_1^m, \ldots , k_{l}^m }. \end{equation}

The residence time elements of ${\boldsymbol{\mathsf{T}}}_{k_1 , \ldots, k_{l}}$ can be assembled in order to determine the moments of transition based on the sequence of visited clusters in $\mathcal {L}_{l}$.

The entering and exiting snapshots defined in (3.15) can be used to describe the vertical transitions. Although they are not necessary to describe the dynamics of the fluidic pinball, the probability of the vertical transitions $Q_{o,j}$ and $Q_{e,j}$, described in Appendix B, completes all possible transitions in our hierarchical structure and makes it consistent with the classic hierarchical Markov model of Fine et al. (Reference Fine, Singer and Tishby1998).

3.2.3. Dynamics reconstruction of the hierarchical network model

The reconstructed flow in (3.17) is a statistical representation of the original snapshot sequence by a few representative centroids, which is a highly discretized description compared with the full dynamics. The approximations in the different layers provide different metrics for the flow dynamics.

The cluster-based hierarchical model uses the centroids $\boldsymbol {c}_{k_1^m, \ldots, k_{l}^m }$ in the original data space together with the time evolution of the cluster index $k_{l}^m$ to simplify the description of the original flow, which is more intuitive and closer to the original flow than the POD reconstruction. For input data with $I$-dimensional state vectors of the velocity field and $M$ snapshots, a POD reconstruction truncated to $R$ modes will lead to a $I \times R$ matrix of POD modes and a $R \times M$ matrix of mode amplitudes. A HiCNM with $K$ centroids leads to a $I \times K$ matrix of centroids and a sequence of the $N$ visited clusters of length $N \ll M$. In this sense, the compressive ability of HiCNM is more powerful for the large amount of continuously sampled data, as $I \times K < (I+M)\times R$. In addition, the hierarchical clustering works as a sparse sampling technique, extracting the representative states according to the clustering subspace.

We use the unbiased auto-correlation function (Protas, Noack & Östh Reference Protas, Noack and Östh2015),

(3.18)\begin{equation} R(\tau)=\frac{1}{T-\tau}\int_{\tau}^{T}(\boldsymbol{u}(\boldsymbol{x}, t-\tau) - \boldsymbol{u}_s(\boldsymbol{x}),\boldsymbol{u}(\boldsymbol{x}, t)- \boldsymbol{u}_s(\boldsymbol{x}))_{\varOmega} \mathrm{d}t, \quad \tau\in [0,T), \end{equation}

after normalization with respect to $R(0)$ to check the accuracy of the dynamics reconstruction of the HiCNM in (3.17). The autocorrelation function without delay $R(0)$ is twice the time-averaging kinetic energy. The modelled autocorrelation function $\hat {R}_{\mathcal {L}_{l}}(\tau )$ in layer $\mathcal {L}_{l}$ is based on the rebuilt flow $\hat {\boldsymbol {u}}_{\mathcal {L}_{l}}$ in (3.17) instead of $\boldsymbol {u}$ in (3.18).

For the discrete snapshots, the root mean-square error (RMSE) of the autocorrelation function $R(\tau )$ of the reference data and that of the model $\hat {R}_{\mathcal {L}_{l}}(\tau )$ is defined as

(3.19)\begin{equation} R_{rms}^l : = \sqrt{\frac{1}{M}\sum_{m=1}^{M}(R(\tau)-\hat{R}_{\mathcal{L}_{l}}(\tau))^2}, \end{equation}

where $M$ is the number of snapshots $\boldsymbol {u}^m$ and $\hat {\boldsymbol {u}}^m_{\mathcal {L}_{l}}$ at $\tau = m \Delta t$.

4.2.1. Hierarchical network model in layer 1

The $K_1=20$ clusters are used to cluster the snapshots from the low-pass filtered data set which removes the coherent structures with frequency $f_c = 0.1074$. To visualize the cluster topology, we apply the classical multidimensional scaling (MDS) to represent the high-dimensional centroids in a two-dimensional subspace $[\gamma _1, \gamma _2]^\textrm {T}$, while the distances between the centroids are preserved (Kaiser et al. Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014). As shown in figure 9, the six exact solutions of the Navier–Stokes equations that organize the state space are well identified. The vorticity field of the resulting centroids can be understood as the slowly varying mean-flow field $\langle \boldsymbol {u}\rangle _T$ along the transient dynamics, with $T\gg 2{\rm \pi} /\omega _c$. The three steady solutions belong respectively to clusters $\mathcal {C}_1$ (symmetric steady solution $\boldsymbol {u}_s$), $\mathcal {C}_{12}$ (asymmetric steady solution $\boldsymbol {u}_s^-$) and $\mathcal {C}_{17}$ (asymmetric steady solution $\boldsymbol {u}_s^+$). The three limit cycles are caught by the time-averaged flow in clusters $\mathcal {C}_7$ (symmetric mean-flow field, centroid $\bar {\boldsymbol {u}}^0$), $\mathcal {C}_9$ (asymmetric upward mean-flow field, centroid $\bar {\boldsymbol {u}}^+$) and $\mathcal {C}_{11}$ (asymmetric downward mean-flow field, centroid $\bar {\boldsymbol {u}}^-$). A network of four transient trajectories connects these clusters as follows:

1. ${Trajectory}$ 1: $\mathcal {C}_1\ (\boldsymbol {u}_s) \rightarrow \cdots \rightarrow \mathcal {C}_7\ (\bar {\boldsymbol {u}}^0)\ \rightarrow \mathcal {C}_8\ \rightarrow \mathcal {C}_9\ (\bar {\boldsymbol {u}}^+)$;

2. ${Trajectory}$ 2: $\mathcal {C}_1\ (\boldsymbol {u}_s) \rightarrow \cdots \rightarrow \mathcal {C}_7\ (\bar {\boldsymbol {u}}^0)\ \rightarrow \mathcal {C}_{10} \rightarrow \mathcal {C}_{11}\ (\bar {\boldsymbol {u}}^-)$;

3. ${Trajectory}$ 3: $\mathcal {C}_{12}\ (\boldsymbol {u}_s^-) \rightarrow \cdots \rightarrow \mathcal {C}_{16} \rightarrow \mathcal {C}_{10} \rightarrow \mathcal {C}_{11}\ (\bar {\boldsymbol {u}}^-)$;

4. ${Trajectory}$ 4: $\mathcal {C}_{17}\ (\boldsymbol {u}_s^+) \rightarrow \cdots \rightarrow \mathcal {C}_{20} \rightarrow \mathcal {C}_8\ \rightarrow \mathcal {C}_9\ (\bar {\boldsymbol {u}}^+)$.

Figure 9. Graph of transitions between clusters in layer 1 at $Re= 80$. Cluster centroids are marked with the coloured squares, with their vorticity fields in colour with $[-1.5, 1.5]$. The snapshots belonging to them are marked as small dots with the same colours. The transitions between clusters are shown with arrows, where the line width presents the probability of transition.

We notice that most of the transitions between clusters are with $100\,\%$ probability, except for the bifurcating cluster $\mathcal {C}_{7}$ with half-probability of transition to clusters $\mathcal {C}_{8}$ or $\mathcal {C}_{10}$. The transition matrix in figure 10(b) illustrates the probability of all the transitions. The probability is $1$ for the surely directed transitions. By contrast, $P_{7\, 8}=0.5$ and $P_{7\, 10}=0.5$ for the bifurcating cluster. After entering into $\mathcal {C}_{8}$ and $\mathcal {C}_{10}$, the flow will surely enter the two clusters $\mathcal {C}_{9}$ and $\mathcal {C}_{11}$, respectively. The two clusters $\mathcal {C}_{9}$ and $\mathcal {C}_{11}$ catch the permanent regimes, from which the dynamics cannot escape, imposing all the terms in the ninth and eleventh columns $P_{j\, 9} = P_{j\, 11} = 0, \forall j$.

Figure 10. Cluster-based analysis at $Re= 80$ in layer 1: (a) transition illustrated with cluster label, (b) transition matrix, (c) residence time matrix. Since subclusters are not considered, the snapshot index in layer 1 $m_1$ is identical to the original index $m$. The colour bar indicates the values of the terms. Residence time larger than 100 is marked with a solid black circle and excluded from the colour bar. An extremely long residence time in a cluster indicates data density, and such a cluster is generally associated with an invariant set.

The time ordering of the cluster transitions is illustrated in figure 10(a). The red vertical line separates the four trajectories. All the trajectories are irreversible transition from one steady solution to one of the two stable periodic solutions. In figure 10(c), the filled black circles emphasize the transitions starting from $\mathcal {C}_{1}, \mathcal {C}_{12}$, $\mathcal {C}_{17}$ and $\mathcal {C}_{7}$. The two attracting clusters $\mathcal {C}_{9}$ and $\mathcal {C}_{11}$ have no transition to any other clusters. Hence, all the terms in the ninth and eleventh columns $T_{j\, 9}, T_{j\, 11}$ are 0, $\forall j$. The residence time associated with clusters $\mathcal {C}_{9}$ and $\mathcal {C}_{11}$ is infinite.

According to the above discussion, the network model in this layer has successfully identified the six invariant sets of the dynamics, four being unstable, the two others being the attractors of the system.

4.2.2. Hierarchical network model in layer 2

We apply Algorithm 2 based on the clustering result for the layer $\mathcal {L}_{1}$ in § 4.2.1. In the second layer $\mathcal {L}_{2}$, we will isolate and analyse the clusters $\mathcal {C}_{k_1}, k_1=9, 7, 1$ associated with three invariant sets of the dynamics.

The permanent regime in cluster $\mathcal {C}_{9}$

Cluster $\mathcal {C}_{9}$ is associated with one of the two asymmetric limit cycles. The $k$-means$++$ algorithm is directly applied to the intra-cluster snapshots $\boldsymbol {u}^m \in \mathcal {C}_{9}$. The limit cycle in figure 11(a) has been divided into $K_2=10$ subclusters according to Algorithm 2. The centroids are distributed on the limit cycle at equal distances. The arrows between the centroids form a closed loop, which results from the periodic nature of the oscillating dynamics. The resulting centroids are phase-averaged flow fields along the complete period of the vortex shedding. The innerjet flows of the centroids in figure 11(a) are all deflected upwards, as expected for the attractor that belongs to cluster $\mathcal {C}_{9}$.

Figure 11. Cluster-based analysis in layer 2 at $Re= 80$ for $\mathcal {C}_{9}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories pass through $\mathcal {C}_{9}$ in the parent layer, one with $m_2=1, \ldots, 5430$ and another with $m_2=5431, \ldots, 9917$.

Figure 11(b) associates the cluster labels with the dynamics. Two different transient trajectories reach $\mathcal {C}_{9}$ with entering snapshots $m = 9571$ and $45\ 514$, one issued from $\mathcal {C}_{8}$, the other from $\mathcal {C}_{16}$, according to the network model in layer $\mathcal {L}_{1}$. According to the entering time, the original snapshot index $m = m_2 + 9571 - 1 = 9571, \ldots, 15\ 000$ and $m = m_2 + 45\ 514 - 5431 = 45\ 514, \ldots, 50\ 000$. The flow periodically travels along the subclusters $\mathcal {C}_{9, k_2}, k_2 = 1, \ldots, 10$ as $\mathcal {C}_{9,\,1} \rightarrow \cdots \rightarrow \mathcal {C}_{9,\,10} \rightarrow \mathcal {C}_{9,\,1}$. The limit cycle has a clear and stable transition matrix, as each cluster only has one possible destination, see figure 11(c). The residence times in each clusters are uniform as shown in figure 11(d). The sum of all residence times is $9.51$, which is close to the real time period $9.50$ of the vortex shedding in the permanent regime computed from the DNS.

The bifurcating state in cluster $\mathcal {C}_{7}$

Cluster $\mathcal {C}_{7}$ is associated with the symmetric limit cycle. $K_2=10$ clusters are used to classify the snapshots in this cluster.

The resulting limit cycle of figure 12(a) differs from the limit cycle of figure 11 by its centroids. As expected with the symmetric limit cycle, the inner jet is not deflected in cluster $\mathcal {C}_{7}$, while it is deflected in cluster $\mathcal {C}_{9}$. In figure 12(b), two different transient trajectories pass through $\mathcal {C}_{7}$ with entering snapshots $m = 6541$ and $21\ 552$, the original snapshot index is $m = 6541, \ldots, 8870$ and $m = 21\ 552, \ldots, 23\ 886$. The limit cycle has a stable transition matrix, as shown in figure 12(c). The sum of the residence times of figure 12(d) is $9.79$, again very close to the period $9.80$ of the symmetric transient vortex shedding computed from the DNS. The bifurcating dynamics cannot be detected with these subclusters but can be captured in the parent layer.

Figure 12. Cluster-based analysis in layer 2 at $Re= 80$ for $\mathcal {C}_{7}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories pass through $\mathcal {C}_{7}$ in the parent layer, one with $m_2=1, \ldots, 2330$ and another with $m_2=2331, \ldots, 4665$.

The destabilizing regime in cluster $\mathcal {C}_{1}$

Cluster $\mathcal {C}_{1}$ is associated with two mirror-conjugated trajectories spiralling out of the symmetric steady solution. $K_2=10$ subclusters are used in the child layer $\mathcal {L}_{2}$. In the $[\gamma _1, \gamma _2]^\textrm {T}$ representation of figure 13(a), the centroids are distributed along two diverging trajectories spiralling out of the fixed point $[0,0]^\textrm {T}$.

Figure 13. Cluster-based analysis in layer 2 at $Re= 80$ for $\mathcal {C}_{1}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories pass through $\mathcal {C}_{1}$ in the parent layer, one with $m_2=1, \ldots, 5319$ and another with $m_2=5320, \ldots, 10\ 640$.

The first three subclusters $\mathcal {C}_{1,\,1}, \mathcal {C}_{1,\,2}$ and $\mathcal {C}_{1,\,3}$ belong to the inner zone of the spirals, where the distribution of snapshots is dense in the $[\gamma _1, \gamma _2]^\textrm {T}$ proximity map of figure 13. As a result, three non-physical closed-loop cycles are formed between these three clusters, namely $\mathcal {C}_{1,\,1} \rightarrow \mathcal {C}_{1,\,2} \rightarrow \mathcal {C}_{1,\,1}, \mathcal {C}_{1,\,1} \rightarrow \mathcal {C}_{1,\,3} \rightarrow \mathcal {C}_{1,\,1}$ and $\mathcal {C}_{1,\,3} \rightarrow \mathcal {C}_{1,\,1} \rightarrow \mathcal {C}_{1,\,2} \rightarrow \mathcal {C}_{1,\,1} \rightarrow \mathcal {C}_{1,\,3}$. The remaining clusters belong to the outer arms of the spirals, with a relatively sparse distribution of the snapshots. The transitions between them form a closed-loop trajectory $\mathcal {C}_{1,\,4} \rightarrow \cdots \rightarrow \mathcal {C}_{1,\,10} \rightarrow \mathcal {C}_{1,\,4}$. The transitions $\mathcal {C}_{1,\,2} \rightarrow \mathcal {C}_{1,\,4}$ and $\mathcal {C}_{1,\,3} \rightarrow \mathcal {C}_{1,\,8}$ correspond to the departing dynamics out of the inner zone, due to the growth of the instability. The varying density of distribution comes from the exponential growth of the instability. The flow perturbations are small in the beginning of the instability while the flow distortions evolve faster in the later stages, leading to a multiscale problem in the transient and post-transient flow dynamics. In this case, the later stages with larger distortions are obviously easier to divide into different clusters.

From figure 13(b), two transient trajectories pass through $\mathcal {C}_{1}$ with entering snapshots $m = 1$ and $15\ 001$. The original snapshot index is $m = 1, \ldots, 5319$ and $m = 15\ 001, \ldots, 20\ 321$, corresponding to the initial stage of destabilization from the symmetric steady solution. In the transition matrix of figure 13(c), the inner and outer portions of the spiral are also apparent. The inner zone is the oscillating dynamics between $\mathcal {C}_{1,\,1}$ and $\mathcal {C}_{1,\,2}$, or $\mathcal {C}_{1,\,1}$ and $\mathcal {C}_{1,\,3}$. The outer zone has a more obvious periodic dynamics through the remaining clusters from $\mathcal {C}_{1,\,4}$ to $\mathcal {C}_{1,\,10}$. In figure 13(d), the residence time in each cluster is very short compared with the residence time in $\mathcal {C}_{1,\,1}$, which the vicinity of the steady solution belongs to.

4.2.3. Dynamics reconstruction of the hierarchical network model at $Re=80$

Figure 14 shows the autocorrelation function of the DNS and the HiCNM in different layers. As the autocorrelation function has been normalized by $R(0)$, the unit one presents the level of kinetic energy of the whole transition. For the transient and post-transient dynamics, the autocorrelation function vanishes with increasingtime shift, as shown in figure 14(b). The autocorrelation function of the DNS identifies the dominant frequency. In layer $\mathcal {L}_{1}$, no oscillation can be identified, due to the centroids by averaging the snapshots within clusters in the state space. The RMSE of the autocorrelation function is $R_{rms}^1= 17.46$. In layer $\mathcal {L}_{2}$, the autocorrelation function of the model matches perfectly over the entire range, with $R_{rms}^2=1.18$, which quantifies the accuracy of the cluster-based model.

Figure 14. Autocorrelation function for $\tau \in [0, 1500)$ from DNS (black solid line) and the hierarchical network model (red dashed line) in the two layers: (a) $\mathcal {L}_{1}$ and (b) $\mathcal {L}_{2}$, at $Re=80$.

4.3.1. Standard CNM at $Re= 80$

We show a standard network model treating the transient dynamics at $Re=80$ with $K=200$ clusters. The directed graph is not shown here as too many clusters overlapped in the two-dimensional subspace $[\gamma _1, \gamma _2]^\textrm {T}$, losing the interpretation of the dynamics. In figure 15(a), the six exact solutions are well classified. The three steady solutions are divided into $\mathcal {C}_1, \mathcal {C}_{192}$ and $\mathcal {C}_{196}$ separately. The three limit cycles are identified with three blocks of oscillating labels, from $\mathcal {C}_{69}$ to $\mathcal {C}_{99}$, from $\mathcal {C}_{100}$ to $\mathcal {C}_{142}$ and from $\mathcal {C}_{148}$ to $\mathcal {C}_{191}$. However, the transient trajectories starting with the asymmetric steady solutions are misidentified, as both travel through the same clusters from $\mathcal {C}_{2}$ to $\mathcal {C}_{68}$ and even the block of the symmetric limit cycle. It indicates that the clustering algorithm failed to distinguish the symmetry breaking in the transient dynamics starting with the different steady solutions. Too many clusters also make the transition matrix hard to read as illustrated in figure 15(c), and also for the residence time matrix not shown here. Ignoring the transitions with low probability in the transition matrix, the transitions in the above-mentioned three blocks are almost definite, indicating the correct identification of the cycles. However, the transient dynamics from $\mathcal {C}_{2}$ to $\mathcal {C}_{68}$ is random. We can hardly find any relevant feature for the building of vortex shedding or symmetric breaking. The cluster distribution directly affects the analysis of the dynamics.

Figure 15. Cluster-based analysis at $Re= 80$ for the fluctuating flow: transition illustrated with cluster label and transition matrix of the standard network model (a,c) and the hierarchical network model in the second layer (b,d), as in figure 10.

In summary, a large number of clusters tends to increase the resolution of the identified network model, but will misidentify the dynamics with the random transitions between clusters during the transient state. The standard network model fails to describe the transient dynamics between different invariant sets. This kind of problem comes from the poor distribution of clusters and can be solved by the hierarchical clustering strategy, as shown in the following sub-section.

4.3.2. HiCNM at $Re= 80$ in layer 2

Based on the hierarchical network model at $Re=80$ recorded in § 4.2, we build a network model ensembling all the clusters in the layer $\mathcal {L}_2$, which also contains $K=200$ clusters. For ease of illustration, the cluster indexes of two layers $\mathcal {C}_{k_1, k_2}$ are denoted with a single index $\mathcal {C}_{k}$, with $k=1,\ldots,200$.

From figure 15(b), we found that the cluster distribution is very uniform during the transient states and the post-transient states. The six exact solutions are well classified, and much fewer clusters are used for the limit cycles. The three steady solutions being divided into $\mathcal {C}_1, \mathcal {C}_{111}$ and $\mathcal {C}_{161}$ separately. The three limit cycles are identified with three blocks of oscillating labels, from $\mathcal {C}_{61}$ to $\mathcal {C}_{70}$, from $\mathcal {C}_{81}$ to $\mathcal {C}_{90}$ and from $\mathcal {C}_{101}$ to $\mathcal {C}_{110}$. Even the transient states between the limit cycles can also be identified, with two blocks from $\mathcal {C}_{71}$ to $\mathcal {C}_{80}$, and from $\mathcal {C}_{91}$ to $\mathcal {C}_{100}$. The transient trajectories starting with three different steady solutions are also well separated. The transition matrix in figure 15(d) are almost full of definite transitions, which shows clearer transient dynamics than figure 15(c). The global matrix keeps the local dynamics for each cluster $\mathcal {C}_{k_1}$ in the first layer $\mathcal {L}_{1}$, as shown in each small block of 10 clusters, and ensembles them together. The transitions from block to block have a much lower probability, comparing with the cycling transitions within the block, which makes them hard to see in the figure. However, this kind of transition between the blocks should be also definite, as shown in § 4.2.1 for the hierarchical model in the first layer. Analogously, the multiscale problem of the dynamics in different layers also exists in the global residence time matrix, which makes the small scale terms unseeable. Therefore, we suggest analysing the slowly varying mean flow and the local dynamics separately in different layers, as in § 4.2, to avoid the influence of different scales.

4.4.1. Hierarchical network model in layer 1

The frequency of the coherent component $\tilde {\boldsymbol {u}}$ is $f_c = 0.1172$. The $K_1=11$ clusters are used on the low-pass filtered data set, following Algorithm 1. The non-trivial transitions are shown in the two-dimensional subspace $[\gamma _1, \gamma _2]^\textrm {T}$ of figure 16. The two attractors belong to clusters $\mathcal {C}_6$ and $\mathcal {C}_{11}$. The symmetric and asymmetric steady solutions belong respectively to clusters $\mathcal {C}_1, \mathcal {C}_2$ and $\mathcal {C}_{7}$. The transient dynamics observed at $Re=105$ is different from that identified in figure 9 for $Re=80$. Figure 16 shows four trajectories, initiated from mirror-conjugated initial conditions close to the symmetric and asymmetric steady solutions:

1. ${Trajectory}$ 1: $\mathcal {C}_1\ (\boldsymbol {u}_s) \rightarrow \mathcal {C}_2 (\boldsymbol {u}_s^+) \rightarrow \cdots \rightarrow \mathcal {C}_6 (\bar {\boldsymbol {u}}^+)$;

2. ${Trajectory}$ 2: $\mathcal {C}_1\ (\boldsymbol {u}_s) \rightarrow \mathcal {C}_7 (\boldsymbol {u}_s^-) \rightarrow \cdots \rightarrow \mathcal {C}_{11} (\bar {\boldsymbol {u}}^-)$;

3. ${Trajectory}$ 3: $\mathcal {C}_2 (\boldsymbol {u}_s^+) \rightarrow \cdots \rightarrow \mathcal {C}_6 (\bar {\boldsymbol {u}}^+)$;

4. ${Trajectory}$ 4: $\mathcal {C}_7 (\boldsymbol {u}_s^-) \rightarrow \cdots \rightarrow \mathcal {C}_{11} (\bar {\boldsymbol {u}}^-)$.

Figure 16. Graph of transitions between clusters in layer 1 at $Re= 105$, displayed as in figure 9.

The state trajectories start from the symmetric steady solution to one of the two asymmetric steady solutions, then converge to the corresponding attracting torus. The two trajectories on the left side have a significant phase delay, while the two on the right side are almost in the same phase. This phase difference is a random function that depends on the initial condition. The two tori also look different because the feature vectors associated with $[\gamma _1, \gamma _2]^\textrm {T}$ are asymmetrical from the MDS.

In figure 17(a), the evolution of the cluster label of the snapshots illustrates four irreversible transient dynamics. The transition matrix in figure 17(b) exhibits two red diagonals associated with trajectories to the final state, and two yellow elements associated with the initial state starting close to the symmetric steady solutions in cluster $\mathcal {C}_1$. The transition from $\mathcal {C}_1$ to clusters $\mathcal {C}_{2}$ and $\mathcal {C}_{7}$ each have 1/2 probability. The two attracting clusters $\mathcal {C}_{6}$ and $\mathcal {C}_{11}$, associated with the attractors of the system, have no transition to any other clusters, as all the terms in the sixth and eleventh columns $P_{j\, 6} = P_{j\, 11} = 0, \forall j$. In the residence time matrix of figure 17(c), the filled black circles indicate three typical clusters $\mathcal {C}_{1}, \mathcal {C}_{2}$ and $\mathcal {C}_{7}$, which correspond to the vicinity of the three steady solutions. All the terms in the sixth and eleventh columns $T_{j\, 6}, T_{j\, 11}$ are 0 $\forall j$, which means that the residence time is infinite, as expected for attractors.

Figure 17. Cluster-based analysis at $Re= 105$ in layer 1: (a) transition illustrated with cluster label, (b) transition matrix, (c) residence time matrix, displayed as in figure 10.

4.4.2. Hierarchical network model in layer 2

We focus on the permanent regime in the cluster $\mathcal {C}_6$, and apply the sub-division clustering Algorithm 2, which results in $K_2=10$ subclusters in layer $\mathcal {L}_{2}$.

As illustrated in figure 18(a), a closed orbit between the 10 clusters is found. Figure 18(b) shows two transient trajectories pass through $\mathcal {C}_{6}$, the original snapshot index is $m = 8124, \ldots, 15\ 000$ and $m = 44\ 196, \ldots, 50\ 000$, involving the asymptotic regime of an attractor. The clusters in this closed orbit have a clear transition rule, as is evidenced by the transition matrix of figure 18(c). The sum of the residence times is $8.45$ for the cycle of $\mathcal {C}_3 \rightarrow \mathcal {C}_4 \rightarrow \cdots \rightarrow \mathcal {C}_{10} \rightarrow \mathcal {C}_{3}$ from the residence time matrix in figure 18(d).

Figure 18. Cluster-based analysis in layer 2 at $Re= 105$ for $\mathcal {C}_{6}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories pass through $\mathcal {C}_{6}$ in the parent layer, one with $m_2=1, \ldots, 6877$ and another with $m_2=6878, \ldots, 12\ 682$.

4.4.3. Hierarchical network model in layer 3

Based on the detected cycle found in layer $\mathcal {L}_{2}$, the entering snapshots of one subcluster can be used to sample the quasi-periodic regime. The snapshots recurrently enter into each cluster of layer $\mathcal {L}_{2}$. The entering states, defined by $T_{o,\,k_1,\,k_2}^{m}$ in (3.15), are the entry in the clusters. The entering snapshots of a given cluster can be considered as hits in a ‘Poincaré section’ (Guckenheimer & Holmes Reference Guckenheimer and Holmes2013).

The clustering algorithm in the layer $\mathcal {L}_{3}$ is applied to the snapshots of all the Poincaré sections. We still apply the sub-division clustering Algorithm 2, with $T_{o,\,k_1,\,k_2}^{m}$ as characteristic function.

We consider the entry in cluster $\mathcal {C}_{6,\,10}$ for illustration, and build a network model in the third layer $\mathcal {L}_{3}$. In the following, the cluster symbol in the first two layers $\mathcal {C}_{k_1=6,\,k_2 = 10, k_3}$ is omitted.

As shown by the graph of non-trivial transitions shown in figure 19(a), there exists a cycle $\mathcal {C}_3 \rightarrow \mathcal {C}_4 \rightarrow \cdots \rightarrow \mathcal {C}_{10} \rightarrow \mathcal {C}_{3}$. The periodically changing cluster label in figure 19(b) indicates the existence of a recurrent dynamics. The original cluster index corresponds to two sets of discrete snapshots in the Poincaré section, with interval approximate to the periodic detected in the $\mathcal {L}_{2}$. From the transition matrix of figure 19(c), the two clusters $\mathcal {C}_1$ and $\mathcal {C}_2$ are not part of the cycle, but they form the transient part of the dynamics, before entering the cycle. This indicates that the low-frequency modulations only start after the high-frequency oscillations have started in the second layer. In other words, the low frequency does not exist during the building process of the vortex shedding, but appears after the vortex shedding has developed to a certain degree.

Figure 19. Cluster-based analysis in layer 3 at $Re= 105$ for $\mathcal {C}_{6,\,10}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories are detected in the ‘Poincaré section’ of $\mathcal {C}_{6,\,10}$, one with $m_3=1, \ldots, 81$ and another with $m_3=82, \ldots, 150$.

From the residence time matrix of figure 19(d), the number of snapshots that belong to the cycle in the Poincaré section, determined by averaging over multiple trajectories, is $11.41$. The clustering analysis in the second layer $\mathcal {L}_{2}$ indicates that the trajectories periodically hit the Poincaré section with $8.45$, by summing up the elements in the blocks from $\mathcal {C}_3$ to $\mathcal {C}_{10}$ in figure 19(d). Hence, the resulting period of the cyclic process, considering both frequencies, is around $96.41$, which is very close to the real period of $97.10$ determined from the DNS.

4.4.4. Dynamics reconstruction of the hierarchical network model at $Re=105$

Figure 20 shows the autocorrelation function of the DNS and the HiCNM in the three layers. The autocorrelation function of the DNS identifies the two dominant frequencies of the dynamics. In layer $\mathcal {L}_{1}$, no oscillation can be identified, and the RMSE of the autocorrelation function is $R_{rms}^1=22.20$. In layer $\mathcal {L}_{2}$, the autocorrelation function of the model matches well with the high-frequency oscillations. The low-frequency oscillations can be also found, but the amplitude does not fit well. The error is $R_{rms}^2=1.52$, which is good enough for the accuracy. In layer $\mathcal {L}_{3}$, the amplitude of the low-frequency oscillations can be better reproduced. The error is further reduced to $R_{rms}^3=0.77$ with higher accuracy.

Figure 20. Autocorrelation function for $\tau \in [0, 1500)$ from DNS (black solid line) and the hierarchical network model (red dashed line) in the three layers: (a) $\mathcal {L}_{1}$, (b) $\mathcal {L}_{2}$ and (c) $\mathcal {L}_{3}$, at $Re=105$.

4.5.1. Hierarchical network model in layer 1

In the first layer, the filtered data set is used to analyse the mean-field dynamics with $K_1=11$ clusters. The detected frequency associated with the coherent component $\tilde {\boldsymbol {u}}$ is $f_c = 0.1225$. The clustering Algorithm 1 is applied to the data set, and the non-trivial transitions are shown in the two-dimensional subspace $[\gamma _1, \gamma _2]^\textrm {T}$ of figure 21. Four trajectories are found, each issued from one of the three steady solutions. The symmetric steady solution and the two asymmetric steady solutions respectively belong to clusters $\mathcal {C}_1, \mathcal {C}_{8}$ and $\mathcal {C}_{10}$. These clusters evolve through $\mathcal {C}_1 \rightarrow \mathcal {C}_2 \rightarrow \mathcal {C}_3, \mathcal {C}_8 \rightarrow \mathcal {C}_9$ and $\mathcal {C}_{10} \rightarrow \mathcal {C}_{11}$ before entering into the same chaotic cloud, consisting of the remaining clusters.

Figure 21. Graph of transitions between clusters in layer 1 at $Re= 130$, displayed as in figure 9.

There is no obvious periodic block of oscillating dynamics in the transient dynamics, as illustrated by the cluster labels of figure 22(a). However, the initial destabilizing process is characterized by a very long residence time in the clusters to which the steady solutions belong. From the transition matrix of figure 22(b), three transitions lead to the chaotic region with $100\,\%$ probability, as $P_{3\,4}, P_{9\,6}$ and $P_{11\, 4}$. Each cluster in the chaotic cloud has at least two possible destinations, with nearly equal probability. The hidden transition dynamics for this chaotic regime will be analysed in the next layer. The residence time matrix of figure 22(c) shows four black filled circles associated with clusters $\mathcal {C}_1, \mathcal {C}_{8}, \mathcal {C}_{10}$ and $\mathcal {C}_2$. The steady solutions belong to the first three, while the symmetric limit cycle belongs to $\mathcal {C}_2$, as it will become clear in the next sections.

Figure 22. Cluster-based analysis at $Re= 130$ in layer 1: (a) transition illustrated with cluster label, (b) transition matrix, (c) residence time matrix, displayed as in figure 10.

4.5.2. Hierarchical network model in layer 2

In the second layer $\mathcal {L}_{2}$, we focus on the clusters associated with four typical states detected in layer $\mathcal {L}_{1}$: the destabilizing state from the symmetric steady solution in $\mathcal {C}_{1}$, the destabilizing state from the (upward) asymmetric steady solution in $\mathcal {C}_{10}$, the transient state with long residence time before chaos in $\mathcal {C}_{2}$ and the chaotic state in the group of clusters $\mathcal {C}_{4},\ldots,\mathcal {C}_{7}$.

The destabilizing state in the cluster $\mathcal {C}_{1}$

Cluster $\mathcal {C}_{1}$ gathers snapshots in the initial stage of the instability starting from the symmetric steady solution. In the second layer $\mathcal {L}_{2}$, the snapshots of $\mathcal {C}_{1}$ are dispatched into $K_2 = 10$ subclusters $\mathcal {C}_{1,\,k_2}$, with $k_2 = 1,\ldots,K_2$. Figure 23 shows the transient trajectories with the centroids in the $[\gamma _1, \gamma _2]^\textrm {T}$ plane.

Figure 23. Cluster-based analysis in layer 2 at $Re= 130$ for $\mathcal {C}_{1}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories pass through $\mathcal {C}_{1}$ in the parent layer, one with $m_2=1, \ldots, 4878$ and another with $m_2=4879, \ldots, 9760$.

Similar to figure 13(a), the snapshots of figure 23(a) form two diverging trajectories spiralling out of the centre $[0,0]$. A loop is formed between the subclusters $\mathcal {C}_{1,\,1}$ and $\mathcal {C}_{1,\,2}$ in the inner zone. In the outer zone, a cycle appears with the periodic trajectory $\mathcal {C}_{1,\,3} \rightarrow \cdots \rightarrow \mathcal {C}_{1,\,10} \rightarrow \mathcal {C}_{1,\,3}$. Figure 23(b) shows two trajectories leaving $\mathcal {C}_{1,\,1}$. The original snapshot index is $m = 1, \ldots, 4878$ and $m = 15\ 001, \ldots, 19\ 882$, corresponding to the initial stage of the instability.

The transition matrix in figure 23(c) corroborates this periodic cycle. The black filled circle in figure 23(d) marks out the transition with a long residence time, due to the unstable centre that belongs to cluster $\mathcal {C}_{1,\,1}$.

The transitions $\mathcal {C}_{1,\,1} \rightarrow \mathcal {C}_{1,\,8}, \mathcal {C}_{1,\,2} \rightarrow \mathcal {C}_{1,\,3}$ and $\mathcal {C}_{1,\,2} \rightarrow \mathcal {C}_{1,\,10}$ correspond to the departing dynamics out of the inner zone, due to the development of the instability. We also note the returning transitions $\mathcal {C}_{1,\,8} \rightarrow \mathcal {C}_{1,\,2}$ and $\mathcal {C}_{1,\,9} \rightarrow \mathcal {C}_{1,\,2}$. However, the latter do not mean that the flow actually returns back to the destabilizing centre, as both trajectories are spiralling out of the centre. This confusing result comes from the clustering process. The edge between the inner and outer zones is not well defined, due to the varying density distribution of snapshots along the arms of the spirals. Cluster $\mathcal {C}_{1,\,2}$ overlaps the outer zone, to the difference of cluster $\mathcal {C}_{1,\,1}$, which fully belong to the inner zone. If we ignore the loop in the inner zone and merge $\mathcal {C}_{1,\,1}$ and $\mathcal {C}_{1,\,2}$, it shows a dynamical evolution from one cluster to the cycle of a group of clusters.

The destabilizing state in cluster $\mathcal {C}_{10}$

Cluster $\mathcal {C}_{10}$ contains the trajectory spiralling out from the upward-deflected asymmetric steady solution, as shown in figure 24(a). The snapshots are dispatched into $K_2 = 10$ of subclusters $\mathcal {C}_{10,\,k_2}$, with $k_2 = 1,\ldots,K_2$. Together with figure 24(b), it shows a one-way transition, departing from the unstable centre with a sparse spiral. The original snapshot index is $m = 40\ 001, \ldots,42\ 782$ for the initial stage of the destabilization from one asymmetric steady solution.

Figure 24. Cluster-based analysis in layer 2 at $Re= 130$ for $\mathcal {C}_{10}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. A sole trajectory passes through $\mathcal {C}_{10}$ in the parent layer with $m_2=1, \ldots, 2782$.

The graph of figure 24(a) is different from the graphs of figures 13(a) and 23(a), the latter being associated with the symmetric steady solution. The flow destabilization from the asymmetric steady solution develops faster than from the symmetric steady solution, as illustrated by the linear growth rates $\sigma _{sym}=0.032$ and $\sigma _{asym}=0.106$ of the respective pairs of unstable eigenmodes. As a result, the distribution of snapshots is sparser in figure 24(a) than in figure 23(a). As indicated by figure 24(b), 93.2 % of the 2782 snapshots of cluster $\mathcal {C}_{10}$ belongs to subcluster $\mathcal {C}_{10,\,1}$. The flow quickly travels through all the remaining clusters in only 188 snapshots.

The centroids of figures 23(a) and 24(a) have two main differences: (i) the inner jet is symmetric in the centroids of figure 23(a) while it is deflected upwards in the centroids of figure 24(a); (ii) the von Kármán street of vortices of figure 23(a) exhibits positive and negative vortices well apart from each other in the $y$-axis. In figure 24(a), the positive and negative vortices are of a larger strength and adjacent to the $x$-axis, together with a longer shear layer.

Transient regime before chaos in cluster $\mathcal {C}_{2}$

Cluster $\mathcal {C}_{2}$ contains two transient trajectories which connect the symmetric steady solution in $\mathcal {C}_{1}$ to the chaotic cloud, with the original snapshot index $m = 4879, \ldots,5966$ and $19\ 883, \ldots,20\ 969$.

Figure 25(b) shows that the flow periodically travels through the ten subclusters in the child layer $\mathcal {L}_{2}$ of cluster $\mathcal {C}_{2}$. As for the centroids of figure 25, the alleys of vortices do not cross the $x$-axis. The closed transition $\mathcal {C}_1 \rightarrow \cdots \rightarrow \mathcal {C}_{10} \rightarrow \mathcal {C}_1$ is further evidenced in the transition matrix of figure 25(c). The residence times of figure 25(d) are rather uniform and the averaged period is 6.62.

Figure 25. Cluster-based analysis in layer 2 at $Re= 130$ for $\mathcal {C}_{2}$: (a) graph of non-trivial transitions between clusters, as in figure 9, (b) transition illustrated with cluster label, (c) transition matrix, (d) residence time matrix, as in figure 10. Two trajectories pass through $\mathcal {C}_{2}$ in the parent layer, one with $m_2=1, \ldots,1088$ and another with $m_2=1089, \ldots,2175$.

The chaotic state in the group of clusters $\mathcal {C}_{4},\ldots,\mathcal {C}_{7}$

When dealing with the chaotic dynamics, considering each cluster $\mathcal {C}_{4},\ldots,\mathcal {C}_{7}$ separately is possible. However, instead of building a network model for each cluster separately, we can build an overall model for these four clusters in the chaotic regime. Therefore, for analysing the chaotic dynamics, we consider all the clusters that belong to the chaotic cloud in layer $\mathcal {L}_{1}$ as a whole.

From figure 22(a), all four trajectories will reach the same chaotic attractor described by the chaotic clusters, with the original snapshot index $m=6281, \ldots,15\ 000, 21\ 301, \ldots,$ $30\ 000, 33\ 365, \ldots, 40\ 000$ and $43\ 068, \ldots,50\ 000$. In figure 26, the closed orbit of the clusters $\mathcal {C}_1 \rightarrow \cdots \rightarrow \mathcal {C}_{6} \rightarrow \mathcal {C}_1$ is formed, with a relatively high probability of transition between the successive clusters. The flow field of the centroids of the first six clusters form a complete cycle of vortex shedding. This is interesting, as it recalls the periodic and quasi-periodic dynamics respectively observed at $Re=80$ and $Re=105$. The periodic block of the transition matrix of figure 27(b), from $\mathcal {C}_1$ to $\mathcal {C}_6$, corroborates the existence of the periodic dynamics. The centroids of the clusters in the cycle present a similar structure of coherence, as can be seen in figure 26. There are also other possible transitions from clusters on this orbit to the remaining clusters with much smaller probability. The residence times shown in figure 27(c) for each cluster on the orbit are uniform, and the averaged period along the complete cycle is $7.78$ by summing up the elements in the block from $\mathcal {C}_1$ to $\mathcal {C}_6$.

Figure 26. Graph of transitions between clusters in layer 2 at $Re= 130$ for the chaotic clusters of $\mathcal {L}_{1}$, displayed as in figure 9.

Figure 27. Cluster-based analysis in layer 2 at $Re= 130$ for the chaotic clusters of $\mathcal {L}_{1}$: (a) transition illustrated with cluster label, (b) transition matrix, (c) residence time matrix, as in figure 10. All four trajectories reach the chaotic clusters, with $m_2=1, \ldots,8720, 8721, \ldots,17\ 420, 17\ 421, \ldots,24\ 056$ and $24\ 057, \ldots,30\ 989$.

The remaining clusters $\mathcal {C}_{7}, \ldots, \mathcal {C}_{10}$ have multiple destinations with quasi-random possibilities. Even though the probabilities of these random transitions are small, they contribute to the chaotic dynamics of the flow field, with recurrent transitions $\mathcal {C}_1 \rightarrow \mathcal {C}_{10}, \mathcal {C}_{9} \rightarrow \mathcal {C}_1$, and so on. The flow fields of the associated centroids are shown in figure 26. Their structure looks like distortions of the vortex shedding cycle formed by the first six clusters. The network model indicates that the fully chaotic state still contains a main cycle of clusters associated with a periodic vortex shedding, together with the random jumping to the clusters associated with a stochastic disorder in the wake. In this case, the transition matrix can be used to build a stochastic model, as shown in Appendix D.

4.5.3. Dynamics reconstruction of the hierarchical network model at $Re=130$

Figure 28 shows the autocorrelation function of the DNS and the HiCNM in the two layers. We stop the hierarchical modelling in layer $\mathcal {L}_{2}$, as both the transient and chaotic dynamics can be fairly reproduced with a limited number of clusters. The autocorrelation function of the DNS shows chaotic oscillations with a dominant frequency. In layer $\mathcal {L}_{1}$, no oscillation can be identified, and the RMSE of the autocorrelation function is $R_{rms}^1 = 6.94$. In layer $\mathcal {L}_{2}, R_{rms}^2 = 3.42$. The autocorrelation function of the model matches well with the dominant frequency of the oscillations, but the $\hat {R}_{\mathcal {L}_{2}}(\tau )$ value can hardly match.

Figure 28. Autocorrelation function for $\tau \in [0, 1500)$ from DNS (black solid line) and the hierarchical network model (red dashed line) in the two layers: (a) $\mathcal {L}_{1}$ and (b) $\mathcal {L}_{2}$, at $Re=130$.