2109.08000.pdf

A Variational Formulation of Resolvent Analysis

Benedikt Barthel,

∗

Salvador Gomez, and Beverley J. McKeon

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

(Dated: September 28, 2021)

The conceptual picture underlying resolvent analysis(RA) is that the nonlinear term in the Navier-

Stokes(NS) 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 trans-

fer 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 advan-

tageous. 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 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 modes. Second, to illustrate the cost

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

a 2-dimensional, 3-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.

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.

[1] and Jovanovi ́c and Bamieh [2] 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 [3] 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 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 non-homogeneous spatial

dimension for which cost of the inversion and SVD of the operator is trivial [2–6]. 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 [7]. In particular, resolvent modes have been found to exhibit self similar behaviour

characteristic of the attached eddy hypothesis proposed by Townsend [5, 8, 9]. Towne

et al.

[10] 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 [11, 12], the flow behind bluff bodies [13, 14], exact

coherent states (ECS) [15], and turbulent jets [16, 17]. 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

∗

bbarthel@caltech.edu.

arXiv:2109.08000v2 [physics.flu-dyn] 25 Sep 2021

[18], and delay flow separation [19]. 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.

[20] 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.

[21] 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 [11, 16]. 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 [5, 22–25].

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 equivalent definition based on an extension of

the Courant–Fischer–Weyl min-max principle (CFL). The CFL principle itself has been used previously by Dawson

and McKeon [26] who formulated a simplified variational problem to estimate the shape of the vorticity component

of the optimal resolvent mode in wall bounded shear flows. We believe the 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, to be novel in the resolvent literature.

This new 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. While in general the resulting

Euler-Lagrange equations remain difficult to solve exactly, this variational formulation allows for the approximation

of resolvent modes as an expansion in any convenient basis, for example the much cheaper one-dimensional resolvent

basis in a two- or three-dimensional problem, an analytical basis such as that described by [26] 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, and to investigate the latter we perform RA around

a 2D/3C exact coherent solution. 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

§

II we derive the proposed variational definition. In

§

III we use the variational

formulation to analyze streamwise constant structures in turbulent channel flow. In

§

IV and

§

V we consider RA applied

to both streamwise periodic and streamwise developing two-dimensional, three velocity component (2D/3C) systems

to illustrate the computational cost and memory saving potential of the proposed VRA formulation. In

§

VI 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

§

VII and conclude in

§

VIII.

II. MATHEMATICAL FORMULATION

Let us consider a general forced linear system

∂

∂t

−

(1)

where

represents a spatial-linear differential operator and

(

)

(

)

∈

∞

. The state variables

and

are

referred to as the ‘response’ and ‘forcing’ respectively. We consider the temporal Fourier transfer of (1) and define

the spatio-temporal linear operator

(

)

≡

iω

−

(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.

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

‖

= 1

(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

§

III. Following the notation of Herrmann

et al.

[21] 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(6), (7) and the definition (9) 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 (10).

The reader is referred to appendix A for a derivation of this property based on the work of [27, 28]. The resolvent

response modes of

−

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

−

= 0

(11)

Equation 11 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 since the matrices in (12) are Hermitian, and the

are orthogonal

w.r.t

since

−

(14)

B. Proof of equivalence

We will now illustrate the equivalence of (12) to the standard SVD-based definition. For simplicity we consider the

case where

‖

. Again, following the notation of Herrmann

et al.

[21], 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 (3), (8), and (15) 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

(

)

QΨΣ

−

(19)

which is equivalent to (12). Again, the resolvent forcing modes are then recovered through

LΨΣ

(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 equations (12) are both analytically intractable and computationally intensive for

complex flows with multiple non-homogeneous spatial dimensions. However, the variational definition introduced

above 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

(

), and let

(

) with (

= 1

...r

) be some known

basis defined on the same domain. We can then write the resolvent response modes as an expansion in this basis.

(21)

Inserting this expansion into (10) transforms the continuous vector field

∈

∞

into a discrete field

∈

, where

is the vector of amplitudes

. The Euler-Lagrange equation (11) takes the form

−

= 0

(22)

where

∈

≡

, and

≡

. The eigenvectors

contain the amplitudes

which

optimally approximate the resolvent response modes and the

are the approximate 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. However, we show in the following examples that for large systems a reduction over the dimension of the

original system by up to two orders of magnitudes is possible due to the lack of matrix inversion. Throughout the

paper we use

to refer to the size of the reduced system (22) and

to refer to the size of the original system.

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 wavenumbers 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 [29]. Therefore, for the remainder of

§

III we assume

= 0.

The forced linearized NSE can be written as

[

][

(

)

(

)

]

[

(

)

(

)

]

(23)

Here

∈

[

−

1] and [

v,u

] 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

(

1) =

(

1) =

(

1) = 0. Note that we write (23) in terms of

instead of the classical formulation in terms of

the wall normal vorticity

−

, since if

= 0 then

∼

. Note that this implies that the off-diagonal term

in (23) does not include the

present in more classical formulations. The Orr-Sommerfeld and Squire operators in

(23) simplify to

−

iω

∇

−

∇

(24)

−

iω

−

∇

(25)

where

∇

≡

∂

−

. The inner product defining the kinetic energy norm is

〈

〉

≡〈

∗

−

∗

i,y

j,y

∗

〉

(26)

where

〈

(

)

〉 ≡

∫

−

(

)

= [

v,u

], and

‖

√

〈

〉

. It is convenient to also define the following norm

associated with the OS operator induced by

〈·〉

≡〈

∗

−

∗

i,y

j,y

〉

(27)

which represents the contribution of

(and thus

) to the kinetic energy and where again the norm is defined as

‖

√

〈

v,v

〉

. Lastly, it is numerically convenient to implement (23) as

[

∇

−

][

]

[

∇

−

]

[

]

(28)

In order to compare our variational results to the direct SVD we use the definition (28) 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 [15]. 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 [22, 30].

∇

−

= ̃

(30)

−

(31)

The SQ family of modes, on the other hand, is the forced response to

, where by construction

= 0.

(32)

Since (32) is a normal scalar operator, the resolvent forcing and response modes are proportional to the eigenmodes

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

(

) =

[

sin

(

jπ

(

+ 1)

)]

(33)

(

) =

[

arctan

(

−

Rω

)

sin

(

jπ

(

+ 1)

)

]

(34)

(

+ 4

)

−

(35)

The problem thus reduces to finding the OS family of modes associated with (29), which in accordance with

§

II,

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. In order to

eliminate the streamwise velocity

we expand the solution to (31) in eigenfunctions of

given by (33).

−

〈

〉

(37)

This allows us to write the kinetic energy constraint as

‖

〈|

−

(

)

〉

‖

(38)

where the third term

(

)

is given by the square of (37). This allows us to rewrite (36) as

(

)

‖∇

−

‖

(39)

with associated Euler-Lagrange equation

δv

(

‖∇

−

‖

)

= 0

(40)

For

= 0, the eigenfunctions of

may also be derived analytically [2, 31]. Using standard methods they are

found to be

(

;

) =

[cos (

(

+ 1))

−

cosh (

(

+ 1))] +

[

sin (

(

+ 1))

−

sinh (

(

+ 1))

]

(41)

(

)

−

iω

(42)

where

and

are defined in appendix B and satisfy

∇

and

〈

〉

. Expanding the

solution to (40) in the basis of OS eigenfunctions (41) such that

(43)

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

∂

∂a

(

‖∇

−

‖

−

(

‖

−

)

∂

∂a

(

−

(

)

= 0

(44)

Here the quantity

represents the projection of the OS eigenfunctions onto the SQ eigenfunctions through (37).

≡−

〈

〉

(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 (29) as a linear combination of the eigenbasis (41).

= [

(

)]

(47)

The singular values

are given by the eigenvalues of (46) and the forcing modes are recovered through

= [

∇

−

0] = [

(48)

Together with Squire family of resolvent modes (33-35) the Orr-Sommerfeld family given by (47) and (48) fully describe

the resolvent basis. In Figures 1 we plot the real part of the variationally reconstructed Orr-Sommerfeld response and

forcing modes along side their numerically computed counterparts for the wave number triplet [

,ω

] = [0

and

= 1000. The singular values plots are plotted in 2a. For this example, the VRA model uses

= 20

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 Figure 2 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

√

∫

−

svd

−

var

(49)

where

u,v

and the subscripts

var

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

1

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

analytical approximation of the leading OS resolvent mode

. Written explicitly, the Lagrangian associated with

(23) is

(

) =

(

∇

)

(

∇

)

(

∇

)

(

∇

)

−

(

)

−

)

(50)

where

is the solution to

−

iωu

−

∇

−

(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

(

∇

)

(

−

)

−

∇

)

= 0

(52)

with boundary conditions

(

1) =

(

1) =

(

1) = 0. Note that we use the the original definition (23) not the

numerical implementation (28) to derive (52). 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.

FIG. 1: Real part of the wall normal component

(a), streamwise component

(b), and forcing

, 3

, and 5

, Orr-Sommerfeld family of resolvent modes. Real part (left column), imaginary part (right column).

Reference modes computed via direct SVD are shown in solid lines, VRA reconstruction using 20 basis eigenmodes

is shown in symbols.

= 6,

= 0

1, and

= 1000.

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

= 6,

= 0

1, and

= 1000 as a function of

the retained basis elements

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 (52) is regularly perturbed problem, the leading order solution may be found by

simply setting

= 0. We may further simplify (52) by considering a high Reynolds number limit

→ ∞

. Analysis

of (29) reveals that for

= 0,

∼

→ ∞

(see appendix C). This allows us introduce the small parameter

≡

−

such that (52) and (51) take the form

∇

(

∇

−

)

∇

−

∇

)

= 0

(53)

∇

(54)

where ̃

(

). We note that (54) implies that

∼

and expand the solution in an asymptotic series.

(

)

(

)

(55)

The leading order solution to (53) and (54) then satisfy

∇

(

∇

−

)

∇

−

)

= 0

(56)

∇

(57)

and the norm constraint takes the form

‖

= 1

(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 eliminate the nonlinearity, the second term in (56) remains prohibitive to analytical

progress. In order to proceed we consider the

→−

symmetry of (23) 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, 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

(

)

(

)

→

(

)

(

). If we formally take the limit as

→ ∞

we may

transform the domain from

∈

[

−

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

∈

∞

] and recover the solution in the

other half through

(

) =

(

−

) =

(

) =

−

(

−

Finally, in order 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

(

κY

)

−

, where

is the Von Karman constant. We note that the mean

shear diverges as like

−

→

0, however, since

(0) =

(0) = 0 we have

(

)

∼

→

0, and thus the

right hand side of (57) remains bounded as

→

Inspection of (57) and (58) reveals that the appropriate scaling of the velocity components is given by

(

) =

(

) and

(

) =

(

). Additionally, we define the scaled Laplacian

∇

≡

∂

Y Y

−

1 such that

∇

→

∇

and note that for

= 0 and

→∞

the singular value scales as ̃

∼

−

(see appendix C). Thus we can write

(56) in our scaled variables as

∇

(

∇

−

)

∇

−

)

= 0

(59)

where

is a constant. We expand

and

in asymptotic series

−

(

−

)

−

(

−

)

(60)

which upon substitution into (59) allows us to eliminate the norm constraint at leading order and reduce the problem

of deriving the leading OS resolvent mode to

∇

= 0

(61)

∇

κY

(62)

where we have again dropped the subscripts to avoid notational clutter. The relevant boundary conditions are

(0) =

(

∞

) =

(

∞

) = 0. The remaining constants of integration are then chosen such that

‖∇

−

‖

is minimized and

‖

= 1. Here we choose to minimize

‖∇

−

‖

instead of

‖

order 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)

(

) =

−

(

+ (4

+ 6

)

+ (6

+ 6

+ 9

) (

+ 1)

)

Y e

−

(64)

The three remaining constants of integration,

a,b,c

, are found by minimizing

‖∇

−

‖

subject to the constraint

‖

= 1. Straight forward integration results in

‖∇

−

‖

‖∇

‖

= 4

(

+ 9

+ 36

+ 72

)

(65)

and

‖

(

128

(112

+ 228

)

351

1089

128

)

= 1

(66)

Minimizing (65) subject to (66) results in the eigenvalue problem





12 0

0 18 36

0 36 144













(

)





8 31

16 351

4 351

64 1089













(67)

Assuming

= 0

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

[

a,b,c

] = [0

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

(

1) = 0. Using (69) and letting

(

)

→

(

) this takes the form

∇

(

) =

−

∇

(

)

(71)

The solution satisfying the boundary condition

(0) = 0 is found to be

(

) =

(

+ (3

−

)

+ (2

−

)

−

(72)

FIG. 3: Optimal resolvent modes:

(a),

(b),

= 0,

= 0

= 10

000 and a

range of

. Numerically calculated modes shown in colored lines, analytically derived modes are shown in the black

open circles. From light to dark, colors indicate increasing

from 6 to 100.

The solutions (63), (64), and (72) with optimal coefficients (68) are plotted in Figure 3 alongside numerically

computed resolvent modes for

= 10

000 and

= 0 over a range of

. Note that for

= 0, and

→

0 the

symmetries of (28) 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

→

1 the

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 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

is smaller by a factor of

= 10

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

§

VI.

Finally, in Figure 3 we also plot the numerically computed leading singular values along side the analytical prediction

(69). While the analytically obtained value of

slightly under-predicts the true singular values for the smaller values

, 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/3C mean flow. We

consider the equilibrium solution EQ1 found in plane Couette flow by Nagata [32]. The data was obtained from the

open-source database

channelflow.org

[33, 34]. In this case the flow has two non-homogeneous spatial dimensions, the

wall normal direction

∈

[

−

1] and the spanwise direction

∈

[

−

2] with

= 0

. The spanwise periodic

EQ1 solution is shown in Figure 4. The 2D/3C resolvent modes computed about this flow are then parameterized

by the streamwise wavenumber 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:

(

y,z

) =

(

;

)

In other words, we seek to approximate the 2D/3C resolvent modes from the 1D resolvent basis as

(

y,z

;

) =

(

y,z

)

(73)

The expansion coefficients

are found by solving the eigenvalue problem

−

= 0

(74)

where

i,j

〈

〉

, and

i,j

〈

〉

. 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 [35]. 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 of

= 2

= 2112 degrees of freedom.

The 1D resolvent modes are computed for the same

as the 2D modes, and a range of

= 3 linearly spaced

wavespeeds 0

≤

where

ω/k

. 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 wavenumbers

= [

−

...

π/L

. 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

SV D

= 8 resolvent modes, resulting in a total of

SV D

= 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 D for an illustration of some representative basis elements). Once

is known,

the construction of the matrices

and

takes approximately 0.5 seconds and the associated eigendecomposition

takes approximately 0.01 seconds on a personal laptop. Meanwhile, the inversion and direct truncated SVD of the

original system takes approximately 5 seconds using the built in Matlab functions

mldivide

() and

svds

().

In Figures 5 and 6 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

= 0

5 and

= 0

75 and

= 400. The

variational approach very accurately reconstructs the true response modes considering the significant reduction in

computational complexity.

In Figures 7 and 8 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 forcing mode is predicted by the VRA model, the

mode is contaminated by higher harmonics. This contamination observed in the VRA reconstruction of the forcing

modes,

, despite the accurate reconstruction of the response modes,

, is due to the directional amplification of

the resolvent operator or equivalently, a sensitivity of the action of the linear operator

, to perturbations in the

input

. This phenomenon is discussed in detail in

§

VI.

Additionally, in Figure 9a we compare the variationally computed singular values with the true values computed via

direct SVD. The singular values are estimated relatively accurately, with our model tending to slightly underestimate

the leading singular values. As before, the true singular values represent the optimal gains and the predicted singular

values are bounded above by the true values. For this example there is no significant separation of singular values, in

other words the resolvent operator is not low rank, and yet our method still accurately predicts the singular values

and resolvent response modes.

In order to quantify the convergence properties of the proposed method, for this example we fix

, include

= [

−

...

π/L

such that

= 1 and

= 11 and compute the error as a function of the number of

retained singular modes

SV D

. The error is based on the kinetic energy norm and is defined as

≡

√

∫

−

svd

−

vra

dydz

(75)

where

= [

u,v,w

]. The error is plotted in Figure 9 alongside the relative error in singular values for two values of

the wave speed,

= 0

75 and

= 0. The former corresponds to the example plotted in Figures 5 through 9a where

there is no significant singular value separation. The latter case corresponds to a case where the 2D/3C resolvent is

FIG. 4: Exact coherent state EQ1 at

= 400 used to compute 2D/3C resolvent modes. Clockwise from top left:

(

y,z

(

y,z

(

y,z

), and spanwise average

(

) used to compute the 1D basis modes.

more low rank, (

/σ

≈

6). In both cases our method is not only able to accurately approximate the leading singular

mode and value but also a large range of suboptimal modes and singular values. Interestingly, we see that our method

is more accurate in the case where there is less singular value separation. Furthermore, for the low rank case, (

= 0)

the largest error in singular value is for

. Again, these findings are a result of the directional nature of the resolvent

operator and are discussed in detail in

§

VI.