Variational formulation of resolvent analysis

PHYSICAL REVIEW FLUIDS

, 013905 (2022)

Benedikt Barthel

Salvador Gomez

, and Beverley J. McKeon

Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, California 91125, USA

(Received 16 September 2021; accepted 13 January 2022; published 28 January 2022)

The conceptual picture underlying resolvent analysis is that the nonlinear term in the

Navier-Stokes equations acts as an intrinsic forcing to the linear dynamics, a description

inspired by control theory. The inverse of the linear operator, defined as the resolvent, is

interpreted as a transfer function between the forcing and the velocity response. From

a theoretical point of view this is an attractive approach since it allows for the vast

mathematical machinery of control theory to be brought to bear on the problem. However,

from a practical point of view, this is not always advantageous. The inversion of the linear

operator inherent in the control theoretic definition obscures the physical interpretation

of the governing equations and is prohibitive to analytical manipulation, and for large

systems it leads to significant computational cost and memory requirements. In this

work we suggest an alternative, inverse-free, definition of the resolvent basis based on

an extension of the Courant–Fischer–Weyl min-max principle in which resolvent modes

are defined as stationary points of a constrained variational problem. This definition leads

to a straightforward approach to approximate the resolvent (response) modes of complex

flows as expansions in any arbitrary basis. The proposed method avoids matrix inversions

and requires only the spectral decomposition of a matrix of significantly reduced size

as compared to the original system. To illustrate this method and the advantages of the

variational formulation we present three examples. First, we consider streamwise constant

fluctuations in turbulent channel flow where an asymptotic analysis allows us to derive

closed form expressions for the optimal resolvent mode. Second, to illustrate the cost-

saving potential and investigate the limits of the proposed method, we apply our method to

both a two-dimensional, three-component equilibrium solution in Couette flow and, finally,

to a streamwise developing turbulent boundary layer. For these larger systems we achieve

a model reduction of up to two orders of magnitude. Such savings have the potential to

open up RA to the investigation of larger domains and more complex flow configurations.

DOI:

10.1103/PhysRevFluids.7.013905

I. INTRODUCTION

Resolvent analysis (RA) can be used to give insight into the forced response of a linearized

dynamical system. This concept was introduced by Trefethen

et al.

[

] and Jovanovi

c and Bamieh

[

] who considered the stability and amplification of linearly stable flows to external forcing. These

ideas were later applied to turbulent flows by McKeon and Sharma [

] who interpreted the nonlinear

term in the Navier-Stokes equations (NSE) as a forcing to the linearized system. The conceptual

framework of RA is inspired by control theory (CT), such that the resolvent operator, the inverse

of the linearized operator, is interpreted as a transfer function from the forcing to the response.

A singular value decomposition (SVD) of the discretized resolvent operator provides two distinct

orthonormal bases (left and right singular modes) for both the response and the forcing, ordered by a

bbarthel@caltech.edu

2469-990X/2022/7(1)/013905(35)

013905-1

BARTHEL, GOMEZ, AND MCKEON

set of gains (singular values) which quantify the linear amplification of the system. This CT-inspired

framework has proven theoretically useful since it is conceptually straightforward and benefits from

years of established mathematical machinery. However, from a practical point of view, the reliance

on the inversion of the linear operator poses some difficulties. It obscures the analytical tractability

of the equations and is computationally costly for all but the simplest systems.

Early studies using RA were largely focused on wall-bounded shear flows with only a single

nonhomogeneous spatial dimension for which the computational cost of the inversion and SVD of

the operator is trivial [

–

]. In these cases, linearly amplified length scales identified by RA were

found to correlate with the energetically active scales observed in experiments and simulations,

and the corresponding resolvent modes capture the qualitative features of the coherent structures

observed in wall turbulence [

]. In particular, resolvent modes have been found to exhibit the

self-similar behavior characteristic of the attached eddy hypothesis proposed by Townsend [

Towne

et al.

[

] have elaborated the assumptions under which resolvent response modes correlate

with spectral proper orthogonal decomposition (SPOD) modes computed from data, illustrating that

RA can predict coherent structure in the full flow field. More recently RA has also been extended to

2D flows such as boundary layers [

], the flow behind bluff bodies [

], exact coherent states

(ECS) [

], and turbulent jets [

]. In particular, modal analysis techniques including RA have

been used by a variety of authors to implement flow control strategies, for example, to suppress

vortex shedding [

], and delay flow separation [

]. For these 2D flows the computational cost

and memory requirements becomes considerable and thus the further extension to 3D flows has

generally remained limited.

The community has endeavoured to address these computational challenges through innovation

in novel methods of estimating resolvent modes. One area of research has been in so called

“matrix-free” methods such as the work of Martini

et al.

[

] who use the transient and steady

state responses of the periodically forced linearized system and its corresponding adjoint system

to estimate the action of the resolvent operator. Another avenue of investigation inspired by the

field of data analysis has been in “equation-free” methods such as Herrmann

et al.

[

] who use

dynamic mode decomposition (DMD) modes to estimate the linear dynamics of a system from a

time series of data. Others have made use of iterative Arnoldi Algorithms that replace the cost of

calculating the SVD and a matrix inverse with the cost of an LU decomposition and a few matrix

multiplications [

]. Furthermore, algorithms such as randomized SVD and others have made it

possible to efficiently and accurately compute singular modes of data sets that would otherwise be

prohibitively expensive [

–

The previously cited research has focused on the CT interpretation of RA and the SVD-based

definition of resolvent modes. In this work we take an alternative approach and propose an equiva-

lent definition based on an extension of the Courant–Fischer–Weyl min-max principle (CFL). The

CFL principle itself and the concept of an optimal forcing and maximum gain are well understood

and have been used extensively by a wide range of authors. For example, Dawson and McKeon [

]

derived analytical models of the optimal resolvent mode in wall-bounded shear flows, Monokrousos

et al.

[

] and Garnaud

et al.

[

] investigated the optimal forcing structure in a Blasius boundary

layer and incompressible jets, respectively, and Towne [

] computed data-driven resolvent modes

in the context of turbulent jets. Here we present an explicit extension from the CFL principle, to

what we coin “variational resolvent analysis” (VRA), which constitutes an alternative definition of

the resolvent basis that includes all modes.

This variational definition is based on the solutions of the Euler-Lagrange equations associated

with the constrained variation of the operator norm of the linearized dynamics. Critically, this

definition does not involve the inversion of any operator, which is useful from both a theoretical

and practical sense. The inversion of large matrices is both costly and obscures the intuitive

interpretation of the underlying linear differential operator. The extension of the CFL principle to

include all resolvent modes was used by Sipp and Marquet [

] who defined the resolvent forcing

modes as the solutions to a generalized Eigenvalue problem, the current work builds on their results

by additionally avoiding the inversion of the linear operator and the need for any adjoint equations.

013905-2

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

While in general the resulting Euler-Lagrange equations remain difficult to solve exactly, this

variational formulation allows for the approximation of two- or three-dimensional resolvent modes

as an expansion in any convenient basis, such as for example a much cheaper one-dimensional

resolvent basis, an analytical basis such as that described by Dawson and McKeon [

], or a

data-driven one. Further, it requires only the eigenvalue decomposition of a matrix of reduced size.

In this paper we illustrate how this variational definition is useful in both gaining physical insights

by allowing for analytical progress in simplified systems, and by reducing computational cost in

complex systems. To illustrate the former we consider the case of streamwise constant fluctuations

in wall-bounded shear flows. To investigate the latter we first perform RA around a 2D

3C exact

coherent solution, where we find that we can accurately approximate the resolvent response modes

and reduce the computational complexity by an order of magnitude. Finally, the VRA formulation

is applied to a streamwise developing turbulent boundary layer, where the near wall modes can be

predicted with a 97% reduction in computational cost using resolvent modes calculated using a 1D

mean flow.

The paper is organized as follows. In Sec.

we derive the proposed variational definition. In

Sec.

III

we use the variational formulation to analyze streamwise constant structures in turbulent

channel flow. In Secs.

and

we consider RA applied to both streamwise periodic and streamwise

developing two-dimensional, three-velocity-component (2D

3C) systems to illustrate the computa-

tional cost and memory saving potential of the proposed VRA formulation. In Sec.

we analyze

the uncertainty and potential sources of error in our method. We provide discussion of the results

and the outlook for future applications in Sec.

VII

and conclude in Sec.

VIII

II. MATHEMATICAL FORMULATION

Let us consider a general forced linear system

∂

−

(1)

where

represents a spatial-linear differential operator and

(

)

(

)

∈

∞

. The state vari-

ables

and

are referred to as the “response” and “forcing,” respectively. We consider the temporal

Fourier transfer of Eq. (

) and define the spatiotemporal linear operator

(

)

≡

−

(2)

as well as the resolvent operator

(

)

≡

(

)

−

(3)

which is classically interpreted as a transfer function from the forcing to the response,

(4)

For readability we have dropped explicit reference to the dependence on

. An SVD of the resolvent

∞

(5)

results in a pair of distinct sets of basis functions for the response (

) and forcing (

) and are

referred to as the resolvent “response modes” and “forcing modes,” respectively. These are ordered

by their gains

that are ordered in descending order, representing the

th largest linear gain

possible. Here and throughout this work superscript

denotes the Hermitian adjoint, or for discrete

matrices the conjugate transpose.

013905-3

BARTHEL, GOMEZ, AND MCKEON

A. A variational definition of resolvent modes

A key contribution of this work is the observation that resolvent response modes may be

equivalently defined as the stationary points,

∗

, of the operator norm of

under the condition

that the argument

∗

satisfies some norm constraint. More explicitly, the resolvent modes of the

linear operator

are defined as the stationary points of the functional

(6)

subject to the constraint

(7)

We note that in general the norms

≡

and

≡

need not be the same, such as

for example in the Orr-Sommerfeld and Squire decomposition discussed in Sec.

III

. Following the

notation of Herrmann

et al.

[

] the Cholesky factorization may be used to decompose the weight

matrix

(8)

This allows a general norm to be related to the Euclidean 2 norm. In other words, we can express

any arbitrary user-defined norm as

(9)

where

is simply a label used to distinguish between different norms.

The method of Lagrange multipliers allows us to combine Eqs. (

), (

), and the definition Eq. (

)

to formulate a constrained variational problem and define a Lagrangian

(

)

−

(10)

Here

and

may be either interpreted as continuous differential operators or discrete matrices.

The vanishing of the variation with respect to the conjugate state

∗

is a necessary and sufficient

condition for the stationarity of Eq. (

). The reader is referred to Appendix

for a derivation of

this property based on the work of Refs. [

]. The resolvent response modes of

−

are

then defined as the solutions to the Euler-Lagrange equations given by

−

(11)

Equation (

) constitutes an eigenvalue problem and thus has a countably infinite set of solutions

which we index by the subscript

−

(12)

We have denoted the eigenvalue

−

and the eigenfunctions

such that the singular values and

resolvent response modes of

are given by

and

, respectively. The resolvent forcing modes

are recovered through

(13)

Note that the

are guaranteed to be orthogonal w.r.t.

since the matrices in Eq. (

)are

Hermitian, and the

are orthogonal w.r.t.

, since

−

(14)

Other researchers have used variational based approaches to the study of forced linear dynamics

and we would like to comment briefly on how Eqs. (

) and (

) differ from these past studies.

Monokrousos

et al.

[

] formulated a variational problem to identify the optimal forcing structure,

which as formulated therein necessitated the introduction of an adjoint equation, which is avoided in

the current formulation. Sipp and Marquet [

] solved a generalized Eigenvalue problem to compute

013905-4

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

the resolvent forcing modes and then recovered the response modes by solving the linear system in

Eq. (

). Here we avoid any operator inversion by first computing the response modes and then

recovering the forcing modes afterwards. More fundamentally, we view the variational definition

proposed herein as a theoretical framework rather than simply a computational strategy. In Sec.

II C

we present one possible technique for how this framework can be used in practice.

B. Proof of equivalence

We will now illustrate the equivalence of Eq. (

) to the standard SVD-based definition. For

simplicity we consider the case where

. Again, following the notation of Herrmann

et al.

[

], the SVD of the properly weighted resolvent operator is given by

≡

FHF

−

(15)

The physical resolvent forcing and response modes are then recovered by left multiplication by

−

, such that

−

and

−

, whose columns give the individual modes

and

respectively. We focus first on the resolvent response modes

. Beginning from the definition of the

weighted resolvent we can write

(16)

Next we use Eqs. (

), (

), and (

) to write the above expression in terms of the linear operator

(

−

)

(17)

Taking the inverse of both sides and noting the unitary nature of

we find

−

(

)

−

(18)

Finally, we right multiply by

and left multiply by

to arrive at

(

)

−

(19)

which is equivalent to Eqs. (

). Again, the resolvent forcing modes are then recovered through

(20)

This establishes the equivalence of the variational and SVD-based definitions of resolvent modes.

We would like to emphasize that a consequence of this equivalence is that the completeness property

of the SVD-based basis also applies to the variational computed basis. The case where

follows similar arguments but for the sake of brevity is not included here.

C. Resolvent mode estimation

In general, the Euler-Lagrange Eqs. (

) are both analytically intractable and computationally

intensive for complex flows with multiple nonhomogeneous spatial dimensions. However, the

variational definition introduced here provides a convenient way to estimate resolvent modes as

an expansion in any convenient basis. Suppose we wish to estimate the resolvent response modes,

(

), of some system defined on a particular domain. Then let

(

) with (

...

)besome

known basis defined on that same domain. We can then write the resolvent response modes as an

expansion in this basis:

(21)

Inserting this expansion into Eq. (

) transforms the continuous vector field

∈

∞

into a discrete

field

∈

, where

is the vector of amplitudes

. The Euler-Lagrange Eq. (

) then takes the

form

−

(22)

013905-5

BARTHEL, GOMEZ, AND MCKEON

where

∈

≡

, and

≡

. The eigenvectors

contain the

amplitudes

which optimally approximate the resolvent response modes in the known basis

and the

are the associated optimally approximated singular values. For

basis elements we will

have

∈

, and thus we will obtain

eigenvalue

eigenvector pairs, representing

singular

mode

singular value pairs. The necessary

depends on both the efficiency of the model basis and

the desired level of accuracy. Since Eq. (

) does not include any inherent approximations, if the

input basis is orthogonal the VRA approximation does converge to the exact solution computed

via SVD as

→∞

. However, any Galerkin type method, such as the one proposed here, is only

valuable as long as the size of the basis,

, is significantly smaller than the size of the discretized

system,

.If

∼

(

), then it would be preferable to compute the SVD directly, since one would

not be restricted to the span of the input basis, which if poorly chosen, may not accurately model

the true resolvent modes. However, we show in the following examples that for large systems a

reduction of order

of two orders of magnitude is possible.

Throughout this work we use lower-dimensional resolvent modes as a modeling basis. In such

cases the input basis elements

are periodic in the wall parallel direction and generally localized

around a critical layer (where the wave speed

is equal to the local mean velocity) in the wall normal

direction. For example, we might have

(

)

∼

(

;

)

ikx

where

is the imposed wave number.

Thus, for this type of basis the number of retained spatial wave numbers determines the wall parallel

resolution of the VRA model and the number and range of retained wave speeds determines the wall

normal span of the VRA reconstruction.

However, we note that other types of modeling basis are possible, for example Towne [

]used

a data driven basis to estimate global resolvent modes in turbulent jets. While the derivation of

the model presented in Towne [

], based on successive applications of the CFL principle, differs

from the direct variational derivation presented here, the final eigenvalue problem being solved is

mathematically equivalent to Eq. (

). His work, coined

empirical resolvent mode decomposition

aims to estimate resolvent modes that are constrained to lie in the span of the input basis, computed

from the full flow field data. This constraint is intended to ensure that the modes are

physical

in the

sense that they reflect structures observed in the real flow [

]. The variational approach described

in this work is, however, an alternative mathematical definition of the resolvent modes without any

reference to external data.

III. 1D RESOLVENT ANALYSIS: TURBULENT CHANNEL FLOW

A. The Orr-Sommerfeld squire system

As a first example we consider the incompressible linearized NSE for streamwise constant

fluctuations about a turbulent mean in a wall-bounded shear flow. This example illustrates the

fundamental theory and highlights the analytical manipulation enabled by the VRA framework.

The equations are nondimensionalized using the channel half-height and friction velocity. A Fourier

transform in the homogeneous spatial directions and time results in a system parametrized by the

Reynolds number,

, and the wave-number triplet,

[

,ω

]. Here

and

denote the wave

numbers in streamwise and spanwise directions, respectively, and

again represents the temporal

frequency. We focus on streamwise constant fluctuations which are useful models of the streamwise

elongated structures known to play a crucial role in the sustenance of turbulence [

]. Therefore,

for the remainder of Sec.

III

we assume

The forced linearized NSE can be written as

(

)

(

)

(

)

(

)

(23)

Here

∈

[

−

1] and [

] are the wall-normal and streamwise velocity fluctuations about the

streamwise, spanwise, and temporal averaged mean velocity

. The spanwise velocity is recovered

through the continuity equation as

−

. The right-hand side [

]

represents an unknown

forcing. The relevant boundary conditions are thus

(

0. Note that we

013905-6

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

write Eq. (

) in terms of

instead of the classical formulation in terms of the wall normal vorticity

≡

−

. This is because if

then

∼

. Note that this implies that the off-diagonal

term in Eq. (

) does not include the

present in more classical formulations. The Orr-Sommerfeld

and Squire operators in Eq. (

) simplify to

=−

∇

−

∇

(24)

=−

−

∇

(25)

where

∇

≡

∂

−

. The inner product defining the kinetic energy norm is

≡

∗

−

∗

(26)

where

(

)

≡

−

(

)

[

], and

√

. It is convenient to also define the

following norm associated with the OS operator induced by

≡

∗

−

∗

(27)

which represents the contribution of

(and thus

) to the kinetic energy and where again the norm

is defined as

√

. Last, it is numerically convenient to implement Eq. (

)as

∇

−

∇

−

(28)

To compare our variational results to the direct SVD we use the definition Eq. (

) going forward.

B. The Orr-Sommerfeld and Squire families

It is instructive to decompose the system into the Orr-Sommerfeld (OS) and Squire (SQ) families

of modes as suggested by Rosenberg and McKeon [

]. The OS family corresponds to the forced

response due to

∇

−

(29)

which upon elimination of ̃

from the equation for

results in a decoupled system reminiscent

of the classical OS

SQ decomposition of linear stability theory [

∇

−

(30)

=−

(31)

The SQ family of modes, however, is the forced response to

, where by construction

(32)

Since Eq. (

) is a normal scalar operator, the resolvent forcing and response modes are proportional

to the eigenmodes of

, and the singular values are equal to the inverse of the norm of the

eigenvalues of

(

)

sin

(

(33)

(

)

arctan

−

sin

(

(34)

−

(35)

013905-7

BARTHEL, GOMEZ, AND MCKEON

The problem thus reduces to finding the OS family of modes associated with Eq. (

), which, in

accordance with Sec.

, are defined as the stationary points of the associated Lagrangian

=∇

−

(36)

where

≡

[

] and we have made use of the fact that

0 to simplify the operator norm.

To eliminate the streamwise velocity

we expand the solution to Eq. (

) in eigenfunctions of

given by Eq. (

=−

(37)

This allows us to write the kinetic energy constraint as

−

(

)

(38)

where the third term

(

)

is given by the square of Eq. (

). This allows us to rewrite Eq. (

)as

=∇

−

(39)

with associated Euler-Lagrange equation

∇

−

(40)

For

0, the eigenfunctions of

may also be derived analytically [

]. Using standard

methods they are found to be

(

;

)

{

cos

[

(

)

]

−

cosh

[

(

)

]

}

sin

[

(

)

]

−

sinh

[

(

)

]

(41)

−

ω,

(42)

where

, and

are defined in Appendix

and satisfy

∇

and

Expanding the solution to Eq. (

) in the basis of OS eigenfunctions Eq. (

) such that

(43)

allows us to transform the variation into an optimization over the coefficients

∂

∇

−

∂

−

(44)

Here the quantity

represents the projection of the OS eigenfunctions onto the SQ eigenfunctions

through Eq. (

) such that

≡−

(45)

Upon carrying out the above differentiation with respect to

we find the eigenvalue problem

−

(

)

(46)

where

. The eigenvectors

correspond to the the coefficients which optimally

represent the resolvent response modes of the system Eq. (

) as a linear combination of the

eigenbasis Eq. (

(47)

013905-8

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

FIG. 1. Real part of the wall normal component

(a), streamwise component

(b), and forcing

(c)ofthe

of the first, third, and fifth Orr-Sommerfeld family of resolvent modes. Reference modes computed via direct

SVD are shown in solid lines, VRA reconstruction using 20 basis eigenmodes is shown in symbols.

1, and

1000.

The singular values

are given by the eigenvalues of Eq. (

) and the forcing modes are recovered

through

∇

−

(48)

Together with the Squire family of resolvent modes Eqs. (

)–(

) the Orr-Sommerfeld family

given by Eqs. (

) and (

) fully describe the resolvent basis. In Fig.

we plot the real part of the

variationally reconstructed Orr-Sommerfeld response and forcing modes along side their numeri-

cally computed counterparts for the wave-number triplet [

,ω

]

1] and

1000. The

singular values plots are plotted in Fig.

2(a)

. For this example, the VRA model uses

basis elements, this value is chosen to show a balance between the accuracy and model reduction

capabilities of the method. Although for this example the computational cost is trivial, the reduction

in size of the relevant matrices and avoiding the need for matrix inversion reduces the computation

time by two orders of magnitude. To quantify the convergence of our method we plot in Fig.

the error in the VRA reconstruction of

, and

as a function of the number of retained OS

eigenfunctions (

) included in the variational reconstruction. The error is defined as

−

svd

−

vra

(49)

where

and the subscripts

vra

and

svd

denote the quantities computed using the VRA model

and direct SVD, respectively. In all cases we observe monotonic convergence. In this this example

the VRA model is extremely effective at reconstructing the results of the direct SVD since our model

basis exactly spans the range of

C. Analytical approximation of

In this section we demonstrate how, under certain assumptions, the variational resolvent for-

mulation allows for the analytical approximation of the leading OS resolvent mode

. Written

explicitly, the Lagrangian associated with Eq. (

)is

(

)

(

∇

)

(

∇

)

(

∇

)

(

∇

)

−

(

)

−

(50)

013905-9

BARTHEL, GOMEZ, AND MCKEON

FIG. 2. Orr-Sommerfeld family of singular values (a) with reference values computed via direct SVD in

red, and variational reconstruction using 20 basis eigenmodes in black. Error in variational reconstruction as a

function of basis elements in

(b),

(c), and

(d) of

for

1, and

1000

as a function of the retained basis elements

where

is the solution to

−

∇

=−

(51)

Here

and

are the streamwise and wall normal components of

and

is the leading OS

singular value. The associated Euler-Lagrange equation written in terms of

is then

∇

−

∇

(52)

with boundary conditions

(

0. Note that we use the the original defini-

tion Eq. (

) not the numerical implementation Eq. (

) to derive Eq. (

). This is done to avoid the

analytically cumbersome treatment of the

∇

−

operator. The problem is now parameterized by

, and

. Our analysis will consider the appropriate limits of each in turn.

It has been shown that for

0 the most linearly amplified frequency is

0, therefore we

will consider the limit as

→

0. Since, in this limit Eq. (

) is regularly perturbed problem, the

leading order solution may be found by simply setting

0. We may further simplify Eq. (

)

by considering a high Reynolds number limit

→∞

. Analysis of Eq. (

) reveals that for

∼

→∞

(see Appendix

). This allows us introduce the small parameter

≡

−

such

that Eqs. (

) and (

) take the form

∇

(

∇

−

)

∇

−

∇

(53)

∇

(54)

where ̃

(

). We note that Eq. (

) implies that

∼

and expand the solution in an asymptotic

series:

(

)

(

)

(55)

013905-10

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

The leading order solution to Eqs. (

) and (

) then satisfy

∇

[(

∇

−

)

∇

−

]

(56)

∇

(57)

and the norm constraint takes the form

(58)

In this work we focus on the leading order solution, and thus to avoid notational clutter we drop the

subscripts

and

moving forward.

While we have managed to simplify the governing equations, the second term in Eq. (

) remains

prohibitive to analytical progress. To proceed we consider the

→−

symmetry of Eq. (

) which

dictates that the resolvent modes come in pairs, one of which is even about the center of the channel,

and one of which is odd. If additionally, the modes have compact support, as is generally the case,

then we have

(

)

(

−

)

(

)

=−

(

−

), and therefore it is sufficient to solve for the

mode shape in one half of the domain.

We assume that

is indeed locally supported and thus introduce the scaling

under

the assumption

1 and make the transformation

(

)

(

)

→

(

)

(

). We note that this

scaling differs from the

and

scaling derived by Arratia and Chomaz [

] in the context

of inviscid transient growth. If we formally take the limit as

→∞

, then we may transform the

domain from

∈

[

−

1], to the “semi-infinite” half channel:

∈

∞

] and recover the solution

in the other half through

(

)

(

−

)

(

)

=−

(

−

Finally, to make progress we require some suitable approximation of the mean velocity profile.

Since we are working within a high Reynolds number limit we choose to assume that the mean

velocity obeys a logarithmic profile over the entirety of the semi-infinite domain. This is a reasonable

assumption since in high Reynolds number channel flow the log-law applies to a large fraction of

the channel. Our approach thus implicitly assumes the support of the resolvent modes is localized

within this region where the log-law approximation is valid. The mean shear is then given in our

scaled variables by

(

)

−

, where

is the Von Karman constant. We note that the mean

shear diverges as like

−

→

0, however, since

(0)

0wehave

(

)

∼

→

0, and thus the right-hand side of Eq. (

) remains bounded as

→

Inspection of Eqs. (

) and (

) reveals that the appropriate scaling of the velocity components

is given by

(

)

(

) and

(

)

(

). Additionally, we define the scaled Laplacian

∇

≡

∂

−

1 such that

∇

→

∇

, and note that for

0 and

→∞

the singular value

scales as ̃

∼

−

(see Appendix

). Thus, we can write Eq. (

) in our scaled variables as

∇

[(

∇

−

)

∇

−

]

(59)

where

is a constant. We expand

and

in asymptotic series

−

(60)

which upon substitution into Eq. (

) allows us to eliminate the norm constraint at leading order

and reduce the problem of deriving the leading OS resolvent mode to

∇

(61)

∇

(62)

where we have again dropped the subscripts to avoid notational clutter. The relevant boundary con-

ditions are

(0)

(

∞

)

(

∞

)

0. The remaining constants of integration

013905-11

BARTHEL, GOMEZ, AND MCKEON

are then chosen such that

∇

−

is minimized and

1. Here we choose to minimize

∇

−

instead of

to facilitate comparison with the numerically computed modes.

However, we note that minimizing the latter functional leads to a very similar solution. Using

standard methods the solutions satisfying the boundary conditions are found to be

(

)

(

)

−

(63)

(

)

=−

)

)(

1)]

−

(64)

The three remaining constants of integration,

, are found by minimizing

∇

−

subject

to the constraint

1. Straightforward integration results in

∇

−

∇

)

(65)

and

128

(112

228

)

351

1089

128

(66)

Minimizing Eq. (

) subject to Eq. (

) results in the eigenvalue problem

⎡

⎣

12 0 0

01836

0 36 144

⎤

⎦

⎡

⎣

⎤

⎦

⎡

⎣

16 7

831

16 351

4 351

64 1089

⎤

⎦

⎡

⎣

⎤

⎦

(67)

Assuming

4 the minimizing solution that satisfies the norm constraint is found to be

[

]

1283

1066

0431]

√

(68)

The leading singular value is

=∇

−

(69)

The wall-normal component

of the optimal resolvent forcing mode

is recovered through

∇

(70)

subject to the boundary conditions

(

0. Using Eq. (

) and letting

(

)

→

(

)this

takes the form

∇

(

)

=−

∇

(

)

(71)

The solution satisfying the boundary condition

(0)

0 is found to be

(

)

−

)

−

)

]

−

(72)

The solutions Eqs. (

), (

), and (

) with optimal coefficients Eq. (

) are plotted in Fig.

alongside numerically computed resolvent modes for

000 and

0 over a range of

Note that for

0, and

→

0 the symmetries of Eq. (

) result in numerical resolvent modes

with constant arbitrary phase, which for ease of comparison we set to zero. With the exception of

the

component for the smallest wave number (

6), the derived scaling laws lead to reasonable

collapse in both the numerically computed resolvent response and forcing modes. As

→

1the

assumption of local support in

is no longer valid. In this limit

tends to have significant support

at the channel center.

For the response modes the analytically-derived mode accurately predicts the shape, amplitude,

and localization of the numerically computed modes. The analytical prediction of the wall normal

013905-12

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

FIG. 3. Optimal resolvent modes:

(a),

(b),

(c), and singular value (d) for

10 000

and a range of

. Numerically, calculated modes

singular values are shown in colored lines

dots, analytically

derived modes

singular values are shown in black open circles

dashed line. From light to dark, colors indicate

increasing

from 6 to 100 (a–c).

velocity is most accurate for the largest wave numbers, tending to slightly over predict the amplitude

of the smaller wave-number modes. This is most likely due to the fact that the amplitude of

smaller by a factor of

000 and is thus susceptible to some numerical uncertainty since

it does not meaningfully contribute to the norm. The streamwise velocity more closely obeys the

derived scaling laws, and thus the analytical model accurately predicts the shape of the numerically

computed modes for all

The prediction of the forcing mode is slightly less accurate. While we capture the location and

amplitude of the peak, the model underpredicts the true mode closer to the wall. The discrepancy in

the forcing despite accurate reconstruction of the response is due to the sensitivity of the action of

linear operator

to perturbations in the argument

. This is discussed in detail in Sec.

Finally, in Fig.

we also plot the numerically computed leading singular values along side the

analytical prediction Eq. (

). While the analytically obtained value of

slightly under-predicts

the true singular values for the smaller values of

, the numerical singular values do converge to

the analytical prediction with increasing

, consistent with the assumption made in our model that

1. This under prediction is consistent with the fact that the true singular value represents the

global maximum gain.

IV. 2D RESOLVENT ANALYSIS: PERIODIC MEAN FLOW

In this section we use VRA to efficiently and accurately compute resolvent modes about a 2D

mean flow. We consider the equilibrium solution EQ1 found in plane Couette flow by Nagata

[

]. The data was obtained from the open-source database

channelflow.org

[

]. In this case

the flow has two nonhomogeneous spatial dimensions, the wall normal direction

∈

[

−

1] and

the spanwise direction

∈

[

−

2] with

. The spanwise periodic EQ1 solution is

shown in Fig.

. The 2D

3C resolvent modes computed about this flow are then parameterized by

the streamwise wave number and frequency pair [

,ω

]. We choose as our modeling basis the local

1D resolvent modes about the mean flow

(

) given by the spanwise average of the EQ1 solution:

(

)

(

;

)

. In other words, we seek to approximate the 2D

3C resolvent modes

013905-13

BARTHEL, GOMEZ, AND MCKEON

FIG. 4. Exact coherent state EQ1 at

400 used to compute 2D

3C resolvent modes:

(

)(a),

(

)

(b),

(

)(c) and spanwise average

(

)(d) used to compute the 1D basis modes.

from the 1D resolvent basis as

(

;

)

(

)

(73)

The expansion coefficients

are found by solving the eigenvalue problem

−

(74)

where

and

. The operator

is the NS operator, in velocity-

vorticity form, linearized about the 2D

3C mean flow, the details of which are discussed in

Rosenberg and McKeon [

]. The operator is discretized in

33 Chebychev points in the

wall normal direction, and

32 linearly spaced points in the spanwise direction, for a total

2112 degrees of freedom.

The 1D resolvent modes are computed for the same

as the 2D modes, and a range of

linearly spaced wave speeds 0

where

ω/

. We use a range of

since the

2D mode is expected to be localized near but not necessarily exactly at the critical layer where

(

). To account for the variation in

we include a range of

11 spanwise wave numbers

[

−

...

π/

. We found that increasing the number of retained harmonics beyond

this range did not meaningfully change the results. At each wave-number triplet [

]we

include

SVD

8 resolvent modes, resulting in a total of

SVD

254 degrees of

freedom. These values were chosen to demonstrate a balance between accuracy and the cost-saving

potential of the proposed method. (The reader is referred to Appendix

for an illustration of some

representative basis elements.) Once

is known, the construction of the matrices

and

takes approximately 0.5 s and the associated eigendecomposition takes approximately 0.01 s on

a personal laptop. Meanwhile, the inversion and direct truncated SVD of the original system takes

approximately 5 s using the built in Matlab functions

mldi

ide

() and

().

In Figs.

and

we compare the real part of the first four resolvent response modes of the

variational reconstruction and the modes computed directly through the SVD of the 2D resolvent

for

5 and

75 and

400. The variational approach very accurately reconstructs the

true response modes considering the significant reduction in computational complexity.

013905-14

VARIATIONAL FORMULATION OF RESOLVENT ANALYSIS

FIG. 5. Real part of the

component of the first 4 resolvent response modes (

)for

75,

and

400. Top row: true modes, bottom row: VRA model with

11,

3, and

SVD

8. From left

to right:

In Figs.

and

we plot resolvent forcing modes computed from the response modes through

. Interestingly we find that while the

component is reproduced accurately the

component shows significant discrepancy. While the qualitative shape of the

component of the

FIG. 6. Real part of the

component of the first four resolvent response modes (

)for

75,

and

400. Top row: true modes, bottom row: VRA model with

11,

3, and

SVD

8. From left

to right:

013905-15