On the shape of resolvent modes in wall-bounded turbulence

TL;DR

This paper introduces a new analytical approach to approximate leading resolvent modes in wall-bounded turbulence, enabling rapid predictions of flow structures based on simplified scalar operators and wavepacket pseudoeigenmodes.

Contribution

It develops a methodology using prescribed analytic functions to approximate resolvent modes, especially wall-normal vorticity modes, in shear-driven turbulent flows, reducing computational complexity.

Findings

Accurate approximation of wall-normal vorticity modes using scalar differential operators.

Mode shape variation with wavenumber and Reynolds number captured by low-dimensional equations.

Framework relates to classical lift-up and Orr mechanisms in shear flows.

Abstract

The resolvent formulation of the Navier$\text{--}$Stokes equations gives a means for the characterization and prediction of features of turbulent flows$\text{---}$such as statistics, structures and their nonlinear interactions$\text{---}$using the singular value decomposition of the resolvent operator based on the appropriate turbulent mean, following the framework developed by McKeon & Sharma (2010). This work will describe a methodology for approximating leading resolvent (i.e., pseudospectral) modes for shear-driven turbulent flows using prescribed analytic functions. We will demonstrate that these functions, which arise from the consideration of wavepacket pseudoeigenmodes of simplified linear operators (Trefethen 2005), in particular give an accurate approximation of the class of nominally wall-detached modes that are centered about the critical layer. Focusing in particular on modeling wall-normal vorticity modes, we present a series of simplifications to the governing equations that result in scalar differential operators that are amenable to such analysis. We demonstrate that the leading wall-normal vorticity response mode for the full Navier$\text{--}$Stokes equations may be accurately approximated by considering a second order scalar operator, equipped with a non-standard inner product. The variation in mode shape as a function of wavenumber and Reynolds number may be captured by evolving a low dimensional differential equation in parameter space. This characterization provides a theoretical framework for understanding the origin of observed structures, and allows for rapid estimation of dominant resolvent mode characteristics without the need for operator discretization or large numerical computations. We relate our findings to classical lift-up and Orr amplification mechanisms in shear-driven flows.

Figures

Equations

$$\dot{\mathbf{u}}=\mathbf{g}(\mathbf{u}).$$

$${\dot{\mathbf{u}}}^{\prime }=\mathbf{g}({\mathbf{u}}_0+{\mathbf{u}}^{\prime })={\frac{\partial \mathbf{g}}{\partial \mathbf{u}}|}_{{\mathbf{u}}_0}{\mathbf{u}}^{\prime }+\mathbf{f}(\mathbf{u}),$$

$$\left(-i\omega -{\frac{\partial \mathbf{g}}{\partial \mathbf{u}}|}_{{\mathbf{u}}_0}\right){\hat{\mathbf{u}}}^{\prime }=\widehat{\mathbf{f}(\mathbf{u})},$$

$${\hat{\mathbf{u}}}^{\prime }={\left(-i\omega -{\frac{\partial \mathbf{g}}{\partial \mathbf{u}}|}_{{\mathbf{u}}_0}\right)}^{-1}\widehat{\mathbf{f}({\mathbf{u}}^{\prime })}={\mathcal{H}}_{\omega }\widehat{\mathbf{f}({\mathbf{u}}^{\prime })},$$

$${\mathcal{H}}_{\omega }=\sum _{j=1}^{\infty }{\sigma }_j{\psi }_j{\varphi }_j^{\ast },$$

$$⟨{\mathbf{x}}_1,{\mathbf{x}}_2⟩={\mathbf{x}}_1^{\bar{T}}\text{mathsfbi}M{\mathbf{x}}_2,$$

$${\sigma }_1$$

$${\psi }_1$$

$${\varphi }_1$$

$${\sigma }_1$$

$${\psi }_1$$

$${\varphi }_1$$

$${\Lambda }_{ϵ}(\mathcal{L})=\left\{z:(\mathcal{L}+\mathcal{E})\mathbf{u}=z\mathbf{u},\text{ for some }\mathbf{u}\text{ and }\mathcal{E},\text{ with }\Vert \mathcal{E}\Vert \le ϵ\right\}.$$

$$z\in {\Lambda }_{ϵ}(\mathcal{L})\backslash \Lambda (\mathcal{L})⟺\Vert (-z+\mathcal{L}{)}^{-1}\Vert \ge {ϵ}^{-1},$$

$${\mathcal{D}}_h^j=(ih{)}^j\frac{d^j}{(dy{)}^j},$$

$$({\mathcal{L}}_{h}u)(y)=\sum _{j=0}^n{c}_j(y){\mathcal{D}}_h^{j}u(y).$$

$$({\mathcal{L}}_h{v}_h)(y)=\sum _{j=0}^n{c}_j(y)k^j{v}_h(y)=f(y,k)v_h(y),$$

$$\text{Im}\left(\frac{{\partial }_{y}f}{{\partial }_{k}f}\right)>0.$$

$$\Vert (-\lambda +{\mathcal{L}}_h)\psi (y;y_{\ast },k_{\ast })\Vert \le M^{-1/h}$$

$$\frac{|\psi (y;y_{\ast },k_{\ast })|}{\max |\psi (y;y_{\ast },k_{\ast })|}\le C\exp \left(-b(y-y_{\ast }{)}^2/h\right),$$

$$\left(\begin{array}{c}\hat{v}\\ \hat{\eta }\end{array}\right)=\underset{{\mathcal{H}}_{\mathbf{k}}}{\underbrace{{\left(\begin{array}{cc}-i\omega +{\Delta }^{-1}{\mathcal{L}}_{os}& 0\\ i{k}_z{U}_y& -i\omega +{\mathcal{L}}_{sq}\end{array}\right)}^{-1}}}\left(\begin{array}{c}{\hat{f}}_v\\ {\hat{f}}_{\eta }\end{array}\right).$$

$${\mathcal{L}}_{os}$$

$${\mathcal{L}}_{sq}$$

$$⟨({\hat{v}}_1,{\hat{\eta }}_1),({\hat{v}}_2,{\hat{\eta }}_2)⟩=\frac{1}{k^2}{\int }_y\left(-{\bar{\hat{v}}}_1\Delta {\hat{v}}_2+{\bar{\hat{\eta }}}_1{\hat{\eta }}_2\right)dy,$$

$$⟨{\hat{z}}_1,{\hat{z}}_2⟩$$

$${⟨{\hat{z}}_1,{\hat{z}}_2⟩}_{\Delta }$$

$${\mathcal{H}}_{vv}$$

$${\mathcal{H}}_{\eta \eta }$$

$$\left(\begin{array}{c}\hat{v}\\ {\hat{\eta }}_{os}\end{array}\right)$$

$$\left(\begin{array}{c}0\\ {\hat{\eta }}_{sq}\end{array}\right)$$

$$\hat{v}$$

$${\hat{\eta }}_{os}$$

$${\hat{\eta }}_{sq}$$

$${\mathcal{H}}_{vv}={\Delta }^{-1}{\mathcal{H}}_{\eta \eta }\Delta .$$

$${\mathcal{H}}_{\eta v}^{Lin}$$

Full text

1 Introduction

The identification of pertinent structures that arise in the transition towards, and as coherent features within, turbulent wall-bounded flows has been the focus of much research over the past several decades. Qualitatively, such analysis includes identification and classification of empirically-observed structures, such as near-wall streaks (Kline et al., 1967), hairpin structures (Theodorsen, 1952; Head & Bandyopadhyay, 1981) and their grouping in large-scale motions (Zhou et al., 1999; Guala et al., 2006), and very large-scale motions or superstructures (Kim & Adrian, 1999; Guala et al., 2006; Hutchins & Marusic, 2007). For more details concerning the (sometimes debated) properties, taxonomy and dynamics of such structures, see reviews such as Robinson (1991), Smits et al. (2011), and Jiménez (2018), and references within.

On a quantitative level, a starting point for the prediction of coherent structure comes from consideration of properties of the governing equations, most typically in linearised form. Features emergent in wall-bounded turbulent flows often bear little resemblance to modes identified from classical stability analysis (Drazin & Reid, 2004), where two-dimensional modes are predicted to be the least stable by Squire’s theorem, and mean-linearised flows often have only stable eigenvalues (Reynolds & Tiederman, 1967; Del Alamo & Jimenez, 2006; Cossu et al., 2009). Perhaps the most important realisation in the study of such linear operators is the fact that their non-normality can result in high amplification (in either the time or frequency domain), which cannot be predicted from their spectra alone (Böberg & Brösa, 1988; Butler & Farrell, 1992; Reddy & Henningson, 1993; Trefethen et al., 1993; Schmid & Henningson, 1994; McKeon & Sharma, 2010; Bamieh & Dahleh, 2001; Jovanović & Bamieh, 2005; Schmid, 2007; Hwang & Cossu, 2010; Schmid & Henningson, 2012). Operator nonnormality is responsible for both finite-time energy growth of the linear system from a from a given initial condition, and the amplification resulting from continual forcing (be it stochastic or harmonic). Indeed, these two notions are mathematically related via the Kriess constant of the operator (Schmid, 2007).

While most initial works considered laminar base flows (with one notable exception being Farrell & Ioannou (1993), who considered mean-linearised equations using stochastic forcing), more recent developments have used similar methods for prediction of turbulent features by linearizing about states computed from turbulent data and/or models (which is most often a mean state). For example, Del Alamo & Jimenez (2006) and Cossu et al. (2009) have demonstrated transient growth of near wall-streaks and large scale motions for turbulent channel and zero-pressure-gradient boundary layers respectively, while Schoppa & Hussain (2002) has showed that a turbulent mean with the addition of low-speed streaks can give rise to the growth of structures in the near-wall region of turbulent channel flow.

The characterisation of the nonlinear Navier–Stokes equations in the frequency domain as the linear resolvent operator acting on the nonlinear “forcing” terms, as developed by McKeon & Sharma (2010) has been particularly fruitful for elucidating operator-based predictions of structure in wall-bounded turbulence, including very large scale motions and their scaling (McKeon & Sharma, 2010), and hairpin structures (Sharma & McKeon, 2013). McKeon (2017) summarizes further developments in this area. Note that this approach has also been applied in other contexts, such as in the study of cavity (Gómez et al., 2016; Qadri & Schmid, 2017), airfoil (Yeh & Taira, 2018), and jet (Garnaud et al., 2013; Jeun et al., 2016; Towne et al., 2018; Schmidt et al., 2018) flows.

While such developments are recent and ongoing, many of the underlying physical mechanisms in shear flows have been understood for at least several decades. These include the Orr mechanism (Orr, 1907; Jiménez, 2013), which amplifies upstream-leaning disturbances while tilting them towards the downstream direction, and the lift-up mechanism (Landahl, 1980, 1975), in which wall-normal disturbances lead to large streamwise responses. Streamwise inhomogeneity provides an additional mechanism for amplification of disturbances (Chomaz, 2005; Hack & Moin, 2017), which, as with the Orr and lift-up mechanisms, arises through different aspects of non normality of the underlying linear operator (Symon et al., 2018).

Recent results concerning the pseudospectral properties of certain matrices (Trefethen & Chapman, 2004) and linear differential operators by (Trefethen, 2005) have revealed criteria by which pseudomodes exist that are localised in both space and spatial frequency, and are highly amplified by the associated resolvent operator. This method of analysis has been used, for example, in the analysis of swept wing flow by Obrist & Schmid (2010, 2011).

An alternative route to quantitatively define coherent structure in turbulent flows is to consider data collected from simulations or experiments. The proper orthogonal decomposition (Lumley, 1967), which computes spatial modes of highest energy in a dataset, is perhaps the most ubiquitous such method. The original formulation computes a set of energetically optimal modes for each temporal frequency. This method was recently shown by Towne et al. (2018) to give modes equivalent to resolvent response modes, under the assumption that the true forcing to the system, arising from the nonlinear term in the Navier–Stokes equations, results in uncorrelated resolvent response mode expansion coefficients. Note that quantitatively correct prediction of characteristics of wall-bounded turbulence using resolvent-based approaches can be improved by considering the properties of the nonlinear forcing terms, which are dependent on second order turbulence statistics (Zare et al., 2017), or through the addition of an eddy viscosity term to the linearised dynamics (Illingworth et al., 2018).

Connections may also be made between resolvent analysis and the dynamic mode decomposition (Schmid & Sesterhenn, 2008; Schmid, 2010; Rowley et al., 2009), which computes spatial modes from the eigendecomposition of an linear operator that best matches the evolution of the data between adjacent snapshots, as described in Sharma et al. (2016) and Towne et al. (2018). A widely used POD variant computes modes without filtering into temporal frequencies, which gives a basis that can be used to project the Navier–Stokes equations onto to obtain a reduced-order model (Berkooz et al., 1993; Holmes et al., 2012). Note in particular that this has been used to describe the dynamics of coherent structures in wall-bounded turbulent flows, such as those associated with the near-wall cycle (Aubry et al., 1988). Reviews of modal decomposition techniques, and their use in reduced-order modelling, may be found respectively in Taira et al. (2017) and Rowley & Dawson (2017).

With the exception of very simple systems, the identification of mode shapes typically requires resorting to numerical methods, applied either to the discretised (and most-often linearised) governing equations, or to data collected from experiments or numerical simulations of their evolution. In this work, we describe methods to approximate mode shapes and amplification mechanisms without requiring the formation of discretised operators, for shear flows with arbitrary mean velocity profiles. In section 2, we provide a review of the mathematical concepts that underpin our analysis. Section 3 presents the formulation and sample results for resolvent analysis of a turbulent boundary layer, and presents a sequence of simplifications to the governing equations that retain the correct features of the leading resolvent response mode. In Section 4, we detail a procedure for estimating mode shapes by solving an optimisation problem using a prescribed template function, which may be derived either from approximations to exact solutions of the given operator, or from the wavepacket pseudomode theory introduced in section 2.

2 Mathematical preliminaries: The resolvent and pseudospectra of a linearised operator

This section presents material on the pseudospectral analysis of a mean-linearised system, which will provide background for the analysis in later sections. Section 2.1 introduces the resolvent formulation of a nonlinear system. The singular value decomposition of the resolvent and its connection with pseudospectra is discussed in Section 2.2. Section 2.3 discusses the underlying theory behind the existence of wavepacket pseudospectral modes, which will be related to our subsequent analysis.

2.1 The resolvent form of a nonlinear dynamical system

We begin by considering a nonlinear dynamical system

[TABLE]

Let ${\bm{u}}_{0}$ denote the temporal mean of the state of the system, where we are assuming that the dynamics are statistically stationary. Expressing the system state as ${\bm{u}}(t)={\bm{u}}_{0}+{\bm{u}}^{\prime}(t)$, we may rewrite equation 1 as

[TABLE]

where we have linearised about the mean state, but retained the full dynamics of the system with the remaining nonlinear dynamics ${\bm{f}}({\bm{u}})$. Taking a Fourier transform in time, equation 2 may be expressed as

[TABLE]

where $\hat{\cdot}$ denotes a Fourier-transformed function. The mean-subtracted state of the system may then be expressed by

[TABLE]

where we refer to $\mathcal{H}_{\omega}$ as the associated resolvent operator for this system for a given frequency $\omega$, where we are assuming here that this inverse exists (i.e., that $i\omega$ is not an eigenvalue of $\left.\frac{\partial{\bf g}}{\partial{\bm{u}}}\right|_{{\bm{u}}_{0}}$). It is important to note that we have not made any approximations to the nonlinear system at this point, and have only made the assumption that the system is statistically stationary in time with a well-defined mean.

2.2 The singular value decomposition of the resolvent operator

From equation 3, the properties of the mean-subtracted state ${\bm{u}}^{\prime}$ will depend both on the nature of the nonlinear term $\hat{\bm{f}}$, and the properties of the linear operator $\mathcal{H}_{\omega}=(-i\omega+\mathcal{L})^{-1}$, where following on from section 2.1, we let $\mathcal{L}=-\left.\frac{\partial{\bf g}}{\partial{\bm{u}}}\right|_{{\bm{u}}_{0}}$. In particular, if $\mathcal{H}_{\omega}$ amplifies a small number of directions, or “modes” to a much larger degree than all others, then so long as these directions are excited to some extent by $\hat{\bm{f}}$, this can allow prediction of the dominant features of ${\bm{u}}^{\prime}$ by studying only $\mathcal{H}_{\omega}$.

More precisely, we consider the singular value decomposition of the resolvent operator

[TABLE]

with $\sigma_{k}>\sigma_{k+1}$ for all $k$. For the remainder of this work, we will consider $\mathcal{H}_{\omega}$ to be a discretised operator, with a corresponding singular value decomposition (SVD) defined in the same manner as 4, but with a finite sum of modes. The SVD of a linear operator requires the definition of an inner product (or more precisely, inner products on both the input and output spaces), which prescribes how the adjoint is prescribed and computed (and also induces the norms used in defining various properties of the SVD). For discrete operators, we may characterise an inner product by a positive definite weighting matrix ${\mathsfbi{M}}$ for which

[TABLE]

where $\cdot^{\bar{T}}$ denotes the conjugate transpose. If ${\mathsfbi{M}}$ characterises the inner product on the spaces containing both the forcing and response functions of a finite-dimensional linear operator ${\mathsfbi{L}}$ (which we consider the discretisation of a linear operator $\mathcal{L}$), the adjoint ${\mathsfbi{L}}^{*}$ satisfies $\left\langle{\mathsfbi{L}}{\bm{x}}_{1},{\bm{x}}_{2}\right\rangle=\left\langle{\bm{x}}_{1},{\mathsfbi{L}}^{*}{\bm{x}}_{2}\right\rangle$, from which one may show that ${\mathsfbi{L}}^{*}={\mathsfbi{M}}^{-1}{\mathsfbi{L}}^{\bar{T}}{\mathsfbi{M}}$. Throughout this paper we will largely work with (infinite-dimensional) operators acting on a continuous domain, with the understanding that numerical computation of the SVD is performed on a finite-dimensional discrete approximation.

For some aspects of this work, it will be convenient to think of the leading singular values and vectors of $\mathcal{H}_{\omega}$ as defined in the following manner:

[TABLE]

where here and throughout we use the spectral norm when taking the norm of an operator. Note that we also have the relationship $\psi_{1}=\sigma_{1}^{-1}\mathcal{H}_{\omega}\phi_{1},$ which may also be rearranged as $\phi_{1}=\sigma_{1}\mathcal{H}^{-1}_{\omega}\psi_{1}$. In other words, $\phi_{1}$ gives the shape of the forcing which gives rise to the largest amplification (a factor of $\sigma_{1}$) upon the application of the operator $\mathcal{H}_{\omega}$, the result of which is the (scaled by $\sigma_{1}$) response mode $\psi_{1}$. It will also be useful to recognise that these singular values and vectors are related to the smallest singular values and vectors of $\mathcal{H}_{\omega}^{-1}=\left(-i\omega+\mathcal{L}\right)$, by

[TABLE]

In this work, we will make use of these definitions, that allow us to define leading singular values and vectors as solutions to an optimisation problem. Roughly speaking, we will search for analytic functions that become as small as possible upon the action of $(-i\omega+\mathcal{L})=\mathcal{H}^{-1}_{\omega}$, and reason from equation 10 they will be close approximations to leading resolvent response modes.

These ideas may be formalised as follows. For a linear operator $\mathcal{L}$ and some $\epsilon>0$, define the $\epsilon$-pseudospectrum as a set $\Lambda_{\epsilon}(\mathcal{L})\subset\mathbb{C}$ satisfying

[TABLE]

Here $\mathcal{E}$ is an operator mapping between the same spaces as $\mathcal{L}$. Note in particular that we have the equivalent definition of pseudospectra based on the norm of the resolvent:

[TABLE]

where here the $\backslash$ operator refers to set exclusion. In particular, for any $z\in\mathbb{C}$, $\epsilon_{z}=\min\{\epsilon:z\in\Lambda_{\epsilon}(\mathcal{L})\}=\sigma_{1}^{-1}(z),$ where $\sigma_{1}(z)$ is the largest singular value of $(-z+\mathcal{L})^{-1}$. Furthermore, the corresponding pseudoeigenvector ${\bm{u}}$ satisfying $(\mathcal{L}+\mathcal{E}){\bm{u}}=z{\bm{u}}$, for some $\mathcal{E}$ with $\|\mathcal{E}\|=\epsilon_{z}$ is the left singular vector of $(-z+\mathcal{L})^{-1}$ corresponding to $\sigma_{1}(z)$. Note that the presence of the operator norm in the definition means that, unlike when looking at the spectrum, the definition of pseudospectra requires that we work in a vector space equipped with a norm, which the $\epsilon$-pseudospectrum associated with an operator is dependent upon.

2.3 Conditions for the existence of wavepacket resolvent modes

Here we briefly describe the conditions under which a linear differential operator permits wavepacket pseudomodes corresponding to small $\epsilon$ (or equivalently, large resolvent norm). More complete description of this phenomenon, as well as the related proofs, may be found in Trefethen & Embree (2005) and Trefethen (2005). For some (small) parameter $h>0$, we may define a family of scaled differential operators

[TABLE]

which we assume to act on a closed finite interval. We may then assemble any arbitrary $n$-th order differential operator $\mathcal{L}_{h}$ acting on a scalar variable $u(y)$ by

[TABLE]

Defining a test function $v_{h}(y)=\exp(-iky/h)$ for a complex scalar $k$, we have

[TABLE]

where we refer to $f(y,k)$ as the symbol corresponding to the family of differential operators $\mathcal{L}_{h}$. We will consider potential pseudoeigenvalues $\lambda=f(y_{*},k_{*})$. The symbol $f(y,k)$ is said to satisfy the twist condition for real $k$ if we have

[TABLE]

When the twist condition is satisfied for some $(y_{*},k_{*})$, one it can be a pseudomode $\psi(y;y_{*},k_{*})$ (of unit norm) that has phase variation matching $v_{h}$ close to $y_{*}$, satisfying both

[TABLE]

for some $M>1$, and also that

[TABLE]

for some $b,C>0$. In other words, we have modes that are spatially localised near $y_{*}$, have localised spatial frequency $k_{*}$, and are “asymptotically good” pseudoeigenfuctions. Note that the region of $\mathbb{C}$ where the twist condition is satisfied is independent of the choice of norm. Note also that this theory developed in Trefethen (2005) builds upon earlier observations of certain classes of equations by Davies (1999a, b). More generally, these ideas are closely related to classical WKBJ expansions, which have also been utilised recently by Leonard (2016) in the context of studying approximate inviscid solutions for channel flow.

In this work, using these ideas as justification and inspiration, for a given temporal frequency (where $\lambda=i\omega$), we will assume that the leading resolvent response mode may be closely approximated by a function that is localised in the wall-normal direction, both spatially (as for functions satisfying equation 18) and in the frequency domain (as is the case for the test functions $v_{h}$). When this assumption is justified, it allows for a reformulation of resolvent analysis in terms of finding the spatial width and frequency of a template function $\psi$ that minimises $\|(-\lambda+\mathcal{L})\psi\|$, and thus are close to the true leading resolvent response mode (from equation 10).

3 The behaviour of leading resolvent modes in wall-bounded turbulence

This section will apply the concepts introduced in section 2 to explore methods by which the shape of resolvent modes may be approximated. In this section and elsewhere, we will largely focus on a boundary layer configuration, which is assumed to be approximately homogeneous in the streamwise direction (as well as being homogeneous in the spanwise direction), though the analysis holds for any (approximately) parallel shear-driven turbulent flow.

In section 3.1 we introduce a resolvent formulation of the Navier–Stokes equations. The setup is similar to that developed in McKeon & Sharma (2010), though following more closely the formulation used in Rosenberg & McKeon (2018). Following this, in section 3.2 we will present a sequence of simplifications to the governing equations, that will render them amenable to application of wavepacket pseudomode theory. We will additionally present some sample results to motivate these developments, as well as those later presented in section 4.

3.1 A resolvent formulation of the Navier–Stokes equations

We will restrict attention to flows which are (approximately) homogenous in the streamwise $(x)$ and spanwise $(z)$ directions, with a mean streamwise velocity $(u,v,w)=(U,0,0)$ that varies in the wall normal $(y)$ direction. We will consider the incompressible Navier–Stokes equations in (wall-normal) velocity-vorticity form, and will take Fourier transforms in the streamwise and spanwise directions, with wavenumbers given by $k_{x}$ and $k_{z}$. Applying the procedure detailed in section 2.1 gives

[TABLE]

where $\hat{v}$ and $\hat{\eta}=ik_{z}\hat{u}-ik_{x}\hat{w}$ are the Fourier-transformed wall-normal velocity and vorticity fields, $U_{y}$ is the wall-normal gradient of the streamwise velocity, and the Laplacian $\Delta=\partial_{yy}-k_{\perp}^{2}$, with $k_{\perp}^{2}=k_{x}^{2}+k_{z}^{2}$. The Orr-Sommerfeld (OS) and Squire (SQ) operators are given respectively by

[TABLE]

Further details concerning the equivalence of this formulation to one using primitive variables (in particular the implicit restriction of the forcing to the solenoidal component) is given in Rosenberg & McKeon (2018). The resolvent operator $\mathcal{H}_{\bm{k}}$ is now parametrised by the spatiotemporal wavevector $\bm{k}=(\omega,k_{x},k_{z})$.

The SVD of the resolvent operator is computed using an inner product that is proportional to kinetic energy for a given set of spatial wavenumbers (Gustavsson, 1986; Butler & Farrell, 1992), defined by

[TABLE]

where the overbar denotes complex conjugation.

In later sections, we will consider the following two scalar inner products arising from the components of 22:

[TABLE]

That is to say, when considering scalar operators, we will assume that the inner product is the “standard” one unless using a $\Delta$-subscript. We will also refer to $\left\langle\cdot,\cdot\right\rangle_{\Delta}$ as the Laplacian inner product.

We will now present some sample results that will motivate many of the developments in the remainder of this paper. We consider a turbulent boundary layer with friction Reynolds number $\Rey_{\tau}=900$. Mean velocity profiles are obtained from the DNS data of Wu et al. (2017), by averaging data at a single streamwise location. Figure 1 shows leading resolvent response mode shapes for this system. The wavenumbers $k_{x}=\pi/2$, and $k_{z}=2\pi$ and wavespeed ($c=0.8U_{\infty}$) have been chosen to be consistent with the typical size of large-scale motions in zero-pressure-gradient boundary layers (e.g., Monty et al. (2009); Kovasznay et al. (1970); Cantwell (1981); Saxton-Fox & McKeon (2017)). Here and throughout, the resolvent operator is discretised using a Chebyshev collocation method, utilising the toolbox of Weideman & Reddy (2000).

We observe that the mode is dominated by the wall-normal vorticity component, which is centred on and localised around the critical layer, and has an approximately linear variation in phase within this region. These observations, which are typical of modes that are “detached” from the wall, suggests that this numerically computed mode resembles a wavepacket mode as described in section 2.3. Here and throughout, we consider a mode to be detached from the wall if its shape is not substantially effected by the wall’s presence. In section 4 we will seek to predict the shape of this mode without explicitly computing an SVD (or indeed, without explicitly forming a discretised resolvent operator). Making this procedure tractable, however, will require simplification of the governing equations, which are described in section 3.2.

3.2 Simplifications to the Navier–Stokes operator for resolvent mode approximation

Following Rosenberg & McKeon (2018), rather than studying the full system we may start by separating the resolvent modes into those being forced by $f_{v}$ and $f_{\eta}$ separately. In particular, defining scalar resolvent operators

[TABLE]

the Orr-Sommerfeld (OS) and Squire (SQ) modes may be computed by taking SVD’s of the reduced resolvent operators given as

[TABLE]

While the modes computed from different subsystems are no longer orthogonal, Rosenberg & McKeon (2018) demonstrated they can provide a basis that better captures dynamical features of the system. Our interest in this decomposition is primarily concerned with using this decomposition to simplify the required analysis. Figure 2 shows the $\eta$-component of the response modes of the OS and SQ subsystems for the same parameters considered in figure 1. We observe in particular that the OS operator gives the same mode shape as the full system, with the SQ mode also sharing the same qualitative characteristics. This suggests that the ability to predict and understand the shape of the vorticity response of the full system can be reduced to studying only the OS subsystem.

In order to make the ensuing analysis more tractable, we will require additional simplifications. Note first that the components of the OS and SQ sub-operators may be written in scalar form as

[TABLE]

where $\mathcal{H}_{\eta v}=(ik_{z}\mathcal{H}_{\eta\eta}U_{y}\mathcal{H}_{vv})$ is the off-diagonal term in $\mathcal{H}_{os,\bm{k}}$ in equation 27.

When considering the SVD of these subsystems, it is important that we use the appropriate inner product on both the input and output spaces. Since the response is energetically dominated by wall-normal vorticity, it is reasonable that a close approximation to the $\eta$ response of the full system (which itself is approximated by the OS response) can be obtained just by considering the scalar operator in equation 30. Figure 3 shows that this assumption is indeed valid for the sample parameters considered in this section. As well as simplifying the analysis, reducing the full resolvent operator to a scalar operator brings us closer to being able to apply the wavepacket pseudomode theory discussed in section 2.3. An additional requirement for the direct application of this theory is that the operator under consideration be representable as a differential operator in the form of equation 14. The operator in equation 30 does not satisfy this property, as the OS component contains the inverse Laplacian.

Motivated by this, we now present an approximating simplification to 14, which we believe to be novel. This development will take advantage of the similarity between the OS and SQ operators, and will consider variations of the “standard” inner products discussed in section 3.1. We start by taking advantage of the localised nature of the resolvent response modes, and assume that the mean velocity profile may be linearised about the critical layer location, $y_{c}$. Under this assumption, we have

[TABLE]

This means that the operator that we are attempting to simplify is

[TABLE]

where $U_{y_{c}}$ is the wall-normal gradient in the mean velocity profile at the critical layer location, $y_{c}$, and the $Lin$ superscript denotes that the mean velocity profile has been linearised about $y_{c}$. Figures 2 and 3 show that the resolvent response modes of $\mathcal{H}_{\eta\eta}$ do not quantitatively match those of $\mathcal{H}_{\eta v}$. This is because the leading forcing mode of $\mathcal{H}_{\eta\eta}$ does not coincide with the mode most amplified by $\mathcal{H}_{vv}$. Given that the leading resolvent modes are dependent on the choice of inner product, it is reasonable to ask if a different choice of inner product could be used in the SVD of $\mathcal{H}_{\eta\eta}$, such that the leading forcing mode approximately coincides with the leading response of $\mathcal{H}_{vv}$. Physically, what we are seeking to do is to find an inner product weighting that will most amplify the lift-up mechanism (at the expense of optimising the Orr mechanism in isolation).

Let $\mathcal{H}^{*}_{\eta\eta}$ denote the adjoint of $\mathcal{H}_{\eta\eta}$ with respect to the standard scalar inner product, equation 23. If we instead consider the “Laplacian” inner product, we obtain

[TABLE]

This has the same form as equation 32, except with the (standard) adjoint of $\mathcal{H}_{\eta\eta}$, which amounts to taking the complex conjugate of the critical layer term. Note also that $\mathcal{H}_{\eta\eta}$ with the regular inner product has leading forcing and response modes (i.e., right and left singular vectors) that differ only in the direction of the phase change. We also then have

[TABLE]

which is identical to equation 32, meaning in particular that both share the same resolvent forcing and response modes. If we are to weight this operator explicitly to compute an SVD in the Laplacian scalar norm, we obtain

[TABLE]

Therefore, computing this SVD should give a leading forcing mode for $\mathcal{H}_{\eta\eta}$ that coincides with the leading response mode of $\mathcal{H}_{vv}$, thus minimising the “projection loss” of the total amplification. In effect, we are modifying the inner product used for $\mathcal{H}_{\eta\eta}$ such that it’s leading forcing mode aligns with the leading response mode of $\mathcal{H}_{vv}$. This analysis therefore relies on the assumption that the total amplification is dominated by the lift-up, rather than the Orr, mechanism. With this assumption, this trick allows us to compute optimal response mode shapes by only considering the Squire operator, which is a differential operator of the general form given in equation 13, rather than the full OS operator, which is not. It will be shown later that this assumption holds except for large $k_{\perp}^{2}$, in which case the Laplacian approaches a constant, and so $\mathcal{H}_{\eta\eta}$ and $\mathcal{H}_{vv}$ converge to the same operator. Figure 3 shows that this modification to the inner product allows us to closely match the wall-normal vorticity component of the response mode of the full Navier–Stokes system, for the sample parameters chosen. We will discuss reasons for the success of this approximation, and in particular the relative balance between optimising the Orr and lift-up mechanisms, in section 4.

4 Predicting the shape of resolvent modes

This section will present a method that allows for the prediction of resolvent mode shapes, focusing in particular on the wall-normal vorticity component. The main idea will be to assume the existence of a mode that is localised in both wall-normal location, and wall-normal spatial frequency, and then find the parameters which result from maximum amplification of the resolvent operator, or equivalently, minimisation of the action of its inverse, on a given function. We will start by predicting the shape of a model operator (which is closely related to the Squire operator) in sections 4.1, which relates mode shapes to asymptotic expansions of Airy functions, and 4.2, which approximates mode shapes under the assumption of the existence of wavepacket pseudomodes. This latter approach is then extended to consider the shape of modes for the full incompressible Navier–Stokes system in section 4.3, using the approximations presented in section 3.2. We will validate our method first on laminar Couette flow with a linear velocity profile in section 4.3, before returning to the turbulent boundary layer configuration in section 4.4.

4.1 Relationship between wavepacket resolvent modes and Airy functions

Before applying wavepacket pseudomode theory more generally, in this section we focus on a model (Airy) operator that is equivalent to the Squire operator with a linearised mean velocity profile. This analysis will show how numerically computed resolvent response modes relate to exact analytical solutions to simplified governing equations. Linearising the mean velocity profile about the critical layer, the Squire equation reduces to

[TABLE]

For clarity, we let $R=k_{x}U_{y_{c}}Re$, and consider the closely-related operator

[TABLE]

where $\lambda=y_{c}+\omega_{c}i$, and $\omega_{c}=R^{-1}{k_{\perp}^{2}}$, and we are assuming that $k_{x}\neq 0$. Note that we are generally interested in performing resolvent analysis along the axis of neutral stability (i.e., for real-valued frequencies with no growth or decay), but here we are incorporating the viscous dissipation term $\omega_{c}$ into $\lambda$, so in general $\lambda$ is complex for this analysis. $\mathcal{T}$ is a complex Airy-type operator, and is identical to the operator considered, for example, in Reddy et al. (1993) as a model for studying the pseudospectra of the Orr-Sommerfeld operator. The equation $\mathcal{T}u=0$ has a general solution that can be expressed by two independent Airy functions, such as

[TABLE]

with

[TABLE]

Note that distinct solutions are obtained through the choice of contour in the complex plane when applying standard Laplace transform methods.

Shown in figure 4 are eigenvalues and selected eigenfuctions of $\mathcal{T}$ with Dirichlet boundary conditions on a finite domain $y\in[-1,1]$, with $R=3000$. In figure 5, these numerical eigenfunctions are compared to $\Ai[z]$, $\Ai[e^{2\pi i/3}z]$, and $\Ai[e^{-2\pi i/3}z]$, and it is seen that $\Ai[e^{2\pi i/3}z]$ (when appropriately-scaled) closely matches the numerically computed eigenfunctions on the finite domain. Note that if we had chosen eigenvalues on the left branch (i.e., with negative real component), then the $\Ai[z]$ solution would have been accurate for all cases. This can be seen more clearly by overlaying the $y$-domain on the amplitude of the Airy functions in the complex plane, as shown in figure 6.

On an infinite domain (and when acting on square-integrable functions) $\mathcal{T}$ has an empty spectrum, and $\epsilon$-pseudospectra that are independent of location on the real axis (that is, the boundaries of the pseudospectrum for various $\epsilon$ are horizontal lines).

As discussed in Reddy et al. (1993), a solution to $\mathcal{T}u=0$ which is within $\epsilon$ of a function that satisfies the appropriate boundary conditions will be an $\epsilon$-pseudoeigenfunction. We consider the pseudospectra of $\mathcal{T}$, again for the finite domain $y\in[-1,1]$. As discussed in section 2.2, this amounts to computing the leading singular values and vectors of the resolvent of $\mathcal{T}$ for each $\lambda\in\mathbb{C}$ of interest. In particular, we may seek analytic approximations of these resolvent (optimal pseudospectral) modes using the same method as for eigenfunctions. In other words, for a given $\lambda\in\mathbb{C}$, we are interested in finding a function which is as close as possible to an eigenmode of $\mathcal{T}$. Noting that the Airy functions are continuous, the fact that there are regions of $\mathbb{C}$ for which they are very close to satisfying the boundary conditions of the finite domain suggest that these regions correspond to very large values of resolvent norm $\sigma_{1}$.

Figure 7 shows $\epsilon$-pseudospectra for $\mathcal{T}$, along with optimal pseudoeigenmodes (i.e., resolvent modes) and Airy function approximations thereof, for specified locations in the centre of the domain. In the central region, the pseudospectral properties of $\mathcal{T}$ resemble that of the operator on an infinite domain, where resolvent modes, which are “detached” from the boundary, are invariant to horizontal translation within this region. Moreover, resolvent modes within this region will resemble translated versions of eigenmodes along the inclined branches at the same vertical position. We observe that, as $\epsilon$ increases and we approach the real axis, the analytical Airy functions become less accurate, owing to their growth in magnitude occurring closer and closer to the critical layer, eventually saturating the component resembling the numerically computed response mode. As discussed in Reddy et al. (1993), these Airy functions can be used to compute constructive lower bounds to the pseudospectrum of $\mathcal{T}$, though these bounds are inaccurate whenever the Airy functions diverge far away from $y=\Real(\lambda)$. The fact that $\omega_{c}=R^{-1}k_{\perp}^{2}$ means that the most physically-relevant region of the complex plane is in the upper half-plane, which means that the Airy functions themselves are not close approximations of physically-relevant resolvent modes.

Away from the origin, one may approximate Airy functions by simpler expressions that do not involve integrals. The derivation of such asymptotic approximations involves finding a location in the complex that has a dominant contribution to the contour integral. We may derive the approximation (e.g., Olver (2014); Vallée & Soares (2010))

[TABLE]

which is accurate for large $|z|$, where here $z=\gamma(y-\lambda)$, where (with reference to equation 38) we have $\gamma=(iR)^{1/3}$, $\lambda=ik_{\perp}^{2}/R=\omega_{c}i$, and we assume for now that $y_{c}=0$. Expanding the exponential term gives

[TABLE]

and thus

[TABLE]

where we must take care when dealing with multi-valued roots. In particular, for a bounded approximation, we require that

[TABLE]

If $\lambda=|\lambda|e^{i\theta_{\lambda}}$, then we find that this condition is only satisfied when $\theta_{\lambda}=\frac{3\pi}{2}.$ Note that this corresponds to the lower half plane region where Airy functions can give close approximations to the numerically-computed solutions. Assuming $y_{c}=0$, so that $\lambda=\omega_{c}i$, for $\omega_{c}<0$ we obtain an approximation of the form

[TABLE]

4.2 Predicting wavepacket modes for a model operator

Section 4.1 demonstrated that, in certain regimes, resolvent modes may be approximated by a function of the general form

[TABLE]

where $a\in\mathbb{R}$ and $b>0$ are functions of $\omega_{c}$ and $R$, and

[TABLE]

is a constant that gives $\psi$ unit norm with respect to the regular scalar inner product.

This section will present an alternative analysis with this template function as a starting point. Applying wavepacket pseudomode theory as described in section 2.3, one may show that $\mathcal{T}$ satisfies the twist condition (equation 16) is satisfied within the half strip

[TABLE]

The $\text{Imag}(\lambda)<0$ condition is consistent with the predicted condition for a spatially localised mode from the truncated asymptotic expansion of an Airy function in equation 41. Note that the region described in equation 43 is outside the region that is most physically relevant for resolvent analysis of the Squire operator, where a positive $k_{\perp}^{2}$ and $R$ give a positive imaginary component of $\omega_{c}$. Despite this, figure 7 shows that even outside this region, the numerically computed pseudomode still maintains the same spatially localised structure. To predict the shape of resolvent modes outside this region, we will seek to directly optimize the shape parameters in the template function, equation 42. In particular, we wish to find the values of $a$ and $b$ for which $\|\mathcal{T}\psi\|$ is minimized. That is, we seek the minimum of the cost function

[TABLE]

where we are assuming that our mode is “detached” and localised, so can integrate over an infinite domain, and are assuming a standard inner product for this model operator (note that including a constant multiplicative factor of $k_{\perp}^{-2}$ as in equation 23 does not affect mode shapes or singular values). For simplicity, we again take $y_{c}=0$. Note that the leading resolvent singular value is related to this cost function by

[TABLE]

Noting that we have

[TABLE]

it can be shown that

[TABLE]

From this, we may compute

[TABLE]

Contours of this cost function for $R=3000$ and several values of $\omega_{c}$ are shown in figure 8.

We may analytically determine the locations of these minima as follows. This cost function is minimised when

[TABLE]

The equations 46 and 47 may be solved to give the parameters $\{a,b\}$ minimising $J_{\mathcal{T}}$, with care taken to select the correct minimising solution. It is also possible to obtain a parametrised family of solutions as follows. Assuming that $R$ is fixed, suppose that we seek the optimal values $\left(a(\omega_{c}),b(\omega_{c})\right)$ for a range of values of $\omega_{c}$. By implicitly differentiating equations 46 and 47 with respect to $\omega_{c}$, the following ordinary differential equations governing the evolution of the optimal mode shape parameters may be found:

[TABLE]

Note that this approach assumes that the global minimiser of $J_{\mathcal{T}}$ stays on the same branch, and is smoothly continuous with $\omega_{c}$. Figure 9 shows that the predicted values of the optimal values of parameters $a$ and $b$ closely match both the asymptotic approximations (for $\omega_{c}<0$), and the values obtained from fitting the numerically computed modes. The width parameter for this fit to the numerically-computed mode data is found by fitting a Gaussian function to the amplitude of the computed mode using MATLAB’s it } command. The phase parameter is ound by considering the gradient of the phase in a small region near the critical layer location. The predicted shape parameters are obtained from evolving the equations 48 and 49 from an initial condition obtained by solving equations 46 and 47 directly (at $\omega_{c}=0$). Here, and in subsequent sections, such solutions (along with the symbolic computation of integrals) are obtained using Mathematica. Note that one could also use the asymptotic approximation 41 for initial conditions, where we must start from a sufficiently large negative value of $\omega_{c}$ to ensure that the initial conditions are accurate. It is additionally shown that the value of the cost function $J_{\mathcal{T}}$ for these optimal shape parameters closely matches the value of the cost function that may be computed directly from the numerically computed singular value. In essence, this shows that the optimal mode and amplification across all functions is closely approximated by the optimal over the class of functions of the form given in equation 42.

We may also consider the behaviour of the optimal shape parameters as $R$ varies. Keeping $\omega_{c}$ constant, from equations 46-47 we obtain

[TABLE]

In the limit of large $R$ and for $\omega_{c}\neq 0$, (and assuming restrictions on the growth of $a$ and $b$ with $R$), we have the approximations

[TABLE]

from which we may infer the scalings $a\to C$, $b\propto R^{1/2}$. For large $R$ but with $R\omega_{c}\ll 1$, we instead obtain

[TABLE]

which can be shown to permit a consistent asymptotic solution with $a\propto R^{1/3}$ and $b\propto R^{2/3}$. This scaling of $b$ corresponds to a mode with that scales with $R^{1/3}$, which agrees with the critical layer scaling (Drazin & Reid, 2004), and can also be inferred from the transformation described in equation 38.

Figure 10 shows that these scalings closely match the shape parameter trends obtained from the numerically computed modes. Moreover, these shape parameters are accurately predicted from finding optimal parameters of cost function $J_{\mathcal{T}}$, which can be obtained either from solving equations 46 and 47 for each $\omega_{c}$ and $R$, or by evolving the differential equations 50 and 51 over $R$. It should be noted that, even for this model system, computing the leading resolvent modes directly can require substantial computational resources for large $R$, highlighting an advantage of approximations using the methods described in this section. We lastly note that one could seek alternative parameterisations of the two-dimensional space of optimal mode shapes, such as by varying $R$ while keeping $\omega_{c}R$ ($=k_{\perp}^{2}$) constant.

4.3 Predicting the shape of resolvent modes for laminar Couette flow

In this section, we will extend the analysis in section 4.2 to consider the Navier–Stokes operator and suboperators considered in sections 3.1 and 3.2. We again will restrict our attention to predicting the vorticity component of the leading response mode of the resolvent operators $\mathcal{H}_{\bm{k}}$, $\mathcal{H}_{os,\bm{k}}$ and $\mathcal{H}_{sq,\bm{k}}$. For the Squire suboperator, the analysis is almost identical to that of the model operator considered in section 4.2. To predict the mode shapes of the full system, we rely on the fact that the response is dominated by the effect of the OS suboperator, which we analyse using the simplifying approximations introduced in section 3.2. In particular, this also results in a methodology similar to that used in section 4.2, but with a modification to use the Laplacian scalar inner product (equation 24.

We begin by again assuming mode shapes of the form

[TABLE]

but now to satisfy the unit-norm requirements of the relevant inner products (equation 23) for the SQ operator, equation 24 for the OS) we have

[TABLE]

where, with reference to the previous section, we have $k_{\perp}^{2}=\omega_{c}R$. Using the approximation to the Orr-Sommerfeld operator introduced in section 3.2, and noting that $U_{y_{c}}=1$, the relevant cost functions for optimizing the shape of the SQ and OS template functions are

[TABLE]

For laminar Couette flow we have $U_{y_{c}}=1$ for all $y_{c}$ in the domain, but we keep this term so that the equations are directly applicable for arbitrary mean velocity profiles. Also note that while the cost function for the Squire modes is “exact” given the assumed template function and a linear mean velocity profile, the OS modes are relying on the accuracy of a simplified operator (being the scalar Squire operator with a Laplacian inner product) capturing the correct behavior of the response mode.

Substituting in the mode template functions (and dropping the subscripts on $a$ and $b$ parameters for brevity), explicit expressions for the cost functions are

[TABLE]

[TABLE]

As before, the minima of these cost functions may be found by selecting the appropriate solution to the equations

[TABLE]

[TABLE]

[TABLE]

Figure 11 compares the predicted shape parameters to those obtained from fitting these parameters to numerically-computed modes, which are computed for $Re=1000$ and various values of $k_{x}$ and $k_{z}/k_{x}$. We observe in subplots (a) and (c) that the Squire mode parameters are accurately predicted from solving equations 60 and 61. Subplots (b) and (d) show that solving equations 62 and 63 predicts the fitted parameters to the leading resolvent response mode of the Squire operator with Laplacian inner product. Furthermore, this prediction also accurately predicts the shape of the wall-normal vorticity component of leading response modes for the full Navier–Stokes system. Note that the scaling laws for small $k_{x}$ are identical to those for small $R$ observed in figure 10. The trends at high $k_{x}$ differ from those in figure 10, due to the fact that here we keep $k_{z}$ proportional to $k_{x}$ (and thus $k_{\perp}\propto R$), whereas in figure 10 constant $\omega_{c}$ resulted in $k_{\perp}^{2}\propto R$. These trends can again be inferred from studying the dominant terms in the evolution equations for the governing parameters, which can be computed from considering the partial derivatives of $J_{sq}$ and $J_{os}$ with respect to $a$ and $b$, which we omit for brevity.

For large $k_{x}$, the mode shapes for both operators converge, which may be explained by the fact that the Laplacian operator is dominated by the constant $k_{\perp}^{2}$ term in this regime. As a consequence, here analysis of the Squire operator gives accurate predictions of the behaviour of the full Navier–Stokes system, particularly for the mode width, which approaches a constant value with increasing $k_{x}$. There is some difference between the phase gradient parameter $a$ for the full system and for the Squire system with both the standard and Laplacian inner product (with the value for the full system lying between those for the two simplified systems), though this difference decreases as $k_{x}$ increases and the phase variation decays towards zero.

Figure 12 plots predicted and true mode amplitudes for several of the cases considered in figure 11. The predicted mode shapes for $\eta$ closely match the computed modes for both the Squire and full Navier–Stokes system, with the modes for the latter being very close to those of the OS subsystem. The largest discrepancy arises in the $k_{x}=1$ case, where the tails of the mode amplitudes are significantly heavier than those of a Gaussian distribution, and extend far from the critical layer towards the wall. This phenomenon is related to a wider distribution in the $v$-component of the mode (not plotted), and gives a larger variation in fitted $b$ values for this $k_{x}$ for the full Navier–Stokes system, as observed in figure 11.

Figure 13 compares the leading resolvent singular values for the Squire and Navier–Stokes systems to those estimated from the minima of the cost functions $J_{sq}$ and $J_{os}$, for the same wavenumbers considered in figure 11. The cost functions are able to accurately predict the singular values of the scalar Squire operator with both the regular and Laplacian inner product, though this is only accurate for the full Navier–Stokes system for large $k_{x}$ (and thus large $k_{\perp}^{2}$). In particular, these model operators are incapable of predicting the increase in singular value for high-aspect-ratio (i.e., large $k_{z}/k_{x}$) modes that is observed, though we have shown that they still accurately predict mode shapes in this regime.

4.4 Predicting mode shapes for a turbulent boundary layer

This section applies the methodology developed in sections 4.2 and 4.3 to a turbulent boundary layer. By linearizing the mean velocity profile about the critical layer, the equations for predicting mode shape parameters are the same as those developed in section 4.3 (i.e., equations 58-63), though now $U_{y_{c}}$ and $R=k_{x}U_{y_{c}}\Rey$ are dependent on the critical layer location. Figure 14 shows the predicted and true response mode shapes for the Squire and Navier–Stokes resolvent operators for two pairs of spatial wavenumbers and various critical layer locations. Note that the case where $\{k_{x},k_{z},c\}=\{\pi/2,2\pi,0.8U_{\infty}\}$ as considered in figures 1-3 as representative of a typical large-scale motion, is included. The location typical of very large-scale motions ($c\approx 0.6U_{\infty}$) is also included, though these structures would correspond to slightly smaller streamwise wavenumbers than $k_{x}=\pi/2$. We observe that mode amplitude and phase variation (in the local region of high amplitude) is accurately estimated, provided the mode is not substantially affected by the presence of the wall. This shows in particular the validity of using a mean velocity profile linearised about the critical layer to estimate mode shapes.

Figure 15 compares predicted and fitted mode shape parameters as a function of $k_{x}$ for various aspect ratios $k_{z}/k_{x}$, for a wavespeed $c=0.8U_{\infty}$, with predicted and true modes shapes plotted in figure 16. We observe the same trends, and similar accuracy in prediction of parameters as for laminar Couette flow (figure 11).

5 Discussion and conclusions

This work has presented a method for approximating leading resolvent response modes for quasi-parallel shear-driven flows. This method relies on the assumption that the true mode may be closely approximated by a simple template function, the general form of which can be reasoned from consideration of wavepacket pseudomode theory. In essence, the method reduces the space of possible mode shapes from an infinite-dimensional space (which in practice is approximated by a high dimensional space defined by the numerical discretisation) to a two-dimensional family of functions. Once this template function is identified, the optimal shape parameters (which govern the width and phase variation of the mode) may be found as the minimisers of a cost function, which is directly related to the resolvent norm of the underlying operator. In practice, this amounts to finding the roots of a pair of coupled equations, which may be arranged to be polynomials in the shape parameters. In addition, it is possible to derive differential equations in parameter space that govern the evolution of these optimal shape parameters. Importantly, this method precludes the need for the formulation and decomposition of discretised linear operators, leading to substantial reduction in computational cost. The extent of the reduction in computational cost is dependent on the size of the discretisation, and on the extent and resolution of the parameter space (e.g., wavenumbers and temporal frequencies) that one wishes to study. The method may be readily applied to a model operator, as considered in section 4.2, and in the analysis of the Squire operator in sections 4.3-4.4. Application for the full Navier–Stokes system relies on additional simplifications to arrive at a scalar differential operator which has a leading response mode (left singular vector) which approximates the wall-normal vorticity component of the response mode of the Navier–Stokes resolvent operator, as detailed in section 3.2. This simplification is made in three steps. Firstly, the observation that the wall-normal vorticity response is dominated by the Orr-Sommerfeld component (equation 27). Secondly, that this operator is in turn dominated by the vorticity response to velocity forcing, governed by the scalar operator (equation 30). Finally, the leading left singular vector of this scalar operator may be approximated by that of the Squire operator furnished with the scalar Laplacian inner product (equation 24). This results in an operator which may be studied by applying the same techniques as those used for the model/Squire system. Note that, even without this analysis, study of the Squire operator (for which the associated cost function has a simpler form) typically gives the same qualitative behaviour observed for the full system. In particular, the Squire system obeys many of the same scaling laws of mode shape parameters with wavenumbers and Reynolds numbers. More detailed analysis of scaling and self-similarity properties of resolvent modes in wall-bounded flows are given in Moarref et al. (2013).

While the vorticity response mode shapes of the Squire and Orr-Sommerfeld resolvent operators are qualitatively similar, they represent quite different physical phenomena. The Squire resolvent operator describes amplification through forcing in the same component, with an upstream-leaning optimal forcing mode giving a downstream-leaning response mode. This is a manifestation of the classical Orr mechanism. Note in particular that the Squire forcing mode has the same amplitude profile as the response mode, but with an opposite phase variation. The Orr-Sommerfeld sub-operator on the other hand typically has a leading response mode that is dominated by the wall-normal vorticity component, but is forced primarily by wall-normal velocity, in a manner resembling the lift-up mechanism. Despite this phenomenological difference, we have shown that the shape of the response may be accurately predicted from the Squire operator with a modified inner product.

This work has focused on characterising the shape of the (dominant) wall-normal vorticity component resolvent modes for parallel wall-bounded flows. It is possible that these methods could be extended for application to more complex geometries, and other mode components. For most of the cases considered in this work, the cost function had only one local minimum corresponding to a wavepacket mode (i.e., with $b>0$), with the exception being for negative $\omega_{c}$ for the model operator considered in section 4.2, as seen in figure 8. More complex geometries might result in more complex cost functions, for which more care must be taken to select the true global optimum. Note that the operators considered here also typically have a large spectral gap between the first and second singular value. Future work could also seek to identify modes corresponding to additional singular values, which would be of particular interest for situations where the operator does not exhibit low-rank behaviour (i.e., there does not exist a large gap between the leading and second singular values). Similar methods could also be applied to study optimal forcing modes in more detail, and to compute nonlinear forcing terms induced from analytic approximations to wavepacket modes. This would provide an alternative route to study phenomena such as the self-similarity of the nonlinear forcing (Sharma et al., 2017).

In terms of the methodology itself, there are a number of possible refinements that could be investigated. For example, nonlinear terms in the expansion of the mean velocity profile about the critical layer could be retained, and additional terms in the template function (for example, a term could be added to allow for the phase variation to be cubic in $y-y_{c}$). The former modification might be particularly useful when dealing with mean velocity profiles that have stationary points, such as in channel flow. Such additions would lead to more complex cost functions, but the same techniques of analysis should be applicable. The accuracy of the assumed form of the modes relies on the mode being far enough away from the wall. Further extensions of the methodology could seek to explicitly model the effect of the wall.

The authors acknowledge support from the Air Force Office of Scientific Research grant FA9550-16-1-0232 (program manager Ivett Leyva). The authors also thank Peter Schmid, Anthony Leonard, and Kevin Rosenberg for valuable discussions, and Xiaohua Wu for allowing us to use his turbulent boundary layer database.