Improved assessment of the statistical stability of turbulent flows using extended Orr–Sommerfeld stability analysis

Abstract

The concept of statistical stability is central to Malkus's 1956 attempt to predict the mean profile in shear flow turbulence. Here we discuss how his original attempt to assess this – an Orr–Sommerfeld (OS) analysis on the mean profile – can be improved by considering a cumulant expansion of the Navier–Stokes equations. Focusing on the simplest non-trivial closure (commonly referred to as CE2) that corresponds to the quasilinearized Navier–Stokes equations, we develop an extended OS analysis that also incorporates information about the fluctuation field. A more practical version of this – minimally extended OS analysis – is identified and tested on a number of statistically steady and, therefore, statistically stable turbulent channel flows. Beyond the concept of statistical stability, this extended stability analysis should also improve the popular approach of mean flow linear analysis in time-dependent shear flows by including more information about the underlying flow in its predictions as well as for other flows with additional physics such as convection.

1. Introduction

This paper revisits the seminal work of Malkus (1956) that attempted to build a theory of shear turbulence. This theory was based upon maximising the momentum transport (or, equivalently, dissipation rate) achieved by the flow amongst all those with a marginally stable mean profile. Malkus clearly had a statistical form of marginal stability in mind but, to make progress, had to resort to specifying marginality with respect to the then 50-year-old Orr–Sommerfeld (OS) equation (Orr 1907; Sommerfeld 1908). So posed, his marginal stability idea was quickly repudiated (Reynolds & Tiederman 1967 and, more recently, Iyer et al. 2019) although further studies showed it could be made to work if an anisotropic eddy viscosity model was used (Reynolds & Hussain 1972; Malkus 1979; Sen & Veeravalli 2000; Sen et al. 2007). The concept of statistical stability was, however, central to Malkus's thinking (remaining so throughout his career, e.g. Malkus 1996 and 2003) and is clearly different from stability as viewed within the context of the governing Navier–Stokes equations (epitomised by the celebrated OS equation). For example, it is fairly uncontentious to assert that the turbulent attractor in, say, pressure-driven channel flow at a high enough Reynolds number and large enough domain has stationary statistics (defined by averaging over one or more homogeneous directions or in an ensemble sense) and so within the partial differential equations (PDEs) that govern how these statistics evolve, the realised turbulence is a stable fixed point – i.e. turbulence is statistically stable to infinitesimal perturbations of the statistics. This is in contrast to the time-dependent turbulent attractor as viewed in the Navier–Stokes equations where adding a small disturbance may well show that disturbance grows and never decays yet the original statistics still recover (an example is shown herein). The difference, of course, is that a flow disturbance can have a component along the turbulent trajectory and, hence, acts as a time shift: this part of the disturbance never decays to zero but does not affect the statistics (see, e.g. Nikitin 2018). For his theory, Malkus wanted a statistical stability criterion based only on the lowest-order statistic – the mean flow. Ideally more statistical information needs to be incorporated and it is our objective here to attempt this. At the very least, doing so should improve the now-standard approach of carrying out linear stability analysis of the mean profile of time-dependent flows in an attempt to understand observed coherent structures (see, e.g. Crighton & Gaster 1976; Gaster, Kit & Wygnanski 1985; Roshko 1993; Barkley 2006; Lesshafft et al. 2006; Sipp & Lebedev 2007; Akervik et al. 2008; Sipp et al. 2010; Mantic-Lugo, Arratia & Gallaire 2014; Beneddine et al. 2016; Lefauve et al. 2018).

The motivation for this work comes from two different directions. The first is the general approach of applying linear analysis around the mean profile of a time varying, possibly turbulent flow to deduce information about the likely fluctuations seen. Initially, this took the form of linear stability analysis stimulated by Malkus's work that tends to work well in free-shear turbulent flows like jets where inviscid (inflectional) instabilities dominate (Crighton & Gaster 1976; Gaster et al. 1985; Roshko 1993), but less well in wall-bounded situations such as channel flow where viscosity can be important (although there have been successes, e.g. Barkley 2006; Sipp & Lebedev 2007; Beneddine et al. 2016). Driven by the fact that shear flow mean profiles tend to be linearly stable, this approach subsequently diversified into non-modal analysis (Butler & Farrell 1993; del Alamo & Jimenez 2006; Cossu, Pujals & Depardon 2009) and input–output or resolvent analyses (Chernyshenko & Baig 2005; Jovanovic & Bamieh 2005; Hwang & Cossu 2010a,b; McKeon & Sharma 2010; Moarref & Jovanovic 2012; Blesbois et al. 2013; Sharma & McKeon 2013; Schmidt et al. 2018). The resolvent approach has been particularly illuminating in showing exactly how far linear analysis can go given the mean flow profile (McKeon 2017; Jovanovic 2021). Predictions can be made of the dominant fluctuation response at a given frequency and wavevector that typically resonates with observations at least up to amplitude and phase. In some sense this solves ‘half’ the problem of turbulence – given a mean flow, linear analysis around this can extract the dominant fluctuation structures – and refocuses attention on the ‘other’ half – predicting the mean profile. Identifying the required amplitudes and phases of all the fluctuation fields to support the observed mean profile, however, may approach the difficulty of ‘just’ solving the Navier–Stokes equations. One alternative is to appeal to some simpler ‘organising’ physical principle as epitomised in Malkus (1956) and then the idea of statistical stability seems a key concept.

The second motivation is the resurgence of interest in dealing with statistical quantities directly through a cumulant expansion of time-varying flows (Hopf 1952; Orszag 1977; Frisch 1995). Generating the evolution equations for cumulants immediately highlights the closure problem of turbulence in which the time derivative of a $n$th-order cumulant requires knowledge of the $(n+1)$th cumulant and so forth. Present-day computational power has, however, renewed interest in the pursuit of simple cumulant-discard closures to see how well they do in modelling flows in a variety of contexts, for example, in atmospheric dynamics (Farrell & Ioannou 2003; Marston, Conover & Schneider 2008; Srinivasan & Young 2012; Tobias & Marston 2013; Parker & Krommes 2014), astrophysics (Tobias, Dagon & Marston 2011), plasmas (Farrell & Ioannou 2009; Parker & Krommes 2013) and wall-bounded shear flows (Farrell & Ioannou 2012; Constantinou et al. 2014; Farrell et al. 2016). The most popular closure ignores third- and higher-order cumulants – commonly referred to as CE2 – and has the nice property of being exact for a quasilinear (or historically a ‘mean field theory’) version of the Navier–Stokes equations (e.g. see the recent review by Marston & Tobias 2022). Significantly for our purposes here, CE2 and higher closures (CE$n$ where $n \geq 3$ is the highest-order cumulant retained) present the most natural framework in which to extend Malkus's idea of statistical stability. In what follows, the focus will be on using CE2 and its relationship to the quasilinear (QL) version of the Navier–Stokes equations to develop an improved version of the usual OS analysis of the mean profile; see figure 1. The overall objective is to develop a way to judge whether a state is statistically stable or not using a subset of its statistics.

Figure 1. Comparison between standard and statistical considerations of the linear stability of turbulent states. On the left are physical space equations (methods are ordered top to bottom by decreasing nonlinearity): full Navier–Stokes (NS) equations, quasilinear (QL) approximation, OS equation around a turbulent mean velocity profile, extended and minimal extended Orr–Sommerfeld equations (EOS and mEOS). On the right are statistical space equations: CE2 equations that include up to second-order statistics. Black arrows indicate the standard path of turbulent flow stability analysis leading to OS; blue arrows indicate the statistical approach to turbulent flow stability analysis leading to EOS and mEOS, emphasizing how the steady statistical state $U,\boldsymbol {C^{mn}}$ can be used to obtain the steady state counterpart $U,\boldsymbol {\tilde {u}_0^{mn}}$ in the physical space.

It is worth emphasizing that the approach taken here is perfectly general and not confined to shear flows or even fluid mechanics. The key idea given a physical system described by a PDE is to generate evolution equations – ‘statistical’ equations – for the first few statistics of the solution. If the solution of the original PDE is spatiotemporally complicated but has steady statistics, the premise pursued here is that a better way to assess the stability of this solution is to look at the linear stability within the statistical equations rather than the original PDEs. The former is an approximation given that a finite number of the system's statistics are considered, but yields a spectral problem built upon a steady state in statistical space. The latter has to rely on a costly ensemble of simulations where a distribution of small perturbations are each added to the solution and their evolution monitored (Iyer et al. 2019; Kashyap, Duguet & Dauchot 2022). A good complementary example of where this approach would be useful is convection (Malkus 1954; Wen et al. 2022).

The plan of this paper is to illustrate the analysis within the context of channel flow described in § 2. The Reynolds number is assumed high enough for the computational box used that the flow can be taken as approximately statistically stationary. An averaging procedure will be assumed below such that the mean flow can be assumed to only depend on the cross-stream variable $y$, i.e. $\boldsymbol {U}=U(y)\boldsymbol {\hat {x}}$. This could be ensemble averaging or averaging over the streamwise ($x$) and spanwise ($z$) directions. Using the latter spatially averaging approach, Malkus pointedly only chose a spanwise average so that his mean profile $\boldsymbol {U}=U(x,y) \boldsymbol {\hat {x}}$ could depend on the cross-stream and streamwise variables. As a result, his OS analysis targeted the stability of a streamwise-independent mean profile $U(y)\boldsymbol {\hat {x}}$ to streamwise-dependent mean flow disturbances $\delta \boldsymbol {U}(y) \exp ({\rm i}k(x-ct))$ rather than fluctuation fields defined as having non-vanishing spanwise dependence. Contrarily, there are growing arguments to only streamwise average to retain the spanwise structure of the mean flow, i.e. $\boldsymbol {U}=U(y,z)\boldsymbol {\hat {x}}$ (e.g. see table 1 of Lozano-Duran et al. (2021) for a sample list of relevant works). Despite this, the focus here is on the simplest mean flow definition for a statistically steady state, $\boldsymbol {U}=U(y)\boldsymbol {\hat {x}}$, given the central role this plays in resolvent flow analysis but, there is no doubt, that extending the mean flow definition is clearly an important direction to extend the approach discussed here. Standard OS analysis is recalled in § 2.1.

Section 2.2 then introduces the cumulant expansion approach (Hopf 1952; Orszag 1977; Frisch 1995) and the hierarchy of evolution equations for these statistical quantities. Solving these equations to (hopefully) reach a steady state is an appealingly direct way to estimate the properties of statistically steady turbulent flows since it avoids having to average across large direct numerical simulation (DNS)-generated datasets. Here, however, the emphasis is on the concept of stability in this statistical framework and its relationship to (dynamic) stability within the Navier–Stokes equations (see § 2.3) not on the accuracy of suitable closures of the cumulant equations in capturing the reference flow state. That is, our strategy is to apply a statistical stability criterion from a statistical closure to a solution of the Navier–Stokes equations to approximate its statistical stability there. Our particular focus will be on the simplest non-trivial cumulant-discard scheme CE2 given the rapidly increasing dimensionality of the approach: CE$n$ works with cumulants up to order $n$ that, before exploiting any symmetries, is typically a $3n$–$2$ rank tensor (three spatial coordinates per field reduced by two averaging directions). Note that CE2 is exact for the quasilinearized Navier–Stokes equations – or QL equations – and translating what a statistically steady state in CE2 means for the QL equations is a crucial step discussed in § 2.4. Appendix A presents an equally important discussion on how the stability predictions within the CE2 and QL systems are related. Quasilinearization (Vedenov, Velikhov & Sagdeev 1961; Herring 1963, 1964) has enjoyed a resurgence of interest recently (see, e.g. Hernandez & Hwang 2020; Skitka, Marston & Fox-Kemper 2020; O'Connor, Lecoanet & Anders 2021; Marston & Tobias 2022) given its accessibility, and together with its ‘generalized’ elaboration (Marston, Chini & Tobias 2016) in which the definition of what constitutes a mean is extended, has the ability to focus on different parts of the nonlinearity in the Navier–Stokes equations (see, e.g. Hernandez, Yang & Hwang 2021).

Section 3.1 then introduces an extended Orr–Sommerfeld analysis – or ‘EOS’ analysis – based on translating the statistical stability problem in CE2 back to the dynamical equations. Intriguingly, this is not the same as just working within the QL approximation (as Appendix A makes clear). Applying EOS analysis carries a substantial overhead so we consider a reduced (practical) version referred to as ‘minimally extended Orr–Sommerfeld analysis’ – or ‘mEOS’ analysis – in § 3.2 that is almost as cheap to apply as OS analysis.

Interestingly, applying the same strategy of mapping a cumulant-based system back to the underlying dynamical equations can only go one level higher in sophistication and requires a jump directly to CE$\infty$. This transforms EOS equations into the familiar linearised Navier–Stokes (LNS) equations albeit based around the steady base state derived by assuming stationary second rank cumulants. This ‘infinitely extended OS (iEOS)’ analysis is described in § 3.2. Section 4 then explores the performance of OS, EOS and mEOS analyses on four different turbulent states realised in two-dimensional (2-D) channel flow that is used as an approximation of a statistically steady flow. The limitations of the analysis are discussed in § 5 followed by a summary and final thoughts in § 6.

2. Formulation: channel flow

For context in this work, we consider a channel flow $\boldsymbol {u}^*(\boldsymbol {x}^*,t^*)$ of a fluid with density $\rho ^*$ and kinematic viscosity $\nu ^*$ between two parallel plates at $y^*=\pm h^*$ across which a time-dependent pressure gradient $9\rho ^* U^{*2}G(t^*)/4h^*\,\boldsymbol {\hat {x}}$ is imposed such that the bulk flow

(2.1)\begin{equation} U^*:= \frac{1}{4h^{*2}L_z} \int^{h^* L_z}_{{-}h^*L_z}\int^{h^*}_{{-}h^*} \,\boldsymbol{u}^*\, {{\rm d} y}^* \,{\rm d}z^* \end{equation}

is kept constant (unstarred/starred quantities are dimensionless/dimensional and periodicity is imposed across the spanwise domain $z^* \in h^*[-L_z,L_z]$). Non- dimensionalizing the Navier–Stokes equations using $h^*$, $3U^*/2$ (so Reynolds numbers based on the bulk speed and the laminar centreline speed $U^{*c}=3U^*/2$ correspond for the laminar parabolic flow) and $\rho ^*$ leads to

(2.2)\begin{equation} \left. \begin{aligned} \boldsymbol{u}_t+ \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} & = G(t)\boldsymbol{\hat{x}}-\boldsymbol{\nabla} p +\frac{1}{Re} \nabla^2 \boldsymbol{u},\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} & = 0, \end{aligned} \right\} \end{equation}

with $({1}/{4L_z})\int ^{L_z}_{-L_z}\int ^1_{-1} u\, \textrm {d}y\,\textrm {d}z=2/3$, where $\boldsymbol {u}=\boldsymbol {u}^*/U^*=u \boldsymbol {\hat {x}}+v \boldsymbol {\hat {y}}+w \boldsymbol {\hat {z}}$, $t=3U^*/(2h^*) t^*$ and $Re:=3U^* h^*/2\nu ^*$. We also impose streamwise periodicity of the flow over $x^* \in h^*[-L_x,L_x]$ so the non-dimensionalised flow domain is $(x,y,z) \in [-L_x,L_x] \times [-1,1] \times [-L_z,L_z]$ with non-slip boundary conditions on the plates at $y=\pm 1$ (fundamental wavenumbers in $x$ and $z$ are labelled $\alpha :={\rm \pi} /L_x$ and $\beta :={\rm \pi} /L_z$, respectively). In terms of an averaging procedure, a number of choices present themselves: ensemble averaging, spatial averaging and time averaging (or even a combination thereof) that all should be equivalent for a statistically stationary system in a large enough domain. However, in what follows we want to treat a numerical experiment in a finite domain with finite $Re$ and then ensemble averaging is the most natural choice as will come clear below. This averaging process is denoted by an overbar, $\overline {(\boldsymbol {\cdot })}$, and then the flow can be decomposed into a mean $U(y,t) \boldsymbol {\hat {x}}:=\overline {\boldsymbol {u}(\boldsymbol {x},t)}$ and fluctuation part $\boldsymbol {\tilde {u}}:=\boldsymbol {u}-\overline {\boldsymbol {u}}=\tilde {u} \boldsymbol {\hat {x}}+\tilde {v} \boldsymbol {\hat {y}} + \tilde {w} \boldsymbol {\hat {z}}$ (due to symmetry, a vanishing mean spanwise component is assumed so $(\overline {\tilde {v} \tilde {w}} )_y=0$).

The Navier–Stokes equations can be similarly decomposed into a mean part,

(2.3)\begin{equation} U_t -\frac{1}{Re}U_{yy} = G - ( \overline{\tilde{u} \tilde{v}} )_y \end{equation}

and a fluctuation part,

(2.4)\begin{equation} \boldsymbol{\tilde{u}}_t = \frac{1}{Re} \nabla^2 \boldsymbol{\tilde{u}} -\boldsymbol{\nabla} \tilde{p} -{U} \boldsymbol{\tilde{u}}_x - \tilde{v} U_y \boldsymbol{\hat{x}} -(\boldsymbol{\tilde{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\tilde{u}} -\overline{\boldsymbol{\tilde{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\tilde{u}}}), \end{equation}

which is incompressible

(2.5)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tilde{u}} = 0, \end{equation}

where subscripts denote derivatives (e.g. $U_y:=\textrm {d}U/{\textrm {d}y}$).

2.1. Orr–Sommerfeld stability analysis

Given a possibly turbulent flow $(U,\boldsymbol {\tilde {u}})$, the ‘standard’ linear stability analysis is to consider small (also known as infinitesimal) perturbations $(0, \delta \boldsymbol {\tilde {u}})$ to a base state $(U,\boldsymbol {0})$ where the fluctuation field is ignored and the mean flow $U$ is assumed steady. As a result, only the (2.4) and (2.5) need be perturbed (and, hence, linearised) and since $U=U(y)$, the ensuing eigenvalue calculation is parameterised by a streamwise and spanwise wavenumber. Squire's theorem Squire (1933) is usually invoked to focus the search for instability to spanwise-independent perturbations and leads to the celebrated OS equation (Orr 1907; Sommerfeld 1908). In primitive variables, $\tilde {w}$ then decouples from $\tilde {u}$ and $\tilde {v}$ and can be ignored, leaving the reduced eigenproblem

(2.6)\begin{gather} -{\rm i}m \alpha c\, \tilde{u} = \frac{1}{Re}(\tilde{u}_{yy}-m^2\alpha^2 \tilde{u})-{\rm i}m \alpha \tilde{p} -{\rm i}m \alpha U\tilde{u}-U_y \tilde{v}, \end{gather}

(2.7)\begin{gather}-{\rm i}m \alpha c \,\tilde{v} = \frac{1}{Re}(\tilde{v}_{yy}-m^2 \alpha^2\tilde{v})-\tilde{p}_y -{\rm i}m \alpha U\tilde{v}, \end{gather}

(2.8)\begin{gather}0 = {\rm i}m \alpha \, \tilde{u}+\tilde{v}_y , \end{gather}

where $(\tilde {u},\tilde {v},\tilde {p}) \propto \exp ({\textrm {i}m \alpha (x-ct)})$ and $c:=c_r+ic_i$ is the (complex) eigenvalue. The OS equation is then reached by defining a streamfunction $\psi$ (so $\tilde {u}=\psi _y$ and $\tilde {v}=-\textrm {i}m \alpha \psi$) and eliminating the pressure so that

(2.9)\begin{equation} (U-c)(\partial_y^2-m^2 \alpha^2)\psi-U_{yy} \psi=\frac{1}{{\rm i} m \alpha Re} (\partial_y^2-m^2 \alpha^2)^2 \psi. \end{equation}

In what follows, we actually work with the primitive variable formulation (2.6)–(2.8) as it is clearer to interpret the origin of new terms added below, the eigenvalue problem is better conditioned (two equations with second-order operators as opposed to one with a fourth-order operator) and it is easily extended to three dimensions if needed (for three-dimensional (3-D) disturbances, the OS equation must be augmented with Squire's equation for the cross-stream vorticity). Nevertheless, we refer to this general approach of doing linear stability around a turbulent mean as ‘OS analysis’ in recognition of its conception in Malkus's work.

The (matrix) size of the eigenproblem (2.6)–(2.8) is only $3N_y \times 3N_y$ for each streamwise wavenumber $m$ and so needs to be repeated $N_x$ times ($N_x,N_y$ represent the streamwise and wall-normal truncations). When the mean profile is known from simulations or experiments, for typical resolutions, this is an easily accessible procedure that, through the structure of the most unstable eigenvectors, can shed some light on the dominant structures of the turbulent flow.

A popular extension to the standard OS analysis involves using an eddy viscosity $E(y)$ instead of the molecular viscosity. An eddy viscosity can be determined self-consistently from the Reynolds stress needed to sustain the turbulent mean profile (e.g. Reynolds & Tiederman 1967) using the steady version of (2.3),

(2.10)\begin{align} \frac{1}{Re} [(1+E(y))U_y]_y&:=\frac{1}{Re}U_{yy}- (\overline{\tilde{u} \tilde{v}} )_y={-}G \Rightarrow E(y) \nonumber\\ &:={-}Re \frac{ \overline{ \tilde{u} \tilde{v}} }{U_y} ={-}\frac{Re}{U_y} \int^y_0 G(\bar{y})\,{\rm d}\bar{y}-1. \end{align}

Then, the eigenproblem can be modified to

(2.11)\begin{gather} -{\rm i}m \alpha c\, \tilde{u} = \frac{1}{Re} \left( \left[(1+E)\tilde{u}_{y} \right]_y-m^2\alpha^2 (1+E) \tilde{u} \right)-{\rm i}m \alpha \tilde{p} -{\rm i}m \alpha U\tilde{u}-U_y \tilde{v}, \end{gather}

(2.12)\begin{gather}-{\rm i}m \alpha c \,\tilde{v} = \frac{1}{Re}\left( \left[(1+E)\tilde{v}_{y} \right]_y-m^2\alpha^2 (1+E) \tilde{v} \right)-\tilde{p}_y -{\rm i}m \alpha U\tilde{v}, \end{gather}

(2.13)\begin{gather}0 = {\rm i}m \alpha \, \tilde{u}+\tilde{v}_y. \end{gather}

We will use this to assess the performance of our EOS analysis in § 4.2.

2.2. Statistics: cumulants

In this section we consider a statistical framework for the flow by working with the equal-time cumulants of the flow (Hopf 1952; Orszag 1977; Frisch 1995). Even within this framework, we specialise further to exclusively consider equal-$x$-and-$z$ cumulants that are the subset of cumulants which influence the mean flow. (In fact, only equal-$y$ cumulants are needed in the mean flow equation – see (2.23) below. However, the evolution equations for these are not available without also solving for ‘non-equal’ $y$ cumulants.) The first cumulant is the mean $U(y,t)$. The second cumulant is the symmetric matrix

(2.14)\begin{equation} \boldsymbol{C}(y_1,y_2,t):= \overline{\boldsymbol{\tilde{u}}(x,y_1,z,t) \otimes \boldsymbol{\tilde{u}}(x,y_2,z,t)} =\left(\begin{array}{@{}ccc@{}} C_{11} & C_{12} & C_{13} \\ C_{21} & C_{22} & C_{23} \\ C_{31} & C_{32} & C_{33} \end{array} \right), \end{equation}

where we introduce the notation $C_{ij}(1,2):= \overline {[\boldsymbol {\tilde {u}}(x,y_1,z,t)]_i [\boldsymbol {\tilde {u}}(x,y_2,z,t)]_j}$ (here $[\boldsymbol {\tilde {u}}]_1=\tilde {u}$, $[\boldsymbol {\tilde {u}}]_2=\tilde {v}$ and $[\boldsymbol {\tilde {u}}]_3=\tilde {w}$) to emphasize the $y$ arguments and de-emphasize the implicit time dependence, and the third cumulant is the third-order tensor

(2.15)\begin{equation} C_{ijk}^{(3)}(1,2,3):=\overline{[\boldsymbol{\tilde{u}}(x,y_1,z,t)]_i [\boldsymbol{\tilde{u}}(x,y_2,z,t)]_j [\boldsymbol{\tilde{u}}(x,y_3,z,t)]_k}. \end{equation}

These correspond to the second and third central moments of the flow, respectively. Cumulants and central moments, however, diverge at fourth order and beyond, e.g.

(2.16)\begin{align} C_{ijkl}^{(4)}(1,2,3,4)&:= \overline{ [\boldsymbol{\tilde{u}}(x,y_1,z,t)]_i [\boldsymbol{\tilde{u}}(x,y_2,z,t)]_j [\boldsymbol{\tilde{u}}(x,y_3,z,t)]_k [\boldsymbol{\tilde{u}}(x,y_4,z,t)]_l } \end{align}

(2.17)\begin{align} & \quad -C_{ij}(1,2)C_{kl}(3,4) -C_{ik}(1,3)C_{jl}(2,4)-C_{il}(1,4)C_{jk}(2,3). \end{align}

We will not go this high in the cumulant expansion used here but just note that the $n$th-order cumulant is $n$-dimensional in space so that storage when doing computations becomes prohibitive very quickly. Hence, the onus is on applying some sort of closure as soon as possible.

To derive evolution equations for the cumulants, we introduce a double Fourier series representation of the flow

(2.18)\begin{equation} \boldsymbol{\tilde{u}}(x,y,z,t):= \sum_m \sum_n \boldsymbol{\tilde{u}}^{mn}(y,t) \exp({{\rm i}m \alpha x+{\rm i}n \beta z}), \end{equation}

where $m,n \in {\mathbb {Z}}$. Clearly $\boldsymbol {\tilde {u}}^{-m-n}=\boldsymbol {\tilde {u}}^{*mn}$ (the complex conjugate of $\boldsymbol {\tilde {u}}^{mn}$) for a real flow but it will be clearer not to build this into the notation in anticipation of deriving perturbation equations later. Hence, we write

(2.19)\begin{align} C_{ij}(1,2) &= \sum_m \sum_n \left\{ C_{ij}^{mn}(1,2):= [\boldsymbol{\tilde{u}}^{mn} (y_1,t)]_i [\boldsymbol{\tilde{u}}^{{-}m-n} (y_2,t)]_j \right\}, \end{align}

(2.20)\begin{align}C^{(3)}_{ijk}(1,2,3) &= \sum_m \sum_n \left\{ C_{ijk}^{(3)mn}(1,2,3) :=\sum_p \sum_q [\boldsymbol{\tilde{u}}^{mn}(y_1,t)]_i [\boldsymbol{\tilde{u}}^{pq}(y_2,t)]_j\right. \nonumber\\ &\quad \left. \times [\boldsymbol{\tilde{u}}^{-(m+p)-(n+q)}(y_3,t)]_k \vphantom{\sum_p}\right\} \end{align}

(note, e.g. $C_{ij}^{mn}(1,2)=C_{ji}^{-m-n}(2,1)$). Equations to evolve the cumulants are obtained by temporally differentiating their definitions in (2.19) and (2.20) and using (2.4). For example, for the second-order cumulant,

(2.21)\begin{align} &\partial_t C_{ij}^{mn}(1,2) = [\boldsymbol{\tilde{u}}^{mn}(y_1,t)]_i \partial_t [\boldsymbol{\tilde{u}}^{{-}m-n}(y_2,t)]_j +\partial_t [\boldsymbol{\tilde{u}}^{mn}(y_1,t)]_i [\boldsymbol{\tilde{u}}^{{-}m-n}(y_2,t)]_j \nonumber\\ &\quad = \frac{1}{Re} \left(\partial^2_1+\partial^2_2 -2m^2\alpha^2-2n^2\beta^2\right) C_{ij}^{mn}(1,2) \nonumber\\ & \qquad -\left[\begin{array}{@{}c@{}} -{\rm i}m\alpha \\ \partial_2 \\ -{\rm i}n\beta \end{array}\right]_j C_{i4}^{mn}(1,2) -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_1 \\ {\rm i}n\beta \end{array}\right]_i C^{mn}_{4j}(1,2) \nonumber\\ &\qquad +{\rm i}m\alpha [U(2)-U(1)]C_{ij}^{mn}(1,2) -U_y(2)C_{i2}^{mn}(1,2) \delta_{1j} \nonumber\\ &\qquad -U_y(1)C_{2j}^{mn}(1,2) \delta_{i1} -\left[\begin{array}{@{}c@{}} -{\rm i}m\alpha \\ \partial_2 \\ -{\rm i}n\beta \end{array}\right]_k C_{ijk}^{(3)mn}(1,2,2) -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_1 \\ {\rm i}n\beta \end{array}\right]_k C_{ikj}^{(3)-m-n}(2,1,1), \end{align}

where $\partial _i:= \partial _{y_i}$, $U(i)=U(y_i)$ for $i=1,2$ and $C_{i4}^{mn}(1,2):= [\boldsymbol {\tilde {u}}^{mn} (y_1,t)]_i \tilde {p}^{-m-n}(y_2,t)=:C_{4i}^{-m-n}(2,1)$ are three extra ‘velocity–pressure’ cumulants that get generated. Incompressibility conditions give the required three extra matrix constraints

(2.22)\begin{equation} \left[\begin{array}{@{}c@{}} {\rm i}m \alpha\\ \partial_1 \\ {\rm i}n \beta \end{array}\right]_i C_{ij}^{mn}(1,2)=0 \quad j \in \{1,2,3\} \end{equation}

along with the equation for the mean equation (2.3)

(2.23)\begin{equation} U_t -\frac{1}{Re}U_{yy} = G - \sum_m \sum_n \partial_y C_{12}^{mn}(y,y,t) \end{equation}

to close the system. The infamous closure problem of the Navier–Stokes equation is immediately evident here in that the evolution equation for the second-order cumulant depends on the third-order cumulant, a pattern that continues for higher-order cumulants so the system never closes. A popular (lowest) closure – commonly called CE2 – is to simply ignore the third-order cumulant which is equivalent to ignoring the fluctuation-fluctuation term (last bracketed term on the right-hand side of (2.4)). This is the QL approximation or sometimes referred to as the mean field theory (e.g. Vedenov et al. 1961; Herring 1963, 1964),

(2.24)\begin{gather} U_t = \frac{1}{Re}U_{yy}+G - ( \overline{\tilde{u} \tilde{v}} )_y, \end{gather}

(2.25)\begin{gather}\boldsymbol{\tilde{u}}_t = \frac{1}{Re} \nabla^2 \boldsymbol{\tilde{u}} -\boldsymbol{\nabla} \tilde{p} -{U} \boldsymbol{\tilde{u}}_x - \tilde{v} U_y \boldsymbol{\hat{x}}, \end{gather}

(2.26)\begin{gather}0 =\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tilde{u}} . \end{gather}

The defining feature of this approximation is that (2.25) is linear in $\boldsymbol {\tilde {u}}$ so that fluctuations with different wavenumbers are only coupled in the mean flow equation (2.24). This linearity also means that any fluctuation field (parametrised by streamwise and spanwise wavenumbers) can not be consistently in the stable manifold of $U$ as it varies with time. In particular, if $U$ is steady, only marginally stable fluctuation fields (typically with a temporal frequency) can be non-vanishing. Malkus argued for this model (and its marginal stability implications) on the basis that the fluctuation-fluctuation nonlinear term was only stabilising. This would be reasonable if bifurcations from unidirectional shear flows were always supercritical but, some decades later, subcriticality was realised the more generic situation (Kerswell 2005; Eckhardt, Schneider & Westerweel 2007; Graham & Floryan 2021; Kawahara, Uhlmann & van Veen 2012).

Applying a closure at next order so $C^{(4)}$ is some assumed function of the lower cumulants or simply ignored (termed CE3) is less straightforward as the ensuing positive definiteness of the second cumulant is not automatic (Marston, Qi & Tobias 2019). This difficulty explains the popularity of the lower-order CE2 approximation where, for example, the existence and stability of steady solutions has recently been investigated for ordinary differential equation (ODE) systems (Li, Marston & Tobias 2021; Li et al. 2022).

2.3. Approximations to statistical stability

The approach here is to consider the stability within the cumulant framework as this presents a natural way to assess stability of the flow statistics. Ideally, this should be attempted for a cumulant expansion which is high enough order to show a robustness against including even higher-order cumulants. However, the rate at which the dimensionality of this procedure explodes means that only second- and perhaps third-order closures are currently practical. As a result, we focus on CE2 here and identify a clear way to progress to higher order (see § 3.2).

Here CE2 is the statistical equations (2.21)–(2.23) with the third-order cumulants in (2.21) set to zero. The corresponding equations for perturbations $(\delta U, \delta C_{ij}^{mn} )$ upon a base statistical state $(U,C_{ij}^{mn})$ – hereafter referred to as the $\delta$CE2 problem – are

(2.27)\begin{align} \partial_t \,\delta C_{ij}^{mn}(1,2) &= \frac{1}{Re} \left(\partial^2_1+\partial^2_2 -2m^2\alpha^2-2n^2\beta^2\right) \delta C_{ij}^{mn}(1,2) \nonumber\\ &\qquad -\left[\begin{array}{@{}c@{}} -{\rm i}m \alpha\\ \partial_2 \\ -{\rm i}n \beta \end{array}\right]_j \delta C_{i4}^{mn}(1,2) -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_1 \\ {\rm i}n\beta \end{array}\right]_i \delta C^{mn}_{4j}(1,2) \nonumber\\ &\qquad +{\rm i}m \alpha [U(2)-U(1)] \delta C_{ij}^{mn}(1,2) -U_y(2) \delta C_{i2}^{mn}(1,2) \delta_{1j} \nonumber\\ &\qquad -U_y(1) \delta C_{2j}^{mn}(1,2) \delta_{i1} \nonumber\\ &\qquad +{\rm i}m \alpha [\delta U(2)- \delta U(1)] C_{ij}^{mn}(1,2) -\delta U_y(2) C_{i2}^{mn}(1,2) \delta_{1j}\nonumber\\ &\qquad -\delta U_y(1) C_{2j}^{mn}(1,2) \delta_{i1}, \end{align}

(2.28)\begin{align} 0 &= \left[\begin{array}{@{}c@{}} {\rm i}m\alpha\\ \partial_1 \\ {\rm i}n \beta \end{array}\right]_i \delta C_{ij}^{mn}(1,2), \end{align}

(2.29)\begin{align} \partial_t \delta U &= \frac{1}{Re}\delta U_{yy} +\delta G - \sum_m \sum_n \partial_y \delta C_{12}^{mn}(y,y,t). \end{align}

An equivalent equation arises in non-modal stability theory (e.g. Farrell & Ioannou 1993; Jovanovic & Bamieh 2005). If the pressure gradient is kept fixed $\delta G=0$, otherwise the constant volume flux condition $\int ^1_{-1} \delta U {\textrm {d}\kern 0.05em y}=0$ is the extra constraint required. Crucially the ansatz $(\delta U, \delta C_{ij}^{mn}) \propto \textrm {e}^{ \lambda t}$ is possible if the base state is independent of time – in other words, the base flow is statistically steady.

Even in the cheapest 2-D situation, the CE2 stability problem requires handling $5 N_x$ correlation matrices of size $N_y^2$ all linked through the mean equation. This leads to a matrix of size $(N_y+5 N_x N_y^2)^2$ or $\approx 25N_x^2 N_y^4$ elements which is out of reach even for modest resolutions. Given this, we explore a route to potentially still capture the essence of the CE2 statistical stability approximation without the considerable cost. To do this, we discuss a connection back to the equations of motion that plausibly retains the steadiness of the stability problem.

2.4. Steady statistics and simplications

A statistically steady base flow has steady mean $U$ and cumulants $\boldsymbol {C}$ with their Fourier components $\boldsymbol {C}^{mn}(1,2)$ also steady under ensemble averaging. This ensures that the associated stability problem (2.27)–(2.29) has temporally constant coefficients and is therefore a (conceptually at least) simple eigenvalue problem. In what follows below, we will not attempt to solve this directly as it is too unwieldy. Instead, a smaller, more practical QL stability problem is sought as a good proxy for it. The key in doing this is identifying a suitable base velocity field around which to develop a QL-type stability problem. A straightforward approach is to find the ‘best’ rank-1 approximation of each $\boldsymbol {C}^{mn}$ and use the associated velocity field, $\boldsymbol {\tilde {u}}^{mn}_0(y)$, as representative of the base flow. This can be accomplished by minimising the Frobenius matrix norm of the difference,

(2.30)\begin{align} &\|\boldsymbol{C}^{mn}-\boldsymbol{\tilde{u}}^{mn}_0 \otimes \boldsymbol{\tilde{u}}^{*mn}_0\|^2_F := \sum_{i,j=1}^3 \sum_{p,q=1}^N (C^{mn}_{ij}(p,q)-[\boldsymbol{\tilde{u}}^{mn}_0(y_p)]_i[\boldsymbol{\tilde{u}}^{*mn}_0(y_q)]_j) \nonumber\\ &\qquad \times (C^{*mn}_{ij}(p,q)-[\boldsymbol{\tilde{u}}^{*mn}_0(y_p)]_i[\boldsymbol{\tilde{u}}^{mn}_0(y_q)]_j) \nonumber\\ &\quad =\left\|\left(\begin{array}{@{}ccc@{}} C^{mn}_{11} & C^{mn}_{12} & C^{mn}_{13} \\ C^{mn}_{21} & C^{mn}_{22} & C^{mn}_{23} \\ C^{mn}_{31} & C^{mn}_{32} & C^{mn}_{33} \end{array} \right) - \left[\begin{array}{@{}c@{}} \tilde{u}^{mn}_{(1:N)}\\ \tilde{v}^{mn}_{(1:N)} \\ \tilde{w}^{mn}_{(1:N)} \end{array} \right] \left[\begin{array}{@{}ccc@{}} \tilde{u}^{*mn}_{(1:N)} & \tilde{v}^{*mn}_{(1:N)} & \tilde{w}^{*mn}_{(1:N)} \end{array} \right] \right\|^2_F, \end{align}

where, e.g. $\tilde {u}^{mn}_{(1:N)}:=[ \tilde {u}^{mn}(y_1)\ \tilde {u}^{mn}(y_2)\ \ldots \tilde {u}^{mn}(y_N)]^\textrm {T}$. Since $\boldsymbol {C}^{mn}=[\boldsymbol {C}^{mn}]^H$ (the Hermitian conjugate) and positive definite, the required $\boldsymbol {\tilde {u}}^{mn}_0$ is the leading right eigenvector of $\boldsymbol {C}^{mn}$ associated with the largest eigenvalue $\sigma _1$ scaled so that $|\boldsymbol {\tilde {u}}^{mn}_0|=\sqrt {\sigma _1}$.

We now discuss how the QL stability problem based upon the $\boldsymbol {\tilde {u}}^{mn}_0$ relates to the CE2 stability problem based on $\boldsymbol {C}^{mn}=\boldsymbol {\tilde {u}}^{mn}_0 \otimes \boldsymbol {\tilde {u}}^{*mn}_0$. Possible disturbances partition into two types: type A where the disturbance has energy in Fourier pairings not excited in the base flow, and type B where the disturbance has energy in Fourier pairings that are a subset of those present in the base flow, i.e. $\boldsymbol {\delta } \boldsymbol {\tilde {u}}^{mn}$ is only non-zero if $\boldsymbol {\tilde {u}}_0^{mn}$ is (see Appendix A for more detail). The former (type A) case is straightforward (since $\boldsymbol {\delta } \boldsymbol {C}=\boldsymbol {0}$) so we focus here on the type B situation. The QL stability problem – hereafter the $\delta$QL problem – is

(2.31)\begin{gather} \partial_t \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} = {\mathbb{L}}^{mn}(U) \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} +{\mathbb{J}}^{mn}(\delta U) \boldsymbol{\tilde{u}}_0^{mn} \quad \forall (m,n) \neq (0,0), \end{gather}

(2.32)\begin{gather}\partial_t \delta U = \frac{1}{Re}\delta U_{yy} +\delta G - \sum_m \sum_n \partial_y \left\{ \tilde{u}_0^{mn} \delta \tilde{v}^{{-}m-n} +\delta \tilde{u}^{mn} \tilde{v}_0^{{-}m-n} \right\}, \end{gather}

where

(2.33)\begin{gather} {\mathbb{L}}^{mn}(U)\, \boldsymbol{\tilde{u}}^{mn} := \frac{1}{Re} [\partial_y^2-m^2\alpha^2 -n^2 \beta^2]\boldsymbol{\tilde{u}}^{mn} -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_y \\ {\rm i}n\beta \end{array} \right] \tilde{p}^{mn} -{\rm i}m \alpha{U} \boldsymbol{\tilde{u}}^{mn} - \tilde{v}^{mn} U_y \boldsymbol{\hat{x}}, \end{gather}

(2.34)\begin{gather}{\mathbb{J}}^{mn}(\delta U) :=\lim_{\epsilon \rightarrow 0} \, \frac{1}{\epsilon}[{\mathbb{L}}^{mn}(U+ \epsilon\delta U)-{\mathbb{L}}^{mn}(U)]={\mathbb{L}}^{mn}(U+\delta U)-{\mathbb{L}}^{mn}(U) \end{gather}

(as ${\mathbb {L}}^{mn}$ is affine in $U$). This has temporally constant coefficients and, therefore, admits an eigenfunction of the form

(2.35)\begin{align} &(\delta U(y,t), \boldsymbol{\delta} \boldsymbol{\tilde{u}}(x,y,z,t)) \nonumber\\ &\quad =\left( \delta \hat{U}(y) {\rm e}^{\lambda t}, \,\sum_m \sum_n \left[ \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn}(y,t):=\boldsymbol{\delta} \widehat{\boldsymbol{u}}^{mn}(y) {\rm e}^{\lambda t} \right]\exp({{\rm i}(m \alpha x+ n \beta z)}) \right). \end{align}

The key point is that this has a direct equivalent in the CE2 problem (2.27)–(2.29) where the same eigenvalue $\lambda$ has the corresponding eigenfunction $\boldsymbol {C}^{mn}:= \boldsymbol {\tilde {u}}_0^{mn} (y_1,t)\otimes \boldsymbol {\delta } \boldsymbol {\tilde {u}}^{-m-n}(y_2,t)+ \boldsymbol {\delta } \boldsymbol {\tilde {u}}^{mn} (y_1,t)\otimes \boldsymbol {\tilde {u}}_0^{-m-n}(y_2,t) =\widehat {\boldsymbol {C}}^{mn}(y_1,y_2)\textrm {e}^{\lambda t}$. Given this, instability within the much more tractable $\delta$QL problem is therefore sufficient to conclude statistical instability within $\delta$CE2. Going further to claim that the stability within the $\delta$QL problem also implies stability within the $\delta$CE2 problem seems very likely but is not assured (see Appendix A). With this one caveat, we nevertheless assume that the stability of the (much) smaller QL system is a proxy for the statistical stability of the CE2 system: in particular, stability in the $\delta$QL system is taken to imply statistical stability within $\delta$CE2. This allows us to generate an extended version of OS, or EOS analysis but, before pursuing this in the next section, we make a remark and an observation.

The remark is that certain notational liberties have been taken here to keep the discussion as clear as possible. For example, the operator ${\mathbb {L}}^{mn}$ strictly maps an incompressible flow field to another that involves a supplementary scalar field (the pressure) entering into the definition. As is well known, this is determined by imposing incompressibility. An implicitly incompressible representation for the velocity field could be used (e.g. reducing the problem down to just using wall-normal velocity and vorticity) to avoid this wrinkle, but staying with primitive variables makes the various manipulations as clear as possible.

The observation is that the special form of the nonlinearity in the QL and CE2 formulations means that different wavenumber pairings only interact through the mean flow equation. As a result, eigenvalues can be sought within any subset of the wavenumbers possible and these are still valid within the full system of wavenumbers, i.e. this is not a truncation merely a subclass of disturbances. In particular, for CE2, the simplest perturbation of the base state $(U,\boldsymbol {C})$ just consists of perturbations in one wavenumber pairing and the mean flow,

(2.36)\begin{equation} (\delta \hat{U}, \,\boldsymbol{\delta} \widehat{\boldsymbol{C}}^{mn}) \,{\rm e}^{\lambda t}. \end{equation}

Even this leads to an eigenvalue matrix calculation of size $(N_y+9N_y^2) \times (N_y+9N_y^2)$ in three dimensions. The equivalent QL perturbation is

(2.37)\begin{equation} (\delta \hat{U}, \,\boldsymbol{\delta} \boldsymbol{\hat{u}}^{mn},\,\boldsymbol{\delta} \boldsymbol{\hat{u}}^{{-}m-n}) \,{\rm e}^{\lambda t} \end{equation}

and requires an eigenvalue matrix calculation of size $(N_y+8N_y) \times (N_y+8N_y)$ that is much smaller (or $(N_y+6N_y)\times (N_y+6N_y)$ in two dimensions – see Appendix B).

3. Extended OS analysis

We refer to the $\delta$QL problem that includes some information on the second-order flow statistics of the base flow as EOS stability analysis. The approach is as follows.

1. (i) Estimate the base flow statistics $(U, \boldsymbol {C})$.

2. (ii) Approximate each base cumulant component tensor $\boldsymbol {C}^{mn}$ as rank 1 to obtain a representative physical base field $\boldsymbol {\tilde {u}}^{mn}_0(y)$.

3. (iii) Solve the EOS eigenvalue problem for small perturbations $(\delta U, \,\boldsymbol {\delta } \boldsymbol {\tilde {u}}, \,\delta \tilde {p})$ that is

   (3.1)\begin{equation} \partial_t\delta U=\frac{1}{Re}\partial^2_y \delta U+\delta G - \sum_{m,n} \partial_y \left(\delta \tilde{u}^{mn} \tilde{v}^{{-}m-n}_0 +\tilde{u}_0^{mn}\delta \tilde{v}^{{-}m-n} \right) \end{equation}
   coupled with the fluctuation equation for every Fourier mode $(m,n)$,
   (3.2)\begin{gather} \partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{mn} = \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn}-\boldsymbol{\nabla} \delta \tilde{p}^{mn} - {\rm i}m U \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} - \delta \tilde{v}^{mn}U_y\boldsymbol{\hat{x}} -{\rm i}m{\delta U} \boldsymbol{\tilde{u}}_0^{mn} -\tilde{v}_0^{mn} \delta U_y \boldsymbol{\hat{x}}, \end{gather}
   (3.3)\begin{gather}{\rm i}m \alpha\, \delta\tilde{u}^{mn}+\partial_y\, \delta \tilde{v}^{mn}+{\rm i}n\beta\, \delta \tilde{w}^{mn} = 0 \end{gather}
   (In (3.2), the last two terms depend on $\delta U$ and will be discussed later).

This eigenvalue problem has size $(N_y+4 N_x N_y N_z)^2 \approx 16N_x^2N_y^2N_z^2$ that is probably impractical for all but the smallest systems since all the wavenumber pairings are coupled through the mean equation. A natural way to simplify the calculation is to only include a targeted subset of the wavenumber pairings and, going further, to only consider one wavenumber pairing. This latter approximation is obviously the most extreme but is also closest in spirit to the original OS analysis – we call this minimally extended OS, or mEOS, analysis.

3.1. Minimally EOS equations

Here, the sum over wavenumber pairings in (3.1) is removed leaving the mean flow equation (3.1) coupled with just one Fourier mode equation (3.2). This eigenvalue problem focuses on the coupling between an individual Fourier mode and the mean flow and its size is now $(5 \times N_y)^2 = 25 N_y^2$, so comparable to OS analysis. Just as for OS analysis, it needs to be repeated for each Fourier wavenumber of interest.

3.2. Infinitely EOS equations

Before going on to test EOS and mEOS analyses on some statistically steady flows, it is worth briefly discussing how far this approach of using the cumulant framework to further motivate more sophisticated versions of OS analysis can go. In fact, it turns out that only one further enhancement is possible and then only for CE$\infty$. This is because any intermediate closure will lead to an equation that when ‘unwrapped’ at the highest level contradicts those obtained at lower cumulant orders. To see this, consider CE3, the next closure after CE2 in which $\boldsymbol {C}^{(4)}$ is ignored. Unwrapping the statistical equation for $\boldsymbol {C}$ will recover the Navier–Stokes equations whereas unravelling the $\boldsymbol {C}^{(3)}$ equations will not as the nonlinear term leading to $\boldsymbol {C}^{(4)}$ has been dropped (any other closure will also suffer this inconsistency). The only way to avoid this is: (1) to only unwrap one cumulant equation – so CE2 which leads to EOS, or (2) ensure that all unwrapped equations are the same, which means no closure and CE$\infty$. In this latter case, the Navier–Stokes equations re-emerge that, under perturbation around a steady base state, lead to the linearized Navier–Stokes equations as the ultimate extension of OS analysis, i.e.

(3.4)\begin{equation} \partial_t\delta U=\frac{1}{Re}\partial^2_y \delta U+\delta G - \sum_{m,n} \partial_y \left(\delta \tilde{u}^{mn} \tilde{v}^{{-}m-n}_0 +\tilde{u}_0^{mn}\delta \tilde{v}^{{-}m-n} \right) \end{equation}

and, for each wavenumber pair $(m,n)$,

(3.5)\begin{gather} \partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{mn} = \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn}-\boldsymbol{\nabla} \delta \tilde{p}^{mn} - {\rm i}m \alpha U \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} - \delta \tilde{v}^{mn}U_y\boldsymbol{\hat{x}} -{\rm i}m \alpha {\delta U} \boldsymbol{\tilde{u}}_0^{mn} -\tilde{v}_0^{mn} \delta U_y \boldsymbol{\hat{x}} \nonumber\\ -[\boldsymbol{\delta} \boldsymbol{\tilde{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\tilde{u}}_0+\boldsymbol{\tilde{u}}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\delta} \boldsymbol{\tilde{u}}]^{mn}, \nonumber\\ {\rm i}m \alpha \, \delta \tilde{u}^{mn}+\partial_y \delta \tilde{v}^{mn}+{\rm i}n \beta \, \delta \tilde{w}^{mn} = 0. \end{gather}

This could be called an iEOS analysis as there is no further extension possible with the cumulant expansion framework considered here. It is worth remarking that the extra term added in the iEOS analysis would get dropped in any minimal version where the individual wavenumber pairings are considered separately, i.e. minimalizing iEOS is just mEOS. There is, however, considerable potential in avoiding this drastic reduction in favour of retaining small subsets of interacting wavenumber pairings such as triads that interact through the newly present term. This or iEOS will not be pursued further here to avoid overcomplicating the discussion.

4. Application

We now test the extended stability approaches, EOS and mEOS, upon turbulent states obtained by DNS of 2-D channel flow. At high enough $Re$, the expectation is that the flow will be statistically stationary (Malkus 1956). The simulations are performed using the open-source PDE solver Dedalus (Burns et al. 2020). The 2-D flow is simulated using a vorticity–streamfunction formulation so that $(u,v):=(\psi _y,-\psi _x)$ and

(4.1)\begin{equation} \omega_t + \psi_y \omega_x - \psi_x \omega_y = \frac{1}{Re} \nabla^2 \omega - \partial_y f(y,t), \quad \omega:={-} \nabla^2 \psi, \end{equation}

where the flow is driven by a streamwise body force $f(y,t)\boldsymbol {\hat {x}}$ so that the volume flux is fixed. Here $Re$ is defined as in § 2.1 so that the conditions

(4.2)\begin{equation} \psi\left({-}1\right)=0, \quad \psi\left(1\right)=\tfrac{4}{3}, \quad \psi_y\left({-}1\right)=0, \quad \psi_y\left(1\right)=0. \end{equation}

are imposed with the latter two reflecting the presence of non-slip walls. Following earlier work (Falkovich & Vladimirova 2018), the length of the channel $2L_x$ is set to four times its height ($L_x=4$), $Re=36\,300$ and, computationally, $1024$ Fourier modes are used to discretize in the $x$ direction and $256$ Chebyshev modes in $y$ direction (see Markeviciute & Kerswell (2021) for details).

The streamwise body force is defined in terms of a profile function $g(y)$ as

(4.3)\begin{equation} f(y,t):=G(t)\left(g(y)-1\right), \end{equation}

where setting $g=0$ recovers the usual $y$-independent applied pressure gradient $-G(t)$. For this situation, it is already known that there is bistability with two statistically steady states possible in a 2-D channel flow at $Re=36\,300$ (Markeviciute & Kerswell 2021): a state that is statistically symmetric about the channel midplane – the symmetric state S – and another that is statistically asymmetric – the asymmetric state A; see figure 2(a). Choosing non-zero $g(y)$ is a device to diversify the test states available with two extra examples generated by using profiles $g_1:=(1-y^2)^6$ (state F1) and $g_2:=\cos \tfrac {3}{2} {\rm \pi}y^2$ (state F2).

Figure 2. Mean velocity profiles (solid lines) and streamwise root-mean-squared velocity profiles (dashed lines) for the four 2-D channel test states at $Re=36\,300$. (a) States with an applied pressure gradient: symmetric S (light blue) and asymmetric A (pink). (b) Body-forced states F1 (blue) and F2 (red). Black solid line shows the laminar parabolic profile in both plots for reference.

4.1. Statistically stable test states

A spatial average over the streamwise direction and temporal averaging was used as a proxy for the ensemble-averaged statistics discussed in deriving EOS equations (these are equivalent for a statistically steady flow over a large enough domain and long enough time). The mean profiles for states S, A, F1 and F2 are shown in figure 2 along with their streamwise root-mean-squared velocity profiles, all obtained by time averaging over $10^4$ time units. Their respective power spectra are shown in figure 3(a). While the dominant wavenumber is the same for all four states ($k_d=2\alpha ={\rm \pi} /2$), significant differences are seen at neighbouring streamwise wavenumbers $2k_d$ and $3k_d$ for example. Typical temporal variations of the fluctuation field energy $E_{fluct}$ compared with the total energy of the flow $E$ – figure 3(b) – show the desired wide variety of mean fluctuation energy and fluctuation amplitudes. In particular, the amplitude of the fluctuations in state F2 are $\approx 60\,\%$ of the mean and only $\approx 30\,\%$ for state A.

Figure 3. Fluctuation energy as a ratio to the total flow energy for the four 2-D channel test states at $Re=36\,300$: symmetric S (light blue), asymmetric A (pink), F1 (blue) and F2 (red). (a) Decomposition by streamwise wavenumber $k$. (b) Temporal fluctuations over a typical time-averaging window.

For each state, the second-order cumulant matrices were computed every unit of time and then averaged over a period of $10^3$. During the process of calculating EOS and mEOS eigenvalues, correlation matrices time averaged over different time windows, $t = 500, 1000, 1500$, were used as well as correlation matrices time averaged with different time steps (0.5 and 1) for one time window $t=500$. The qualitative results were all the same, with only minor quantitative differences that did not affect the relative positions of eigenvalues obtained by different methods. The time-averaged cumulant matrix was then diagonalised and the flow field corresponding to the leading (real) eigenvalue used to typify the base fluctuation field. The implications of the time averaging to the rank of the matrix are discussed in § 5.1.

The EOS problem – (3.1) and (3.2) – and the mEOS problem were posed as generalized matrix eigenvalue problems of matrix size $\approx ( 4 N_x N_y)^2$ and $(5N_y)^2$, respectively. For EOS analysis, $N_x=20$ streamwise wavenumbers were used as a balance between accurate eigenvalue convergence and computational accessibility; see, e.g. figure 4. Wall-normal resolution of the eigenvalue problem was also varied to ensure eigenvalue convergence in the mEOS case. Doubling the wall-normal resolution from the DNS resolution of $256$ Chebyshev modes to $512$ produced less than a $10^{-8}$ relative error in the eigenvalues, so all stability calculations were done using the DNS wall-normal resolution.

Figure 4. Typical eigenvalue convergence with respect to the total number of the streamwise wavenumbers $N_x$ included in the EOS model, shown for $k=2$ (a) and $k=3$ (b) for the F2 state.

4.2. Comparisons

Here we compare the eigenvalues obtained by standard (OS), extended (EOS) and mEOS stability analyses on the four different turbulent states. Since all our test cases are presumed statistically steady, for a model to be good at predicting statistical stability of the turbulent state, we expect all the eigenvalues to be stable, i.e. to have a negative real part. For the symmetric (figure 5a) and asymmetric (figure 5b) states, we observe no significant difference in the leading eigenvalues between the three stability models. While all three models predict the asymmetric state to be statistically stable, they also all predict statistical instability for the symmetric state. For these test states, the leading eigenvalue is associated with the first streamwise wavenumber (marked blue in the plot). We observe some changes for the other eigenvalues which we do not examine further as they are all stable and do not affect characterization of the statistical stability of the turbulent states.

Figure 5. Comparison of the eigenvalues $\sigma =-\textrm {i}m\alpha (c_r+\textrm {i}\,c_i)$ – so growth rate is $m \alpha c_i$ – for OS ($\bullet$), OS with eddy viscosity ($\square$), EOS ($\Diamond$) and mEOS ($\bigcirc$) analyses. Only the leading eigenvalues for $m=1$ (blue), $m=2$ (red) and $m=3$ (green) are shown. (a) Symmetric S, (b) asymmetric A, (c) F1, (d) F2 states.

Since all streamwise wavenumbers are coupled in the extended model, the leading wavenumber for the eigenvalue is assigned by examining the power spectrum of the corresponding eigenvector. Example power spectra for leading eigenvalues are shown in figure 6 where it is apparent that the mean velocity component of an eigenvector is at least an order of magnitude smaller than the leading wavenumber contribution.

Figure 6. Power spectra for the EOS eigenvectors for the leading eigenvalues. (a) Symmetric S state ($\times$) and asymmetric A state ($\bullet$, red) with $m=1$; (b) states F1 ($\blacksquare$, cyan) and F2 ($+$, blue) with $m=2$. The mean velocity perturbation ($m=0$) is at least an order of magnitude smaller than the leading wavenumber contribution to the power spectrum in all 4 cases.

The leading eigenvalue for the body-forced states states F1 (figure 5c) and F2 (figure 5d) is associated with the second streamwise wavenumber (marked red in the plot). For both F1 and F2 states, this eigenvalue is unstable in the standard stability analysis but becomes stable in EOS analysis. The stabilisation effect is also captured by mEOS analysis but is not so strong. While the statistical stability of the F1 state is predicted correctly by the extended model, the improvement is only partial for the F2 state. In addition to the unstable leading eigenvalue that is stabilised by the extended model, there are two subsequent unstable eigenvalues corresponding to the first and third streamwise wavenumbers (marked blue and green on the plot). Similar to the symmetric state case, these two unstable eigenvalues see no improvement towards stability when the extended stability models are applied.

We also repeat the standard OS analysis with eddy viscosity as described in § 2.1 (marked using squares in figure 5). While eddy viscosity has a general stabilising effect on the OS eigenvalues, this effect is rarely enough to stabilise the unstable OS eigenvalues. Significant improvement is only seen for the symmetric S state (top left plot) where the unstable eigenvalue corresponding to $m=1$ is stabilised, while the unstable eigenvalue corresponding to $m=3$ for the F2 state (bottom right plot) moves to just ${O} (10^{-4})$ under the instability line. However, the most unstable eigenvalues corresponding to $m=2$ for the states F1 and F2 where significant improvement is seen using the EOS analyses, are not notably affected by the eddy viscosity. Note, to avoid the singularity in the expression for eddy viscosity (2.10) for the asymmetric state, the mean velocity profile $U$ was symmetrised by $\tfrac {1}{2}(U(y)+U(-y))$ and the Reynolds stress $\overline { \tilde {u} \tilde {v}}$ was anti-symmetrised by $\tfrac {1}{2}( \overline { \tilde {u} \tilde {v}}(y) \unicode{x2013} \overline { \tilde {u} \tilde {v}}(-y) )$.

So, in summary, for each unstable eigenvalue observed in our test cases, the standard and extended models either agree on their statistical stability prediction or the extended model shifts the prediction towards statistical stability. The stabilising effect is qualitatively captured by both EOS and mEOS analyses with the mean velocity components making a non-zero power contribution to the eigenvectors of the extended model. We now examine why the extended stability models sometimes perform better than OS analysis and sometimes make very little improvement.

4.3. Eigenvalue perturbation analysis

In an attempt to rationalise the improvement (stabilisation) or not in mEOS of the leading eigenvalues, we carry out an eigenvalue perturbation analysis. While this neglects the interactions between different streamwise wavenumbers, it includes the minimal improvement that our models offer: adding the mean velocity perturbation into the stability consideration. This leads to an additional equation that governs the evolution of the mean velocity perturbation, and two extra terms – advection and shear – which affect the fluctuation perturbations. While formally limited to only small effects, the perturbation analysis provides a framework to study the structures of the mean velocity perturbations, base fluctuations and fluctuation perturbations and how their interactions affect the eigenvalues.

The unperturbed mEOS eigenvalue problem is the OS problem ((3.2) for one Fourier pairing without the $\delta U$ dependent terms) and (3.1) (we are imagining a small number $\varepsilon$ inserted in front of the $\delta U$ dependent terms and developing a perturbation expansion in this but actually $\varepsilon =1$). The latter is passive: it sets $\delta U$ but there is no feedback to the OS equation. Discretized, it reads

(4.4)\begin{equation} A \boldsymbol{y_R} = \sigma B \boldsymbol{y_R} \end{equation}

with matrices $A,B$, eigenvalue $\sigma$ and right eigenvector $\boldsymbol {y_R}:=[\delta U\ \boldsymbol {\delta } \boldsymbol {\tilde {u}}^{m0}]^\textrm {T}$. The the last two $\delta U$ terms in (3.2) are then treated as perturbations of $A$ that causes $(\sigma,\boldsymbol {y_R}) \rightarrow (\sigma +\delta \sigma, \,\boldsymbol {y_R}+\boldsymbol {\delta y_R})$ so that

(4.5)\begin{equation} \left(A + {\delta A_A} +\delta A_S \right)\left(\boldsymbol{y_R} + \boldsymbol{\delta y_R } \right) = \left( \sigma +\delta \sigma \right) B \left( \boldsymbol{y_R} +\boldsymbol{\delta y_R} \right), \end{equation}

where ${\delta A_A\boldsymbol {y_R}}$ corresponds to the advection term ${{-\textrm {i}m\alpha \boldsymbol {\tilde {u}}_0^{m0}\delta U_R }}$ and $\delta A_S \boldsymbol{y_R}$ to the shear term $-\tilde {v}_0^{m0} \partial _y \delta U_R \boldsymbol {\hat {x}}$. By taking the inner product using the corresponding left eigenvector of the unperturbed system $\boldsymbol {y_L}$ (${\boldsymbol y_L}^{{\dagger} } A=\sigma {\boldsymbol y_L}^{{\dagger} } B$), and considering first-order terms only, an expression then follows for the first-order perturbation to the eigenvalue $\delta \sigma$,

(4.6)\begin{equation} {\delta \sigma} = {\delta \sigma_A} + {\delta \sigma_S}:= \frac{\boldsymbol{y_L}^{{\dagger}}{\delta A_A\boldsymbol{y_R}}}{\boldsymbol{y_L}^{{\dagger}}B\boldsymbol{y_R}}+\frac{\boldsymbol{y_L}^{{\dagger}}{\delta A_S\boldsymbol{y_R}}}{\boldsymbol{y_L}^{{\dagger}}B\boldsymbol{y_R}}. \end{equation}

We present the two most interesting cases of the eigenvalue perturbation analysis. In the first case, we perform eigenvalue perturbation analysis on the leading eigenvalue of the standard OS analysis for the F2 state, which is seen to stabilise in mEOS and EOS analysis. The second case corresponds to the leading eigenvalue of the symmetric state, which shows no significant difference between the standard and extended approaches. The eigenvalues as well as perturbation analysis predictions are summarised in table 1. The left and right eigenvectors and perturbed terms are shown in figures 7 and 8.

Figure 7. Perturbation analysis results for F2 state. (a–c) Streamwise component of the unperturbed left eigenvector $\delta u_{sL}$, streamwise component of the perturbed advection term $-\textrm {i}m\alpha u_{0s}^{m0} \delta U_{sR}$ and streamwise component of the perturbed shear term $-v_{0s}^{m0} \partial _y \delta U_{sR}$. (d–f) Wall-normal component of the unperturbed left eigenvector $\delta v_{sL}$, wall-normal component of the perturbed advection term $-\textrm {i}m\alpha v_{0s}^{m0} \delta U_R$ and mean profile component of the unperturbed right eigenvector $-\delta U_{sR}$. Solid (dashed) lines indicate real (imaginary) parts of the vectors. Subscript $(\boldsymbol {\cdot })_s$ refers to the sine components of the eigenfunctions as explained in Appendix B.

Figure 8. Perturbation analysis results for S state. (a–c) Streamwise component of the unperturbed left eigenvector $\delta u_{cL}$, streamwise component of the perturbed advection term $-\textrm {i}m\alpha u_{0c}^{m0} \delta U_{cR}$ and streamwise component of the perturbed shear term $-v_{0c}^{m0} \partial _y \delta U_{cR}$. (d–f) Wall-normal component of the unperturbed left eigenvector $\delta v_{cL}$, wall-normal component of the perturbed advection term $-\textrm {i}m\alpha v_{0c}^{m0} \delta U_R$ and mean profile component of the unperturbed right eigenvector $-\delta U_{cR}$. Solid (dashed) lines indicate real (imaginary) parts of the vectors. Subscript $(\boldsymbol {\cdot })_c$ refers to the cosine components of the eigenfunctions as explained in Appendix B.

Table 1. Perturbation analysis results for the leading eigenvalue in the F2 state (left) and the symmetric S state (right). The first-order change in the eigenvalue caused by the advection ($\sigma_A$) and shear ($\sigma_S$) terms in the mEOS equations. Unperturbed (standard OS) eigenvalues are given by $\sigma _{OS}$ and mEOS eigenvalues are given by $\sigma _{mEOS}$. The table shows that F2 changes to being stable when using mEOS analysis, whereas S remains unstable.

For the body-force-driven F2 case, the advection term has a small and positive contribution to the real part of $\delta \sigma$ while the shear term has a large negative contribution to the real part of $\delta \sigma$, as summarised in table 1. This latter prediction is only qualitatively correct as the shear term is not a small perturbation giving a much larger predicted decay rate of $-0.169$ than the actual value of $-0.023$ obtained from mEOS analysis. The perturbation analysis, however, indicates that the shear in the perturbed mean field is causing the stabilisation. Looking in more detail at the fields in figure 7, it is clear that the active regions of the left and right eigenvectors coincide. In particular, the shear in $\delta U_R$ coincides with where the left eigenvector is significant and is of opposite phase giving the large stabilising effect. In contrast, the resulting advection term is not only an order of magnitude smaller than the shear term, but also has the same sign contribution as the left eigenvector, yielding a small destabilising contribution to the eigenvalue.

For the symmetric S case, the perturbation analysis shows both advection and shear terms make small destabilising contributions to the eigenvalue commensurate with mEOS analysis although there is still an error due to the finiteness of the perturbation. Examining the structure of the left and right eigenvectors shown in figure 8 indicates that they are ill-matched. The left eigenvector is concentrated much closer to the channel walls than the right eigenvector. Moreover, even though the right eigenvector component $\delta U$ has some structure near the wall, the interaction of this component with the perturbed operators moves this structure towards the centre of the wall for the shear term, avoiding any significant interaction with the left eigenvector. The advection term has some overlap with the left eigenvector, but it is also an order of magnitude smaller than the shear term, overall resulting in small destabilising contributions. It is also worth remarking that even at a resolution of $2048$ wall-normal Chebyshev modes, the eigenvectors possess the same small-scale structure seen in the figures.

The ‘take home’ message from this section is that the mean velocity perturbation, even if an order of magnitude smaller than the other components of the eigenvector (see figure 6), can interact with the base flow fluctuation fields in a way that produces a strong stabilising influence.

5. Model limitations

In this section we discuss the various approximations made in deriving EOS and mEOS which include time averaging the statistics, approximating correlation matrices as rank 1 and neglecting higher-order statistics. As already mentioned, the channel flow states considered above were assumed statistically steady. Making this assumption is required to apply Malkus's OS analysis in the first place and so subsequently EOS and mEOS analysis, but is only an approximation. Time averaging the statistics in a small computational box is a commonly used approach to get a better estimate of the steady statistics hypothetically generated in an infinite computational box.

5.1. Rank 1 approximation of correlation matrices

Time averaging the correlation matrix increases its rank from 1 up to potentially its full dimension; see figure 9 for singular values $\sigma _i$ of correlation matrices time averaged over 1000 time units. This is not necessarily a reflection that the statistics are non-stationary as a fluctuation field with two frequencies will give a rank-2 correlation matrix under averaging. For $m=2$, which is the dominant streamwise wavenumber in all of the test states, a gap between the leading and the second singular values is larger than an order of magnitude indicating that a rank-1 approximation of the time-averaged correlation matrix may be reasonable. On the other hand, time-averaged correlation matrices for $m \in \{1,3\}$ do not show a significant gap between the leading and second singular values. In this case, it is much more likely that some important fluctuation field features might not be accurately represented by the rank-1 representation of the correlation matrix and, as a result, the approach is less justified.

Figure 9. Leading singular values $\sigma _i$ of the correlation matrices $\boldsymbol {C}^{m0}$ time averaged over 1000 time units. Shown for $m=1$ (blue), $m=2$ (red) and $m=3$ (green). (a) Symmetric S, (b) asymmetric A, (c) F1, (d) F2 states. Note the gap after the leading singular value for $m=2$ (red) across all states confirming rank 1 approximation of the correlation matrix for the cases when mEOS and EOS analyses show stabilisation.

In an attempt to explore a more accurate representation of $\boldsymbol {C}^{mn}=\sum _{i=1} \sigma _i \,\boldsymbol {\hat {e}_i} \otimes \, \boldsymbol {\hat {e}_i}$ (where $\boldsymbol {\hat {e}_i}$ is the $i$th normalised eigenvector and $\sigma _i$ the largest eigenvalue or singular value of $\boldsymbol {C}^{mn}$ as it is Hermitian) beyond just the rank-1 approximation $\boldsymbol {C}^{mn} \approx \sigma _1 \,\boldsymbol {\hat {e}_1} \otimes \, \boldsymbol {\hat {e}_1}$ ($\boldsymbol {\tilde {u}}^{mn}_0=\sqrt {\sigma _1} \boldsymbol {\hat {e}_1}$), a weighted average of the leading $N$ fields,

(5.1)\begin{equation} \boldsymbol{\tilde{u}}^{mn}_0=\sum_{i=1}^{N} \frac{\sigma_i}{\sum_{i=1}^{N} \sigma_i} \,\sqrt{\sigma_i} \boldsymbol{\hat{e}_i}, \end{equation}

was also considered for the F2 state where both $m=1$ and $m=2$ wavenumbers have a non-rank-1 time-averaged correlation matrix $C^{m0}$. However, including $N=1,2,5$ or even $10$ eigenvectors into the expansion showed no significant difference in the leading eigenvalues (less than $1\,\%$ difference in the absolute value of the eigenvalue). Understanding the effect of a rank-1 approximation really requires time stepping the full statistical equations which is beyond the scope of this work.

5.2. Time-averaged statistics

To understand the effect of time averaging the statistics in the EOS stability models, numerical experiments were performed in which the perturbations to the turbulent base flow were exposed to a time-varying base flow. This experiment was first performed for the symmetric (S) state where the unstable OS eigenvalue remains unstable in EOS and mEOS analyses. From the DNS data, 10 turbulent flow histories were taken of length $\Delta t=500$ separated by at least 1000 time units to avoid statistical dependence. Each was taken as a turbulent, time-dependent base flow. Two random perturbations with energy $E_{pert} = 10^{-5} E_{base}$ of the base flow were then added to each of the 10 base flow histories (making a total of 20 cases) and then their evolution computed over a $\Delta t=500$ time window using the QL equations. For $(\delta U, \boldsymbol {\delta }\boldsymbol {\tilde {u}}^{m0})$, these read

(5.2)\begin{gather} \partial_t\delta U=\frac{1}{Re}\partial^2_y \delta U+\delta G -[\boldsymbol{\tilde{u}}_0(t).\boldsymbol{\nabla} \boldsymbol{\delta} \boldsymbol{\tilde{u}}+ \boldsymbol{\delta}\boldsymbol{\tilde{u}}.\boldsymbol{\nabla} \boldsymbol{\tilde{u}}_0(t)]^{00}, \end{gather}

(5.3)\begin{align}\partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{m0} &= \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0}-\boldsymbol{\nabla} \delta \tilde{p}^{m0} - {\rm i}m U \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0} - \delta \tilde{v}^{m0}U_y\boldsymbol{\hat{x}}\nonumber\\ &\quad -{\rm i}m{\delta U} \boldsymbol{\tilde{u}}_0^{m0} -\tilde{v}_0^{m0} \delta U_y \boldsymbol{\hat{x}}, \quad m>0, \end{align}

so that the perturbation-base flow fluctuation interactions are not included in the perturbation fluctuation equation (compare with (5.4) below where they are retained). Averaging the perturbation energy growth over the 20 different simulations revealed an unstable mode with exponential growth after initial transients decayed. The growth rate of this unstable mode was found to be $\sigma =0.00263$ which is within $5\,\%$ of the growth rate predicted by the EOS analysis; see figure 10. The most unstable eigenvectors of the EOS and mEOS analyses agree and are observed as the growing structure in this QL experiment; see figure 11. We conclude that time dependence of the base flow is not enough to recover the statistical stability of the symmetric state S.

Figure 10. Results of various numerical experiments: time evolution of the perturbation energy $E_{pert}$ as a ratio to the base flow energy $E_{base}$ averaged over multiple DNS runs. Perturbations were time stepped using linearised QL (blue), LNS (orange) or full nonlinear Navier–Stokes (green) equations. For each case, base flow was time dependent and generated using full Navier–Stokes DNS. Results shown for the symmetric S (a) and body-force-driven BF2 (b) states. Black dashed lines are provided as guides to emphasize the zero gradient of the saturated states and to indicate the growth rate ${\rm Re}(\sigma )$ of the unstable states. For comparison, the real part of the leading mEOS eigenvalue ${\rm Re}(\sigma _{mEOS})$ is given for each of the states.

Figure 11. Comparison of the leading eigenvectors from the EOS (a), mEOS $m=1$ (b) analyses and an instantaneous snapshot of the velocity perturbation field in the QL experiment at $t=100$ (c) for the symmetric (S) state. Only the streamwise velocity component is shown. There is a good correspondence in the leading eigenvector across EOS and mEOS analyses and the QL experiment using a time-dependent base flow.

5.3. Importance of higher-order statistics

To evaluate the need of flow statistics beyond the second order to recover statistical stability of the symmetric (S) state, we repeated the numerical experiment described above but this time the perturbation fields were time stepped using the linearised Navier Stokes (LNS) equations where perturbation-base flow fluctuation interactions are now included and the base flow is time varying,

(5.4)\begin{gather} \partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{m0} = \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0}-\boldsymbol{\nabla} \delta \tilde{p}^{m0} - {\rm i}m U(t) \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0} - \delta \tilde{v}^{m0}U_y(t)\boldsymbol{\hat{x}} -{\rm i}m{\delta U} \boldsymbol{\tilde{u}}_0^{m0}(t)\nonumber\\ -\tilde{v}_0^{m0}(t) \delta U_y \boldsymbol{\hat{x}} -[\boldsymbol{\tilde{u}}_0(t).\boldsymbol{\nabla} \boldsymbol{\delta} \boldsymbol{\tilde{u}}+ \boldsymbol{\delta}\boldsymbol{\tilde{u}}.\boldsymbol{\nabla} \boldsymbol{\tilde{u}}_0(t)]^{m0}. \end{gather}

We repeat the experiment with initial $E_{pert} = 10^{-5} E_{base}$ and $E_{pert} = 10^{-7} E_{base}$ obtaining the same qualitative results. When the energy growth of the perturbation field is time averaged over 20 different runs, the perturbation growth is seen to saturate after an initial transient growth; see figure 10(a). This suggests that the higher-order statistics ignored in QL equations are important here for recovering statistical stability of the symmetric (S) turbulent state.

5.4. Nonlinear effects

Extended OS stability analysis applied to the body-force F2 state shows a significant stabilisation of the eigenvalue corresponding to the second streamwise wavenumber, but shows no significant difference for the other two unstable eigenvalues of the standard OS approach. Repeating the numerical experiments described above applied to the F2 state shows that both the QL equations and the LNS equations yield rapid exponential growth of the perturbations applied to the time-dependent base flow; see figure 10(b). Periodically renormalising the growing perturbation in the latter LNS case and continuing the simulation shows sustained growth over a period of $O(1000)$ time units (not shown). This behaviour is in contrast to what was found for state S where the perturbation saturated (see figure 10a). On the face of it, this experiment seems to show that state F2 is linearly unstable but this does not necessarily follow. Time stepping a perturbation using the LNS equations (or in the tangent space of a turbulent attractor) can continually sample the stretching or contracting dynamics along the attractor which can dominate the anticipated decaying behaviour perpendicular to the attractor. It seems stretching along the attractor dominates for state F2 but not for state S. A comparison of the leading growing structures in figure 12 is consistent with this: the QL and LNS flows bear no similarity with the predictions of the mEOS equations. These LNS calculations for the state F2 highlight perfectly the difficulty in assessing the (linear) statistical stability using the Navier–Stokes equations touched on in the introduction.

Figure 12. Comparison of the leading eigenvectors from the mEOS $m=1$ (a) and mEOS $m=3$ (b) analyses with the instantaneous snapshots of the perturbation fields obtained from the QL experiment at $t=100$ (c) and linear NS experiment at $t=100$ (d) for the F2 state. The structures of the leading mEOS eigenvectors seem unrelated to the growing perturbations in the QL or LNS regimes.

Repeating the experiment (with $E_{pert} = 10^{-5} E_{base}$) but now time stepping the perturbation field with the full Navier–Stokes equations, we see the perturbation energy saturate after an initial transient; see figure 10(b). Importantly, and as anticipated, the statistics of the turbulent state F2 recover albeit slowly. Figure 13 shows the relative difference between the correlation matrices averaged over a rolling time window between the perturbed and unperturbed F2 state. The statistical differences ebb away after an initial transient of growth. Since F2 is statistically linearly stable, this initial growth indicates that the statistical evolutionary operator linearised around the state F2 fixed point has to be non-normal. Finally, figure 14 shows the evolution of the streamwise perturbation velocity field. We see that the highly localised structure visible at $t=10$ spreads in the channel at $t=20$ and reaches its saturated state at $t=30$ that does not change qualitatively even at much longer times $t=200$. It takes much longer for the statistics to recover suggesting the adjustment along the attractor happens much quicker than normal to it; see figure 13.

Figure 13. Time evolution of the relative difference between correlation matrices of the full flow $C^{m0}$ and base flow $C^{m0}_0$ shown for $m=1$ (blue), $m=2$ (red) and $m=3$ (green). The correlation matrices were time averaged over a rolling time window to display the statistical behaviour. Shown for the full nonlinear NS experiment where the difference slowly approaches $0$ signifying the recovery of the same statistical state.

Figure 14. Snapshots of the streamwise velocity component of the perturbation field taken from the full nonlinear Navier–Stokes experiment with the F2 state. At time $t = 10$ the unstable structure is still growing until time $t = 20$ when it reaches its maximum size and starts to spread out through the channel, at time $t = 30$ the perturbation saturates and shows no qualitative change at a much later time $t = 200$. Results are shown for (a) $t = 10$, (b) $t = 20$, (c) $t = 30$, (d) $t = 200$.

6. Discussion

6.1. Summary

We first summarise what has been done in the paper. The motivating objective has been to develop a theoretical approach which can assess the statistical stability of a turbulent state and goes beyond just examining the mean flow as Malkus did in the 1950s (Malkus 1956). That such an approach should exist in some form is predicated on the fact that the statistics of a disturbed turbulent flow always return to their undisturbed values if the disturbance is not too large regardless of whether it decays or not. That is, a turbulent flow is assumed linearly statistically stable (note this does not preclude the possibility of the statistical disturbance initially increasing before ultimately decaying away as the linear statistical operator may be non-normal). The way forward is to include more statistical information than just the mean of the flow in the analysis and the simplest such is the second rank cumulant or two-point equal-time correlations of the flow fluctuations. To do this, we have proposed to examine the statistical stability of the turbulent flow from within a cumulant equation framework truncated to just consider these statistical quantities (called CE2). Even in CE2, however, the ensuing statistical spectral problem is so unwieldy that a substitute has had to be sought. Here the well-known connection between CE2 and the QL equations can be exploited and we have argued that a suitably designed QL spectral problem defined in § 3 can act as a good proxy to detect statistical instability. The key in designing this QL problem is to translate a statistically steady base flow in CE2 into a representative steady base flow in the QL problem with no pretence that this base flow satisfies the QL equations itself (the QL equations offer an approximation of uncertain accuracy to the Navier–Stokes equations after all).

The issue of how the statistical spectral problem in CE2 is related to the QL spectral problem is an intriguing one that seems not to have been discussed before. It is important because it underpins the tacit assumption that a CE2 simulation starting from given initial conditions should exactly shadow (statistically) the equivalent QL simulation given one is derivable from the other (e.g. Marston et al. 2008; Tobias et al. 2011; Tobias & Marston 2013). However, this ignores the difference in dimensionality between the two formulations and the possibility of triggering unstable ‘non-physical’ perturbations in the larger CE2 problem (meaning perturbations that cannot exist in the QL setting). In Appendix A we have tried to clarify the different types of instability that exist in these formulations – type A and type B – and how they are related. In particular, for type A instabilities, the spectral QL problem captures the most unstable disturbance possible in the CE2 system and so is an accurate proxy for CE2: CE2 can only be type A unstable if the QL problem is. For type B, the situation is less certain. While it remains the case that any type B instability in the QL setting is mirrored in CE2 and CE2 will likely have further type B instabilities, it is not clear whether CE2 can have a type B instability and the QL problem not. Given the as-yet unreported consequence of this – the CE2 equations could suffer a bifurcation whereas the QL equations not – it is tempting to assume that this is not possible and then that the QL analysis is also a good proxy for type B perturbations. With this one caveat, probing the QL spectral problem – named here EOS analysis – then appears the best way to examine the statistical stability of a flow where the two-point correlations as well as the mean profile of the flow are incorporated.

Despite improving upon the CE2 spectral problem, EOS analysis is still unfortunately quite costly to implement so a reduced version was proposed in § 3.1 – mEOS analysis – which has a similar cost to the usual OS analysis of the flow mean and is the most straightforward extension of doing linear stability analysis only around the mean flow. Just as in OS analysis, one chooses a particular Fourier mode $(m,n)$ to analyse but now extra information about the corresponding base Fourier mode is included. Minimally EOS analysis was found to capture the qualitative effect of EOS analysis with either both producing an improved estimate of statistical stability (e.g. the $m=2$ mode for state F1; see figure 5c) or offering no correction of the OS eigenvalue (e.g. the $m=1$ mode for state F1). Reassuringly, both approaches never made a poorer statistical stability prediction than OS analysis at least over the turbulent channel flow states tested. Minimally EOS analysis works comparatively well because the structure of the key eigenmodes in the EOS problem tend to consist of one dominant (primary) fluctuation field interacting with a mean flow perturbation (see figure 6). This is presumably because the comparatively small size of the cumulants or base fluctuation fields relative to the mean profile (e.g. indicated by the singular values in figure 9) means that only the few fluctuation fields close to marginality on the mean flow play a role. This may well change with increasing $Re$ but the relative size of the cumulants should also increase possibly acting to deprioritise this near-marginal stability property. Clearly higher $Re$ flows warrant further study.

6.2. Implications

By incorporating two-point fluctuation correlations, EOS analysis and its slimmed down version mEOS analysis, represent an attempt to generate a better way to assess the statistical stability of a turbulent flow than just considering the mean profile. Malkus (1956) wanted to use this to constrain the set of realisable turbulent mean profiles and it still remains an interesting requirement given the paucity of predictive alternatives. The concept of marginal stability still periodically rears its head in the literature to post-rationalise turbulent mean profiles (see, e.g. Brauckmann & Eckhardt (2017) and references therein) and EOS analysis represents a finer tool for this. As noted in the introduction, the mean flow is increasingly taken as the starting point for resolvent flow analysis. To move this towards a full wall-bounded turbulence theory, a way has to be found to close the cycle between the fluctuation field and the mean. Just maybe, requiring statistical stability with respect to the dominant fluctuation fields could help constrain this.

A byproduct of this work is that it has also suggested an improved approach for predicting coherent structures in turbulent flows. Coherence here means weakly damped or growing fluctuation and mean disturbances that emerge from the linear operator built around the base flow rather than, say, proper orthogonal decomposition modes derived from the two-point correlations that have no direct relation to the underlying equations. Initial work (Crighton & Gaster 1976; Gaster et al. 1985; Roshko 1993) linearising around the turbulent mean proved successful in free-shear flows but less so when viscosity plays a leading role. Admittedly, EOS or mEOS analysis needs extra (second cumulant) information about the base flow that adds cost to the analysis but the potential for better predictions is clear. A popular (cheap) alternative approach has been to apply OS analysis around the mean profile that incorporates an eddy viscosity to acknowledge the presence of the base fluctuations (Reynolds & Hussain 1972; Sen & Veeravalli 2000; del Alamo & Jimenez 2006; Sen et al. 2007; Cossu et al. 2009; Hwang & Cossu 2010a,b). This eddy viscosity is taken to be that which is needed to sustain the turbulent mean profile and so makes no attempt to incorporate fluctuation information for each excited wavenumber of the base flow as in EOS analysis. In the well-known triple expansion of Hussain & Reynolds (Hussain & Reynolds 1970), the coherent wave part ($\tilde {f}$ in their equation (1.1)) is the equivalent of $\boldsymbol {\delta } \boldsymbol {\tilde {u}}^{mn}$ that there is excited by a wavemaker and then studied as it travels downstream in a turbulent flow. Orr–Sommerfeld analysis was found unable to reconcile what was seen experimentally but EOS analysis should do better.

The cumulant framework examined here to motivate EOS analysis was used only to the first non-trivial level where only the first and second cumulants are retained. Whether this represents a good approximation or not to modelling the turbulent flow itself is a question not considered here. Rather the focus is on the associated spectral problem characterising the statistical stability within this framework taking the turbulent base flow as a given from DNS or experiment. Clearly considering higher cumulants beyond the second would be better but the computational cost of handling the ensuing spectral problem becomes prohibitive. Interestingly, ‘projecting down’ to an underlying Navier–Stokes-type problem fails immediately unless no closure at all is performed, which means CE$\infty$ is considered. This gives what we have called the iEOS analysis that is as costly as EOS analysis and unfortunately does not have a ‘minimal’ version. Again, only targeting certain subsets of fluctuation fields – e.g. key triads suggested by resolvent flow analysis – may yet make this tractable.

6.3. Future perspectives

In terms of future work, an obvious focus is to improve the main approximation in EOS analysis that is deriving a representative base flow. Here, we have used the leading singular vector of the second cumulant for each excited (vector) wavenumber but other strategies involving information from higher-order cumulants are certainly possible (leaving aside the cost of computing these$\ldots$). This would seem sensible in the case of iEOS analysis that is a proxy for CE$\infty$ where all the cumulants are retained and the only other cumulant system beyond CE2 that allows a Navier–Stokes-type proxy. Pragmatically, a steady base flow has been sought that leads to a (conceptually) simple eigenvalue problem. Trying to go beyond this profoundly complicates matters.

One interesting future direction hinted at in the introduction is to relax the averaging procedure to just being over the streamwise direction or, more generally, redefine what constitutes the ‘mean’. A key feature we have been able to exploit is the connection between CE2 and the QL equations that still holds if the mean is more flexibly defined to include some wavenumber subset of the fluctuation field. The then ‘generalized’ QL equations (Marston et al. 2016) have the advantage of including more of the nonlinearity in the Navier–Stokes equations and can interpolate between the QL equations and the full Navier–Stokes equations depending on how the mean is defined. A minimal mean extension would be to include one spanwise wavenumber field whereas including all fluctuation fields with no streamwise variation is equivalent to a mean defined by only streamwise averaging. There are many theoretical reasons for considering the latter, where the mean is $U(y,z)\boldsymbol {\hat {x}}$ rather than just $U(y)\boldsymbol {\hat {x}}$, in wall-bounded flows (e.g. the SSP/VWI mechanism Waleffe 1997; Hall & Sherwin 2010), but capturing it experimentally is a very time-consuming task and using it reduces the predictive power of resolvent flow analysis. Even in computations, the spanwise homogeneity needs to be broken otherwise there is no reason to suppose the mean flow should not be spanwise invariant (and then the spanwise structure of the mean profile needs to be divorced from the source of spanwise inhomogeneity). Nevertheless, this choice is starting to receive attention at least at the level of doing OS analysis on $U(y,z)\boldsymbol {\hat {x}}$ (see Lozano-Duran et al. (2021) and references herein) and this analysis can be extended as envisaged here.

Many challenges remain. Statistical stability is an important fundamental feature of turbulent flows yet is difficult to assess and, hence, to exploit. Hopefully, the work reported here represents a useful but admittedly small step forward.