2.1. PCF

We consider a PCF of an incompressible fluid, where the two walls move in opposite directions with streamwise velocity, ${\pm }U_w$. The wall-normal distance between the two walls is given by $2h$. The kinematic viscosity of the fluid is $\nu$, and the Reynolds number is defined by $Re = U_w h/\nu$. The streamwise, wall-normal and spanwise coordinates are made dimensionless with $h$, and they are denoted by $\boldsymbol {x} = (x,y,z)$. The two walls are set to be located at $y=\pm 1$. The equations of motion are then given by

(2.1)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} ={-}(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} - y\frac{\partial \boldsymbol{u}}{\partial x} - {v} \boldsymbol{e}_x - \boldsymbol{\nabla} p + \frac{1}{Re}\nabla^{2}\boldsymbol{u}, \end{equation}

where $t$ is the time, $\boldsymbol {u}=({u},{v},{w})$ is the perturbation velocity from the laminar base flow, $\boldsymbol {U}=(y,0,0)$, $p$ the pressure and $\boldsymbol {e}_x$ the unit vector in the streamwise direction.

To build a reduced-order model, we first perform a direct numerical simulation (DNS) confined to a minimal flow unit (Jiménez & Moin Reference Jiménez and Moin1991; Hamilton et al. Reference Hamilton, Kim and Waleffe1995). Following Hamilton et al. (Reference Hamilton, Kim and Waleffe1995), the computational domain of $(L_x,L_y,L_z)=(1.75{\rm \pi},2,1.2{\rm \pi} )$ is considered at $Re=400$. The simulation was performed using channelflow 2.0, an open-source DNS code (https://www.channelflow.ch/). The code uses the Fourier–Galerkin discretisation in the streamwise and spanwise directions and the Chebyshev-tau discretisation in the wall-normal direction. The number of grid points in each spatial direction is given by ($N_x,N_y,N_z) = (16,33,16$) (after dealiasing). A third-order semi-implicit backward difference formula (SDBF3) is used for the time-integration scheme. The simulation has been performed by setting a zero pressure gradient. The domain size normalised by viscous inner units (denoted by the superscript $(\cdot )^{+}$) is obtained as $(L_x^{+},{L_y^{+}},L_z^{+}) \simeq (186,68,127)$, in good agreement with that in Hamilton et al. (Reference Hamilton, Kim and Waleffe1995).

Figure 1 shows the DNS results, which exhibit the SSP. To examine the time evolution of the flow fields, the square root of energy of each plane Fourier mode is introduced,

(2.2)\begin{equation} M(n_x,n_z) = \left(\int_{{-}1}^{1} [\hat{{u}}^{2}(n_x,n_z,y) + \hat{{v}}^{2}(n_x,n_z,y) + \hat{{w}}^{2}(n_x,n_z,y)] \,\mathrm{d} y\right)^{{1}/{2}}, \end{equation}

where $(\hat {\cdot })$ denotes the Fourier coefficients for the perturbation velocity, and $n_x$ and $n_z$ define the streamwise and spanwise wavenumbers such that $n_x \alpha$ and $n_z \beta$ ($\alpha = 2{\rm \pi} /L_x$ and $\beta = 2{\rm \pi} /L_z$). Figure 1(a) shows the time trace featuring the quasi-periodic oscillation at the SSP time scale, $T_{SSP} \approx 80\text {--}90$. A sequence of flow field snapshots for a single SSP cycle, which correspond to the time trace in figure 1(b), are shown in figure 1(c). The initial flow field is featured with an amplified state of the high- and low-speed streamwise velocity streaks ($t=1768$) (note that the time evolution of the streaks is depicted by $M(0,1)$ in figure 1b). The streaks subsequently become unstable (Hamilton et al. Reference Hamilton, Kim and Waleffe1995) or experience the related transient growth (Schoppa & Hussain Reference Schoppa and Hussain2002), leading to their breakdown in a sinuously meandering manner. The streak breakdown emerges at $t=1809$ where the streaks evidently disappear. Nonlinear processes subsequently regenerate streamwise vortices, leading to an increase in $M(1,1)$. Finally, for $t=1819\text {--}1840$, the regenerated streamwise vortices ($y$–$z$ cut planes) redistribute the momentum from mean shear, resulting in the formation of new streaks especially near the lower wall region. This mechanism is known as the ‘lift-up’ effect.

Figure 1.

Time trace of $M(0,1)$, $M(0,2)$, $M(1,0)$, $M(1,1)$ (a) for $t=1000\text {--}3000$ and (b) for $t=1768\text {--}1840$. (c) Flow snapshots at $t=1768, 1778, 1789, 1799, 1809, 1819, 1830, 1840$, where the blue and red isosurfaces indicate $u={\pm }0.65$, respectively.

2.2. POD

The POD seeks a set of orthonormal functions that maximises the ensemble-averaged projection of the velocity perturbation, $\boldsymbol {u}$ (e.g. Holmes et al. Reference Holmes, Lumley and Berkooz1996). The optimisation is performed by solving the following eigenvalue problem:

(2.3)\begin{equation} \int_{\varOmega} \left\langle{\boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{u}^{H}(\boldsymbol{x'},t)}\right\rangle \boldsymbol{\varPhi}_{n_x,n_z}^{(n)}(\boldsymbol{x'})\, {\rm d}^{3}\boldsymbol{x'} = \lambda_{n_x,n_z}^{(n)}\boldsymbol{\varPhi}_{n_x,n_z}^{(n)}(\boldsymbol{x}), \end{equation}

where $(\cdot )^{H}$ is the complex conjugate transpose, $\left\langle {\cdot }\right\rangle $ an ensemble average, $\lambda _{n_x,n_z}^{(n)}$, the eigenvalue representing the average kinetic energy contained in each POD mode and $\boldsymbol {\varPhi }_{n_x,n_z}^{(n)}$ are the POD modes. Here, the eigenvalue and the corresponding POD modes are indexed by a positive integer $n$ in decreasing order of the eigenvalue for each pair of the streamwise and spanwise wavenumber indices, $(n_x, n_z)$. The perturbation velocities are then represented in terms of a linear superposition of the POD modes,

(2.4a)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \sum_{n_x={-}\infty}^{\infty}\sum_{n_z={-}\infty}^{\infty}\sum_{n=1}^{\infty} a_{n_x,n_z}^{(n)} (t) \boldsymbol{\varPhi}_{n_x,n_z}^{(n)}(\boldsymbol{x}), \end{equation}

where $a_{n_x,n_z}^{(n)}(t)$ is the time-dependent amplitude of each POD mode. Given the translational invariance in the streamwise and spanwise directions, Fourier expansions are optimal and the POD mode is further written as

(2.4b)\begin{equation} \boldsymbol{\varPhi}_{n_x,n_z}^{(n)}(\boldsymbol{x}) = \frac{1}{\sqrt{L_xL_z}} \exp({{\rm i}(\alpha n_x x+\beta n_z z)}) \boldsymbol{\phi}_{n_x,n_z}^{(n)}(y), \end{equation}

where $\boldsymbol {\phi }_{n_x,n_z}^{(n)}(y)$ describes the wall-normal structure of each POD mode. Since the velocities in physical space are real valued, the following conjugate symmetries are also satisfied:

(2.4c)\begin{equation} a_{n_x,n_z}^{(n)}(t) = a_{{-}n_x,-n_z}^{(n)*}(t) \quad \textrm{and} \quad \boldsymbol{\varPhi}_{n_x,n_z}^{(n)} = \boldsymbol{\varPhi}_{{-}n_x,-n_z}^{(n)*}, \end{equation}

where the superscript $(\cdot )^{*}$ indicates the complex conjugate. Substituting (2.4a) and (2.4b) into (2.3), the eigenvalue problem is simplified to

(2.5)\begin{equation} \int_{{-}1}^{1} \underbrace{\left\langle{\hat{\boldsymbol{u}}(n_x,n_z,y,t)\hat{\boldsymbol{u}}^{H}(n_x,n_z,y',t)}\right\rangle}_{{\equiv} \hat{\boldsymbol{R}}}\boldsymbol{\phi}_{n_x,n_z}^{(n)}(y') \,\mathrm{d} y' = \lambda_{n_x,n_z}^{(n)}\boldsymbol{\phi}_{n_x,n_z}^{(n)}(y). \end{equation}

The ensemble-averaged covariance $\hat {\boldsymbol {R}}$ is computed by enforcing the discrete symmetries of the PCF,

(2.6)\begin{equation} \left.\begin{gathered} \mathcal{I}:[{u},{v},{w},p](x,y,z,t) = [{u},{v},{w} ,p](x,y,z,t),\\ \mathcal{P}:[{u},{v},{w},p](x,y,z,t) = [{-u},{-v},{-w}, p]({-}x,-y,-z,t),\\ \mathcal{R}:[{u},{v},{w},p](x,y,z,t) = [{u},{v},{-w}, p](x,y,-z,t),\\ \mathcal{R} \mathcal{P}:[{u},{v},{w},p](x,y,z,t) = [{-u},{-v},{w}, p]({-}x,-y,z,t), \end{gathered}\right\} \end{equation}

where $\mathcal {I}$ is the identify transformation, $\mathcal {P}$ is a point reflection about the origin, $\mathcal {R}$ is a reflection about the $z$-plane and $\mathcal {RP}$ is a $180^{\circ }$ rotation about the $z$-axis. Given the statistically stationary nature of the turbulent state and the invariance of the PCF under the discrete group transformation in (2.6), the covariance operator is obtained as

(2.7)\begin{equation} \hat{\boldsymbol{R}}(n_x,n_z,y,y')=\frac{1}{4T} \int_{t_0}^{t_0+T} \sum_{\mathcal{T}\in \mathcal{D}_2}(\widehat{\mathcal{T}:{\boldsymbol{u}}})(n_x,n_z,y,t)(\widehat{\mathcal{T}:{\boldsymbol{u}}})^{H}(n_x,n_z,y',t) \,\mathrm{d}t, \end{equation}

where $\mathcal {D}_2 = \{\mathcal {I},\mathcal {P},\mathcal {R},\mathcal {RP}\}$.

For the computation of the POD modes, the ensemble-averaged covariance $\hat {\boldsymbol {R}}$ is first obtained and the eigenvalue problem (2.5) is subsequently solved. Since the ensemble average is equivalent to the time average for a statistically stationary flow, $\hat {\boldsymbol {R}}$ is obtained by averaging in time over an interval $t \in [-10\,000,0]$ with a sampling time interval $\Delta t=1$ (i.e. $t_0=-10\,000$ and $T=10\,000$). Within this time interval, the turbulent state is chosen to be statistically stationary. We also note that the typical time period of the SSP is approximately $T_{SSP} \approx 80\text{--} 90$, implying that more than 100 cycles are considered for the construction of the POD modes. Table 1 shows the leading eigenvalues obtained in the present study, and they are found to match closely with those reported in Smith et al. (Reference Smith, Moehlis and Holmes2005). The structures of the nine leading POD modes are also visualised in figure 2.

Table 1.

Eigenvalues of the first seven POD modes, ranked in terms of the eigenvalue of the POD mode $\lambda$. Here, $E \, [\%]$ is the total energy content of both $(n_x,n_z,n)$ and $(-n_x,-n_z,n)$ wavenumber pairs.

Figure 2.

Visualisation of nine leading POD modes $(n_x,n_z,n)$: (a) $(0,0,1)$; (b) $(0,1,1)$; (c) $(1,0,1)$; (d)  $(1,1,1)$; (e) $(1,-1,1)$; ( f) $(0,2,1)$; (g) $(1,2,1)$; (h) $(1,-2,1)$; (i) $(0,3,1)$. Isosurfaces of (a,b,d–i) denote $u = -0.15$ (red) and $u= 0.15$ (blue). Isosurfaces of (c) denote $w = 0.15$ (teal) and $w=-0.15$ (purple).

3.1. POD-Galerkin projection

To build a reduced-order model, we consider the velocity given by (2.4a) with a finite number of POD modes,

(3.1)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \sum_{n_x={-}M_x}^{M_x}\sum_{n_z={-}M_z}^{M_z}\sum_{n=1}^{N_p} a_{n_x,n_z}^{(n)} (t) \boldsymbol{\varPhi}_{n_x,n_z}^{(n)}(\boldsymbol{x}) + \boldsymbol{u}_R(\boldsymbol{x},t), \end{equation}

where $M_x$, $M_z$ and $N_p$ are the numbers of streamwise Fourier modes, spanwise Fourier modes and $\boldsymbol {\phi }_{n_x,n_z}^{(n)}$, respectively. Here $\boldsymbol {u}_R(\boldsymbol {x},t)$ is the residual velocity field that will not be resolved by the reduced-order model. After substituting (3.1) into (2.1), the projection onto each POD basis yields the following system of ordinary different equations:

(3.2)\begin{equation} \dot{\boldsymbol{a}} = \boldsymbol{L} \boldsymbol{a} +\boldsymbol{N}(\boldsymbol{a},\boldsymbol{a})+\boldsymbol{T}, \end{equation}

where $\dot {(\cdot )}\equiv \textrm {d}/\textrm {d}t$, $\boldsymbol {a} \in \mathbb {C}^{r}$ with $r = [(2M_x + 1)(2M_z+1)N_P+1]/2$ as the column vector, the element of which is given by $a_{n_x,n_z}^{(n)}$ for $-M_x \leqslant n_x \leqslant M_x$, $-M_z \leqslant n_z \leqslant M_z$ and $1 \leqslant n \leqslant N_p$, $\boldsymbol {L}$ and $\boldsymbol {N}$ are from the projection of the linear and quadratic nonlinear parts of (2.1) onto the finite number of POD modes (for further details on their definitions, see Appendix A), and $\boldsymbol {T}$ is the residual term originating from $\boldsymbol {u}_R(\boldsymbol {x},t)$. From (3.2), the simplest low-dimensional model is obtained by ignoring the residual term, i.e.

(3.3)\begin{equation} \boldsymbol{T}=\boldsymbol{0}. \end{equation}

This case shall be referred to as the POD-Galerkin model.

3.2. Sparse POD-Galerkin regression

Several recent studies have been proposed to identify $\boldsymbol {L}$ and $\boldsymbol {N}$ directly from the data (Brunton et al. Reference Brunton, Proctor and Kutz2016; Loiseau & Brunton Reference Loiseau and Brunton2018). If the dimension of the given nonlinear oscillation is sufficiently low (e.g. two-dimensional laminar cylinder wake), $\boldsymbol {T}$ can be ignored with the use of a small number of POD modes. In this case, $\boldsymbol {L}$ and $\boldsymbol {N}$ can directly be obtained from the snapshots of $\boldsymbol {a}$ using a sparse regression technique (i.e. SINDy; e.g. Brunton et al. (Reference Brunton, Proctor and Kutz2016)). Now, SINDy is formulated by first collecting a set of time snapshots of the POD amplitudes from DNS (e.g. $\boldsymbol {a}_{dns}(t)$) into the following data matrices:

(3.4a)\begin{equation} \boldsymbol{X} = \left[\boldsymbol{a}_{dns}(t_1) \quad \boldsymbol{a}_{dns}(t_2) \quad \ldots \quad \boldsymbol{a}_{dns}(t_m) \right]^{{\rm T}} \end{equation}

and

(3.4b)\begin{equation} \dot{\boldsymbol{X}} = \left[ \dot{\boldsymbol{a}}_{dns}(t_1) \quad \dot{\boldsymbol{a}}_{dns}(t_2) \quad \ldots \quad \dot{\boldsymbol{a}}_{dns}(t_m) \right]^{{\rm T}}. \end{equation}

A set of candidate library functions, $\varTheta (\boldsymbol {X})$, is subsequently constructed. In the present study, we restrict the library functions to being

(3.5a)\begin{equation} \varTheta(\boldsymbol{X}) = \left[\mathbb{P}_{\boldsymbol{L}}(\boldsymbol{X})\quad \mathbb{P}_{\boldsymbol{N}}(\boldsymbol{X})\right] =\begin{bmatrix} \mathbb{P}_{\boldsymbol{L}}(\boldsymbol{a}_{dns}(t_1)) & \mathbb{P}_{\boldsymbol{N}}(\boldsymbol{a}_{dns}(t_1))\\ {}\mathbb{P}_{\boldsymbol{L}}(\boldsymbol{a}_{dns}(t_2)) & \mathbb{P}_{\boldsymbol{N}}(\boldsymbol{a}_{dns}(t_2))\\ \vdots & \vdots\\ \mathbb{P}_{\boldsymbol{L}}(\boldsymbol{a}_{dns}(t_m)) & \mathbb{P}_{\boldsymbol{N}}(\boldsymbol{a}_{dns}(t_m)) \end{bmatrix}, \end{equation}

where $\mathbb {P}_{\boldsymbol {L}}(\boldsymbol {X})$ and $\mathbb {P}_{\boldsymbol {N}}(\boldsymbol {X})$ denote the linear and quadratic combinations of the state vector admitted by the form of $\boldsymbol {L}$ and $\boldsymbol {N}$ in (3.2), respectively. We also introduce a coefficient matrix

(3.5b)\begin{equation} \boldsymbol{\varXi} = \left[{\boldsymbol{\varXi}}_{\boldsymbol{L}} \quad {\boldsymbol{\varXi}}_{\boldsymbol{N}} \right]^{\rm T}, \end{equation}

where ${\boldsymbol {\varXi }}_{\boldsymbol {L}}$ and ${\boldsymbol {\varXi }}_{\boldsymbol {N}}$ contain the (unknown) coefficients for $\mathbb {P}_{\boldsymbol {L}}(\boldsymbol {X})$ and $\mathbb {P}_{\boldsymbol {N}}(\boldsymbol {X})$. The coefficient matrix, $\boldsymbol {\varXi }$, is subsequently determined by solving the following convex least squares regression problem:

(3.6)\begin{equation} \min_{\boldsymbol{\varXi}} \| \dot{\boldsymbol{X}} - \varTheta(\boldsymbol{X}) \, \boldsymbol{\varXi} \|_2 + \gamma \| \boldsymbol{\varXi} \|_1, \end{equation}

where $\| \cdot \|_2$ and $\| \cdot \|_1$ denote the standard ${\ell }_2$ and ${\ell }_1$ norms, respectively, and $\gamma$ is the penalty introduced for the sparsity promoting ${\ell }_1$-regulariser. An advantage of using an ${\ell }_1$-regularisation as opposed to an ${\ell }_2$-regularisation is that they tend to prevent data overfit by promoting a model with the least complexity (sparse) required to model the dynamics (Brunton & Kutz Reference Brunton and Kutz2019). The optimisation problem in (3.6) can be solved with the well known LASSO (least absolute shrinkage and selection operator) algorithm. For large datasets, an alternative based on sequential threshold least squares was recommended instead – i.e. the SINDy approach (Brunton et al. Reference Brunton, Proctor and Kutz2016). This approach is used in the present study. We note that if the given dataset is even larger, the approach proposed by Gelß et al. (Reference Gelß, Klus, Eisert and Schütte2019) may further be considered (i.e. multidimensional approximation of nonlinear dynamical systems (known as MANDy)). The sparse regression using a template given by the Galerkin projection will be referred to as the POD-SINDy model.

It is worth mentioning that the optimisation (3.6) was recently proposed to be solved with an equality constraint which explicitly enforces the energy conservation in the nonlinear operator $\boldsymbol {N}$ (Loiseau & Brunton Reference Loiseau and Brunton2018): i.e. $\boldsymbol {a}^{H}\boldsymbol {N(a,a)}=0$. However, this approach will not be considered in the present study, where the number of POD modes for the construction of a reduced-order model will not necessarily be large enough to fully cover the dimension of the given chaotic state. In such a case, the sparse regression in (3.6) implies that the residual term $\boldsymbol {T}$ in (3.2) is modelled as

(3.7)\begin{equation} \boldsymbol{T}=\boldsymbol{L}_{\boldsymbol{\varDelta}}\boldsymbol{a}+\boldsymbol{N}_{\boldsymbol{\varDelta}}(\boldsymbol{a},\boldsymbol{a}), \end{equation}

where $\boldsymbol {L}_\varDelta$ and $\boldsymbol {N}_{\boldsymbol {\varDelta }}$ are the differences in $\boldsymbol {L}$ and $\boldsymbol {N}$ obtained by the Galerkin projection and the sparse regression, respectively. Therefore, the energy-conservation constraint for the nonlinear operator proposed by Loiseau & Brunton (Reference Loiseau and Brunton2018) would enforce $\boldsymbol {a}^{H}(\boldsymbol {N(a,a)}+\boldsymbol {N}_{\boldsymbol {\varDelta }}(\boldsymbol {a},\boldsymbol {a}))=0$. However, this is not necessarily desirable in the present case, as will be discussed in the following subsection. Hence in § 4.2, the library terms for the set of POD modes from the same wavenumber space, $\boldsymbol {a}_{n_x,n_z}$, are made dependent on the wavenumber of the POD modes as is suggestive from a Galerkin projection. The library function will consist of first-order polynomials of POD modes from the same wavenumber space. It attempts to model the $\boldsymbol {L_{\varDelta }}$ terms of the dynamics. The modelling of the $\boldsymbol {N_{\varDelta }}$ terms are accounted for by including the combination of POD modes that form a set of second-order polynomials which reflect the triadic wavenumber interactions, where $(n_x,n_z) = (k_x,k_z) + (m_x,m_z)$, from a POD-Galerkin model.

3.3. The needs for an eddy-viscosity closure

Now, we consider a physical model for $\boldsymbol {T}$. To rationalise the need of such a model, we first define $\boldsymbol {a}_{\infty }$ as the solution to (3.2) for infinitely large $M_x$, $M_z$ and $N_p$. In this case, the residual $\boldsymbol {T}$ in (3.2) should vanish, thus $\boldsymbol {a}_{\infty }$ would be identical to those obtained by projecting $\boldsymbol {u}$ from DNS onto the POD modes. Given the energy-conserving nature of the nonlinear term in (2.1) (e.g. Joseph Reference Joseph1976), the contribution of the resulting nonlinear term to the change rate of the perturbation kinetic energy ($\boldsymbol {a}_\infty ^{H} \boldsymbol {a}_\infty$) should be zero for every time instance, $t$, i.e.

(3.8)\begin{equation} \boldsymbol{a}_{\infty}^{H}\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty},\boldsymbol{a}_{\infty}) = 0, \end{equation}

where $\boldsymbol {N}_\infty$ is the quadratic nonlinear term obtained by considering infinitely large $M_x$, $M_z$ and $N_p$. We note that (3.8) must also be true even if $\boldsymbol {a}_{\infty }$ is replaced by any arbitrary vector. This observation motivated Loiseau & Brunton (Reference Loiseau and Brunton2018) to impose $\boldsymbol {a}^{H}\boldsymbol {N}(\boldsymbol {a},\boldsymbol {a})=0$ as an equality constraint into the optimisation problem in (3.6).

Let us now consider small values of $M_x$, $M_z$ and $N_p$ which define the size of the reduced-order model. In particular, we will assume that $M_x$, $M_z$ and $N_p$ are not large enough to cover the full energy cascade dynamics of the given turbulent state. We define a projection operator $\mathcal {P}_l$ onto the subspace defined by the small values of $M_x$, $M_z$ and $N_p$. Then, $\boldsymbol {a}_{\infty }$ can be decomposed into $\boldsymbol {a}_{\infty }=\boldsymbol {a}_{\infty,l}+\boldsymbol {a}_{\infty,h}$, where $\boldsymbol {a}_{\infty,l}=\mathcal {P}_l[\boldsymbol {a}_\infty ]$ and $\boldsymbol {a}_{\infty,h}=\mathcal {P}_h[\boldsymbol {a}_\infty ]$ with $\mathcal {I}_0 [\cdot ]=\mathcal {P}_l[\cdot ]+\mathcal {P}_h[\cdot ]$ (where $\mathcal {I}_0 [\cdot ]$ is the identity operator). Using this decomposition and the quadratic nature of $\boldsymbol {N}_\infty$, (3.8) can be written as

(3.9)\begin{align} &\boldsymbol{a}_{\infty,l}^{H}[\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,l},\boldsymbol{a}_{\infty,h})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,l})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,h})] \nonumber\\ &\quad +\boldsymbol{a}_{\infty,h}^{H}[\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,l},\boldsymbol{a}_{\infty,l})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,l},\boldsymbol{a}_{\infty,h})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,l})]=0, \end{align}

where the top and bottom lines describe the nonlinear energy transport of the perturbation kinetic energy in the $\mathcal {P}_l$ and $\mathcal {P}_h$ subspaces, respectively. Here, we note that the term $\boldsymbol {a}_{\infty,l}^{H}\boldsymbol {N}_\infty (\boldsymbol {a}_{\infty,l},\boldsymbol {a}_{\infty,l})$, equivalent to $\boldsymbol {a}^{H}\boldsymbol {N}(\boldsymbol {a},\boldsymbol {a})$ expected from (3.2), vanishes due to (3.8), and so does $\boldsymbol {a}_{\infty,h}^{H}\boldsymbol {N}_\infty (\boldsymbol {a}_{\infty,h},\boldsymbol {a}_{\infty,h})$ (note that the only difference between $\boldsymbol {a}_{\infty,l}^{H}\boldsymbol {N}_\infty (\boldsymbol {a}_{\infty,l},\boldsymbol {a}_{\infty,l})$ and $\boldsymbol {a}^{H}\boldsymbol {N}(\boldsymbol {a},\boldsymbol {a})$ is their dimension). The top line in (3.9) should indicate the rate of perturbation energy transferred from the $\mathcal {P}_l$ to the $\mathcal {P}_h$ subspace. Importantly, the energy cascade from large to small scales in the three-dimensional Navier–Stokes equations implies

(3.10a)\begin{equation} \langle \boldsymbol{a}_{\infty,l}^{H}[\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,l},\boldsymbol{a}_{\infty,h})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,l})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,h})]\rangle<0 \end{equation}

and

(3.10b)\begin{equation} \langle \boldsymbol{a}_{\infty,h}^{H}[\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,h})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,l})+\boldsymbol{N}_\infty(\boldsymbol{a}_{\infty,h},\boldsymbol{a}_{\infty,h})]\rangle>0 \end{equation}

for a statistically stationary flow where the ensemble average is equivalent to a time average. Indeed, this has been observed in a number of previous studies where the interscale energy transfer is analysed in detail (e.g. Cho et al. Reference Cho, Hwang and Choi2018; Lee & Moser Reference Lee and Moser2019; Hwang & Lee Reference Hwang and Lee2020; Doohan et al. Reference Doohan, Willis and Hwang2021). Given the equivalence of $\boldsymbol {a}_{\infty,l}^{H}\boldsymbol {N}_\infty (\boldsymbol {a}_{\infty,l},\boldsymbol {a}_{\infty,l})$ to $\boldsymbol {a}^{H}\boldsymbol {N}(\boldsymbol {a},\boldsymbol {a})$, a condition for $\boldsymbol {T}$ to meet from a physical viewpoint would then be

(3.11)\begin{equation} \langle {\boldsymbol{a}^{H}\boldsymbol{T}} \rangle<0, \end{equation}

indicating that the residual term, $\boldsymbol {T}$, in (3.2) must contain an energy removal mechanism.

The discussion above evidently justifies the use of an eddy viscosity model in many previous studies (e.g. Aubry et al. Reference Aubry, Holmes, Lumley and Stone1988; Rempfer & Fasel Reference Rempfer and Fasel1994a,Reference Rempfer and Faselb; Smith et al. Reference Smith, Moehlis and Holmes2005; Östh et al. Reference Östh, Noack, Krajnović, Barros and Borée2014; Protas et al. Reference Protas, Noack and Östh2015), i.e.

(3.12)\begin{equation} \boldsymbol{T}={\nu}_t\boldsymbol{Da}, \end{equation}

where $\nu _t$ is a scalar-valued eddy viscosity and $\boldsymbol {D}$ is the Laplacian operator for the reduced-order model defined in Appendix A. For $\nu _t>0$, (3.12) yields a sufficient condition for (3.11), as it ensures $\boldsymbol {a}^{H}\boldsymbol {T}(=\nu _t \boldsymbol {a}^{H}\boldsymbol {D})\boldsymbol {a<0}$ at every time for non-zero $\boldsymbol {a}$ with $(n_x,n_z) \neq (0,0)$ due to the negative semidefinite nature of $\boldsymbol {D}$ (see Appendix A). In the previous studies, various form of eddy viscosity have been introduced and examined. In the present study, we will first consider a simple spectral eddy-viscosity model similar to the one in Smith et al. (Reference Smith, Moehlis and Holmes2005) where an empirical real-valued constant of $\nu _t$ is employed. It is evident that the performance of this simple single-valued eddy-viscosity model is expected to be limited, as there is evidence supporting the complex nature of the intermodal energy transfer (e.g. Couplet, Sagaut & Basdevant Reference Couplet, Sagaut and Basdevant2003; Podvin Reference Podvin2009, see also § 4.3). Therefore, the approach here is considered for the purpose of comparing with the other models. This approach shall be referred to as the POD-Galerkin-E model (the ‘E’ stands for an ‘empirical’ eddy viscosity).

3.4. Sparse optimal eddy-viscosity closure

Although the eddy-viscosity model in (3.12) ensures the physical property that originates from the energy cascade of turbulent state, (3.11), it is evidently too crude. Indeed, $\boldsymbol {a}^{H}\boldsymbol {T}$ does not have to be negative for every point in time like the one ensured by (3.12) – only its time average needs to be negative. Furthermore, there is no physical reason why different POD modes should feel the ‘same’ eddy viscosity; for example, the POD modes for large $n_x$, $n_z$ and $n$ are not expected to experience a large amount of energy removal from their dynamics by $\boldsymbol {T}$. More flexible forms of eddy viscosity have therefore been proposed previously (for a review, see Östh et al. (Reference Östh, Noack, Krajnović, Barros and Borée2014), where various forms of eddy viscosity have been examined for a three-dimensional turbulent bluff-body wake). In particular, Protas et al. (Reference Protas, Noack and Östh2015) introduced the concept of ‘optimal’ eddy viscosity by formulating an adjoint-based optimisation problem that minimises the difference between the POD amplitudes from the measurement and the reduced-order model. However, as discussed in § 1, the application of the adjoint-based optimisation to a turbulent flow does not always allow for a sufficiently long optimisation time interval due to the unstable nature of the adjoint system around a chaotic state. Also, the approach of Protas et al. (Reference Protas, Noack and Östh2015) is still based on a scalar-valued eddy viscosity, although its generalisation to a sophisticated form of eddy viscosity is easily possible.

In the present study, we take an alternative formulation which enables us to consider a long time horizon for a similar optimisation problem. In particular, we determine an optimal eddy viscosity with the sparse regression in § 3.2. For the purpose of demonstration, we consider a nonlinear closure model for $\boldsymbol {T}$ (Östh et al. Reference Östh, Noack, Krajnović, Barros and Borée2014; Protas et al. Reference Protas, Noack and Östh2015), in which a complex-valued eddy viscosity is set to vary with POD modes,

(3.13a)\begin{equation} \boldsymbol{T}=\boldsymbol{V}_t\boldsymbol{Da} \end{equation}

with

(3.13b)\begin{equation} \boldsymbol{V}_t(t)=e(t) \,{\rm diag}[\boldsymbol{c}], \end{equation}

where $e(t)$ is a scalar-valued function of $\boldsymbol {a}$ that can be chosen for a nonlinear eddy-viscosity model and $\boldsymbol {c}$ the constant vector to be determined for each $n_x$, $n_z$ and $n$. We note that the choice of $e$ usually ensures $\boldsymbol {V}_t$ will vanish as $\boldsymbol {a}$ becomes zero. In previous studies (e.g. Östh et al. Reference Östh, Noack, Krajnović, Barros and Borée2014; Protas et al. Reference Protas, Noack and Östh2015), the perturbation kinetic energy, $e(t)=(\boldsymbol {a}^{H}\boldsymbol {a})^{1/2}$, has often been considered. In the present study, a simpler form, $e(t)=|a_{0,0}^{(1)}|$, is chosen, given that the corresponding mode contains approximately $70\,\%$ of the total perturbation kinetic energy (see table 1). A preliminary test also reveals that this choice makes the optimised reduced-order model perform slightly better than that of $e(t)=(\boldsymbol {a}^{H}\boldsymbol {a})^{1/2}$ (see Appendix B).

Now, we formulate an optimisation problem determining $\boldsymbol {V}_t$ in (3.13b). Using the POD mode amplitudes obtained from DNS (i.e. $\boldsymbol {a}_{dns}$) and (3.2), the desired residual term of the low-dimensional system is given by

(3.14)\begin{equation} \boldsymbol{T}_{dns}(t) \equiv \dot{\boldsymbol{a}}_{dns}(t)-\boldsymbol{La}_{dns}-\boldsymbol{N}(\boldsymbol{a}_{dns},\boldsymbol{a}_{dns}), \end{equation}

where $\boldsymbol {T}_{dns}$ indicates the residual term calculated with $\boldsymbol {a}_{dns}$. Similarly to the sparse regression in § 3.2, a set of time snapshots of $\boldsymbol {T}_{dns}(t)$ is introduced into a data matrix using (3.14):

(3.15a)\begin{equation} \boldsymbol{Y} = \left[ \boldsymbol{T}_{dns}(t_1) \quad \boldsymbol{T}_{dns}(t_2) \quad \ldots \quad \boldsymbol{T}_{dns}(t_m) \right]^{{\rm T}}. \end{equation}

Given (3.13a), the related library function for the optimisation is subsequently formed to be

(3.15b)\begin{equation} \varTheta_e(\boldsymbol{Y}) = \left[ e_{dns}(t_1)\boldsymbol{a}_{dns}(t_1)\boldsymbol{D} \quad e_{dns}(t_2)\boldsymbol{a}_{dns}(t_2)\boldsymbol{D} \quad \ldots \quad e_{dns}(t_m)\boldsymbol{a}_{dns}(t_m)\boldsymbol{D} \right]^{{\rm T}}, \end{equation}

where $e_{dns}(t)$ is $e$ obtained from DNS. Since we seek a complex-valued constant vector $\boldsymbol {c}$ that minimises the difference between $\boldsymbol {T}_{dns}$ given by (3.14) and $\boldsymbol {T}$ from the residual model in (3.13a), the following optimisation is defined:

(3.16)\begin{equation} \min_{\boldsymbol{c}} \| \boldsymbol{\dot{Y}} - \varTheta_e(\boldsymbol{Y}) \, {{\rm diag}[\boldsymbol{c}]} \|_2 + \gamma_e \| \boldsymbol{c} \|_1, \end{equation}

where $\gamma _e$ is the parameter for the sparsity-promoting ${\ell }_1$-regulariser. Like the optimisation problem in (3.6), (3.16) is solved by applying the SINDy approach (Brunton et al. Reference Brunton, Proctor and Kutz2016). This approach will be referred to as POD-Galerkin-R model (the ‘R’ stands the determination of an eddy viscosity with a sparse ‘regression’). As mentioned in § 1, several previous studies (e.g. Couplet et al. Reference Couplet, Basdevant and Sagaut2005; Cordier et al. Reference Cordier, El Majd and Favier2010) proposed a similar idea of calibrating the residual term, $\boldsymbol {T}$, using Tikhonov regularisation, an optimisation based on the ${\ell _2}$-regularisation. Here, the ${\ell }_1$-regularisation has a benefit over the ${\ell }_2$-regularisation, as it is designed to prevent data overfit, offering a parsimonious low-dimensional model. Lastly, we note that setting $\gamma _e \rightarrow \infty$ results in the POD-Galerkin model, while $\gamma _e \rightarrow 0$ yields the least-square eddy viscosity (for a further discussion on the effect of $\gamma _e$, see Appendix C).

It is also useful to make some remarks on the least-square sparse regression in (3.16). First, the regression (3.16) can now consider a very long sampling time, as it simply relies on the POD mode amplitudes, $\boldsymbol {a}_{dns}$, taken from DNS. Therefore, it no longer suffers from the finite optimisation time-interval issue that one may face in the conventional adjoint-based formulation. Second, the regression problem (3.16) can flexibly be formulated by accounting for various forms of eddy viscosity, and this can be achieved by adding more library functions in (3.15b). Third, the regression (3.16) allows for negative elements of $\boldsymbol {V}_t$ (or $\boldsymbol {c}$), indicating that the ‘backwards scattering’ in the energy cascade can be taken into account. Although this issue may cause a potential numerical instability of the resulting reduced-order model, it may well be fixed by imposing an additional inequality constraint (e.g. all the element of $\boldsymbol {c}$ is greater than or equal to zero). Finally, $\boldsymbol {V}_t$ (or $\boldsymbol {c}$) in (3.13b) can also be a complex vector, given that $\boldsymbol {a}$ is complex. In other words, some ‘dispersive’ effect in the dynamics of the reduced-order model can also be added with non-zero imaginary part of $\boldsymbol {V}_t$.