Dynamics-augmented cluster-based network model

Abstract

In this study we propose a novel data-driven reduced-order model for complex dynamics, including nonlinear, multi-attractor, multi-frequency and multiscale behaviours. The starting point is a fully automatable cluster-based network model (CNM) (Li et al., J. Fluid Mech., vol. 906, 2021, A21) that kinematically coarse grains the state with clusters and dynamically predicts the transitions in a network model. In the proposed dynamics-augmented CNM (dCNM) the prediction error is reduced with trajectory-based clustering using the same number of centroids. The dCNM is first exemplified for the Lorenz system and then demonstrated for the three-dimensional sphere wake featuring periodic, quasi-periodic and chaotic flow regimes. For both plants, the dCNM significantly outperforms the CNM in resolving the multi-frequency and multiscale dynamics. This increased prediction accuracy is obtained by stratification of the state space aligned with the direction of the trajectories. Thus, the dCNM has numerous potential applications to a large spectrum of shear flows, even for complex dynamics.

1. Introduction

Advancements in computational capabilities and flow measurement technologies are producing a rapidly increasing amount of high-fidelity flow data. The coherent spatio-temporal structures of the flow data enable data-driven reduced-order models (ROMs). In terms of kinematics, ROMs furnish simplified descriptions that enrich our understanding of fundamental flow processes (Holmes, Lumley & Berkooz 1996), facilitated by increasingly powerful machine learning methods (Brunton, Noack & Koumoutsakos 2020). The ROMs may also allow the prediction of future states with acceptable accuracy. In the context of flow control, ROMs are serving as efficient tools for designing and testing control strategies, replacing costly high-fidelity simulations with an acceptable trade-off in accuracy (Bergmann & Cordier 2008).

First-principle-based ROMs have historically been the foundation of the ROM community, as only a limited number of large data sets were available at that time. The Galerkin framework is one of the most classic methods in this category. By projecting the Navier–Stokes equations onto a low-dimensional subspace, the Galerkin model elegantly describes the original dynamics, exhibiting self-amplified amplitude-limited dynamics (Stuart 1971; Busse 1991; Noack & Eckelmann 1994). Landau (1944) and Stuart (1958) pioneered the mean-field model, a major progress in first-principle-based ROMs that provides insight into flow instabilities and bifurcation theory. For instance, in the case of a supercritical Hopf bifurcation, mean-field models have been applied to the vortex shedding behind a cylinder (Strykowski & Sreenivasan 1990; Schumm, Berger & Monkewitz 1994; Noack et al. 2003) and high-Reynolds-number turbulent wake flow (Bourgeois, Noack & Martinuzzi 2013). For more complex flows undergoing successive bifurcations, including both Pitchfork and Hopf bifurcations, weakly nonlinear mean-field analysis is applied to the wake of axisymmetric bodies (Fabre, Auguste & Magnaudet 2008), the wake of a disk (Meliga, Chomaz & Sipp 2009) and the fluidic pinball (Deng et al. 2020). Furthermore, in the field of resolvent analysis, the mean-field theory also contributes by decomposing the system into time-resolved linear dynamics and a feedback term involving quadratic nonlinearity (McKeon et al. 2004; Gómez et al. 2016; Rigas et al. 2017).

In contrast to a first-principle ROM, a data-driven version is based on a low-dimensional representation of flow snapshots. Proper orthogonal decomposition (POD) is a commonly used example. Proper orthogonal decomposition begins with the eigenvalue or singular value decomposition of the correlation matrix, yielding a low-dimensional subspace comprising leading orthogonal eigenvectors. This subspace provides an ‘optimal’ Galerkin expansion with minimal average residual in the energy norm. Since Aubry et al. (1988) introduced the groundbreaking POD–Galerkin model for unforced turbulent boundary layers, numerous POD models have emerged for various configurations. Examples include POD models for channel flow (Podvin & Lumley 1998; Podvin 2009), the wake of a two-dimensional square cylinder (Bergmann, Bruneau & Iollo 2009), laminar and turbulent vortex shedding (Iollo, Lanteri & Désidéri 2000) and flow past a circular cylinder with a dynamic subgrid-scale model and variational multiscale model (Wang et al. 2012). There are also various variations of the POD model, e.g. integrating the actuation terms into the projection system for control design (Bergmann & Cordier 2008; Luchtenburg et al. 2009) and balanced POD (Rowley 2005), which is derived from a POD approximation to the product of controllability and observability Gramians to obtain an approximately balanced truncation (Moore 1981). Increasingly powerful machine learning methods can make data-driven ROMs more automated. Examples include the sparse identification of nonlinear dynamics (SINDy) aim at human interpretable models (Brunton, Proctor & Kutz 2016), ROMs with artificial neural networks (San & Maulik 2018; San, Maulik & Ahmed 2019; Zhu et al. 2019; Kou & Zhang 2021), turbulence modelling and flow estimation with multi-input multi-output by deep neural networks (Kutz 2017; Li, Li & Noack 2022) and manifold learning methods (Farzamnik et al. 2023).

In this work we focus on automated data-driven modelling. The starting point is cluster-based ROMs (CROMs), pioneered by Burkardt, Gunzburger & Lee (2006). Clustering is an unsupervised classification of patterns into groups commonly used in data science (Jain & Dubes 1988; Jain, Murty & Flynn 1999; Jain 2010), it is popular in data mining, document retrieval, image segmentation and feature detection (Kim et al. 2022). The foundation of the CROM lies in the cluster-based Markov model (CMM) proposed by Kaiser et al. (2014), which combines a cluster analysis of an ensemble of snapshots and a Markov model for transitions between different flow states reduced by clustering. The CMM has provided a valuable physical understanding of the mixing layer, Ahmed body wakes (Kaiser et al. 2014), combustion-related mixing (Cao et al. 2014) and the supersonic mixing layer (Li & Tan 2020). Nair et al. (2019) applied the cluster-based model to feedback control for drag reduction and first introduced a directed network for dynamical modelling. Building on this concept, Fernex, Noack & Semaan (2021) and Li et al. (2021) further proposed the cluster-based network model (CNM) with improved long-time-scale resolution. Instead of the ‘stroboscopic’ view of the CMM, the CNM focuses on non-trivial transitions. The dynamics is restricted to a simple network model between the cluster centroids, like a deterministic–stochastic flight schedule that allows only a few possible flights with corresponding probabilities and flight times consistent with the data set. Networks of complex dynamic systems have gained great interest, forming an increasingly important interdisciplinary field known as network science (Watts & Strogatz 1998; Albert & Barabási 2002; Börner, Sanyal & Vespignani 2007; Barabási 2013). Network-based approaches are often used in fluid flows (Nair & Taira 2015; Hadjighasem et al. 2016; Taira, Nair & Brunton 2016; Yeh, Gopalakrishnan Meena & Taira 2021; Taira & Nair 2022), in conjunction with clustering analysis (Bollt 2001; Schlueter-Kuck & Dabiri 2017; Murayama et al. 2018; Krueger et al. 2019). The critical structures that modify the dynamical system can be identified by the intra- and inter-cluster interactions using community detection (Gopalakrishnan Meena, Nair & Taira 2018; Gopalakrishnan Meena & Taira 2021).

The CROMs are fully automated, robust and physically interpretable, while the model accuracy is strongly related to the clustering process. The state space is equivalently discretized in the above-mentioned CROMs, leading to a lack of dynamic coverage. For example, the CNM struggles to capture multiscale behaviours such as the oscillations near attractors and the amplitude variations between trajectories. To address this issue, an effective solution is to employ dynamics-augmented clustering to determine the centroid distribution. Inspired by the hierarchical clustering (Deng et al. 2022) and the network sparsification (Nair & Taira 2015), we propose a dynamics-augmented cluster-based network model (dCNM) with an improved resolution of complex dynamics. In this case, the time-resolved dynamics is reflected by the evolution of trajectory segments after the state space is clustered. These segments are automatically identified from cluster transitions and are represented by centroids obtained through segment averaging. A second-stage clustering further refines the centroids, eliminating the network redundancy and also deepening the comprehension of underlying physical mechanisms. The proposed dCNM can systematically identify complex dynamics involved in the case of multi-attractor, multi-frequency and multiscale dynamic systems. Figure 1 provides a comparative illustration of CNM and dCNM in terms of kinematics and dynamics, exemplified by an inward spiral trajectory in a two-dimensional state space.

Figure 1. Principle sketches: the CNM and the dCNM are illustrated using an inward spiral trajectory in a two-dimensional state space with the same number of centroids. The thick solid lines denote cluster divisions and the thin solid lines represent sub-cluster divisions. The centroids are represented by coloured dots and their colours represent their cluster affiliations. The CNM centroids are derived from snapshot averages within each cluster and show uniform geometric coverage, whereas the dCNM centroids incorporate dynamic features and exhibit a weighted distribution. Consequently, dCNM accurately reconstructs the cycle-to-cycle variations and also ensures precise transition sequencing.

The dCNM is first applied to the Lorenz (1963) system as an illustrative example. The Lorenz attractor is notable for the ‘butterfly effect’, showcasing the chaotic dynamics governed by only three ordinary differential equations. Subsequently, we demonstrate the dCNM on the sphere wake of the periodic, quasi-periodic and chaotic flow regimes. The sphere wake is a well-investigated benchmark configuration, serving as a prototype flow of bluff body wakes commonly encountered in many modern applications, for instance, the design of drones and air taxis. Despite the simple geometry, the sphere wake can experience a series of bifurcations with increasing Reynolds number. Along the route to turbulence, the flow system exhibits steady, periodic, quasi-periodic and chaotic flow regimes. The transient and post-transient flow dynamics, characterised by multi-frequency and multiscale behaviours, provide a challenging testing ground for reduced-order modelling.

This paper is organised as follows. In § 2 the clustering algorithm and the different perspectives on the dCNM strategy are described. In § 3 the dCNM is illustrated on the Lorenz system and in § 4 it is demonstrated on the sphere wake of three flow regimes: the periodic flow, the quasi-periodic flow and the chaotic flow. In § 5 the main findings and improvements are summarised and future directions are suggested.

2. Dynamics-augmented CNM

In this section we detail the process of the dCNM. In § 2.1 the $k$-means++ clustering algorithm and its demonstration on the state space are introduced. The second-stage clustering on the trajectory segments is further discussed in § 2.2. In § 2.3 the transition characteristics are described and in § 2.4 different criteria are introduced to evaluate the performance of the proposed model. The variables used in this section are listed in table 1.

Table 1. Table of variables. Subscripts $k$ and $i$ are related to the level of clusters from the state space clustering, and subscripts $l$ and $j$ are related to the level of trajectory segments.

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 1967; Lloyd 1982; Arthur & Vassilvitskii 2007) 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. (2014), 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. 2021; Li et al. 2021). 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. 2021) or add the trajectory supporting points (Hou, Deng & Noack 2022) 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. 2021). 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. 2021). 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 2015) 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. 2022). 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}

3. Lorenz system as an illustrative example

In this section we apply the dCNM to the Lorenz (1963) system to illustrate its superior spatial resolution in handling multiscale dynamics. We also compare it with the CNM (Fernex et al. 2021; Li et al. 2021) of the same rank as a reference.

The Lorenz system is a three-dimensional autonomous system with non-periodic, deterministic and dissipative dynamics that exhibit exponential divergence and convergence to strange fractal attractors. The system is governed by three coupled nonlinear differential equations:

(3.1)\begin{equation} \left. \begin{aligned} \mathrm{d} x/\mathrm{d}t & =\sigma(y-x),\\ \mathrm{d} y/\mathrm{d}t & =x(\rho-z)-y,\\ \mathrm{d}z/\mathrm{d}t & =x y-\beta z. \end{aligned} \right\} \end{equation}

The system parameters are set as $\sigma = 10$, $\rho = 28$ and $\beta = 8/3$. These equations emulate the Rayleigh–Bénard convection. The trajectory of the system revolves around two weakly unstable oscillatory fixed points, forming two sets of attractors, that are loosely called ‘ears.’ These two ears have similar but not identical shapes, with the left ear being rounder and thicker in the toroidal region. The region where the ears overlap is called the branching region. The Lorenz system has two main types of dynamics. One is that the inner loop in each ear varies and oscillates for several cycles. The other is that the inner loop may randomly switch from one ear to another in the branching region and resume oscillatory motion.

We numerically integrate the system using the fourth-order explicit Runge–Kutta method. The time-resolved $10\,000$ snapshots data with $\boldsymbol {u}^m=[x, y, z]^{\intercal }$ are collected at a sampling time step of $\Delta t = 0.015$ with an initial condition of $[-3, 0, 31]^{\intercal }$ (Fernex et al. 2021). This time step corresponds to approximately one-fiftieth of a typical cycle period. The first $5\,\%$ of the snapshots are neglected to reserve only the post-transient dynamics.

Figure 4 shows the phase portrait of the clustered Lorenz system from the CNM and dCNM.

Figure 4. Phase portrait of the clustered Lorenz system from the CNM and dCNM. 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. As a 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.90$.

We set $K =10$ for the state space clustering of the dCNM, which is consistent with previous studies (Kaiser et al. 2014; Li et al. 2021). This number is large enough for the further subdivision of transition dynamics and is also small enough to obtain a simple structure for understanding. The sparsification index $\beta$ is chosen with large numbers as $\beta =0.90$ to allow for a distinct visualisation of the centroids. In addition, since the trajectory in each ‘ear’ is confined to a two-dimensional surface, a high value of $\beta$ is deemed suitable. The normalised transverse cluster size vector $\boldsymbol {\hat {R}}^{\mathcal {T}} = [0.1163, 0.1262,0.1164,0.0921,0.0943,0.0908,0.1116,0.0840,0.0866,0.0817]^{\intercal}$ corresponds to the number of sub-clusters $\boldsymbol {L} = [13,15,13,11,11,11,13,10,10,10]^{\intercal }$.

The two models exhibit notable differences in centroid distribution. The CNM clustering relies solely on the spatial topology in the phase space, evenly dividing the entire attractor and dispersing centroids uniformly throughout the phase portrait. It can be inferred that increasing the number of centroids under this uniform distribution does not lead to substantial changes, merely resulting in a denser centroid distribution. This uniform distribution possesses certain disadvantages regarding the dynamics. First, it unnecessarily complicates the transition rhythm as the deterministic large-scale transition may be fragmented into several stochastic transitions. Second, even with many centroids, it fails to capture the increasing oscillation amplitude between the loops in one ear, as the uniform distribution provides only a limited number of centroid orbits. The same result occurs for the branching region where these limited numbers of centroids usually oversimplify the switch between ears. In contrast, the distribution of the dCNM centroids resembles a weighted reallocation. For the Lorenz system, the state space is stratified along the trajectory direction, leading to a concentrated distribution of the dCNM centroids in the radial direction of the attractor and the branching region, which correspond to the system's primary dynamics. Additionally, varying quantities of the centroids can be observed in the radial direction in the toroidal region, depending on its thickness. In thinner toroidal regions with smaller variations between trajectory segments, the second-stage clustering assigns fewer sub-clusters and, consequently, builds fewer centroids.

The cluster transition matrices, which are a distinctive feature of cluster modelling, are preserved because the dCNM maintains the coarse-grained transitions at the cluster level. Figure 5 illustrates the cluster transition probability matrix $\boldsymbol {Q}$ and the corresponding transition time matrix $\boldsymbol {T}$ to illustrate the significant dynamics of the Lorenz system. It is worth noting that in the case of the CNM with an equivalent number of centroids, the matrices become considerably larger, which diminishes their readability and interpretability. The matrices reveal three distinct cluster groups. The first group comprises clusters $\mathcal {C}_{1}$ and $\mathcal {C}_{2}$, which resolve the branching region and exhibit similar transition probabilities to clusters $\mathcal {C}_{3}$ and $\mathcal {C}_{7}$. The branching region is further linked to different ears and is crucial to the attractor oscillation. Clusters $\mathcal {C}_{1}$ and $\mathcal {C}_{2}$ can be referred to as flipper clusters (Kaiser et al. 2014), representing a switch between the different groups. The equivalent transition probability from $\mathcal {C}_{2}$ is consistent with the random jumping behaviour of the two ears. The other two groups demonstrate an inner-group circulation corresponding to the main components of the two ears, exemplified by the cluster chains $\mathcal {C}_{3} \to \mathcal {C}_{4} \to \mathcal {C}_{5} \to \mathcal {C}_{6}$ and $\mathcal {C}_{7} \to \mathcal {C}_{8} \to \mathcal {C}_{9} \to \mathcal {C}_{10}$. These chains exhibit deterministic transition probabilities that resolve the cyclic behaviour. In the second-stage clustering these two groups are further categorised into numerous centroid orbits. Moreover, the transition time matrix resolves the variance in the transition times, with significantly shorter transition times observed in the cyclic groups compared with transitions involving the flipper clusters.

Figure 5. Transition matrices of the Lorenz system. The colour bar indicates the values of the terms. (a) Transition probability matrix $\boldsymbol {Q}$. (b) Transition time matrix $\boldsymbol {T}$.

The original and reconstructed trajectories in the phase space are directly compared. We focus solely on the spatial resolution, disregarding phase mismatches during temporal evolution. Figure 6 shows the original Lorenz system and the reconstruction by the CNM and dCNM with the same parameters as in figure 4. To ensure clarity, we select a time window from $t=0$ to $t=30$ for the trajectories and employ spline interpolation for a smooth reconstruction. Inaccurate or non-physical centroid transitions, along with incomplete dynamic coverage, can lead to substantial deformations in the reconstructed trajectory. As expected, the dCNM provides a more accurate reconstruction than the CNM. The CNM uses a finite number of centroid orbits to represent oscillating attractors, converting slow and continuous amplitude growth into limited and abrupt amplitude jumps. Furthermore, the CNM may group one continuous snapshot loop into clusters belonging to different centroid orbits, often when these clusters are adjacent to each other. This can lead to unnecessary orbit-crossing centroid transitions and result in non-physical radial jumps in the reconstructed trajectory. In contrast, the dCNM provides more comprehensive dynamic coverage, resolving more cyclic behaviour with additional centroid orbits. Dual indexing also guarantees accurate centroid transitions. The radial jumps are eliminated, as departing centroids can only transition to destination centroids within the same centroid orbits. Consequently, oscillations are effectively resolved by the centroid orbits, and transitions between them are constrained by densely distributed centroids in the branching region, ensuring a smoothly varied oscillation.

Figure 6. Trajectory of the Lorenz system. 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 performed with the same parameters as in figure 4.

The auto-correlation function is computed to reflect the model accuracy, as shown in figure 7. In the original data set, the normalised auto-correlation function $R(\tau ) / R(0)$ vanishes smoothly as $\tau$ increases, and the variance between the periodic behaviour can be clearly observed. However, the CNM reconstruction captures only the first four periods of oscillation dynamics. As $\tau$ increases, there is a sudden amplitude decay accompanied by a phase mismatch. This can be attributed to amplitude jumps between the centroid loops and commonly occurring orbit-crossing transitions. In contrast, the dCNM reconstruction accurately captures both the amplitude and frequency of the oscillation dynamics, demonstrating robust and precise long-time-scale behaviours.

Figure 7. Auto-correlation function for $\tau \in [0, 30)$ of the Lorenz system. The thin black curves represent the original data set and the thick red curves represent the models: (a) CNM and (b) dCNM.

4. Dynamics-augmented modelling of the sphere wake

In this section we demonstrate the dCNM for the transient and post-transient flow dynamics of the sphere wake. The numerical method for obtaining the flow field data set and the flow characteristics is presented in § 4.1. The performance of the dCNM for the periodic, quasi-periodic and chaotic flow regimes is evaluated in §§ 4.2, 4.3 and 4.4, respectively. The physical interpretation of the modelling strategy is discussed in § 4.5.

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 2018; Lorite-Díez & Jiménez-González 2020). 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 (1999) and Rajamuni, Thompson & Hourigan (2018). 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 2003; Eshbal et al. 2019). The length of this wake grows linearly with $\ln ({{Re}})$. When ${{Re}}$ surpasses $130$ (Taneda 1956), 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. 2008), 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 1999). At a subsequent Hopf bifurcation around ${{Re}} = 270 \sim 272$ (Johnson & Patel 1999; Fabre et al. 2008), 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. 2019), 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 1999; Pan et al. 2018; Eshbal et al. 2019), 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 1948)

(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. 2021).

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.

5. Conclusions and outlook

We propose an automatable data-driven ROM for nonlinear dynamics. This model can resolve the periodic, quasi-periodic and chaotic dynamics of the sphere wake featuring multi-frequency and multiscale behaviours. The starting point is the CNM (Fernex et al. 2021; Li et al. 2021), which is an automated framework employing clustering and network science. The dynamics within the CNM are described using a deterministic–stochastic approach on a network, where centroids act as nodes, and transitions serve as edges. However, the clustering process in the CNM relies on a uniform geometric coverage of the snapshot data, agnostic of the temporal dynamic relevance. For multi-frequency dynamics, this can result in large prediction errors. One example is the long transient to a limit cycle. Here, the slow increase in the radius requires a finer resolution than the robust angular dynamics. Hence, the CNM can be expected to be more accurate if the centroids are much denser in the radial direction than in the angular motion. This idea is incorporated in the proposed dynamics-augmented CNM (dCNM). The model can automatically stratify the state space along the trajectory direction.

The dCNM was applied to the Lorenz system (in § 3) and the three-dimensional sphere wake (in § 4), with $K=10$ clusters for the coarse-graining of the state space. The Lorenz system features oscillatory dynamics, presented as two ‘ears’ consisting of many unstable orbits, and stochastic dynamics, presented as random switching between the ‘ears.’ For the future state, the phase can be accurately predicted, but the amplitude requires a higher resolution. The CNM is only capable of reconstructing limited loops of the cyclic behaviours and their related transitions in the branching area. Non-physical radial jumps also occur due to transition errors. On the other hand, the dCNM coarsely resolves the deterministic phases but accurately resolves the slowly varying amplitude. The attractor oscillations are distinctly defined, and the transitions in the branching region are subsequently constrained. For the transient and post-transient dynamics of the periodic sphere wake, the dCNM accurately resolves the slowly growing amplitude between the cyclic behaviours. Regarding the quasi-periodic sphere wake, the dCNM successfully captures both the periodic behaviour and cycle-to-cycle variations. Notably, it discerns intrinsic deterministic transition behaviours, which are often misinterpreted as stochastic transitions by the CNM. For the chaotic flow dominated by unstable periodic orbits with varying scales, the dCNM accurately distinguishes between these orbits and captures their transitions. Even after sparsification, chaotic features remain preserved, with transition dynamics demonstrating stochastic characteristics. Overall, these findings underscore the notable improvement of the dCNM in capturing and accurately representing multi-frequency and multiscale dynamics.

The dCNM offers several advantages over other reduced-order modelling strategies. It preserves the advantages of previous cluster-based approaches and adds the following new noteworthy features.

1. (i) The prediction error is minimised. The slow evolution of amplitude oscillations, the deterministic quasi-periodic dynamics and the stochastic chaotic dynamics can be automatically resolved without any prior knowledge.

2. (ii) The model complexity is significantly reduced as the number of non-trivial transitions is mitigated by design. The CNM often requires more clusters and a higher order to achieve similar accuracy, with more complex cluster transition relationships.

3. (iii) The physical interpretability of the model is enhanced.

Our results suggest entropy as a guiding principle of future cluster-based models. We characterise the prediction accuracy of the CNM with the Kullback–Leibler entropy of the transition matrix, called transition entropy for brevity. This transition entropy is significantly reduced for the dCNM as compared with the CNM for the same number of centroids. Further improvements may be expected by optimizing the $\beta$ parameter. Thus, the dCNM development from snapshots can be fully automated. The results even inspire a new clustering based on the prediction uncertainty expressed with the transition entropy. An intrusive framework with the Navier–Stokes propagator may be a further avenue for improvement.

The dCNM may be compared with the POD-based Galerkin method. By construction, centroids are physically interpretable as coherent structures. In contrast, POD models have no intrinsic meaning and typically mix different frequencies. However, in select cases, the Galerkin method may yield deep insights into linear and nonlinear dynamics. Examples are the Galerkin mean-field models for the effect of forcing on a vortex shedding (Semaan et al. 2016), for a single oscillator (Noack et al. 2003) and for frequency cross-talk (Luchtenburg et al. 2009). The authors work on combining the advantages of clustering and POD in human-interpretable dynamic models.