2.1. Clustering the state space

The dynamics-augmented clustering procedure is divided into two steps. Initially, the state space is clustered, yielding coarse-grained state transition dynamics with trajectory segments composed of time-continuous snapshots within each cluster. Subsequently, we cluster these trajectory segments, utilising centroids derived from the average of each segment. This step optimises the centroid distribution and eliminates the redundancy of the trajectory segments.

The first-stage clustering discretizes the high-dimensional state space by grouping the snapshots. We first define a Hilbert space $\mathscr {L}^{2}(\varOmega )$, in which the inner product of vector fields in the domain $\varOmega$ is given by a square-integrable function:

(2.1)\begin{equation} (\boldsymbol{u}, \boldsymbol{v})_{\varOmega} = \int_{\varOmega} \,\mathrm{d}\kern0.7pt \boldsymbol{x} \, \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{v}. \end{equation}

Here $\boldsymbol {u}$ and $\boldsymbol {v}$ represent snapshots of this vector field, also known as observations in the machine learning context. The corresponding norm is defined as

(2.2)\begin{equation} \|\boldsymbol{u}\|_{\varOmega}:=\sqrt{(\boldsymbol{u}, \boldsymbol{u})_{\varOmega}}. \end{equation}

The distance $D$ between two snapshots can be calculated as

(2.3)\begin{equation} D(\boldsymbol{u},\boldsymbol{v}) = \| \boldsymbol{u} - \boldsymbol{v} \|_{\varOmega}. \end{equation}

The unsupervised $k$-means++ algorithm (MacQueen Reference MacQueen1967; Lloyd Reference Lloyd1982; Arthur & Vassilvitskii Reference Arthur and Vassilvitskii2007) is used for clustering. It operates automatically, devoid of assumptions or data categorisation. Serving as the foundation of cluster analysis, this algorithm partitions a set of $M$ time-resolved snapshots $\boldsymbol {u}^m$, where $m=1 \ldots M$, into $K$ clusters $\mathcal {C}_k$, where $k = 1 \dots K$. Each cluster corresponds to a centroidal Voronoi cell, with the centroid defined as the average of the snapshots within the same cluster. The algorithm comprises the following steps.

1. (i) Initialisation: $K$ centroids $\boldsymbol {c}_{k}$, where $k = 1 \ldots K$, are randomly selected. In contrast to the $k$-means algorithm, $k$-means++ optimises the placement of these centroids to prevent sensitivity to initial conditions.

2. (ii) Assignment: each snapshot $\boldsymbol {u}^{m}$ is allocated to the nearest centroid by $\underset {k}{\rm arg\min } D(\boldsymbol {u}^{m},\boldsymbol {c}_{k})$. The characteristic function is used to mark their affiliation, and it is defined as follows:

   (2.4)\begin{equation} \chi_{k}^{m}:=\left\{\begin{array}{@{}ll@{}} 1, & \text{if } \boldsymbol{u}^{m} \in \mathcal{C}_{k}, \\ 0, & \text{otherwise}. \end{array}\right. \end{equation}

3. (iii) Update: each centroid is recalculated by averaging all the snapshots belonging to the corresponding cluster as

   (2.5)\begin{equation} \boldsymbol{c}_{k}=\frac{1}{M_{k}} \sum_{\boldsymbol{u}^{m} \in \mathcal{C}_{k}} \boldsymbol{u}^{m}=\frac{1}{M_{k}} \sum_{m=1}^{M} \chi_{k}^{m} \boldsymbol{u}^{m} , \end{equation}
   where
   (2.6)\begin{equation} M_{k}=\sum_{m=1}^{M} \chi_{k}^{m}. \end{equation}

4. (iv) Iteration: the assignment and update steps are repeated until convergence is reached. Convergence means that the centroids do not move or stabilise below a certain threshold. The algorithm minimises the intra-cluster variance and maximises the inter-cluster variance. The intra-cluster variance is computed as follows:

   (2.7)\begin{equation} J\left(\boldsymbol{c}_{1}, \ldots, \boldsymbol{c}_{K}\right) = \sum_{k=1}^{K} \sum_{m=1}^{M} \chi_{k}^{m} \| \boldsymbol{u}^{m}-\boldsymbol{c}_{k} \|_{\varOmega}^{2} . \end{equation}
   Each iteration reduces the value of the criterion $J$ until convergence is reached.

The cluster probability distribution $\boldsymbol {P} = [P_1, \ldots, P_K]$ is determined by $P_k = M_{k}/M$ for each cluster $\mathcal {C}_k$, and satisfies the normalisation condition ${\sum }_{k=1}^{K}P_k = 1$.

The geometric properties of the clusters are quantified for further analysis. The cluster standard deviation on the snapshots $R_{k}^{u}$ measures the cluster size, following Kaiser et al. (Reference Kaiser2014), as

(2.8)\begin{equation} R_{k}^{u} = \sqrt{ \frac{1}{M_{k}} \sum_{m=1}^{M} \chi_{k}^{m} \left\| \boldsymbol{u}^{m}-\boldsymbol{c}_{k} \right\|_{\varOmega}^{2}} . \end{equation}

The time-resolved snapshots should be equidistantly sampled and cover a statistically representative time window of the coherent structure evolution. As a rule of thumb, at least ten periods of the dominant frequency are needed to obtain reasonably accurate statistical moments and at least $K$ snapshots per characteristic period to capture an accurate temporal evolution.

2.2. Clustering the trajectory segments

After the state space is discretized, the trajectory is also divided into segments. We use the cluster transition information to identify the trajectory segments that pass through a cluster.

Based on the temporal information from the given data set, the nonlinear dynamics between snapshots are modelled as linear transitions between clusters, known as the classic CNM (Fernex et al. Reference Fernex, Noack and Semaan2021; Li et al. Reference Li, Fernex, Semaan, Tan, Morzyński and Noack2021). We infer the probability of cluster transition from the data as

(2.9)\begin{equation} Q_{ik}=\frac{n_{ik}}{n_{k}}, \quad i, k=1, \ldots, K, \end{equation}

where $Q_{ik}$ is the direct cluster transition probability from cluster $\mathcal {C}_{k}$ to $\mathcal {C}_{i}$ and $n_{ik}$ is the number of transitions from $\mathcal {C}_{k}$ only to $\mathcal {C}_{i}$:

(2.10)\begin{equation} n_{ik}=\sum_{m=1}^{M} \chi_{ik}^m. \end{equation}

Here

(2.11)\begin{equation} \chi_{ik}^m=\left\{\begin{array}{@{}ll@{}} 1, & \text{if } \boldsymbol{u}^{m} \in \mathcal{C}_{k}\ \text{and} \ \boldsymbol{u}^{m+1} \in \mathcal{C}_{i},\\ 0, & \text{otherwise}, \end{array}\right. \end{equation}

$n_{k}$ is the total number of transitions from $\mathcal {C}_{k}$ regardless of the destination cluster:

(2.12)\begin{equation} n_{k}=\sum_{i=1}^{K} n_{ik}, \quad i, k=1, \ldots, K. \end{equation}

If $Q_{ik} \neq 0$, it can be inferred that in cluster $\mathcal {C}_{k}$ there exists at least one trajectory segment that is bound for cluster $\mathcal {C}_{i}$. We assign distinct labels to each trajectory segment corresponding to all destination clusters, denoted as $\mathcal {T}_{(kl)}$, where $k$ and $l$ represent the $l$th segment in $\mathcal {C}_{k}$. Therefore, the snapshots are marked according to their trajectory affiliations by a characteristic function:

(2.13)\begin{equation} \chi_{(kl)}^{m} = \left\{\begin{array}{@{}ll@{}} 1, & \text{if } \boldsymbol{u}^{m} \in \mathcal{T}_{(kl)} ,\\ 0, & \text{otherwise}. \end{array}\right. \end{equation}

Here $k$ represents the cluster affiliation and $l$ represents the trajectory segment affiliation. The total number of trajectory segments in $\mathcal {C}_{k}$ equals $n_{k}$. Note that the final trajectory segment of the data set will not be considered as it will not lead to any destination cluster and is usually incomplete. The total number of trajectory segments in the data set can be obtained by the sum of $n_{k}$ as follows:

(2.14)\begin{equation} n_{traj} = \sum_{k=1}^{K} n_{k}. \end{equation}

Analogous trajectory segments within the same cluster will be merged in the subsequent clustering stage. Operations on the trajectories can often be costly. Efficiency in clustering can be achieved by mapping the operations performed on trajectory segments to their corresponding averages, i.e. the trajectory segment centroids, given their topological relationship. Additionally, the propagation of our model relies on centroids, rendering the trajectory information essentially unnecessary. We define the centroids $\boldsymbol {c}_{(kl)}$ as the average of snapshots belonging to the same trajectory segment:

(2.15)\begin{equation} \boldsymbol{c}_{(kl)}=\frac{1}{M_{(kl)}} \sum_{\boldsymbol{u}^{m} \in \mathcal{T}_{(kl)}}\boldsymbol{u}^{m}=\frac{1}{M_{(kl)}} \sum_{m=1}^{M} \chi_{(kl)}^{m} \boldsymbol{u}^{m}. \end{equation}

Here

(2.16)\begin{equation} M_{(kl)}= \sum_{m=1}^{M} \chi_{(kl)}^{m}. \end{equation}

The subsequent question pertains to how to determine the number of sub-clusters. The allocation of sub-clusters within each cluster can be automatically learnt from the data. To maintain the spatial resolution, more sub-clusters should be assigned to clusters with a larger transverse size. We first introduce a transverse cluster size vector $\boldsymbol {R}^{\mathcal {T}}$, which is defined by the standard deviation of the $n_{k}$ centroids $\boldsymbol {c}_{(kl)}$ with respect to the cluster centroid $\boldsymbol {c}_{k}$ as follows:

(2.17)\begin{equation} R_{k}^{\mathcal{T}} = \sqrt{\frac{1}{n_{k}} \sum_{l=1}^{n_{k}}\|\boldsymbol{c}_{(kl)}-\boldsymbol{c}_{k}\|_{\varOmega}^{2}}. \end{equation}

Next, we denote the number of sub-clusters as $L_{k}$ for clustering the centroids in cluster $\mathcal {C}_{k}$. A $K$-dimensional vector $\boldsymbol {L} = [L_{1}, \ldots, L_{K}]^{\intercal }$ records the numbers of sub-clusters in each cluster, with $L_{k}$ determined by

(2.18)\begin{equation} L_{k} = \min (\lfloor \hat{R}^{\mathcal{T}}_{k} n_{traj} (1 - \beta) \rfloor +1, n_{k}). \end{equation}

Here, the vector $\boldsymbol {R}^{\mathcal {T}}$ is normalised with the sum ${\sum }_{k=1}^{K} R_{k}^{\mathcal {T}}$, denoted as $\boldsymbol {\hat {R}}^{\mathcal {T}}$, which ensures a suitable distribution of sub-clusters for the ensemble of $n_{traj}$ trajectories. To increase the flexibility of the model, we introduce a sparsification controller $\beta \in [0, 1]$ in this clustering process. For the extreme value of $\beta =1$, all the centroids are merged into one centroid, and the dCNM is identical to a classic CNM, with the maximum sparsification. For the other extreme $\beta =0$, the dCNM is minimally sparsified according to the transverse cluster size. For periodic or quasi-periodic systems, the dCNM with a large $\beta$ can capture most of the dynamics, while for complex systems such as chaotic systems, a small $\beta$ may be needed. In addition, the minimum function prevents the possibility that the left-hand side of the equation exceeds the number of centroids $n_{k}$ when $\beta$ is too small, causing the second-stage clustering to not be performed. The choice of $\beta$ is discussed in Appendix C.

The refined centroids are obtained by averaging a series of centroids related to analogous trajectory segments. The redundancy of the $n_{traj}$ centroids is mitigated, and the corresponding transition network becomes sparse. The $k$-means++ algorithm is also used in the second-stage clustering. It will iteratively update the centroids $\boldsymbol {c}_{(kl)}$ and the characteristic function $\chi _{(kl)}^{m}$ until convergence or the maximum number of iterations is reached. The overall clustering process of the dCNM is summarised in Algorithm 1.

Algorithm 1: Pseudocode for the dynamics-augmented clustering procedure

2.3. Characterising the transition dynamics

We use the centroids obtained from § 2.2 as the nodes of the network and the linear transitions between these centroids as the edges of the network. First, we introduce two transition properties: the centroid transition probability $Q_{(ij) (kl)}$ and the transition time $T_{ik}$.

Figure 2 illustrates the definition of the subscripts in the centroid transition probability $Q_{(ij) (kl)}$, which can contain all possible transitions between the refined centroids of clusters $\mathcal {C}_{k}$ and $\mathcal {C}_{i}$. Considering the transitions between these centroids, we define $Q_{(ij) (kl)}$ as

(2.19)\begin{equation} Q_{(ij) (kl)} = \frac{n_{(ij) (kl)}}{n_{k}}, \quad i, k =1, \ldots, K, \ j = 1, \ldots, L_{i}, \ l = 1, \ldots, L_{k}, \end{equation}

where $n_{(ij) (kl)}$ is the number of transitions from $\boldsymbol {c}_{(kl)}$ only to $\boldsymbol {c}_{(ij)}$. This definition differs from that of the CNM, which uses the cluster transition $Q_{ik}$ in (2.9) to define the probability. In fact, we can compute $Q_{ik}$ by summing up $Q_{(ij) (kl)}$ as follows:

(2.20)\begin{equation} Q_{ik}=\sum_{j=1}^{L_{i}} \sum_{l=1}^{L_{k}} Q_{(ij) (kl)}. \end{equation}

The definition of the transition time $T_{ik}$ is identical to the CNM, as shown in figure 3. This property is not further investigated in the present work, as the transition time crossing the same clusters varies little in most dynamic systems.

Figure 2. Illustration of the subscripts in the refined centroid transitions. After the state space is clustered, only one subscript is needed to distinguish the different clusters, such as $\mathcal {C}_{k}$ and $\mathcal {C}_{i}$. After the trajectory segments are clustered, two subscripts are needed to represent the refined centroids, such as $\boldsymbol {c}_{(kl)}$ in $\mathcal {C}_{k}$ and $\boldsymbol {c}_{(ij)}$ in $\mathcal {C}_{i}$.

Figure 3. Individual transition time $\tau ^{n}_{ik}$ for the transition from cluster $\mathcal {C}_{k}$ to $\mathcal {C}_{i}$.

Let $t^{n}$ be the instant when the first snapshot enters, and $t^{n+1}$ be the instant when the last snapshot leaves on one trajectory segment passing through cluster $\mathcal {C}_{k}$. The residence time $\tau _{k}^{n}$ is the duration of staying in cluster $\mathcal {C}_{k}$ on this segment, which is given by

(2.21)\begin{equation} \tau_{k}^{n}=t^{n+1}-t^{n}. \end{equation}

For an individual transition from $\mathcal {C}_{k}$ to $\mathcal {C}_{i}$, the transition time is defined as $\tau _{ik}^{n}$, which can be obtained by the average of the residence times from both clusters:

(2.22)\begin{equation} \tau_{ik}^{n}=(\tau_{k}^n + \tau_{i}^n)/2. \end{equation}

By averaging $\tau _{ik}^{n}$ for all the individual transitions from $\mathcal {C}_{k}$ to $\mathcal {C}_{i}$, the transition time can be expressed as

(2.23)\begin{equation} T_{ik}=\frac{\sum_{n=1}^{n_{ik}} \tau_{ik}^{n}}{n_{ik}}. \end{equation}

The essential dynamics can also be summarised into single entities as in the CNM, since the cluster-level information is still retained in the current model. For completeness, we introduce the cluster transition probability matrix $\boldsymbol {Q}$ and the cluster transition time matrix $\boldsymbol {T}$ as

(2.24)\begin{equation} \left. \begin{aligned} \boldsymbol{Q} & =Q_{ik} \in \mathbb{R}^{K \times K}, \quad i, k=1, \ldots, K\\ \boldsymbol{T} & =T_{ik} \in \mathbb{R}^{K \times K}, \quad i, k=1, \ldots, K. \end{aligned} \right\} \end{equation}

The cluster indices are reordered in both matrices to enhance readability. Cluster $\mathcal {C}_1$ is the cluster with the highest distribution probability, $\mathcal {C}_2$ is the cluster with the highest transition probability leaving from $\mathcal {C}_1$ and $\mathcal {C}_3$ is the cluster with the highest transition probability leaving from $\mathcal {C}_2$, so on and so forth. If the cluster with the highest probability is already assigned, we choose the cluster with the second highest probability. If all the clusters with non-zero transition probabilities are already assigned, we choose the next cluster with the highest distribution probability among the rest.

By analogy with $\boldsymbol {Q}$, the centroid transition probability $Q_{(ij) (kl)}$ for given affiliations of the departure cluster $k$ and destination cluster $i$ can form a centroid transition matrix $\boldsymbol {Q}_{ik}$ that captures all possible centroid dynamics between the two clusters:

(2.25)\begin{equation} \boldsymbol{Q}_{ik} = Q_{(ij) (kl)} \in \mathbb{R}^{L_{i} \times L_{k}}, \quad j = 1, \ldots, L_{i}, \ l = 1, \ldots, L_{k}. \end{equation}

Moreover, to summarise the centroid transition dynamics, the centroid transition probability $Q_{(ij) (kl)}$ for a given affiliation $k$ of only the departure cluster can form a centroid transition tensor $\mathcal {Q}_k$ that captures all the possible centroid dynamics from this cluster, as

(2.26)\begin{equation} \mathcal{Q}_{k} = Q_{(ij) (kl)} \in \mathbb{R}^{K \times L_{i} \times L_{k}}, \quad i =1, \ldots, K, \ j = 1, \ldots, L_{i}, \ l = 1, \ldots, L_{k}. \end{equation}

The dCNM propagates the state motion based on the centroids $\boldsymbol {c}_{(kl)}$ for the reconstruction. To determine the transition dynamics, we first use $\mathcal {Q}_{k}$ to find the centroid transitions from the initial centroid $\boldsymbol {c}_{(kl)}$ to the destination centroid $\boldsymbol {c}_{(ij)}$. As the destination centroids are determined, the cluster-level dynamics is determined correspondingly. Then, $\boldsymbol {T}$ is used to identify the related transition time.

We assume a linear state propagation between the two centroids $\boldsymbol {c}_{(kl)}$ and $\boldsymbol {c}_{(ij)}$ obtained from the tensors, as follows:

(2.27)\begin{equation} \boldsymbol{u}^m(t) = \alpha_{ik}(t) \boldsymbol{c}_{(ij)} + \left[ 1-\alpha_{ik}(t) \right] \boldsymbol{c}_{(kl)}, \quad \alpha_{ik}=\frac{t-t_{k}}{T_{ik}}. \end{equation}

Here $t_{k}$ is the time when the centroid $\boldsymbol {c}_{(kl)}$ is left. Note that we can use splines (Fernex et al. Reference Fernex, Noack and Semaan2021) or add the trajectory supporting points (Hou, Deng & Noack Reference Hou, Deng and Noack2022) to interpolate the motion between the centroids for smoother trajectories.

Intriguingly, we observe that the trajectory-based clustering of the dCNM enhances the resolution of the cluster transitions. Now each centroid only has a limited number of destination centroids, often within the same cluster. This minimises the likelihood of selecting the wrong destination cluster based solely on the cluster transition probability matrix, as is the case in classic CNM. Consequently, it becomes feasible to accurately resolve long-term cluster transitions without the need for historical information. It can be argued that dCNM effectively constrains cluster transitions, leading to outcomes similar to those obtained with the higher-order CNM (Fernex et al. Reference Fernex, Noack and Semaan2021). This improvement is attained by replacing higher-order indexing with higher-dimensionality dual indexing. Specifically, the dual indexing also results in a substantial reduction in the model complexity. While the complexity of the high-order CNM is defined as $K^{\tilde {L}}$, where $K$ is the number of clusters and $\tilde {L}$ is the order, the model complexity of the dCNM is expressed as ${\sum }_{{k}=1}^{K} {L}_{k}$, which is a significantly lower value, particularly when $\tilde {L}$ is relatively large. In terms of computational efficiency, dCNM with $\beta = 0.80$ reduces the computational time by $40\,\%$ as compared with CNM with the same number of centroids. This improvement is primarily attributed to the hierarchical clustering. The computational load of the first-stage clustering on the state space is reduced by a small number of clusters $K$. The second-stage clustering on the trajectory segments accounts only for $20\,\%$ of the total computation time.

2.4. Validation

The auto-correlation function and the representation error are used for validation. We examine the prediction errors for cluster-based models considering both spatial and temporal perspectives. The spatial error arises from the inadequate representation by cluster centroids, as evidenced by the representation error and the auto-correlation function. The temporal error arises due to the imprecise reconstruction of intricate snapshot transition dynamics. This can be observed directly through the temporal evolution of snapshot affiliations and, to some extent, through the auto-correlation function.

The auto-correlation function is a practical tool for evaluating ROMs, as it can statistically reflect the prediction errors. Additionally, the auto-correlation function circumvents the problem of directly comparing two trajectories with finite prediction horizons, which may suffer from phase mismatch (Fernex et al. Reference Fernex, Noack and Semaan2021). This is particularly relevant for chaotic dynamics, whereby minor differences in initial conditions can lead to divergent trajectories, making the direct comparison of time series meaningless. The unbiased auto-correlation function of the state vector (Protas, Noack & Östh Reference Protas, Noack and Östh2015) is given by

(2.28)\begin{equation} R(\tau )=\frac{1}{T-\tau } \int_{0}^{T-\tau} (\boldsymbol{u}(x,t) , \boldsymbol{u}(x,t+\tau ))_\varOmega\,\mathrm{d}t, \quad \tau\in \left [ 0,T \right ]. \end{equation}

In this study, $R(\tau )$ will be normalised by $R(0)$ (Deng et al. Reference Deng, Noack, Morzyński and Pastur2022). This function can also infer the spectral behaviour by computing the fluctuation energy at the vanishing delay.

The representation error can be numerically computed as

(2.29)\begin{equation} E_{r}=\frac{1}{M} \sum_{m=1}^{M} D_{\mathcal{T}}^{m}, \end{equation}

where $D_{\mathcal {T}}^{m}$ is the minimum distance from the snapshot $\boldsymbol {u}^m$ to the states on the reconstructed trajectory $\mathcal {T}$:

(2.30)\begin{equation} D_{\mathcal{T}}^{m}=\min _{\boldsymbol{u}^{n} \in \mathcal{T}}\|\boldsymbol{u}^{m}-\boldsymbol{u}^{n}\|_{\varOmega}. \end{equation}

4.1. Numerical methods and flow features

Numerical simulation is performed to obtain the data set, as shown in figure 8. A sphere with a diameter $D$ is placed in a uniform flow with a streamwise velocity $U_\infty$. The computational domain takes the form of a cylindrical tube, with its origin at the centre of the sphere and its axial direction along the streamwise direction ($x$ axis). The dimensions of the domain in the $x$, $y$ and $z$ directions are $80D$, $10D$ and $10D$, respectively. The inlet is located $20D$ upstream from the sphere. These specific domain parameters are chosen to minimise any potential distortion arising from the outer boundary conditions while also mitigating computational costs (Pan, Zhang & Ni Reference Pan, Zhang and Ni2018; Lorite-Díez & Jiménez-González Reference Lorite-Díez and Jiménez-González2020). The fluid flow is governed by the incompressible Navier–Stokes equations:

(4.1)\begin{equation} \left. \begin{gathered} \partial{\boldsymbol{u}} / \partial{t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} p- \nabla^{2} \boldsymbol{u} / {{Re}} = 0,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0. \end{gathered} \right\} \end{equation}

Here $\boldsymbol {u}$ denotes the velocity vector $(u_x, u_y, u_z)$, $p$ is the static pressure and ${{Re}}$ is the Reynolds number, which is defined by

(4.2)\begin{equation} {{Re}} = U_{\infty} D / \nu, \end{equation}

where $\nu$ is the kinematic viscosity.

Figure 8. Numerical sketch of the sphere wake.

The net forces on the sphere have three components $F_\alpha$, $\alpha = x, y, z$, and the corresponding force coefficients $C_{\alpha }$ are defined as

(4.3)\begin{equation} C_{\alpha}=\frac{2 F_{\alpha}}{\rho U_\infty^2 S}, \end{equation}

where $S={\rm \pi} D^{2}/4$ is the projected surface area of the sphere in the streamwise direction. The total drag force coefficient is $C_D = C_x$. Since the lift coefficient can have any direction in the $yz$ plane on the axisymmetric sphere, the total lift force coefficient $C_L$ is given by

(4.4)\begin{equation} C_L=\sqrt{C_{y}^{2}+C_{z}^{2}}. \end{equation}

The flow parameters are non-dimensionalised based on the characteristic length $D$ and the free-stream velocity $U_\infty$. This implies that the time unit scales are $D/U_\infty$ and the pressure scales are $\rho U_\infty ^2$, where $\rho$ is the density. The Strouhal number $St$ is correspondingly expressed as

(4.5)\begin{equation} St = f, \end{equation}

where $f$ is the characteristic frequency.

The ANSYS Fluent $15.0$ software is used as the computational fluid dynamics (CFD) solver for the governing equations with the cell-centred finite volume method. We impose a uniform streamwise velocity $\boldsymbol {u} = [U_\infty, 0, 0]$ at the inlet boundary and an outflow condition at the outlet boundary. The outflow condition is set as a Neumann condition for the velocity, $\partial _x \boldsymbol {u} = [0, 0, 0]$, and a Dirichlet condition for the pressure, $p_{out} = 0$. We apply a no-slip boundary condition on the sphere surface and a slip boundary condition on the cylindrical tube walls to prevent wake-wall interpolations. The pressure-implicit split-operator algorithm is chosen for pressure–velocity coupling. For the governing equations, the second-order scheme is used for the spatial discretization, and the first-order implicit scheme is used for the temporal term. To satisfy the Courant–Friedrichs–Levy condition, a small integration time step is set as $\Delta t = 0.01$ non-dimensional time unit, such that the Courant number is below $1$ for all simulations. For the periodic flow at ${{Re}} = 300$, the simulation starts in the vicinity of the steady solution and runs for $t = 200$ time units, incorporating the transient and post-transient dynamics. For the quasi-periodic flow, the simulations are performed for $t = 500$ time units and for the chaotic flow for $t = 700$ time units. The snapshots are collected at a sampling time step of $\Delta t_s = 0.2$ time units for all the test cases. Moreover, we discard the first $200$ time units to eliminate any transient phases for the quasi-periodic and chaotic cases. The relevant numerical investigation approach can be found in Johnson & Patel (Reference Johnson and Patel1999) and Rajamuni, Thompson & Hourigan (Reference Rajamuni, Thompson and Hourigan2018). For the convergence and validation studies, see Appendix A.

The wake of a sphere exhibits different flow regimes as ${{Re}}$ increases, ultimately transitioning to a chaotic state. At ${{Re}} = 20 \sim 24$, flow separation occurs, forming a steady recirculating bubble, as observed in previous studies (Sheard, Thompson & Hourigan Reference Sheard, Thompson and Hourigan2003; Eshbal et al. Reference Eshbal, Rinsky, David, Greenblatt and van Hout2019). The length of this wake grows linearly with $\ln ({{Re}})$. When ${{Re}}$ surpasses $130$ (Taneda Reference Taneda1956), the wake bubble starts oscillating in a wave-like manner, while the flow maintains axisymmetry. The first Hopf bifurcation takes place at approximately ${{Re}} \approx 212$ (Fabre et al. Reference Fabre, Auguste and Magnaudet2008), leading to a loss of axisymmetry and the emergence of a planar-symmetric double-thread wake with two stable and symmetric vortices. The orientation of the symmetry plane can vary (Johnson & Patel Reference Johnson and Patel1999). At a subsequent Hopf bifurcation around ${{Re}} = 270 \sim 272$ (Johnson & Patel Reference Johnson and Patel1999; Fabre et al. Reference Fabre, Auguste and Magnaudet2008), the flow becomes time dependent, initiating periodic vortex shedding with the same symmetry plane as before. In the range $272 < {{Re}} < 420$ (Eshbal et al. Reference Eshbal, Rinsky, David, Greenblatt and van Hout2019), periodicity and the symmetry plane diminish, with the vortex shedding becoming quasi-periodic and then fully three dimensional. Beyond ${{Re}} = 420$, shedding becomes irregular and chaotic (Ormières & Provansal Reference Ormières and Provansal1999; Pan et al. Reference Pan, Zhang and Ni2018; Eshbal et al. Reference Eshbal, Rinsky, David, Greenblatt and van Hout2019), due to the azimuthal rotation of the separation point and lateral oscillations of the shedding.

In this study we examine three baseline flow regimes of the sphere wake: periodic flow at ${{Re}} = 300$, quasi-periodic flow at ${{Re}} = 330$ and chaotic flow at ${{Re}} = 450$. Figure 9 illustrates the flow characteristics of these regimes. Figure 9(a,b) depicts the lift coefficient $C_L$ of the periodic flow. The oscillation of $C_L$ starts after a few time units, with the amplitude gradually increasing in the cycle-to-cycle evolution, eventually forming the limit cycle. Figure 9(c) shows the instantaneous vortex structures of the quasi-periodic flow, identified by the $Q$ criterion and colour coded by the non-dimensional velocity $U_{\infty }$. The vortex shedding forms hairpin vortices with slight variations between successive shedding events, signifying the absence of short-term periodicity while retaining long-term periodic behaviour. Figure 9(d,e) shows the temporal evolution of $C_L$ and its phase portrait, respectively. The amplitude of $C_L$ is strongly associated with the quasi-periodic dynamic, and the modulation also thickens the limit cycle of the oscillator on the phase portrait. The power spectral density in figure 9(f) indicates two dominant frequencies: a higher frequency linked to natural shedding and a lower frequency associated with amplitude modulation resulting from variations between shedding events. For chaotic flow, periodicity entirely vanishes, and the flow regime displays the typical features of a chaotic system. The hairpin vortexes in figure 9(g) shed irregularly, with varying separation angles and even double spirals. The temporal evolution of $C_L$ in figure 9(h) exhibits more complex dynamics, with the phase diagram in figure 9(i) depicting many random loops that no longer exhibit circular patterns. Furthermore, the power spectral density of $C_L$ in figure 9(j) shows a broad peak, also indicating chaotic features.

Figure 9. Flow characteristics of the sphere wake. The periodic flow at ${{Re}}=300$ including the transient and post-transient dynamics is displayed by (a) the phase portrait of the lift coefficient $C_L$ and (b) the temporal evolution of $C_L$. The quasi-periodic flow at ${{Re}}=330$ is displayed by (c) vortex structures, where the vortexes are identified by the $Q$ criteria and are colour coded by the non-dimensional velocity $U_{\infty }$, (d) temporal evolution of $C_L$, (e) phase portrait of $C_L$ and (f) power spectral density of $C_L$ on a time series of length $T_{traj} = 100$ (red curve) and $T_{traj} = 300$ (black curve). The chaotic flow at ${{Re}}=450$ is displayed by (g) vortex structures, (h) temporal evolution of $C_L$, (i) phase portrait of $C_L$ and (j) power spectral density of $C_L$ on a time series of length $T_{traj} = 500$.

We performed a lossless POD preprocessing on the snapshots to reduce the computational cost of clustering the three-dimensional flow field data set, as described in Appendix B. This preprocessing is optional and does not affect the distance measure in the clustering algorithm. For consistency, the notation snapshot is maintained in the following for the preprocessed data.

4.2. The periodic flow regime at ${{Re}} = 300$

We compare the dCNM to the CNM for the periodic flow regime of the sphere wake at ${{Re}} = 300$. The transient and post-transient dynamics are considered, providing insights into the mechanisms for the instability and nonlinear saturation.

The comparison between the CNM and the dCNM clustering for the periodic flow is presented in figure 10. In the dCNM we set $K=10$ for state space clustering and $\beta = 0.50$ for the sub-clustering. The value of $\beta$ characterises the trade-off between a small number of sub-clusters and the model accuracy. The normalised transverse cluster size vector $\boldsymbol {\hat {R}}^{\mathcal {T}} = [0.1510,0.1101,0.0855,0.1120,0.1069,0.1031,0.0650,0.0934,0.0883,0.0847]^{\intercal }$ corresponds to the number of sub-clusters $\boldsymbol {L} = [3,3,3,4,4,4,3,3,3,3]^{\intercal }$. Classical multidimensional scaling is applied to project the high-dimensional snapshots and centroids into a three-dimensional subspace $[a_1, a_2, a_3]^\intercal$ for visualisation. The snapshots form a conical surface in the three-dimensional subspace, where the trajectory spirals up from a fixed point to a periodic motion. This behaviour is indicative of a Hopf bifurcation, which involves an unstable steady solution and nonlinear saturation to a periodic limit cycle. Most of the CNM centroids are located on the limit cycle and only a few resolve the transient phase. In contrast, the dCNM centroids offer a finer resolution of the amplitude growth.

Figure 10. Three-dimensional visualisation of the clustered periodic flow regime of the sphere wake at ${{Re}} = 300$. Classical multidimensional scaling is applied to the data set to visualise the high-dimensional snapshots and centroids in the subspace. The small dots represent the snapshots and the large dots represent the centroids. Snapshots and centroids with the same colour belong to the same cluster. For comparison, the CNM result in (a) is shown with the same number of centroids as the corresponding dCNM result. The dCNM result in (b) is shown with $K = 10$ and $\beta = 0.50$.

The original and reconstructed trajectories of the CNM and dCNM for the periodic flow regime are shown in figure 11. The CNM fails to resolve the transient dynamics and only captures the stable limit cycle. In contrast, the dCNM can effectively resolve both the cyclic behaviour and the growing oscillation amplitude. Given the deterministic nature of these transient and post-transient dynamics, the centroid transition should yield a permutation matrix – each column and each row of the matrix should only have one element being unity. However, in the CNM, under-resolved transient dynamics result in a stochastic transition matrix. The transition uncertainty expressed within a column of the matrix often leads to prediction errors. Examples are the unphysical jumps between erratic cycles in the transient stage, as displayed in figure 10(a). The sub-clusters in dCNM are designed to capture a one-way forward transition from the starting point to the limit cycle, ensuring an accurate reconstruction of the deterministic dynamics.

Figure 11. Trajectory of the periodic flow at ${{Re}} = 300$. The thin grey curve represents the original trajectory, the thick red curve represents the reconstructed trajectory and the red dots represent the centroids. (a) The CNM reconstruction and (b) the dCNM reconstruction are obtained with the same parameters as in figure 10.

The other extreme is the maximum transition uncertainty, which can be represented by the least informative transition matrix – a perfectly mixing matrix with elements $Q_{ik} \equiv 1/K$. The information entropy (Shannon Reference Shannon1948)

(4.6)\begin{equation} S(\boldsymbol{Q})={-}\sum_{i=1}^{K} \sum_{k=1}^{K} Q_{i k} \ln Q_{i k} \end{equation}

of the permutation matrix vanishes. In contrast, the maximum information entropy $K\ln K$ is obtained from the perfectly mixing transition matrix with equal elements $Q_{ik} \equiv 1/K$. Here, the knowledge of the current state has no predictive value for the future population. For the current case, the reference model has an entropy $S_{CNM} = 7.0586$, which is much smaller than the upper bound $33 \ln 33 \approx 115.38$. The proposed model minimises entropy, $S_{dCNM} = 0$. Thus, our novel clustering measurably increases the prediction accuracy.

4.3. The quasi-periodic flow regime at ${{Re}} = 330$

The clustered quasi-periodic flow of the CNM and dCNM are shown in figure 12. We set $K=10$ for the state space clustering and $\beta = 0.80$ for the subsequent clustering, which proves adequate for capturing the quasi-periodic dynamics and ensuring clarity in visualisation. The choice of $\beta$ and other results with different values of $\beta$ are discussed in Appendix C. In this case, the normalised transverse cluster size vector $\boldsymbol {\hat {R}}^{\mathcal {T}} = [0.1043, 0.1046,0.1056,0.1068,0.0889,0.1073,0.1061,0.1050,0.1047,0.0668]^{\intercal }$ corresponds to the number of sub-clusters $\boldsymbol {L} = [7,7,7,7,6,7,7,7,7,5]^{\intercal }$. In the three-dimensional subspace, the snapshots collectively form a hollow cylinder. The system's dynamics is chiefly governed by two underlying physical phenomena: a cyclic behaviour that synchronises with natural vortex shedding and a quasi-stochastic component responsible for introducing variations between cycles, which is, in turn, synchronised with the oscillator amplitude.

Figure 12. Same as figure 10 but for the quasi-periodic flow at ${{Re}} = 330$. (a) The CNM result with the same number of centroids as the dCNM. (b) The dCNM result with $K = 10$ and $\beta = 0.80$.

The centroid distribution of the CNM reveals that the clustering algorithm fails to distinguish between the shedding dynamics and inter-cycle variations. It uniformly groups them based solely on spatial topology. Nonetheless, the CNM centroids effectively capture the cyclic behaviour, as there exist deterministic transitions between adjacent centroids within an orbit, forming a limit cycle structure akin to the ‘ear’ of the Lorenz system. However, this centroid distribution inadequately models the quasi-stochastic component, as it overlooks the inter-cycle transitions. To comprehensively represent this dynamic, clear transitions between the limit cycles are essential. The clustering process obscures these transitions, causing the quasi-stochastic behaviour to resemble a random walk governed by a fully stochastic process. In essence, the clustering process cannot differentiate between the random jumps in the Lorenz system and the quasi-stochastic behaviour in this flow regime. This explains why the CNM often struggles with multifrequency problems. In contrast, the dCNM centroids automatically align along the axial direction of the cylinder with equidistant circumferential spacing, resulting in a greater number of centroid orbits compared with the CNM. This enhancement enables the accurate resolution of inter-cycle variations. For the quasi-stochastic behaviour, the denser and occasionally overlapping centroids in the axial direction ensure precise spatial representation of the transitions between the limit cycles. Additionally, this behaviour can be further constrained by the dual indexing approach for long-time-scale periodicity, eliminating random jumps and ensuring accurate transitions between limit cycles.

The cluster transition matrices of the quasi-periodic flow regime are illustrated in figure 13. The quasi-periodic dynamics is evident from $\boldsymbol {Q}$, which displays dominant transition probabilities corresponding to deterministic cyclic behaviour and minor wandering transitions signifying inter-cycle variations. Cluster $\mathcal {C}_{4}$ serves as the transition cluster with two destination clusters, $\mathcal {C}_{5}$ and $\mathcal {C}_{10}$. The two destination clusters have similar transition probabilities since they are visited for comparable times during the quasi-periodic transitions. The transition cluster bridges the deterministic cluster chains $\mathcal {C}_{1} \to \mathcal {C}_{2} \to \mathcal {C}_{3} \to \mathcal {C}_{4}$ and $\mathcal {C}_{6} \to \mathcal {C}_{7} \to \mathcal {C}_{9} \to \mathcal {C}_{1}$ as two different limit cycles through two short cluster chains: $\mathcal {C}_{4} \to \mathcal {C}_{5} \to \mathcal {C}_{6}$ and $\mathcal {C}_{4} \to \mathcal {C}_{10} \to \mathcal {C}_{6}$. These two limit cycles alternate with a fixed order, ultimately forming an extended cluster chain that constitutes the fundamental elements of the long-term periodicity. However, this characteristic is not effectively portrayed in the transition matrix. The purely probabilistic transitions from this matrix can result in arbitrary cluster transitions within the network model, introducing additional transition errors. Since the CNM relies on this cluster-level matrix, these transition errors present a notable challenge. While the transition tensors $\mathcal {Q}$, which resolve the refined centroid transitions, mitigate this issue, we further discuss the transition tensors and the corresponding centroid transition matrices in Appendix D. The time matrix $\boldsymbol {T}$ reveals that the transitions within a cyclic behaviour possess a generally similar time scale, with residence times in adjacent clusters changing smoothly, indicating the presence of a gradually evolving limit cycle.

Figure 13. Same as figure 5 but for the quasi-periodic flow at ${{Re}} = 330$. (a) Transition probability matrix $\boldsymbol {Q}$. (b) Transition time matrix $\boldsymbol {T}$.

The original and reconstructed trajectories using the CNM and dCNM for the quasi-periodic flow regime are displayed in figure 14. The reconstruction is achieved with the same parameters as in figure 12. As anticipated, the trajectory reconstructed by the CNM undergoes substantial deformation, featuring discontinuous cyclic behaviours and a serrated trajectory. Conversely, the dCNM produces cleaner cyclic behaviours with more noticeable variations. The reconstructed trajectory accurately replicates the intersecting limit cycles and guides the inter-cycle transition with reduced spatial errors. These observations highlight the ability of the dCNM centroids to capture significant dynamics without assuming any prior knowledge of the data set. A kinematic comparison between the POD reconstruction and the dCNM reconstruction is presented in Appendix E.

Figure 14. Same as figure 11 but for the quasi-periodic flow at ${{Re}} = 330$. (a) The CNM reconstruction and (b) the dCNM reconstruction are obtained with the same parameters as in figure 12.

In the following sections we shift our focus to the temporal aspects. The CNM uses the transition matrices to predict the next destination state for each step. The quasi-periodic feature will be obscured by the stochastic transition probability matrix and the missing historical information. In contrast, the dCNM preserves the transition sequence by embedding the sub-clusters (see Appendix D). Initially, we explore the cluster and trajectory segment affiliation for each snapshot in both the original data set and the dCNM reconstruction to illustrate the accuracy of transition dynamics, as depicted in figure 15. We maintain the same parameters as those used in figure 12 for the reconstruction. The affiliation of the original data reveals that the dual clustering effectively represents the quasi-periodic dynamics. The transition dynamics exhibit significant regularity, with centroids being sequentially and periodically visited, confirming deterministic transitions. Each period of centroid visits corresponds to an extended cluster chain, encompassing multiple centroid orbits and capturing cycle-to-cycle variations. The periodic visits of these extended cluster chains are instrumental in determining the long-time-scale periodicity. These transition characteristics are fully preserved by the dCNM due to the dual indexing constraint. In this case, each departure centroid corresponds to only one destination centroid, eliminating the stochastic transition in the model and mitigating the transition errors.

Figure 15. Transition illustrated with the temporal evolution of the cluster and trajectory segment affiliation of the quasi-periodic flow at ${{Re}} = 330$. The vertical direction represents the cluster-level transition and the horizontal direction indicates the trajectory segments inside this cluster. The transition with black markers represents the CFD data and the transition with red markers represents the reconstructed dynamics by the dCNM. The $x$ axis is the non-dimensionalised time $t$, the $y$ axis is the trajectory segment affiliation $l$ and the $z$ axis is the cluster affiliation $k$. The reconstruction is achieved from the same parameters as in figure 12.

Envelope demodulation can clearly reveal the long-time-scale behaviour and is more efficient in reflecting the quasi-periodic dynamics. We analyse the envelope spectrum of the streamwise fluctuation velocity $u'_{x}$ from the data set, CNM, high-order CNM and dCNM, as depicted in figure 16. The spectrum of the data set exhibits a dominant frequency $f=0.05$, representing long-time-scale periodicity. However, the CNM spectrum shows significant noise and lacks a clear dominant frequency due to frequent transition errors. This observation supports the CNM's limitation in capturing multifrequency dynamics effectively. By incorporating the historical information, the high-order CNM demonstrates superior performance by producing a cleaner spectrum closely aligned with the CFD data's dominant frequency. Remarkably, the dCNM outperforms all other models by precisely reconstructing both the frequency and amplitude while minimising noise.

Figure 16. The envelope spectrum of the streamwise fluctuation velocity $u'_{x}$, obtained by the surface average in $x=5D$. The dCNM reconstruction is achieved with the same parameters as in figure 12 and the CNM reconstruction is achieved with the same number of centroids as the dCNM.

4.4. The chaotic flow regime at ${{Re}} = 450$

The comparison between the CNM and dCNM clustering with $K = 10$ and $\beta = 0.40$ for the chaotic flow are illustrated in figure 17. This value of $\beta$ is the sweet point between model complexity and model accuracy for this test case. Discussions with different values of $\beta$ for this flow regime are presented in Appendix C. The normalised transverse cluster size vector $\boldsymbol {\hat {R}}^{\mathcal {T}} = [0.1068,0.1073,0.1063,0.1068,0.1061, 0.0966,0.1001,0.0954,0.0984,0.0763]^{\intercal }$ corresponds to the number of sub-clusters $\boldsymbol {L} = [22,22,22,22,22,20,21,20,20,8]^{\intercal }$. As the dynamics become more complex, the snapshots form a chaotic cloud, which is driven by numerous cyclic behaviours of different scales and indicates irregular three-dimensional vortex shedding. The CNM continues to cluster the data set primarily based on spatial properties, essentially dividing the chaotic cloud into different segments in an evenly distributed manner. Figure 17(a) illustrates this process, with the uniformly spread centroids capturing only part of one whole cyclic behaviour, limiting their ability to resolve the multiscale dynamics. The dCNM centroids concentrate in regions of rich dynamics, enabling a more comprehensive resolution of the cyclic behaviours. These centroids, in various combinations, form the basis of multi-frequency and multiscale cyclic behaviour. Even after sparsification, the dCNM centroids can encompass a significant amount of scale diversity by merging only those that are spatially close to each other.

Figure 17. Same as figure 10 but for the chaotic flow at ${{Re}} = 450$. (a) The CNM result with the same number of centroids as the dCNM. (b) The dCNM result with $K = 10$ and $\beta = 0.40$.

The cluster transition matrices of the chaotic flow are illustrated in figure 18. The probability matrix $\boldsymbol {Q}$ in figure 18(a) shows that most of the clusters have three or more destination clusters, indicating complex transition dynamics among them. Several dominant transition loops are identifiable, such as the large-size cluster chain $\mathcal {C}_{1} \to \mathcal {C}_{2} \to \mathcal {C}_{3} \to \mathcal {C}_{4} \to \mathcal {C}_{5} \to \mathcal {C}_{6} \to \mathcal {C}_{1}$, the mid-size cluster chains $\mathcal {C}_{1} \to \mathcal {C}_{2} \to \mathcal {C}_{3} \to \mathcal {C}_{4} \to \mathcal {C}_{5} \to \mathcal {C}_{1}$ and $\mathcal {C}_{3} \to \mathcal {C}_{4} \to \mathcal {C}_{5} \to \mathcal {C}_{6} \to \mathcal {C}_{7} \to \mathcal {C}_{3}$, and the small-size cluster chain $\mathcal {C}_{6} \to \mathcal {C}_{7} \to \mathcal {C}_{8} \to \mathcal {C}_{6}$. These cluster chains with different lengths represent cyclic loops at different scales. The small number of chains facilitates human understanding of the transition dynamics but is insufficient for accurately capturing the dynamics. The dominant loops have their key transition clusters inside, from which they can randomly jump into each other by choosing periodic or stochastic routes. This is where the transition error often occurs. The time matrix $\boldsymbol {T}$ in figure 18(b) shows the difference in the transition times between different types of transitions. Regarding the main loops, the time scale changes smoothly within its transitions. However, for the jumps between the loops, the time scale fluctuates considerably, and some transitions can be very large, showing the diversity of the dynamics. Moreover, this observation implies that the distribution density of snapshots differs among clusters. In other words, the distribution of the trajectory segments in different clusters also exhibits significant variations. This explains the necessity for determining the number of sub-clusters in the second-stage clustering based on the deviation $\boldsymbol {R}^{\mathcal {T}}$. The refined transition matrices between the centroids are shown in Appendix D.

Figure 18. Same as figure 5 but for the chaotic flow at ${{Re}} = 450$. (a) Transition probability matrix $\boldsymbol {Q}$. (b) Transition time matrix $\boldsymbol {T}$.

The original and reconstructed trajectories by the CNM and dCNM for this flow regime are shown in figure 19. The reconstruction is achieved with the same parameters as in figure 17. For a clear visualisation, only the trajectories from the first half of the entire time window are plotted. This selection suffices to analyse the precision of the current trajectory, as it contains ample dynamics. We exclude trajectory discrepancies triggered by phase mismatch and focus exclusively on the accuracy of the present trajectory. In the case of the CNM, noticeable disparities exist between the original trajectory and the reconstructed trajectory. These differences include variations in the shape, spatial location and inclination angle of the cyclic loops. These disparities can be attributed to the elimination of small-scale structures and the blending of certain large-scale structures due to the uniform distribution of centroids. Regarding the dCNM, the reconstructed trajectory nearly occupies the entire chaotic cloud, closely resembling the original trajectory. The external and internal geometries are accurately reproduced, capturing both large-scale and small-scale structures. However, despite the improved accuracy, some deformations persist. These deformations arise from the interpolations between the limited centroids during one single cyclic loop. Notably, due to its complexity, achieving a superior reconstruction of a chaotic system often requires more refined centroids compared with a quasi-periodic system. The kinematic comparison with the POD reconstruction is also introduced in Appendix E.

Figure 19. Same as figure 11 but for the chaotic flow at ${{Re}} = 450$. (a) The CNM reconstruction and (b) the dCNM reconstruction are obtained from the same parameters as in figure 17.

Figure 20 shows the auto-correlation function of the CNM, high-order CNM and dCNM. We still normalise this function by $R(0)$, and the time window is chosen from $t = 0$ to $t = 400$, which is sufficient for comparison. For the chaotic flow regime, $R(\tau ) / R(0)$ denotes the kinetic energy level of the time window. Nonetheless, a notable discrepancy arises in the CNM, where the amplitude experiences a distinct decay after the initial few periods. It eventually stabilises with minimal variation, primarily due to the distorted reconstructed trajectory and transition errors. This limited variance is indicative of inaccuracies in capturing short-term dynamics, consistent with the absence of historical information. The high-order CNM, which incorporates this historical information, outperforms the CNM in this regard. Its amplitude decays gradually and exhibits variance akin to that of the data set. Additionally, it reveals some peaks with similar time delays, due to the potential introduction of unnecessary long-time-scale periodicity into the reconstruction via the high-order cluster chain. The dCNM also surpasses the CNM with regard to accuracy. Both the amplitude and phase are faithfully retained, with a gradual amplitude decay and more pronounced variation. Eventually, the amplitude diminishes, similar to the original data set. As $\tau$ increases, all three models exhibit some degree of phase delay or lead. This is a consequence of averaged transition times introducing some errors into the model (Li et al. Reference Li, Fernex, Semaan, Tan, Morzyński and Noack2021).

Figure 20. Auto-correlation function of the chaotic flow at ${{Re}} = 450$; here $R$ is normalised by $R(0)$. The thin black curve represents the CFD data and the thick red curve represents the reconstruction from different models. (a) The CNM reconstruction with the same number of centroids as the dCNM. (b) The high-order CNM reconstruction with $L = 10$ and the same number of centroids as the dCNM. (c) The dCNM reconstruction with the same parameters as in figure 17.

4.5. Physical interpretation

One of the major advances of the cluster-based model is its strong physical interpretability. The dCNM also maintains and even enhances this nature while improving the model accuracy. In this section we discuss the physical interpretation of the CROM exemplified for the sphere wake, with particular emphasis on the dCNM.

The cluster-based model spatially coarse grains the snapshots into groups and represents them by centroids to reduce dimensionality. In contrast to the projection-based methodology, such as the POD–Galerkin model, the cluster-based model uses cluster centroids that are linear combinations of several snapshots, and thus, reflects the representative patterns. This feature contributes to its high interpretability. The snapshot dynamics is mapped into the pattern dynamics, followed by the construction of a probabilistic mechanism to reduce temporal dimensionality. The network model, with centroids as nodes and centroid transitions as edges, converts the complex dynamics into pure data analysis. The centroids act as a bridge between the data-driven model and its underlying physical background. Furthermore, it is conceivable that the same model can be easily transferred to analogous pattern dynamics, even with distinct backgrounds, through adjustments to the centroids’ backstory.

The sphere wake offers a concise physical interpretation based on the centroids. The coherent structure evolution governs the flow field and manifests as vortex-shedding events with diverse dynamics. These shedding events can be captured well by a limit cycle, with a set of centroids representing flow patterns at different shedding phases as foundational elements. The cyclic transitions between these specific flow patterns collectively characterise the entire shedding process. The deterministic–stochastic transitions between different shedding events contribute to the overall periodic–chaotic dynamics.

To explain the physical mechanisms of the flow regimes, we propose a chord transition diagram for the cluster transitions and sub-cluster transitions along with centroid visualisation, which provides a comprehensive view of the flow regime. We start with the periodic flow, as shown in figure 21. The cluster probability distribution $P_{k}$ and the cluster size $\boldsymbol {R^u}$ used for visualising the blocks in figure 21(a) are shown in figure 22(a,b). The blocks in figure 21(b) are split based on the sub-clusters, the transverse cluster size is shown in figure 22(c) and the sub-cluster size is shown in figure 22(d). The cluster transition diagram is capable of clearly distinguishing the dynamic behaviour categories. The circumferential arrows along the boundary represent the cyclic behaviours, with centroids transitioning to adjacent destination centroids. The radial arrows crossing the graph signify cycle-to-cycle transitions, with the centroids transitioning to non-adjacent destinations. These arrows usually originate from the transition clusters. The number of radial arrows indicates the dynamic characteristics, with more arrows indicating more chaotic features.

Figure 21. Transition diagram of the periodic flow at ${{Re}} = 300$. The centroids are depicted by the vortex distribution. The vortices are identified by the isosurfaces of $z$ vorticity, with $-1$ for the negative vortices coloured in blue and $1$ for the positive vortices coloured in red. The transition dynamics is depicted by the directed arrows, the size of the arrow tail represents the transition probability and the colour is consistent with the departure block. (a) Cluster transitions. Different blocks represent different clusters, the colour of the block represents the corresponding cluster probability distribution $P_{k}$, and the size of the block represents the cluster size $\boldsymbol {R^u}$. (b) Sub-cluster transitions of $\beta = 0.50$, with transitions specifically departing from $\mathcal {C}_{1}$, corresponding to the red-bordered cluster transition in (a). Blocks with the same colour belong to the same cluster, the colour still represents the cluster probability distribution $P_{k}$, and the size of each block represents the sub-cluster size $\boldsymbol {R}^{\boldsymbol {u}}_{sub}$.

Figure 22. Cluster and centroid properties of the periodic flow at ${{Re}} = 300$. (a) Cluster probability distribution. (b) Normalised cluster size. (c) Normalised transverse cluster size. (d) Normalised sub-cluster size, where the elements from the same cluster sum to unity.

The cluster dynamics is illustrated in figure 21(a). The limit cycle is captured by the clusters $\mathcal {C}_4$ to $\mathcal {C}_{10}$, shown as the deterministic transitions between adjacent clusters. The transient phase is resolved by clusters from $\mathcal {C}_1$ to $\mathcal {C}_3$. However, the stochastic transitions between the first three clusters are in contrast to the slowly varying amplitude in the transient state and are insufficient to represent this deterministic process. The distinct vortex structures of the three centroids also indicate a need for higher resolution. The sub-cluster transitions leaving from $\mathcal {C}_1$ are shown in figure 21(b). The sub-cluster centroids manifest as varying vortex structures, corresponding to the growing amplitudes. Each sub-cluster has only one destination, resulting in deterministic transitions. The stochastic cluster transitions from $\mathcal {C}_1$ to $\mathcal {C}_2$ and $\mathcal {C}_3$ are now terminated into a chain of deterministic sub-cluster transitions, effectively reducing the prediction error.

The transition graph of the quasi-periodic flow regime is shown in figure 23, with the corresponding cluster and centroid properties given in figure 24. In figure 23(a) the flow regime is characterised by the cyclic cluster transitions, with only three non-adjacent transitions, i.e. $\mathcal {C}_{4}$ to $\mathcal {C}_{10}$, $\mathcal {C}_{10}$ to $\mathcal {C}_{6}$ and $\mathcal {C}_{9}$ to $\mathcal {C}_{1}$. The clusters involved in the cyclic transitions exhibit varying vortex structures within one shedding period. Their relatively higher probability distribution, as shown in figure 24(a), suggests dominant flow patterns. For the bifurcating cluster $\mathcal {C}_{4}$, its destination clusters $\mathcal {C}_{5}$ and $\mathcal {C}_{10}$ manifest visible differences in the far wake. Further distinction of transitions from $\mathcal {C}_{4}$ is provided by the sub-clusters, as shown in figure 23(b). The centroids belonging to the same cluster are roughly at the same shedding phase, but exhibit different vortex structures, exemplified by $\mathcal {C}_{4\, 1}$ and $\mathcal {C}_{4\, 2}$. This difference leads to distinct sheddings, such as $\mathcal {C}_{10\, 1}$ and $\mathcal {C}_{5\, 1}$. Consequently, finer dynamic resolution originating from $\mathcal {C}_{4}$ is captured, enabling the deterministic transitions to $\mathcal {C}_{5}$ and $\mathcal {C}_{10}$, respectively.

Figure 23. Same as figure 21 but for the quasi-periodic flow at ${{Re}} = 330$. (a) Cluster transitions. (b) Sub-cluster transitions of $\beta = 0.95$, with transitions specifically departing from $\mathcal {C}_{4}$, corresponding to the marked cluster transition in (a).

Figure 24. Same as figure 22 but for the quasi-periodic flow at ${{Re}} = 330$. (a) Cluster probability distribution. (b) Normalised cluster size. (c) Normalised transverse cluster size. (d) Normalised sub-cluster size, where the elements from the same cluster sum to unity.

The chaotic flow regime exhibits a more complex transition graph, as shown in figure 25. The relative information is illustrated in figure 26. In figure 25(a), similar to other flow regimes, adjacent cluster transitions continue to dominate the flow field, reflecting cyclic behaviours. However, the increasing number of radial arrows with varying transition probabilities indicates chaotic features. Each centroid represents a distinct flow field with different scales of vortex structures, even within the same cyclic cluster transition. This discrepancy indicates that the current flow patterns are inadequate for capturing the entire shedding dynamic. Concerning the sub-cluster transitions in figure 25(b), the increased arrows maintain the transition rhythm but offer more specificity. The flow states can be inferred from the vortex structures surrounding the cyclic diagram. The centroids within the same cluster also represent vortex structures sharing the same shedding phase but exhibiting different scales. Only centroids with similar scales of vortex structures are connected by sub-cluster transitions. The departing sub-clusters are thus restricted, for instance, $\mathcal {C}_{(2\,3)}$ and $\mathcal {C}_{(10\,1)}$ each have only one departing sub-cluster. The diversity of the centroid transitions guarantees diverse flow scales, while simultaneously maintaining a consistent scale within the same cluster loop. Therefore, the dCNM facilitates multiscale fidelity and smoother cyclic behaviours, significantly enhancing the representation capacity of the model.

Figure 25. Same as figure 21 but for the chaotic flow at ${{Re}} = 450$. (a) Cluster transitions. (b) Sub-cluster transitions of $\beta = 0.95$, with transitions specifically departing from $\mathcal {C}_{1}$.

Figure 26. Same as figure 22 but for the chaotic flow at ${{Re}} = 450$. (a) Cluster probability distribution. (b) Normalised cluster size. (c) Normalised transverse cluster size. (d) Normalised sub-cluster size, where the elements from the same cluster sum to unity.

When comparing the cluster transition and the sub-cluster transition in dCNM, it is evident that the state space clustering can be seen to automatically introduce prior knowledge into the model. This prior knowledge includes the inner-state kinematic information, as resolved by the trajectory segments, and the inter-state dynamic information, as resolved by the transitions between the trajectory segments. The incorporation aids in the automatic assignment of refined centroids within each cluster and constrains the probabilistic transition dynamics. In essence, the dCNM can be regarded as having a built-in unsupervised physics-informing process, which results in superior model accuracy.