# Trajectory-optimized cluster-based network model for the sphere wake

Chang Hou (侯昶)  School of Mechanical Engineering and Automation, Harbin
Institute of Technology, Shenzhen 518055, Peoples’ Republic of China  Nan Deng
(邓楠)  School of Mechanical Engineering and Automation, Harbin Institute of
Technology, Shenzhen 518055, Peoples’ Republic of China<EMAIL_ADDRESS>Bernd R. Noack  School of Mechanical Engineering and Automation, Harbin
Institute of Technology, Shenzhen 518055, Peoples’ Republic of China  Hermann-
Föttinger-Institut, Technische Universität Berlin, Müller-Breslau-Straße 8,
D-10623 Berlin, Germany<EMAIL_ADDRESS>
###### Abstract

We propose a novel trajectory-optimized Cluster-based Network Model (tCNM) for
nonlinear model order reduction from time-resolved data following Li _et al._
["Cluster-based network model," J. Fluid Mech. 906, A21 (2021)] and improving
the accuracy for a given number of centroids. The starting point is
$k$-means++ clustering which minimizes the representation error of the
snapshots by their closest centroids. The dynamics is presented by ‘flights’
between the centroids. The proposed trajectory-optimized clustering aims to
reduce the kinematic representation error further by shifting the centroids
closer to the snapshot trajectory and refining state propagation with
trajectory support points. Thus, curved trajectories are better resolved. The
resulting tCNM is demonstrated for the sphere wake for three flow regimes,
including the periodic, quasi-periodic, and chaotic dynamics. The
representation error of tCNM is 5 times smaller as compared to the
approximation by the closest centroid. Thus, the error is at the same level as
Proper Orthogonal Decomposition (POD) of same order. Yet, tCNM has distinct
advantages over POD modeling: it is human interpretable by representing
dynamics by a handful of coherent structures and their transitions; it shows
robust dynamics by design, i.e., stable long-time behavior; and its
development is fully automatable, i.e., it does not require tuneable auxiliary
closure and other models.

††preprint: AIP/123-QED

##  I Introduction

Reduced-order modeling is a cornerstone of theoretical fluid mechanics and a
critical enabler for understanding, estimation, dynamic modeling, control, and
optimization. For over one century, a myriad of vortex models explain
fundamental aspects of coherent structure dynamics (Pullin and Saffman, 1998;
Cottet and Koumoutsakos, 2000), e.g., the vortex pair behind a cylinder
(Föppl, 1913) and the formation of a vortex street (von Kármán, 1911). Over
hundred years ago, the Galerkin method (Galerkin, 1915) served as a foundation
for computational fluid mechanics (Fletcher, 1984) and modern reduced-order
models (Holmes  _et al._ , 2012). In the 50s and 60s, Stuart (1971) employed
stability modes in the Galerkin framework for an elegant description of self-
amplified amplitude-limited dynamics. Another wave of physical Galerkin models
starts with this stability theory framework. Mathematical Galerkin models
(Busse, 1991; Noack and Eckelmann, 1994) aim at more general dynamics by
employing complete Hilbert space bases. These modes are typically only
available for simple geometries.

The Proper Orthogonal Decomposition (POD), particularly the snapshot method
(Sirovich, 1987), opened an increasingly popular avenue to data-driven
reduced-order modeling (Taira  _et al._ , 2017), fueled by the increasing
availability of high-quality numerical and experimental data. Aubry  _et al._
(1988) pioneered POD modeling for the unforced turbulent boundary layer, one
of the most complex flows to start with. An avalanche of POD models followed
for a myriad of configurations. Podvin and Lumley (1998) proposed a low-
dimensional model for the minimal channel flow unit for the purpose of
physical understanding. She continued to develop an accurate high-dimensional
POD model for the wall region of a turbulent channel flow (Podvin, 2009).
Bergmann, Bruneau, and Iollo (2009) further proposed a pressure extended POD
model with no additional cost and tested that on the two-dimensional square
cylinder wake.

The first principles based Galerkin model has been completed with a pressure-
term representation (Noack, Papas, and Monkewitz, 2005) and modal energy-flow
analysis. The unresolved fluctuations may be incorporated with closures as
recently reviewed by Ahmed  _et al._ (2021). Non-intrusive Galerkin models
avoid any first-principle information. In this category, Sparse Identification
of Nonlinear Dynamics (SINDy) aim at human interpretable models (Brunton,
Proctor, and Kutz, 2016). In addition to the spectrum of dynamical system
identification, also different Galerkin expansions have been developed.
Examples are the balanced POD by Rowley (2005), the multiscale-POD by Mendez,
Balabane, and Buchlin (2019), the Dynamic Mode Decomposition (DMD) by Rowley
_et al._ (2009); Schmid (2010a); Kutz  _et al._ (2016), and the recursive DMD
by Noack  _et al._ (2016). All Galerkin models tend to be fragile, inaccurate
for transient dynamics, and restricted to a narrow range of operating
conditions. A key reason is the elliptical Galerkin ansatz with global modes
neglecting the convective dynamics(Noack, 2016). The advances in machine
learning algorithms also attracted widespread concern for data-driven
modeling, e.g., with Artificial Neural Networks (ANN) (San, Maulik, and Ahmed,
2019; Zhu  _et al._ , 2019), Deep Neural Networks (DNN) (Kutz, 2017).

A fundamentally different data-centric avenue, cluster-based modeling, rapidly
gains traction. While the Galerkin method is based on a linear superposition
of modes, the cluster-based approach conceptualizes the dynamics as ‘flights’
in a state-space which is discretized by select centroids. This avenue was
pioneered by Burkardt, Gunzburger, and Lee (2006). Kaiser  _et al._ (2014)
proposed the Cluster-based Reduced-Order Modeling (CROM) method. CROM
partitions the flow data into clusters and describes the flow dynamics with a
Cluster-based Markov Model (CMM). The temporal evolution is modeled as a
probabilistic Markov model of the transition dynamics. The state vector of
cluster probabilities may initially start in a single cluster but eventually
diffuses to a fixed point representing the post-transient attractor. CMM has
provided valuable physical insights for the Ahmed body wake (Kaiser  _et al._
, 2014), combustion-related mixing (Cao  _et al._ , 2014), and the supersonic
mixing layer (Li and Tan, 2020).

A challenge for CMM is the temporal evolution: the state quickly diffuses over
the whole attractor, often within one typical time period. This loss of
dynamic information is alleviated by Cluster-based Network Model (CNM) from
time-resolved data (Li  _et al._ , 2021; Fernex, Noack, and Semaan, 2021). In
CNM, the dynamics is modeled on a directed network(Barabási, 2013; Arenas  _et
al._ , 2008) based on motions between centroids. This motion may be pictured
in an airport analogy. The centroids constitute the nodes, like airports, and
the flights represent the edges, like a deterministic-stochastic flight
schedule. The schedule allows only a few possible flights with corresponding
probabilities and flight times consistent with the data.

The CNM can be seen as an extension of the CMM using the network model instead
of the standard Markov model to describe the transient dynamics. Deng  _et
al._ (2022) later applied a scale-dependent hierarchical clustering to the CNM
under a mean-field consideration and generalized the CNM in the case of
multiple invariant sets. Nair  _et al._ (2019) applied clustering as the basis
for coarse-graining the control law before optimization.

The accuracy of the CNM depends on the centroids. Clustering minimizes the
representation error of the snapshots to the closest centroids, but ignores
the possibility of linear interpolations (flights). Evidently, the averaging
of the snapshots in a cluster leads to centroids where distinct vortex
structures may be smoothed out. In addition, the flights between centroids may
not accurately track the real trajectories. In this study, centroids are
shifted closer to the snapshot trajectories, and constraints between centroids
are allowed to get an accurate state propagation. The resulting trajectory-
optimized Cluster-based Network Model (tCNM) shows better consistency with the
original dataset as compared to the CNM, and leads to a better understanding
of the physical mechanisms involved in the flow dynamics.

Table 1 compares the discussed data-driven reduced-order models. The starting
point is the POD model with a deterministic description in the corresponding
subspace, followed by the CROM (CMM) with a stroboscopic temporal prediction
of discrete states. The CNM with time-continuous uniform motion on a network
of routes between the nodes, and the tCNM with the trajectory optimization
based on centroid shift and inter-cluster constraint.

Table 1: Comparison between POD, CROM, CNM and tCNM. For details see text. Model | POD | CROM | CNM | tCNM
---|---|---|---|---
Data driven | ✓ | ✓ | ✓ | ✓
Time resolved data |  | ✓ | ✓ | ✓
Intrusive formulation exists | ✓ |  |  |
Robust long time behavior |  | ✓ | ✓ | ✓
Accurate long time behavior |  |  | ✓ | ✓
Trajectory tracking |  |  |  | ✓

In analogy to CFD models, the POD model can be conceptualized as the Large
Eddy Simulations (LES) describing the coherent structures with the
statistically rare events being underrepresented, and the CROM is similar to
the unsteady Reynolds-averaged Navier-Stokes equations describing the
transient mean flow. The CNM is seen as the LES with the subgrid-scale
vortices being simply modeled, while the tCNM can be regarded as the LES with
a more accurate subgrid-scale model.

As demonstration, we consider the well-investigated three-dimensional
incompressible flow around a sphere. Despite the simple geometry, the wake is
fully three-dimensional and capable of featuring complicated vortical
structures. Many studies have investigated the flow structure in the wake of
spheres using dye or smoke visualizations (Möller, 1937; Achenbach, 1974;
Sakamoto and Haniu, 1990; Leweke  _et al._ , 1999; Chrust, Goujon-Durand, and
Wesfreid, 2013) and numerical simulations (Johnson and Patel, 1999; Kim, Kim,
and Choi, 2001; Ploumhans  _et al._ , 2002; Yun, Kim, and Choi, 2006;
Giacobello, Ooi, and Balachandar, 2009; Rajamuni, Thompson, and Hourigan,
2018). As the Reynolds number $Re$ increases, the sphere wake undergoes a
series of bifurcations and transitions on the route to turbulence. This
transition scenario is an exquisite testing ground for reduced-order modeling.
In this study, we focus on three representative flow regimes of the sphere
wake, including periodic, quasi-periodic, and chaotic flow.

The manuscript is organized as follows. Section II describes the numerical
plant and the computational validation of the simulation. Section III proposes
the trajectory-optimized Cluster-based Network Model. In section IV, the
detailed flow characteristics of sphere wake at different $Re$ are introduced.
In section V, tCNM is performed on the sphere wake featuring three dynamics,
including periodic flow, quasi-periodic flow, and chaotic flow, and compared
to the CNM. Section VI summarizes the main findings and gives some suggestions
for future directions.

##  II Flow configuration and numerical method

###  II.1 Flow configuration

The three-dimensional sketch of the flow configuration is presented in figure
1 (a).

Figure 1: (a) Three-dimensional sketch of the flow configuration. $D$ is the
diameter of the sphere. (b) Grid expanded view of plane $z=0$. (c) Grid
expanded view of plane $x=0$.

A sphere of diameter $D$ is placed in a uniform flow with streamwise velocity
$U_{\infty}$. The computational domain is described in a Cartesian coordinate
system where the origin is placed in the center of the sphere, and the
$x$-axis is aligned with the streamwise direction. The whole domain is a
cylinder tube, and its sizes in the $x$, $y$, $z$ directions are 80$D$, 10$D$,
and 10$D$, respectively. The inlet is 20$D$ upstream from the center of the
sphere. The choice of these domain parameters ensures a balance of minimal
distortion of the flow structures due to the outer boundary conditions and a
reasonable computational cost. The computational domain is discretized into a
number of hexahedral grids, figure 1 (b) and (c) illustrate the representative
expanded view of the grids around the sphere.

Based on the Cartesian coordinates, the net forces acting on the sphere can be
decomposed into three components $F_{\alpha}$, where $\alpha=x,y$ and $z$.
These forces are non-dimensionalized as

$C_{\alpha}=\frac{F_{\alpha}}{1/2\rho U_{\infty}^{2}S},$ (1)

where $C_{\alpha}$ represents the streamwise force coefficient $C_{x}$,
transverse force coefficients $C_{y}$ of the $y$-axis and $C_{z}$ of the
$z$-axis, $\rho$ is the flow density, $S=\pi D^{2}/4$ is the projected surface
area of sphere in streamwise direction. The total drag force coefficient
$C_{d}=C_{x}$. For an axisymmetric sphere, the lift coefficient can be in any
direction in the $yz$ plane (Zhao, 2021), and the total lift force coefficient
$C_{l}$ can be defined as

$C_{l}=\sqrt{C_{y}^{2}+C_{z}^{2}}.$ (2)

Furthermore, the Strouhal number $St$ can be written as

$St=\frac{Df}{U_{\infty}},$ (3)

where $f$ is the characteristic frequency.

###  II.2 Direct Numerical Simulation

The fluid flow is simulated by solving the non-dimensionalized incompressible
Navier-Stokes equations,

$\displaystyle\frac{\partial\bm{u}}{\partial t}+\bm{u}\cdot\nabla\bm{u}+\nabla
p-\frac{1}{Re}\nabla^{2}\bm{u}=0,$ (4) $\displaystyle\nabla\cdot\bm{u}=0,$ (5)

where $\bm{u}$ denotes the velocity vector $(u,v,w)$ referred to a cylindrical
$(x,r,\theta)$ coordinates system whose origin is placed at the center of the
sphere, $p$ is the static pressure, $Re=U_{\infty}D/\nu$ and $\nu$ is the
kinematic viscosity. The scales for the non-dimensionalization are
characteristic length $D$, free stream velocity $U_{\infty}$, time units
$D/U_{\infty}$ and pressure scales $\rho U_{\infty}^{2}$.

The commercially available CFD package ANSYS Fluent 15.0 is employed in the
present study to discretize the governing equations using the cell-centered
Finite Volume Method (FVM). The Pressure-Implicit Split-Operator (PISO)
algorithm is used to achieve the pressure-velocity coupling. The spatial
discretization in the governing equations is in the second order. The first-
order implicit scheme is used to discretize the temporal term. In addition,
small integration temporal steps are selected in order to satisfy the Courant-
Friedrichs-Levy (CFL) condition. Thus the Courant number is below 1 for all
simulations undertaken. Note that all the simulations included herein are run
at least throughout $t=400$ non-dimensional time units to ensure enough
characteristics of the quasi-periodic or chaotic flow. Moreover, the data
processing is mainly performed using the whole converged range, discarding the
first 200 units to avoid any transient stage.

Regarding the boundary conditions, a uniform streamwise velocity is imposed at
the inlet boundary, with $\bm{u}=(U_{\infty},0,0)$. For the outlet, an outflow
condition is imposed which implements a Neumann condition for the velocity,
$\partial_{x}\bm{u}=(0,0,0)$, and a Dirichlet condition for the pressure,
$p=0$. A no-slip boundary condition is imposed at the surface of the sphere.
Additionally, slip boundary conditions are set on the cylinder tube walls
therefore any wake-wall interpolations are excluded.

###  II.3 Validation

Grid independence tests are performed at $Re=300$ to select the appropriate
grid size for the numerical study. The relative errors of the typical
parameters using three grids with different node numbers are compared before
conducting extensive simulations.

Table 2: Numerical results at $Re=300$ using three grids with different node numbers (million). $\overline{C_{d}}$, $\overline{C_{l}}$ are the mean values of drag and lift coefficients, respectively. $C_{d}^{{}^{\prime}}$, $C_{l}^{{}^{\prime}}$ are the standard deviations. $St$ is the Strouhal number. The corresponding error deviations in brackets are calculated with respect to the values of the dense grid. Grids | $\overline{C_{d}}$ | $C_{d}^{{}^{\prime}}$ | $\overline{C_{l}}$ | $C_{l}^{{}^{\prime}}$ | $St$
---|---|---|---|---|---
| Coarse
---
3.6
| 0.6653
---
(0.42$\%$)
| 0.001902
---
(2.51$\%$)
| 0.07045
---
(1.37$\%$)
| 0.009839
---
(3.63$\%$)
0.1363
| Moderate
---
5.1
| 0.6631
---
(0.10$\%$)
| 0.001946
---
(0.25$\%$)
| 0.06959
---
(0.13$\%$)
| 0.01012
---
(0.89$\%$)
0.1363
| Dense
---
8.1
0.6625 | 0.001951 | 0.06950 | 0.01021 | 0.1363

Table 2 shows the grid independence analysis results with three sets of grids,
the coarse grid, the moderate grid, and the dense grid. For all three sets,
the grids are refined near the surface of the sphere. The spacing is increased
with a ratio $1.1$, and the first layer is small enough to ensure $y^{+}$ on
the sphere surface is less than $1$. The difference between these three sets
mainly comes from the difference in node numbers on the sphere surface and the
wake region. The mean and standard deviation of the typical parameters of the
sphere wake, e.g., lift and drag coefficients and Strouhal number, are
compared. For this set of variables, there is less than a $1\%$ variation
between the moderate and the dense grid, which means further refinement of the
grid will only affect the results weakly. Thus, the moderate grid confirm the
adequacy of the numerical results and are adopted in the present study.

The results are compared against the available data from other literature to
validate the numerical method.

Table 3: Validation for the numerical method at $Re=300$, with compared with the listed literature. Study | $\overline{C_{d}}$ | $\overline{C_{l}}$ | $St$
---|---|---|---
Present Study | 0.663 | 0.069 | 0.136
Johnson and Patel (1999) | 0.656 | 0.069 | 0.137
Kim, Kim, and Choi (2001) | 0.657 | 0.067 | 0.137
Giacobello, Ooi, and Balachandar (2009) | 0.658 | 0.067 | 0.134
Rajamuni, Thompson, and Hourigan (2018) | 0.665 | 0.070 | 0.137
Chizfahm and Jaiman (2021) | 0.662 | 0.070 | 0.137

Table 3 illustrates the mean drag coefficient, mean lift coefficient, and the
Strouhal number, calculated in this study and other literature for $Re=300$.
The results from different studies are similar to each other. Since the value
of these variables is small and sensitive, it can be concluded that the
present study has a satisfactory match with them. In addition, the result from
the present study is in agreement with the most recent numerical results by
Chizfahm and Jaiman (2021).

Overall, the comparative validations presented here provide confidence that
our computational grid and selected schemes are adequate to resolve the flow
characteristics of the sphere wake appropriately, and sufficient for testing
the reduced-order modeling method.

##  III Trajectory-optimized Cluster-based Network Model

Figure 2: Overview of the presented trajectory-optimized Cluster-based Network
Model framework, exemplified for three flow regimes of the sphere wake
involving periodic, quasi-periodic, and chaotic dynamics. After clustering,
the datasets are grouped into 10 clusters. Trajectory optimization procedure
is used to shift the centroids and build the trajectory support points, which
are defined as nodes of the network. The dynamics is presented by ‘flights’
between the nodes based on two matrixes, including transition probability
$\mathbf{Q}$ and transition time $\mathbf{T}$ inferred from data.

In this section, we propose a novel data-driven reduced-order modeling
methodology, namely the trajectory-optimized Cluster-based Network Model
(tCNM). The tCNM is based on the Cluster-based Network Model (CNM) (Fernex,
Noack, and Semaan, 2021; Li  _et al._ , 2021) which constitutes a
deterministic-stochastic grey-box model for the coherent structure dynamics of
complex flow. The general framework of tCNM exemplified for the sphere wake is
illustrated in figure 2. The first step is the data collection of baseline
flow, where a set of snapshots are collected from experiments or numerical
simulations. The second step is the identification of the cluster distribution
using an unsupervised clustering algorithm to achieve the coarse-graining.
Following is the automated optimization of the cluster centroids and the
building of trajectory support points, which correspond to the nodes of the
network. This step is a critical improvement of tCNM compared with CNM, with
the representation error being reduced. The last step is to characterize the
motion along with the nodes and defines a smoother, more accurate state
propagation trajectory than CNM. We discuss the details of the methodology in
the following subsections. The cluster analysis is described in § III.1. The
characterization of the transition dynamics is illustrated in § III.2. The
trajectory optimization strategy used in the network model is discussed in §
III.3. In the end, § III.4 defines the validation of the model.

###  III.1 Coarse-graining with clustering

In this study, clustering is achieved with the unsupervised $k$-means++
algorithm, which groups the data into clusters such that the intra-cluster
similarity is maximized and the inter-cluster similarity is minimized (Kaiser
_et al._ , 2014). For an ensemble of $M$ statistically representative, time-
resolved snapshots $\bm{u}^{m}$, where $m=1,\ldots,M$, the $k$-means++
algorithm partitions the snapshots into $K$ clusters, defined as
$\mathcal{C}_{k}$, where $k=1,\ldots,K$. The snapshots are equidistantly
sampled, resolving the coherent structure evolution on a time window which is
statistically representative. In particular, the time step determines the
temporal resolution of the snapshots and the accuracy of the transition
dynamics by the cluster-based model. The resolved time range affects the
accuracy of the first and second moments. In our experience, at least 10
periods of the dominant frequency are required for reasonably accurate
statistical moments, and at least 10 snapshots per characteristic period are
needed to resolve the evolution of the dominant coherent structures. In this
study, snapshots are sampled at a time step of $\Delta t=0.2$ over a time
window more than 20 periods, with approximately 40 snapshots per period. The
clustering algorithm can also be applied to the low-dimensional
representations or observables of the full state snapshots.

The algorithm iteratively updates the optimal location of the cluster
centroids according to the data distribution in the state space. Firstly, the
$K$ centroids, namely $\bm{c}_{k}^{0}$ in this study, are chosen from the
snapshots. Compared to the $k$-means algorithm, the $k$-means++ algorithm
optimizes the initial centroids to avoid local optimum and speed up
convergence. Next, each snapshot is affiliated to the closest centroid,
following the cluster affiliation function $f$:

$f\left(\bm{u}^{m}\right)=\arg\underset{k}{\min}\left\|\bm{u}^{m}-\bm{c}_{k}^{0}\right\|_{\Omega},$
(6)

where $\left\|\ \right\|_{\Omega}$ denotes the standard norm on a Hilbert
space of the domain $\Omega$, which is defined as:

$\|\bm{u}\|_{\Omega}:=\sqrt{(\bm{u},\bm{u})_{\Omega}},$ (7)

where the inner product of two square-integrable velocity fields $\bm{u}(x)$
and $\bm{v}(x)$ in domain $\Omega$ reads:

$(\bm{u},\bm{v})_{\Omega}=\int_{\Omega}\bm{u}(x)\cdot\bm{v}(x)\ \mathrm{d}x.$
(8)

After all the snapshots are allocated, the intra-cluster variance $J$ of this
cluster partition is computed as:

$J\left(\bm{c}_{1}^{0},\ldots,\bm{c}_{K}^{0}\right)=\sum_{k=1}^{K}\sum_{\bm{u}^{m}\in\mathcal{C}_{k}}\left\|\bm{u}^{m}-\bm{c}_{k}^{0}\right\|_{\Omega}^{2}.$
(9)

Next, the original cluster centroids $\bm{c}_{k}^{0}$ are updated by averaging
all the snapshots within the cluster $\mathcal{C}_{k}$:

$\bm{c}_{k}^{0}=\frac{1}{n_{k}}\sum\bm{u}^{m},\quad\bm{u}^{m}\in\mathcal{C}_{k},$
(10)

where $n_{k}$ is the number of snapshots in the cluster $\mathcal{C}_{k}$. The
snapshots are allocated to clusters again, according to the new $K$ centroids.
These processes will be repeated until the intra-cluster variance $J$ is
minimized below a specified tolerance. As the initial cluster centroids are
generated randomly, $k$-means++ algorithm is repeated 100 times in this study
to ensure an optimal partition. The averaged centroids $\bm{c}_{k}^{0}$ and
the cluster affiliation $k(\bm{u}^{m})$ of snapshots within the corresponding
clusters $\mathcal{C}_{k}$ are taken as the final results.

The computational cost of clustering algorithms is also considered for large
high-dimensional datasets. An optional pre-processing can be done with a
lossless POD method to conduct data compression (Kaiser  _et al._ , 2014; Li
_et al._ , 2021; Fernex, Noack, and Semaan, 2021; Deng  _et al._ , 2022). The
snapshot $\bm{u}^{m}$ can be exactly expressed by the POD expansions as:

$\bm{u}^{m}(\bm{x})-\bm{u}_{0}(\bm{x})\approx\sum_{i=1}^{N}a_{i}^{m}\bm{u}_{i}(\bm{x}),$
(11)

where $\bm{u}_{0}$ is the constant mean flow field, $\bm{u}_{i}(\bm{x})$ are
the spatial POD modes, with the corresponding mode coefficients
$\bm{a}^{m}=\left[a_{1}^{m},\ldots,a_{N}^{m}\right]$ for temporal evolution. A
complete basis for the decomposition can be expected when $N=M$(Berkooz,
Holmes, and Lumley, 1993). Thus, the distance can be computed with
$\bm{a}^{m}$ instead of the full-dimensional snapshots, which reads as:

$d\left(\bm{u}^{m},\bm{u}^{n}\right)=\left\|\bm{u}^{m}-\bm{u}^{n}\right\|_{\Omega}=\left\|\bm{a}^{m}-\bm{a}^{n}\right\|.$
(12)

The computational load of clustering is significantly reduced with
dimensionality reduction, e.g., an optional POD pre-processing. The clustering
algorithm is applied to the $N$-dimensional state vector $\bm{a}^{m}$ instead
of performing an expensive computation on the full domain with millions of
grid nodes. This computational saving of POD pre-processing can be further
amplified according to the number of iterations required for a $k$-means
algorithm until convergence (Fernex, Noack, and Semaan, 2021). We emphasize
that the POD is only an optional pre-processing to make clustering
computationally feasible and is not necessary for the clustering algorithm
itself. In this work, we compress the time-resolved snapshots with a
_lossless_ POD pre-processing, where the leading $500$ POD modes can resolve
more than $99.9\%$ of the fluctuation energy. Thus, the model performance is
not affected by POD pre-processing. For the small chosen time step, CNM and
tCNM resolve only the non-trivial transitions between the centroids. Thus, the
rapid diffusion of states in the original cluster-based Markov models(Kaiser
_et al._ , 2014) is avoided.

###  III.2 Characterizing the transition dynamics

Since tCNM and CNM share similar transition dynamics, before we detail the
trajectory optimization procedure, CNM is briefly reviewed.

In CNM, the datasets are grouped into clusters by the clustering algorithm,
and cluster centroids are considered as network nodes. Moreover, the motion
between these nodes is characterized by a network model with transition
probabilities and transition times.

The nonlinear dynamics is reconstructed as a linear transition between
centroids. The transitions between clusters rely on the direct transition
matrix, which is typically modeled with a first-order Markov model. The
probability of moving from cluster $\mathcal{C}_{j}$ to $\mathcal{C}_{k}$ is
described by the direct transition matrix $\mathbf{Q}=Q_{kj}\in R^{K\times
K}$, which is deducted from the snapshots as:

$Q_{kj}=\frac{n_{kj}}{n_{j}},$ (13)

where $n_{kj}$ is the number of snapshots moving from cluster
$\mathcal{C}_{j}$ to $\mathcal{C}_{k}$, $n_{j}$ is the number of transitions
departing from $\mathcal{C}_{j}$ regardless of the destination cluster.
Moreover, the direct transition ignores inner-cluster residence probability
with $Q_{kj}=0$ when $k=j$. The clusters are renumbered so that the first
cluster has the highest state probability, the second cluster is chosen with
the highest transition probability from the first one, and so on until all are
sorted (Kaiser  _et al._ , 2014). Note that the matrix is uniformed such that
all elements of this matrix are non-negative, and the elements of each column
sum up to unity.

The transition times are the other transition property, let $t^{n}$ be the
time of the first snapshots to enter and $t^{n+1}$ be the last snapshots to
leave during one individual transition in cluster $\mathcal{C}_{j}$. The
residence time $\tau_{j}^{n}=t^{n+1}-t^{n}$ corresponds to the duration of
staying in cluster $\mathcal{C}_{j}$ at this transition. For the transition
time, we further define the individual transition time $\tau_{kj}^{n}$ from
cluster $\mathcal{C}_{j}$ to $\mathcal{C}_{k}$ as half the residence time of
both clusters, as shown in figure 3.

Figure 3: Transition from cluster $\mathcal{C}_{j}$ to $\mathcal{C}_{k}$.
Different colors correspond to different clusters. $\tau_{j}^{n}$ and
$\tau_{k}^{n}$ are the residence time in cluster $\mathcal{C}_{j}$ and
$\mathcal{C}_{k}$, respectively. $\tau_{kj}^{n}$ is the individual transition
time.

Averaging $\tau_{kj}^{n}$ from all individual transitions, the transition time
is expressed as:

$T_{kj}=\frac{\sum_{n=1}^{n_{kj}}\tau_{kj}^{n}}{n_{kj}},$ (14)

where $n_{kj}$ is the number of individual transitions from $\mathcal{C}_{j}$
to $\mathcal{C}_{k}$ from all the transitions. Similar to the direct
transition matrix $\mathbf{Q}$, we can also build the transition time matrix
$\mathbf{T}=T_{kj}\in R^{K\times K}$.

Based on the transition probabilities and the transition times, we assume a
uniform state propagation between two centroids $\bm{c}_{j}^{0}$ and
$\bm{c}_{k}^{0}$ as:

$\bm{u}^{m}(\bm{x},t)=\alpha_{kj}(t)\,\bm{c}_{k}^{0}(\bm{x})+\left(1-\alpha_{kj}(t)\right)\,\bm{c}_{j}^{0}(\bm{x}),$
(15)

where $\alpha_{kj}(t)=\left(t-t_{j}\right)/T_{kj}$ , $t_{j}$ is the time when
the centroids $\bm{c}_{j}^{0}$ is left.

###  III.3 Trajectory optimization

The proposed improvement relates to the uniform straight-line flight (15) of
CNM. First, the trajectory from centroid $\bm{c}_{j}^{0}$ to $\bm{c}_{k}^{0}$
will be curved, i.e., enter cluster $\mathcal{C}_{k}$ on a different point
than predicted by (15). Second, a typical trajectory may not even be near
$\bm{c}_{k}^{0}$. In the tCNM, we replace the straight-line flights from
cluster $\mathcal{C}_{j}$ to $\mathcal{C}_{k}$ by a curved trajectory through
a shifted centroid $\bm{c}_{j}$, an averaged entry point $\bm{r}_{kj}$, and a
shifted centroid $\bm{c}_{k}$. The shifted centroid $\bm{c}_{k}$ is on the
line connecting $\bm{c}_{k}^{0}$ and the average of all entry points and all
exit points to and from cluster $\mathcal{C}_{k}$. The shift amplitude is
chosen to reduce the representation error of all trajectories through cluster
$\mathcal{C}_{k}$. In other words, the shifted centroid $\bm{c}_{k}$ is
employed for all trajectories regardless of the departure cluster
$\mathcal{C}_{j}$ while the support point $\bm{r}_{kj}$ is chosen for a
specific cluster transition from $\mathcal{C}_{j}$ to $\mathcal{C}_{k}$. We
emphasize a shift in the clustering paradigm: The original centroids optimize
the average squared distance of the snapshots to their nearest centroids. In
contrast, the shifted centroids reduce the distance from the snapshots to the
corresponding approximation of the trajectory. In the following, the
construction is detailed.

For simplicity, we first explain the centroid shift for a single trajectory,
as illustrated in figure 4.

Figure 4: Centroid shift. The cluster affiliation is displayed as in figure 3.
For cluster $\mathcal{C}_{k}$, a straight line is generated based on
$\bm{r}_{k}^{\mathrm{avg}}$ and $\bm{c}_{k}^{0}$, where
$\bm{r}_{k}^{\mathrm{avg}}$ is the average of all entry and exit snapshots.
The averaged centroid $\bm{c}_{k}^{0}$ is shifted to $\bm{c}_{k}$, which is
defined by the closest snapshot from the straight line. The black dashed line
is the propagation with $\bm{c}_{k}^{0}$ of CNM, the red dashed line is the
propagation with the shifted centroid $\bm{c}_{k}$ of tCNM.

A curved trajectory can be found for a transition from cluster
$\mathcal{C}_{j}$ to $\mathcal{C}_{k}$ and then to $\mathcal{C}_{l}$. CNM
uniformly propagates the state through the motion between the three original
centroids $\bm{c}_{j}^{0}$, $\bm{c}_{k}^{0}$ and $\bm{c}_{l}^{0}$, as shown by
the black dashed line. The state propagation shows significant deviation from
the curved trajectory, because the centroid $\bm{c}_{k}^{0}$ from the mean
operator erases lots of critical information and is not suitable to represent
the states in $\mathcal{C}_{k}$. The centroid is shifted to be as close to the
trajectory as possible, ensuring that the error between the state propagation
and the snapshots is minimized. This optimization directly benefits kinematic
state propagation by reducing trajectory variance rather than intra-cluster
variance. We use $\bm{c}_{k}$ to denote the shifted centroids to distinguish
them from the original centroids $\bm{c}_{k}^{0}$.

The centroid shift is inferred from the the snapshots. The entry and exit
points in cluster $\mathcal{C}_{k}$ are defined as:

$\displaystyle\overline{\bm{r}_{k}^{\mathrm{in}}}=\frac{1}{n_{c}}\sum\bm{u}^{m},\quad\bm{u}^{m}\in\mathcal{C}_{k}\
\&\ \bm{u}^{m-1}\notin\mathcal{C}_{k},$ (16)
$\displaystyle\overline{\bm{r}_{k}^{\mathrm{out}}}=\frac{1}{n_{c}}\sum\bm{u}^{m},\quad\bm{u}^{m}\in\mathcal{C}_{k}\
\&\ \bm{u}^{m+1}\notin\mathcal{C}_{k},$

where $n_{c}$ is the number of trajectory segments in $\mathcal{C}_{k}$.
Hence, $\overline{\bm{r}_{k}^{\mathrm{in}}}$ and
$\overline{\bm{r}_{k}^{\mathrm{out}}}$ are the averages of inlet and outlet
snapshots of $\mathcal{C}_{k}$, respectively. The midpoint of two boundary
averages of cluster $\mathcal{C}_{k}$ reads:

$\bm{r}_{k}^{\mathrm{avg}}=\frac{\overline{\bm{r}_{k}^{\mathrm{in}}}+\overline{\bm{r}_{k}^{\mathrm{out}}}}{2}.$
(17)

A straight line $\bm{l}_{k}$ passing through $\bm{r}_{k}^{\mathrm{avg}}$ and
the original centroid $\bm{c}_{k}^{0}$ in Euclidean geometry is given by the
two point form of line equation in Euclidean geometry:

$\bm{l}_{k}(\alpha)=\alpha\,\bm{r}_{k}^{\mathrm{avg}}+\left(1-\alpha\right)\,\bm{c}_{k}^{0},\quad\alpha\in\mathbb{R}.$
(18)

The most representative state for a curved trajectory can be identified at the
junction with the line $\bm{l}_{k}$. For cluster with multiple trajectories,
the shifted centroid is defined as the average of all their junctions. With
the above geometric method, the shifted centroid is close to the extreme in a
curved trajectory and the average in a straight trajectory. The trajectory
optimization by shifting centroids to more representative states is data-
driven and generalizable, achieving the minimum deviation for state motion on
trajectories.

For a trajectory consisting of discretized snapshots, the junction is the
nearest snapshot to the straight line, with the distance read as:

$D_{\bm{l}_{k}}^{m}=\min_{\bm{u}^{m}\in\mathcal{C}_{k}}\left\|\bm{u}^{m}-\bm{l}_{k}(\alpha^{\prime})\right\|_{\Omega},$
(19)

where
$\alpha^{\prime}=(\bm{u}^{m}-\bm{c}_{k}^{0},\bm{r}_{k}^{\text{avg}}-\bm{c}_{k}^{0})_{\Omega}/\left\|\bm{r}_{k}^{\text{avg}}-\bm{c}_{k}^{0}\right\|_{\Omega}^{2}$
consider the foot of the perpendicular on the line $\bm{l}_{k}$ from
$\bm{u}^{m}$. The shifted centroid can hence be obtained by averaging these
chosen snapshots from all the trajectories.

The second trajectory optimization is introducing inter-cluster constraint for
the transitions between two clusters. Since a cluster can have multiple
destination clusters, a refined state motion according to the cluster
transitions will provide more accurate state propagation.

Figure 5: Inter-cluster constraint. The cluster affiliation is displayed as in
figure 3. $\bm{r}_{kj}$, $\bm{r}_{lk}$ are the trajectory support points of
the adjacent clusters. The red dashed line is the direct propagation between
cluster centroids, black dashed line is the optimized propagation considering
the trajectory support points.

Figure 5 illustrates the sketch of inter-cluster constraint for a single
trajectory. Trajectory support points are introduced to ensure the proximity
of an accurate enter to the destination cluster.

The average entry point from cluster $\mathcal{C}_{j}$ to $\mathcal{C}_{k}$ is
noted as trajectory support points $\bm{r}_{kj}$ and defined by

$\bm{r}_{kj}=\left\langle\bm{u}^{m}\right\rangle_{j\mapsto k}.$ (20)

Here, $\left\langle\bm{u}^{m}\right\rangle_{j\mapsto k}$ averages over all
snapshots on segment crossing the boundary between cluster $j$ and cluster
$k$. These snapshots satisfy
$\bm{u}^{m}\in\mathcal{C}_{k}\quad\bm{u}^{m-1}\in\mathcal{C}_{j}$ or
$\bm{u}^{m}\in\mathcal{C}_{j}\quad\bm{u}^{m+1}\in\mathcal{C}_{k}$. The model
trajectory connects $\bm{c}_{j}$, the entry point $\bm{r}_{kj}$ and
$\bm{c}_{k}$ and can be expected to be more accurate than the original
straight-line segment.

In CNM, the transition between cluster $\mathcal{C}_{j}$ and $\mathcal{C}_{k}$
is based on a uniforn motion between two centroids. In tCNM, the model
trajectory is a spline through centroids and ‘border-crossing’ trajectory
support points. The times from centroid $\bm{c}_{j}$ to $\bm{r}_{kj}$ from
$\bm{r}_{kj}$ to $\bm{c}_{k}$ are estimated to be $T_{kj}^{-}$ and
$T_{kj}^{+}$ respectively. In other words, the transition time is decomposed
in two parts, $T_{kj}=T_{kj}^{-}+T_{kj}^{+}$.

An efficient realization of CNM requires POD of the snapshot data. The
clustering adds a small computational cost in comparison to building POD,
about 10% for this study. The computational load of the transition model,
i.e., building the matrixes $\mathbf{Q}$ and $\mathbf{T}$, and integrating the
model is small even compared to the cost of clustering. Details on the
computational load of these algorithms are provided in the previous study by
Fernex, Noack, and Semaan (2021). In this study, the computational costs of
CNM and tCNM add 11% and 30%, respectively, to the cost of POD, taken as 100%.
tCNM is more costly than CNM because of the matrix operations to identify
$\bm{c}_{k}$ as well as to get $\bm{r}_{kj}$.

###  III.4 Validation

There are usually two kinds of prediction errors for cluster based models:
from the spatial perspective, the representation error due to the insufficient
representation by cluster centroids, and from the temporal perspective, the
transition error due to the inaccurate reconstruction of complex transition
relationships. Considering CNM, for complex phase-space trajectories, more
accuracy will be achieved by the high-order Markov chains (Fernex, Noack, and
Semaan, 2021). The high-order CNM requires more previous states with time
delays to decide the destination of the current state, enabling finer
differentiation of transition dynamics and lower transition error. However,
our focus is on decreasing representation error with a trajectory-optimized
extension, thus, the traditional first-order model is chosen in this study.
The trajectory optimization by centroid shift and inter-cluster constraint can
be seen as a procedure to reduce the representation error without producing
extra transition error. It is worth noting that there will be no coupling
relationship between these two kinds of errors once the number of clusters is
determined.

The chosen number of clusters $K$ significantly affects the prediction error.
The network model with a small number of clusters can capture the dominant
behavior of transition dynamics. Fewer transitions between clusters will be
accounted for, enabling easier interpretation of the results with physical
insights (Fernex, Noack, and Semaan, 2021). Therefore the transition error is
relatively small. However, a larger representation error is induced due to
representing the dynamics by fewer clusters. The details of snapshots are
eliminated as the centroids average the snapshots in each cluster. In
contrast, a large $K$ will allow CNM to model the transition dynamics with
more detail, which means a more accurate representation of snapshots. However,
considering more clusters, the transition relationship will be more complex.
Especially with an excessively large $K$, the centroids are too dense so that
the nearly identical motions are clustered separately into adjacent
trajectories, and even random switching between trajectories is introduced. In
this case, the resulting CNM exhibits non-physical transitions, leading to a
significant increase in transition error. Thus, we expect a sweet spot with
optimal prediction error based on good representation error and an accurate
prediction for the number of clusters. Nair  _et al._ (2019) determined an
optimal choice of cluster number by the F-test considering the ratio of inter-
cluster variance to the total variance. A similar balancing strategy, like the
elbow method (Kaiser  _et al._ , 2014), can also be used to get an appropriate
cluster number. In this study, the snapshots are coarse-grained into $K=10$
clusters, which is consistent with many previous studies(Kaiser  _et al._ ,
2014; Nair  _et al._ , 2019; Li  _et al._ , 2021; Li and Tan, 2020; Tan, Li,
and Noack, 2020). This value is in good agreement with the F-test and elbow
method, and is large enough to resolve the main transition mechanism.

To statistically evaluate the performance of the reduced-order models,
including the representation error and the transition error. Several metrics
are presented: the autocorrelation function, the cluster probability vector
and the function of representation error.

The autocorrelation function avoids the problem of comparing two trajectories
directly with finite prediction horizons due to phase mismatch (Fernex, Noack,
and Semaan, 2021). Especially for chaotic flow, where a slight change in
initial conditions leads to completely different trajectories, thus, the
direct comparison of time series is often pointless. Note that for pure
periodic flow, the direct comparison of time series often makes sense since
there is no phase mismatch. The function is defined as(Protas, Noack, and
Östh, 2015):

$R(\tau)=\frac{1}{T-\tau}\int_{0}^{T-\tau}(\bm{u}(\bm{x},t)\cdot\bm{u}(\bm{x},t+\tau))_{\Omega}\
\mathrm{d}t,\tau\in[0,T].$ (21)

The cluster probability vector is defined as:

$\bm{p}=\left[p_{1},\ldots,p_{k}\right],$ (22)

where $p_{k}$ represents the probability of being in cluster
$\mathcal{C}_{k}$, which is the cumulative residence time normalized by the
simulation time horizon, as:

$p_{k}=\frac{\sum\tau_{k}^{n}}{T}.$ (23)

The vector $\bm{p}$ indicates whether the predicted trajectories could
populate the phase space similar to the original data. The comparison, to some
extent, reveals the transition error of the modeling method.

The representation error of the reconstructed dynamics is defined as:

$E_{r}=\frac{1}{M}\sum_{m=1}^{M}D_{\mathcal{T}}^{m},$ (24)

where $D_{\mathcal{T}}^{m}$ is defined as the minimal distance from the
snapshot $\bm{u}^{m}$ to the states on the reconstructed trajectory
$\mathcal{T}$:

$D_{\mathcal{T}}^{m}=\min_{\bm{u}^{n}\in\mathcal{T}}\left\|\bm{u}^{m}-\bm{u}^{n}\right\|_{\Omega}.$
(25)

##  IV Flow features at different Reynolds number

Flow past a sphere is a prototype wake of bluff bodies commonly encountered in
many engineering applications, for instance, the design of drones, air taxis,
or micro-robots. Due to the geometric simplicity and symmetry, the sphere wake
has been well investigated for its three-dimensional flow characteristics,
which exhibits different patterns at different $Re$, such as steady, periodic,
quasi-periodic, and chaotic flow. These different kinds of baseline flow are
ideal for testing the performance of reduced-order modeling. We can get the
dataset of general flow dynamics without changing the numerical configuration
but simply adjusting $Re$.

As $Re$ increases, the wake behind a sphere undergoes a series of bifurcations
and transitions on its way to chaos. The flow separation around a sphere is
widely accepted to occur at $Re\approx 20$ (Chrust, Goujon-Durand, and
Wesfreid, 2013).

Figure 6: Representative vortex structures and temporal evolution of $C_{d}$,
$C_{l}$ at different $Re$. The vortices are visualized using iso-surfaces of
the $Q$-criterion ($Q$ = 0.001) color-coded with non-dimensional velocity
$U/U_{\infty}$, and the visualization of all baseline flows uses the same
legend and view. (a), (b) Planar symmetry steady flow at $Re=250$, with
$C_{d}$ and $C_{l}$ being constants. (c), (d) Periodic flow at $Re=300$, with
periodic shedding of hairpin vortices and a symmetry plane ($xy$ plane). (e),
(f) Quasi-periodic flow at $Re=350$, with a newly appeared lower frequency
modulates the shedding. (g), (h) Chaotic flow at $Re=450$, with the irregular
shedding and the azimuthal rotation of the separation point, as characterized
by chaotic features of $C_{d}$ and $C_{l}$.

The resulting wake remains axisymmetric until $Re\approx 210$, where a
pitchfork bifurcation occurs, which breaks the axial symmetry and leads to a
steady non-axisymmetric wake consisting of a pair of steady streamwise
vortices that are symmetrical about a longitudinal plane (Fabre, Auguste, and
Magnaudet, 2008), like in figure 6 (a). For better observation, a rotation
along the $x$-axis is applied to keep the symmetry plane (if present)
coincident with the $xy$ plane. A Hopf bifurcation is observed at $Re\approx
270$, corresponding to the appearance of vortex streets or the shedding of
hairpin vortices, like in figure 6 (c). The symmetry of the wake with respect
to the longitudinal plane still preserves (Fabre, Auguste, and Magnaudet,
2008). At higher $Re\approx 350$, the planar symmetry of vortex shedding
starts to lose its trace, and the hairpin vortices exhibits differences among
each shedding period (Zhao, 2021). A lower frequency appears and modulates the
shedding, e.g., the flow exhibits a longer periodicity of about seven shedding
periods at $Re=350$, which results from a secondary Hopf bifurcation, as shown
in figure 6 (e). The successive bifurcations eventually lead to a fully
chaotic flow at $Re>420$ (Pan, Zhang, and Ni, 2018), revealing more complex
vortical structures with three-dimensional properties, as shown in figure 6
(g).

In this study, the four representative baseline flow regimes mentioned above
are investigated, including the planar symmetry steady state at $Re=250$, the
periodic, quasi-periodic, and chaotic flow at $Re=300$, $350$, and $450$,
respectively. The three-dimensional vortex structures and $C_{d}$, $C_{l}$ of
the sphere are shown in the right column of figure 6. For the planar symmetry
steady flow, the streamwise vortices remain stationary, with two threads of
opposite axial vorticity, and $C_{d}$ and $C_{l}$ are all constants. Since the
vortex intensity of the corresponding two vortices is not the same, $C_{l}$
would not be zero but is one order of magnitude smaller than $C_{d}$. For
periodic flow, the shedding of hairpin vortices dominates the whole flow
field, which remains identical among each shedding period with $C_{d}$ being
periodic, while the vortex intensity of the corresponding vortex pair still
exhibits a slight difference from each other. Thus, the average $C_{l}$ would
not be zero, but its amplitude will remain periodic. For quasi-periodic flow,
the vortex intensity of the corresponding vortex pair exhibits a larger
difference. The similarity of the hairpin vortices among each shedding period
is lost, and a newly appeared lower frequency modulates the shedding, with
$C_{d}$ and $C_{l}$ becoming quasi-periodic. For chaotic flow, the shedding of
hairpin vortices becomes fully irregular, and displays different orientations,
with an azimuthal rotation of the separation point at the rear of the
sphere(Lorite-Díez and Jiménez-González, 2020). Meanwhile, the wake exhibits
obvious lateral oscillations and secondary vortices(Eshbal  _et al._ , 2019),
as characterized by chaotic features of $C_{d}$ and $C_{l}$.

In the following section, tCNM will be applied to the nonlinear dynamics of
the periodic, quasi-periodic, and chaotic regimes.

##  V Modeling the sphere wake

In this section, tCNM is applied to periodic (§ V.1), quasi-periodic (§ V.2)
and chaotic (§ V.3) flow regimes. We evaluate the performance of tCNM by
comparing the temporal evolution, autocorrelation functions, and
representative vortex distribution profiles against CNM.

At higher $Re$, coherent structures can still dominate the flow and be
resolved by centroids. The small-scale stochastic component is filtered out in
the clustering process(Fernex, Noack, and Semaan, 2021). In contrast, cluster-
based modeling of transient flows from the steady solution to post-transient
dynamics requires careful preparation of the simulation data and thorough
post-processing(Deng  _et al._ , 2022).

###  V.1 Periodic flow

Periodic unsteady flow is characterized by a single frequency and spatial-
temporal symmetry. The flow regime involves only a Hopf bifurcation, and the
dynamic reconstruction is simple to implement. The dataset is firstly
compressed with a lossless POD pre-processing. The temporal evolution of the
POD mode coefficients is visualized in figure 7 (a). Only the first, third and
fifth POD mode coefficients $a_{1}$, $a_{3}$ and $a_{5}$ are presented, as the
POD modes usually appear in pairs, and the leading modes can be enough to
represent most of the flow characteristics. The first pair of POD modes
$a_{1}$, $a_{2}$ are an order of magnitude larger than all the other
coefficients, which is common in periodic flows because the coherent structure
dominates and holds most of the kinematic energy. The phase diagram based on
$a_{1}$, $a_{2}$ exhibits a standard cycle because of the first harmonic, as
viewed in figure 7 (b). Higher-order modes present higher harmonics and even a
mixture of different frequencies, which usually do not affect the macro-scale
flow characteristics but have an impact on the micro-scale characteristics.

Figure 7: Temporal evolution (a) and phase diagram (b) of the leading POD
coefficients for periodic flow at $Re=300$.

The clustering is performed on the vectors composed of the leading 500 POD
coefficients. The assignment of each snapshot to the clusters is shown in
figure 8. The cluster indexes are reordered as § III.2 illustrates. For
periodic flow, the vortex shedding is resolved by the ten clusters, and a
clear periodicity can be seen from the cluster affiliation. The whole shedding
period is divided by ten different vortex patterns corresponding to different
clusters, these vortex patterns are further represented by the cluster
centroids of CNM and the shifted centroids by tCNM, which are severed as the
nodes of the network. The almost uniform residence time in each cluster shows
a strict periodicity and indicates that the size of each cluster is relatively
homogeneous.

Figure 8: Cluster affiliation for periodic flow at $Re=300$. Figure 9: Three-
dimensional proximity map for periodic flow at $Re=300$. Gray balls are the
original CNM centroids, orange balls are the shifted tCNM centroids, solid
circle are the trajectory support points. $a_{1}$, $a_{2}$, $a_{3}$ are the
corresponding features of the three-dimensional space by MDS, and a projection
are conducted into the $a_{1}a_{2}$ plane. The difference between the shifted
centroids and the original centroids is comparatively small for the periodic
regime with respect to the other two cases in the following study. Figure 10:
Transition matrixes for periodic flow at $Re=300$. (a) Direct transition
matrix $\mathbf{Q}$. (b) Transition time matrix $\mathbf{T}$. Figure 11:
Cluster probability distribution for periodic flow at $Re=300$. Only tCNM is
exhibited as CNM and tCNM share the same probability distribution.

Figure 12: Temporal evolution of POD mode coefficients for periodic flow at
$Re=300$. (a) Coefficient $a_{1}$. (b) Coefficient $a_{3}$. (c) Coefficient
$a_{5}$. The black solid curve is CFD data, the green dashed curve is the
reconstructed temporal dynamics by CNM, the red dotted curve is the
reconstructed temporal dynamics by tCNM. Figure 13: Vortex distribution in one
typical vortex shedding period of $xy$ plane at $z=0$ (a) and $xz$ plane at
$y=0$ (b) for periodic flow at $Re=300$. For the $xy$ plane, the distribution
is color-coded with $z$ vorticity, while for the $xz$ plane with $y$
vorticity. The first row is the CFD data, the second row is the reconstructed
dynamics by CNM, and the third row is by tCNM.

The high-dimensional dataset and their cluster affiliation can be easily
interpreted by multi-dimensional scaling (MDS), which can project the high-
dimensional feature space into a low-dimensional one and then build a
proximity map (Li and Tan, 2020; Tan, Li, and Noack, 2020). In this study, we
project all the POD coefficient vectors of the snapshots and the corresponding
nodes, including the original cluster centroids $\bm{c}_{k}^{0}$ of CNM, the
shifted cluster centroids $\bm{c}_{k}$ and the trajectory support points
$\bm{r}_{kj}$ of tCNM, into a three-dimensional feature, as shown in figure 9.
The gray balls represent $\bm{c}_{k}^{0}$, the orange balls represent
$\bm{c}_{k}$, and the solid circle represent $\bm{r}_{kj}$. The distance
between any two dots or balls in the feature space simply represents their
similarity level. These snapshots, represented by dots and colored by their
cluster affiliation, form a clear orbit. A clear difference between CNM and
tCNM in the building of nodes can be seen from the proximity map. For typical
clusters, $\bm{c}_{k}^{0}$ usually exhibits certain deviation from the
original trajectory, while $\bm{c}_{k}$ and the trajectory support points
$\bm{r}_{kj}$ are exactly on the trajectory, making the propagation optimized
substantially. Thus, the propagation of tCNM based on both $\bm{c}_{k}$ and
$\bm{r}_{kj}$ is more accurate with less prediction error, as compared to CNM
based on $\bm{c}_{k}^{0}$.

The direct transition matrix $\mathbf{Q}$ and the transition time matrix
$\mathbf{T}$ can reveal the temporal evolution of all the possible
transitions. For periodic flow, the flow characteristics are fully resolved,
as these matrixes briefly illustrate the limit cycle in one shedding period,
as shown in figure 10. The transition matrix shows a strictly periodic
transition route, which indicates that the transition dynamics is
deterministic and periodic, as also revealed by the cluster affiliation shown
in figure 8. Each entry of the direct transition matrix $\mathbf{Q}$ is unity,
from cluster $\mathcal{C}_{k}=1$ to $\mathcal{C}_{k}=2$ then gradually to
$\mathcal{C}_{k}=10$ and finally return to $\mathcal{C}_{k}=1$. The transition
times of all the transitions are similar, revealing that the idling time of
each cluster is also similar. Thus, the population probability of these
clusters will be uniform.

The temporal dynamics of CNM and tCNM are built based on the transition
matrixes $\mathbf{Q}$ and $\mathbf{T}$, as shown in figure 11. We choose the
same initial cluster $\mathcal{C}_{k}$ as the CFD data and integrate the
dynamics over 200 time units. As the two models share the same transition
matrixes, the same temporal evolution are used to ensure a clear comparison
for their representation error. Thus, the cluster population distribution of
CNM and tCNM will be the same, we further compare that with the CFD data, as
shown in figure 11. For periodic flow, the probability distribution of tCNM
(or CNM) is similar to the CFD data, which means the baseline flow is
probabilistically well reconstructed with minor transition error. The
distribution of each cluster is nearly $10\%$, which is consistent with the
cluster affiliation in figure 8.

Figure 12 shows the temporal evolution of POD mode coefficients $a_{1}$,
$a_{3}$ and $a_{5}$. The black solid curve represents the original CFD data,
the green dashed curve represents CNM and the red dotted curve represents
tCNM. Following Fernex, Noack, and Semaan (2021) and Li  _et al._ (2021), the
temporal evolution is smoothed by spline interpolation. For $a_{1}$, the two
modes all track well the phase, while as for the amplitude, tCNM performs
slightly better near the extreme, practically the same as the CFD data. For
the other two modes associated with higher harmonics, the information around
the extremes cannot be fully resolved by the averaged centroids
$\bm{c}_{k}^{0}$ in CNM, leading to some distortion in the reconstructed
dynamics. However, the distortion can be easily avoided in tCNM using the
optimized centroids.

The evolution of vortical structures in one vortex shedding period is
illustrated through the slices of 3D flow field, which brings a more intuitive
view of the difference between the two models. Here, the vortex distribution
on the middle $xy$ plane and $xz$ plane are visualized in figure 13. The top
row is the CFD data, the middle row is the reconstructed CNM results, and the
bottom row is the tCNM results. The results of CNM and tCNM are equidistantly
sampled with a time step $T/4$ from $t_{0}=100$, where $T$ is one shedding
period detected from the CFD data. We do these in a similar way for quasi-
periodic and chaotic flow. The corresponding shedding period of the CFD data
is chosen for comparison. Since there is no phase mismatch in periodic flow,
the CFD data is also chosen from $t_{0}=100$ with the same time step. The two
models perform well on the first pair of modes for periodic flow, which catch
most of the flow characteristics. Therefore, the contour of the two models is
all extremely close to the CFD data for most of the snapshots, exhibiting
strong applicability of cluster-based reduced-order modeling methods for
periodic flow. tCNM shows its advantages for the snapshots of $t=T/2$, where
the snapshots are at the critical flow state of one shedding period, or in
other words, at the moment when one vortex is about to separate from the
sphere. CNM lost certain micro-scale vortex structures due to its inaccuracy
in high-order modes, while tCNM can better represent the original data and
retain most of the details of the vortex. As this moment is the critical time
point of the shedding period, the importance of accurate reconstruction near
this moment is self-evident. Thus, the improved strategy with trajectory
optimization in tCNM can be essential. In addition, for the case where both
models can reconstruct the original data accurately, tCNM performs slightly
better, resolving more details than CNM due to the impact of inter-cluster
constraint, which provides a refined transition reconstruction on the
trajectory.

###  V.2 Quasi-periodic flow

Figure 14: Same as figure 7, but for quasi-periodic flow at $Re=350$. (a)
Phase diagram. (b) The leading POD coefficients.

For quasi-periodic flow, a newly appeared lower frequency modulates the
shedding of hairpin vortices. Figure 14 illustrates the temporal evolution of
the POD mode coefficients $a_{1}$, $a_{3}$, $a_{5}$ and the phase diagram of
$a_{1}$, $a_{2}$. The mode coefficients exhibit a long periodicity in the time
evolution with around seven shedding periods and the phase diagram shows a
torus. Moreover, the similar order of magnitude of the three coefficients
indicates that the leading pair of modes cannot completely resolve the hairpin
vortices in the quasi-periodic regime. These differences bring challenges to
the accuracy of reduced-order modeling approaches, especially for the
representation error of cluster-based methods.

Figure 15: Same as figure 8, but for quasi-periodic flow at $Re=350$.

The assignment of each snapshot to the clusters is shown in figure 15. A
longer period is observed, which contains two kinds of transition evolution.
In the first part ($t\approx 0$ to $25$), the ten clusters are visited
sequentially, and their residence times are almost uniform. In the second part
($t\approx 25$ to $60$), the transition dynamics is more complex, regularly
returning to the visited clusters. This long-period dynamic includes around
seven shedding periods, and one shedding period is resolved by a combination
of only three or four clusters rather than all the clusters. The different
periodic sheddings within a complete loop are resolved by different sorts of
clusters. Consequently, more snapshots are involved in one cluster than
periodic flow, which means the representation error will be larger for CNM.

Figure 16: Same as figure 9, but for quasi-periodic flow at $Re=350$.

Figure 16 shows the three-dimensional view of all the snapshots and
propagation nodes of quasi-periodic flow. The distribution of snapshots
resembles the side surface of a cone. Different parts of the side surface
correspond to different clusters depicted in different colors. The dynamics
appear to be driven by two physical phenomena: a cyclic behavior and a quasi-
periodic component that forces the limit cycle to experience cycle-to-cycle
variations (Cao  _et al._ , 2014). The averaged cluster centroids
$\bm{c}_{k}^{0}$ (gray balls) are inside the surface, while the shifted
cluster centroids $\bm{c}_{k}$ (orange balls) and the trajectory support
points $\bm{r}_{kj}$ (solid circle) are right on the surface. This indicates
that the averaged centroids are out of the data manifold, while the optimized
centroids are closer to the data, thus having a smaller representation error.

Figure 17: Same as figure 10, but for quasi-periodic flow at $Re=350$. (a)
Direct transition matrix $\mathbf{Q}$. (b) Transition time matrix
$\mathbf{T}$.

The direct transition matrix $\mathbf{Q}$ and the transition time matrix
$\mathbf{T}$ of quasi-periodic flow are illustrated in figure 17. An
inspection of $\mathbf{Q}$ reveals that most of the clusters have at least two
destination clusters. However, the cluster transition with higher
probabilities form a main loop from cluster $\mathcal{C}_{k}=1$ to
$\mathcal{C}_{k}=10$, which indicates a relatively deterministic transition
appeared in figure 15. In addition, there are some transitions with minor
probabilities, which will lead to random transitions appearing in the network
model and introduce the additional transitions errors. The high-order network
model is considered as an effective way to solve this problem (Fernex, Noack,
and Semaan, 2021). From $\mathbf{T}$, the transition time of each cluster in
the main loop is different but presents a smooth variation. That is, the
residence time in the adjacent clusters along the trajectory varies
continuously, as well as the corresponding vortex structure.

Figure 18: Same as figure 11, but for quasi-periodic flow at $Re=350$.

The probability distribution of quasi-periodic flow is well reconstructed by
tCNM and CNM, as shown in figure 18. The probability distribution of each
cluster is not uniform, which is reasonably larger for the cluster with higher
transition probability and longer transition time. The difference between the
distribution of the two models and the CFD data become larger, as compared to
the periodic case. It is reasonable from the following two aspects. One is the
transition error becoming larger for quasi-periodic flow. The other is the
assumed constant transition time $T_{kj}$ for all trajectories from cluster
$\mathcal{C}_{j}$ to $\mathcal{C}_{k}$ is a crude assumption due to the
coarse-graining of the clustering. The transition times can vary by a large
factor and can thus give rise to significant systematic error (Li  _et al._ ,
2021).

Figure 19: Same as figure 12, but for quasi-periodic flow at $Re=350$. (a)
Coefficient $a_{1}$. (b) Coefficient $a_{3}$. (c) Coefficient $a_{5}$.

Figure 20: Autocorrelation function of POD mode amplitudes for quasi-periodic
flow at $Re=350$. (a) Coefficient $a_{1}$. (b) Coefficient $a_{3}$. (c)
Coefficient $a_{5}$. (d) Coefficient $a_{7}$. The black solid curve is CFD
data, the green dashed curve is the CNM data, the red dotted curve is the tCNM
data. Figure 21: Same as figure 13, but for quasi-periodic flow at $Re=350$.
(a) $xy$ plane at $z=0$. (b) $xz$ plane at $y=0$.

The temporal evolution of POD mode coefficients $a_{1}$, $a_{3}$ and $a_{5}$
is illustrated in figure 19. For CNM, as the vortex shedding period is highly
coarse grained, the loss of the critical flow states by the averaged centroids
correspondingly gets more visible than the periodic flow. Therefore the
amplitude of the restored coefficients is much smaller. Despite the phase
mismatch due to the complex transition relationship of several clusters, tCNM
is shown to be more accurate for restoring the critical dynamics of the quasi-
periodic flow

Considering the impact of phase mismatch, the autocorrelation function of POD
mode coefficients $a_{1}$, $a_{3}$, $a_{5}$ and $a_{7}$ is calculated and
shown by figure 20. Similar to the temporal evolution, the black solid curve
represents the original CFD data, the green dashed curve represents CNM and
the red dotted curve represents tCNM. Note that a small range of $t$ by
definition also makes sense for the comparison. To show the difference between
the two modes more clearly, we limit the range of $t$ to 20 time units. For
larger $t$, situations are similar. The fluctuation level of CNM is much lower
than tCNM, due to the information loss. Moreover, the higher-order mode is,
the lower the fluctuation level becomes, which is consistent with the periodic
flow. The spatial-temporal resolution of cluster-based models depends on the
number of clusters. A handful of clusters can well resolve the dominant
frequency, but a significantly larger number of clusters is needed to resolve
the higher harmonics with lower amplitudes. This is apparent in figure 20,
where the dominant dynamics of leading modes is well captured by only ten
centroids, but the high-frequency dynamics of higher-order modes is
compromised.

The vortex distribution on the middle $xy$ and $xz$ plane are illustrated in
figure 21. As expected, the most visible difference still occurs at the
snapshots of $t=T/2$, as the contour of CNM is more like an averaged flow
state rather than an instantaneous one. Several coherent structures disappear,
and slight deformation can be seen on the vortex pair. While for tCNM, many of
the details are still well retained. As the propagation of CNM is generally
based on only three cluster centroids for one vortex shedding period for
quasi-periodic flow, extensive information loss of one centroid will distort
the whole trajectory. Thus, for the other snapshots, the inaccuracy of CNM can
also be seen. Note that the first and last snapshots of tCNM also exhibit
certain differences with CFD data. We infer it is mainly due to the phase
mismatch of the nearby shedding periods.

###  V.3 Chaotic flow

For chaotic flow, the vortex shedding becomes irregular and fully three-
dimensional chaotic, which is explained as the consequence of the azimuthal
rotation of the separation point and lateral oscillations of the vortex
shedding. Figure 22 shows the temporal evolution of the POD mode coefficients
$a_{1}$, $a_{3}$, $a_{5}$ and the phase diagram of $a_{1}$, $a_{2}$ for the
fully chaotic three-dimensional flow at $Re=450$. The three coefficients are
in the same order of magnitude. However, the temporal evolution becomes fully
irregular for each coefficient, and no obvious periodic pattern can be found.
The phase diagram also manifests chaotic characteristics. As the temporal
evolution of these coefficients gets more complex in chaotic cases, the inter-
cluster constraint of tCNM will exhibit more obvious effects in the dynamic
representation.

Figure 22: Same as figure 7, but for chaotic flow at $Re=450$. (a) Phase
diagram. (b) The leading POD coefficients. Figure 23: Same as figure 8, but
for chaotic flow at $Re=450$.

Figure 23 illustrates the assignment of each snapshot to the clusters. It is
hard to find stable cycling from the cluster affiliation for chaotic flow. One
vortex shedding period is typically resolved by the first six clusters, which
is unstable and quickly jumps into the rest four clusters randomly. The
chaotic shedding still retains some characteristics of periodic flow but the
incoherent fluctuations and stochastic dynamics dominates. However, when
clustering the chaotic flow, slightly different vortex structures derived from
lateral oscillations and secondary vortexes are grouped into the same
clusters, resulting in an increase in the intra-cluster variance and a
decrease in the representation ability of the cluster centroids. This
introduces a more significant extra representation error. This kind of
representation error is relatively small in periodic or quasi-periodic flow,
as the coherent structures dominate in these flows.

Figure 24: Same as figure 9, but for chaotic flow at $Re=450$. Figure 25: Same
as figure 10, but for chaotic flow at $Re=450$. (a) Direct transition matrix
$\mathbf{Q}$. (b) Transition time matrix $\mathbf{T}$. Figure 26: Same as
figure 11, but for chaotic flow at $Re=450$.

Figure 27: Same as figure 12, but for chaotic flow at $Re=450$. (a)
Coefficient $a_{1}$. (b) Coefficient $a_{3}$. (c) Coefficient $a_{5}$.

Figure 28: Same as figure 20, but for chaotic flow at $Re=450$. (a)
Coefficient $a_{1}$. (b) Coefficient $a_{3}$. (c) Coefficient $a_{5}$. (d)
Coefficient $a_{7}$. Figure 29: Same as figure 13, but for chaotic flow at
$Re=450$. (a) $xy$ plane at $z=0$. (b) $xz$ plane at $y=0$.

A three-dimensional view of all the snapshots and propagation nodes is shown
in figure 24. The distribution of snapshots in the state space is rather
scattered and irregular, but a ring can be seen in the top view. The
clustering algorithm divides them into different clusters, like cutting a
chaotic cloud area into different colored pieces. CNM and tCNM try to figure
out a clear pattern of trajectory motion statistically. From the top view, the
averaged centroids $\bm{c}_{k}^{0}$ locates deep in each piece, while the
shifted cluster centroids $\bm{c}_{k}$ and the trajectory support points
$\bm{r}_{kj}$ are near the physical trajectories. Most of the original
centroids $\bm{c}_{k}^{0}$ are away from the shifted centroids $\bm{c}_{k}$,
and the trajectory support points will further help to refine the cluster
transitions, which means that tCNM can significantly reduce the representation
error of CNM for the chaotic flow.

The transition matrix $\mathbf{Q}$ and the transition time matrix $\mathbf{T}$
of chaotic flow are illustrated in figure 25. The transition matrix do not
exist a main loop traveling through all the clusters. The cluster transitions
with high probabilities form a cycling transition from cluster
$\mathcal{C}_{k}=1$ to $\mathcal{C}_{k}=6$, which corresponds to a vortex
shedding period resolved by six clusters, appearing most frequently in this
flow regime. This also reveals the advantage of the cluster-based reduced-
order modeling method in resolving the dynamics of complex dynamic systems, as
a clear limit cycle can be recognized even for chaotic flow. The transitions
from the rest four clusters show to be highly random. The cycling block from
cluster $\mathcal{C}_{k}=1$ to $\mathcal{C}_{k}=6$ can jump into the random
block from cluster $\mathcal{C}_{k}=7$ to $\mathcal{C}_{k}=10$ through
$\mathcal{C}_{k}=1,4,5,6$, Conversely, the random clusters can enter the
circular orbit from $\mathcal{C}_{k}=7,8,10$. $\mathcal{C}_{k}=6,8$ are the
key transient points for choosing the periodic and stochastic routes,
corresponding to the different shedding features. From the time matrix
$\mathbf{T}$, we can also find the dissimilarity of transition times in
different transitions, as in the quasi-periodic flow. However, no smooth
variation can be found for the random clusters. For the clusters affiliated
with the limit cycle, the transition times are similar. For the random
clusters, some transition times are much longer, indicating dense distribution
of data in these clusters.

The cluster probability distribution is shown in figure 26. The clusters
affiliated with the limit cycle have a relatively higher distribution due to
their higher transition possibility, and a uniform distribution due to their
analogous transition times. Even the transition times are longer for the other
clusters, they still have a lower distribution due to the lower transition
probability. Note that the lower distribution for cluster $\mathcal{C}_{k}=7$
attributed to its low transition probability as destination clusters from
cluster $\mathcal{C}_{k}=6$. Despite the dataset being probabilistically well
reconstructed, certain differences can be found between the distribution of
two models and the CFD data, and even more significant than the quasi-periodic
flow. As the distribution tendency is still similar and these differences only
affect the transition error, we deem these differences can have little impact
in this study to compare two models with representation error.

The temporal evolution of $a_{1}$, $a_{3}$ and $a_{5}$ is illustrated in
figure 27, the autocorrelation function of $a_{1}$, $a_{3}$, $a_{5}$ and
$a_{7}$ is shown by figure 28. From the temporal evolution, it is clear that
the amplitude of the restored coefficients by tCNM is much larger than that of
CNM due to the lesser information loss. Most of the amplitudes resolved by
tCNM show a higher accuracy, similarly with the autocorrelation function.
Additionally, since the shifted centroids $\bm{c}_{k}$ are the averages of all
the representative flow states, the trajectory of tCNM through one cluster can
be seen as the average of all the optimized trajectories, which corresponds to
all the slightly dissimilar vortex structures in this cluster. Thus, the extra
representation error can be substantially moderated by tCNM, as demonstrated
in the autocorrelation function of $a_{1}$ in figure 28. For the higher-order
mode coefficients, the fluctuation level of the autocorrelation function for
both tCNM and CNM is lower than the CFD data, which is mainly caused by the
transition error.

The vortex distribution on the middle $xy$ and $xz$ plane are illustrated in
figure 29. For chaotic flow, as the separation point is rotating around the
sphere and the symmetric plane is lost, we did not do the rotation to the
numerical domain. The contours of CNM are like the averaged flow state nearly
in every snapshot, and the coherent structures and secondary vortexes cannot
be easily recognized compared with tCNM. Meanwhile, slight deformation also
occurs on the vortex pair. For tCNM, the largest differences compared with the
CFD data still occurs on the first and last snapshots, the reason is the same
as the quasi-periodic flow. While for the other snapshots, most of the
delicate vortex structures are well preserved, which is explained as the
consequence of the combined improvement by centroid shift and inter-cluster
constraint.

Figure 30: Representation error of three cluster-based models for the three
flow regimes and the corresponding average. All the levels are normalized with
respect to the representation error of CROM.

In figure 30, to illustrate the performance improvement of tCNM, a statistical
comparison based on the representation error between three cluster-based
models is conducted for all the three flow regimes. The representation error
evaluates the reconstructed dynamics to the raw snapshots. For tCNM and CNM,
this error is defined as the average of the distances from the snapshot to the
reconstructed trajectory. Whereas for CROM, the dynamics is expressed as
direct jumps between discrete cluster centroids. Hence this error is defined
as the average of the distances from the snapshots to their nearest centroids.
The results of the three models are normalized by the error of CROM for
intuitive comparison. CNM and tCNM achieve a significant improvement in the
representation accuracy by counting the state propagation of the network
model, especially for the periodic flow. The improvement for the other flow
regimes is relatively unapparent, mainly due to the increasing intra-cluster
variance for the more complex flow features. Compared with CNM, the
representation errors of tCNM on periodic flow, quasi-periodic flow, and
chaotic flow are reduced by $60\%$, $31\%$, and $18\%$, respectively, with an
average of $27\%$.

##  VI Conclusions

In this work, we propose a novel data-driven reduced-order modeling
methodology for time-resolved data. The starting point is the Cluster-based
Network Model (CNM) (Fernex, Noack, and Semaan, 2021; Li  _et al._ , 2021).
The $k$-means++ algorithm optimizes the representation error of the snapshots
by the closest centroids. Kaiser  _et al._ (2014) proposed a Cluster-based
Markov Model to analyze the flow dynamics between these representative
centroids. Yet, the CNM includes time-continuous ‘flights’ between neighboring
centroids and builds a directed graph emphasizing the non-trivial transitions
between nodes in the network. The key innovation of the proposed trajectory-
optimized Cluster-based Network Model (tCNM) is to shift the centroids and
allow constraints between centroids in order to decrease the kinematic
representation error of interpolated states and thus to increase the accuracy
of the reconstructed dynamics. In particular, the snapshots are coarse-grained
into a few clusters following the $k$-means++ algorithm. Then, the centroids
$\bm{c}_{k}$ are shifted closer to the snapshot trajectory, meanwhile, the
trajectory support points $\bm{r}_{kj}$ are built in order to refine the
transitions between two cluster centroids.

The resulting tCNM has several desirable features:

  1. 1.

tCNM retains all the advantages of CNM. The framework does not require any
first-principle knowledge. It can be extended to incorporate multiple
operating conditions or control-oriented predictions beyond the training data
.

  2. 2.

The representation error is reduced, without losing simplicity and robustness.
The propagation trajectory is more accurate, especially near the critical
moment when the vortex is about to generate or separate.

  3. 3.

The off-line computational load is only slightly larger than CNM. The extra
computations are mainly due to the building of the new nodes, which requires
only a tiny fraction of total computational operations.

  4. 4.

The methodology is data-driven and fully automatable, including the centroid
shift and the inter-cluster constraint.

However, tCNM, or any other data-driven model order reduction, requires a
careful calibration of the observation domain. If the domain is too small,
future states are difficult to predict and challenge the transition model. If,
on the other hand, the domain comprises too many uncorrelated flow structures,
the number of necessary centroids becomes excessive. Also, the required
snapshot data will become impractically large, as each centroid represents a
coarse-grained recurrent state which should be found in multiple snapshots. We
note that similar observations also hold for POD.

The three-dimensional wake behind a sphere of three different dynamics is
considered for the evaluation of tCNM: (1) The periodic flow features a
periodic shedding of hairpin vortices at $Re=300$. (2) The quasi-periodic flow
exhibits an irregular shedding with a long-time periodicity at $Re=350$. (3)
The chaotic flow shows a fully three-dimensional irregular shedding with the
separation point slowly rotating around the sphere at $Re=450$. For all three
cases, the snapshots are coarse-grained into $K=10$ clusters yielding
satisfactory dynamics predictions.

On a kinematic level, clustering discretizes the flow states. Each centroid
$\bm{c_{k}^{0}}$ represents realizable states modulo a coarse graining. In
fact, the $K$ centroids constitute a complete set of probable and recurrent
states with distinct vortex formations. The proposed tCNM shifts of the
centroids improve the representation error to POD level. In contrast, linear
combinations of POD modes are needed to approximate flows, and the individual
modes typically consist of a mixture of frequencies. Both features mitigate
their interpretability. On a dynamic level, the time-relevant dynamics between
the snapshots are represented by transitions between different clusters, or,
in other words between different centroids. These transitions illustrated by
the transition matrixes serve as a deterministic–stochastic gray-box model
resolving the coherent-structure evolution, e.g., the deterministic limit
cycle for a periodic flow regime; the stochastic model for a chaotic flow
regime, which chooses the next destination according to the probability of
transitions. Hence, CNM and tCNM can automatically capture the dynamics of the
unsteady flow field, while POD models can only capture the energy-dominate
modes and require additional nonlinear modeling for the dynamics. Even though
CNM well resolves the temporal evolution of the flow dynamics, a noteworthy
deficiency can be seen in dealing with the representation error. However,
compared with CNM, tCNM achieved a reduction of $27\%$ for the representation
error on average of the three flow regimes.

tCNM opens a novel automatable avenue for data-driven nonlinear reduced-order
modeling. tCNM enables cluster-based methods with optimized trajectory to
enhance linear interpolation, thereby reducing representation error and
improving the accuracy of the dynamical model. The methodology can be applied
to numerous other applications and extended to parametric and control-oriented
reduced-order models.

As outlook, tCNM foreshadows the advantages of a new optimization paradigm of
clustering: the centroids may be chosen to minimize representation errors with
respect to local interpolations. In addition, the trajectory models may be
further refined with respect to historical data, like for the higher-order
Markov chain model (Fernex, Noack, and Semaan, 2021). The authors are actively
pursuing these avenues.

###### Acknowledgements.

Chang Hou appreciates the support from the School of Mechanical Engineering
and Automation during his PhD thesis at Harbin Institute of Technology
(Shenzhen). This work is supported by the National Natural Science Foundation
of China (NSFC) under grant 12172109 and 12172111, and by the Natural Science
and Engineering grant 2022A1515011492 of Guangdong province, China. We
appreciate valuable discussions with Guy Y. Cornejo-Maceda, Hao Li, Richard
Semaan and the HIT-Hanghua team. Last but not least, we thank the referees for
their many constructive suggestions. The authors declare that there is no
conflict of interest regarding the publication of this paper.

## Appendix: Cluster centroids versus POD modes

Figure 31: Comparison between the tCNM centroids and the first POD modes for
the periodic flow at $Re=300$. (a) Graph of transitions between clusters. The
centroids are marked with large gray dots in the proximity map, and their
vortex structures are visualized like in figure figure 6. The snapshots
belonging to different clusters are marked with small dots in different
colors. (b) The leading 6 POD modes. The vortex structures are visualized by
the $Q$-criterion ($Q=0.001$). Figure 32: Same as figure 31, but for chaotic
flow at $Re=450$. (a) Graph of transitions between clusters. (b) The leading 6
POD modes.

The tCNM centroids and the leading 6 POD modes are compared in figure 31 for
periodic flow and in figure 32 for the chaotic flow of the sphere wake. After
clustering, the time-resolved snapshots of the sphere wake are grouped into
different clusters depending on their flow characteristics. Each cluster is
represented by its centroid. For tCNM, these snapshots are shifted towards the
trajectory and can therefore be seen as more realistic flow states. In
contrast, POD modes come with a mathematically optimal representation error
but are typically not physically interpretable flow states: First, POD needs
the superposition of individual modes to approximate the snapshot data.
Second, POD modes typically comprise multiple frequencies, making the physical
interpretation difficult. Third, POD modes represent a first and second-moment
statistical analysis ignoring the temporal dynamics. The last two features
might be improved with other empirical Galerkin expansions like DMD (Rowley
_et al._ , 2009; Schmid, 2010b), recursive DMD (Noack  _et al._ , 2016),
spectral (Sieber, Paschereit, and Oberleithner, 2016), etc. by sacrificing the
optimal representation error.

## Data Availability Statement

The data that support the findings of this study are available from the first
author upon reasonable request.

## References

  * Pullin and Saffman (1998) D. Pullin and P. Saffman, “Vortex dynamics in turbulence,” Annual Review of Fluid Mechanics 30, 31–51 (1998).
  * Cottet and Koumoutsakos (2000) G. H. Cottet and P. Koumoutsakos,  _Vortex Methods — Theory and Practice_ (Cambridge University Press, Cambridge, 2000).
  * Föppl (1913) L. Föppl, “Wirbelbewegung hinter einen Kreiszylinder,” Sitzb. d. k. bayr. Akad. d. Wiss. 1, 1–18 (1913).
  * von Kármán (1911) T. von Kármán, “über den Mechanismus des Widerstandes, den ein bewegter Körper in einer Flüssigkeit erfährt,” Nachrichten der Kaiserlichen Gesellschaft der Wissenschaften zu Göttingen , 509–517 (1911).
  * Galerkin (1915) B. G. Galerkin, “Rods and plates: series occuring in various questions regarding the elastic equilibrium of rods and plates (translated),” Vestn. Inzhen. 19, 897–908 (1915).
  * Fletcher (1984) C. A. J. Fletcher,  _Computational Galerkin Methods_ , 1st ed. (Springer, New York, 1984).
  * Holmes  _et al._ (2012) P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley,  _Turbulence, Coherent Structures, Dynamical Systems and Symmetry_ , 2nd ed. (Cambridge University Press, Cambridge, 2012).
  * Stuart (1971) J. T. Stuart, “Nonlinear stability theory,” Annual Review of Fluid Mechanics 3, 347–370 (1971).
  * Busse (1991) F. H. Busse, “Numerical analysis of secondary and tertiary states of fluid flow and their stability properties,” Applied Scientific Research 48, 341–351 (1991).
  * Noack and Eckelmann (1994) B. R. Noack and H. Eckelmann, “A global stability analysis of the steady and periodic cylinder wake,” Journal of Fluid Mechanics 270, 297–330 (1994).
  * Sirovich (1987) L. Sirovich, “Turbulence and the dynamics of coherent structures, Part I: Coherent structures,” Quarterly of Applied Mathematics XLV, 561–571 (1987).
  * Taira  _et al._ (2017) K. Taira, S. L. Brunton, S. T. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: An overview,” AIAA Journal 55, 4013–4041 (2017).
  * Aubry  _et al._ (1988) N. Aubry, P. Holmes, J. L. Lumley, and E. Stone, “The dynamics of coherent structures in the wall region of a turbulent boundary layer,” Journal of Fluid Mechanics 192, 115–173 (1988).
  * Podvin and Lumley (1998) B. Podvin and J. Lumley, “A low-dimensional approach for the minimal flow unit,” Journal of Fluid Mechanics 362, 121–155 (1998).
  * Podvin (2009) B. Podvin, “A proper-orthogonal-decomposition–based model for the wall layer of a turbulent channel flow,” Physics of Fluids 21, 015111 (2009).
  * Bergmann, Bruneau, and Iollo (2009) M. Bergmann, C.-H. Bruneau, and A. Iollo, “Enablers for robust pod models,” Journal of Computational Physics 228, 516–538 (2009).
  * Noack, Papas, and Monkewitz (2005) B. R. Noack, P. Papas, and P. A. Monkewitz, “The need for a pressure-term representation in empirical galerkin models of incompressible shear flows,” Journal of Fluid Mechanics 523, 339–365 (2005).
  * Ahmed  _et al._ (2021) S. E. Ahmed, S. Pawar, O. San, A. Rasheed, T. Iliescu, and B. R. Noack, “On closures for reduced order models—a spectrum of first-principle to machine-learned avenues,” Physics of Fluids 33, 091301:1–33 (2021).
  * Brunton, Proctor, and Kutz (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Proceedings of the National Academy of Sciences 113, 3932–3937 (2016).
  * Rowley (2005) C. W. Rowley, “Model reduction for fluids, using balanced proper orthogonal decomposition,” International Journal of Bifurcation and Chaos 15, 997–1013 (2005).
  * Mendez, Balabane, and Buchlin (2019) M. Mendez, M. Balabane, and J.-M. Buchlin, “Multi-scale proper orthogonal decomposition of complex fluid flows,” Journal of Fluid Mechanics 870, 988–1036 (2019).
  * Rowley  _et al._ (2009) C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. Henningson, “Spectral analysis of nonlinear flows,” Journal of Fluid Mechanics 645, 115–127 (2009).
  * Schmid (2010a) P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” Journal of Fluid Mechanics 656, 5–28 (2010a).
  * Kutz  _et al._ (2016) J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor,  _Dynamic mode decomposition: data-driven modeling of complex systems_ (SIAM, 2016).
  * Noack  _et al._ (2016) B. R. Noack, W. Stankiewicz, M. Morzyński, and P. J. Schmid, “Recursive dynamic mode decomposition of transient and post-transient wake flows,” Journal of Fluid Mechanics 809, 843–872 (2016).
  * Noack (2016) B. R. Noack, “From snapshots to modal expansions–bridging low residuals and pure frequencies,” Journal of Fluid Mechanics 802, 1–4 (2016).
  * San, Maulik, and Ahmed (2019) O. San, R. Maulik, and M. Ahmed, “An artificial neural network framework for reduced order modeling of transient flows,” Communications in Nonlinear Science and Numerical Simulation 77, 271–287 (2019).
  * Zhu  _et al._ (2019) L. Zhu, W. Zhang, J. Kou, and Y. Liu, “Machine learning methods for turbulence modeling in subsonic flows around airfoils,” Physics of Fluids 31, 015105 (2019).
  * Kutz (2017) J. N. Kutz, “Deep learning in fluid dynamics,” Journal of Fluid Mechanics 814, 1–4 (2017).
  * Burkardt, Gunzburger, and Lee (2006) J. Burkardt, M. Gunzburger, and H.-C. Lee, “POD and CVT-based reduced-order modeling of navier–stokes flows,” Computer Methods in Applied Mechanics and Engineering 196, 337–355 (2006).
  * Kaiser  _et al._ (2014) E. Kaiser, B. R. Noack, L. Cordier, A. Spohn, M. Segond, M. Abel, G. Daviller, J. Östh, S. Krajnović, and R. K. Niven, “Cluster-based reduced-order modelling of a mixing layer,” Journal of Fluid Mechanics 754, 365–414 (2014).
  * Cao  _et al._ (2014) Y. Cao, E. Kaiser, J. Borée, B. R. Noack, L. Thomas, and S. Guilain, “Cluster-based analysis of cycle-to-cycle variations: application to internal combustion engines,” Experiments in Fluids 55, 1–8 (2014).
  * Li and Tan (2020) H. Li and J. Tan, “Cluster-based markov model to understand the transition dynamics of a supersonic mixing layer,” Physics of Fluids 32, 056104 (2020).
  * Li  _et al._ (2021) H. Li, D. Fernex, R. Semaan, J. Tan, M. Morzyński, and B. R. Noack, “Cluster-based network model,” Journal of Fluid Mechanics 906 (2021).
  * Fernex, Noack, and Semaan (2021) D. Fernex, B. R. Noack, and R. Semaan, “Cluster-based network modeling—from snapshots to complex dynamical systems,” Science Advances 7, eabf5006 (2021).
  * Barabási (2013) A.-L. Barabási, “Network science,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120375 (2013).
  * Arenas  _et al._ (2008) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, “Synchronization in complex networks,” Physics Reports 469, 93–153 (2008).
  * Deng  _et al._ (2022) N. Deng, B. R. Noack, M. Morzyński, and L. R. Pastur, “Cluster-based hierarchical network model of the fluidic pinball–cartographing transient and post-transient, multi-frequency, multi-attractor behaviour,” Journal of Fluid Mechanics 934 (2022).
  * Nair  _et al._ (2019) A. G. Nair, C.-A. Yeh, E. Kaiser, B. R. Noack, S. L. Brunton, and K. Taira, “Cluster-based feedback control of turbulent post-stall separated flows,” Journal of Fluid Mechanics 875, 345–375 (2019).
  * Möller (1937) W. Möller,  _Experimentelle Untersuchungen Zur Hydrodynamik Der Kugel_ (Hirzel, 1937).
  * Achenbach (1974) E. Achenbach, “Vortex shedding from spheres,” Journal of Fluid Mechanics 62, 209–221 (1974).
  * Sakamoto and Haniu (1990) H. Sakamoto and H. Haniu, “A study on vortex shedding from spheres in a uniform flow,” ASME Journal of Fluids Engineering 112, 386–392 (1990).
  * Leweke  _et al._ (1999) T. Leweke, M. Provansal, D. Ormieres, and R. Lebescond, “Vortex dynamics in the wake of a sphere,” Physics of Fluids 11, S12–S12 (1999).
  * Chrust, Goujon-Durand, and Wesfreid (2013) M. Chrust, S. Goujon-Durand, and J. Wesfreid, “Loss of a fixed plane of symmetry in the wake of a sphere,” Journal of Fluids and Structures 41, 51–56 (2013).
  * Johnson and Patel (1999) T. Johnson and V. Patel, “Flow past a sphere up to a reynolds number of 300,” Journal of Fluid Mechanics 378, 19–70 (1999).
  * Kim, Kim, and Choi (2001) J. Kim, D. Kim, and H. Choi, “An immersed-boundary finite-volume method for simulations of flow in complex geometries,” Journal of Computational Physics 171, 132–150 (2001).
  * Ploumhans  _et al._ (2002) P. Ploumhans, G. Winckelmans, J. K. Salmon, A. Leonard, and M. Warren, “Vortex methods for direct numerical simulation of three-dimensional bluff body flows: application to the sphere at Re= 300, 500, and 1000,” Journal of Computational Physics 178, 427–463 (2002).
  * Yun, Kim, and Choi (2006) G. Yun, D. Kim, and H. Choi, “Vortical structures behind a sphere at subcritical reynolds numbers,” Physics of Fluids 18, 015102 (2006).
  * Giacobello, Ooi, and Balachandar (2009) M. Giacobello, A. Ooi, and S. Balachandar, “Wake structure of a transversely rotating sphere at moderate reynolds numbers,” Journal of Fluid Mechanics 621, 103–130 (2009).
  * Rajamuni, Thompson, and Hourigan (2018) M. M. Rajamuni, M. C. Thompson, and K. Hourigan, “Transverse flow-induced vibrations of a sphere,” Journal of Fluid Mechanics 837, 931–966 (2018).
  * Zhao (2021) M. Zhao, “Flow past a circular cylinder and a downstream sphere for re< 300,” Journal of Fluid Mechanics 913 (2021).
  * Chizfahm and Jaiman (2021) A. Chizfahm and R. Jaiman, “Data-driven stability analysis and near-wake jet control for the vortex-induced vibration of a sphere,” Physics of Fluids 33, 044104 (2021).
  * Berkooz, Holmes, and Lumley (1993) G. Berkooz, P. Holmes, and J. L. Lumley, “The proper orthogonal decomposition in the analysis of turbulent flows,” Annual Review of Fluid Mechanics 25, 539–575 (1993).
  * Tan, Li, and Noack (2020) J. Tan, H. Li, and B. R. Noack, “On the cavity-actuated supersonic mixing layer downstream a thick splitter plate,” Physics of Fluids 32, 096102 (2020).
  * Protas, Noack, and Östh (2015) B. Protas, B. R. Noack, and J. Östh, “Optimal nonlinear eddy viscosity in galerkin models of turbulent flows,” Journal of Fluid Mechanics 766, 337–367 (2015).
  * Fabre, Auguste, and Magnaudet (2008) D. Fabre, F. Auguste, and J. Magnaudet, “Bifurcations and symmetry breaking in the wake of axisymmetric bodies,” Physics of Fluids 20, 051702 (2008).
  * Pan, Zhang, and Ni (2018) J.-H. Pan, N.-M. Zhang, and M.-J. Ni, “The wake structure and transition process of a flow past a sphere affected by a streamwise magnetic field,” Journal of Fluid Mechanics 842, 248–272 (2018).
  * Lorite-Díez and Jiménez-González (2020) M. Lorite-Díez and J. Jiménez-González, “Description of the transitional wake behind a strongly streamwise rotating sphere,” Journal of Fluid Mechanics 896 (2020).
  * Eshbal  _et al._ (2019) L. Eshbal, V. Rinsky, T. David, D. Greenblatt, and R. van Hout, “Measurement of vortex shedding in the wake of a sphere at Re = 465,” Journal of Fluid Mechanics 870, 290–315 (2019).
  * Schmid (2010b) P. J. Schmid, “Dynamic mode decomposition for numerical and experimental data,” J. Fluid. Mech 656, 5–28 (2010b).
  * Sieber, Paschereit, and Oberleithner (2016) M. Sieber, C. O. Paschereit, and K. Oberleithner, “Spectral proper orthogonal decomposition,” Journal of Fluid Mechanics 792, 798–828 (2016).