1. Introduction
Modelling a turbulent flow often involves a large number of degrees of freedom, especially at high Reynolds numbers. One of the most fundamental challenges in this modelling process is to have an accurate description of the spatio-temporal evolution of flow fields while reducing the degrees of freedom (DoFs) as much as possible. An approach that has been pursued is the development of reduced-order models (ROMs) based on the concepts of dynamical systems. In this approach, the development of a finite-dimensional dynamical system with much reduced DoFs is aimed at (see review by Rowley & Dawson Reference Rowley and Dawson2017). It has seen success in flow modelling especially for identification of the dominant coherent structures and related physical processes (e.g. Sirovich Reference Sirovich1987; Aubry et al. Reference Aubry, Holmes, Lumley and Stone1988; Holmes, Lumley & Berkooz Reference Holmes, Lumley and Berkooz1996; Waleffe Reference Waleffe1997; Moehlis, Faisst & Eckhardt Reference Moehlis, Faisst and Eckhardt2004; Cavalieri Reference Cavalieri2021; Cavalieri, Rempel & Nogueira Reference Cavalieri, Rempel and Nogueira2022; McCormack, Cavalieri & Hwang Reference Maia and Cavalieri2024) and flow control applications in transitional flows (e.g Rowley, Colonius & Murray Reference Rowley, Colonius and Murray2004; Rowley Reference Rowley2005; Ilak & Rowley Reference Ilak and Rowley2008; Dergham et al. 2011; Rowley & Dawson Reference Rowley and Dawson2017; Maia & Cavalieri Reference Maia and Cavalieri2025, and many others).
Perhaps, one of the most well-studied approaches in reduced-order modelling is based on a Galerkin projection of the equations of motion using the proper orthogonal decomposition (POD) modes (Sirovich Reference Sirovich1987; Aubry et al. Reference Aubry, Holmes, Lumley and Stone1988; Holmes et al. Reference Holmes, Lumley and Berkooz1996; Moehlis et al. Reference Moehlis, Smith, Holmes and Faisst2002; Gibson Reference Gibson2002; Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003; Smith, Moehlis & Holmes Reference Smith, Moehlis and Holmes2005). This method first obtains a set of orthogonal basis functions that best represent the given flow in terms of its energy content (POD modes) and subsequently projects the equations of motion onto the space spanned by these basis functions. With suitable truncation in the number of POD modes, a finite-dimensional representation of the original fluid flow is obtained in the form of a dynamical system. An important challenge in this type of approach is often associated with truncation in the number of POD modes. The exclusion of higher-order POD modes, which often represent small-scale motions with a lower energy content, causes an issue in modelling dynamically important, but not energetic, flow structures involved in complex nonlinear processes, such as turbulent energy cascade and dissipation. This greatly affects the accuracy of the ROMs and sometimes causes numerical instabilities. To address this limitation, additional modelling efforts have been made, such as incorporation of eddy viscosity, to improve the accuracy and stability of ROMs (see e.g. Rempfer & Fasel Reference Rempfer and Fasel1994a , Reference Reddy, Schmid, Baggett and Henningsonb ; Moehlis et al. Reference Moehlis, Smith, Holmes and Faisst2002; Smith et al. Reference Smith, Moehlis and Holmes2005; Östh et al. Reference Östh, Noack, Krajnović, Barros and Borée2014; Chua Khoo et al. 2022).
POD modes are, however, not the only possible basis functions that can be used for the design of ROMs. In principle, any set of orthogonal function satisfying the boundary conditions may be used to approximate the projection of the Navier–Stokes equations into a finite dimensional dynamical system due to its fast (often exponential) converging nature, and this is the theoretical basis underpinning the spectral discretisation (e.g. Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988). Indeed, some ROMs that aim to model the key underlying physical processes are based on simple orthogonal functions, such as trigonometric functions (e.g. Waleffe Reference Waleffe1997; Moehlis et al. Reference Moehlis, Faisst and Eckhardt2004; Cavalieri Reference Cavalieri2021). Furthermore, some recent studies have used eigenfunctions of the controllability Gramian operator of the linearised Navier–Stokes equations: i.e. controllability modes (Cavalieri & Nogueira Reference Cavalieri and Nogueira2022; McCormack et al. Reference Maia and Cavalieri2024; Maia & Cavalieri Reference Maia and Cavalieri2024). The controllability Gramian operator is the covariance of the stochastic response of a given linear dynamical system. Therefore, the resulting modes are the POD modes of the flow fields from the linearised Navier–Stokes equations subject to a stochastic forcing. If the equations are used to model a turbulent flow, the corresponding controllability modes may well be a good candidate for the basis functions to design a ROM; such operator-based modes do not require large data sets and may thus lead to ROMs directly from governing equation.
POD modes and the controllability modes are the basis functions that represent the most energetic structures in the given flow fields, but they do not necessarily capture all dynamically important structures effectively. This is particularly true when energy-containing flow structures are initiated by some non-energetic structures with very different shapes. Balanced truncation is a model reduction method that accounts for this issue from a viewpoint of control theory for linear dynamical systems (Moore Reference Moore1981; Zhou, Doyle & Glover Reference Zhou, Doyle and Glover1996). Furthermore, an approach of approximating the associated modes (i.e. balanced truncation modes) has been proposed to deal with the high state-space dimension of the linearised Navier–Stokes equations (Rowley Reference Rowley2005): balanced POD modes. Application of this type of approach has been shown to be particularly useful for modelling the flow dynamics associated with non-normality of the linearised Navier–Stokes equations (Ilak & Rowley Reference Ilak and Rowley2008; Bagheri et al. Reference Bagheri, Henningson, Hœpffner and Schmid2009).
The scope of the present study is to systematically develop a ROM based on Galerkin projections for turbulent flows, with the ultimate goal of reaching high Reynolds numbers. Recently, there have been an increasing number of studies to develop ROMs for turbulent flows using data-driven approaches (e.g. Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016; Kaszás et al. Reference Kaszás, Cenedese and Haller2022; Linot & Graham Reference Linot and Graham2023, and many others), but here we shall focus on the construction of ROMs based on Galerkin projections. Although this approach is classical and its capability in reducing the DoFs may not be as good as that of the data-driven counterparts, it offers some important advantages from a physical viewpoint. First, it preserves the form of the original Navier–Stokes equations, allowing one to interpret and assess its complex chaotic dynamics in terms of the turbulent energy budget (i.e. turbulence production, inter-scale and space turbulent transport, dissipation, etc.; see McCormack et al. Reference Maia and Cavalieri2024, for example). Second, the Galerkin projection is robust compared with data-driven approaches, as it neither has any overfitting issues nor relies on the geometrical structure of state space. This is important, especially at high Reynolds numbers, where the number of degrees of freedom of a suitable ROM is generally expected to be large and the geometry of state space would be too complicated to obtain useful physical insights. Finally, when properly implemented, the Galerkin projection ensures that the ROM converges to the corresponding DNS, providing a way to properly assess its performance.
Of particular interest in the present study is understanding how the form of the basis functions affects the convergence and performance (especially with severe truncation) of a ROM around both laminar and turbulent state, in the hope of finding a suitable set of basis functions that are efficient and computationally economical to generate for a given modelling purpose. To this end, we will consider various basis functions discussed previously: (i) POD; (ii) controllability; and (iii) balanced truncation modes. Furthermore, the controllability and balanced modes will be generated by considering a few different forms of the linearised Navier–Stokes equations. We note that the POD modes are generated by the flow field, i.e. the modelling target itself. Therefore, numerical simulations or experiments need to be carried out before reduced-order modelling and they are often expensive. In contrast, the controllability and balanced truncation modes are generated by computing controllability and observability Gramian operators of the linearised Navier–Stokes equations. In parallel shear flows, the computational cost of these Gramian operators is modest because Fourier decomposition in the wall-parallel directions allows their size to be relatively small (e.g. Farrell & Ioannou Reference Farrell and Ioannou1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; Hwang & Cossu Reference Hwang and Cossu2010a , Reference Hwang and Cossub ). Even in complex flows, there are a few suitable computational approaches to approximating them (e.g. Rowley Reference Rowley2005; Ilak & Rowley Reference Ilak and Rowley2008). However, in this case, each Gramian operator may well be formulated to be a function of the flow about which the linearised Navier–Stokes equations are obtained. For modelling of transition, the laminar base flow may be the obvious choice (Jovanović & Bamieh Reference Jovanović and Bamieh2005; Rowley Reference Rowley2005; Ilak & Rowley Reference Ilak and Rowley2008), but this is not necessarily the case for turbulent flows. Although the laminar base flow can still be used for the generation of Gramian operators (Cavalieri et al. Reference Cavalieri, Rempel and Nogueira2022; McCormack et al. Reference Maia and Cavalieri2024), turbulent mean flow, which may be approximated by considering an RANS model in practice, can also be used (Hwang & Cossu Reference Hwang and Cossu2010a ; McKeon & Sharma Reference McKeon and Sharma2010). More recently, there has been growing evidence that the linearised Navier–Stokes equations enhanced by a suitable eddy viscosity (Reynolds & Hussain Reference Reynolds and Hussain1972) provide a better description of coherent structures in turbulent flows (e.g. Hwang Reference Hwang2016; Illingworth, Monty & Marusic Reference Illingworth, Monty and Marusic2018; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021; Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023).
The emphasis of the present study will therefore be on understanding the effectiveness of different basis functions in approximating each state space (around laminar or turbulent state) using a Galerkin-based ROM. In addition to POD modes obtained from turbulent flow, several sets of controllability and balanced truncation modes with various forms of the linearised Navier–Stokes operator will be considered: for example, variations in base flow profile (laminar base flow or turbulent mean flow) and presence of the eddy viscosity model. To identify the strength and weakness of each of the basis functions in different parts of the state space, we evaluate the linear stability and transient growth of the laminar state while examining statistics (mean and fluctuations) and dynamics (time correlations) of the turbulent state. The benchmark case is chosen as plane Couette flow at a relatively low Reynolds number within a minimal flow unit (Jiménez & Moin Reference Jiménez and Moin1991; Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995), a canonical parallel shear flow where the dynamics and coherent structures are well understood.
This paper is organised as follows. In § 2, we provide a general, but concise introduction to the basis functions to be examined: POD, controllability and balanced truncation modes. POD modes and the linear operators for the generation of controllability and balanced truncation modes will then be introduced in § 3. In § 4, the performance and convergence of each ROM will be presented around both laminar and turbulent state, followed by detailed discussions in § 5. The paper will conclude in § 6, with some perspectives.
2. Basis functions
We start by introducing the different sets of basis functions to be used in the present study. Throughout this section, these basis functions are denoted in matrix form,
$\varPhi (\boldsymbol{x}) =[ \phi _1(\boldsymbol{x})\,\phi _2(\boldsymbol{x})\,\ldots ]$
, so that each column represents the individual basis function. If the basis functions for incompressible flows are sufficiently smooth, satisfy homogenous or periodic boundary condition and are divergence-free, a velocity fluctuation field
$\boldsymbol{u}(\boldsymbol{x},t)$
(
$\boldsymbol{x}$
and
$t$
are space and time, respectively) may be written in terms of their linear combination:
\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \sum _{i=1}^\infty a_i(t){\phi }_i(\boldsymbol{x}), \end{equation}
where
$\phi _i$
is the
$i$
th basis function. The reconstructed velocity field in (2.1) satisfies the continuity equation and the boundary conditions by construction. If the basis functions
$\phi _i(\boldsymbol{x})$
are orthogonal, the series expansion in (2.1) would often converge quickly for sufficiently smooth
$\boldsymbol{u}(\boldsymbol{x},t)$
. If not (e.g. balanced modes), the convergence is not necessarily guaranteed in general, although such a behaviour has not been observed in the present study.
2.1. POD modes
The POD seeks the basis functions
$\phi$
by formulating an optimisation problem that maximises the projection of
$\boldsymbol{u}$
onto
$\phi$
, such that
\begin{equation} \max _\phi \frac { \overline {| \langle \boldsymbol{u}, \phi \rangle |^2} }{ \langle \phi ,\phi \rangle }, \end{equation}
where
$\overline {(\boldsymbol{\cdot })}$
denotes ensemble average and
$\langle \boldsymbol{\cdot }, \boldsymbol{\cdot }\rangle$
represents the inner product defined by
$\langle \boldsymbol{f},\boldsymbol{g}\rangle = \int _\varOmega \boldsymbol{g}^H\boldsymbol{f}\,\textrm {d}\boldsymbol{x}$
in the flow domain
$\varOmega$
, with the superscript
$(\boldsymbol{\cdot })^H$
indicating complex conjugate transpose. The first-order necessary condition for the solution to (2.2) yields an eigenvalue problem that takes the form of
\begin{equation} \int _\varOmega \overline {\boldsymbol{u}(\boldsymbol{x})\boldsymbol{u}^H(\boldsymbol{x}')} \varPhi (\boldsymbol{x}')\, \textrm {d}\boldsymbol{x}' = \varPhi (\boldsymbol{x}) \varLambda , \end{equation}
where
\begin{equation} \varLambda = \textrm {diag}[\lambda _1, \lambda _2,\ldots ]. \end{equation}
Here, the eigenvector
$\phi _i$
in each column of
$\varPhi (\boldsymbol{x})$
forms the individual POD mode, and the eigenvalue
$\lambda _i$
, which is positive and real, indicates the content of the averaged kinetic energy in each POD mode. The POD modes are also orthogonal: i.e.
$\langle \phi _i, \phi _j \rangle = \delta _{ij}$
, where
$\delta _{ij}$
is the Kronecker delta and
$\phi _i$
is normalised accordingly. Consequently, the POD modes form an optimal set of basis functions for reconstruction of the flow field in terms of kinetic energy.
2.2. Controllability modes
We now assume that the velocity fluctuation field
$\boldsymbol{u}$
is obtained from the following linear dynamical system driven by a stochastic body forcing:
\begin{align} \frac {\partial {\boldsymbol{q}}}{\partial t} &= \boldsymbol{A}\boldsymbol{q} + \boldsymbol{B}\boldsymbol{f}, \end{align}
\begin{align} \boldsymbol{u} &= \boldsymbol{C}\boldsymbol{q}, \end{align}
where
$\boldsymbol{q}$
,
$\boldsymbol{A}$
,
$\boldsymbol{B}$
and
$\boldsymbol{C}$
are defined later in § 3.2 for plane Couette flow. For linear stable
$\boldsymbol{A}$
and zero initial condition,
$\boldsymbol{u}$
is given as
\begin{equation} \boldsymbol{u}(t)=\int _0^t \boldsymbol{C}e^{(t-s)\boldsymbol{A}}\boldsymbol{B}\boldsymbol{f}(s)\,{\rm d}s \end{equation}
for
$t\rightarrow \infty$
. Using (2.5), the covariance of the velocity fluctuation can subsequently be formulated. If the body forcing is white in time such that
$\overline {\boldsymbol{f}(t_1) \otimes \boldsymbol{f}(t_2)}=\boldsymbol{I}\delta (t_1-t_2)$
, it can be shown that the covariance operator of the velocity fluctuation is computed from
\begin{equation} \boldsymbol{W}_c\boldsymbol{v}=\int _\varOmega \overline {\boldsymbol{u}(\boldsymbol{x})\boldsymbol{u}^H(\boldsymbol{x}')} \boldsymbol{v}(\boldsymbol{x}')\, \textrm {d}\boldsymbol{x}', \end{equation}
where
$v(\boldsymbol{x})$
is a test function sharing the same spatial dimension with
$\boldsymbol{u}$
, and
\begin{equation} \boldsymbol{W}_c=\boldsymbol{C}\boldsymbol{W}_c^q\boldsymbol{C}^{\dagger } \end{equation}
with
\begin{equation} \boldsymbol{A}\boldsymbol{W}_c^q + \boldsymbol{W}_c^q\boldsymbol{A}^{\dagger } + \boldsymbol{BB}^\dagger = 0. \\ \end{equation}
Here, the superscript
$(\boldsymbol{\cdot })^\dagger$
denotes adjoint with respect to the inner product defined previously and
$\boldsymbol{W}_c^q$
is the controllability Gramian that characterises the reachability of
$\boldsymbol{q}$
in state space by a suitable action of forcing
$\boldsymbol{f}$
. For detailed derivations, the reader may refer to Farrell & Ioannou (Reference Farrell and Ioannou1993) and Jovanović & Bamieh (Reference Jovanović and Bamieh2005). Similarly to POD modes, the following eigenvalue problem produces a set of orthogonal basis functions that are ranked in terms of the most responsive velocity field
$\boldsymbol{u}$
to the white-in-time forcing:
\begin{equation} \boldsymbol{W}_c \varPhi = \varPhi \varLambda . \end{equation}
In the present study, we shall refer to these orthogonal basis functions in
$\varPhi$
as ‘controllability modes’. We note that if (2.4) can generate
$\boldsymbol{u}$
whose covariance operator is identical to that from numerical simulations or experiments, the POD and controllability modes become identical.
2.3. Balanced truncation modes
Instead of computing the most energetic structure in the output, balanced truncation produces a set of biorthgonal modes that are ranked in terms of dynamical importance (Moore Reference Moore1981; Zhou et al. Reference Zhou, Doyle and Glover1996; Willcox & Peraire Reference Willcox and Peraire2002). To this end, we first introduce the observability Gramian
$\boldsymbol{W}_o^q$
. Using (2.5), it can be shown that the response of
$\boldsymbol{q}$
driven by an impulse of
$\boldsymbol{f}$
applied at
$t=0$
is characterised by the observability Gramian
$\boldsymbol{W}_o^q$
, which is obtained by solving the following equation:
\begin{equation} \boldsymbol{A}^\dagger \boldsymbol{W}_o^q + \boldsymbol{W}_o^q \boldsymbol{A}+ \boldsymbol{C}^\dagger \boldsymbol{C} = 0, \end{equation}
where
$\boldsymbol{W}_o^q$
characterises the forcing structure
$\boldsymbol{B}\boldsymbol{f}$
for the most energetic (or responsive)
$\boldsymbol{u}$
. In terms of forcing
$\boldsymbol{f}$
,
$\boldsymbol{W}_o = \boldsymbol{B}^\dagger \boldsymbol{W}_o^q \boldsymbol{B}$
is obtained instead.
The controllability Gramian characterises the velocity structure of the most responsive
$\boldsymbol{u}$
, whereas the observability Gramian is for the structure of the forcing
$\boldsymbol{f}$
associated with the most responsive
$\boldsymbol{u}$
. They typically appear in different forms, especially when the Reynolds number is high (see the discussion in § 5.3). To account for this difference, the balanced truncation introduces a coordinate transformation
$\boldsymbol{T}$
, defined by
$\boldsymbol{u}=\boldsymbol{T}\boldsymbol{z}$
, and subsequently seeks the
$\boldsymbol{T}$
that makes
$\boldsymbol{W}_c$
and
$\boldsymbol{W}_o$
are identical (i.e. balancing):
\begin{align} \boldsymbol{W}_c &\mapsto \boldsymbol{T}^{-1} \boldsymbol{W}_c \left(\boldsymbol{T}^{-1}\right)^H, \end{align}
\begin{align} \boldsymbol{W}_o &\mapsto \boldsymbol{T}^H\boldsymbol{W}_o \boldsymbol{T}. \end{align}
It can be shown that such a coordinate exists, and the controllability and observability Gramian operators in the transformed coordinate are diagonal. In practice, the transformation
$\boldsymbol{T}$
is obtained by computing the right eigenvector of
$\boldsymbol{W}_c\boldsymbol{W}_o$
:
\begin{equation} \boldsymbol{W}_c\boldsymbol{W}_o\boldsymbol{T} = \boldsymbol{T}\varLambda , \end{equation}
where the transformation
$\boldsymbol{T}$
provides a set of basis functions (i.e.
$\varPhi =\boldsymbol{T}$
) known as balanced truncation modes and the diagonal entries in
$\varLambda$
are the eigenvalues that provide their dynamical importance. We direct the reader to Moore (Reference Moore1981) and Willcox & Peraire (Reference Willcox and Peraire2002) for detailed derivation. Unlike the POD and controllability modes, the balanced truncation modes are not orthogonal to each other. Therefore, to apply a Galerkin projection with them, another set of functions that are orthogonal to the balanced truncation modes is required. This can be obtained by calculating the left eigenvector of
$\boldsymbol{W}_c\boldsymbol{W}_o$
, denoted by
$\boldsymbol{S}(\equiv \varPsi =[ \psi _1(\boldsymbol{x})\,\psi _2(\boldsymbol{x})\, \ldots ])$
. After a suitable normalisation of
$\phi _i$
and
$\psi _i$
, the biorthogonality between
$\varPhi$
and
$\varPsi$
can be set to be satisfied:
$\langle \phi _i, \psi _j \rangle = \delta _{ij}$
.
3. Reduced-order models for plane Couette flow
As a benchmark case, we consider plane Couette flow of incompressible fluid with density
$\rho$
and kinematic viscosity
$\nu$
. The two parallel walls move in a direction opposite to each other with a streamwise velocity of
$\pm U_w$
and are separated by a distance
$2h$
. All variables are dimensionless with
$U_w$
,
$h$
and
$\nu$
, such that the resulting Reynolds number is defined to be
$Re=U_wh/\nu$
. The dimensionless streamwise, wall-normal and spanwise directions are denoted by
$x$
,
$y$
and
$z$
, respectively.
To compute POD modes and assess the performance of the ROMs, a direct numerical simulation (DNS) is performed using an open source code, channelflow (Gibson, Halcrow & Cvitanović Reference Gibson, Halcrow and Cvitanović2008). The code employs the Fourier–Galerkin method for the discretisation in the streamwise and spanwise directions, and the Chebyshev-tau method in the wall-normal direction. The computational domain is defined as
$(x,y,z) \in [0, 1.75\pi )\times [-1,1]\times [0, 1.2\pi )$
and is close to the minimal flow unit of Hamilton et al. (Reference Hamilton, Kim and Waleffe1995). The resolution of the grid is
$(N_x, N_y, N_z) = (32, 65, 32)$
after dealiasing. The Reynolds number considered is
$Re=500$
and the corresponding friction Reynolds number is obtained as
(
$u_\tau$
is the friction velocity). The DNS is integrated forward in time, covering a time span of
$t\in [0,5000]$
, with 10 000 snapshots taken evenly in time.
3.1. POD modes
The POD modes in plane Couette flow are computed following the procedure described by Smith et al. (Reference Smith, Moehlis and Holmes2005) and Chua Khoo et al. (2022). Since the flow is translationally invariant in the streamwise and spanwise directions, the POD modes are given in the form of the Fourier modes in those directions:
\begin{equation} \phi _i(x,y,z)=\frac {1}{\sqrt {L_x L_z}}\hat {\phi }_{n_x,n_z}^{(n)}(y)\textrm {exp}(\textrm {i}( n_x \alpha _0 x+n_z\beta _0 )), \end{equation}
where
$n=1,2,\ldots ,N_y$
,
$\alpha _0=2\pi /L_x$
,
$\beta _0=2\pi /L_z$
and (
$n_x$
,
$n_z$
) is a pair of integers indicating each Fourier mode. Here, the index
$i$
for the
$i$
th POD mode on the left-hand side can now be suitably defined in terms of
$(n_x,n_z,n)$
. The covariance operator in (2.3a
) can then be formulated in terms of individual Fourier modes, significantly reducing its size. This covariance operator is also computed by imposing the following discrete symmetries in plane Couette flow:
\begin{align} \mathcal{I}(u,v,w,p)[x,y,z,p] &= (u,v,w,p)[x,y,z,p], \end{align}
\begin{align} \mathcal{P}(u,v,w,p)[x,y,z,p] &= (-u,-v,-w,p)[-x,-y,-z,p], \end{align}
\begin{align} \mathcal{R}(u,v,w,p)[x,y,z,p] &= (u,v,-w,p)[x,y,-z,p], \end{align}
\begin{align} \mathcal{RP}(u,v,w,p)[x,y,z,p] &= (-u,-v,w,p)[-x,-y,z,p], \end{align}
where
$p$
is pressure. Here, the group
$\mathcal{I}$
indicates the identity transformation,
$\mathcal{P}$
is a point reflection about the origin,
$\mathcal{R}$
is the rotation about the plane of
$z=0$
and
$\mathcal{RP}$
is a 180
$^\circ$
rotation about the
$z$
-axis. For further details, the reader may refer to Smith et al. (Reference Smith, Moehlis and Holmes2005) and Chua Khoo et al. (2022).
3.2. Linear operators
For the computation of the controllability and balanced modes, here we define the linear operators introduced in § 2.2. We first decompose the flow field into a base flow
$\boldsymbol{U} = (U(y),0,0)$
and perturbation/fluctuation
$\boldsymbol{u} = (u,v,w)$
. The dimensionless momentum equations for
$\boldsymbol{u}$
are written as
\begin{equation} \frac {\partial \boldsymbol{u} }{\partial t} + (\boldsymbol{U} \boldsymbol{\cdot } \boldsymbol{\nabla }) \boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla }) \boldsymbol{U} = -\boldsymbol{\nabla }p + \frac {1}{\textit{Re}}{\nabla} ^2 \boldsymbol{u}+\boldsymbol{N}, \end{equation}
where
$\boldsymbol{N}$
is the nonlinear term given depending on the definition of
$\boldsymbol{U}$
. Here, the ensemble average
$\overline {(\boldsymbol{\cdot })}$
is equivalent to an average in time, and in the streamwise and spanwise directions, as the given flow is statistically stationary and homogeneous in these wall-parallel directions.
To model (3.3) in terms of a linear dynamical system, we consider various combinations of different base flows
$\boldsymbol{U}$
and models for the nonlinear term
$\boldsymbol{N}$
such that
\begin{equation} U(y)=\left \{ \!\!\begin{array}{ll} U_0(y) \\[5pt] \overline {U}(y)\\ \end{array}\right . \quad \textrm {and} \quad \boldsymbol{N}=\left \{\!\! \begin{array}{ll} \boldsymbol{f} \\[5pt] \boldsymbol{\nabla }\boldsymbol{\cdot }\left(\nu _t\left(\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^T\right)\right)+\boldsymbol{f},\\ \end{array}\right . \end{equation}
where
$U_0(y)(=y)$
is the laminar base flow,
$\overline {U}$
the turbulent mean velocity obtained from DNS and
$\boldsymbol{f} = [f_u,f_v,f_w]$
the forcing term providing a white noise in time. Here, the dimensionless eddy viscosity
$\nu _t$
in
$\boldsymbol{N}$
is given by a simple mixing length model (Reynolds & Hussain Reference Reynolds and Hussain1972; Hwang & Cossu Reference Hwang and Cossu2010a
):
\begin{equation} -\overline {u v} = \nu _t \frac {{\rm d}\overline {{U}}(y)}{{\rm d}y}. \end{equation}
Given the homogeneity in the streamwise and spanwise directions, we consider a Fourier mode representation of the velocity fluctuation/perturbation and forcing:
$\hat {\boldsymbol{u}}(y, t; \alpha , \beta ) \text{e}^{\text{i}(\alpha x + \beta z)}$
and
$\hat {\boldsymbol{f}}(y, t; \alpha , \beta ) \text{e}^{\text{i}(\alpha x + \beta z)}$
, where
$ \alpha = n_x\alpha _0$
and
$\beta = n_z\beta _0$
. Substituting these Fourier modes into (3.3) and eliminating the pressure fluctuation
$p$
with the standard method, we obtain the following Orr–Sommerfeld–Squire system in terms of the Fourier modes of the wall-normal velocity and vorticity,
$\hat {v}$
and
$\hat {\eta }$
:
\begin{equation} \frac {\partial }{\partial t}\underbrace {\begin{bmatrix} \hat {v}\\ \hat {\eta } \end{bmatrix}}_{\hat {\boldsymbol{q}}} = \underbrace { \begin{bmatrix} \Delta ^{-1}\mathcal{L}_{OS} & 0 \\ -i\beta U^\prime & \mathcal{L}_{SQ} \end{bmatrix}}_{\boldsymbol{A}} \begin{bmatrix} \hat {v} \\ \hat {\eta } \end{bmatrix} + \underbrace { \begin{bmatrix} -i\alpha \Delta ^{-1} \mathcal{D} & -k^2\Delta ^{-1} & -i\beta \Delta ^{-1}\mathcal{D} \\ i\beta & 0 & -i\alpha \\ \end{bmatrix} }_{\boldsymbol{B}} \underbrace {\begin{bmatrix} \hat {f_u} \\[3pt] \hat {f_v} \\[3pt] \hat {f_w} \\ \end{bmatrix}}_{\hat {\boldsymbol{f}}}, \end{equation}
where
\begin{align} \mathcal{L}_{OS} &= -i\alpha (U\Delta - \mathcal{D}^2 U) + \nu _T \Delta ^2 + 2\mathcal{D}\nu _t\Delta \mathcal{D} + \mathcal{D}^2\nu _t(\mathcal{D}^2 + k^2), \end{align}
\begin{align} \mathcal{L}_{SQ} &= -i\alpha U + \nu _T \Delta + \mathcal{D}\nu _t \mathcal{D} \end{align}
with the boundary conditions,
$\hat {v}(y=\pm 1) =\mathcal{D}\hat {v}(y = \pm 1) = 0$
and
$\hat {\eta }(y=\pm 1) = 0$
. Here,
$\mathcal{D}\equiv \partial / \partial y$
,
$\Delta = \mathcal{D}^2 - k^2$
and
$k^2 = \alpha ^2 + \beta ^2$
. Also,
$\nu _T$
represents the total dimensionless eddy viscosity defined as
$\nu _T=1/Re+\nu _t$
and the terms involving
$\nu _t$
in (3.5) appear if
$\boldsymbol{N}$
contains the related terms (see (3.4a)). Using
$\hat {\eta } = \partial \hat {{u}} / \partial z - \partial \hat {{w}}/{\partial x}$
, the primitive variables are recovered with
\begin{equation} \underbrace {\begin{bmatrix} \hat {u} \\ \hat {v} \\ \hat {w} \\ \end{bmatrix}}_{\hat {\boldsymbol{u}}} = \underbrace {\frac {1}{k^2} \begin{bmatrix} i\alpha \mathcal{D} & -i\beta \\ k^2 & 0 \\ i\beta \mathcal{D} & i\alpha \\ \end{bmatrix} }_{\boldsymbol{C}} \begin{bmatrix} \hat {v} \\ \hat {\eta } \end{bmatrix} . \end{equation}
By replacing
$\boldsymbol{u}$
and
$\boldsymbol{q}$
in (2.4) with
$\hat {\boldsymbol{u}}$
and
$\hat {\boldsymbol{q}}$
, it is now possible to generate the controllability and balanced modes in the form of (3.1).
The eddy viscosity
$\nu _t(y)$
is obtained using (3.4b
), but the mean velocity profile there was obtained from DNS with a slightly different domain size given by
$(x,y,z) \in [0, 2\pi )\times [-1,1]\times [0, \pi )$
. This is to avoid some unphysical behaviour in the mean velocity profile of the minimal flow unit, where
${\rm d}U/{\rm d}y\lt 0$
was found around
$y=0$
due to the highly restricted spatial domain. We note that in such a case, the eddy viscosity calculated by (3.4b
) yields singularities at the wall-normal locations where
${\rm d}\overline {U}/{\rm d}y=0$
, an undesirable feature in terms of smoothness of
$\nu _t(y)$
to be modelled. For the computation of the controllability and balanced modes, the wall-normal direction of (3.5) is discretised using a Chebyshev collocation method described by Weideman & Reddy (Reference Weideman and Reddy2000) with
$N_y=65$
. The Lyapunov equations in (2.8) and (2.10) are solved using the lyap function in MATLAB. Lastly, for
$\alpha = \beta = 0$
, which corresponds to the flow averaged in the streamwise and spanwise directions, the linear operator
$\boldsymbol{A}$
in (3.5) is singular. In this case, we instead used a set of the Stokes modes, the orthogonal basis functions obtained from the viscous diffusion operator, as was done by Cavalieri & Nogueira (Reference Cavalieri and Nogueira2022) and McCormack et al. (Reference Maia and Cavalieri2024).
List of basis functions used for ROMs in the present study.
3.3. Galerkin projection
Using each set of POD modes from the turbulent state and controllability/balanced modes from various forms of (3.5) depending on the base flow
$U$
and the nonlinear term model
$\boldsymbol{N}$
, we apply a Galerkin projection following the procedure described by Cavalieri & Nogueira (Reference Cavalieri and Nogueira2022) based on the Navier–Stokes equations for flow fluctuations around the laminar solution
$U_0$
. The list of basis functions and their construction approaches are summrised in table 1. Notice that although different base flows were considered to obtain the basis functions, the Galerkin projection is constructed here based on fluctuations around the laminar solution
$U_0$
; thus, the Galerkin procedure is the same for all cases studied in table 1 and the only change between ROMs is the choice of basis functions. For the Galerkin projection, the test functions used for the POD and controllability modes were themselves, resulting in
\begin{align} \underbrace {\frac {{\rm d}{a}_i}{{\rm d}t}}_{{\rm d}\boldsymbol{a}/{\rm d}t} & = \underbrace { \frac {1}{\textit{Re}} \sum _j \langle {\nabla} ^2 \phi _j, \phi _i \rangle a_j- \sum _j \langle [(\phi _j \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U_0} + (\boldsymbol{U_0} \boldsymbol{\cdot }\boldsymbol{\nabla })\phi _j] , \phi _i \rangle a_j}_{\boldsymbol{La}}\nonumber \\ &\quad - \underbrace { \sum _j \sum _k \langle (\phi _j \boldsymbol{\cdot }\boldsymbol{\nabla })\phi _k,\phi _i \rangle a_j a_k}_{\boldsymbol{N}(\boldsymbol{a})}, \end{align}
where
$\boldsymbol{a}=[a_1\,a_2\,a_3\,\ldots ]^T$
and
$\boldsymbol{U}_0=(U_0(y),0,0)$
. Here, similarly to (3.1), the index
$i$
of
$a_i$
can be expressed in terms of
$(n_x,n_z,n)$
for the
$n$
th modes obtained from the Fourier component of
$(n_x,n_z)$
: i.e.
$a_i=a_{n_x,n_z}^{(n)}$
. Note that (3.7) is formulated using the Navier–Stokes equations in a perturbation form about the laminar base flow. Alternatively, it may be formulated with turbulent mean flow instead of the laminar base flow. In such a case, the turbulent mean flow must be known in prior. However, not all the Galerkin models in this study require such a requirement and we do not assume either that the very precise mean flow is available. Therefore, this is not considered. Finally, for the ROM based on balanced truncation, the test functions were replaced with the left eigenvectors of
$\boldsymbol{W}_c\boldsymbol{W}_o$
,
$\boldsymbol{S}$
(or
$\boldsymbol{\varPsi }$
), defined in § 2.3. For further details on the Galerkin projection, the reader is referred to Cavalieri & Nogueira (Reference Cavalieri and Nogueira2022).
4. Results
We present this section by dividing it into two subsections: one is on the convergence and performance of the ROMs around the laminar state, in which the linear stability and optimal transient growth computed with the ROMs are compared with those obtained with the full equations of motion, and the other is on their convergence and performance for the turbulent state, where the turbulence statistics and the coherent structure dynamics, evaluated by temporal correlations, are examined against the corresponding DNS results.
The computational domain of ROMs in the present study is the same as that of DNS described in § 3:
$(x,y,z) \in [0, 1.75\pi )\times [-1,1]\times [0, 1.2\pi )$
. Throughout the present study, we consider a single non-zero Fourier mode in the streamwise diection and two modes in the spanwise directions for the given computational domain: i.e.
$n_x=0,\pm 1$
and
$n_z=0, \pm 1, \pm 2$
in (3.1). The number of DoFs of each of the ROMs is varied by changing the number of modes in the wall-normal direction,
$N_y$
, for each plane Fourier component. Here, we consider
$N_y = 5,10,20,30$
. The corresponding total DoFs are
$N_{{tot}}=75,150,300,450$
, respectively.
4.1. Laminar state
We first examine the ROMs around the laminar base state, focussing on two pairs of wavenumbers,
$(n_x,n_z) = (0,1)$
and
$(n_x,n_z)= (1,1)$
. These two pairs of wavenumbers were chosen due to their physical importance in the dynamics of the transition to turbulence (Reddy et al. Reference Reddy, Schmid, Baggett and Henningson1998). The
$(0,1)$
mode provides the maximum (linear) transient growth associated with streaks, a key physical element of transition, and the excitation of the
$(1,1)$
mode leads to an oblique transition, which is energetically more efficient than the transition through the excitation of the
$(0,1)$
mode. We also note that the oblique transition scenario has recently been shown to be initiated by the Orr mechanism (Jiao, Hwang & Chernyshenko Reference Jiao, Hwang and Chernyshenko2021) and is physically very similar to the transition triggered by the so-called ‘minimal seed’ (i.e. the minimal energy initial condition for transition; for a review, see Kerswell Reference Kerswell2018).
Figure 1 shows the eigenvalues of the linearised operator for the two pairs of wavenumbers, on increasing the number of degrees of freedom in the wall-normal direction,
$N_y$
. Specifically, the linear operator for each pair of wavenumbers is extracted from
$\boldsymbol{L}$
in (3.7) to study its linear stability and here we note that the laminar base state in plane Couette flow is linearly stable. In all cases, the eigenvalues from the ROMs gradually converge to those obtained directly from the linearised Navier–Stokes (or the Orr–Sommerfeld–Squire) operator, as
$N_y$
increases, and in particular, the first couple of least stable eigenvalues show excellent convergence. In the cases of C-LNL and BT-LNL, the degrees to which the eigenvalues converge as
$N_y$
increases are also seen to be consistent with those reported by Ilak & Rowley (Reference Ilak and Rowley2008): for example, compare figure 1 with their figure 5.
Eigenvalues of the linearised operator about the laminar base state: (a–d)
$(n_x,n_z) = (0,1)$
; (e–h)
$(n_x,n_z) = (1,1)$
. Here, (a,e)
$N_y=5$
; (b, f)
$N_y=10$
; (c,g)
$N_y=20$
; (d,h)
$N_y=30$
. Only top half of the complex plane is shown, as the eigenvalues are symmetric about
$\textrm {Imag}(\lambda )=0$
. In each case, dashed line indicates
$\text{Real}(\lambda ) = 0$
. Here, the symbols denoted by OSQ represent the eigenvalues directly computed from the Orr–Sommerfeld–Squire system operator.
Of all the cases tested, it was found that the cases C-LNL and BT-LNL, where the basis functions are obtained from the linearised operator with the laminar base flow and molecular viscosity, do not generate unstable eigenvalues for all examined
$N_y$
, indicating the preservation of the linear stability of the laminar base state even with very severe truncations. However, this was not necessarily the case for the ROMs with all the other basis functions, which are obtained with some information from turbulent flows. In this case, at
$N_y \lesssim 20$
, the laminar base state of the ROMs often produced linear instabilities, which are not desirable for modelling the state space around the laminar base flow: for example,
$(n_x,n_z)=(0,1)$
with
$N_y=5$
(figure 1
a) and
$(n_x,n_z)=(1,1)$
with
$N_y=5,10$
(figures 1
e and 1
f). When increasing
$N_y$
further, the first few eigenvalues are gradually captured well by the ROMs (figure 1
d, f). Eventually, many leading eigenvalues of all ROMs converged to those of the linearised Navier–Stokes (or the Orr–Sommerfeld–Squire) operator with excellent accuracy when
$N_y\gtrsim 100$
(not shown).
Lastly, table 2 summarises the minimum
$N_y$
required to retain the linear stability of the laminar base state for all
$(n_x,n_z)$
incorporated in the present ROMs. In all wavenumber pairs examined, cases C-LNL and BT-LNL consistently reproduced the linear stability of the laminar base state only with
$N_y=1$
. However, it was found that other ROMs require a substantially large value of
$N_y$
for that.
The minimum
$N_y$
greater than or equal to which laminar base state becomes stable.
We now consider the optimal transient growth computed from the ROMs for each wavenumber (e.g. Schmid & Henningson Reference Schmid and Henningson2001). The optimal transient growth is defined as
\begin{equation} G(t{;}n_x,n_z) \equiv \lVert \boldsymbol{C}e^{\boldsymbol{L}t}\boldsymbol{B} \rVert ^2=\max _{\boldsymbol{u}_0\not = 0} \frac {\lVert \boldsymbol{u} (t{;}n_x, n_z) \rVert ^2}{\lVert \boldsymbol{u}_0 (n_x, n_z) \rVert ^2}, \end{equation}
where
$\boldsymbol{u}_0$
is the initial condition. Here,
$\|(\boldsymbol{\cdot })\|$
is the standard energy norm (Schmid & Henningson Reference Schmid and Henningson2002). We note that optimal transient growth
$G$
is also the square of the energy norm of the state transition operator evaluated at every
$t$
, as shown previously.
The optimal transient growth for the linearised Navier–Stokes operator and the ROMs was calculated using the standard method (Schmid & Henningson Reference Schmid and Henningson2002): for the linearised operator, the eigenfunction expansion technique was used and, for the ROMs, the projected linear operator
$\boldsymbol{L}$
was used. The transient growth is reported in figure 2 for the two wavenumber pairs considered:
$(n_x,n_z) = (0,1)$
and
$(n_x,n_z)= (1,1)$
. It is evident that the best performing ROM is the case BT-LNL (the ROM using balanced truncation modes obtained with laminar base flow and molecular viscosity), consistent with the previous finding of Ilak & Rowley (Reference Ilak and Rowley2008). The C-LNL case also performs fairly well, as it only needs
$N_y \approx 10$
to capture the transient growth reasonably well. Similarly to the eigenvalue assessment in figure 1, the ROMs, whose basis functions are constructed using the information from turbulent flows, do not perform well and this is evident from the results for
$N_y\lesssim 10$
. In all cases, the ROMs reproduced the optimal transient growth very well when
$N_y=30$
, although the C-LNTe and BT-LNTe cases still showed a slight deviation for
$(n_x,n_z)= (1,1)$
(figure 2
h).
Optimal transient growth around laminar base flow for (a–d)
$(n_x,n_z)=(0,1)$
and (e–h)
$(n_x,n_z)=(1,1)$
: (a,e)
$N_y=5$
; (b, f)
$N_y=10$
; (e,g)
$N_y=20$
; (d,h)
$N_y=30$
. Here, the lines missing indicate the cases that the corresponding ROM exhibit either rapid exponential growth (i.e. instability) or decay (not sufficient transient growth). In panels (c,d), the lines overlap with each other.
We also consider the response of the ROM to harmonic forcing for each wavenumber (Schmid & Henningson Reference Schmid and Henningson2002). The maximum response is measured by computing the square of the resolvent norm, such that
\begin{equation} R(\omega {;}n_x,n_z) \equiv \rVert \boldsymbol{C}(-i\omega \boldsymbol{I}-\boldsymbol{A})^{-1}\boldsymbol{B}\lVert ^2 = \max _{\hat {\boldsymbol{f}}\neq 0} \frac {\rVert \boldsymbol{C}(-i\omega \boldsymbol{I}-\boldsymbol{A})^{-1}\boldsymbol{B} \tilde {\boldsymbol{f}} \lVert ^2}{\lVert \tilde {\boldsymbol{f}}\lVert ^2}, \end{equation}
where the harmonic forcing is defined as
$\hat {\boldsymbol{f}} = \tilde {\boldsymbol{f}}e^{-i\omega t}$
. The same method described for transient growth was used for the computation of
$R(\omega {;}n_x,n_z)$
. We present the computed
$R(\omega {;}n_x,n_z)$
in figure 3 with the same sets of wavenumbers considered in figure 2. The trends in the maximum response to harmonic forcing closely resemble those observed in the transient growth results: the BT-LNL model exhibits excellent convergence at
$N_y=5$
, followed by the C-LNL model, which also performs reasonably well at the same
$N_y$
(figure 3
b, f). Consistent with the transient growth behaviour at
$(n_x, n_z) = (1,1)$
, the C-LNTe and BT-LNTe cases display slight deviations even at
$N_y = 30$
(figure 3
h).
Optimal harmonic response for (a–d)
$(n_x,n_z)=(0,1)$
and (e–h)
$(n_x,n_z)=(1,1)$
: (a,e)
$N_y=5$
; (b, f)
$N_y=10$
; (e,g)
$N_y=20$
; (d,h)
$N_y=30$
. In panels (b,c,d,h), the lines overlap with each other.
4.2. Turbulent state
We now examine the modelling capability of the ROMs for the turbulent state. The ROMs were integrated in time using the function ode45 in MATLAB. The initial condition for the ROMs in (3.7) is first generated as a random sequence with
$a_i \in (0,0.1)$
. Equation (3.7) is subsequently integrated to see if it reaches a statistically stationary state. When
$N_y=5$
, not all ROMs were able to generate a statistically stationary turbulent state. Only the cases C-LNL, C-LNT and POD can generate a statistically stationary turbulent state. The case BT-LNL, which was able to most effectively capture the stability of the laminar state, showed rapid relaminarisation, and BT-LNT, C-LNTe and BT-LNTe generated solutions that were rapidly blown up. Regardless of the stability of the laminar base state, when
$N_y=10$
, all ROMs were able to generate a stationary turbulent state. Finally, when
$N_y\gtrsim 20$
, all ROMs were able to retrieve the linearly stable laminar state (table 2) and a stationary turbulent state.
4.2.1. Turbulence statistics
Figure 4 reports the mean velocity profiles obtained by averaging over a period of time
$T_{{avg}}=10\,000$
. At the heaviest truncation (
$N_y=5$
; figure 4
a), the mean velocity profiles from the C-LNL, C-LNT and POD cases, which were able to generate a statistically stationary turbulent state, are considerably different from the one obtained from DNS. In particular, C-LNL and C-LNT show negative values of
${\rm d}U/{\rm d}y$
and the related overshoot in the mean velocity profile
$U(y)$
around
$y=0$
. We note that the mean profile of the C-LNL case appears to be better than that of the C-LNT case. However, this is likely a coincidence due to the non-monotonic convergence of the ROMs at very low values of
$N_y$
: for example, velocity fluctuations of the C-LNL case are seen to be not better than those of the C-LNT case for
$N_y=10$
(see figure 5
d–f). When
$N_y$
is doubled (
$N_y=10$
; figure 4
b), all ROMs generate a turbulent state, the mean of which appears to be very difficult to distinguish from each other. In this case, the wall-normal locations, at which
${\rm d}U/{\rm d}y\lt 0$
, are very small or do not exist. Therefore, there is no visible overshoot in
$U(y)$
. Finally, for
$N_y\gtrsim 20$
, the mean velocity profiles of all the ROMs appear to converge into a single one, which shows only a small difference from the mean velocity profile obtained with DNS. This difference originates from the modelling setting of the ROMs, where we only considered a limited number of Fourier modes in the wall-parallel directions. In fact, when the number of Fourier modes is further increased, the mean velocity profile of ROMs converges to that of DNS, as shown by Cavalieri & Nogueira (Reference Cavalieri and Nogueira2022) for the C-LNL case with the basis functions obtained at
$Re=100$
.
Mean streamwise velocity of turbulent state: (a)
$N_y=5$
; (b)
$N_y=10$
; (c)
$N_y=20$
; (d)
$N_y=30$
. In panel (a), the cases missing from the plot either show the solution rapidly blowing-up or relaminarisation (see text). In panels (c, d), the ROM profiles are hardly discernible from each other due to their convergence.
Root-mean-square (a,d,g,j) streamwise, (b,e,h,k) wall-normal and (c, f,i,l) spanwise velocity of turbulent state: (a–c)
$N_y=5$
; (d–f)
$N_y=10$
; (g–i)
$N_y=20$
; (j–l)
$N_y=30$
. In panel (a), the cases missing from the plot either show the solution rapidly blowing-up or relaminarisation (see text), except the C-LNT case that is not shown due to huge deviations from the DNS profiles.
Root-mean-square (r.m.s.) velocity profiles are also shown in figure 5. For
$N_y=5$
(figure 5
a–c), the case POD shows remarkably good performance even with such a severe truncation. This is particularly true, as the profiles from the cases C-LNL and C-LNT significantly deviate from those obtained with DNS (the case C-LNT is not shown in figure 5(a) due to its huge deviation). When
$N_y=10$
(figure 5
d–f), the r.m.s. profiles of all ROMs begin to predict reasonably well compared with those of DNS. The POD case still predicts the r.m.s. profiles well, and the ROMs constructed using the mean velocity and the corresponding eddy viscosity (i.e. the cases C-LNTe and BT-LNTe) also perform almost equally well. All the other ROMs show a non-negligible difference in the r.m.s. profiles from those obtained with DNS, except the BT-LNT case, the prediction of which is only slightly worse than the POD, C-LNTe and BT-LNTe cases. Finally, at
$N_y\gtrsim 20$
(figure 5
g–l), the velocity fluctuations of all ROMs converge to a profile that deviates slightly from the DNS statistics. The same trend was observed in Reynolds shear stress, as reported in figure 6: at
$N_y = 10$
, the eddy viscosity significantly improves statistics (figure 6
b) and at
$N_y \gtrsim 20$
, the convergence is achieved. Incorporating additional streamwise and spanwise Fourier harmonics in (3.1) improves the convergence of rms velocity profiles, as already discussed for mean velocity. The convergence of mean and r.m.s. velocity profiles is shown in Appendix A.
Reynold stress of turbulent state: (a)
$N_y=5$
; (b)
$N_y=10$
; (c)
$N_y=20$
; (d)
$N_y=30$
. In panel (a), the cases missing from the plot either show the solution rapidly blowing-up or relaminarisation (see text), except the C-LNT case that is not shown due to huge deviation from the DNS profile. In panels (c, d), the ROM profiles are hardly discernible from each other due to their convergence.
Overall, the observations here highlight that the ROMs using the basis functions from the linearised Navier–Stokes operator with turbulent mean velocity and the eddy viscosity model (C-LNTe and BT-LNTe) generally perform well in predicting turbulence statistics. Their performance is comparable to that of the ROM with POD modes and superior to those using the modes obtained from the linearised Navier–Stokes operator without the eddy viscosity model. This suggests that adding the eddy viscosity model to the linearised Navier–Stokes operator would be useful for the generation of basis functions, especially at high Reynolds numbers, as they would potentially be able to keep the dimension of Galerkin-projection-based ROMs relatively low.
4.2.2. Dynamics of coherent structures
Coherent structure dynamics in the minimal flow units of plane Couette flow has been well understood, since the pioneering analysis by Hamilton et al. (Reference Hamilton, Kim and Waleffe1995). Figure 7 explains the coherent structure dynamics, often referred to as the ‘self-sustaining process’. There, the time evolution of the key observables is presented in figure 7(a,b) and the corresponding flow visualisations are shown in figure 7(c–g). The flow field is dominated by elongated streamwise velocity structures called streaks (figure 7
c), the spanwise spacing of which is captured by
$L_z$
. The streaks are therefore well characterised by the energy of the plane Fourier mode, defined by
\begin{equation} M(t{;}n_x,n_z) = \int _{-1}^{1} \tilde {\boldsymbol{u}}^H(t,y;n_x,n_z)\tilde {\boldsymbol{u}}(t,y;n_x,y,n_z)\,{\rm d}y \end{equation}
for
$(n_x,n_z)=(0,1)$
. The streaks gain their energy from mean shear and are amplified. When their energy is high enough (top one in figure 7
c), they experience instability and/or transient growth, both typically involving sinuous meandering motion with the fundamental streamwise wavelength in the minimal flow unit (figure 7
d). This meandering motion can be characterised by
$M(t,n_x,n_z)$
with
$(n_x,n_z)=(1,1)$
. After experiencing the meandering motion, the streaks break down and reach a low-energy state, while generating streamwise vortices (figure 7
e). The streamwise vortices regenerate the streaks (figure 7
f,g), and the physical processes described here repeat with the typical time scale of
$T=80 {-} 100$
in the minimal flow unit of plane Couette flow. To quantitatively assess these temporal dynamics of the coherent structures in the ROMs, we introduce the temporal auto- and cross-correlations between
$M(t;n_x,n_z)$
:
\begin{equation} C_{(m_x,m_z)}^{(n_x,n_z)}(\tau ) = \frac {\overline { \tilde {M}(t+\tau {;}n_x,n_z) \tilde {M}(t{;}m_x,m_z)}}{\sqrt {\tilde {M}^2(t{;}n_x,n_z)}\sqrt {\tilde {M}^2(t{;}m_x,m_z)}}, \end{equation}
where
$\tilde {M}(t{;}n_x,n_z) = M(t{;}n_x,n_z) - \overline {M(t{;}n_x,n_z)}$
is the fluctuation around the mean.
(a) Time trajectories of
$M(0,1)$
and
$M(1,1)$
and (b) their magnified view for
$t \in [3880,3980]$
in DNS: blue,
$M(0,1)$
; orange,
$M(1,1)$
. (c–g) Flow field snapshots at different stages of SSP, with blue/red isosurfaces indicating
$u = \pm 0.65$
, respectively. Here, the streamlines on
$x = 0$
and
$x=1.75\pi$
indicate the spanwise and wall-normal velocity fields.
Figure 8 shows the correlations evaluated from the ROMs and subsequently compares them with those from DNS using
$M(t;0,1)$
(streaks) and
$M(t;1,1)$
(streak meandering motions). At the heaviest truncation (
$N_y=5$
; figure 8
a–c), among the three ROMs that exhibit no numerical instability (POD, C-LNL, C-LNT), the POD case was found to best capture the time scale of correlations obtained with DNS. In particular,
$C_{(0,1)}^{(1,1)}$
was well reproduced with such a heavy truncation and the negative values of
$C_{(0,1)}^{(1,1)}$
for
$\tau \lesssim 0$
, which represent the phase difference between
$M(0,1)$
and
$M(1,1)$
, were also observed. In contrast, the other two cases show auto-correlations with very short time scales compared with DNS and their cross-correlations appear to be largely unphysical (e.g. the behaviour of correlation tails around
$\tau \simeq \pm 40$
).
Auto- and cross-correlations of
$M(0,1)$
and
$M(1,1)$
for different ROMs: (a–c)
$N_y=5$
; (b–f)
$N_y=10$
; (g–i)
$N_y=20$
; (j–l)
$N_y=30$
.
As the number of degree of freedom increases to
$N_y=10$
(figure 8
d–f), the auto-correlations of the ROMs appear to improve, in addition to the fact that all ROMs can now generate a turbulent solution. In particular, the cases POD, C-LNTe and BT-LNTe show relatively well-converged auto-correlations (figure 8
d,e) compared with those for
$N_y=30$
(figure 8
j,k). However, the negative values of
$C_{(0,1)}^{(1,1)}$
for
$\tau \lesssim 0$
are still not seen for all ROMs, except the POD case, suggesting that the ROM based on POD modes performs best in terms of capturing the coherent structure dynamics of turbulent flow.
Finally, when
$N_y\gtrsim 20$
, all the correlation functions of the ROMs appear to be reasonably well converged and their
$C_{(0,1)}^{(1,1)}$
now also shows negative correlations for
$\tau \lesssim 0$
, consistent with DNS (figure 8
g–l).
4.2.3. Short-time tracking
We also consider the short time tracking of modal energy of the
$x$
-independent streak mode (0,1). We take the instantaneous velocity field in figure 7(c) (t = 3890) as an initial condition and the ROMs are integrated forward in time. To obtain the initial conditions in terms of mode coefficients, we project the DNS velocity field into the basis function of interest. Therefore, the mode-coefficient trajectories in figure 9 are different because the basis functions for each projection are not the same. We do not expect long-time tracking of the solution, as the system is chaotic, with exponential divergence of neighbouring trajectories in state space. However, short-time tracking of DNS trajectories is an indication of the accuracy of the ROMs.
Trajectories of ROMs compared to DNS. The initial conditions are taken from point (c) in figure 7(b). Here, we omit
$N_y=5$
, as some ROMs show numerical instability or do not generate a turbulent state.
Figure 9 shows the time trajectories of the leading mode coefficient of
$(n_x,n_z)=(0,1)$
(i.e.
$a^1_{(0,1)}$
) for DNS and ROMs with different values of
$N_y$
, when integrated from the same initial condition. At
$N_y = 10$
, the POD and C-LNTe cases are found to be the best in the tracking of the trajectory (figure 9
a, f), where their leading mode coefficient follows reasonably well that of DNS for
$t\lesssim 15$
. Notably, the basis generated around the turbulent state with eddy viscosity gives overall a better performance at
$N_y = 10$
(figure 9
f,g). At
$N_y = 20$
, the other operator-based ROMs, which do not use the eddy viscosity, begin to show some improvements; although the C-LNT case still shows an immediate deviation from DNS trajectory, the BT-LNL and BT-LNT cases follow the fluctuation longer to
$t\sim 5$
. Finally, at
$N_y = 30$
, except for the BT-LNL and BT-LNT cases, all ROMs are able to track the trajectory of DNS well up to
$t\approx 15$
.
5. Discussions
Thus far, we have explored the performance of several ROMs constructed with various forms of basis functions: POD modes from turbulent flows, and controllability and balanced truncation modes of the linearised Navier–Stokes operator considering different base flows (laminar base flow or turbulent mean flow) and incorporation of an eddy viscosity model. It has become evident that the performance and convergence of these ROMs around both laminar and turbulent state depends on the state space, the information of which is used to generate the different sets of the basis functions. Around the laminar state, the cases C-LNL and BT-LNL, which use the controllability and balanced truncation modes of the Navier–Stokes operator linearised about the laminar base flow, do not exhibit any linear instabilities even with only a single mode in the wall-normal direction (
$N_y=1$
). In particular, the BT-LNL case rapidly captures the transient growth and harmonic forcing response, the linear route of the transition to turbulence, as
$N_y$
increases (
$N_y \lesssim 10$
), consistent with previous findings of Ilak & Rowley (Reference Ilak and Rowley2008). In contrast, all other ROMs require at least
$N_y \gtrsim 20$
to achieve the linear stability and transient growth of the laminar base state (table 2 and figure 2).
The turbulent state has been found to be best modelled with the ROM based on POD modes, which have often been criticised due to the nature accounting for only the energy content without considering the dynamics. It is also worth noting that the cases C-LNTe and BT-LNTe model turbulence statistics fairly well at relatively low
$N_y$
(e.g.
$N_y=10$
; figures 4, 5, 6 and 8), indicating that the modes from the eddy-viscosity-enhanced linearised Navier–Stokes operator with turbulent mean flow have some non-negligible advantages over the other operator-based modes. An advantage of the cases C-LNTe and BT-LNTe over the POD-based ROMs is that they do not require prior simulation or experimental data, as they are derived directly from the linearised Navier–Stokes system, with some elements from the turbulent state via the use of mean flow and eddy viscosity. However, unlike the laminar state, the advantage of using the balanced truncation modes does not appear to be evident, given the comparable performance of the C-LNTe and BT-LNTe cases in modelling the turbulent state. The advantage of using the balanced truncation modes appears more evident in the case BT-LNT compared with case C-LNT, but both of the cases, obtained with the linearised Navier–Stokes operator with turbulent mean flow, do not perform better than their counterparts with the eddy viscosity model (cases C-LNTe and BT-LNTe).
To quantify the above-mentioned discussion, we report the errors from different observables studied around both laminar and turbulent states in tables 3 and 4. Overall, it is seen that there is a good correlation between the states around which the ROMs perform well and those whose information was used to create the basis functions. The remaining part of the present study will be spent understanding where this correlation originates from a physical viewpoint. In particular, we shall focus on the dynamics around the laminar and turbulent state in relation to the flow structures of each mode, visualised in figures 10 and 11 for
$(n_x,n_z)=(0,1)$
and
$(n_x,n_z)=(1,1)$
, respectively.
Error measures between ROM and DNS around laminar state at
$N_y = 10$
(see Appendix B for their definitions). Here, N/A in
$\epsilon _G$
indicates the case with linear instability or rapid relaminarisation, thus the error in transient growth is not reported.
Error measures between ROM and DNS at
$N_y = 10$
for turbulent statistics and correlations (see Appendix B for their definitions).
Visualisation of the leading mode of different basis functions for (
$n_x,n_z) = (0,1)$
on the
$y$
–
$z$
plane. The positive/negative contours indicate positive/negative streamwise velocity, and the vectors represents the wall-normal and spanwise velocity components. Each mode is normalised by its energy and the size of the vectors are scaled to be the same for all the cases: (a) POD; (b) C-LNL; (c) C-LNT; (d) C-LNTe; (e) BT-LNL direct mode; ( f) BT-LNL adjoint mode; (g) BT-LNT direct mode; (h) BT-LNT adjoint mode; (i) BT-LNTe direct mode; (j) BT-LNTe adjoint mode.
Cross-streamwise visualisation of the leading mode of different basis functions for (
$n_x,n_z) = (1,1)$
on the
$x$
–
$y$
plane. The positive/negative contours indicate positive/negative spanwise velocity, and the vectors represents the streamwise and wall-normal velocity components. Each mode is normalised by its energy and the size of the vectors are scaled to be the same for all the cases: (a) POD; (b) C-LNL; (c) C-LNT; (d) C-LNTe; (e) BT-LNL direct; ( f) BT-LNL adjoint; (g) BT-LNT direct; (h) BT-LNT adjoint; (i) BT-LNTe direct; (j) BT-LNTe adjoint.
5.1. Laminar state
In a previous study by Smith et al. (Reference Smith, Moehlis and Holmes2005), where a ROM was built using a Galerkin projection of POD modes obtained from turbulent state in plane Couette flow, it was reported that a sufficiently large number of
$N_y$
is necessary for the ROM to have linear stability of the laminar base state. Appropriate capturing of the linear stability of the laminar state is crucial, especially from a dynamical systems viewpoint, because it is the most fundamental feature required to depict the correct state-space dynamics for the transition to turbulence in plane Couette flow: e.g. saddle-node bifurcation for the first non-trivial invariant solutions to the system (Nagata Reference Nagata1990), the edge of turbulence (Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006) and so on. For this reason, Smith et al. (Reference Smith, Moehlis and Holmes2005) subsequently introduced a decoupling of each of the POD modes for
$n_x=0$
or
$n_z=0$
, so that their streamwise and cross-streamwise components are assigned into two different orthogonal modes. This was to better capture the so-called ‘lift-up’ effect (especially for
$n_x=0$
and
$n_z\neq 0$
; see also Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1991), by which the streamwise component of velocity is amplified due to the energy from an interaction between mean shear and the cross-streamwise components. In doing so, they were able to restore the linear stability of the laminar state without considering large values of
$N_y$
.
In the present study, we have shown that the C-LNL and BT-LNL cases retrieve the linear stability of the laminar state only with a single mode in the wall-normal direction (
$N_y=1$
). This was possible without introducing any decoupling procedure used by Smith et al. (Reference Smith, Moehlis and Holmes2005). We note that the basis functions for the C-LNL and BT-LNL cases are directly associated with the response of the given system around the laminar state to a small white-noise forcing: the C-LNL modes are physically the POD modes of the linear white-noise response of the laminar base state, while the BT-LNL modes combine the C-LNL modes with the observability modes that capture the flow structures responsible for the C-LNL modes. Therefore, they are expected to capture well the flow structures associated with the perturbation dynamics and the related stable eigenmodes of the laminar state. Despite this expectation, the related performance appears to be remarkably superior to that of the other basis functions, which require
$N_y \gtrsim 20$
to achieve the linear stability of the laminar state.
This is also evident in figure 10, where the leading modes used for the ROMs in the present study are visualised for
$(n_x,n_z)=(0,1)$
. The leading C-LNL (figure 10
b) and BT-LNL modes (figure 10
e) are structured with highly amplified streamwise velocity alternating in the spanwise direction (i.e. streaks), and they are highly energetic around the channel centre. However, the streaks captured in the leading POD mode have their largest amplitude near the wall and they involve counter-rotating streamise vortices, the amplitude of which is considerably larger than that of the leading C-LNL (figure 10
b) and BT-LNL modes. This implies that the most energetic modes relevant to turbulent state are quite different from the most energetic or dynamically relevant modes associated with the laminar base state even at the present Reynolds number (
$Re=500$
or
$Re_\tau = 37$
), which may well be viewed as low. Furthermore, the fact that the leading C-LNL and BT-LNL modes contain considerably smaller cross-streamwise velocity components relative to the POD mode supports the decoupling approach of Smith et al. (Reference Smith, Moehlis and Holmes2005) employed to restore the linear stability of the laminar base state: the mode containing only the streamwise component due to their decoupling would presumably be more appropriate to model the perturbation dynamics energetically dominated by the streamwise velocity like the leading C-LNL and BT-LNL modes.
5.2. Balanced truncation modes revisited
Controllability modes capture the most energetic flow structures from the response of the given linear system to a white-noise forcing, whereas observability modes capture the structures leading to the most energetic flow response. As discussed in § 2.3, the balanced truncation approach considers a trade-off between controllability and observability of the given (linear) system, providing a ROM that can better describe the dynamical behaviour of the full system only with small DoF (Rowley Reference Rowley2005). This is also very well supported by the BT-LNL case, which performs better than any other ROMs around the laminar state (see § 4.1). Visualisation of the balanced truncation modes for (
$n_x,n_z$
) = (0,1) and
$(n_x,n_z)=(1,1)$
in figures 10 and 11 provides further physical insight into why the balanced truncation is so effective around the laminar case. For
$(n_x,n_z)=(0,1)$
, all the direct balanced truncation modes (the right eigenvector of
$\boldsymbol{W}_c\boldsymbol{W}_o$
), which were used as basis functions for the corresponding ROMs, are dominated by streamwise streaks (figure 10
e,g,i), while the adjoint balanced truncation modes (the left eigenvector of
$\boldsymbol{W}_c\boldsymbol{W}_o$
) are mainly composed of counter-rotating streamwise vortices (figure 10
f,h,j). This is reminiscent of previous studies where controllability and observability modes were shown to be dominated by streaks and counter-rotating streamwise vortices, respectively (e.g. Farrell & Ioannou Reference Farrell and Ioannou1993; Hwang & Cossu Reference Hwang and Cossu2010a
,
Reference Hwang and Cossub
), indicating that these balanced truncation modes contain physical elements to suitably capture the lift-up effect. It also explains well why the BT-LNL case captures the optimal transient growth so well compared with all the other ROMS in addition to the fact that the related modes are sampled from the dynamics around the laminar state. A similar conclusion can be drawn for
$(n_x,n_z)=(1,1)$
, illustrated in figure 11. In this case, direct balanced truncation modes contain flow structures tilted downstream, while the corresponding adjoint modes have flow structures tilted upstream, indicating that they inherently contain suitable physical elements for the description of the Orr mechanism (Orr Reference Orr1907; Jiao et al. Reference Jiao, Hwang and Chernyshenko2021), the energy production mechanism of spanwise vortical structures due to tilting by base flow shear.
Despite the substantial advantages in the description on the dynamics around the laminar base state, the ROMs based on balanced truncation (cases BT-LNL, BT-LNT and BT-LNTe) were found to not perform substantially better than those based on controllability modes (cases C-LNL, C-LNT and C-LNTe) for concerning the description of a turbulent state. In fact, for turbulent states, it was found that the POD case consistently outperforms all the other ROMs. As discussed previously, the balanced truncation is able to account for the physical processes associated with the non-normality of the linearised Navier–Stokes equations, such as the lift-up effect and the Orr mechanism. It has been well understood that these processes are also highly relevant for the dynamics of coherent structures in turbulent flows: for example, the lift-up effect has been known to be the key process for the amplification of streaks (e.g. Butler & Farrell Reference Butler and Farrell1992; Hamilton et al. Reference Hamilton, Kim and Waleffe1995; del Á lamo & Jiménez 2006; Hwang & Cossu Reference Hwang and Cossu2010a
) and there has been considerable evidence on the importance of the Orr mechanism for coherent structures in turbulent flows (Jiménez Reference Jiménez2013; Enciner & Jiménez Reference Enciner and Jiménez2020). Furthermore, in the case of BT-LNTe, the leading (direct) balanced truncation mode for the most energetic plane Fourier mode
$(n_x,n_z)=(0,1)$
, which represents the streaks, compares fairly well with the leading POD mode of turbulent flow (compare figure 10
a with figure 10
i) – the spatial structure of the streaks and the relative strength cross-streamwise velocity components compare much better than those of the BT-LNL and BT-LNT cases (for a further discussion, see also § 5.3). However, for
$(n_x,n_z)=(1,1)$
, the leading direct balanced truncation mode of the BT-LNTe case is seen to be considerably different from the leading POD mode of turbulent flow (figure 11): the downstream tilt angle and the relative intensity of the spanwise velocity (figure 11
i), and the vector field indicating the wall-normal and spanwise velocity components of the BT-LNTe case are all seen to be considerably different from those of the POD mode (figure 11
a). In fact, in this case, the modes of the BT-LNL and BT-LNT cases are seen to compare better at least for the spanwise velocity, although their wall-normal velocity component also appears to be non-negligibly different from that of the POD mode. We note that the Fourier mode of
$(n_x,n_z)=(1,1)$
has been understood as a consequence of streak instability and/or transient growth (Hamilton et al. Reference Hamilton, Kim and Waleffe1995; Schoppa & Hussain Reference Schmid and Henningson2002). In other words, the relevant physical process for this mode would be depicted by the linearised Navier–Stokes equations about a suitable streaky base flow,
$\boldsymbol{U}_b(y,z)=\boldsymbol{U}(y)+\boldsymbol{u}_s(y,z)$
, where
$\boldsymbol{u}_s(y,z)$
represents the streaks (e.g. Figure 10
a). The leading POD mode for
$(n_x,n_z)=(1,1)$
would then reflect the flow structures associated with the streak instability and/or transient growth, a nonlinear process that cannot be modelled by the linearised Navier–Stokes operator about turbulent mean flow
$\boldsymbol{U}(y)$
with a simple white-noise forcing regardless of the presence of the eddy viscosity. This also explains why the POD case reproduces the qualitative behaviour of the correlation
$C_{(0,1)}^{(1,1)}$
at low values of
$N_y$
(
$N_y=5$
; figure 8
c), whereas all other ROMs are unable to do so even at
$N_y=10$
(figure 8
f) – this correlation is associated with the phase relation between the streaks and the streak instability/transition growth.
5.3. Linear operator modes with turbulent mean velocity and eddy viscosity
Although the ROM based on POD modes shows the best performance in modelling the turbulent state, the ROMs accounting for the information from the turbulent state perform substantially better than those without it. This is particularly true for the C-LNTe and BT-LNTe cases, where the linear operator employs both the turbulent mean flow and the corresponding eddy viscosity. As briefly discussed earlier, the benefit of using such a linear operator for the generation of basis functions is most apparent for
$(n_x,n_z)=(0,1)$
. This Fourier mode depicts the streamwise elongated structures that have been understood to be well modelled by the linearised Navier–Stokes equations with turbulent mean flow and eddy viscosity (Hwang & Cossu Reference Hwang and Cossu2010a
). In this case, the inertial part of the streamwise momentum equation is given by
\begin{equation} \frac {\partial \hat {u}}{\partial t} \sim -\hat {v}\frac {{\rm d}U}{{\rm d}y}, \end{equation}
implying that, for a given
$\hat {v}$
,
$\hat {u}$
would be more energetic at which
${\rm d}U/{\rm d}y$
is large. This explains why the streamwise velocity of the leading modes from the linear operators with turbulent mean velocity tends to be relatively strong in the region close to the wall (C-LNT, C-LNTe, BT-LNT and BT-LNTe in figure 10).
Furthermore, enhancing the linearised equations with the eddy viscosity model is seen to further improve the structure of the leading modes, making it more favourably comparable to the leading POD mode. The underlying reason can also be explained by the diffusion term of the eddy viscosity. In particular, the wall-normal transport of the streamwise Fourier component by this term is associated with the following part of the linearised Navier–Stokes equations:
\begin{equation} \frac {\partial \hat {u}}{\partial t} \sim \nu _t\frac {\partial ^2 \hat {u}}{\partial y^2}+\frac {{\rm d} \nu _t}{{\rm d}y}\frac {\partial \hat {u}}{\partial y}. \end{equation}
Here, the first term on the right-hand side represents the diffusion by eddy viscosity and the second term arises from the nature of the eddy viscosity varying in the wall-normal direction. Given the nature of the first term that would enhance the diffusion and weaken the non-normality of the linearised Navier–Stokes operator, the related leading modes are expected to be more isotropic. This is consistent with the structure of the leading C-LNTe and BT-LNTe (direct) modes for
$(n_x,n_z)=(0,1)$
, which contain relatively large cross-streamwise velocity components (figure 10
d,i) with amplitudes comparable to those of the leading POD modes (figure 10
a). The second term on the right-hand side of (5.2) is given in the form of a wall-normal advection with a local speed of
${{\rm d} \nu _t}/{{\rm d}y}$
, and this may be interpreted as a model for wall-normal turbulent (nonlinear) transport (see also Hwang Reference Hwang2016; Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023; Holford, Lee & Hwang Reference Holford, Lee and Hwang2024). Since
$\nu _t$
vanishes on the wall,
${{\rm d} \nu _t}/{{\rm d}y}\gt 0$
around the lower wall and
${{\rm d} \nu _t}/{{\rm d}y}\lt 0$
around the upper wall. This term would then create an advection flux towards the wall, which is why the streamwise velocity of the leading C-LNTe and BT-LNTe (direct) modes is more energetic than that of the other cases near the wall, similar to that of the leading POD modes.
Finally, it is worth mentioning that the balanced truncation modes of the linearised Navier–Stokes operator are non-orthogonal due to the advection term. Indeed, given the inner product defined for any solenoidal three-dimensional vector field in § 2,
$\boldsymbol{B}^\dagger =\boldsymbol{C}$
(
$\boldsymbol{C}^\dagger =\boldsymbol{B}$
, equivalently; see Holford Reference Holford2024). Also, consider that
$\boldsymbol{A}$
is the Stokes operator where the advection term is absent,
$\boldsymbol{A}^\dagger =\boldsymbol{A}$
. In such a case,
$\boldsymbol{W}_c=\boldsymbol{W}_o$
, which makes the controllability and observability modes identical. Consequently, for the Stokes operator,
$\boldsymbol{W}_c\boldsymbol{W}_o$
becomes self-adjoint and the resulting balanced truncation modes are orthogonal. From this observation, in general, the linearised Navier–Stokes operator is expected to yield highly non-orthogonal balanced modes on increasing the Reynolds number. The use of these modes as the basis functions is then not necessarily robust for the construction of a ROM, as they generally do not guarantee rapid convergence to any arbitrary functions with increasing number of degrees of freedom.
To see this, the squared values of the mode coefficients for all ROMs are shown in figure 12. Indeed, except for the BT-LNTe case, where the corresponding linear operator employed the eddy viscosity model, the ROMs with balanced truncation modes (BT-LNL and BT-LNT) have mode coefficients an order of magnitude higher than those with orthogonal modes (i.e. POD and controllability modes) especially for
$(n_x,n_z)=(1,1)$
. It is expected that these balanced-truncation-based ROMs will have even higher mode coefficients at higher Reynolds numbers, because of their increasing non-orthogonality, a potential issue in their applications to such flows. In this respect, adding an eddy viscosity model to the linearised Navier–Stokes operator appears to be useful, as it would reduce the non-orthogonal nature of the balanced truncation modes. This is well observed from the mode coefficients of the BT-LNTe case, where their values are comparable to those of the ROMs based on orthogonal modes.
Time-averaged mode coefficient fluctuations: (a)
$(n_x,n_z)=(0,1)$
; (b)
$(n_x,n_z)=(1,1)$
.
6. Concluding remarks
In the present study, we have investigated how the choice of basis functions affects the performance and convergence of Galerkin projection-based ROMs for Couette flow. To this end, POD modes from turbulent flow and controllability and balanced truncation modes of the linearised Navier–Stokes operator considering different base flows (laminar base flow and turbulent mean flow) and incorporation of an eddy viscosity model were considered. With Galerkin projection-based ROMs, these basis functions can ensure an eventual convergence with no concern in terms of overfitting, with which modern data-driven approaches may struggle. It was found that the performance and convergence of the ROMs constructed with various basis functions are significantly dependent on the state and are highly correlated with the information that each set of the basis functions contains. In the neighbourhood of the laminar base state, the balanced truncation modes from the Navier–Stokes operator linearised about the laminar base flow (BT-LNL) were found to be the best basis functions, as the resulting ROM well captures the linear stability of the laminar base flow and transient growth only with a very small number of degrees of freedom (
$N_y \simeq 5$
). The controllability modes based on the same linear operator (C-LNL) were also effective in retrieving the linear stability of the laminar state and the related transient growth. It is presumable that the good performance of these models would deteriorate if their solution trajectory reaches around the edge of turbulence (i.e. the state-space boundary between laminar and turbulent states), as it is physically associated with breakdown of streaks just before transition to turbulence, a nonlinear process not captured by the given linearised Navier–Stokes operator (e.g. Skufca et al. Reference Skufca, Yorke and Eckhardt2006). In contrast, POD modes were found to be the best basis functions for describing a turbulent state, and balanced and controllability modes based on the linearised Navier–Stokes operator employing turbulent mean flow and an eddy viscosity model (BT-LNTe and C-LNTe) were also effective for the construction of ROMs compared with other operator-based basis functions. We also expect that these models would begin to perform poorly around the edge of turbulence, where the related flow structures are significantly different from those of fully developed turbulence (e.g. Gibson et al. Reference Gibson, Halcrow and Cvitanović2008; Doohan, Willis & Hwang Reference Doohan, Willis and Hwang2019).
In a way, the conclusion here might have been expected somehow: the ROM using information from a state would converge quickly to model the dynamics around the same state. However, it is worth recalling that the degree to which each of the ROMs shows the state-dependent convergence appears to be remarkably high. For example, the cases BT-LNL and C-LNL only need a single mode in the
$y$
-direction (
$N_y=1$
) to retrieve the linear stability of the laminar base state and BT-LNL can reliably reproduce the optimal transient growth only with
$N_y=5$
. Similarly, the POD case only needs
$N_y=5$
to retrieve most of the statistical and dynamical characteristics of turbulent state, whereas all the other cases require at least
$N_y=20$
to do the same.
Nevertheless, it is certainly possible that there may be another set of basis which model the dynamics of the entire state space better than any basis considered here with even smaller
$N_y$
. Furthermore, approximation of the turbulent state here with linear subspaces may have certain limitations. However, it has been understood that even a very simple attractor, such as a limit cycle, may not necessarily be contained in any finite-dimensional subspace (Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012). Indeed, some recent data-driven approaches, which try to model a turbulent state on a nonlinear subspace (e.g. latent space created by an autoencorder), have demonstrated that neural-network-based ROMs can model the same flow with fewer than 20 degrees of freedom (e.g. Linot & Graham Reference Linot and Graham2023). However, unlike such a model, which requires an extensive set of training dataset and the related computational cost, the ROMs considered here, except the POD case, are operator-based and obtained with a minimal amount of data, such as mean flow, which may be obtained with a very cheap RANS model. Furthermore, it is robust to extrapolate the model to different Reynolds numbers, given that it is a direct approximation of the given physical system, guaranteeing convergence to the full Navier–Stokes equations. It also allows to interpret turbulent dynamics in a physically meaningful way, especially at high Reynolds numbers, due to the form identical to the Navier–Stokes equations (e.g. McCormack et al. Reference Maia and Cavalieri2024). Such benefits should not be ignored, when one considers modelling a flow with a huge number of degrees of freedom, such as turbulent flows at high Reynolds numbers.
The ultimate goal of any ROM lies in the reduction of dimensionality of the original Navier–Stokes equations. However, in practice, it is very difficult to achieve this goal without losing any of the modelling capacity of the ROM. At high Reynolds numbers, the expected dimension of any ROM would have to be non-negligibly large, as it would also increase with the Reynolds number: for example, the Lyapunov (Kaplan–Yorke) dimension of the present minimal Couette flow at a slightly lower friction Reynolds number (
) is
$D_L=14.8$
(Inubushi, Takehiro & Yamada Reference Inubushi, Takehiro and Yamada2015) (the smallest number of degrees of freedom of ROMs that capture the mean and fluctuations of turbulent state in the present study is 75 (C-LNL/LNT and POD), consistent with this number), but a similar setting in Poiseulle flow at higher Reynolds number at
$Re_\tau =80$
is
$D_L \sim 780$
(Keefe, Moin & Kim Reference Keefe, Moin and Kim1992), indicating a very rapid increase of the Lyapunov dimension with the Reynolds number. This also implies that the dimension of the ROM would also need to increase very rapidly, especially if one wishes to capture the turbulent dynamics very accurately. It is then practical and probably necessary that a ROM for high Reynolds numbers would need to make some sacrifice in its modelling capability to serve its original purpose. This rationale is similar to that of large eddy simulation (LES), where the dynamics of small scales are replaced with a simple subgrid-scale model or are not even modelled (e.g. implicit LES). The observation in the present study is useful in this regard. For example, at high Reynolds numbers, the laminar base flow is practically unseen because of its very small basin of attraction or its strongly unstable nature. Then, as long as a ROM is able to describe the turbulent state well, there would probably be no need for this ROM to capture the dynamics around the laminar base flow accurately, unless one is interested in modelling the full bifurcation scenario. From this viewpoint, the present study highlights the usefulness of POD modes from turbulent flows. These POD modes also need to be obtained from simulation and experiment data. However, obtaining a large number of converged modes requires, apart from a significant amount of temporal data, a significantly large discretisation. This can be particularly challenging for experimental data. The balanced truncation and controllability modes from the linearised Navier–Stokes equations with turbulent mean flow and an eddy viscosity model may be a suitable alternative to the POD modes (especially at high Reynolds number), despite their limitations. Developing a ROM with these modes at high Reynolds numbers is our on-going task combining with a suitable data-driven augmentation (e.g. Chua Khoo et al. 2022).
Acknowledgements
This work was initiated at the Isaac Newton Institute for Mathematical Sciences at the University of Cambridge through the programme ‘Mathematical aspects of turbulence: where do we stand?’ supported by the Engineering and Physical Sciences Research Council (EPSRC; EP/R014604/1).
Funding
I.A.M. acknowledges funding from the São Paulo Research Foundation-FAPESP through Grant No. 2022/06824-4. Y.H. gratefully acknowledges a financial support from CAPES (Brazilian Federal Agency for Support and Evaluation of Graduate Education) for the visit to Instituto Tecnológico de Aeronáutica.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Convergence of ROMs
Here, we demonstrate the convergence of C-LNTe in figure 13. This benchmark ROM contains three non-zero Fourier modes in the streamwise and spanwise directions (
$n_x = n_z = 0,\pm 1, \pm 2, \pm 3$
) and each Fourier modes contains 30 controllability modes, totalling 1470 degrees of freedom. This ROM is plotted against the C-LNTe ROM with 450 DoFs. We can notice a significant improvement in turbulent mean flow and all fluctuating quantities, where the mean flow,
$\sqrt {(u^\prime )^2}$
and
$-\overline {u^\prime v^\prime }$
, shows a nearly overlapping profile with DNS result. However, only subtle improvements are observed in terms of cross-correlations, hinting more wavenumbers in the streamwise and spanwise directions are required for C-LNTe.
Appendix B. Error measures for ROMs
In this section, we describe the formulations used in tables 3 and 4 to compute the errors between DNS and ROM quantities. For eigenvalues, the following is used:
\begin{equation} \epsilon _\lambda (n_x,n_z) = \frac {\sqrt {|\lambda _{\textit{DNS}}(n_x,n_z)-\lambda _{\textit{ROM}}(n_x,n_z)|^2}}{\sqrt {|\lambda _{\textit{DNS}}(n_x,n_z)|^2}}, \end{equation}
where
$\lambda _{\textit{DNS}}$
is the least stable eigenvalue and
$\lambda _{\textit{ROM}}$
is the eigenvalue of a ROM closest to
$\lambda _{\textit{DNS}}$
. For transient growth and harmonic response, the error is computed as follows over a given time/frequency interval:
\begin{equation} \epsilon _G (n_x,n_z)= \frac {\frac {1}{T} \int _0^{T} \sqrt {(G_{\textit{DNS}}(t;n_x,n_z)-G_{\textit{ROM}}(t;n_x,n_z))^2} \,\textrm {d}t}{\frac {1}{T} \int _0^{T} \sqrt {(G_{\textit{DNS}}(t;n_x,n_z)^2} \,\textrm {d}t}, \end{equation}
\begin{equation} \epsilon _R (n_x,n_z) = \frac {\frac {1}{\varOmega } \int _0^{\varOmega } \sqrt {(R_{\textit{DNS}}(\omega ;n_x,n_z)-R_{\textit{ROM}}(\omega ;n_x,n_z))^2}\, \textrm {d}\omega }{\frac {1}{\varOmega }\int _0^{\varOmega } \sqrt {R_{\textit{DNS}}(\omega ;n_x,n_z)^2} \,\textrm {d}\omega }. \end{equation}
For
$(n_x,n_z) = (0,1)$
,
$T=500$
and
$\varOmega =1$
, and for
$(n_x,n_z) = (1,1)$
,
$T=50$
and
$\varOmega =1.5$
. For the errors in turbulent statistics, difference integrated over the physical domain
$y$
is considered. For example, the error in the mean velocity is defined as
\begin{equation} \epsilon _{\overline {U}} = \frac {\frac {1}{2h} \int _{-y}^{y} \sqrt {({\overline {U}}_{\textit{DNS}}(y)-{\overline {U}}_{\textit{ROM}}(y))^2}\, \textrm {d}y}{\frac {1}{2h} \int _{-y}^{y}\sqrt {({\overline {U}}_{\textit{DNS}}(y))^2} }. \end{equation}
The same is applied to the r.m.s. profiles of the streamwise, wall-normal and spanwise velocities, and that of the Reynolds shear stress, denoted by
$\epsilon _{u'}$
,
$\epsilon _{v'}$
,
$\epsilon _{w'}$
and
$\epsilon _{u'v'}$
, respectively. Finally, for all the auto- and cross-correlations, the error was defined by considering a correlation time interval,
$\tau \in [-40,40]$
,
\begin{equation} \epsilon _{(m_x,m_z)}^{(n_x,n_z)}= \frac {\frac {1}{2T} \int _{-T}^{T} \sqrt {\left (\Delta C_{(m_x,m_z)}^{(n_x,n_z)}\right )^2}\, \textrm {d}\tau }{\frac {1}{2T} \int _{-T}^{T} \sqrt {\big ( C_{(m_x,m_z)}^{(n_x,n_z)}\Big |_{\textit{DNS}}\big )^2} \,\textrm {d}\tau }, \end{equation}
where
$\Delta C_{(m_x,m_z)}^{(n_x,n_z)}=C_{(m_x,m_z)}^{(n_x,n_z)}\Big |_{\textit{DNS}}-C_{(m_x,m_z)}^{(n_x,n_z)}\Big |_{\textit{ROM}}$
is difference in the correlation functions from DNS and ROM.
