Multi-scale invariant solutions in plane Couette flow: a reduced-order model approach

Abstract

Plane Couette flow at Reynolds number $Re=1200$ (based on the channel half-height and half the velocity difference between the top and bottom plates) is investigated with a spatial domain designed to retain only two spanwise integral length scales. In this system, the computation of invariant solutions that are physically representative of the turbulent state has been understood to be challenging. To address this challenge, our approach is to employ an accurate reduced-order model with 600 degrees of freedom (Cavalieri & Nogueira, Phys. Rev. Fluids, vol. 7, 2022, L102601). Using the two-scale energy budget and the temporal cross-correlation of key observables, it is first demonstrated that the model contains most of the multi-scale physical processes identified recently (Doohan et al., J. Fluid Mech., vol. 913, 2021, A8); i.e. the large- and small-scale self-sustaining processes, the energy cascade for turbulent dissipation, and an energy-cascade mediated small-scale production mechanism. Invariant solutions of the reduced-order model are subsequently computed, including 96 equilibria and 43 periodic orbits. It is found that none of the computed equilibrium solutions are able to reproduce an accurate energy balance associated with the multi-scale dynamics of the turbulent state. Incorporation of unsteadiness into invariant solutions is seen to be essential for a sensible description of the multi-scale turbulent dynamics and the related energetics, at least in this type of flow, as periodic orbits with a sufficiently long period are mainly able to describe the complex spatio-temporal dynamics associated with the known multi-scale phenomena.

1. Introduction

The theory of chaotic dynamical systems has played an increasing role in the study of turbulence in recent years, stemming from the work of Hopf (1948), who envisaged the state space of a fluid flow as that of the infinite-dimensional function space of the Navier–Stokes equations. In this view, each unique flow field is represented by a single point in the state space, with a time evolving flow tracing a continuous curve. As in traditional finite-dimensional dynamical systems theory, this trajectory through the state space is largely influenced and organised by invariant solutions such as equilibria, periodic orbits and higher-dimensional tori, which attract and repel trajectories along their stable and unstable manifolds.

This approach has seen notable success in describing the laminar–turbulent transition in parallel wall-bounded shear flows that exhibit a sub-critical transition (e.g. Kerswell 2005; Eckhardt et al. 2007; Kawahara, Uhlmann & Van Veen 2012; Graham & Floryan 2021), especially when the computational domain of interest is sufficiently small so that the concept of ‘space’ in transition can be ignored – note that when the domain size is very large, the spatial evolution of disturbances becomes important and transition has been understood to be better described by the framework of statistical mechanics, such as directed percolation (for this issue, see Barkley 2016). In this case, the transition, induced by finite-amplitude perturbations, can be viewed to be initiated by a saddle-node bifurcation, through which a pair of non-trivial exact solutions to the Navier–Stokes equations emerge. As the Reynolds number is increased, the solution with higher skin friction (upper branch) subsequently experiences a bifurcation cascade involving homoclinic tangency (e.g. Kreilos & Eckhardt 2012; Lustro et al. 2019) and evolves into a chaotic saddle (i.e. a chaotic state with finite lifetime). On the contrary, the solution with lower skin friction (lower branch) often becomes part of the laminar–turbulent separatrix, referred to as edge of turbulence, a high-dimensional surface in the state space, where a relative attractor, known as the edge state, lives (Itano & Toh 2001; Skufca, Yorke & Eckhardt 2006; Schneider, Eckhardt & Yorke 2007).

The application of dynamical systems concepts to wall-bounded shear flows has been largely restricted to the minimal flow unit of Jiménez & Moin (1991). This flow unit, periodic in both the streamwise and spanwise directions, represents the smallest near-wall flow domain that can sustain turbulence, in units based on the viscous inner length scale ($\delta _\nu = \nu /u_\tau$, where $\nu$ is the kinematic viscosity and $u_\tau$ the friction velocity), hereafter referred to with a ${\cdot }^+$ superscript. This box size is directly related to the spanwise spacing of streaks ($\lambda ^+_z \approx 100$) found in the near-wall region of turbulent flows, first identified experimentally by Kline et al. (1967). Therefore, by construction, the minimal flow unit isolates a single pair of low- and high-speed streaks and captures the related minimal dynamics required to sustain turbulence, often referred to as the self-sustaining process (Hamilton, Kim & Waleffe 1995).

This represents a quasi-cyclic process whereby streaks are first generated by streamwise vortices that take energy from the mean shear through the lift-up effect (Moffatt 1965; Ellingsen & Palm 1975; Landahl 1980, 1990; Butler & Farrell 1993; del Alamo & Jimenez 2006; Cossu, Pujals & Depardon 2009; Pujals et al. 2009; Hwang & Cossu 2010a; McKeon & Sharma 2010). The streaks then undergo a normal-mode instability and/or transient growth, forming downstream undulation of streaks, which subsequently promote the generation of streamwise vortices as a result of nonlinear self-interactions (Hamilton et al. 1995; Schoppa & Hussain 2002; Park, Hwang & Cossu 2011; Hwang & Bengana 2016; Cassinelli, de Giovanetti & Hwang 2017; Lozano-Durán et al. 2021).

From a dynamical systems perspective, the self-sustaining process has firmly been associated with a variety of non-trivial invariant solutions that have been computed in plane Couette flow (Nagata 1990; Clever & Busse 1997; Kawahara & Kida 2001; Viswanath 2007; Gibson, Halcrow & Cvitanović 2009; Cvitanović & Gibson 2010; Hall & Sherwin 2010; Deguchi, Hall & Walton 2013), as well as in other canonical shear flows such as pipe flow (Faisst & Eckhardt 2003; Wedin & Kerswell 2004; Willis, Short & Cvitanović 2016), plane Poiseuille flow (Waleffe 1998, 2001; Park & Graham 2015; Hwang, Willis & Cossu 2016; Yang, Willis & Hwang 2019), stress-driven shear flows (Doohan, Willis & Hwang 2019; Doohan et al. 2022) and boundary layer flow (Khapko et al. 2013; Kreilos et al. 2013; Deguchi & Hall 2014). In the minimal flow unit, these invariant solutions have been shown to form a skeleton for the chaotic dynamics of the flows (Gibson, Halcrow & Cvitanović 2008; Willis et al. 2016; Doohan et al. 2019), and are typically seen to capture the underlying mathematical structure and dynamics of the coherent structures. These solutions typically have both stable and unstable manifolds which attract and repel chaotic trajectories throughout the state space, with turbulence being viewed as a chaotic walk between these solutions (Kawahara et al. 2012; Graham & Floryan 2021). In this sense, periodic orbits have been understood to provide a low-dimensional representation of a smooth dynamical system, and suitable averaging over these solutions has the ability to yield a meaningful statistical description of the chaotic state (Artuso, Aurell & Cvitanovic 1990; Chandler & Kerswell 2013; Yalnız, Hof & Budanur 2021; Page et al. 2022).

As the Reynolds number is increased, it has been shown that the self-sustaining process possibly exists at each of the integral length scales varying from the viscous inner to the large-scale outer units in wall-bounded turbulence (Flores & Jiménez 2010; Hwang & Cossu 2010b, 2011; Hwang 2015; Hwang & Bengana 2016). At each integral length scale, the self-sustaining process supports energy-containing coherent structures, composed of streaks and streamwise vortices, and they are statistically and dynamically self-similar with respect to the integral length scale, consistent with the attached eddy hypothesis of Townsend (1956). There have been efforts to characterise the possible interactions between these self-similar coherent structures and the related turbulent energy transfer (Cho, Hwang & Choi 2018; Lee & Moser 2019), but the dynamical processes underpinning these scale interactions currently remain poorly understood.

The simplest system containing the dynamics of these interactions was recently studied by Doohan, Willis & Hwang (2021), who considered a computational domain twice greater than the minimal flow unit, referred to as the minimal multi-scale flow unit. By construction, this domain only admits the coherent structures at two integral length scales (i.e. $\lambda _z^+\simeq 100, 200$). Several important multi-scale processes have been numerically identified: (i) large- and small-scale self-sustaining processes; (ii) energy cascade via the streak instability (or transient growth) and breakdown; (iii) energy transport from the large to small scales which drives small-scale turbulent production; (iv) the feeding of energy from small to large scales (Doohan et al. 2021). Furthermore, a number of equilibria and two periodic orbits in this system were computed (Doohan et al. 2022). However, none of them were representative of the full system's chaotic dynamics, as had been seen previously in the minimal flow unit – most of them were found to contain the single-scale dynamics only associated with the corresponding self-sustaining process, and they did not appear to be key players forming the state space structure of turbulence, since they lacked some of the other physically relevant processes.

The observations of Doohan et al. (2022) raise some questions regarding the structure of the state space of turbulence at high Reynolds numbers. It has often been believed that the closure of a set of relevant invariant solutions of the Navier–Stokes equations form the ergodic subspace of the turbulent state, in the sense that the solution trajectory visits or shadows many of these solutions, at least for some amount of time. This is an important requirement for a meaningful statistical description of the turbulent state in terms of the invariant solutions. However, as the Reynolds number is increased, many of these solutions, especially those which are relatively easily accessible through the numerical continuation of their counterpart computed at lower Reynolds numbers, are no more directly relevant to turbulence, as the solution trajectory appears to very rarely visit them. This now emphasises the necessity of capturing the invariant solutions buried around the turbulent state at high Reynolds numbers. Such solutions are expected to contain the full multi-scale features of the turbulent state like the one recently obtained by Motoki, Kawahara & Shimizu (2021) for thermal convection.

Computation of multi-scale invariant solutions organised around the turbulent state, however, appears to be increasingly challenging as the Reynolds number is increased. This is directly connected to the extremely high computational expense of finding these solutions, as the computation typically involves a large number of iterative direct numerical simulations (DNS) combined with a suitable searching algorithm, such as the Newton–Krylov method (Viswanath 2007; Willis, Cvitanović & Avila 2013) and adjoint method (Farazmand 2016). Furthermore, any recurrence measure, monitored to find a suitable initial condition for the searching algorithm, becomes increasingly large with increased Reynolds number, as turbulence contains increasingly many different scales. Therefore, a few new approaches to bypass this difficulty have been proposed (Page & Kerswell 2020; Azimi, Ashtari & Schneider 2022; Page et al. 2022; Parker & Schneider 2022). Finally, the leading Lyapunov exponent of turbulence grows faster than the inverse of the Kolmogorov time scale (Mohan, Fitzsimmons & Moser 2017), indicating that finding good initial guesses for a periodic orbit becomes increasingly difficult at higher Reynolds numbers, at least with some of the most popular techniques currently available (Viswanath 2007; Willis et al. 2013).

In earlier studies, these difficulties have been bypassed by modelling some physical processes of turbulence instead of using full DNS. A popular approach taken was to model the small-scale eddies around the large-scale coherent structures of interest using an eddy viscosity implemented by an over-damped large-eddy simulation (Rawat et al. 2015; Hwang et al. 2016; van Veen, Vela-Martín & Kawahara 2019; Yang et al. 2019; Azimi, Cossu & Schneider 2020; Sasaki, Kawahara & Jiménez 2021). However, this approach eliminates the complex nonlinear interactions between the large-scale coherent structures and background turbulence, often offering only limited information on the physical process involved (e.g. self-sustaining process). To this end, in the present study, a step forward from the early approach will be taken. In particular, the objective of the present study is to explore the state space of a reduced-order model (ROM) for the minimal multi-scale (i.e. two-scale) flow unit in plane Couette flow with the computation of its invariant solutions. We will see that the ROM-based approach enables us to delicately control the degrees of freedom of the flow to be sufficiently low (only 600 degrees of freedom), while retaining all the multi-scale physical processes of interest. This significantly relieves the difficulties in the computation of invariant solutions, allowing us to compute a large number of equilibria and unstable periodic orbits.

The ROM used in the present study was recently proposed by Cavalieri & Nogueira (2022), who demonstrated that it exhibits long-time numerical stability with reasonable statistical agreement to reference DNS. A classical approach to formulating a ROM is based on the Galerkin projection of the proper orthogonal decomposition (POD) modes, first introduced in the context of turbulence by Lumley (1967). In this approach, the POD basis functions, representing coherent structures, are ranked by their turbulent kinetic energy, allowing the highest energy-containing structures to be included in the model. The approach has been used extensively in a wide variety of applications in fluid dynamics using both experimental and numerical data (e.g. Aubry et al. 1988; Noack et al. 2003; Smith, Moehlis & Holmes 2005; Khoo, Chan & Hwang 2022 and many others). Such an approach, however, requires a priori DNS data to construct the model, and thus it can be computationally demanding. A benefit of the approach by Cavalieri & Nogueira (2022) is that the orthogonal basis functions considered are the POD modes of the stochastic response of the linearised Navier–Stokes equations, as such basis functions can be obtained a priori by solving the related Lyapunov equation (Farrell & Ioannou 1993; Bamieh & Dahleh 1999; Jovanović & Bamieh 2005) at minimal computational cost.

The ROM of Cavalieri & Nogueira (2022) is designed for plane Couette flow, but the multi-scale nonlinear processes of interest in this study were discovered from the shear stress-driven flow model of Doohan et al. (2021). Given that the only difference between the two systems is the top boundary condition, we expect that the multi-scale processes of interest will be well captured by the ROM of Cavalieri & Nogueira (2022) – indeed, we shall see that this ROM retains all the multi-scale physical processes discovered by Doohan et al. (2021), thereby making the dynamics of its invariant solutions physically meaningful. Furthermore, plane Couette flow is a flow configuration where a large number of equilibria and periodic orbits previously computed are well documented (e.g. Gibson et al. 2009; Cvitanović & Gibson 2010). As will be evident (see also § 5), the ROM of Cavalieri & Nogueira (2022) also has the potential to be used for the improvement of existing search methods for equilibria and periodic orbits. From this perspective, retaining the plane Couette flow setting in the ROM would also be beneficial.

This paper is organised as follows: § 2 summarises the formulation of the ROM, § 3 discusses its validation with respect to a time-averaged energy budget as well as the temporal dynamics of spatially averaged relevant observables, § 4 discusses the computation of the invariant solutions and the corresponding results, followed by a conclusion in § 5. Additional details about the terms in the energy budget equations are contained in Appendix A, and further results discussing the validity of the ROM are contained in the supplementary material available at https://doi.org/10.1017/jfm.2024.108.

2. Reduced-order model

2.1. Model formulation

Incompressible plane Couette flow is considered, where two infinitely wide and long parallel plates, separated by a wall-normal distance $2h$, are set to move in opposite directions with velocity $\pm U_0$. All variables are made dimensionless with length and velocity scales $h$ and $U_0$, respectively. The dimensionless flow domain of interest is considered to be $\boldsymbol {x}(=(x,y,z)) \in [0,2{\rm \pi} ] \times [-1,1] \times [0,{\rm \pi} ]$, equipped with periodic boundary conditions in the streamwise $(x)$ and spanwise $(z)$ directions and the velocity $\tilde {\boldsymbol {u}}|_{y=\pm 1} = (\pm 1,0,0)$ specified on the upper and lower wall-normal $(y)$ boundaries. The velocity field is first decomposed into the laminar state, $\boldsymbol {u}_0=(y,0,0)$, and the fluctuations about this state, $\boldsymbol {u}'$

(2.1)\begin{equation} \tilde{\boldsymbol{u}}(\boldsymbol{x},t) = \boldsymbol{u}_0(y) + \boldsymbol{u}'(\boldsymbol{x},t). \end{equation}

The evolution of this fluctuating velocity field is described by the following equations:

(2.2)\begin{gather} \partial_t \boldsymbol{u}' + (\boldsymbol{u}_0\boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}' + (\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}_0 + (\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}'={-}\boldsymbol{\nabla} p + \frac{1}{Re}\nabla^2\boldsymbol{u}', \end{gather}

(2.3)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}' = 0, \end{gather}

with boundary conditions given by $\boldsymbol {u}'|_{y=\pm 1} = (0,0,0)$, where $Re=U_0 h/ \nu$ is the Reynolds number, and $\nu$ is the kinematic viscosity.

Following the formulation of Cavalieri & Nogueira (2022), we define an inner product

(2.4)\begin{equation} \langle \,\boldsymbol{f},\boldsymbol{g} \rangle = \frac{1}{4{\rm \pi}^2} \int_0^{\rm \pi} \int_{{-}1}^1 \int_0^{2{\rm \pi}} \boldsymbol{f}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{g}(\boldsymbol{x}) \,\mathrm{d}\kern0.06em x \, \mathrm{d}\kern0.05em y \, \mathrm{d}z, \end{equation}

and introduce a modal decomposition of the fluctuating velocity field

(2.5)\begin{equation} \boldsymbol{u}'(\boldsymbol{x},t) = \sum_{j=1}^N a_j(t) \boldsymbol{\varPhi}_j(\boldsymbol{x}), \end{equation}

where the modes $\boldsymbol {\varPhi }_j$ are defined to be divergence free $\boldsymbol {\nabla } \boldsymbol {{\cdot }}\boldsymbol {\varPhi }_j =0$, and orthonormal with respect to the defined inner product, i.e. $\langle \boldsymbol {\varPhi }_i,\boldsymbol {\varPhi }_j \rangle = \delta _{ij}$, where $\delta _{ij}$ is the Kronecker delta. Hence, by construction, any linear superposition of modes satisfies both the incompressibility condition and the boundary conditions of the system.

Substituting the decomposition (2.5) into (2.2) and taking an inner product with the mode $\boldsymbol {\varPhi }_i$ leads to a system of ordinary differential equations describing the evolution of the mode coefficients $a_i(t)$

(2.6)\begin{equation} \frac{{\rm d}a_i}{{\rm d}t} = \frac{1}{Re} \sum_{j=1}^N L_{ij} a_j + \sum_{j=1}^N \tilde{L}_{ij} a_j + \sum_{j=1}^N \sum_{k=1}^N Q_{ijk} a_j a_k, \end{equation}

where

(2.7a)\begin{gather} L_{ij} = \langle \nabla^2 \boldsymbol{\varPhi}_j,\boldsymbol{\varPhi}_i\rangle, \end{gather}

(2.7b)\begin{gather}\tilde{L}_{ij} ={-}\langle (\boldsymbol{\varPhi}_j \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_0 + \boldsymbol{u}_0\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{\varPhi}_j), \boldsymbol{\varPhi}_i \rangle, \end{gather}

(2.7c)\begin{gather}Q_{ijk} ={-}\langle (\boldsymbol{\varPhi}_j \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\varPhi}_k), \boldsymbol{\varPhi}_i \rangle. \end{gather}

Here, $L_{ij}$ and $\tilde {L}_{ij}$ represent the viscous term and the linear interaction with the laminar state, respectively, whereas $Q_{ijk}$ represents the quadratic nonlinear interaction between the different modes. The pressure dependence is seen to vanish due to the divergence-free property of the modes.

This system of ordinary differential equations still describes the exact dynamics in the limit as $N\rightarrow \infty$ by the orthogonality of the modes, although here we attempt to choose the modes in an optimal way such that we can model the dynamics of the system using a greatly reduced number of degrees of freedom.

Many choices exist for the modes $\boldsymbol {\varPhi }_j(\boldsymbol {x})$. In this study, POD modes of the stochastic response of the linearised Navier–Stokes equations are chosen, and the detailed computational procedure is described in Cavalieri & Nogueira (2022), who followed the method of Jovanović & Bamieh (2005). The velocity covariance resulting from the stochastic response is obtained by formulating the controllability Gramian of the Navier–Stokes equations linearised about the laminar basic state, $\boldsymbol {u}_0$, and is computed by solving the associated Lyapunov equation. The POD modes are obtained by computing the eigenfunctions of the velocity covariance that are ranked in terms of their energy (i.e. eigenvalue).

Compared with a typical POD approach, the use of POD modes from the stochastic response of the linearised Navier–Stokes equations removes the necessity to form these modes using large numerical or experimental datasets, as they can be directly computed from the forced linear system. The results of Cavalieri & Nogueira (2022) suggest that these modes provide an efficient set of basis functions for the construction of the ROM, as the resulting ROM provides equal or more accurate results than a Galerkin–POD approach constructed from DNS data, at least in terms of the second-order velocity statistics.

2.2. Numerical implementation

To implement this model numerically, the POD modes from the controllability Gramian are represented using a Fourier discretisation in the streamwise and spanwise directions, and a Chebyshev discretisation in the wall-normal direction

(2.8)\begin{equation} \boldsymbol{\varPhi}_j(\boldsymbol{x}) = \hat{\boldsymbol{\varPhi}}_j(y)\,{\rm e}^{\mathrm{i}(m k_{x,0}x + n k_{z,0}z)}, \end{equation}

where $\widehat {{\cdot }}$ denotes the Fourier transform, and $k_{x,0} = {2{\rm \pi} }/{L_x} = 1$ and $k_{z,0} = {2{\rm \pi} }/{L_z} = 2$ denote the fundamental wavenumbers in the streamwise and spanwise directions, respectively. In particular, we use the modes constructed in Cavalieri & Nogueira (2022) with the number of modes $N=600$: $m=\{0,1,2\}, n=\{-2,-1,0,1,2\}$ with $N_{y}=24$, where $N_{y}$ is the number of POD modes used in the expansion for a given $m$ and $n$. Here, we note that only positive streamwise wavenumbers need be considered since negative wavenumber modes can be reconstructed by complex conjugation. The number of degrees of freedom here was chosen by running a preliminary test, and $N=600$ was found to be close to the smallest possible dimension for a system retaining self-sustaining processes at two scales. Two wavenumbers must be considered for each of the two wall-parallel directions to retain the ‘two’ scales, and $N_y=24$ was chosen to ensure the linear stability of the laminar solution. A bifurcation analysis of the ROM was subsequently carried out, revealing that its behaviour is remarkably similar to that of plane Couette flow; i.e. saddle-node bifurcation and the resulting bifurcation cascade leading to turbulence (e.g. Kreilos & Eckhardt 2012). These modes can then be used to construct the system of ordinary differential equations (2.6) which are integrated in time using a standard fourth- or fifth-order Runge–Kutta scheme at $Re = 1200$ using the ode45 function in MATLAB.

The time evolution of this system starting from a random initial condition displays chaotic behaviour consistent with turbulence, as shown by Cavalieri & Nogueira (2022). The long trajectory was initialised from a random perturbation of the mode coefficients $a_i$ with the amplitude between 0 and 0.1, and was time integrated using ode45 with a relative error tolerance of $10^{-8}$. Temporal averaging of the flow snapshots from the model over a sufficiently long time leads to the mean and root mean square (r.m.s.) velocity profiles shown in figure 1. These profiles have been compared with reference DNS statistics at a grid resolution of $(N_x, N_y, N_z)=(64, 65, 64)$ computed using the open-source code Channelflow (Gibson 2014, http://channelflow.org/), showing reasonable agreement relative to the greatly reduced computational cost. Note that the number of degrees of freedom of the DNS is $N_{DNS} \approx 10^6$. Computation of the velocity gradient at the wall leads to a friction Reynolds number $Re_\tau \approx 80$, giving a domain size of $(L_x^+, L_y^+, L_z^+) = (504,160,250)$ based on inner units. These dimensions are bigger than the minimal flow unit configurations of Hamilton et al. (1995) and Kawahara & Kida (2001), and are of similar size to that of the minimal multi-scale flow unit formulated by Doohan et al. (2019). In this regard, we expect similar two-scale behaviour as observed by Doohan et al. (2019) to be captured by the ROM.

Figure 1. Mean and r.m.s. velocity profiles: (a) $U^+$; (b) $u_{rms}^+$; (c) $v_{rms}^+$; (d) $w_{rms}^+$. The superscript $({\cdot })^+$ denotes the viscous inner-scaled variables. Here, the black and blue lines indicate those from the ROM and DNS, respectively.

3. Energy budget and dynamics

To ensure that any invariant solutions found for the ROM are physically relevant to the corresponding exact flow, we first examine the two-scale time-averaged energy budget of the flow as a function of the wall-normal direction. The results of this energy budget can then be qualitatively compared with the results of the energy budget analysis of the similar two-scale system studied by Doohan et al. (2021).

3.1. Two-scale energy budget

To formulate a relevant energy budget for a statistically stationary turbulent state, the full velocity field is decomposed as

(3.1)\begin{equation} \tilde{\boldsymbol{u}}(\boldsymbol{x},t) = \boldsymbol{U}(y) + \boldsymbol{u}(\boldsymbol{x},t), \end{equation}

where $\boldsymbol {U}$ describes the mean flow, defined as $\boldsymbol {U}(\boldsymbol {x}) = ({1}/{T})\int _0^T \tilde {\boldsymbol {u}}(\boldsymbol {x},t)\,\mathrm {d}t$ for some suitably long time interval $T$, and $\boldsymbol {u}(\boldsymbol {x},t)$ is the velocity fluctuation about the mean flow.

As discussed in Doohan et al. (2019), the inner-scaled flow domain at $Re_\tau (\simeq 80)$ considered here is approximately $L_z^+\simeq 250$. Given that the smallest spanwise size of energy-containing structures is $\lambda _z^+\simeq 100$, the computational domain would admit energy-containing structures mainly at two different spanwise length scales: $\lambda _z^+\simeq 125, 250$. To understand the interactions of the two different scales within the flow, we further decompose the fluctuating velocity field into a large- and small-scale fluctuating velocity field

(3.2)\begin{equation} \boldsymbol{{u}}(\boldsymbol{x},t) = \boldsymbol{{u}}_l(\boldsymbol{x},t) + \boldsymbol{u}_s(\boldsymbol{x},t), \end{equation}

where the large- and small-scale velocity fields are defined through the large- and small-scale projection operators, $\mathcal {P}_l\{{\cdot }\}$ and $\mathcal {P}_s\{{\cdot }\}$, which split the fields by spanwise wavenumber in the following way:

(3.3a)\begin{gather} \boldsymbol{u}_l(\boldsymbol{x},t) = \mathcal{P}_l\{\boldsymbol{u}(\boldsymbol{x},t)\} = \sum_{|m|\le m_x} \sum_{|n|\le 1} \hat{\boldsymbol{u}}(y)\, {\rm e}^{\mathrm{i}(m k_{x,0}x + n k_{z,0}\,z)}, \end{gather}

(3.3b)\begin{gather}\boldsymbol{u}_s(\boldsymbol{x},t) = \mathcal{P}_s\{\boldsymbol{u}(\boldsymbol{x},t)\} =\sum_{|m|\le m_x} \sum_{2\le|n|\le n_z} \hat{\boldsymbol{u}}(y)\, {\rm e}^{\mathrm{i}(m k_{x,0} x + n k_{z,0}\,z)}. \end{gather}

The split by spanwise wavenumber here is justified by the observation that the statistical structure of eddies in wall-bounded flows is self-similar with respect to the spanwise wavelength and proportional to the distance from the wall (Hwang 2015). Additionally, both large and small scales are expected to contain a range of streamwise wavenumbers due to the presence of both long streaks and short streamwise vortices as part of the self-sustaining process (Hwang 2015).

The evolution equations for turbulent fluctuations are subsequently projected onto the large- and small-scale velocity subspaces using the definition of $\mathcal {P}_l\{{\cdot }\}$ and $\mathcal {P}_s\{{\cdot }\}$, yielding two fluctuation equations, one equation each for the large and small scales, respectively. Multiplying each of these equations by the corresponding large- or small-scale velocity and averaging in the streamwise and spanwise direction yields energy budget equations for each component of the velocity at large

(3.4a)\begin{gather} \frac{\partial E_{ul}}{\partial t} = {P}_{ul} + T_{ul} + \varPi_{ul} + T_{\nu,ul} + \varepsilon_{ul}, \end{gather}

(3.4b)\begin{gather}\frac{\partial E_{vl}}{\partial t} = T_{vl} + \varPi_{vl} + T_{\nu,vl} + T_{p,vl} + \varepsilon_{vl}, \end{gather}

(3.4c)\begin{gather}\frac{\partial E_{wl}}{\partial t} = T_{wl} + \varPi_{wl} + T_{\nu,wl} + \varepsilon_{wl}, \end{gather}

and small scales

(3.5a)\begin{gather} \frac{\partial E_{us}}{\partial t} = {P}_{us} + T_{us} + \varPi_{us} + T_{\nu,us} + \varepsilon_{us}, \end{gather}

(3.5b)\begin{gather}\frac{\partial E_{vs}}{\partial t} = T_{vs} + \varPi_{vs} + T_{\nu,vs} + T_{p,vs} + \varepsilon_{vs}, \end{gather}

(3.5c)\begin{gather}\frac{\partial E_{ws}}{\partial t} = T_{ws} + \varPi_{ws} + T_{\nu,ws} + \varepsilon_{ws}, \end{gather}

where $P$ denotes the energy production, $T$ the turbulent transport, $\varPi$ the pressure strain, $T_\nu$ the viscous transport, $T_p$ the pressure transport and $\varepsilon$ the dissipation at each scale and component, where the subscripts $l$ and $s$ denote large and small scales, and $u$, $v$ and $w$ the streamwise, wall-normal and spanwise components. Full details of each term are included in Appendix A. Since only the velocity field is computed directly by the ROM, the fluctuating pressure field must be computed using the Poisson equation for pressure projected onto the large and small scales

(3.6a)\begin{gather} \nabla^2 p_l ={-} U\frac{\partial^2 u_l}{\partial x^2} - \frac{\partial v_l}{\partial x} \frac{\partial U}{\partial y} - \mathcal{P}_l\{\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u})\}, \end{gather}

(3.6b)\begin{gather}\nabla^2 p_s ={-} U\frac{\partial^2 u_s}{\partial x^2} - \frac{\partial v_s}{\partial x} \frac{\partial U}{\partial y} - \mathcal{P}_s\{\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u})\}, \end{gather}

with boundary conditions specified as $\partial _y p_l |_{y=\pm 1}= 0$, and $\partial _y p_s |_{y=\pm 1}= 0$ (Kim 1989).

The various multi-scale energy budget terms have been computed from the ROM. The temporal averaging has been taken over a dataset which is roughly 18 000 dimensionless time units in length $(T^+ = T u_\tau ^2 /\nu \approx 96\,500)$. They have been plotted over a channel half-height as a function of the wall-normal direction in figure 2. The overall statistical behaviour of the two-scale energy budget terms are qualitatively very similar to those found by Doohan et al. (2021), who used the DNS results of their two-scale shear stress-driven flow model.

Figure 2. Wall-normal variation of the time-averaged terms of the energy budget equations (3.4) and (3.5) for (a,b) streamwise, (c,d) wall-normal and (e,f) spanwise components: (a,c,e) large scale; (b,d,f) small scale.

A large amount of energy production is seen in both the large and small scales (blue lines in figure 2a,b), featuring peaks in production in the near-wall region at approximately $y^+\approx 10$, the location of which is consistent with that reported by Bech et al. (1995) and Afzal, Seena & Bushra (2018) as well as many others. The magnitude of the peak of the small-scale production is seen to be larger than that of the large scales, similar to the results reported by Doohan et al. (2021), although slightly less production is observed overall in the logarithmic region of the flow compared with this work.

The pressure strain terms (yellow lines in figure 2) sum to zero at each scale as a consequence of continuity, i.e. $\varPi _{ul}+\varPi _{vl}+\varPi _{wl}=0$ and $\varPi _{us}+\varPi _{vs}+\varPi _{ws}=0$. The streamwise components of the pressure strain are negative across the whole wall-normal domain, whereas the wall-normal and spanwise pressure strains are seen to be positive, except at $y^+ \lesssim 20$. This is indicative of the pressure strain terms being responsible for the redistribution of the turbulent kinetic energy from the streamwise to the wall-normal and spanwise components, also observed in Doohan et al. (2021). However, the pressure strain terms have a relative magnitude lower than that of Doohan et al. (2021), presumably due to the low number of streamwise and spanwise Fourier modes used for the present ROM. For $y^+ \lesssim 20$, large peaks in the spanwise and wall-normal components are observed. The near-wall negative wall-normal pressure strain is attributed to the well-known ‘splat effect’ (Mansour, Kim & Moin 1988; Perot & Moin 1995; Lee & Moser 2019), whereby fluid motions moving towards the wall must be turned to move parallel to the wall due to its impermeability.

The streamwise turbulent transport term at the large scales is seen to be positive, especially in the near-wall region (green line in figure 2$a$), implying that the turbulent transport from the small scales injects energy into the large scales. This is the inverse energy transfer (or energy feeding from the small scales), previously observed (Cho et al. 2018; Doohan et al. 2021), and it gradually disappears on increasing $y^+$. In general, the turbulent transport terms follow similar trends to those of Doohan et al. (2021) but are more oscillatory (green lines in figure 2). This lack of smoothness is likely a result of the nonlinear interactions between the truncated modes not being fully captured by the model, although the consequences of this are not immediately clear. For this reason, the turbulent transport terms have been further decomposed and studied in § 3.2 as their role is of significant interest in relation to the accuracy of the multi-scale behaviour in ROM.

The largest amounts of dissipation are seen in the near-wall region through the streamwise component (orange line in figure 2$a$), similar to Doohan et al. (2021). There is slightly greater large-scale dissipation than small-scale dissipation for all components (orange lines in figure 2), suggesting that the cascade of energy from large to small scales (Kolmogorov 1941) is not fully resolved. This is likely due to the truncated nature of the model, whereby energy cannot be transferred to the smaller scales neglected by the model, as a result of spectral blocking. Nevertheless, the amounts of dissipation at each scale are still comparable due to the low $Re$ considered here, as in Doohan et al. (2021). However, at higher Reynolds numbers, the energy cascade is expected to be more poorly resolved due to the increasing number of small-scale eddies for turbulence dissipation. As such, the number of degrees of freedom of the ROM would need to be increased further.

Lastly, the pressure and viscous transport terms appropriately balance the other terms, and their behaviours are also consistent with those found in Doohan et al. (2019). The integrated energy is seen to be balanced over all terms to within 2 % relative to the magnitude of the turbulent production.

3.2. Inter- and intra-scale energy transport

To access the role of the turbulent transport in greater detail, we decompose these terms further, following the methods of Doohan et al. (2021) and Kawata & Alfredsson (2019)

(3.7a)\begin{gather} T_{ul} = T_{ul,-} +T_{ul,\#} -T_{u,\updownarrow}, \end{gather}

(3.7b)\begin{gather}T_{vl} = T_{vl,-} +T_{vl,\#} -T_{v,\updownarrow}, \end{gather}

(3.7c)\begin{gather}T_{wl} = T_{wl,-} +T_{wl,\#} -T_{w,\updownarrow}, \end{gather}

(3.7d)\begin{gather}T_{us} = T_{us,-} +T_{us,\#} +T_{u,\updownarrow}, \end{gather}

(3.7e)\begin{gather}T_{vs} = T_{vs,-} +T_{vs,\#} +T_{v,\updownarrow}, \end{gather}

(3.7f)\begin{gather}T_{ws} = T_{ws,-} +T_{ws,\#} +T_{w,\updownarrow}. \end{gather}

Here, $T_-$ denotes the intra-scale spatial turbulent transport, $T_\#$ the inter-scale spatial turbulent transport and $T_\updownarrow$ the inter-scale turbulent transport. Here, the use of intra refers to the involvement of velocity fields of only the same scale, whereas inter refers to terms that include the nonlinear coupling of velocity fields of different scales. Spatial transport refers to terms that can be rewritten as the divergence of a vector field, meaning that they do not contribute to the loss or gain of energy when integrated over the spatial domain. On the contrary, the inter-scale turbulent transport cannot be written as the divergence of a vector field and is thus directly responsible for energy exchanges between the large and small scales. Importantly, these inter-scale turbulent transport terms appear in both the large- and small-scale energy transport equations with opposing signs, thus cancelling each other out when considering the total energy budget of both the large and small scales. The full details of these equations are included in Appendix A.

The time-averaged energy transport terms in (3.7) are plotted as a function of the wall-normal direction in figure 3. The inter-scale turbulent transport terms, $T_{\updownarrow }$, are responsible for the direct transport of energy between scales. In the streamwise component, only inter-scale turbulent transport, $T_{\updownarrow }$, from the small to large scales is observed across the entire domain in the mean (figure 3a,b). This is contrary to the DNS results of Doohan et al. (2021), which show a contribution to the energy cascade in the outer region of the flow in the streamwise direction. With that being said, the spanwise and wall-normal components are seen to contribute to the energy cascade from large to small scale in a similar way to that of Doohan et al. (2021), suggesting that the cascade is still being captured, at least to some degree, in the outer region of the flow (figure 3c–f). The streamwise and spanwise components of the inter-scale turbulent transport, $T_{\updownarrow }$, show a large amount of energy transport from the small to large scales in the near-wall region (figure 3). This is representative of the near-wall energy feeding process from small to large scale, consistent with both Cho et al. (2018) and Doohan et al. (2021). Note that this inverse energy transfer from small to large scale (e.g. Cho et al. 2018; Andreolli, Quadrio & Gatti 2021) has been known to break the self-similarity of the energy-containing motions in the region close to the wall, where the attached eddy hypothesis is not supposed to be valid (for a further discussion, see Hwang 2016).

Figure 3. Wall-normal variation of the time-averaged intra- and inter-scale energy transport terms in (3.7) for (a,b) streamwise, (c,d) wall-normal and (e,f) spanwise components: (a,c,e) large scale; (b,d,f) small scale.

The inter-scale spatial transport, $T_{\#}$, describes an energy transport that results from interactions between the different scales. The streamwise inter-scale spatial transport terms show local maxima at a wall-normal location similar to that of the peaks of the inter-scale turbulent transport and production, suggesting that these transport terms also play a role in transporting the turbulent kinetic energy produced by the flow between the scales. However, all of these inter-scale spatial transport terms exhibit similar trends at both the large and small scales (at least with respect to their signs) to those in Doohan et al. (2021), suggesting that these terms act in a similar manner at both scales.

The intra-scale spatial transport, $T_{-}$, represents energy transport within each scale. At the large scales, these terms follow reasonably similar trends to Doohan et al. (2021). However, no intra-scale transport is exhibited by the ROM at the small scales, a feature that is not in line with the DNS results of Doohan et al. (2021). This is a consequence of the spanwise averaging over the small-scale intra-scale transport terms equalling zero. For example, the streamwise intra-scale spatial turbulent transport is given by

(3.8)\begin{equation} T_{us,-} ={-}\langle u_s (\boldsymbol{u}_s\boldsymbol{\cdot} \boldsymbol{\nabla} u_s) \rangle_{x,z}. \end{equation}

We note that $T_{us,-}$ takes the form $\langle\, f(x,y,n,t) \,{\rm e}^{\mathrm {i}2(n_1+n_2+n_3)z}\rangle _{x,z}$, where $n_i$ refers to the spanwise wavenumber related to the $i$th term in the triple correlation of (3.8). Since $n_1+n_2+n_3 \neq 0$ at the small scales on account of the model's truncation (i.e. $n_i=\pm 2$ from (3.3)), this spanwise average will always be zero at the small scales unless additional wavenumbers are considered in the present ROM.

The immediate effect of this on the results of the model is not apparent, but with that being said, intra-scale processes that are responsible for the transport of energy from the production to the dissipation are seen to persist, with the pressure strain being largely responsible for this redistribution. Additionally, the viscous transport, $T_\nu$, which itself only involves intra-scale terms, is seen to be active in relation to the production and dissipation terms at both scales.

3.3. Temporal dynamics

Having studied the time-averaged variation of the different terms in the energy equations, it is also of interest to study the temporal relationship between these quantities and to confirm the presence of known dynamic processes found in the similar system of Doohan et al. (2021). These include the large- and small-scale self-sustaining process, the energy cascade, the driving mechanism which transports energy from large-scale processes and injects it into the small scales through the small-scale production and the feeding mechanism which is responsible for transporting energy from small to large scales. To study the temporal correlation of such quantities, we introduce the temporal cross-correlation of two time-varying functions $f(t)$ and $g(t)$

(3.9)\begin{equation} \{\,f\star g\}(\tau) = \frac{\displaystyle\int_{\Delta T} f(t+\tau)g(t)\,\mathrm{d}t}{\left( \displaystyle\int_{\Delta T} f(t)^2\,\mathrm{d}t \right) ^{1/2} \left( \displaystyle\int_{\Delta T} g(t)^2\,\mathrm{d}t \right) ^{1/2}}, \end{equation}

where $\Delta T$ is a suitable time interval of interest.

Of particular interest is the presence of the large- and small-scale self-sustaining processes (Hwang & Bengana 2016; Doohan et al. 2021). With this in mind, the energy related to the different velocity components is separated into streamwise-dependent and -independent parts (Lucas & Kerswell 2017; Doohan et al. 2019, 2021). Specifically, these provide observables based on the kinetic energy for the structural components of the self-sustaining process; i.e. straight streaks $(E_{ss})$, wavy streaks $(E_{ws})$, straight rolls $(E_{sr})$ and wavy rolls $(E_{wr})$ at both the large

(3.10a)\begin{gather} E_{ss,l}(t) = \tfrac{1}{2} \langle \langle u_l \rangle_x^2 \rangle _z \big|_0^{80} , \end{gather}

(3.10b)\begin{gather}E_{ws,l}(t) = \tfrac{1}{2} \langle (u_l - \langle u_l \rangle_x)^2 \rangle _{x,z} |_0^{80}, \end{gather}

(3.10c)\begin{gather}E_{sr,l}(t) = \tfrac{1}{2} \langle \langle v_l \rangle_x^2 + \langle w_l \rangle_x^2 \rangle _z |_0^{80}, \end{gather}

(3.10d)\begin{gather}E_{wr,l}(t) = \tfrac{1}{2} \langle (v_l - \langle v_l \rangle_x)^2 + (w_l - \langle w_l \rangle_x)^2 \rangle _{x,z} |_0^{80}, \end{gather}

and small scales

(3.11a)\begin{gather} E_{ss,s}(t)= \tfrac{1}{2} \langle \langle u_s \rangle_x^2 \rangle _z |_0^{40}, \end{gather}

(3.11b)\begin{gather}E_{ws,s}(t)= \tfrac{1}{2} \langle (u_s - \langle u_s \rangle_x)^2 \rangle _{x,z} |_0^{40}, \end{gather}

(3.11c)\begin{gather}E_{sr,s}(t)= \tfrac{1}{2} \langle \langle v_s \rangle_x^2 + \langle w_s \rangle_x^2 \rangle_z|_0^{40}, \end{gather}

(3.11d)\begin{gather}E_{wr,s}(t)= \tfrac{1}{2} \langle (v_s - \langle v_s \rangle_x)^2 + (w_s - \langle w_s \rangle_x)^2 \rangle _{x,z} |_0^{40}, \end{gather}

where the shorthand $f(t)|_{y_1}^{y_2}$ denotes averaging of quantity $f(t)$ in the wall-normal direction over the interval $[y_1, y_2]$ given in wall units. The choice of $y^+=40$ to separate the boundary of the near-wall region (viscous sublayer and buffer layer) from the logarithmic region is chosen to approximately separate phenomenological differences observed in these regions based on figures 1–3. In practice, small changes in the choice of this separating $y^+$ value was seen to have little impact on the results of the following sections.

The cross-correlation of these energies at both scales has been computed and is shown in figure 4. Indeed, correlation does exist between the different kinetic energies, with the shift of the peak of the correlations revealing the order of the processes: $E_{sr} \rightarrow E_{ss} \rightarrow E_{ws} \rightarrow E_{wr} \rightarrow E_{sr}$, consistent with the results of Hamilton et al. (1995), Hwang & Bengana (2016) and Doohan et al. (2021), indicating that the self-sustaining process at both scales is captured by the ROM. Further, the time scale of these processes both at the large and small scales is also similar to that found by Doohan et al. (2021). The small-scale process occurs on a time scale that is significantly shorter than that of its large-scale counterpart, as seen by the width of the cross-correlation functions.

Figure 4. Temporal cross-correlation between the energy observables associated with the structural components of the self-sustaining process at (a) large and (b) small scales. Markers show the peak of each cross-correlation.

Similar correlations have been produced for the other energy transport processes identified by Doohan et al. (2021), namely for the energy cascade, driving and feeding processes, all showing agreement with those found in the original DNS. The details of these correlations can be found in the supplemental material of this paper. In summary, this system contains an energy cascade which transports energy from the large-scale production to small-scale wall-normal and spanwise dissipations mainly through a particular part of the wall-normal inter-scale turbulent transport

(3.12)\begin{equation} T_{v,\updownarrow c} = \langle v_l(\boldsymbol{u}_l\boldsymbol{\cdot} \boldsymbol{\nabla} v_s) \rangle_{x,z}. \end{equation}

Energy transported from large to small scales also involves the driving mechanism, the activation process of small-scale production via energy cascade. The related energy transport is well characterised by the other part of the wall-normal inter-scale turbulent transport

(3.13)\begin{equation} T_{v,\updownarrow d} ={-}\langle v_s(\boldsymbol{u}_s\boldsymbol{\cdot} \boldsymbol{\nabla} v_l) \rangle_{x,z}. \end{equation}

It results in subsequent activation of the small-scale production via the Orr mechanism (Orr 1907), promoting a range of other processes such as an increase in small-scale dissipation in response to the breakdown of small-scale wavy streaks as well as a feeding process which transports energy from the small to large scales. Finally, the feeding process is also seen to transport energy from small-scale production mainly to large-scale dissipation terms through the streamwise and spanwise inter-scale turbulent transport, and it involves the activation of small-scale wavy streaks and small-scale wavy rolls (through the spanwise pressure strain).

4. Equilibria and periodic orbits

Having validated the ROM by ensuring that it captures the time-averaged energy budget both at large and small scales, as well as capturing the known processes associated with the temporal energy dynamics of the flow, we now search for invariant solutions using the model.

Equilibria and periodic orbits were searched using a trust-region algorithm, implemented in the fsolve function in MATLAB. Equilibria are steady states of (2.6), such that

(4.1)\begin{equation} \frac{1}{Re} \sum_{j=1}^N L_{ij} a_j + \sum_{j=1}^N \tilde{L}_{ij} a_j + \sum_{j=1}^N \sum_{k=1}^N Q_{ijk} a_j a_k = 0, \end{equation}

whereas periodic orbits are solutions of (2.6) satisfying $a_j(t+T) = a_j(t)$, $j=1, \ldots, N$. Equilibria may be sought by looking directly for the zeros of (4.1). Periodic orbits are searched as zeros of the normalised residual

(4.2)\begin{equation} r(\boldsymbol{a}) = \frac{\| \boldsymbol{a}(t+T) - \boldsymbol{a}(t) \|}{\| \boldsymbol{a}(t)\| }, \end{equation}

where the $\boldsymbol {a}$ vector collects the temporal coefficients $a_j(t)$ of the ROM.

In order to avoid travelling waves and relative periodic orbits, and thus simplify the search for invariant solutions, the reflection symmetry $[u,v,w](x,y,z) = [-u,-v,-w](-x,-y,-z)$ was imposed in the ROM; this is accomplished by truncating the ROM to the modes that satisfy the said symmetry, resulting in 300 degrees of freedom. This lower dimension simplifies the search for invariant solutions. When looking for equilibria, the analytical Jacobian of (4.1) was supplied to the Newton solver, which greatly accelerates the search by reducing the number of function calls. Periodic orbits require time integrations to obtain the Jacobian, and using the direct linearisation of the ROM did not lead to clear benefits in our tests; we have thus resorted to a finite-difference approximation of the Jacobian of the normalised error of (4.2). The tolerance of the solver was specified such that (4.1) is correct to within $10^{-8}$ for equilibria, and that the residual is lower than $10^{-5}$ for periodic orbits – the final residuals of many periodic orbits were of $O(10^{-6})-O(10^{-7})$. Lastly, the convergence of time period, which was updated in the Newton iterations, was not explicitly monitored, although an approach for this was recently proposed (Page et al. 2022).

Using this method, starting from initial guesses taken from a recurrence analysis of chaotic trajectories, 96 equilibria, named EQ1 to EQ96, and 43 periodic orbits, named PO1 to PO43, have been computed. The period of the periodic orbits found varies in the range $T^+ \in [26.53, 141.57]$. In the following sections, the inner units are formulated with respect to the friction velocity of the ROM for each invariant solution and not the friction velocity of each invariant solution unless specified otherwise.

Equilibria and periodic orbits computed in the similar two-scale system of Doohan et al. (2022) through DNS were not seen to capture the full multi-scale behaviour of the system. They typically captured the dynamics of only a single integral length scale. The solutions that appeared to capture any of the multi-scale dynamics were from the upper branch, i.e. the solution family characterised by a large wall shear stress (i.e. high $Re_\tau$). Even still, these upper branch solutions failed to capture the full multi-scale dynamics of the long-time chaotic state.

To compare the solutions of Doohan et al. (2022) with those found in this study using the ROM, we have constructed phase portraits of similar observables, as well as the mean and r.m.s. statistics of the invariant solutions, which are presented in the following sections, first for the equilibria, followed by the periodic orbits. This is followed by a discussion of periodic orbits identified to contain essential dynamics relating to key physical processes in the flow.

4.1. Equilibria

The mean and r.m.s. velocity profiles of the computed equilibrium solutions as functions of the wall-normal direction are shown in figure 5, compared with a long-time trajectory from the ROM shown in black and the DNS in blue. Four of these solutions have been highlighted in figure 5 which closely capture the wall shear stress of the long-time trajectory. Phase portraits showing similar energy observables to those of Doohan et al. (2022) are shown in figure 6, also highlighting the same equilibria as in the velocity profiles. Here, we note that many of the computed equilibria have non-negligible amounts of energy and production at both large and small scales (figure 6a–d) and, in this sense, they may be viewed to be multi-scale equilibrium solutions.

Figure 5. Mean and r.m.s. velocity profiles for the equilibrium solutions (grey) compared with the long-time trajectory generated by the ROM (black) and DNS (blue): (a) $U^+$; (b) $u_{rms}^+$; (c) $v_{rms}^+$; (d) $w_{rms}^+$. Highlighted equilibria are EQ8 (orange); EQ28 (cyan); EQ48 (purple); EQ52 (green).

Figure 6. Phase portraits of the large- and small-scale energy observables of the (a–d) self-sustaining process (equations (3.10), (3.11)) and (e–i) the two scale energy budget (equations (3.4), (3.5), (3.7)) for the equilibria (blue), EQ8 (orange), EQ28 (cyan), EQ48 (purple), EQ52 (green), compared with linearly spaced contours of the respective bivariate probability density function (dark–light grey) generated from a long-time trajectory of the ROM.

A large number of the equilibria are found to capture the mean flow reasonably well, especially in the near-wall region (figure 5a). However, many did not accurately capture the velocity fluctuations of the flow and are not seen to be highly representative of the chaotic state with respect to observables other than the mean flow. The equilibria found are typically seen to overestimate the near-wall streamwise velocity fluctuations. This is clearly observed in the maxima of the streamwise r.m.s. profiles (figure 5b), which appear to be related to the large amount of energy associated with straight streaks in figure 6(a). Equilibria overestimating the amount of energy related to straight streaks compared with the long-time trajectory were also observed by Doohan et al. (2022) with the relative magnitude of the overestimation being similar, especially for the small-scale straight streaks.

Wall-normal fluctuations are seen to be of generally lower magnitude than that of the long-time trajectory, although some of the highlighted equilibria with high friction Reynolds numbers are seen to capture the wall-normal r.m.s. velocity profiles quite well (EQ8, EQ28 and EQ48; figure 5c). With an accurate reproduction of the mean wall shear stress, some of the computed equilibria capture the energy production of the large and small scales reasonably well (figure 6e). The spanwise r.m.s. velocities of the computed equilibria are also typically lower than that of the chaotic trajectory (figure 5d). These observations are consistent with some of the phase portraits, in which straight rolls, wavy streaks and rolls are not well captured by the equilibria (figure 6b–d).

Finally, the phase portraits associated with energy cascade and nonlinear scale interactions show that most of the equilibria computed do not properly capture the dynamics of the chaotic trajectory. The streamwise dissipation of many equilibria does not appear near the chaotic trajectory (figure 6f). Furthermore, the phase portraits of the wall-normal and spanwise dissipation (figure 6g), inter-scale energy transport (figure 6h) and pressure strain (figure 6i) all show that most of the computed equilibria have much smaller amplitudes than those of the instantaneous chaotic state. This is consistent with the finding of the DNS study by Doohan et al. (2021). However, it should be pointed out that most of their equilibria feature the self-sustaining process either at large or small scale, whereas many equilibria in this study have non-negligible amounts of energy and production at both of the scales (figure 6a,b,e) – note that the equilibria of Doohan et al. (2021) were obtained by continuing the existing solutions obtained in a smaller computational domain (equivalent to lower Reynolds number), whereas those in our study are obtained by directly looking for the zeros of (4.1). This observation suggests that the relatively low amplitude of the observables associated with energy cascade and nonlinear scale interactions is possibly more associated with the lack of unsteadiness in the equilibria, as we shall indeed see in § 4.2.

4.2. Periodic orbits

Similar data as shown for the equilibria are now shown for the periodic orbits computed using the ROM, with mean and r.m.s. velocity profiles shown in figure 7, and phase portraits of large- and small-scale energy observables in figure 9. Highlighted periodic orbits in figures 7 and 9 are chosen as they show reasonable agreement with the mean and r.m.s. statistics, and their dynamics in the state space will be further explored with respect to figure 9.

Figure 7. Mean and r.m.s. velocity profiles for the periodic orbits (grey) compared with the long-time trajectory generated by the ROM (black) and DNS (blue): (a) $U^+$; (b) $u_{rms}^+$; (c) $v_{rms}^+$; (d) $w_{rms}^+$. Highlighted periodic orbits: PO17 (cyan); PO20 (purple); PO26 (green); PO42 (orange).

In general, the periodic orbits characterise the chaotic state more accurately than the equilibria, showing closer agreement to the long-time trajectory in many cases. The mean profiles displayed by the periodic orbits appear to be similar to those shown by the equilibria, with the wall shear stress being captured well (figure 7a). Large maxima seen in the near-wall region of the streamwise r.m.s. velocity profiles in the equilibria are less apparent in the periodic orbits with a multitude of periodic orbits capturing this profile well (compare figure 7b with 5b). Similar to the equilibria, the wall-normal fluctuations displayed by the periodic orbits appear to be mostly smaller than that of the long-time trajectory or DNS, although wall-normal fluctuations in the near-wall region are seen to be a little closer to the long-time trajectory than those of the equilibria (compare figure 7c with 5c). The spanwise r.m.s. profiles also generally exhibit values less than that of the long-time trajectory, although some periodic orbits replicate the near-wall peaks better than those of the equilibria: PO17 and PO26, in particular (compare figure 7c with 5c).

Some recent studies have suggested that taking simple averages over computed invariant solutions (Chandler & Kerswell 2013; Page et al. 2022) is often found to be as good as other more sophisticated approaches (cyclic expansion; Cvitanovic et al. 2005). In figure 8, turbulence statistics obtained from a simple average of equilibria and periodic orbits are further compared with those of the long-time trajectory of DNS and ROM. While there is very little difference in the mean profiles from the simple average of equilibria and periodic orbits (figure 8a), the velocity fluctuations from periodic orbits show better agreements with those from the long-time trajectory of DNS and ROM (figure 8b–d). However, the velocity fluctuations from periodic orbits still show considerable differences from those from the long-time trajectory of ROM, indicating that there are a number of orbits which may not be near the turbulent solution trajectory in the state space.

Figure 8. Mean and r.m.s. velocity profiles from simple average over the equilibria (red) and the periodic orbits (green) compared with those of the long-time trajectory generated by the ROM (black) and DNS (blue): (a) $U^+$; (b) $u_{rms}^+$; (c) $v_{rms}^+$; (d) $w_{rms}^+$.

The periodic orbits generally appear to capture the cross-streamwise velocity fluctuations better than the equilibria, as seen in figure 7. The observables associated with energy transport between different velocity components would therefore be better captured by them, and pressure strain has been understood to be mainly involved in this process (e.g. Tennekes & Lumley 1967; Cho et al. 2018). Given the form of pressure strain shown in (A4) and (A11), the state-space projections of the periodic orbits with the observables involving $\partial u/\partial x$ are expected to be improved compared with that of the equilibria. Indeed, an examination of the phase portraits in figure 9 reveals that such observables are more convincingly captured in the periodic orbits at both large and small scales than in the equilibria: for example, wavy streaks (figure 9c), cross-streamwise dissipation (figure 9g) and pressure strain (figure 9i), in particular.

Figure 9. Phase portraits of the large- and small-scale energy observables of the (a–d) self-sustaining process (equations (3.10), (3.11)) and (e–i) the two scale energy budget (equations (3.4), (3.5), (3.7)) for the periodic orbits (blue), PO17 (cyan), PO20 (purple), PO26 (green), PO42 (orange), compared with linearly spaced contours of the respective bivariate probability density function (dark–light grey) generated from a long-time trajectory of the ROM.

Overall, the periodic orbits describe the chaotic trajectory both statistically and dynamically better than the equilibria. The periodic orbits are seen to be more representative with respect to the energy of straight streaks in the flow, with less high-energy small-scale outliers, although these high-energy outliers continue to be seen at the large scales (figure 9a). There is a modest improvement in the description of rolls by the periodic orbits (figure 9b,d). Production and dissipation appear to be more consistently captured by the periodic orbits than the equilibria (figure 9e–g). The inter-scale energy transport terms from the periodic obits are also able to explore more of the state space than those from the equilibria.

Despite the improved state-space dynamics offered by the periodic orbits compared with that by the equilibria, the periodic orbits computed in this study do not sufficiently cover the state space visited by the chaotic trajectory, as is evident from figure 9. Imposition of the reflection symmetry might be a potential reason, but the phase portrait and turbulence statistics from the chaotic trajectory subject to the symmetry did not yield any significant difference to those without the symmetry (see Appendix B). It is, however, important to mention that this does not necessarily imply a ‘theoretical’ failure of the attempt to construct the chaotic dynamics in terms of periodic orbits. As mentioned earlier, the periods of the computed orbits are less than $T^+ \leq 141.57$. However, the time scale of the self-sustaining process at large scale is expected to be $T^+ \sim 400\unicode{x2013}500$ (Hwang & Bengana 2016; Doohan et al. 2021), considerably longer than the periods of the computed orbits in this study. Therefore, it is presumable that the lack of state-space coverage by the computed orbits originates from their relatively short periods, and the computation of the periodic orbits with longer periods is expected to improve the state-space description by the periodic orbits. Nevertheless, the computation of orbits with longer periods has been found to be challenging even with the present ROM. In particular, we have observed the vanishing recurrence of the chaotic trajectory on increasing Reynolds numbers, causing trouble in identifying suitable initial conditions for the Newton algorithm. This also suggests that a numerical approach with a fundamentally different concept would probably be required to overcome this challenge. The Newton method employed for the computation of periodic orbits in this study is expected to have limited applicability, when the leading Lyapunov exponent of the given system is too large – the leading Lyapunov exponent of the present ROM computed using the algorithm in Parker & Chua (1989) was $\lambda _1=0.107$ ($\lambda _1^+=0.58$), comparable to the leading eigenvalues of the periodic orbits in table 1. This is also a common feature of some of the most popular algorithms currently available based on the same type of Newton iteration (e.g. Viswanath 2007; Willis et al. 2016). There have been a few recent works addressing this challenge, and the proposed algorithms do not appear to depend on the magnitude of the leading Lyapunov exponent (e.g. Azimi et al. 2022; Page et al. 2022; Parker & Schneider 2022). Employing such a technique may be a way to overcome this challenge in the future.

Table 1. Friction Reynolds number $Re_\tau$ (with respect to the friction velocity of each orbit), period $T^+$ (with respect to the friction velocity of the ROM), maximum streamwise r.m.s. velocity $\max \{u_{rms}^+\}$, the dimension of the unstable manifold $\dim (\mathcal {W}^u)$, the largest Floquet multiplier $\mu _1$ and corresponding Lyapunov exponent $\lambda _1$ of selected periodic orbits, compared with the ROM whenever possible.

4.3. Multi-scale dynamics in periodic orbits

Thus far, we have seen that a few periodic orbits computed appear to provide some multi-scale dynamics representative of the chaotic state, consistent with the results of figure 9, which show contributions in energy terms from both the large and small scales. Now, we will consider some periodic orbits, highlighted with distinct colours in figures 7 and 9 (PO17, PO20, PO26, PO42, in particular), to explore the exact dynamics that each orbit offers. We shall see that these orbits describe: (i) the self-sustaining process at small and large scales (Hwang & Bengana 2016); (ii) the driving mechanism of small-scale production via energy (Doohan et al. 2021); (iii) energy cascade from large to small scales. An overview of these chosen orbits is tabulated in table 1. In particular, their stability properties show a large number of unstable modes with a high leading Lyapunov exponent, with values much larger than those reported at low Reynolds numbers (Viswanath 2007). We also note the large dimension of the unstable manifold of these orbits, comparable to or larger than that of the multi-scale equilibrium solution found in Doohan et al. (2022) (i.e. $L_{L2b}$). Invariant solutions with a high-dimensional unstable manifold have often been observed to more closely resemble the turbulent state (see also Doohan et al. 2019).

Figure 10 shows the large- and small-scale velocity fields at an instant in time in PO42. It is clear that a large pair of streaks that extends from the wall into the centre of the domain dominates the large-scale field, with twice as many small-scale streaks dominating the near wall of the small-scale velocity field, in agreement with the phenomenology of the large- and small-scale self-sustaining processes. This structure in the large- and small-scale velocity fields, common to a number of the computed periodic orbits, combines to exhibit a rich variety of dynamics.

Figure 10. (a) Instantaneous large-scale streamwise velocity isosurfaces $u_l^+= \pm 1.5$ (red/blue) and (b) instantaneous small-scale streamwise velocity isosurfaces $u_s^+= \pm 0.75$ (red/blue) shown at the same instant in time in PO42.

A number of the spatio-temporal processes statistically identified in the long-time chaotic trajectory have also been found in the periodic orbits. Evidence of the small-scale self-sustaining process has been seen in a large number of the periodic orbits, with an example (PO20) shown in figure 11. The time trace of the observables of the small-scale self-sustaining process shows a consistent interchange of dominance between its structural components: $E_{sr,s} \rightarrow E_{ss,s} \rightarrow E_{ws,s} \rightarrow E_{wr,s} \rightarrow E_{sr,s}$. The flow field reconstructions of figure 11 show the evolution of two pairs of streaks (red/blue $u_s^+$ contours) and vortices (black $\omega _{x,s}$ contours), showing straight streaks undergoing an instability, resulting in a wavy streak. This wavy streak is accompanied by a growth of intensity of the neighbouring streamwise vortices which slowly begin to grow between approximately $t^+=60$ and $t^+=10$ in the absence of streaks. During this time, straight streaks are seen to regenerate and the cycle continues.

Figure 11. (a) The evolution of the small-scale self-sustaining process observables normalised by their maximum values over a single period of PO20. (b–d) Flow field visualisations at selected times show small-scale streamwise velocity isosurfaces $u_s^+=\pm 1.5$ (red/blue), and small-scale streamwise vorticity isosurfaces $\omega _{x,s}=\pm 0.55$ (black).

Due to the increased difficulty and expense of converging longer periodic orbits, evidence of the large-scale self-sustaining process has only been found in a single longer orbit (PO42), shown in figure 12. As in the case of the small-scale self-sustaining process, we see a consistent interchange of dominance between the structural components of the large-scale self-sustaining process: $E_{sr,l} \rightarrow E_{ss,l} \rightarrow E_{ws,l} \rightarrow E_{wr,l} \rightarrow E_{sr,l}$. A single pair of large-scale streaks (red/blue $u_l^+$ contours) are seen to dominate the flow field, starting as straight streaks, slowly undergoing an instability to form a wavy streak, resulting in a breakdown in the streak along its length. From this point, a large intensity of large-scale streamwise vorticity is formed which sits alongside the remnants of the broken streak. This vortex is then seen to stretch in the streamwise direction (Schoppa & Hussain 2002) accompanying the regeneration of the original straight streak occupying the full streamwise length of the domain.

Figure 12. (a) The evolution of the large-scale self-sustaining process observables normalised by their maximum values over a single period of PO42. (b–d) Flow field visualisations at selected times show large-scale streamwise velocity isosurfaces $u_l^+=\pm 1.5$ (red/blue), and large-scale streamwise vorticity isosurfaces $\omega _{x,l}=\pm 1$ (black).

Other processes, such as the energy cascade and the driving mechanism, identified statistically in long-time chaotic trajectories of the ROM (§ 3.3) as well as in the DNS of Doohan et al. (2021), have also been identified in a number of the periodic orbits. An example of the driving mechanism is observed in PO26, shown in figure 13. The first snapshot ($t^+=10$) shows a point of inactivity in the cycle featuring small amounts of small-scale production accompanying some high $u_s^+$ (streaky) structures (figure 13b). By the second snapshot at $t^+=40$, a large peak in the driving wall-normal inter-scale turbulent transport (see (3.13)) from large to small scale (i.e. positive $T_{v,\updownarrow d}^+$) is observed both in the centre of the channel and around the streaks, transporting energy from large to small scales. This is accompanied by a rise in turbulent production (figure 13c), due to an injection of energy into the wall-normal small-scale velocity $v_s^+$ (see also figure 13a). The result of this is an amplification of the small-scale streaky motions observed at $t^+=60$ (figure 13d).

Figure 13. (a) The evolution of the small-scale streak energy $E_{ss,s}^+$, driving component of wall-normal inter-scale turbulent transport $T_{v,\updownarrow d}^+$ and small-scale production $P_{us}^+$ normalised by their maximum values over a single period of PO26, showing the driving mechanism. (b–d) The flow field visualisations at selected times show small-scale streamwise velocity isosurfaces $u_s^+=\pm 2.25$ (red/blue), driving wall-normal inter-scale turbulent transport isosurfaces $T_{v,\updownarrow d}^+=0.11$ (black) and small-scale production isosurfaces $P_{us}^+=0.67$ (grey).

An example of the energy cascade is observed in PO17, shown in figure 14. In the statistical analysis of the long-time chaotic trajectory, energy transport primarily through a particular part of the wall-normal inter-scale turbulent transport associated with the cascade $T_{v,\updownarrow c}$ (see (3.12)) was seen to transport energy from the large-scale production to the wall-normal small-scale dissipation near $y \pm 0.5$. Thus, figure 14 shows snapshots of the flow in the central region of the channel. Large-scale production isosurfaces reveal long streak-like structures in the centre of the channel at the peak of the large-scale production at around $t^+=12$ (figure 14b). Between $t^+=12$ and $t^+=30$ these structures shorten slightly in length accompanied by a growth in regions of high wall-normal inter-scale turbulent transport (figure 14c), which are seen to move down and intersect with the small-scale dissipation contours at $t^+=35$, continuing until roughly $t^+=50$ (figure 14d,e). At this point, an increase in the wall-normal small-scale dissipation is observed, with a drift of the small-scale dissipation isosurfaces towards the central region of the channel. It is also interesting to note that $T_{v,\updownarrow c}^+<0$ during this time interval, indicating the possible existence of an inverse energy cascade, although the time-averaged value of $T_{v,\updownarrow c}^+$ is positive. From $t^+=50$ to $t^+=60$ (figure 14f), the large-scale production increases again, correlated to the intensification of large-scale straight streaks, completing the cycle.

Figure 14. (a) The evolution of the large-scale production $P_{ul}^+$, cascade component of the wall-normal inter-scale turbulent transport $T_{v,\updownarrow c}^+$ and wall-normal small-scale dissipation $\varepsilon _{vs}^+$ normalised by their maximum values over a single period of PO17, showing the energy cascade. (b–d) The flow field visualisations at selected times show large-scale production isosurfaces $P_{ul}^+=0.22$ (grey), wall-normal inter-scale turbulent transport isosurfaces $T_{v,\updownarrow c}^+=0.22$ (black) and wall-normal and spanwise small-scale dissipation $\varepsilon _{vs}^++\varepsilon _{ws}^+=-0.3$ (cyan).

5. Conclusions

In this paper, we have validated the model of Cavalieri & Nogueira (2022) for plane Couette flow retaining only two integral length scales through a two-part process. The first involved the examination of the time-averaged two-scale energy budget based on a split by spanwise wavenumber, including accessing the exchange of energy within and between scales. Overall, a qualitative agreement was seen between the ROM and the DNS results of Doohan et al. (2021) in nearly every respect. The second part of the model validation process involved examining these processes further by studying the temporal cross-correlations of various energy measures integrated over the spatial domain. Identified in the model was both a large- and small-scale self-sustaining process, an energy cascade from large to small scales, the driving of small-scale production via energy transport from the large scales in the near-wall region and a feeding process transporting energy from the small to large scales. The underlying mechanisms of these processes, at least in relation to the various energy budget terms and the structural components of the self-sustaining process at each scale, have been analysed and compared with the results of Doohan et al. (2021), with good qualitative agreement.

A total of 96 equilibria and 43 periodic orbits have subsequently been computed for the ROM, the statistics and phase portraits of which have been directly compared with those of a long-time chaotic trajectory generated using the ROM. In general, the invariant solutions of the ROM reproduce its multi-scale dynamics more convincingly than those found by Doohan et al. (2022) for the full Navier–Stokes equations. Equilibrium solutions were largely seen to capture straight streaks and rolls, resulting in large near-wall peaks in the streamwise r.m.s. velocity profiles, with high $Re_\tau$ equilibria seen to exhibit energy production values very similar to that of the mean long-time trajectory. In general, equilibria failed to accurately capture wavy streaks and rolls, as well as the wall-normal and spanwise energy dissipation. The periodic orbits were seen to characterise the chaotic state more accurately than the equilibria, showing closer agreement to the long-time trajectory in nearly all cases, exhibiting a rich range of multi-scale dynamics and spatio-temporal evolution.

The periodic orbits in particular have been seen to contain useful dynamical information about the known processes in the flow, which have been more easily identified than by using a statistical approach applied to the long-time chaotic trajectory. In particular, various periodic orbits have contained valuable information, including the large- and small-scale self-sustaining process, the driving mechanism which transports energy from the large scales to the small-scale production and the energy cascade, allowing for a more clear physical interpretation of these processes. It might be seen as evident that periodic orbits would be more useful and relevant than equilibrium solutions for the statistical and dynamical description of turbulence. However, it is probably worth pointing out that, at lower Reynolds numbers, where the flow is characterised by a single integral length scale (i.e. the minimal flow unit), many equilibrium solutions in plane Couette flow have also previously been found to be organised around the turbulent trajectory in the state space (e.g. Gibson et al. 2008; Doohan et al. 2019). These equilibrium solutions often contain the most fundamental physical ingredients for the underlying physical mechanism of sustaining turbulence (i.e. the self-sustaining process). Similarly, in Rayleigh–Bénard convection, the recent study by Motoki et al. (2021) showed that there exists a multi-scale equilibrium solution which resembles turbulence. Even from a dynamical systems viewpoint, the equilibrium solutions can still act as a signpost for the chaotic solution trajectory through the formation of homoclinic/heteroclinic orbits. From this perspective, it is interesting to observe that a large number of the equilibrium solutions in a wall-bounded shear flow lose their statistical and dynamical importance on increasing the Reynolds number. It remains to be seen if this is a general feature in other flows.

As a final remark, the invariant solutions computed using the ROM are not invariant solutions of the exact system, and this is an important limitation of this study. However, the computation of invariant solutions to the full Navier–Stokes equations is prohibitively expensive and very challenging for high Reynolds numbers. The present study has demonstrated that the significantly reduced cost of generating invariant solutions for a carefully crafted ROM would provide a gateway to applying dynamical systems approaches to flows at high Reynolds numbers and exploring the increasingly complex multi-scale turbulent dynamics identified in the full-scale dynamics – indeed, the dramatically reduced computational cost has enabled us to carry out an extensive search of unstable periodic orbits buried in the turbulent state, a task that is likely impossible to tackle using DNS at the current computational cost for the given Reynolds number. In this respect, the key to the application of this approach is in the construction of a successful ROM, which is able to capture the turbulent dynamics of interest. Thankfully, there are an increasing number of studies towards the reliable construction of ROMs for turbulent flows (e.g. see the review by Rowley & Dawson 2017). For the purpose of extracting the key turbulent dynamics through invariant solutions, the ROM must be in the form of a ‘nonlinear’ dynamical system, and it is highly desirable to preserve the mathematical structure of the Navier–Stokes equations so that the dynamics described by invariant solutions of the ROM can be analysed in comparison with that of the real flow. By doing so, the invariant solutions would be able to describe the exact ‘reduced-order’ dynamics of the given turbulent flow and the underlying complex physical processes can be studied by formulating a suitable energy budget analysis, a key tool to explore energy cascade and scale interactions in turbulence. From this perspective, it is also important to carefully control the construction of a ROM by recognising its exact capability in the description of the turbulent dynamics. Lastly, it is worth noting that the invariant solutions from a ROM could be continued with respect to the neglected modes of the exact system, provided that the ROM has a sufficient number of degrees of freedom. This is a direction we currently pursue towards a cost-effective method for the computation of physically meaningful invariant solutions to the Navier–Stokes equations.