2.1. Clustering as coarse-graining

We consider velocity fields in a steady domain $\varOmega$ which may be obtained from experiments or from numerical simulations. Starting point is an ensemble of $M$ statistically representative, time-resolved snapshots as employed for cluster-based models (Kaiser et al. Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014), Dynamic mode decomposition (DMD) (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Schmid Reference Schmid2010) or POD (see e.g. Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012). The velocity field is equidistantly sampled with time step $\Delta t$, i.e. the $m$th instant reads $t^{m}=m \Delta t$. The corresponding snapshot velocity field is denoted by $\boldsymbol {u}^{m}(\boldsymbol {x} ) :=\boldsymbol {u}(\boldsymbol {x}, t^{m})$, $m=1,\ldots ,M$.

Cluster analysis lumps similar objects into clusters as illustrated in figure 2. This lumping of data is performed in an unsupervised manner, i.e. no advance labelling or grouping of the data has been performed. In cluster-based models, the $M$ snapshots $\boldsymbol {u}^{m} (\boldsymbol {x} )$ are coarse-grained into $K$ clusters represented by the centroids $\boldsymbol {c}_k (\boldsymbol {x})$, $k=1,\ldots ,K$, using the unsupervised k-means++ algorithm (Steinhaus Reference Steinhaus1956; MacQueen Reference MacQueen1967; Lloyd Reference Lloyd1982). The centroids characterize the typical flow patterns of each cluster, also called modes in the ROM community. The corresponding cluster-affiliation function maps a velocity field $\boldsymbol {u}$ to the index of the closest centroid,

(2.1)\begin{equation} k (\boldsymbol{u}) = \arg\min_i \Vert \boldsymbol{u} - \boldsymbol{c}_i \Vert_{\varOmega}, \end{equation}

where $\Vert \cdot \Vert _{\varOmega }$ denotes the standard Hilbert space norm in the domain $\varOmega$ (see appendix A). This function defines cluster regions as Voronoi cells around the centroids

(2.2)\begin{equation} \mathcal{C}_i = \left \{ \boldsymbol{u} \in \mathcal{L}^2 ( \varOmega ) \colon k (\boldsymbol{u}) = i \right \}. \end{equation}

This function can also be employed to map the snapshot index $m$ to the representative cluster index $k(m) := k ( \boldsymbol {u}^m )$. Alternatively, the characteristic function

(2.3)\begin{equation} \chi_i^m := \left\{\begin{array}{ll} 1, & \text{if} \ i=k(m), \\ 0, & \text{otherwise}, \end{array}\right. \end{equation}

describes whether the $m$th snapshot is affiliated with the $l$th centroid. The latter two quantities are equivalent.

Figure 2. Clustering exemplified in a two-dimensional state space. The snapshots $\{\boldsymbol {u}^m\}_{m=1}^M$ are coarse-grained into three clusters. Each centroid $\boldsymbol {c}_{k}, k=1, 2 ,3$, is the averaged field of all snapshots belonging to this cluster. Every centroid $\boldsymbol {c}_k$ is associated with a Voronoi cell $\mathcal {C}_k$, i.e. a region in which the points are closer to $\boldsymbol {c}_k$ than any other centroid. The cluster affiliation for a snapshot is the cluster index of the closest centroid which has been indicated by an arrow. By definition, the cluster affiliation is the index of the corresponding Voronoi cell.

The performance of a set of centroids $\{ \boldsymbol {c}_k \}_{k=1}^K$ with respect to a given set of snapshots $\{ \boldsymbol {u}^m \}_{m=1}^M$ is measured by the average variance of the snapshots with respect to their closest centroid. The corresponding inner-cluster variance reads

(2.4)\begin{equation} J \left( \boldsymbol{c}_1,\ldots,\boldsymbol{c}_K \right ) = \frac{1}{M} \sum_{m=1}^{M} \left \Vert \boldsymbol{u}^{m}-\boldsymbol{c}_{k(m)} \right \Vert^2_{\varOmega} . \end{equation}

The optimal centroids $\{ \boldsymbol {c}_k^{\star } \}_{k=1}^K$ minimize this inner-cluster variance:

(2.5)\begin{equation} \left ( \boldsymbol{c}_1^{\star}, \ldots, \boldsymbol{c}_K^{\star} \right ) =\underset{\boldsymbol{c}_1,\ldots,\boldsymbol{c}_K}{\arg\min} J \left( \boldsymbol{c}_1,\ldots,\boldsymbol{c}_K \right ). \end{equation}

The argument is indeterminate with respect to a re-ordering. For CMM, we chose the first cluster as the one with the highest population, i.e. the largest number of associated snapshots. The $(k+1)$th cluster, $k>1$, has the largest transition probability from the $k$th one.

The optimization problem (2.5) is solved by the k-means++ algorithm. The $K$ centroids are initialized randomly and then iterated until convergence is reached or when the variance $J$ is small enough; k-means++ repeats the clustering process 30 times and take the best set of centroids.

The number of snapshots $n_k$ in cluster $k$ is given by

(2.6)\begin{equation} n_k=\sum_{m=1}^{M} \chi_k^{m}. \end{equation}

The centroids are the mean velocity field of all snapshots in the corresponding cluster. In other words,

(2.7)\begin{equation} \boldsymbol{c}_k = \frac{1}{n_k} \sum_{\boldsymbol{u}^{m} \in \mathcal{C}_k} \boldsymbol{u}^{m}=\frac{1}{n_k} \sum_{m=1}^{M} \chi_k^{m} \boldsymbol{u}^{m}. \end{equation}

In the following centroid visualizations, we accentuate the vortical structures by displaying the fluctuations $\boldsymbol {c}_k - \bar {\boldsymbol {u}}$ around the snapshot mean $\bar {\boldsymbol {u}}$ and not the full velocity field $\boldsymbol {c}_k$.

2.2. Cluster-based Markov model

We briefly recapitulate CMM by Kaiser et al. (Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014) as our benchmark cluster-based reduced-order model. In CMM, the state variable is the cluster population $\boldsymbol {p}=[ p_1,\ldots ,p_K ]^{\mathrm {T}}$, where $p_i$ represents the probability of being in cluster $i$ and the superscript $\mathrm {T}$ denotes the transpose. The transition between clusters in a given time step $\Delta t^c$ is described by the transition matrix ${\boldsymbol{\mathsf{P}}} = (P_{ij}) \in \mathcal {R}^{K \times K}$. The superscript ‘$c$’ refers to cluster-based model. Here, $P_{ij}$ is the transition probability of moving from cluster $j$ to cluster $i$. Let $\boldsymbol {p}^l$ be the probability vector at time $t^l = l \Delta t^c$, then the change in one time step is described by

(2.8)\begin{equation} \boldsymbol{p}^{l+1} = {\boldsymbol{\mathsf{P}}} \boldsymbol{p}^l. \end{equation}

With increasing iterations, the iteration (2.8) converges to the asymptotic probability $\boldsymbol {p}^{\infty } := \lim _{l\to \infty } \boldsymbol {p}^l$. In a typical case, (2.8) has a single fixed point $\boldsymbol {p}^{\infty }$. For completeness, a continuous form of Markov models with a new transition matrix ${{\boldsymbol{\mathsf{P}}}}^c$ is mentioned

(2.9)\begin{equation} \frac{\textrm{d}\boldsymbol{p}}{\textrm{d}t}={\boldsymbol{\mathsf{P}}}^{c} \boldsymbol{p}. \end{equation}

From the time-continuous form (2.9), the time-discrete one (2.8) can be derived. The opposite is not generally true. In the following, no continuous Markov models are employed.

A CMM of the time-resolved snapshots starts with cluster affiliation (2.1) which can also be considered function of time $k(t)$. We refer to the original paper for the determination of $P_{ij}$ from $k(t)$. The time step $\Delta t^c$ is a critical design parameter for CMM. A good choice is a value where the transition from one cluster to the next is likely. If the time step is too small, the Markov model idles many times in each cluster for a stochastic number of times before transitioning to the next cluster. The model-based transition time may thus significantly deviate from the deterministic data-driven trajectories through the clusters. If the time step is too large, one may miss intermediate clusters. We choose $\Delta t^c= T/10$, where $T$ is the characteristic period of the flow. On a circular limit cycle with uniform rotation, this value is optimal for $K=10$ clusters, enforcing the transition from one cluster to the next in each time step.

In figure 3(a,c,e), the effect of the suboptimal time step is illustrated for the CMM of a uniform rotation $u_1 = \cos ( 2{\rm \pi} t)$, $u_2 = \sin (2 {\rm \pi}t)$. Here, four clusters and a time step $\Delta t^c = 1/16$ are chosen. The probability of staying in the cluster during one time step is $P_{11}=P_{22}=P_{33}=P_{44}= 3/4$ and the transition probability to the next counter-clockwise neighbour is $P_{14}=P_{21}=P_{32}=P_{43}= 1/4$. Thus, the probability of staying in one cluster for $l$ steps exponentially decays, $P_{11}^l$. In contrast, the uniform rotation commands that the state is exactly three time steps in one cluster before it leaves in the fourth step. This example motivates the proposed cluster-based reduced-order model, foreshadowed in figure 3(b,d,f) and explained in the following section.

Figure 3. Introduction of cluster-based models for a limit-cycle example. The CMM and CNM are displayed in (a,c,e) and (b,d,f), respectively. The uniform rotation $u_{1}=\cos (2{\rm \pi} t)$, $u_{2}=\sin (2{\rm \pi} t)$ in a two-dimensional plane is discretized by four centroids $\boldsymbol {c}_{k}$, $k=1,\ldots ,4$. (a) Phase portrait of CMM with time step $\Delta t=1/16$. The centroids are near the limit cycle (red dashed line). The state vector residing in centroid $\boldsymbol {c}_{i}$ has the probability $P_{ii}$ of staying in its state and $P_{ji}$ of transitioning to centroid $\boldsymbol {c}_{j}$ in the considered time step. (b) Phase portrait of the CNM. The state in centroid $\boldsymbol {c}_{i}$ moves uniformly to its counter-clockwise neighbour taking a quarter period $T_{14}=T_{21}=T_{32}=T_{43}=1/4$. Here, $Q_{14}=Q_{21}=Q_{32}=Q_{43}=1$ and $Q_{ij}=0$ otherwise. The estimated probability evolution starting in cluster $i=1$ at $t=0$ is illustrated for CMM (c) and CNM (d). Panels (e) and (f) present the model-based evolution of the first coordinate $u_{1}$ for CMM and CNM, respectively. In (c–f), the solid black lines correspond to the uniform rotation and the dashed red line to the model.

2.3. Cluster-based network model

For CMM, the time step $\Delta t^c$ is, as mentioned, an important design parameter. This design parameter can be avoided by the new proposed CNM. The key idea is to abandon the ‘stroboscopic’ view of CMM and focus on non-trivial transitions from cluster $j$ to cluster $i$. These transitions are characterized by two parameters: the probability $Q_{ij}$ and a time scale $T_{ij}$. Evidently, no time step is needed for the description and the assumption of a constant transition time is found to be much more aligned with shear flow modes than assuming an exponential decay of residence time. Moreover, it could be relaxed by assuming a probability distribution of transition times.

In the following, the transition probability and transition time are inferred from the cluster affiliation function $k(t)$. The continuous form is convenient for discussion. The time-discrete affiliation function $k(m)$ can be made continuous by taking the cluster of the snapshot which is closest in time:

(2.10)\begin{equation} k(t) = k \left ( \arg\min_m \left\vert t - t_m \right\vert \right). \end{equation}

The $n$th transition time $t_n$ of the cluster affiliation is recursively defined as the first discontinuity of $k(t)$ for $t>t_{n-1}$. Here, $t_0 = 0$. The transition time $t_n$ satisfies

(2.11)\begin{equation} k \left ( t_{n}-\varepsilon \right ) \neq k \left ( t_{n}+\varepsilon \right ) \end{equation}

for any sufficiently small positive $\varepsilon$. For $t \in ( t_n, t_{n+1} )$, the data-based trajectory is assumed to stay in cluster $k$ at the averaged time $(t_{n+1}+t_n)/2$ (see figure 4). The residence time in this cluster is defined by

(2.12)\begin{equation} \tau_{n} :=t_{n+1}-t_{n}. \end{equation}

Let $j$ and $i$ be the indices of the clusters after $t_n$ and $t_{n+1}$ respectively. Then the transition time from $j$ to $i$ is defined as half of the residence time of both clusters:

(2.13)\begin{equation} \tau_{ij} := \frac{\tau_{n}+\tau_{n+1}}{2} = \frac{t_{n+2}-t_n}{2}. \end{equation}

This definition may appear arbitrary but is the least-biased guess consistent with the available data. The sum of all residence times from a given data set add up to the total investigated time period $T_0$.

Figure 4. Sketch of times and periods employed in the cluster-based network model; $\times$ marks the centre of the cluster residence time, while $\bullet$ denotes the transition between clusters.

The direct transition probability $Q_{ij}$ and transition time $T_{ij}$ can be inferred from the data. Then,

(2.14)\begin{equation} Q_{ij}=\frac{n_{ij}} {n_{j}},\quad i,j=1,\ldots,K, \end{equation}

where $n_{ij}$ is the number of transitions from $\boldsymbol {c}_{j}$ to $\boldsymbol {c}_{i}$ and $n_j$ the number of transitions departing from $\boldsymbol {c}_{j}$ regardless of the destination point

(2.15)\begin{equation} n_j=\sum_{i=1}^K n_{i j},\quad i,j=1,\ldots,K. \end{equation}

We emphasize that $n_{ii}=0$ for $i=1,\ldots ,K$ by very definition of a direct transition. The direct transition matrix (DTM) ${\boldsymbol{\mathsf{Q}}} = (Q_{ij} ) \in \mathcal {R}^{K \times K}$ lumps these probabilities into a single entity.

Similarly, the direct transition time $T_{ij}$ from cluster $j$ to cluster $i$ is taken to be the average of all values. This average is symbolically denoted by

(2.16)\begin{equation} T_{i j}= \langle \tau_{ij} \rangle. \end{equation}

These values are lumped into the matrix ${\boldsymbol{\mathsf{T}}} = (T_{ij} ) \in \mathcal {R}^{K \times K}$.

It should be noted that a given trajectory may repeatedly pass through the same clusters (Voronoi cells) with different residence and transition times. With enough data, this variability may be incorporated into the model. Our goal is to compare the Markov model with the most simple network model where constant (averaged) transition times are assumed.

The CNM predicts the asymptotic cluster probability $p^{\infty }_i$ in cluster $i$. Let $[0,T_0]$ be a sufficiently long time horizon simulated by the model. Then, the probability of staying in cluster $i$ is the cumulative residence time normalized by the simulation time

(2.17)\begin{equation} p^{\infty}_i= \frac{\sum \tau_i}{T_0}. \end{equation}

We return to the introductory example depicted in figure 3. The CNM is seen to accurately describe the uniform rotation (figure 3f) and correctly yields the asymptotic cluster probability $p_i^{\infty } =1/4$, $i=1,\ldots ,4$. In contrast, the prediction horizon of the CMM is limited to approximately one period. After this time, the initial condition is forgotten and the asymptotic distribution is reached – rendering CMM unsuitable for dynamic prediction. However, CMM predicts the asymptotic state faster than the CNM. For this particular example, the CMM could be made equivalent to the CNM by choosing $\Delta t^c=1/4$. However, the Markov model will inevitably diffuse the state with a range of cluster transition times, e.g. for non-uniform rotation or more complex dynamics.

2.4. Velocity fields associated with the cluster-based reduced-order models

The CMM describes the cluster population

(2.18)\begin{equation} \boldsymbol{p}=\left[ p_{1},\ldots,p_{K} \right]^{\textrm{T}} \end{equation}

at discrete times $t=l \Delta t^c$. In the following, this population is considered to be continuous in time, e.g. by using linear or higher-order interpolation. The corresponding velocity field $\boldsymbol u(\boldsymbol x ,t)$ at time $t$ is defined as the expectation value,

(2.19)\begin{equation} \boldsymbol u(\boldsymbol x ,t) =\sum_{i=1}^{K} p_{i}(t) \boldsymbol c_{i}(\boldsymbol x), \end{equation}

where $\boldsymbol c_{i}$ is the $i$th centroid.

The CNM is based on centroid visits at discrete times. The clusters $k_0, k_1, k_2, \ldots$ are visited at times

(2.20)\begin{equation} t_0=0, \quad t_1=T_{k_1 k_0}, \quad t_2= t_1 + T_{k_2 k_1}, \ldots \end{equation}

consistent with the direct transition matrix $(Q_{ij})$ and the transition times $T_{ij}$. A uniform motion is assumed between these visits. In other words, for $t \in [t_n,t_{n+1}]$ the velocity field reads

(2.21)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \alpha_{n} (t) \boldsymbol{c}_{k_n} (\boldsymbol{x}) + \left [ 1-\alpha_n(t) \right ] \boldsymbol{c}_{k_{n+1}} (\boldsymbol{x}), \quad \alpha_n = \frac{t_{n+1}-t}{t_{n+1}-t_n} . \end{equation}

We note that a smoother motion may be achieved with splines.

The actual flow computations are based on a lossless POD, as elaborated in appendix A. The interpolations are performed with the mode amplitudes $\boldsymbol {a}= [ a_1,\ldots ,a_N ]^{\textrm {T}}$ before being transcribed into velocity fields via the POD expansion.

Figure 5 compares the classical CMM with stroboscopic temporal prediction of discrete states and the proposed CNM with time-continuous uniform motion on a network of routes between two centroids. In the top row, the possible states are illustrated. In case of the CMM, the states (2.19), denoted by red dots, quickly converge to the mean flow, like Reynolds-averaged Navier–Stokes simulations. The CNM-predicted state (2.21) moves on the directed network, marked by red arrows, and is reminiscent of large-eddy simulations. As displayed in the middle row, the CMM is discrete in time while the CNM dynamics is time continuous with cluster visits after pre-specified transition times $T_{jk}$. The bottom row shows another difference: CMM describes averages over all centroids while the CNM only allows for linear interpolations between two neighbouring centroids. This interpolation is consistent with the purpose of accurately resolving evolving coherent structures. Averaging over many centroids acts like a low-pass filter mitigating the fluctuation level and thus underresolving the coherent structures.

Figure 5. Comparison of the CMM and CNM. For reasons of simplicity, an example with four centroids $i$, $j$, $k$ and $l$ is shown. The CNM exemplifies the evolution of the weight of centroid $k$. For the CMM, the transition from cluster $k$ to $j$ is shown. The predicted state (2.21) is determined by the weights $w_k(t) = \alpha (t)$ and $w_j = 1-\alpha (t)$, $\alpha =(t_{n+1}-t)/T_{jk}$ in the time interval $[ t_n,t_{n+1} ]$. For details, see text.

2.5. Validation of the cluster-based reduced-order models

Following Protas, Noack & Östh (Reference Protas, Noack and Östh2015), the cluster-based model is validated based on the computed and predicted autocorrelation function of the velocity field. The unbiased autocorrelation function reads

(2.22)\begin{equation} R(\tau) :=\frac{1}{T-\tau} \int_{\tau}^{\textrm{T}} \int_{\varOmega} \boldsymbol{u}(\boldsymbol{x}, t-\tau) \boldsymbol{\cdot} \boldsymbol{u}(\boldsymbol{x}, t)\,\textrm{d}\boldsymbol{x}\,\textrm{d}t, \quad \tau \in[0, T). \end{equation}

This function reveals the turbulent fluctuation level $R(0)$ and the frequency spectrum. Moreover, the problem of comparing two trajectories with finite dynamic prediction horizons due to the increasing phase mismatch is avoided (Pastoor et al. Reference Pastoor, King, Noack, King and Tadmor2005).

In case of the CNM, the modelled autocorrelation function $\hat {R}$ is based on the modelled velocity field (2.21). In case of the CMM, the time integration quickly leads to the average flow and is not indicative of the range of possible initial conditions. Hence, $K$ trajectories are considered starting with $p_k=1$ for each cluster $k$, or, equivalently, $\boldsymbol {p}^{\circ k} (t=0)= [ \delta _{1k}, \ldots , \delta _{Kk} ]^{\mathrm {T}}$. These cluster-specific autocorrelation functions are weighted with the cluster probability $p_i^{\infty }$

(2.23)\begin{equation} \hat{R} (\tau) :=\sum_{k=1}^K p_i^{\infty} \int_{0}^T \int_{\varOmega} \boldsymbol{u}^{\circ k} (\boldsymbol{x}, t) \boldsymbol{\cdot} \boldsymbol{u}^{\circ k} (\boldsymbol{x}, t+\tau ) \,\textrm{d}\boldsymbol{x}\,\textrm{d}t, \quad \tau \in[0, T), \end{equation}

where $\boldsymbol {u}^{\circ k}$ denotes the CMM-predicted velocity field starting in cluster $k$.

2.6. Lorenz system as an illustrating example

Following the original CMM paper by Kaiser et al. (Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014), the CNM is illustrated for the celebrated Lorenz (Reference Lorenz1963) system, arguably the first demonstration of a chaotic dynamics in low-dimensional dynamics. The Lorenz system is a three-dimensional autonomous system of nonlinear ordinary differential equations. The derivation was inspired by a Galerkin model of Rayleigh–Bénard convection, but typically selected parameters clearly exceed the range of model validity (Sparrow Reference Sparrow1982). The system features non-periodic, deterministic, dissipative dynamics associated with exponential divergence and convergence to a fractal strange attractor. The three coupled nonlinear differential equations read

(2.24a)\begin{gather} \frac{\mathrm{d} x}{\mathrm{d} t} =\sigma(y-x), \end{gather}

(2.24b)\begin{gather}\frac{\mathrm{d} y}{\mathrm{d} t} =x(r-z)-y, \end{gather}

(2.24c)\begin{gather}\frac{\mathrm{d} z}{\mathrm{d} t} =x y-b z, \end{gather}

with the system parameters $\sigma =10$, $b=8/3$ and $r=28$. For these parameters, there are three unstable fixed points at $(0,0,0)$ and $(\pm \sqrt {72},\pm \sqrt {72},27)$, denoted by $F^+$ and $F^-$, respectively. The attractor of the Lorenz system resembles two butterfly wings around $F^+$ and $F^-$ in phase space. The trajectory typically oscillates for several periods with increasing amplitude around a fixed point ($F^+$ or $F^-$) before it moves to the other wing. The number of revolutions made on either side varies unpredictably from one cycle to the next.

The Lorenz equations (2.24) are solved employing an explicit fourth-order Runge–Kutta scheme with an initial condition on the attractor. The time-resolved snapshot data $\boldsymbol {x}(t_{m})$ with $\boldsymbol x=[x,y,z]^{\textrm {T}}$ are collected at a sampling time step $\Delta t=0.005$ corresponding roughly to one thousandth of a typical oscillation period. The k-means++ algorithm partitions $M=1\,000\,000$ snapshots into $K=10$ clusters. Figure 6(a) displays a phase portrait of the corresponding clusters. The snapshots associated with one cluster are highlighted by the same colour. The $10$ clusters feature three different subsets: the transition clusters $k=1,2$ between two butterfly wings, the $F^-$ wing related cluster $k=3,4,5,6$ and clusters $k=7,8,9,10$ associated with the $F^+$ wing. The wing-related cluster groups represent approximately $90^{\circ }$ phase bins and do not resolve the amplitude. Evidently, the 10 clusters are coarse representations of the state.

Figure 6. Cluster-based network model for the Lorenz attractor. (a) Cluster partitioning: the centroids are displayed as coloured solid circles. A trajectory is illustrated by a black curve. The dots on this trajectory are coloured according to their cluster affiliation. The clusters $k=3,4,5,6$ oscillate around the fixed point $F^-$ and clusters $k=7,8,9,10$ around $F^+$. The clusters $k=1,2$ connect both ‘ears’ of the Lorenz attractor. (b) The trajectory of the CNM (red dashed line). The centroids represent the network nodes and edges represent possible transitions. Here, trajectory of CNM is obtained by a spline interpolation through the visited centroids.

In the following, the dynamics is resolved by the network model of § 2.3. The $10$ centroids are considered as nodes in the network. The transition between these centroids define directed edges characterized by direct transition matrix ${\boldsymbol{\mathsf{Q}}}$ and the flight times ${\boldsymbol{\mathsf{T}}}$. The connectivity is described by the adjacency matrix $H({\boldsymbol{\mathsf{Q}}})$ where $H$ denotes the Heaviside function: non-vanishing elements of ${\boldsymbol{\mathsf{Q}}}$ are replaced by unity (Newman Reference Newman2010). Figure 7, displays the DTM ${\boldsymbol{\mathsf{Q}}}$ (figure 7a) and associated transition time matrix ${\boldsymbol{\mathsf{T}}}$ (figure 7b). The matrices reveal three distinct cluster groups consistent with the phase diagram of figure 6. Clusters $1$ and $2$ allow transitions to $3$ and $10$, i.e. the $F^-$ and $F^+$ wing, respectively, and have been called flipper clusters by Kaiser et al. (Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014). The cluster transition sequence $3 \to 4 \to 5 \to 6 \to 1 \to 2 \to 3$ is the dominant cyclic group associated with the $F^-$ wing. Another cyclic group skips cluster $2$: $3 \to 4 \to 5 \to 6 \to 1 \to 3$. A cyclic group through the $F^+$ wing reads $10 \to 9 \to 8 \to 7 \to 1 \to 10$. A longer sequence includes the second cluster: $10 \to 9 \to 8 \to 7 \to 1 \to 2 \to 10$. The transition times in the wing centroids are noticeably smaller than the passage through the flipper clusters.

Figure 7. Cluster-based network model of the Lorenz attractor. (a) Direct transition matrix $(Q_{ij})$. (b) Transition time $(T_{ij})$.

Figure 8 compares the asymptotic population $\boldsymbol {p}^{\infty }$ predicted by CNM with the population from a long-term simulation. The CNM statistics are based on $20\,000$ transitions while the integration of the Lorenz equations is performed over 5000 time units. Both statistics correspond to the roughly 800 periods found to be sufficient for an accurate statistics. The relative error of the CNM is no more than approximately 10 %. This error does not decrease with much larger integration times, but is linked to the coarse-graining of the state to Voronoi cells around the centroids. The assumed constant transition time $\tau _{ij}$ for all trajectories from cluster $j$ to cluster $i$ is a crude assumption. In fact, the transition times can vary by a large factor and can thus give rise to significant systematic errors. A more accurate transition model may, for instance, include earlier transitions for a more realistic representation of the trajectory. Intriguingly, the CMM has an error of only 0.5 % which is one order of magnitude lower. Due to the stroboscopic monitoring of the CMM states, no estimates of the transition times are required and one source of systematic errors is excluded by construction.

Figure 8. Cluster probability distribution of the Lorenz system (solid rectangles) and the corresponding cluster-based network model (open rectangles). The model results are based on $20\,000$ transitions.

A distinguishing feature of the CNM is the resolution of the temporal dynamics illustrated in figure 9. The evolution of the model-based trajectory is hardly distinguishable from the one obtained by numerical integration. Smoothness of the CNM trajectory has been achieved by splines connecting the states between two consecutive centroid visits. Yet, the oscillatory amplitude growth in both wings cannot be resolved with this low cluster-based resolution. The CNM can only follow the simulations for a short time period, as nearby trajectories exponentially diverge with Lyapunov exponent 2.16 (Wolf et al. Reference Wolf, Swift, Swinney and Vastano1985) and the initial separation in each cluster is already large. Yet, the fluctuation amplitude, frequency content and bi-modality are well reproduced. A CNM with $K=100$ clusters yields more realistic dynamics but require orders of magnitude more simulation data. At the least-order extreme, a CNM with two or three clusters only coarsely resolves the transitions between both ears of the Lorenz attractor, not the growing oscillations in each ear.

Figure 9. Evolution of the Lorenz system $x$, $y$, $z$, $t\in [0,20]$ from integrating the dynamical system (black solid line) and from the prediction of the cluster-based network model (red dashed line).

Finally, the autocorrelation of the simulation (black solid curve) and the CNM (red dashed curve) are presented for aggregate comparison in figure 10. CNM roughly reproduces the fast oscillatory decay of the autocorrelation function in the first five periods.

Figure 10. Autocorrelation function of the Lorenz system (black solid line) and the cluster-based network model (red dashed line).

3.1. Flow configuration and direct numerical simulation

The two-dimensional incompressible mixing layer and a velocity ratio of $3:1$ is considered as the test plant in this paper. The velocity ratio is a common choice in the literature (Comte, Silvestrini & Bégou Reference Comte, Silvestrini and Bégou1998; Noack et al. Reference Noack, Papas and Monkewitz2005; Kaiser et al. Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014). The high- and low-speed streams have velocities $U_1$ and $U_2$, respectively. The convection velocity $U_c$ of coherent structure is well approximated by the average velocity (Monkewitz Reference Monkewitz1988)

(3.1)\begin{equation} U_{c}=\frac{U_{1}+U_{2}}{2}. \end{equation}

The initial vorticity thickness is denoted by $\delta _0$. The Newtonian fluid is characterized by the density $\rho$ and kinematic viscosity $\nu$. The flow characteristics are described by the Reynolds number based on the convection velocity $Re = U_c \delta _0 / \nu$ and velocity ratio. In the sequel, all quantities are assumed to be non-dimensionalized with the initial vorticity thickness $\delta _0$, the high-speed velocity $U_1$ and the density $\rho$. The Reynolds number is set to 200.

The flow is described in a Cartesian coordinate system $(x,y)$ with the origin at the maximum gradient location of the inlet profile. The $x$-axis points in the streamwise direction and the $y$-axis points in the direction of the high-speed stream. The velocity components in the $x$- and $y$-directions are denoted by $u$ and $v$, respectively.

Figure 11 illustrates the rectangular computational domain

(3.2)\begin{equation} \varOmega : = \{ ( x , y ) \in \mathcal{R}^2 : 0 \leq x \leq 80 \wedge \vert y \vert \leq 15 \} \end{equation}

with 10 237 nodes and 2248 triangular elements. The location vector is denoted by $\boldsymbol {x}=(x,y)$. Similarly, the velocity vector is denoted by $\boldsymbol {u}=(u,v)$. The inlet velocity profile reads

(3.3)\begin{equation} u=2+\tanh \left( \frac{2y}{\delta_{0}} \right), \quad v=0, \quad \text{where} \ \delta_{0}=1. \end{equation}

The Kelvin–Helmholtz vortices are triggered at the inlet by a stochastic perturbation of the $u$-component for $y \in [-2,2]$ with a standard deviation of $0.01 U_c$.

Figure 11. Numerical simulation sketch of incompressible mixing layer. An unperturbed tanh velocity profile $u(y) =2+\tanh (2y)$ is patched in the inlet of a rectangular domain. Time-averaged streamwise velocity profiles separated by $\Delta x$=10 are visualized by blue lines. The red dashed curves mark the mixing-layer thickness. We chose the 90 % thickness of the profile starting with the average velocity, i.e. $u=2.9$ and $u=1.1$.

The streamwise extent of the domain is 80. This corresponds to a downwash time of 40 given the convection velocity of 2. This is the minimum time for the transient time as all initial interior vortices will leave the domain. A simulation over 400 convective units corresponds to 10 downwash times. This period is found to be sufficient for a good statistics of the mean value and fluctuation level. One simulation yields $M=10\,000$ velocity snapshots $\boldsymbol {u}^{m}(\boldsymbol {x})=\boldsymbol {u}(\boldsymbol {x},t^m)$, where the sampling times $t^m = 0.04\ \textrm {m}$ start with $t=0$ in the converged post-transient phase. The sampling frequency $25$ is two orders of magnitude larger than the dominant shear-layer frequency of $f=0.1075$ in the most active downstream region.

An in-house direct numerical simulation solver was employed to simulate the incompressible mixing layer. This solver is based on the finite-element method with third-order Taylor–Hood elements and implicit third-order time integration. The solver has been used for numerous configurations, like the cylinder wake (Noack et al. Reference Noack, Stankiewicz, Morzyński and Schmid2016), the mixing layer (Shaqarin, Noack & Morzyński Reference Shaqarin, Noack and Morzyński2018) and the fluidic pinball (Ishar et al. Reference Ishar, Kaiser, Morzynski, Albers, Meysonnat, Schröder and Noack2019), to name only a few.

3.2. Flow features

The incompressible mixing layer exhibits two typical behaviours. First, the initial dynamics is characterized by the roll-up of vorticity originating from the Kelvin–Helmholtz (K–H) instability (see figure 12b i). Second, these vortices pair further downstream, as can be seen at the outlet region of figure 12(b ii). This vortex pairing contributes to the mixing layer growth. The location of vortex pairing may change significantly in time. Upstream (downstream) vortex pairing is associated with high (low) fluctuation energy.

Figure 12. Mixing-layer simulation. (a) Energy fluctuation over time with the maximum marked by a blue square and minimum by a red square. (b) Vorticity fields associated with the maximum and minimum fluctuation energies. The minimum (b i) corresponds to a K–H vortex at $t=167.60$, The maximum (b ii) features vortex pairing at $t=195.92$. The curves represent the isolines of vorticity. Higher values corresponds to darker green areas.

The time-averaged velocity field in figure 11 shows the mixing layer growth. The velocity thickness is visualized by a red dashed line and is defined as the distance between transverse locations where the mean streamwise velocity was equal to $U_{1}-0.05 \Delta U$ and $U_{2}+0.05 \Delta U$. The mixing-layer thickness increases significantly between $x=30$ and $x=60$. Here, vortex pairing leads to this thickness increase.

The temporal dynamics may be inferred from the evolution of the fluctuation energy in figure 12(a). The fluctuations indicate narrow bandwidth oscillatory behaviour. More refined insights may be gained from the correlation function between the flows at time $t$ and $t^{\prime }$,

(3.4)\begin{equation} C(t,t^{\prime}) = \int_{\varOmega} \boldsymbol{u}^{\prime} \left ( \boldsymbol{x},t \right ) \boldsymbol{\cdot} \boldsymbol{u}^{\prime} \left ( \boldsymbol{x},t^{\prime} \right )\, \textrm{d}\boldsymbol{x}. \end{equation}

Figure 13 illustrates the autocorrelation matrix for $t,t^{\prime } \in [0,400]$. The fluctuation energy of figure 12(a) is quantified in the diagonal, $\mathcal {K}(t) = C(t,t)/2$. The wavy pattern indicates oscillatory coherent structures. The changes from pure periodicity are caused by vortex pairing at a large range of streamwise locations.

Figure 13. Autocorrelation matrix (3.4) of the mixing layer for $t ,t^{\prime } \in [0,400]$. The value is presented in the colour bar. The plot is based on 401 snapshots collected at uniform time steps $\Delta t=1$.

3.3. Clustering

Both considered reduced-order models are based on the direct numerical simulation of the two-dimensional incompressible mixing layer described in § 3.1; $M=10\,000$ velocity field snapshots of the post-transient phase are sampled with a time step $\Delta t = 0.04$.

The computational load of clustering is significantly reduced by an effectively lossless POD compression detailed in appendix A. In fact, all operations are performed on the POD amplitude vector $\boldsymbol {a}= [a_{1},a_{2},\ldots ,a_{N}]^{\textrm {T}}$ instead of the snapshots.

The $M$ snapshots are clustered with the k-means++ algorithm into $K=10$ centroids. This number is small enough to allow for the physical interpretation of all centroids and all transitions but large enough for a meaningful reduced-order model. Figure 14 illustrates the transverse velocity fluctuation of the centroids. The first six centroids show the streamwise convection of K–H vortices, while the next four centroids resolve a vortex pairing (VP) event. In centroid $7$, two vortices merge at the beginning of the vortex chain. In the following three centroids, the merging is completed and leads to a large vortex. Note that the VP centroids $k=7,8,9,10$ have pronounced vortices at a similar position as the K–H centroids $k=4,5,6$, respectively. The structures of the K–H and VP centroids are noticeably different. The main vortices of the K–H centroids are elliptical and the major axis is rotated in clockwise direction, i.e. the upper part of the vortices follow the faster stream. In contrast, the main elliptical vortices of VP centroids are rotated in mathematically positive direction, i.e. the upper part of these vortices move upstream with respect to their centre.

Figure 14. The cluster centroids $\boldsymbol {c}_{k}$, $k=1, \ldots , 10$, of the mixing layer. The transverse velocity fluctuation is depicted with contour lines. Red and blue regions mark positive and negative values.

The centroids represent characteristic stages in the mixing-layer dynamics as can be elucidated in a proximity map. This map reflects the configuration matrix ${\boldsymbol{\mathsf{D}}} = (D_{ij}) \in $ $\mathcal {R}^{K \times K}$ comprising the distance between two centroids

(3.5)\begin{equation} D_{i j}^{c} :=\left\|\boldsymbol{c}_{i}-\boldsymbol{c}_{j}\right\|_{\varOmega},\quad i,j=1,2,\ldots,K. \end{equation}

Following Kaiser et al. (Reference Kaiser, Noack, Cordier, Spohn, Segond, Abel, Daviller, Östh, Krajnović and Niven2014), the proximity map is used to represent the configuration matrix ${\boldsymbol{\mathsf{D}}}$ in a two-dimensional feature space $\boldsymbol {\gamma } \in \mathcal {R}^2$ optimally preserving the relative distances. The proximity map employs classical multidimensional scaling (Mardia, Kent & Bibby Reference Mardia, Kent and Bibby1979). Figure 15(a) displays centroids close to a circle which is characteristic of vortex shedding.

Figure 15. Cluster-based Markov model of the mixing layer. (a) Proximity map of centroids. Each centroid is marked by a solid coloured circle. The colour denotes the relative energy content (see colour bar on top). Unity corresponds to the average value. (b) Transition matrix. The probability value is displayed by the background colour and the radius of the corresponding circle.

3.4. Markov model

The temporal mixing-layer evolution is characterized by the cluster transition matrix ${\boldsymbol{\mathsf{P}}}$ illustrated in figure 15(b); $P_{ij}$ represents the probability of moving from cluster $j$ to $i$ in one forward time step. Here, we choose a time step $\Delta t^c=T/10=1$ where the $T=10$ is the dominant period of the evolved mixing layer.

The cluster transition matrix reveals two cyclic groups. The first group $1 \to 2 \to 3 \to 4 \to 5 \to 6 \to 7 \to 1$ is consistent with the convection process of the K–H vortex shedding observed in the centroid visualization. This periodic process corresponds to a nearly uniform clockwise rotation in the proximity map. The second cyclic group $8 \to 9 \to 10 \to 7 \to 1 \to 2\to 8$ comprises VP centroids $k=8,9,10$ and shares two centroids with the K–H regime. This dynamics also leads to a nearly uniform clockwise rotation in the feature space. There are also transitions from the VP to K–H regime, e.g. $8 \to 4$, $9 \to 5$, $10 \to 6$ and $10 \to 7$ and in the opposite direction. All these transitions are between similar centroids of both groups. From the cluster index the orientation of the main elliptical vortices can be inferred. For $k \le 6$ ($k\ge 7$), the upper part of the vortices are displaced in (against) the direction of the flow with respect to their centres.

The evolution of the cluster population vector $\boldsymbol {p}^l$ at $t=l \Delta t^c$ is investigated by iterating equation (2.8). Figure 16 compares the probability distribution of DNS data and the model-based asymptotic vector $\boldsymbol {p}^{\infty }$. The agreement is astonishingly good for such a low-order model. The probability vector converges quickly to a unique, stationary probability distribution near $t=20$.

Figure 16. Cluster probability distribution of the mixing layer from the DNS data and the cluster-based Markov model. Solid rectangles denote the probability $q_{k}$ from the DNS data. Open rectangles represent asymptotic values from the CMM after $l=35$ iterations.

In figure 17, the dynamics of CMM is illustrated for the first cluster probability $p_1$ and the first POD mode amplitude $a_1$ inferred from the flow state (2.19). Starting point is direct numerical simulation starting at $t=0$ close to the first cluster $\boldsymbol {c}_1$ which corresponds to the probability vector $\boldsymbol {p}=[1,0,0,0,0,0,0,0,0,0]^{\textrm {T}}$. The probability and POD mode amplitude of CMM show a convergence after around $l=35$ iterations or, equivalently, $t\approx 35$. The solid horizon line denotes $q_1$, i.e. the population of the first cluster from DNS data. The POD mode amplitude $a_1$ performs three strongly damped oscillations before vanishing.

Figure 17. Dynamics of the mixing layer from DNS and the cluster-based Markov model. $(a)$ Probability evolution for DNS (black solid line) and for CMM (red dashed line). The probability of the first cluster $p_{1}$ quickly converges to $p^{\infty }_i$ around $t=35$. The corresponding DNS value is $0.1325$ and is represented by a horizontal line. $(b)$ The evolution of the first POD mode amplitude $a_{1}$ for DNS (black solid line) and CMM (red dashed line). $(c)$ The autocorrelation function for DNS (black solid line) and for CMM (red dashed line).

Figure 17(c) shows an oscillating quickly decaying autocorrelation function of CMM which is consistent with the observations for $a_1$ and $p_1$. In contrast, the autocorrelation function associated with the DNS keeps oscillating around with an amplitude around 50 % of the average fluctuation level. This level indicates that half of the fluctuation energy resides in repeating oscillatory flow structures while the other half is of non-repeating stochastic nature.

3.5. Network model

In this section, a CNM is developed using the same snapshot data and same centroids. Starting point for the dynamic network is the cluster-affiliation function $k(t)$. Following § 2.3, the direct cluster transition matrix ${\boldsymbol{\mathsf{Q}}}$ with associated average transition times ${\boldsymbol{\mathsf{T}}}$ are derived. Figure 18 illustrates both matrices. These matrices have the almost same structure as the Markov model except for the diagonal elements which are vanishing by design. In other words,

(3.6a)\begin{gather} Q_{ii}=T_{ii}=0, \quad \forall i \in \{1,\ldots,K\} , \end{gather}

(3.6b)\begin{gather}H[P_{ij}]=H[Q_{ij}]= H[T_{ij}], \quad \forall i, j \in \{1,\ldots,K\} \wedge i \not = j , \end{gather}

with $H$ being again the Heaviside function. Vanishing diagonal elements (3.6a) arise from the requirement of non-trivial transitions. Theoretically, the trajectory may terminate in a cluster, like in a stable fixed point of a linear dynamical system. This case is not compatible with the goal to model a well-resolved non-trivial attractor and shall be ignored in this study. Equation (3.6b) requires a sufficiently small time step of the CMM. Otherwise, the stroboscopic view on the trajectory may miss a crossing of an intermediate cluster. This happens with the transition from $1 \to 2 \to 3$ in one CMM time step $\Delta t^c$. Hence, $P_{31} \not = 0$ while $Q_{31} = 0$. However, this is a rare event, as indicated by the small value of $P_{31}$.

Figure 18. Dynamics of the cluster-based network model for the mixing layer. (a) Direct transition matrix. (b) Averaged transition time. The non-vanishing values are denoted by the circle radius and the colour code from the bottom caption.

An inspection of ${\boldsymbol{\mathsf{T}}}$ reveals that the transition time between K–H and VP centroids is relatively small. This is consistent with the closeness of the corresponding centroids in the proximity map (figure 15a). An exception is the transition between K–H centroid 2 to VP centroid 8 which are well separated in the proximity map. Intriguingly, the transitions within the K–H and VP regimes are also strongly correlated with the distances depicted in the proximity map. For instance, the smallest (largest) inner-regime transition from centroid $6$ to $7$ ($8$ to $9$) is associated with a small (large) distance in the proximity map. The physical interpretation of the cycle-to-cycle variations of the CMM persist for CNM.

In the following, the temporal dynamics of CNM is investigated based on the identified centroids $\boldsymbol {c}_{k}$, the description of their connectivity DTM ${\boldsymbol{\mathsf{Q}}}$ and their flight times ${\boldsymbol{\mathsf{T}}}$. Like a POD model, CNM is a grey-box model resolving the temporal dynamics and the associated coherent structures. We choose cluster $k=1$ as initial condition for DNS and for the CNM and integrate over $l=20\,000$ transitions. In figure 19, the asymptotic cluster population $\boldsymbol {p}^{\infty }$ from (2.17) is compared with $\boldsymbol {q}$ from the DNS. The discrepancies of few per cent seem expectable and tolerable for a 10-cluster model. This difference is not cured by increasing the amount of transition data in CNM. Intriguingly, the probability distribution of the CMM displayed in figure 16 is significantly more accurate. This behaviour can be linked to the simple transition time estimate which employs one single average value for a large range of observed transition times. We have developed more refined and more accurate transition time estimates leading to much better agreements of the cluster probability distributions. The price is increased complexity of the CNM which we deemed not helpful for our first publication.

Figure 19. Cluster probability distribution of the mixing layer from DNS (solid rectangles) and cluster-based network model (open rectangles). The modelled values are obtained from simulating $20\,000$ clusters transitions.

Figure 20 shows the evolution of the first two POD mode amplitudes (red dashed curve). The CNM tracks well the amplitude and phase of the DNS over $100$ time units. As for the Lorenz system, the temporal evolution is smoothed by a spline and does not use the non-smooth uniform motion between two consecutive centroid visits.

Figure 20. Evolution of mode amplitudes $a_{1}$, $a_{2}$ for the mixing layer $t\in [0,100]$. The curves correspond to DNS (black solid line) and the cluster-based network model (red dashed line).

Figure 21 compares the autocorrelation function of the CNM and the DNS. We intentionally do not normalize this function to reveal the resolved fluctuation level at vanishing time delay. As expected, the model-based fluctuation level is significantly lower than the DNS value. This difference is quantified by the unresolved inner-cluster variance. Intriguingly, CNM and DNS functions become already similar after half a period. The asymptotic fluctuation level represents coherent structures which are well resolved by the chosen centroids and serve as coarse-grained recurrence points of the DNS. Due to the dominant oscillatory dynamics, the autocorrelation does not vanish with increasing time. The good reproduction autocorrelation function is a posteriori justification for the chosen cluster number.

Figure 21. Autocorrelation function of the mixing layer for DNS (black solid line) and cluster-based network model (red dashed line).

4.1. Flow configuration and large-eddy simulation

In this section, the actuated turbulent boundary layer configuration for skin-friction reduction is detailed. In particular, the actuation mechanism is presented, and the numerical set-up is described. For more details, the reader is referred to Albers et al. (Reference Albers, Meysonnat and Schröder2019) and Fernex et al. (Reference Fernex, Semaan, Albers, Meysonnat, Schröder and Noack2020). The fluid flow is described in a Cartesian frame of reference where the streamwise, wall-normal and spanwise coordinates are denoted by $\boldsymbol {x} = (x,y,z)$ and the velocity components by $\boldsymbol {u} = (u,v,w)$. The Mach number is set to Ma $= 0.1$ corresponding to a nearly incompressible flow. An illustration of the rectangular physical domain is shown in figure 22. A momentum thickness of $\theta = 1$ at $x_0$ is achieved such that the momentum thickness based Reynolds number is $Re_{\theta } = 1000$ at $x_0$. The domain length and height in the streamwise and wall-normal directions are $L_x = 190\,\theta$ and $L_y = 105\,\theta$. In the spanwise direction, different domain widths $L_z \in [21.65\,\theta , 108.25\,\theta ]$ are used to simulate different actuation wavelengths.

Figure 22. Overview of the physical domain of the actuated turbulent boundary layer flow, where $L_x, L_y,$ and $L_z$ are the domain dimensions in the Cartesian directions, $\lambda$ is the wavelength of the spanwise travelling wave and $x_0$ marks the actuation onset. The shaded red surface $A_{surf}$ marks the integration area of the wall-shear stress $\tau _w$.

At the domain inlet, a synthetic turbulence generation method is applied to generate a natural turbulent boundary layer flow after a transition length of 2–4 boundary layer thicknesses (Roidl, Meinke & Schröder Reference Roidl, Meinke and Schröder2013). Characteristic boundary conditions are used at the domain exit and a no-slip wall boundary condition is enforced at the lower domain boundary for the unactuated and actuated wall.

The actuation is performed by a transverse travelling wave on the surface. The corresponding wall motion is prescribed by the space- and time-dependent function

(4.1)\begin{equation} y_{{wall}}^+(z^+,t^+) = A^+ \cos\left( \frac{2{\rm \pi}}{\lambda^+}z^+ - \frac{2{\rm \pi}}{T^+}t^+ \right), \end{equation}

in the interval $-5 \leq x/\theta \leq 140$. The quantities $\lambda ^+$, $T^+$ and $A^+$ denote the wavelength, period and amplitude in inner coordinates, i.e. the parameters are scaled by the viscosity $\nu$ and the friction velocity of the unactuated reference case $u^n_{\tau }$. In the area just upstream and downstream of the wave actuation region, a spatial transition is used from a flat plate to an actuated plate and vice versa (Albers et al. Reference Albers, Meysonnat and Schröder2019). In total, 38 actuation configurations with wavelength $\lambda ^+ \in [200,500,3000]$, period $T^+ \in [20,120]$ and amplitude $A^+ \in [10,78]$ are simulated. In the current study, we model one test case with $\lambda ^+=1000$, $T^+=120$ and $A^+ =40$ which yields the largest drag reduction of 3 % found at that wavelength. These actuation parameters correspond to case N36 in table 3 of Ishar et al. (Reference Ishar, Kaiser, Morzynski, Albers, Meysonnat, Schröder and Noack2019) and in table 2 of Albers et al. (Reference Albers, Meysonnat, Fernex, Semann, Noack and Schröder2020).

The physical domain is discretized by a structured block-type mesh with a resolution of $\Delta x^+ = 12.0$ in the streamwise and $\Delta z^+ = 4.0$ in the spanwise direction. In the wall-normal direction, a resolution of $\Delta _y^+|_{{wall}} = 1.0$ at the wall is used with gradual coarsening away from the wall. Depending on the domain width, the meshes consist of 24–120 million cells.

The actuated flat plate turbulent boundary layer flow is governed by the unsteady compressible Navier–Stokes equations in the arbitrary Lagrangian–Eulerian formulation for time-dependent domains. A second-order accurate finite-volume approximation of the governing equations is used in which the convective fluxes are computed by the advection upstream splitting method (AUSM) and time integration is performed via a 5-stage Runge–Kutta scheme. The smallest dissipative scales are implicitly modelled through the numerical dissipation of the AUSM scheme. This monotonically integrated large-eddy simulation approach (Boris et al. Reference Boris, Grinstein, Oran and Kolbe1992) is capable of accurately capturing all physics of the resolved scales (Meinke et al. Reference Meinke, Schröder, Krause and Rister2002).

The actuated simulations are initialized by the solution from the unactuated reference case and the temporal transition from the flat plate to the actuated wall is initiated. When a converged state of the friction drag is obtained, statistics are collected for $t U_{\infty } / \theta = 1250$ convective times.

The actuation effects on the near-wall flow features are illustrated in figure 23, which shows contours of the streamwise velocity fluctuation of a reference natural (figure 23a) and the actuated case (figure 23b). The intensity of the near-wall streaks, which are known to contribute to skin friction, are observed to diminish with the actuation.

Figure 23. Contours of the random streamwise velocity fluctuations in the near-wall region of (a) a non-actuated reference case and (b) the actuated case. The actuation diminishes the near-wall streak intensity. The surface deformation is indicated at the bottom.

4.2. Clustering

Similar to the mixing layer, the clustering of the actuated boundary layer large-eddy simulation (LES) snapshots is performed using a lossless POD compression. Again, this compression dramatically reduces the computational load of clustering. Here, we perform the POD and the clustering on all 38 test cases simultaneously. Employing this enlarged set of POD modes yields a richer, and thus more accurate, dynamical representation of the individual test cases and allows for a direct comparison of different actuations (Ishar et al. Reference Ishar, Kaiser, Morzynski, Albers, Meysonnat, Schröder and Noack2019). Concatenating all configurations results in $M=15\,873$ snapshots sampled at $\Delta t = 0.94$ time units.

Following Ishar et al. (Reference Ishar, Kaiser, Morzynski, Albers, Meysonnat, Schröder and Noack2019), the $M$ snapshots are clustered with the k-means++ algorithm into 50 centroids, corresponding to $K=10$ centroids populated by the investigated actuation. It is worth noting that increasing $K$ significantly, say $K=100$, uncovers centroids with smaller length scale features associated with broadband turbulence of the boundary layer. In this study, we purposely choose to focus on the main energy-containing dynamics and thus limit the number of centroids to $K=10$. Figure 24 presents four centroid distributions of the test case with $\lambda ^+=1000$, $T^+ =120$ and $A^+ =40$. As the figure shows, the centroids have similar spatial distributions and are phase shifted with respect to one another. Such behaviour is consistent with a limit-cycle dynamics, indicating partial lock-on of the boundary layer dynamics to the periodic surface actuation. This lock-on phenomenon is sometimes associated with aerodynamic gains or losses depending on the targeted flow instability. It is synonymous with synchronization, and has been repeatedly investigated for drag reduction problems (Barros et al. Reference Barros, Borée, Noack and Spohn2016; Taira & Nakao Reference Taira and Nakao2018; Herrmann et al. Reference Herrmann, Philipp, Richard and Brunton2020). Similar to these studies, a lower actuation threshold with sufficient authority is required to synchronize the flow.

Figure 24. The cluster centroids $\boldsymbol {c}_{k}, k= 2, 4, 9, 10$ of the actuated boundary layer. The cluster numbers are denoted in the state space of figure 25. The iso-surfaces correspond to constant wall-normal velocity of $V^+={\pm }0.04$. Red and blue regions mark positive and negative values.

The dynamics is well represented in the state space illustrated in figure 25, which is spanned by the first three POD mode coefficients. The cluster centroids are displayed as black solid circles and their indices are labelled. The snapshots are coloured according to their cluster affiliation. Similar to the shear layer, the dynamics of the actuated boundary layer appears to be driven by two physical phenomena: a cyclic behaviour synchronized with the surface actuation, and a quasi-stochastic component that forces the limit cycle to experience cycle-to-cycle variations (Cao et al. Reference Cao, Kaiser, Boreé, Noack, Thomas and Guilan2014). The latter phenomenon is associated with broadband turbulence of the boundary layer.

Figure 25. The state space is spanned by the first three mode coefficients of the lossless proper orthogonal decomposition of the LES data. The cluster centroids are displayed as black solid circles and the snapshots are coloured according to their cluster affiliation.

4.3. Network model

The CNM is generated based on the direct transition matrix and the averaged transition time matrix, which are illustrated in figure 26. We reiterate the vanishing diagonal elements of both matrices, i.e. $Q_{ii}=T_{ii}=0$, which are a result of enforcing non-trivial transitions. The direct transition matrix (cf. figure 26a) shows both the dominant transition probability of subsequent centroids being associated with the limit-cycle behaviour, and the wandering dynamics from the remaining transitions. The transition time matrix between the centroids (cf. figure 26b) reflects the same behaviour, and exhibits a quasi-constant transition time for limit-cycle ‘flight times’ and diverse transition times for the wandering effect.

Figure 26. Dynamics of the cluster-based network model for the actuated boundary layer. (a) Direct transition matrix. (b) Averaged transition time. The non-vanishing values of the matrix elements are proportional to the circle radius and can be inferred from the colour code from the bottom caption.

Figure 27 compares the probability distribution of LES data and the model-based asymptotic vector $\boldsymbol {p}^{\infty }$. Again, we choose cluster $k = 1$ as the initial condition for LES and for the CNM and integrate over $l = 106$ transitions, which corresponds to a similar time range as that of the snapshots. The agreement between the two distributions is good.

Figure 27. Probability distribution of the actuated boundary layer from large-eddy simulation (solid rectangle) and cluster-based network model (open rectangle). The CNM values are obtained from simulating $l=106$ cluster transitions.

The model performance is assessed against the reference LES results. Figure 28 shows the evolution of the first four POD mode amplitudes (red dashed curve). The dominance of the first two POD modes compared to the subsequent modes is expected for the current quasi-synchronous actuated flow. Similar to the previously presented results, the temporal evolution is smoothed with a spline. As the figure shows, CNM agrees very well with the amplitude and phase of the LES reference data over the entire approximately $400$ time units.

Figure 28. Evolution of mode amplitudes $a_{1}$–$a_{4}$ for the actuated boundary layer $t\in [0,400]$. The curves correspond to LES (black solid line) and cluster-based network model (red dashed line).

The agreement between the model and the reference data is further corroborated by comparing the autocorrelation function. Figure 29 displays the autocorrelation function of the CNM and the LES. As with the mixing layer, the model-based fluctuation level at vanishing time delay is lower than the LES value but becomes similar to oscillation level for an arbitrary larger time horizon. This large representation error at $\tau =0$ relates to the unresolved inner-cluster variance. Yet, the centroids adequately resolve the periodic flow response of the flow to the periodic surface actuation.

Figure 29. Autocorrelation function of the actuated turbulent boundary layer for LES (black solid line) and cluster-based network model (red dashed line).