Simple aerodynamic configurations under even modest conditions can exhibit complex flows with a wide range of temporal and spatial features. It has become common practice in the analysis of these flows to look for and extract physically important features, or modes, as a first step in the analysis. This step typically starts with a modal decomposition of an experimental or numerical dataset of the flowfield, or of an operator relevant to the system. We describe herein some of the dominant techniques for accomplishing these modal decompositions and analyses that have seen a surge of activity in recent decades [1 –8]. For a nonexpert, keeping track of recent developments can be daunting, and the intent of this document is to provide an introduction to modal analysis that is accessible to the larger fluid dynamics community. In particular, we present a brief overview of several of the well-established techniques and clearly lay the framework of these methods using familiar linear algebra. The modal analysis techniques covered in this paper include the proper orthogonal decomposition (POD), balanced proper orthogonal decomposition (balanced POD), dynamic mode decomposition (DMD), Koopman analysis, global linear stability analysis, and resolvent analysis.

In the study of fluid mechanics, there can be distinct physical features that are shared across a variety of flows and even over a wide range of parameters such as the Reynolds number and Mach number [9,10]. Examples of common flow features and phenomena include von Kármán shedding [11 –17], Kelvin–Helmholtz instability [18–20], and vortex pairing/merging [21–23]. The fact that these features are often easily recognized through simple visual inspections of the flow, even under the presence of perturbations or variations, provides us with the expectation that the features can be extracted through some mathematical procedure [24]. We can further anticipate that these dominant features provide a means to describe in a low-dimensional form that appears as a complex high-dimensional flow. Moreover, as computational techniques and experimental measurements are advancing their ability in providing large-scale high-fidelity data, the compression of a vast amount of flowfield data to a low-dimensional form is ever more important in studying complex fluid flows and developing models for understanding and modeling their dynamical behavior.

To briefly illustrate these ideas, let us provide a preview of modal decomposition. In Fig. 1, we present a modal decomposition analysis of two-dimensional laminar separated flow over a flat-plate wing [25,26]. By inspecting the flowfield, we clearly observe the formation of a von Kármán vortex street in the wake as the dominant unsteady feature. A modal decomposition method discussed later (proper orthogonal decomposition [1,27,28]; see Sec. III) can extract the important oscillatory modes of this flow. Moreover, two of these most dominant modes and the mean represent (reconstruct) the flowfield very effectively, as shown in the bottom figure. Additional modes can be included to reconstruct the original flow more accurately, but their contributions are much smaller in comparison to the two unsteady modes shown in this example. What is also encouraging is that the modes seen here share a striking resemblance to the dominant modes for three-dimensional turbulent flow at a much higher Reynolds number of 23,000 with a different airfoil and angle of attack (see Sec. III.B.1).

We refer to modal decomposition as a mathematical technique to extract energetically and dynamically important features of fluid flows. The spatial features of the flow are called (spatial) modes, and they are accompanied by characteristic values, representing either the energy content levels or growth rates and frequencies. These modes can be determined from the flowfield data or from the governing equations. We will refer to modal decomposition techniques that take flowfield data as input to the analysis as data-based techniques. This paper will also present modal analysis methods that require a more theoretical framework or discrete operators from the Navier–Stokes equations, and we will refer to them as operator-based techniques.

The origin of this document lies with an AIAA Discussion Group, titled “Modal Decomposition of Aerodynamic Flows,” formed under the auspices of the Fluid Dynamics Technical Committee. One of the initial charters for this group was to organize an invited session where experts in the areas of modal decomposition methods would provide an introductory crash course on the methods. The intended audience for these talks was the nonspecialist, e.g., a new graduate student or early-career researcher who, in one afternoon, could acquire a compact yet intensive introduction to the modal analysis methods. This session (121-FC-5) appeared at the 2016 AIAA Aviation Conference* (13–17 June 2016 in Washington, D.C.) and provided the foundation for the present overview paper.

In this overview paper, we present key modal decomposition and analysis techniques that can be used to study a range of fluid flows. We start by reprising the basics of eigenvalue and singular value decompositions as well as pseudospectral analysis in Sec. II, which serve as the backbone for all decomposition and analysis techniques discussed here. We then present data-based modal decomposition techniques: proper orthogonal decomposition in Sec. III, balanced POD in Sec. IV, and dynamic mode decomposition in Sec. V. These sections are then followed by discussions on operator-based modal analysis techniques. The Koopman analysis is briefly discussed in Sec. VI as a generalization of the DMD analysis to encapsulate nonlinear dynamics using a linear (but infinite-dimensional) operator-based framework. The global linear stability analysis and resolvent analysis are presented in Secs. VII and VIII, respectively. Table 1 provides a brief summary of the techniques to facilitate comparison of the methods before engaging in the details of each method.

For each of the methods presented, we provide subsections on overview, description, illustrative examples, and future outlook. We offer in the Appendix an example of how the flowfield data can be arranged into vector and matrix forms in preparation for performing the (data-based) modal decomposition techniques presented here. At the end of the paper, in Sec. IX, we provide concluding remarks on the modal decomposition and analysis methods.

II. Eigenvalue and Singular Value Decompositions

The decomposition methods presented in this paper are founded on the eigenvalue and singular value decompositions of matrices or operators. In this section, we briefly present some important fundamental properties of the eigenvalue and singular value decomposition techniques. We also briefly discuss the concepts of pseudospectra and nonnormality.

Eigenvalue decomposition is performed on a square matrix, whereas singular value decomposition can be applied on a rectangular matrix. Analyses based on the eigenvalue decomposition are usually employed when the range and domain of the matrix or operator are the same [29]. That is, the operator of interest can take a vector and map it into the same space. Hence, eigenvalue decomposition can help examine the iterative effects of the operator (e.g, $A^{k}$ and $\exp (A t) = I + A t + \frac{1}{2} A^{2} t^{2} + \dots$ ).

The singular value decomposition, on the other hand, is performed on a rectangular matrix, which means that the domain and range spaces are not necessarily the same. As a consequence, singular value decomposition is not associated with analyzing iterative operators. That is, rectangular matrices cannot serve as propagators. However, singular value decomposition can be applied on rectangular data matrices compiled from dynamical processes (see Sec. II.C and the Appendix for details).

The theories and numerical algorithms for eigenvalue and singular value decompositions are not provided here but are discussed extensively in textbooks by Horn and Johnson [30], Golub and Loan [31], Trefethen and Bau [29], and Saad [32]. Numerical programs and libraries to perform eigenvalue and singular value decompositions are listed in Sec. II.D.

A. Eigenvalue Decomposition

The eigenvalues and eigenvectors of a matrix (linear operator) capture the directions in which vectors can grow or shrink. For a given matrix $A \in C^{n \times n}$ , a vector $v \in C^{n}$ and a scalar $λ \in C$ are called an eigenvector and an eigenvalue, respectively, of $A$ if they satisfy

A v = λ v

(1)

Note that the eigenvectors are unique only up to a complex scalar. That is, if

v

is an eigenvector,

α v

is also an eigenvector (where

α \in C

). The eigenvectors obtained from computer programs are commonly normalized such that they have unit magnitude. The set of all eigenvalues† of

A

is called a spectrum of

A

Although the preceding expression in Eq. (1) appears simple, the concept of an eigenvector has great significance in describing the effect of premultiplying $A$ on a vector. The aforementioned expression states that, if an operator $A$ is applied to its eigenvector (eigendirection), the operation can be captured solely by the multiplication of a scalar $λ$ , which is the eigenvalue associated with that direction. The magnitude of the eigenvalue tells us whether the operator $A$ will increase or decrease the size of the original vector in that particular direction. If multiplication by $A$ is performed in an iterative manner, the resulting vector from the compound operations can be predominantly described by the eigenvector having the eigenvalue with the largest magnitude, as shown by the illustration in Fig. 2.

Fig. 2
Collection of random points (vectors $x$ ) stretched in the direction of the dominant eigenvector $v_{1}$ with iterative operations $A^{k}$ for matrix $A$ , which has eigenvalues of $λ_{1} = 1.2$ and $λ_{2} = 0.5$ .

If $A$ has $n$ linearly independent eigenvectors $v_{j}$ with corresponding eigenvalues $λ_{j}$ ( $j = 1, \dots, n$ ), then we have

A V = V Λ

(2)

where

V = [\begin{matrix} v_{1} & v_{2} & \dots & v_{n} \end{matrix}] \in C^{n \times n}

and

Λ = diag (λ_{1}, λ_{2}, \dots, λ_{n}) \in C^{n \times n}

. Postmultiplying

V^{- 1}

to the preceding equation, we have

A = V Λ V^{- 1}

(3)

This is called the eigenvalue decomposition. For the eigenvalue decomposition to hold,

A

needs to have a full set of

n

linearly independent eigenvectors.‡

For linear dynamical systems, we often encounter systems for some state variable $x (t) \in C^{n}$ described by

\dot{x} (t) = A x (t)

(4)

with the solution of

x (t) = \exp (A t) x (0) = V \exp (Λ t) V^{- 1} x (0)

(5)

where

x (0)

denotes the initial condition. Here, the eigenvalues characterize the long-term behavior of linear dynamical systems [6,34] for

x (t)

, as illustrated in Fig. 3. The real and imaginary parts of

λ_{j}

represent the growth (decay) rate and the frequency at which the state variable evolves in the direction of the eigenvector

v_{j}

. For a linear system to be stable, all eigenvalues need to be on the left-hand side of the complex plane, i.e.,

R e (λ_{j}) \leq 0

for all

j

Fig. 3
Dynamic response of a linear system characterized by the eigenvalues (stable: $R e (λ) < 0$ , and unstable: $R e (λ) > 0$ ). Location of example eigenvalues $λ$ are shown by the symbols with corresponding sample solutions $\exp (λ t)$ in inset plots.

For intermediate dynamics, the pseudospectra [33,35,36] can provide insights. The concept of pseudospectra is associated with nonnormality of operators and the sensitivity of the eigenvalues to perturbations. We briefly discuss the pseudospectra in Sec. II.E.

For some problems, there can be a mass matrix $B \in C^{n \times n}$ that appears on the left-hand side of Eq. (4):

B \dot{x} = A x

(6)

In such a case, we are led to a generalized eigenvalue problem of the form

A v = λ B v

(7)

B

is invertible, we can rewrite the preceding equation as

B^{- 1} A v = λ v

(8)

and treat the generalized eigenvalue problem as a standard eigenvalue problem. However, it may not be desirable to consider this reformulation if

B

is not invertible§ or if the inversion of

B

results in ill conditioning (worsening of scaling) of the problem. Note that generalized eigenvalue problems can also be solved with many numerical libraries, which are similar to the standard eigenvalue problem [Eq. (1)]. See the work of Trefethen and Embree [33] and Golub and Loan [31] for additional details on the generalized eigenvalue problems.

B. Singular Value Decomposition

The singular value decomposition (SVD) is one of the most important matrix factorizations, generalizing the eigendecomposition to rectangular matrices. The SVD has many uses and interpretations, especially for dimensionality reduction, where it is possible to use the SVD to obtain optimal low-rank matrix approximations [37]. The singular value decomposition also reveals how a rectangular matrix or operator stretches and rotates a vector. As an illustrative example, consider a set of vectors $v_{j} \in R^{n}$ of unit length that describe a sphere. We can premultiply these unit vectors $v_{j}$ with a rectangular matrix $A \in R^{m \times n}$ as shown in Fig. 4. The semiaxes of the resulting ellipse (ellipsoid) are represented by the unit vectors $u_{j}$ and magnitudes $σ_{j}$ . Hence, we can view the singular values to capture the amount of stretching imposed by matrix $A$ in the directions of the axes of the ellipse.

Fig. 4
Graphical representation of singular value decomposition transforming a unit radius sphere, described by right singular vectors $v_{j}$ , to an ellipse (ellipsoid) with semiaxes characterized by the left singular vectors $u_{j}$ and magnitude captured by the singular values $σ_{j}$ . In this graphical example, we take $A \in R^{3 \times 3}$ .

Generalizing this concept for complex $A \in C^{m \times n}$ , $v_{j} \in C^{n}$ , and $u_{j} \in C^{m}$ , we have

A v_{j} = σ_{j} u_{j}

(9)

In matrix form, the aforementioned relationship can be expressed as

A V = U Σ

(10)

where

U = [\begin{matrix} u_{1} & u_{2} & \dots & u_{m} \end{matrix}] \in C^{m \times m}

and

V = [\begin{matrix} v_{1} & v_{2} & \dots & v_{n} \end{matrix}] \in C^{n \times n}

are unitary matrices¶ and

Σ \in R^{m \times n}

is a diagonal matrix with

σ_{1} \geq σ_{2} \geq \dots \geq σ_{p} \geq 0

along its diagonal, where

p = \min (m, n)

. Now, multiplying

V^{- 1} = V^{*}

from the right side of the preceding equation, we arrive at

A = U Σ V^{*}

(11)

which is referred to as the singular value decomposition. In the preceding equation,

*

denotes the conjugate transpose. The column vectors

u_{j}

and

v_{j}

U

and

V

are called the left and right singular vectors, respectively. Both of the singular vectors can be determined up to a complex scalar of magnitude one (i.e.,

e^{i θ}

, where

θ \in [0, 2 π]

Given a rectangular matrix $A$ , we can decompose the matrix with the SVD in the following graphical manner:

((12)) where we have taken

m > n

in this example. Sometimes, the components in

U

enclosed by the broken lines are omitted from the decomposition, as they are multiplied by zeros in

Σ

. The decomposition that disregards the submatrices in the broken-line boxes are called the reduced SVD (economy-sized SVD), as opposed to the full SVD.

In a manner similar to the eigenvalue decomposition, we can interpret the SVD as a means to represent the effect of matrix operation merely through the multiplication by scalars (singular values) given the appropriate directions. Because the SVD is applied to a rectangular matrix, we need two sets of basis vectors to span the domain and range of the matrix. Hence, we have the right singular vectors $V$ that span the domain of $A$ and the left singular vectors $U$ that span the range of $A$ , as illustrated in Fig. 4. This is different from the eigenvalue decomposition of a square matrix: in which case, the domain and the range are (generally) the same. Although the eigenvalue decomposition requires the square matrix to be diagonalizable, the SVD (on the other hand) can be performed on any rectangular matrix.

C. Relationship Between Eigenvalue and Singular Value Decompositions

The eigenvalue and singular value decompositions are closely related. In fact, the left and right singular vectors of $A \in C^{m \times n}$ are also the orthonormal eigenvectors of $A A^{*}$ and $A^{*} A$ , respectively. Furthermore, the nonzero singular values of $A$ are the square roots of the nonzero eigenvalues of $A A^{*}$ and $A^{*} A$ . Therefore, instead of the SVD, the eigenvalue decomposition can be performed on $A A^{*}$ or $A^{*} A$ to solve for the singular vectors and singular values of $A$ . For these reasons, the smaller of the square matrices of $A A^{*}$ and $A^{*} A$ are often chosen to perform the decomposition in a computationally inexpensive manner as compared to the full SVD. This property is taken advantage of in some of the decomposition methods discussed in the following because flowfield data usually yield a rectangular data matrix that can be very high-dimensional in one direction (e.g., the snapshot POD method [28] in Sec. III).

D. Numerical Libraries for Eigenvalue and Singular Value Decompositions

Eigenvalue and singular value decompositions can be performed with codes that are readily available. We list a few standard numerical libraries to execute eigenvalue and singular value decompositions.

MATLAB: In MATLAB®, the command eig finds the eigenvalues and eigenvectors for standard eigenvalue problems as well as generalized eigenvalue problems. The command svd outputs the singular values and the left and right singular vectors. It can also perform the economy-sized SVD. For small- to moderate-sized problems, MATLAB can offer a user-friendly environment to perform modal decompositions. We provide in Table 2 some common examples of eig and svd in use for canonical decompositions.**

LAPACK: LAPACK (linear algebra package) offers standard numerical library routines for a variety of basic linear algebra problems, including eigenvalue and singular value decompositions. The routines are written in Fortran 90. See the users’ guide [38].††

ScaLAPACK: ScaLAPACK (scalable LAPACK) comprises high-performance linear algebra routines for parallel distributed memory machines. ScaLAPACK solves dense and banded eigenvalue and singular value problems. See the users’ guide [39].‡‡

ARPACK: ARPACK (Arnoldi package) is a numerical library, written in FORTRAN 77, that is specialized to handle large-scale eigenvalue problems as well as generalized eigenvalue problems. It can also perform singular value decompositions. The library is available for both serial and parallel computations. See the users’ guide [40].§§

E. Pseudospectra

Before we transition our discussion to the coverage of modal analysis techniques, let us consider the pseudospectral analysis [33,35], which reveals the sensitivity of the eigenvalue spectra with respect to perturbations to the operator. This is also an important concept in studying transient and input–output dynamics, complementing the stability analysis based on eigenvalues. Concepts from the pseudospectral analysis appear later in the resolvent analysis (Sec. VIII).

For a linear system described by Eq. (4) to exhibit stable dynamics, we require all eigenvalues of its operator $A$ to satisfy $R e (λ_{j} (A)) < 0$ , as illustrated in Fig. 3. Although this criterion guarantees the solution $x (t)$ to be stable for large $t$ , it does not provide insights into the transient behavior of $x (t)$ . To illustrate this point, let us consider an example of $A = V Λ V^{- 1}$ with stable eigenvalues of

λ_{1} = - 0.1, λ_{2} = - 0.2

(13)

and eigenvectors of

v_{1} = {[\cos (\frac{π}{4} - δ), \sin (\frac{π}{4} - δ)]}^{T}, v_{2} = {[\cos (\frac{π}{4} + δ), \sin (\frac{π}{4} + δ)]}^{T}

(14)

where

δ

is a free parameter to choose. Observe that, as

δ

becomes small, the eigenvectors become nearly linearly dependent, which makes the matrix

A

ill conditioned.

Providing an initial condition of $x (t_{0}) = {[1, 0.1]}^{T}$ , we can solve Eq. (5) for different values of $δ$ , as shown in Figs. 5a and 5b. Although all solutions decay to zero due to the stable eigenvalues, the transient growths of $x_{1} (t)$ and $x_{2} (t)$ become noticeable as $δ \to 0$ . The large transient for small $δ$ is caused by the eigenvectors becoming nearly parallel, which necessitate large coefficients to represent the solution [i.e., $x (t) = α_{1} (t) v_{1} + α_{2} (t) v_{2}$ , where $| α_{1} |$ and $| α_{2} | ≫ 1$ during the transient]. As such, the solution grows significantly during the transient before the decay from the negative eigenvalues starts to take over the solution behavior as $t$ becomes large. Thus, we observe that the transient behavior of the solution is not controlled by the eigenvalues of $A$ . Nonnormal operators (i.e., operators for which $A A^{*} \neq A^{*} A$ ) have nonorthogonal eigenvectors and can exhibit this type of transient behavior. Thus, it is important that care is taken when we examine transient dynamics caused by nonnormal operators. In fluid mechanics, the dynamics of shear-dominant flows often exhibit nonnormality.

Fig. 5
Representations of a) example of transient growth caused by increasing level of nonnormality from decreasing $δ$ , b) trajectories of $x_{1} (t)$ vs $x_{2} (t)$ exhibiting transient growth, and c–e) pseudospectra expanding for different values of $δ$ . $ε$ pseudospectra are shown with values of $\log_{10} (ε)$ placed on the contours, and stable eigenvalues are depicted with $\times$ .

To further assess the influence of $A$ on the transient dynamics, let us examine here how the eigenvalues are influenced by perturbations on $A$ . That is, we consider

Λ_{ε} (A) = {z \in C : z \in Λ (A + Δ A), w h e r e ‖ Δ A ‖ \leq ε}

(15)

This subset of perturbed eigenvalues is known as the

ε

pseudospectrum of

A

. It is also commonly known with the following equivalent definition:

Λ_{ε} (A) = {z \in C : {‖ z I - A ‖}^{- 1} \geq ε^{- 1}}

(16)

Note that, as

ε \to 0

, we recover the eigenvalues (0 pseudospectrum); and as

ε \to \infty

, the subset

Λ_{\infty} (A)

occupies the entire complex domain. To numerically determine the pseudospectra, we can use the following definition based on the minimum singular value of

(z I - A)

Λ_{ε} (A) = {z \in C : σ_{\min} (z I - A) \leq ε}

(17)

which is equivalent to

Λ_{ε} (A)

, described by Eqs. (15) and (16). If

A

is normal, the pseudospectrum

Λ_{ε} (A)

is the set of points away from

Λ_{0} (A)

(eigenvalues) by only less than or equal to

ε

on the complex plane. However, as

A

becomes nonnormal, the distance between

Λ_{0} (A)

and

Λ_{ε} (A)

may become much larger. As will be discussed later, the resolvent analysis in Sec. VIII considers the pseudospectra along the imaginary axis [6] (i.e.,

z \to ω i

, where

ω \in R

Let us return to the example given by Eqs. (13) and (14) and compute the pseudospectra for decreasing $δ$ of 0.01, 0.001, and 0.0001, as shown in Figs. 5c–5e, respectively. Here, the contours of the $ε$ pseudospectra are drawn for the same values of $ε$ . With decreasing $δ$ , the matrix $A$ becomes increasingly nonnormal and susceptible to perturbations. The influence of nonnormality on the spectra is clearly visible with the expanding $ε$ pseudospectra. It should be noticed that some of the pseudospectra contours penetrate into the right-hand side of the complex plane, suggesting that perturbations of such magnitude may thrust the system to become unstable even with stable eigenvalues. This nonnormal feature can play a role in destabilizing the dynamics with perturbations or nonlinearity.

The transient dynamics of $\dot{x} = A x$ can be related to how the $ε$ pseudospectrum of $A$ expands from the eigenvalues as parameter $ε$ is varied. The pseudospectra of $A$ can provide a lower bound on the amount of transient amplification by $\exp (A t)$ . If $Λ_{ε} (A)$ extends a distance $η$ into the right half-plane for a given $ε$ , it can be shown through Laplace transform that $‖ \exp (A t) ‖$ must be as large as $η / ε$ for some $t > 0$ . If we let a constant $κ$ for $A$ be defined as the supremum of this ratio over all $ε$ , the lower bound for the solution can then be shown to take the form of [41]

\sup_{t \geq 0} ‖ \exp (A t) ‖ \geq κ

(18)

This constant

κ

is referred to as the Kreiss constant, which provides an estimate of how the solution [Eq. (5)] behaves during the transient. This estimate is not obtained from the eigenanalysis but from the pseudospectral analysis. The same concept applies to time-discretized linear dynamics [42]. Readers can find applications of the pseudospectral analysis to fluid mechanics in the works of Trefethen et al. [35], Trefethen and Embree [33], and Schmid [36].

III. Proper Orthogonal Decomposition

The proper orthogonal decomposition is a modal decomposition technique that extracts modes based on optimizing the mean square of the field variable being examined. It was introduced to the fluid dynamics/turbulence community by Lumley [27] as a mathematical technique to extract coherent structures from turbulent flowfields. The POD technique, also known as the Karhunen–Loève procedure [43,44], provides an objective algorithm to decompose a set of data into a minimal number of basis functions or modes to capture as much energy as possible. The method itself is known under a variety of names in different fields: POD, principal component analysis (PCA), Hotelling analysis, empirical component analysis, quasi-harmonic modes, empirical eigenfunction decomposition, and others. Closely related to this technique is factor analysis, which is used in psychology and economics. Roots of the POD can be traced back to the middle of the 19th century to the matrix diagonalization technique, which is ultimately related to the SVD (Sec. II). Excellent reviews on the POD can be found in [1,45], and chapter 3 of [46].

In applications of POD to a fluid flow, we start with a vector field $q (ξ, t)$ (e.g., velocity) with its temporal mean $\bar{q} (ξ)$ subtracted and assume that the unsteady component of the vector field can be decomposed in the following manner:

q (ξ, t) - \bar{q} (ξ) = \sum_{j} a_{j} ϕ_{j} (ξ, t)

(19)

where

ϕ_{j} (ξ, t)

and

a_{j}

represent the modes and expansion coefficients, respectively. Here,

ξ

denotes the spatial vector.¶¶ This expression represents the flowfield in terms of a generalized Fourier series for some set of basis functions

ϕ_{j} (ξ, t)

. In the framework of the POD, we seek the optimal set of basis functions for given flowfield data. In early applications of the POD, this typically led to modes that were functions of space and time/frequency [47

–

51], as will also be discussed in the following.

Modern applications of modal decompositions have further sought to split space and time, hence only needing spatial modes. In that context, the preceding equation can be written as

q (ξ, t) - \bar{q} (ξ) = \sum_{j} a_{j} (t) ϕ_{j} (ξ)

(20)

where the expansion coefficients

a_{j}

are now time dependent. Note that Eq. (20) explicitly employs a separation of variables, which may not be appropriate for all problems. The application of the two forms listed previously should depend on the properties of the flow and the information one wishes to extract, as discussed by Holmes et al. [52]. In what follows, we will discuss the properties of the POD by assuming that the desire is to extract a spatially dependent set of modes.

The POD is one of the most widely used techniques in analyzing fluid flows. There are a large number of variations of the POD technique, with applications including fundamental analysis of fluids flows, reduced-order modeling, data compression/reconstruction, flow control, and aerodynamic design optimization. Because the POD serves as the basis and motivation for the development of other modal decomposition techniques, we provide a somewhat detailed overview of the POD in the following.

A. Description

1. Algorithm

The inputs are snapshots of any scalar (e.g., pressure, temperature) or vector (e.g., velocity, vorticity) field $q (ξ, t)$ over one-, two-, or three-dimensional discrete spatial points $ξ$ at discrete times $t_{i}$ .

The outputs are a set of orthogonal modes $ϕ_{j} (ξ)$ with their corresponding temporal coefficients $a_{j} (t)$ and energy levels $λ_{j}$ arranged in the order of their relative amount of energy. The fluctuations in the original field are expressed as a linear combination of the modes and their corresponding temporal coefficients:

q (ξ, t) - \bar{q} (ξ) = \sum_{j} a_{j} (t) ϕ_{j} (ξ)

(21)

We discuss three main approaches to perform the POD of the flowfield data: namely, the spatial (classical) POD method, the snapshot POD method, and the SVD. In the following, we briefly describe these three methods and discuss how they are related to each other.

Spatial (Classical) POD Method. With the POD, we determine the set of basis functions that optimally represents the given flowfield data. First, given the flowfield $q (ξ, t)$ , we prepare snapshots of the flowfield stacked in terms of a collection of column vectors $x (t)$ . That is, we consider a collection of finite-dimensional data vectors that represents the flowfield:

x (t) = q (ξ, t) - \bar{q} (ξ) \in R^{n}, t = t_{1}, t_{2}, \dots, t_{m}

(22)

Here,

x (t)

is taken to be the fluctuating component of the data vector with its time-averaged value

\bar{q} (ξ)

removed. Although the data vector can be written as

x (ξ, t)

, we simply write

x (t)

to emphasize that it is being considered as a snapshot at time

t

. An example of forming the data vector

x (t)

for a given flowfield is provided in the Appendix.

The objective of the POD analysis is to find the optimal basis vectors that can best represent the given data. In other words, we seek the vectors $ϕ_{j} (ξ)$ in Eq. (20) that can represent $q (ξ)$ in an optimal manner and with the least number of modes. The solution to this problem [37] can be determined by finding the eigenvectors $ϕ_{j}$ and the eigenvalues $λ_{j}$ from

R ϕ_{j} = λ_{j} ϕ_{j}, ϕ_{j} \in R^{n}, λ_{1} \geq \dots \geq λ_{n} \geq 0

(23)

where

R

is the covariance matrix of vector

x (t)

.***

R = \sum_{i = 1}^{m} x (t_{i}) x^{T} (t_{i}) = X X^{T} \in R^{n \times n}

(24)

where the matrix

X

represents the

m

snapshot data being stacked into a matrix form of

X = [\begin{matrix} x (t_{1}) & x (t_{2}) & \dots & x (t_{m}) \end{matrix}] \in R^{n \times m}

(25)

The size of the covariance matrix

n

is based on the spatial degrees of freedom of the data. For fluid flow data,

n

is generally large and is equal to the number of grid points times the number of variables to be considered in the data, as illustrated in Eq. (A2) of the Appendix. See the Appendix for an example of preparing the data matrix from the velocity field data.

The eigenvectors found from Eq. (37) are called the POD modes. It should be noted that the POD modes are orthonormal. That means that the inner product††† between the modes satisfies

⟨ ϕ_{j}, ϕ_{k} ⟩ \equiv \int_{V} ϕ_{j} \cdot ϕ_{k} d V = δ_{j k}, j, k = 1, \dots, n

(26)

Consequently, the eigenvalues

λ_{k}

convey how well each eigenvector

ϕ_{k}

captures the original data in the

L_{2}

sense (scaled by

m

). When the velocity vector is used for

x (t)

, the eigenvalues correspond to the kinetic energy captured by the respective POD modes. If the eigenvalues are arranged from the largest to the smallest, in decreasing order, the POD modes are arranged in the order of importance in terms of capturing the kinetic energy of the flowfield.

We can use the eigenvalues to determine the number of modes needed to represent the fluctuations in the flowfield data. Generally, we retain only $r$ number of modes to express the flow such that

\sum_{j = 1}^{r} λ_{j} / \sum_{j = 1}^{n} λ_{j} \approx 1

(27)

With the determination of the important POD modes, we can represent the flowfield only in terms of the finite or truncated series

q (ξ, t) - \bar{q} (ξ) \approx \sum_{j = 1}^{r} a_{j} (t) ϕ_{j} (ξ)

(28)

in an optimal manner, effectively reducing the high-dimensional (

n

) flowfield to be represented only with

r

modes. The temporal coefficients are determined accordingly by

a_{j} (t) = ⟨ q (ξ, t) - \bar{q} (ξ), ϕ_{j} (ξ) ⟩ = ⟨ x (t), ϕ_{j} ⟩

(29)

Method of Snapshots. When the spatial size of the data $n$ is very large, the size of the correlation matrix $R = X X^{T}$ becomes very large ( $n \times n$ ), making the use of the classical spatial POD method for finding the eigenfunctions practically impossible. Sirovich [28] pointed out that the temporal correlation matrix will yield the same dominant spatial modes while giving rise to a much smaller and computationally more tractable eigenvalue problem. This alternative approach, called the method of snapshots, takes a collection of snapshots $x (t_{i})$ at discrete time levels $t_{i}$ , $i = 1, 2, \dots, m$ , with $m ≪ n$ ; and it solves an eigenvalue problem of a smaller size ( $m \times m$ ) to find the POD modes. The number of snapshots $m$ should be chosen such that important fluctuations in the flowfield are well resolved in time.

The method of snapshots relies on solving an eigenvalue problem of a much smaller size:

X^{T} X ψ_{j} = λ_{j} ψ_{j}, ψ_{j} \in R^{m}, m ≪ n

(30)

where

X^{T} X

is of size

m \times m

instead of the original eigenvalue problem of size

n \times n

[Eq. (37)]. Although we are analyzing the smaller eigenvalue problem, the same nonzero eigenvalues are shared by

X^{T} X

and

X X^{T}

and the eigenvectors of these matrices can be related to each other (see Sec. II.C). With the eigenvectors

ψ_{j}

of the aforementioned smaller eigenvalue problem determined, we can recover the POD modes through

ϕ_{j} = X ψ_{j} \frac{1}{\sqrt{λ_{j}}} \in R^{n}, j = 1, 2, \dots, m

(31)

which can equivalently be written in matrix form as

ϕ = X Ψ Λ^{- 1 / 2}

(32)

where

ϕ = [\begin{matrix} ϕ_{1} & ϕ_{2} & \dots & ϕ_{m} \end{matrix}] \in R^{n \times m}

and

Ψ = [\begin{matrix} ψ_{1} & ψ_{2} & \dots & ψ_{m} \end{matrix}] \in R^{m \times m}

Due to the significant reduction in the required computation and memory resources, the method of snapshots has been widely used to determine POD modes from high-dimensional fluid flow data. In fact, this snapshot-based approach is presently the most widely used POD method in fluid mechanics.

SVD and POD. Let us consider the relation between the POD and SVD, as discussed in Sec. II.C. Recall that the SVD [3,29,31] can be applied to a rectangular matrix to find the left and right singular vectors. In matrix form, a data matrix $X$ can be decomposed directly with the SVD as

X = ϕ Σ Ψ^{T}

(33)

where

ϕ \in R^{n \times n}

Ψ \in R^{m \times m}

, and

Σ \in R^{n \times m}

, with

m < n

. The matrices

ϕ

and

Ψ

contain the left and right singular vectors‡‡‡ of

X

, and matrix

Σ

holds the singular values (

σ_{1}, σ_{2}, \dots, σ_{m}

) along its diagonal. These singular vectors held in

ϕ

and

Ψ

are identical to the eigenvectors of

X X^{T}

and

X^{T} X

, respectively. Moreover, the singular values and the eigenvalues are related by

σ_{j}^{2} = λ_{j}

. This means that the SVD can be directly performed on

X

to determine the POD modes

ϕ

. Note that the SVD finds

ϕ

from Eq. (33) with a full size of

n \times n

, as compared to

ϕ

from the snapshot approach with a size of

n \times m

, holding only the leading

m

modes.

The terms POD and SVD are often used interchangeably in the literature. However, the SVD is a decomposition technique for rectangular matrices and POD can be seen as a decomposition formalism for which the SVD can be one of the approaches to determine its solution. Although the method of snapshot is preferred for handling large datasets, the SVD-based technique to determine the POD modes is known to be robust against roundoff errors [3].

2. Notes

Optimality. The POD modes are computed in the optimal manner in the $L_{2}$ sense [1]. If the velocity or vorticity field is used to determine the POD modes, the modes are optimal to capture the kinetic energy or enstrophy, respectively, of the flowfield. Moreover, POD decomposition is optimal, not only in terms of minimizing the mean-square error between the signal and its truncated representation but also minimizing the number of modes required to describe the signal for a given error [53].

The optimality (the fastest-convergent property) of the POD reduces the amount of information required to represent statistically dependent data to a minimum. This crucial feature explains the wide usage of the POD in a process of analyzing data. For this reason, the POD is used extensively in the fields of detection, estimation, pattern recognition, and image processing.

Reduced-Order Modeling. The orthogonality of POD modes [i.e., $⟨ ϕ_{j}, ϕ_{k} ⟩ = δ_{j k}$ ] is an attractive property for constructing reduced-order models [1,54]. Galerkin projection can be used to reduce high-dimensional discretizations of partial differential equations into reduced-order ordinary differential equation models for the temporal coefficients $a_{j} (t)$ . POD modes have been used to construct Galerkin projection-based reduced-order models for incompressible [48,55–57] and compressible [58] flows.

Traveling Structures. With real-valued POD modes, traveling structures cannot be represented as a single mode. In general, traveling structures are represented by a pair of stationary POD modes, which are similar but appear shifted in the advection direction. See, for instance, modes 1 and 2 in Figs. 1 and 6. One way to understand the emergence of POD mode pairs for traveling structures is to consider the following traveling sine wave example:

\sin (ξ - c t) = \cos (c t) \sin (ξ) - \sin (c t) \cos (ξ)

(34)

Here, we have a pair of spatial modes,

\sin (ξ)

and

\cos (ξ)

, which are shifted by a phase of

π / 2

. Note that we can combine the pair of modes into a single mode if a complex representation of POD modes is considered [59]. There are variants of the POD analysis specialized for traveling structures [60,61].

Constraints. With linear superposition of the POD modes in representing the flowfield, each and every POD mode also satisfies linear constraints, such as the incompressibility constraint and the no-slip boundary condition. This statement assumes that the given data also satisfy these constraints.

Homogeneous Directions. For homogeneous, periodic, or stationary (translationally invariant) directions, POD modes reduce to Fourier modes [62].

Spectral POD. It is also possible to consider the use of the POD in the frequency domain. Spectral POD provides time-harmonic modes at discrete frequencies from a set of realizations of the temporal Fourier transform of the flowfield. This application of the POD provides an orthogonal basis of modes at discrete frequencies that are optimally ranked in terms of energy because the POD reduces to a harmonic analysis over directions that are stationary or periodic [62,63]. Hence, the method is effective at extracting coherent structures from statistically stationary flows and has been successfully applied in the early applications of the POD, as well as other flows including turbulent jets [64 –68]. This frequency-based approach to the POD overcomes some of the weaknesses, and the discussions in the strengths and weaknesses do not directly apply here [69].

Spectral POD can be estimated from a time series of snapshots in the form of Eq. (25) using Welch’s method [70]. First, the data are segmented into a number of $n_{b}$ (of potentially overlapping) blocks or realizations, consisting of $m_{FFT}$ snapshots each. Under the ergodicity hypothesis, each block can be regarded as statistically independent realization of the flow. We proceed by calculating the temporal Fourier transform

{\hat{X}}^{(l)} = [\begin{matrix} \hat{x} {(ω_{1})}^{(l)} & \hat{x} {(ω_{2})}^{(l)} & \dots & \hat{x} {(ω_{m_{FFT}})}^{(l)} \end{matrix}] \in R^{n \times m_{FFT}}

(35)

of each block, where superscript

l

denotes the

l

th block. All realizations of the Fourier transform at a specific frequency

ω_{k}

are now collected into a new data matrix:

{\hat{X}}_{ω_{k}} = [\begin{matrix} \hat{x} {(ω_{k})}^{(1)} & \hat{x} {(ω_{k})}^{(2)} & \dots & \hat{x} {(ω_{k})}^{(n_{b})} \end{matrix}]

(36)

The product

{\hat{X}}_{ω_{k}} {\hat{X}}_{ω_{k}}^{T}

forms the cross-spectral density matrix, and its eigenvalue decomposition

{\hat{X}}_{ω_{k}} {\hat{X}}_{ω_{k}}^{T} ϕ_{ω_{k, j}} = λ_{ω_{k, j}} ϕ_{ω_{k, j}}, ϕ_{ω_{k, j}} \in R^{n}, λ_{ω_{k}, 1} \geq \dots \geq λ_{ω_{k}, n_{b}} \geq 0

(37)

yields the spectral POD modes

ϕ_{ω_{k, j}}

and corresponding modal energies

λ_{ω_{k, j}}

, respectively. As an example, the most energetic spectral POD mode of a turbulent jet at one frequency is shown in Fig. 16d. The spectral POD as described here dates back to the early work of Lumley [62] and is not related to the method under the same name proposed by Sieber et al. [71], whose approach blended between the POD and discrete Fourier transform modes by filtering the temporal correlation matrix.

3. Strengths and Weaknesses

Strengths:

The POD gives an orthogonal set of basis vectors with the minimal dimension. This property is useful in constructing a reduced-order model of the flowfield.
POD modes are simple to compute using either the (classical) spatial or snapshot methods. The method of snapshots is especially attractive for high-dimensional spatial datasets.
Incoherent noise in the data generally appears as high-order POD modes, provided that the noise level is lower than the signal level. The POD analysis can be used to practically remove the incoherent noise from the dataset by simply removing high-order modes from the expansion.
The POD (PCA) analysis is very widely used in a broad spectrum of studies. It is used for pattern recognition, image processing, and data compression.

Weaknesses:

As the POD is based on second-order correlation, higher-order correlations are ignored.
The temporal coefficients of spatial POD modes generally contain a mix of frequencies. The spectral POD discussed previously addresses this issue.
The POD arranges modes in the order of energy contents and not in the order of the dynamical importance. This point is addressed by balanced POD and DMD analyses.
It is not always clear how many POD modes should be kept, and there are many different truncation criteria.

B. Illustrative Examples

1. Turbulent Separated Flow over an Airfoil

We present an example of applying the POD analysis on the velocity field obtained from three-dimensional large-eddy simulation (LES) of turbulent separated flow over a NACA 0012 airfoil [72]. The flow is incompressible with spanwise periodicity at $R e = 23, 000$ and $α = 9 \deg$ . Visualized in Fig. 6 (left) are the instantaneous and time-averaged streamwise velocities on a spanwise slice. We can observe that there are large-scale vortical structures in the wake from von Kármán shedding, yielding spatial and temporal fluctuations about the mean flow. Also present are the finer-scale turbulent structures in the flow. Performing the POD on the flowfield data, we can find the dominant modes [73,74]. Here, the first four dominant POD modes with the percentage of kinetic energy held by the modes are shown in Fig. 6 (middle and right). The four modes shown together capture approximately 19% of the unsteady fluctuations over the examined domain. Modes 1 and 2 (first pair) represent the most dominant fluctuations in the flowfield, possessing equal levels of kinetic energy, and amounting to oscillatory (periodic) modes. Modes 3 and 4 (second pair) represent the subharmonic spatial structures of modes 1 and 2 in this example. Compared to the laminar case shown in Fig. 1, the number of modes required to reconstruct this turbulent flow is increased due to the emergence of multiple spatial scales and higher-dimensionality of the turbulent flow. The dominant features of the shown POD modes share similarities with the laminar flow example shown in Fig. 1, despite the large difference in the Reynolds numbers.

Fig. 6
POD analysis of turbulent flow over a NACA0012 airfoil at $R e = 23,000$ and $α = 9 \deg$ . Shown are the instantaneous and time-averaged streamwise velocity fields and the associated four most dominant POD modes [73,74]. Reprinted with permission from Springer.

2. Compressible Open-Cavity Flows

As the second sample application of the POD, let us look at an analysis of the velocity field from open-cavity flow experiments [75]. In this study, the snapshot POD was applied to two-component particle image velocimetry data acquired with several different freestream Mach numbers from 0.2 through 0.73. Here, the dominant mode (mode 1) contains between 15 and 20% of the energy, with 50% represented by the first seven modes. This study reveals the similarity of the modes among four of the freestream Mach numbers investigated, even though there were differences in the mean flow patterns. The first five POD modes associated with the vertical velocity component are presented in Fig. 7. In this figure, we can observe a representation of the vortical structures in the cavity shear layer with similar wavelengths, regardless of the freestream Mach number. Further quantitative analysis of the similarity of the modes was verified by checking the orthogonality between the modes for the various applications. The similarity among the modes implies that the underlying turbulence had the same structure, at least over the range of freestream Mach numbers investigated, regardless of the mean flow differences.

Fig. 7
POD modes of the vertical fluctuating velocity for flows over a rectangular cavity at Mach numbers of $M_{\infty} = 0.19$ and 0.29 [75]. Reprinted with permission from AIP Publishing.

C. Outlook

The POD has been the bedrock of modal decomposition techniques to extract coherent structures for unsteady fluid flows. To address some of the shortcomings of the standard POD analysis, many variations have emerged: namely, the balanced POD [76] (see Sec. IV), the split POD [77], the sequential POD [78,79], The temporal POD [80], and the joint POD [81], among others. A number of overarching studies have emerged to bridge the gap between the POD and other decomposition methods, revisiting some of the early POD discussions by Lumley [27] and George [63]. Recently, theoretical connections have been made between the spectral POD and several other methods, including resolvent analysis (Sec. VIII) [69,82,83] and other data-based methods, including the spatial POD described earlier and dynamic mode decomposition (Sec. V) [69].

One of the most attractive properties of POD modes is orthogonality. This feature allows us to develop models that are low in order and sparse. Taking advantage of such a property, there have been efforts to construct reduced-order models based on Galerkin projection to capture the essential flow physics [1,48,55–57] and to implement model-based closed-loop flow control [54,84]. More recently, there are emerging approaches that take advantage of the POD modes to model and control fluid flows, leveraging cluster-based analysis [85] and networked-oscillator representation [59] of complex unsteady fluid flows.

IV. Balanced Proper Orthogonal Decomposition

Balanced proper orthogonal decomposition is a modal decomposition technique that can extract two sets of modes for specified inputs and outputs. Here, the inputs are typically external disturbances or actuation used for flow control. The outputs are typically the available sensor measurements or the quantities we want to capture with a model (for instance, they could be amplitudes of the POD modes).

This method is an approximation of a technique called balanced truncation [86], which is a standard method used in control theory that balances the properties of controllability and observability. The most controllable states correspond to those that are most easily excited by the inputs, and the most observable states correspond to those that excite large future outputs. In a reduced-order model, we wish to retain both the most controllable modes and the most observable modes; but, the difficulty is that, for some systems (particularly systems that are nonnormal, which arise in many shear flows), states that have very small controllability might have very large observability, and vice versa. Balancing involves determining a coordinate system in which the most controllable directions in state space are also the most observable directions. We then truncate the states that are the least controllable/observable.

Balanced POD is closely related to the POD: both procedures produce a set of modes that describe the coherent structures in a given fluid flow, and the computations required are similar (both involve the SVD). However, there are some important differences. The POD provides a single set of modes that is orthogonal and ranked by energy content. In contrast, the balanced POD provides two sets of modes (balancing modes and adjoint modes), which form a biorthogonal set and are ranked by controllability/observability (which we can think of as the importance to the input–output dynamics). With both the POD and balanced POD, a quantity $q (ξ, t)$ is expanded as

q (ξ, t) = \sum_{j = 1}^{n} a_{j} (t) ϕ_{j} (ξ)

(38)

where

ϕ_{j}

are the POD modes or direct balancing modes, and

a_{j} (t)

are scalar temporal coefficients. For the POD, the modes

ϕ_{j}

are orthonormal (which means

⟨ ϕ_{j}, ϕ_{k} ⟩ = δ_{j k}

), so the coefficients

a_{j}

are computed by

a_{j} (t) = ⟨ q, ϕ_{j} ⟩

. With balanced POD, biorthogonality of the balancing modes

ϕ_{j}

and adjoint modes

ψ_{k}

means that these satisfy

⟨ ϕ_{j}, ψ_{k} ⟩ = δ_{j k}

, and the coefficients

a_{j}

are then determined by

a_{j} (t) = ⟨ q, ψ_{j} ⟩

The dataset used for the balanced POD is also quite specific: it consists of the linear responses of the system to impulsive inputs (one time series for each input), as well as impulse responses of an adjoint system (one adjoint response for each output). It is these adjoint simulations that enable the balanced POD to determine the observability (or sensitivity) of different states, which makes the procedure so effective for nonnormal systems. However, because adjoint information is required, it is usually not possible to apply the balanced POD to experimental data. It has been shown that a system identification method called the eigensystem realization algorithm (ERA) [87] produces reduced-order models that are equivalent to balanced-POD-based models, without the need for adjoint responses, and can therefore be used on experimental data [88]. For the full details of the balanced POD, see the work of Rowley [76] or the second edition of the book by Holmes et al. ([1] Chap. 5). A description of a related method was also given in the work of Willcox and Peraire [89].

A. Description

1. Algorithm

As the balanced POD analysis is founded on linear state-space systems, the inputs to the balanced POD should be obtained from linear dynamics.

The inputs are two sets of snapshots from a linearized forward simulation and a companion adjoint simulation.

The outputs are sets of balancing modes and adjoint modes ranked in the order of the Hankel singular values. These modes comprise a coordinate transform that balances the controllability and observability of the system.

The balanced POD is based on the concept of balanced truncation that provides a balancing measure between controllability and observability in the transformed coordinate. This approach seeks for the balancing transform $ϕ$ and its inverse transform $Ψ$ that can diagonalize and equate the (empirical) controllability and observability Gramians, $W_{c}$ and $W_{o}$ , respectively, such that $Ψ^{*} W_{c} Ψ = ϕ W_{o} ϕ^{*} = Σ$ is a diagonal matrix. Although the controllability Gramian can be determined from the forward simulation, the observability Gramian requires results from the adjoint simulation.

Now, let us introduce the forward and adjoint linear systems and present the snapshot-based balanced POD technique [76] that is analogous to the snapshot-based POD method. Following standard state-space notation from control theory, the forward and adjoint simulations solve

{\begin{cases} \dot{x} = A x + B u \\ y = C x \end{cases}

(39)

and

{\begin{cases} \dot{z} = A^{T} z + C^{T} v \\ w = B^{T} z \end{cases}

(40)

For the forward dynamics [Eq. (39)],

x (t) \in R^{n}

u (t) \in R^{p}

, and

y (t) \in R^{q}

are the state, input, and output vectors; and

A \in R^{n \times n}

B \in R^{n \times p}

, and

C \in R^{q \times n}

are the state, input, and output matrices, respectively. For the adjoint system [Eq. (40)],

z (t) \in R^{n}

v (t) \in R^{q}

, and

w (t) \in R^{p}

are the adjoint state, input, and output vectors.

Based on the solution of linear and adjoint simulations, we can construct the data matrices $X$ and $Z$ (also see the Appendix). For simplicity, let us consider in the following a single-input single-output system with $p = q = 1$ , for which we can construct the data matrices as

X = [\begin{matrix} x (t_{1}) & x (t_{2}) & \dots & x (t_{m}) \end{matrix}] \in R^{n \times m}

(41)

and

Z = [\begin{matrix} z (t_{1}) & z (t_{2}) & \dots & z (t_{m}) \end{matrix}] \in R^{n \times m}

(42)

For a multi-input multioutput system, the solutions from the linear and adjoint simulations can be stacked in an analogous manner to construct

X

and

Z

, as discussed by Rowley [76]. Consequently, the empirical controllability and observability Gramians are given by

W_{c} \approx X X^{T}

and

W_{o} \approx Z Z^{T}

, respectively. The balancing transform

ϕ

and its inverse

Ψ

can simply be found as

ϕ = X V Σ^{- 1 / 2} and Ψ = Z U Σ^{- 1 / 2}

(43)

where

U

Σ

, and

V

are determined from the SVD of the matrix product

Z^{T} X

Z^{T} X = U Σ V^{T}

(44)

The columns

ϕ_{j}

and

ψ_{j}

ϕ

and

Ψ

correspond to the balancing and adjoint modes ranked in the order of

Σ = diag (σ_{1}, σ_{2}, \dots, σ_{m})

, which are referred to as the Hankel singular values.

2. Notes

Biorthogonality. The property of biorthogonality ( $⟨ ϕ_{j}, ψ_{k} ⟩ = δ_{j k}$ ) can be used as a projection to derive a reduced-order model known as the Petrov–Galerkin model [90].

Eigenvalue Realization Algorithm. If we are only interested in deriving a reduced-order model based on balanced truncation without the need to access the balancing and adjoint modes, we can use the ERA and remove the requirement to perform the adjoint simulation. For additional details on the derivation of the models and the ERA, the readers can refer to the work of Juang and Pappa [87] and Ma et al. [88].

3. Strengths and Weaknesses

Strengths.

The balanced POD is particularly attractive for capturing the dynamics of nonnormal systems with large transient growth [90]. Because the POD ranks the modes based on energy content, nonnormal characteristics of the flow may not be captured. In contrast, with the balanced POD models capture small energy perturbations that are highly observable (typically through the adjoint modes). For this reason, balanced POD models typically perform better than POD models for nonnormal systems.
The balanced POD provides an input–output model suitable for feedback control [90].

Weaknesses.

Snapshots from adjoint simulations are needed, which makes the balanced POD analysis difficult or impossible to perform with experimental measurements. However, the ERA can be used instead if only the balanced-POD-based model is sought for without access to the balancing and adjoint modes [87,88].
Both forward and adjoint simulations should be based on linear dynamics (although various extensions to nonlinear systems have been introduced [91,92]).

B. Illustrative Example

Control of Wake Behind a Flat-Plate Wing. The balanced POD has been applied to analyze and control the unsteady wake behind a flat-plate wing at $R e = 100$ . In the work of Ahuja and Rowley [93], they performed linearized forward and adjoint simulations of flow over a flat-plate wing to determine balancing and adjoint modes, as shown in Fig. 8 (left). The balancing modes resemble those of the traditional POD modes, but the adjoint modes highlight regions of the flow that can trigger large perturbations downstream. These modes are used to develop models and closed-loop controllers to stabilize the naturally unstable fluid flows, as depicted in Fig. 8 (right). Although their work necessitated snapshots from the forward and adjoint simulations, the requirement for adjoint simulations was later removed by the use of the ERA in the work by Ma et al. [88]. Although the balancing and adjoint modes were not revealed from the ERA, the resulting reduced-order model was shown to be identical to the balanced-POD-based model. This was numerically demonstrated using the same flat-plate wing problem along with a successful implementation of observer-based feedback control.

Fig. 8
Use of balanced POD analysis for feedback stabilization of the unstable wake behind a flat plate ( $R e = 100$ , $α = 35 \deg$ ) [93]. First and third balancing and adjoint modes are shown on the left. Baseline lift history is shown with a dashed line, and controlled cases with different initiation times of feedback control using the balanced-POD-based reduced-order model are plotted in solid lines on the right. Reprinted with permission from Cambridge University Press.

C. Outlook

Although the balanced POD has been successfully applied to model and control a number of fluid systems [90,93 –96], further research could validate and extend its use to a wider range of fluids applications. Work toward this goal currently includes application to the control of nonlinear systems [97], direct application to unstable systems [98], and use with harmonically forced data [99]. Further work could also seek algorithmic variants that improve efficiency. Examples of work in this direction include analytic treatment of impulse response tails [100] and the use of randomized methods [101,102]. Although, at present, the need for adjoint simulations makes the balanced POD unsuitable for experimental data, future work could possibly remove this restriction by making use of the connections that the balanced POD shares with methods such as the ERA and DMD.

V. Dynamic Mode Decomposition

Dynamic mode decomposition provides a means to decompose time-resolved data into modes, with each mode having a single characteristic frequency of oscillation and growth/decay rate. DMD is based on the eigendecomposition of a best-fit linear operator that approximates the dynamics present in the data. This technique was first introduced to the fluids community in an American Physical Society talk [103], which was subsequently followed by an archival paper by Schmid [104]. The connections with the Koopman operator (see Sec. VI) were given by Rowley et al. [105], which explained the meaning of DMD for a nonlinear system (see also the review articles by Mezić [106] and Tu et al. [107]). There have been a number of different formulations and interpretations of DMD since then [4], which are mentioned in this section.

In many ways, DMD may be viewed as combining favorable aspects of both the POD and the discrete Fourier transform [105,108], resulting in spatiotemporal coherent structures identified purely from data. Because DMD is rooted firmly in linear algebra, the method is highly extensible, spurring considerable algorithmic developments. Moreover, as DMD is purely a data-driven algorithm without the requirement for governing equations, it has been widely applied beyond fluid dynamics: in finance [109], video processing [110–112], epidemiology [113], robotics [114], and neuroscience [115]. As with many modal decomposition techniques, DMD is most often applied as a diagnostic to provide physical insight into a system. The use of DMD for future-state prediction, estimation, and control is generally more challenging and less common in the literature [4].

A. Description

1. Algorithm

The inputs are a set of snapshot pairs from fluids experiments or simulations, where the two snapshots in each pair are separated by a constant interval of time. Often, this will come from a time series of data.

The outputs are DMD eigenvalues and modes. The modes are spatial structures that oscillate and/or grow/decay at rates given by the corresponding eigenvalues. These come from the eigendecomposition of a best-fit linear operator that approximates the dynamics present in the data.

We begin by collecting snapshots of data and arranging them as columns of matrices $X$ and $X^{#}$ , such that

X = [\begin{matrix} x (t_{1}) & x (t_{2}) & \dots & x (t_{m}) \end{matrix}] \in R^{n \times m} and X^{#} = [\begin{matrix} x (t_{2}) & x (t_{3}) & \dots & x (t_{m + 1}) \end{matrix}] \in R^{n \times m}

(45)

In DMD, we approximate the relationship between the snapshots in a linear manner such that

X^{#} = A X

(46)

The matrix

A

may be defined by

A = X^{#} X^{+}

, where

X^{+}

denotes the pseudoinverse of

X

. The DMD eigenvalues and modes are then defined as the eigenvalues and eigenvectors of

A

[107]. It is common that the number of snapshots is smaller than the number of components of each snapshot (

m ≪ n

). In this case, it is not efficient to compute

A

explicitly, so we normally use an algorithm such as the following, which is similar to algorithm 2 in the work of Tu et al. [107]:

Perform the reduced SVD (Sec. II.B) of $X$ , letting $X = U Σ V^{T}$ .
Truncate the SVD by considering only the first $r$ columns of $U$ and $V$ , as well as the first $r$ rows and columns of $Σ$ (with the singular values ordered by size), to obtain $U_{r}$ , $Σ_{r}$ , and $V_{r}$ . Note that this is optional.
Let
$\tilde{A} = U_{r}^{T} A U_{r} = U_{r}^{T} X^{#} V_{r} Σ_{r}^{- 1} \in R^{r \times r}$ (47)
and find the eigenvalues $μ_{j}$ and eigenvectors ${\tilde{v}}_{j}$ of $\tilde{A}$ , with $\tilde{A} {\tilde{v}}_{j} = μ_{j} {\tilde{v}}_{j}$ .
Every nonzero $μ_{i}$ is a DMD eigenvalue, with the corresponding DMD mode given by
$v_{i} = μ_{i}^{- 1} X^{#} V_{r} Σ_{r}^{- 1} {\tilde{v}}_{i}$ (48)

It is common to compute the projected DMD modes, which are simply

P v_{i} = U_{r} {\tilde{v}}_{i}

, where

P = U_{r} U_{r}^{T}

is the orthogonal projection onto the first

r

POD modes of the data in

X

Note that matrix $A$ in Eq. (46) is related to operator $\exp (A Δ t)$ in Eq. (5), with $Δ t = t_{i + 1} - t_{i}$ . Hence, the eigenvalues are related by

λ_{j} = \frac{1}{Δ t} \log (μ_{j})

(49)

Using this relationship, the growth/decay rates and frequencies of the DMD modes can be inferred by examining the real and imaginary components of

λ_{j}

In addition to the DMD algorithms presented, there are a number of variants. Some of them are discussed below in the Notes. For further details, refer to the work of Rowley and Dawson [2], Tu et al. [107], and Kutz et al. [4].

2. Notes

Computationally Efficient Algorithms. A fast method to perform DMD in real time on large datasets was recently proposed by Hemati et al. [116]. A library of tools for computing variants of DMD is available.§§§ A parallelized implementation of DMD (as well as other system identification/modal decomposition techniques) was described by Belson et al. [117].¶¶¶

Sparsity. It is often desirable to represent a dataset sparsely, i.e., in terms of a small number of DMD modes. Because DMD modes are not orthogonal and have no objective ranking (as POD modes have), this is not an easy task. A number of variants of DMD has been proposed to provide such a sparse representation: for instance, the optimized DMD [108] and the optimal mode decomposition [118] are two such variants. An algorithm called sparsity-promoting DMD has also been proposed [119].**** Note that finding a sparse representation of a dataset is different from applying DMD to sparse data. On this latter problem, it has been shown that compressed sensing can be used to apply DMD to temporally [120] or spatially [121,122] sparse data, and that randomized methods can be used for enhanced computational efficiency [123].

Systems with Inputs/Control. Some systems have external inputs that are known (unlike process noise): for instance, the input would be the signal to the actuator in a flow control situation. DMD has recently been extended to handle such systems [124].

Nonlinear Systems. As mentioned previously, the connections with the Koopman operator give meaning to DMD when applied to a nonlinear system; however, when nonlinearity is important, DMD often gives unreliable results, unless a sufficiently large set of measurements is used in the snapshots. For instance, it can help to include the square of the velocities in each snapshot, in addition to the velocities themselves. Recent work [125] clarifies these issues and extends DMD to allow for a better global approximation of Koopman modes and eigenvectors, as well as Koopman eigenfunctions. Data-driven approximation of Koopman eigenfunctions can be used as a set of universal coordinates, through which disparate datasets may be related [126].

Connections with Other Methods. In many situations, DMD is equivalent to a discrete Fourier transform [105,108]. In addition, DMD shares algorithmic similarities with a number of other techniques, such as the eigensystem realization algorithm [107] and linear inverse modeling [104,107] (a technique used in climate science). If DMD is applied to linearized flow about its steady state, the extracted DMD modes capture the global modes (also see Sec. VII).

3. Strengths and Weaknesses

Strengths.

DMD does not require any a priori assumptions or knowledge of the underlying dynamics. DMD is an entirely data-driven analysis.
DMD can be applied to many types of data, or even concatenations of disparate data sources.
Under certain conditions, DMD gives a finite-dimensional approximation to the Koopman operator, which is an infinite-dimensional linear operator that can be used to describe nonlinear dynamics (see Sec. VI).
DMD modes can isolate specific dynamic structures (associated with a particular frequency).
DMD has proven to be quite customizable, in the sense that a number of proposed modifications address the weaknesses outlined in the following.

Weaknesses.

It can be difficult (or at least subjective) to determine which modes are the most physically relevant (i.e., there is no single correct way to rank eigenvalue importance, unlike other methods such as POD).
DMD typically requires time-resolved data to identify dynamics, although extensions exist [107,120].
If DMD is used for system identification (without any modifications), the resulting model will be linear.
DMD can be unreliable for nonlinear systems. In particular, for a nonlinear system, we must be careful to choose a sufficiently rich set of measurements (in each snapshot). Without care, the connection with the Koopman operator and the underlying dynamical system may be lost. Furthermore, for nonlinear systems with complex (e.g., chaotic) dynamics, there are further complications that could limit the applicability of DMD and related algorithms.
The outputs of DMD can be sensitive to noisy data, which were shown empirically [127] and analytically [128]. The effect of process noise (that is, a disturbance that affects the dynamics of the system) has also been investigated [129]. However, there are algorithms that are robust to sensor noise [128,130].
DMD should generally be used only for autonomous systems (i.e., the governing equations should have no time dependence or external inputs, unless these are explicitly accounted for [124]).
DMD modes are not orthogonal. This has a number of drawbacks: for instance, if the modes are used as a basis/coordinate system for a reduced-order model, the model will have additional terms due to the spatial inner product between different modes being nonzero. Note that a recent variant, recursive DMD [131], considers orthogonalized DMD modes.
DMD relies fundamentally on the separation of variables, as does POD, and hence does not readily extend to traveling wave problems.
DMD does not typically work well for systems with highly intermittent dynamics. However, multiresolution [132] and time-delay [133] variants show promise for overcoming this weakness.

B. Illustrative Example

1. Jet in Crossflow

We show an example from Rowley et al. [105], where DMD is applied to three-dimensional jet-in-crossflow direct numerical simulation (DNS) data. A typical snapshot of the complex flow is visualized in Fig. 9 (top) using the $λ_{2}$ criterion. The results of applying DMD to a sequence of 251 snapshots are shown in Fig. 9 (bottom). The bottom-left plot shows the frequency and amplitude of each DMD mode, with two modes having large amplitudes visualized in the bottom plots. Note that these modes capture very different flow structures, with each having a different characteristic frequency of oscillation identified by DMD.

Fig. 9
Snapshot of jet-in-crossflow data used for DMD [105], shown with $λ_{2}$ contour in the top plot. The frequency and amplitude of the DMD modes are shown in the bottom left plot, whereas the bottom right figures show the modes corresponding to (dimensionless) frequencies of 0.0175 and 0.141, visualized with contours of streamwise velocity. Reprinted with permission from Cambridge University Press.

2. Canonical Separated Flow with Control

Let us present a second example [134] of DMD (and POD) analysis performed on three-dimensional separated turbulent flow over a finite-thickness plate with an elliptical leading edge at $R e = 100, 000$ . The plate is aligned with the freestream, and separation is induced by imposing a steady blowing/suction boundary condition above the plate. The flowfield is obtained from LES, and the wake dynamics in the example is modified by a synthetic jet actuator placed on the top surface of the wing. First, we consider the dominant and secondary POD modes, shown in Fig. 10 (left). We observe that some interactions between the actuation input (at the 0.6 chord location) and the wake are captured by these energetically dominant modes. Also shown on the right are the DMD modes corresponding to the actuation frequency and its superharmonic. Although the fundamental DMD mode is similar to the dominant POD mode, the secondary mode is able to clearly identify the synchronization of the actuator input with the downstream wake. In contrast, POD modes comprise spatial structures having a distribution of temporal frequencies that makes pinpointing the POD mode to a specific frequency difficult.

Fig. 10
Comparison of POD and DMD modes for separated three-dimensional turbulent flow over a flat-plate wing with synthetic-jet actuation on the top surface [134]. The fundamental and secondary DMD modes correspond to the actuation frequencies of 4.40 (blue) and its superharmonic (red). Reprinted with permission from J. Tu.

C. Outlook

Although DMD has quickly become a widely used method for analyzing fluid flow data, there remain challenges and applications that have yet to be fully addressed. Connections between DMD and the Koopman operator indicate its potential to model (and control) nonlinear systems. However, choosing suitable observables to give an accurate finite-dimensional approximation to the Koopman operator generally remains an open question. Algorithmic improvements that have been and should continue to be made will allow DMD to remain a practical tool to analyze increasingly large fluid flow datasets.

VI. Koopman Analysis

Koopman analysis provides an important alternative perspective to classical dynamical systems theory [135] for the description of complex systems. The Koopman operator was introduced in the early 1930s [136] to show how the dynamics of Hamiltonian systems could be described by an infinite-dimensional linear operator on the space of observable functions of the state of the system. Recently, this theory has been at the center of efforts for the data-driven characterization of complex systems [137]. There is particular interest in obtaining finite-rank approximations to the linear Koopman operator. Practically, dynamic mode decomposition (Sec. V) is the most implemented numerical framework for Koopman mode decomposition in fluids [103–105]. There have been a number of excellent in-depth reviews on Koopman analysis recently [106,138].

Koopman analysis not only provides a set of modes, of which DMD modes are a subset, but also a set of eigenvalues that determine the modal dynamics and a set of Koopman eigenfunctions that serve as intrinsic observable functions. In many contexts, Koopman analysis provides an equation-free method [139] to extract coherent structures and dynamics from data measurements of a complex system; these coherent structures are related to POD modes in fluids [1].

A. Description

1. Algorithm

The inputs are a nonlinear dynamical system:

x_{i + 1} = f (x_{i})

(50)

The state

x

typically lives in a vector space, such as

R^{n}

C^{n}

, although the theory is defined more generally on curved manifolds. Although the aforementioned dynamical system is written in discrete time, the theory also holds for continuous dynamical systems.

The output is the infinite-dimensional linear operator $U_{t}$ that describes the evolution of scalar observables $g (x)$ on state space:

U_{t} g = g \circ f

(51)

Here,

\circ

denotes function composition.†††† This operator advances all functions in the Hilbert space of measurement functions, and it holds for all states

x

. An implication is that measurements are advanced in time:

U_{t} g (x_{i}) = g \circ f (x_{i}) = g (x_{i + 1})

(52)

An alternative description, related to DMD, is as follows:

The inputs are data snapshots of observables on the state space of a dynamical system.

The outputs are the modal decomposition and linear dynamical system describing modal evolution (see Sec. V; DMD).

2. Notes

Approximation by DMD. A major applied goal of the Koopman analysis is the identification of eigenfunctions from data. These eigenfunctions provide an intrinsic coordinate system, along which the dynamics appear linear. The Koopman mode decomposition is usually approximated via the dynamic mode decomposition, described in Sec. V. DMD is essentially a linear regression of data onto dynamics, and augmenting the measurement data with nonlinear functions of the state may enable a more accurate approximation of the Koopman operator [125] and resulting eigenfunctions. Another promising algorithm computes Koopman eigenfunctions based on the diffusion operator [140].

Hamiltonian Systems. For Hamiltonian systems, the Koopman operator is unitary, meaning that the inner product of any two observable functions remains the same before and after the operator. Unitarity is a familiar concept, as the discrete Fourier transform and the POD basis both provide unitary coordinate transformations. In the original paper of Koopman [136], connections were drawn between the Koopman eigenvalue spectrum and both conserved quantities and integrability.

3. Strengths and Weaknesses

Strengths:

The Koopman analysis provides an alternative operator-theoretic perspective to dynamical systems and allows nonlinear systems to be represented and analyzed using linear techniques.
The Koopman operator is amenable to standard spectral decomposition in terms of eigenvalues and eigenvectors.
Powerful techniques in control theory apply to linear systems and may improve nonlinear control performance via Koopman linear systems [141].

Weaknesses:

By introducing the Koopman operator, we trade finite-dimensional nonlinear dynamics for infinite-dimensional linear dynamics. Dealing with infinite-dimensional operators and obtaining low-order approximations is challenging, although there are a larger number of computational techniques to analyze linear operators than there are for differential geometry on manifolds.
It may be difficult to find Koopman eigenfunctions, which define an intrinsic observable measurement coordinate system. Without Koopman eigenfunctions, it may be difficult or impossible to obtain a finite-dimensional subspace of the Hilbert space of measurement functions that remains closed under the Koopman operator. However, if discovered, a Koopman invariant subspace defines a measurement system where measurements are propagated by a finite-dimensional linear dynamical system.
The integration of control theory and Koopman analysis remains to be fully developed.

B. Illustrative Example

Simple Nonlinear Dynamical System. As an example, consider the following simple nonlinear dynamical system in two variables with a single fixed point at the origin [107,142]:

{\dot{x}}_{1} = μ x_{1}, {\dot{x}}_{2} = λ (x_{2} - x_{1}^{2})

(53)

For

λ ≪ μ < 0

, the fixed point at the origin is stable and there is a slow manifold given by

x_{2} = x_{1}^{2}

; trajectories quickly attract onto this manifold before converging to the origin (Fig. 11a). Introducing a nonlinear change of coordinates, given by

(y_{1}, y_{2}, y_{3}) = (x_{1}, x_{2}, x_{1}^{2})

, the nonlinear dynamics become linear:

[\begin{matrix} {\dot{y}}_{1} \\ {\dot{y}}_{2} \\ {\dot{y}}_{3} \end{matrix}] = \underset{K}{\underset{⏟}{[\begin{matrix} μ & 0 & 0 \\ 0 & λ & - λ \\ 0 & 0 & 2 μ \end{matrix}]}} [\begin{matrix} y_{1} \\ y_{2} \\ y_{3} \end{matrix}]

(54)

The nonlinear system in Eq. (53) and the linear Koopman system in Eq. (54) are shown in Fig. 11. The system in Eq. (54) is defined on a vector of observables $y = {(y_{1}, y_{2}, y_{3})}^{T}$ , which defines a Koopman-invariant subspace. The infinite-dimensional Koopman operator restricted to this subspace defines a $3 \times 3$ matrix $K$ . With this linear representation, it is possible to analytically predict the state of the system at any future time [i.e., $y (t) = \exp (K t) y (0)$ ] or extend textbook linear estimation and control techniques to nonlinear systems.

The preceding example may appear exceedingly simple, but there are a few direct applications of the Koopman analysis for more complex systems because it is challenging to discover the coordinate transformations that linearize the problem. In fluid dynamics, Bagheri [143] analyzed the Koopman mode decomposition of the flow past a cylinder, relating the decomposition to a POD–Galerkin with a shift mode [56]. Related work investigated the Liouville equation [144], which resulted in cluster-based reduced-order models in fluid systems [85]. In power electronics, Koopman mode decomposition has been used to model and predict instabilities [145,146]. However, DMD is used in most fluid applications to approximate the Koopman mode decomposition, although the quality of this approximation depends on the measurements used for DMD [4]. True Koopman linearization of complex systems in fluid dynamics relies on the choice of good observable functions that provide linear embeddings for the nonlinear dynamics, as in the preceding example. The extended DMD [125] provides one approach to augment DMD with nonlinear measurements of the state, potentially improving the approximation of Koopman eigenfunctions.

C. Outlook

The promise of Koopman analysis hinges on the discovery of good observable functions that provide a coordinate transformation in which measurements behave linearly. Identifying these Koopman eigenfunctions from data is both a major goal and central challenge moving forward, making it a focus of research efforts. Dynamic mode decomposition, which is a workhorse of Koopman analysis, implicitly uses linear observable functions. In fluid flows, these often take the form of direct velocity field measurements from numerical simulations or particle image velocimetry in experiments. In other words, the observable function is an identity map on the fluid flow state. This set of linear observables is often too limited to describe the rich dynamics observed in fluid systems. Recently, DMD has been extended to include a richer set of nonlinear observable functions, providing the ability to effectively analyze nonlinear systems [125]. To avoid overfitting, sparse regression is emerging as a principled technique to select active terms in the dynamical system [147], and similar techniques may be used to improve the identification of observable functions. Recent techniques also construct Koopman eigenfunctions using a regularized advection–diffusion operator [140] or approximate the Koopman operator in delay coordinates [133,148]. Discovering Koopman eigenfunctions is a central computational issue to obtain closed systems that may be used for nonlinear control and estimation, and the ultimate success of this analysis will depend on our ability to accurately and efficiently approximate these eigenfunctions from data. It is likely that developments in machine learning, such as deep learning, will continue to advance these capabilities.

VII. Global Linear Stability Analysis

Global linear stability analysis solves the linearized Navier–Stokes equations pertinent to an underlying base flow with multiple inhomogeneous spatial directions. The introduction of a modal formulation for the linear perturbations converts the linearized equations of motion (initial-value problem) to an eigenvalue problem with respect to the base flow, which needs to be an exact steady (or unsteady) solution of the equations of fluid motion. The term global is used to distinguish the analysis from the classic local linear stability theory [6,34,149], where the base flow is independent of two coordinate directions. The eigenvalue problem delivers information on the spectrum of the linearized Navier–Stokes operator at asymptotically large times. Much like in classic local linear stability theory, the eigenvalue determines whether the corresponding eigenvector grows (exponentially) in time or space and delivers information on the frequency of the linear perturbation (Sec. II and Fig. 3). The eigenvector describes the spatial structure of the linear perturbation and may be used to reconstruct the full flowfield at conditions consistent with the linear approximation.

Global modes are solutions to the two- or three-dimensional partial derivative eigenvalue problems to which the linearized Navier–Stokes equations may be recast. They describe disturbances developing upon base flows that vary inhomogeneously in two (biglobal) or three (triglobal) spatial directions [150]. Thus, global modes are the counterparts, in base flows that vary in multiple spatial dimensions, of well-known linear instabilities developing in one-dimensional base flows such as the inviscid Kelvin–Helmholtz instabilities, viscous Tollmien–Schlichting, and crossflow modes. When the two- or three-dimensional base flow contains portions of flow that could be analyzed with local theory, the spatial distribution of the global eigenfunction is seen to contain structures related to the eigenfunction delivered by local analysis, e.g., Tollmien–Schlichting waves in the boundary layer on a flat plate [151] and on the downstream wall of an open cavity [152]. The applications of global stability analysis to study linear stability of two- and three-dimensional base flow have been reviewed recently [5] and are continuously expanding.

A. Description

1. Algorithm

The input is the steady or time-periodic laminar base flow, inhomogeneous in two or all three spatial directions, at a given Reynolds number and Mach number. The base flow can be a stable or an unstable steady state.

For the outputs, the global linear modal stability theory delivers global modes (spatial distribution pattern) and the associated growth rates and frequencies of small-amplitude perturbations about the base state.

Unlike the aforementioned POD, balanced POD, or DMD techniques, global stability analysis is not based on the snapshots of the flowfield. Global modes are found numerically by discretizing the Navier–Stokes equations, which result from the decomposition of any flow quantity $q$ into a base flow $q_{0}$ and small-amplitude perturbations $q^{'}$ (i.e., $‖ q^{'} ‖ / ‖ q_{0} ‖ ≪ 1$ ):

q (ξ, t) = q_{0} (ξ, t) + q^{'} (ξ, t)

(55)

We can then examine the growth or decay of the perturbations with respect to the base state

q_{0}

. Substitution of the flow variable in the form of Eq. (55) into the Navier–Stokes equations and neglecting the quadratic terms of small perturbation

q^{'}

yields the linearized Navier–Stokes equations.

Now, as an example, let us assume the base flow $q_{0} = q_{0} (ξ)$ to be steady and homogeneous (e.g., periodic in space) in the $ζ$ direction and perform a biglobal stability analysis. In this case, the perturbation can assume the form of

q^{'} (ξ, η, ζ, t) = \hat{q} (ξ, η) e^{i (β ζ - ω t)} + complex conjugate

(56)

where

β

is a real wave number in the

ζ

direction. The real and imaginary components of the complex number

ω = ω_{r} + i ω_{i}

, respectively, correspond to the frequency and the growth/decay rate of the amplitude function of the global mode

\hat{q} (ξ, η)

. There are two main approaches to find the global stability modes and the associated eigenvalues. Namely, they are the matrix-based approach and the time-stepping approach, which are described in the following. Various techniques have also been summarized in the review by Theofilis [5], where the relative merits of matrix-based and time-stepping techniques were discussed.

Matrix-Based Approach:

The matrix-based approach determines the global modes by solving a generalized eigenvalue problem. The substitution of the assumed form of Eq. (56) into the linearized Navier–Stokes equations results in

A (q_{0}; β) \hat{q} = ω B \hat{q}

(57)

which is written in terms of a generalized eigenvalue problem where

ω

is the eigenvalue and

\hat{q}

is the global stability mode (eigenvector). The matrix

A

is dependent on the base state

q_{0}

as well as the wave number

β

and flow parameters, such as the Reynolds number and Mach number. In general,

B

is invertible for compressible flow but is not invertible for incompressible flow due to the incompressibility constraint. For detailed discussions on the linearized Navier–Stokes equations, see the work of Schmid and Henningson [6] and Theofilis [5].

The matrices $A$ and $B$ are approximately of size $n_{grid} n_{var} \times n_{grid} n_{var}$ , where $n_{grid}$ and $n_{var}$ are the number of grid points and number of variables, respectively. If the operators are of modest size, the matrices may be formed and stored on memory. For larger-sized operators, the eigenvalues can be determined without storing the operators but by only relying on matrix–vector operations (matrix-free approach) [5,153–155]. We can also distribute the matrix entries over several processors on a cluster and use appropriate libraries (e.g., ScaLAPACK) for the linear algebra operations [156].

Time-Stepping Approach:

In the time-stepping approach [5,157–159], a code is employed to compute the linear operators of the Navier–Stokes equations (Fréchet derivative [160]), without the need to store the matrix. Because, typically, the numerical discretization of the operators does not expand the solution on bases that automatically satisfy the boundary conditions for the perturbations, spurious modes can be introduced. We need to extract the physically meaningful modes from the numerically obtained part of the eigenspectrum.

2. Notes

Computational Effort. With either the matrix-based approach or the time-stepping approach, the computational effort required to find eigenvalues and eigenfunctions can be very high. Hence, the use of sparse yet high-order spatial discretization methods can be beneficial [161,162]. If the flow of interest has symmetry, it can be exploited to reduce the computational cost of solving two- and three-dimensional eigenvalue problems [159,163,164].

Choice of Base Flow. For stable flows, the steady-state solution to the Navier–Stokes equations can be found by integrating the solution in time until all unsteadiness is eliminated. If the unstable steady state is needed, techniques such as the selective frequency damping method [165] or a Newton–Krylov-type iterative solver [166,167] can be used. We note that time-periodic flow can also be used as the base flow [168], as we will see in the first illustrative example in the following [169].

Occasionally, researchers compute global modes for base flows that are not themselves solutions of the governing equations (e.g., the mean of a turbulent flow). In these cases, the resulting modes are not associated with a question of stability of the base flow but may be useful in modeling and identifying the frequency content of large-scale coherent structures [154,170]. Additional details with a cautionary note were provided by Sipp and Lebedev [171].

Global Nonmodal/Transient-Growth Analysis. Global modes may represent either stable or unstable disturbances with respect to the base flow. When unstable modes are identified, the base flow itself is said to be asymptotically unstable, meaning that at least one perturbation will grow exponentially in time and the equilibrium is not expected to be observed in nature. As discussed in Sec. II.E, systems that are nonnormal can exhibit significant growth of linear combinations of modal perturbations, known as transient growth, which is a phenomenon that may occur even when the global modes are all stable [33,36,172]. Nonnormality of the governing linear operator implies that global transient-growth analysis may deliver qualitatively different short-time behavior of linear combinations of perturbations when modal analysis predicts only asymptotically stable perturbations. Global transient-growth analysis has unraveled initial optimal perturbations that are qualitatively different from their corresponding global eigenmodes on the cylinder [173], low-pressure turbine blades [174,175], and several stalled NACA airfoils [169].

3. Strengths and Weaknesses

Strengths:

A spectrum of (discrete and continuous) eigenmodes can be determined, especially with the matrix-based approach.
Modal analysis, based on the solution of the eigenvalue problem (see Sec. II), determines whether unstable modes exist. Exponential growth implies that flow dynamics will be asymptotically dominated by the characteristics of the most unstable eigenmode (or the least stable mode if stable).
If only stable modes are found, global nonmodal/transient-growth analysis, using a linear combination of the leading members of the eigenvalue spectrum or a solution of a related SVD problem (see Sec. II), can determine the level of energy growth of small-amplitude perturbations over a short time horizon.

Weaknesses:

Global linear stability analysis is inherently linear.
The base flow provides the spatially variable coefficients of the underlying partial derivative eigenvalue/initial-value problem. Any algorithm for the solution of the latter problems can be influenced by the quality (accuracy) of the underlying base flow.
The analysis requires a base flow that is an exact solution of the equations of motion. Although counterexamples exist in both classic [176] and global linear stability theories [177], an analysis of mean turbulent flow requires validation of the turbulence closures employed and typically aims at the prediction of the frequencies and spatial structure of coherent turbulent structures.

B. Illustrative Examples

1. Large-Scale Separation Cells on Stalled Airfoils

He et al. [169] employed a suite of matrix-forming and time-stepping techniques to reexamine large-scale separation patterns on the suction side of wings in near-stall flight conditions. The origins of these structures were explained through systematic application of primary and secondary linear global eigenvalue and transient-growth stability analyses to massively separated spanwise homogeneous laminar flows over spanwise-periodic wings of different thicknesses and cambers. At low chord Reynolds numbers, the dominant flow structure arising from either primary modal or nonmodal mechanisms is associated with two- or three-dimensional Kelvin–Helmholtz eigenmodes. As the wing aspect ratio is shortened, the stationary three-dimensional global eigenmode discovered by Theofilis et al. [178] is amplified more than the Kelvin–Helmholtz mode and stall cells may arise through this primary modal linear amplification mechanism [179]. As the Reynolds number increases, linear amplification of the two-dimensional Kelvin–Helmholtz global eigenmode leads to a time-periodic wake that, in turn, is linearly unstable with respect to two distinct classes of three-dimensional secondary (Floquet) eigenmodes, which peak at short- and long-spanwise wavelengths, respectively. The short-wavelength Floquet eigenmode is the stronger of the two at moderate and high Reynolds numbers. The three-dimensional spatial structure of the short-wavelength mode can be seen in Fig. 12a, whereas the wall-streamline pattern resulting from its linear superposition upon the underlying time-periodic base flow gives rise to stall-cell-like patterns on the wing surface, as can be seen in Fig. 12b.

Fig. 12
Three-dimensional short-wavelength instability mode determined from a global stability analysis for a NACA 4415 airfoil at $R e = 500$ and $α = 20 \deg$ : a) streamwise vorticity $ω_{ζ}$ isosurface of the short-wavelength eigenmode, and b) the corresponding wall streamlines, superposed upon the time-periodic base state [169]. Reprinted with permission from Cambridge University Press.

2. Finite-Span Open-Cavity Flow

The second example leverages the time-stepper-based approach to gain insights into the critical conditions and spatial characteristics of distinct classes of eigenmodes of incompressible flow over a rectangular three-dimensional lateral-wall-bounded open cavity. The analysis has been performed in a triglobal setting [159]. The leading traveling shear-layer mode responsible for transition in this geometry, as well as the next in significant stationary and traveling centrifugal modes, respectively, are shown in Fig. 13 at the slightly subcritical Reynolds number of $R e = 1050$ . A detailed temporal biglobal analysis of the related spanwise homogeneous open cavity has also been performed [155,161], and the effect of the presence of lateral walls on the linear instability mechanisms in this class of flows can now be quantified on the basis of global stability theory.

Fig. 13
Amplitude functions of the components of the three-dimensional perturbation velocity components of the leading eigenvectors in lateral wall-bounded open rectangular cavity flow at $R e = 1050$ [159]. Reprinted with permission from Cambridge University Press.

C. Outlook

With the availability of enhanced computational resources, the use of temporal and spatial biglobal stability analyses is becoming prevalent to examine a range of fluid flows. An area in which the theory is expected to become particularly useful is laminar–turbulent transition prediction in hypersonic flow, in support of increasing efforts to understand and control flow phenomena critical to the design of next-generation vehicles. Unlike typical incompressible or supersonic conditions at which turbulent flow prevails in flight, hypersonic flow is mostly laminar and the underlying base flows are exact solutions of the equations of motion, which can be computed either using a continuum assumption [180] or by employing direct simulation Monte Carlo methods [181 –184]. As an example, a linear global instability analysis has recently addressed hypersonic flow over the HIFiRE-5 (Hypersonic International Flight Research Experimentation) elliptic cone [164,185]. The spatial biglobal eigenvalue problem has been solved, and four classes of (leading) centerline, attachment-line, crossflow, and second Mack modes have been recovered as two-dimensional global mode amplitude functions, underlying the potential of the theory to provide predictions of laminar–turbulent transition at hypersonic speeds.

In addition to taking advantage of the symmetry present in flows of interest, we can incorporate the parabolized assumption [186] into a global stability analysis. For instance, three-dimensional inhomogeneous base flows that depend strongly on two directions and weakly on the third spatial direction are an extension of the nonlocal stability analysis, which is based on the parabolized stability equations. Examples of successful application include the system of trailing vortices behind an aircraft wing [187–189], the wake of an isolated roughness element in supersonic flow [190], streaks in a boundary layer [191], corner flows [192], and models of duct intakes [193]. A computationally expensive triglobal stability analysis is emerging to provide insights into the influence of three-dimensionality on the stability of flows [157–159]. It is anticipated that there will be growing use of these approaches to analyze increasingly complex flows.

As a companion to the global stability analysis, the adjoint global stability analysis [194] can provide efficient responses to questions on flow receptivity, sensitivity, and control. Global modes are the right eigenvectors of the eigenvalue problem, whereas the left eigenvectors are the adjoint global modes. The adjoint linearized Navier–Stokes equations, which are at the heart of the linear instability analysis and linear control methodologies [194], can be solved with some added effort to the global eigenvalue problem. At its simplest implementation, a single solution of the global eigenvalue problem delivers the eigenvalues, as well as linear perturbations and their adjoint variables as the right and left eigenvectors, respectively.

VIII. Resolvent Analysis

The resolvent is defined in the context of fluid mechanics as the linear operator relating an input forcing to a linear system to the output. For the complete incompressible Navier–Stokes equations, the output is the (divergence-free) velocity field and the input to the resolvent (linear dynamics) is provided by the nonlinear term, as shown in Fig. 14. The resolvent analysis is therefore a complement to eigenvalue and global stability analyses, which are linear analyses with perturbations assumed to be sufficiently small such that the nonlinear advective forcing term (second order in the perturbation) can be neglected. More generally, inputs to the resolvent analysis can include actuation and other kinds of forcing terms added to the governing equations, and outputs can be any observable.

Fig. 14
Schematic of a divergence-free projection of the resolvent formulation of the incompressible Navier–Stokes equations, showing the nonlinear term $f$ providing input forcing to the linear dynamics (described by the resolvent operator), with velocity output, $u$ [i.e., $q^{'}$ in Eq. (58)] [195]. Reprinted with permission from AIP Publishing.

The resolvent analysis relies on the pseudospectrum of the linearized Navier–Stokes operator rather than the spectrum itself (see Sec. II.E). For stable base flows, the damped modes of the linear system can be superposed to construct the response to stationary (real-frequency) inputs. Even for these stable flows, when the operator is nonnormal (most fluid systems linearized about spatially nonuniform flows give rise to nonnormal operators), stationary inputs can be amplified. For these systems, there are typically a limited number of highly amplified inputs, implying that the resolvent can be effectively approximated at low rank. An SVD of the resolvent (in the discrete setting) can be used to identify these inputs and their corresponding outputs, as well as to provide a low-rank approximation of the input–output dynamics of the full system. High-gain inputs are useful for flow control efforts by showing where actuation will produce the largest effect.

A. Description

1. Algorithm

The inputs are a real-valued frequency $ω$ and a linear time-invariant input–output system (operator) including its boundary conditions: for example, the linearized Navier–Stokes equations with an appropriate assumed base or mean velocity profile $\bar{q}$ . These must typically have been discretized in space using a finite difference, finite element, spectral, or other numerical method.

The outputs are a ranked set of modes (total number equal to the number of degrees of freedom of the discretized system) that represent the spatial pattern of inputs, outputs, and the positive gain (amplification) from input to output. High-gain modes represent stationary (statistically steady) inputs that are most amplified at the specified frequency.

Consider $q = \bar{q} + q^{'}$ and the spatially discretized governing equations for the perturbation variable expressed in the compact form

M \frac{\partial q^{'}}{\partial t} + A q^{'} = f (q^{'})

(58)

where the vector

q^{'}

represents the fluctuating discrete flow variables; the (typically sparse) matrices

M

and

A = A (\bar{q})

represent the linearization about a time-invariant flow

\bar{q}

; and

f

denotes the remaining nonlinear terms, as well as any additional forcing terms added to the equations. Such terms can include, for example, stochastic forcing (synthetic turbulence) added to inflow regions of spatially developing flows. For statistically steady (stationary) flows, we may work in the frequency domain. By expressing

q^{'} = \hat{q} e^{i ω t} and f = \hat{f} e^{i ω t}

(59)

for temporal frequency

ω \in R

, Eq. (58) can be written as

(i ω M + A) \hat{q} = \hat{f}

(60)

Proceeding formally, we write

\hat{q} = R \hat{f}

, where

R = {(i ω M + A)}^{- 1}

is the resolvent operator. For a particular frequency

ω

, the maximum amplification

G

(defined in a chosen norm

{\hat{q}}^{*} Q \hat{q}

, where

Q

is a positive definite matrix) of any forcing may be found:

G = \max_{\hat{f}} \frac{{\hat{q}}^{*} Q \hat{q}}{{\hat{f}}^{*} Q \hat{f}} = \max_{\hat{f}} \frac{{\hat{f}}^{*} R^{*} Q R \hat{f}}{{\hat{f}}^{*} Q \hat{f}}

(61)

The solution to this maximization problem can be found through the SVD of

R

. The largest singular value represents the (squared) gain, and the left and right singular vectors represent the optimal response and forcing, respectively.

In the context of the incompressible Navier–Stokes equations, the resolvent $R$ is the linear operator that relates velocity perturbation $u^{'}$ and pressure perturbation $p^{'}$ to the nonlinear term $f$ :

q^{'} = [\begin{matrix} u^{'} \\ p^{'} \end{matrix}] = R f

(62)

Note that a base velocity profile is required to formulate the resolvent and to define the perturbations.

Spatial homogeneity can be exploited to reduce the computational burden by first decomposing the linear system into the Fourier basis in the homogeneous spatial directions. This is also useful for identifying dynamics associated with individual wave number components of interest in the data. For simplicity, we consider here an example flow that is homogeneous in the streamwise and spanwise directions and statistically stationary in time, such that the wave number/frequency triplet $k = (k_{x}, k_{z}, ω)$ associated with oblique downstream traveling waves can be defined. Accordingly, we have

R_{k} = {(- i ω [\begin{matrix} I \\ 0 \end{matrix}] - [\begin{matrix} L_{k} & - \nabla \\ \nabla^{T} & 0 \end{matrix}])}^{- 1} [\begin{matrix} I \\ 0 \end{matrix}]

(63)

where

L_{k}

is the linear Navier–Stokes (advection and diffusion) operator evaluated at

k

\nabla

and

\nabla^{T}

are the gradient and divergence operators, and

I

is the identity operator. As such, the resolvent operator can be obtained by a relatively minor extension of code for the more common Orr–Sommerfeld–Squire analysis. Of course, the system may be formulated in terms of either velocity–vorticity or primitive variables. For an analogous derivation of the resolvent for compressible flow, see the work of Jeun et al. [196].

Two closely related formulations of the resolvent analysis are common in the literature. In the input–output formulation, an output quantity is defined as the state vector premultiplied by a matrix that is chosen to extract an observable of interest [196]. Similarly, the forcing is premultiplied by a (not necessarily) different matrix to restrict the input–output relation to forcings of interest. In the linear frequency response approach [197], the gain is defined as the quotient of the input and output in terms of their energy defined through (not necessarily different) norms.

For a given $k$ , the SVD procedure described previously in Sec. II returns right singular vectors corresponding to the most dangerous (most amplified) inputs or forcing modes, and left singular vectors corresponding to the associated response modes, with gains given by the corresponding singular values. The example code for a resolvent analysis of turbulent pipe flow in primitive variables is available online for general use [198]. The analysis proceeds, depending on the specific objective. For example, the globally most amplified input can be sought, the forcing and response modes can be compared with observations, or data can be projected onto the resolvent basis [195,199].

2. Notes

Origin. The resolvent is a familiar construction in the study of forced, ordinary differential equations and, more generally, linear operators and their spectra. It arises in control theory and in eigenvalue/eigenvector perturbation analysis. Identified by Schmid and Henningson [6], its use as a tool to study structures in transitional and turbulent flows seems to have been stimulated by Farrell and Ioannou [200], who considered a stochastically forced Navier–Stokes system, and follows contemporaneous work on understanding transient growth in stable base flows with pseudospectra [35].

Choice of Base Flow. Similar to the linear global stability analysis discussed in the previous section, the choice of the base flow $q_{0}$ is important in resolvent analysis because operator $A$ and forcing input $f$ are established with respect to the base state $q_{0}$ . If the (stable or unstable) steady-state solution to the Navier–Stokes equations is available, it can be used as the base flow [6,35] and the nonlinear term can be assumed to be negligible within the context of linear analysis. In this case, the input forcing must be provided from an external source. We also note that a resolvent analysis has been performed about the time-averaged flow $\bar{q}$ . In such a case, the nonlinear term is not negligible but is considered as an internal source of the forcing mechanism [195], as illustrated in Fig. 14. For this reason, a resolvent analysis has been used to examine how fluctuations including those from nonlinear effects are amplified or attenuated with respect to the time-averaged flow. When applying a resolvent analysis to the temporal mean of inherently unsteady or turbulent flows, it is advisable to examine whether the operator $A = A (q_{0})$ has stable eigenvalues to separate the forced response from the unforced response.

Implementation in Different Geometries. The use of a resolvent analysis is increasing in fluid dynamics to study both transitional and fully turbulent flows. Flows that have been analyzed in terms of the most amplified (forced) modes include, among others, canonical flows such as boundary layers, Couette and Poiseuille flows [199,201], jets [196,197,202,203], backward-facing steps [204], and cavity flows [205].

Computational Costs. Resolvent modes can be routinely and efficiently computed for base flows that are homogeneous in two spatial dimensions (one-dimensional base flow). A resolvent analysis for flows that are homogeneous in one spatial dimension (two-dimensional base flow) incurs significantly more computing costs [206]. With current computational resources, fully inhomogeneous flows represent an unmet challenge, although flows with narrowband perturbations can be tackled [205].

Wall Boundary Conditions and Control. A resolvent analysis admits any linear control formulation through the boundary conditions. For example, passive (a compliant surface [207]), open-loop (dynamic roughness [208]), and active (opposition control [209]) control inputs at the wall have been investigated in wall turbulence.

Wave Packets and Coherent Structures. Although formally a tool for linear systems, a resolvent analysis can provide insights into coherent structures in turbulent flow by analyzing the response of the system associated with the linear resolvent operator formed using the turbulent mean flow. This is not necessarily equivalent to linearization because the nonlinear forcing is required to sustain an otherwise stable system. Note that the mean velocity itself can be recovered via an extension of the analysis shown in Fig. 14. Resolvent modes can be identified, and the corresponding inputs can be thought to represent how such coherent structures are forced through nonlinear (triadic) interactions among modes at other frequencies, or through other stochastic inputs to the system (for example, noise added at an inflow in a DNS to excite turbulence). A particular success has been the recovery of packets of hairpin vortices in wall turbulence from resolvent modes [210].

The resolvent-based approach is also particularly useful to analyze convectively unstable configurations that exhibit no intrinsic dynamics. Despite the absence of an unstable linear global mode (see Sec. VII) in such cases, coherent structures are often observed in such flows. A typical example is a turbulent jet under most operating conditions. Here, continuous forcing through the background turbulence sustains large-scale coherent structures in the jet shear layer over a wide range of frequencies. The resolvent analysis allows identification of these wave packets as optimal responses, and it associates each response mode to its corresponding optimal input via the gain of the input–output relation (see second example in the following).

3. Strengths and Weaknesses

Strengths.

A resolvent analysis can use the mean flow as the base state (instead of the exact solution to the Navier–Stokes equations), even in the case of turbulent flows.
When resolvent analysis is applied to turbulent mean flows, large-scale coherent structures can be interpreted as modal solutions that are sustained through forcing by the turbulent background.
The resolvent operator can be formulated to describe the linear dynamics associated with the governing Navier–Stokes equations.
The resolvent analysis identifies the form of the most amplified inputs, the corresponding output, and the gain.
The resolvent operator is often low rank, with an underlying mathematical structure that can be used to understand the coherent flow structure and permit exploitation of state-of-the-art matrix approximation techniques.

Weaknesses. A resolvent analysis identifies the most amplified inputs and outputs, whereas the product of the amplification (singular value) and the weight of the nonlinear forcing determines the observations. Information on the nonlinear forcing term is therefore required to predict the most energetic observed modes and synthesize nonlinear models.

B. Illustrative Examples

1. Turbulent Channel Flow

We show in Fig. 15 an example of the typical form of the first (i.e., most amplified) resolvent response modes in turbulent channel flow. Here, a left- and right-going pair of modes (i.e., $\pm k_{ζ}$ ) are considered with wave numbers and wave speeds representative of the very large-scale motions in wall turbulence. For this mode, the amplitudes of the three velocity components have relative magnitudes $(| u |, | v |, | w |) = (1, 0.05, 0.17)$ , i.e., the turbulent kinetic energy of the mode is mostly associated with the streamwise velocity fluctuations. The isocontours of velocity correspond to a quasi-streamwise vortex (inclined roll) structure in the cross-stream plane that is familiar from observations of fully developed turbulence at all scales, but it is especially reminiscent of the very large-scale motions, suggesting that the resolvent mode is a useful simplified model for that structure at least.

2. Turbulent Jet

In Fig. 16a, the gain of first five resolvent modes associated with a turbulent jet is given as a function of the Strouhal number $S t$ (dimensionless frequency in acoustic units). The gain of the leading resolvent mode clearly dominates over a large frequency interval. Two features of the gain curve are of special physical interest. The sharp peak at $S t \approx 0.4$ is associated with an acoustic resonance in the potential core [68] and the underlying broadband response peak with the preferred mode of forced jets [197]. The forcing structure at $S t = 0.5$ in Fig. 16b reveals that the response depicted in Fig. 16c is most efficiently sustained by an upstream forcing distribution that uses the Orr mechanism in the shear layer and direct acoustic forcing in the potential core to achieve optimality. The optimal response closely resembles the coherent structure shown in Fig. 16d educed from a LES database using spectral POD (see Sec. III) based on the spectral cross-correlation matrix.

C. Outlook

The use of a resolvent analysis to develop parametric sensitivity studies is an especially promising direction for future research. The utility of the resolvent approach in developing dynamical models rests on identifying important nonlinear interactions determining the correct weighting of the nonlinear forcing. These weights can be obtained or approximated in various ways: for example, by direct calculation of the forcing from full-field snapshots [202,211]; projection of the snapshots onto resolvent response modes [212–214]; scaling arguments [215]; or direct solution of the quadratic programming problem defining the coupling between response modes and wave numbers [195], which is the subject of ongoing work.

Resolvent modes have recently been shown to provide an efficient basis to represent exact coherent solutions of the Navier–Stokes equations [214], suggesting that the analysis may both provide an inexpensive approximate continuation method within families of solutions and aid the search for new exact solutions. Formal connections between the resolvent analysis, Koopman analysis, DMD, and exact nonlinear traveling wave solutions have recently been outlined [216]. This is out of the scope of this paper, but it may be of interest to expert practitioners.

IX. Conclusions

An overview of the modal decomposition/analysis techniques that are widely used to examine a variety of fluid flows was provided. Specifically, the POD, balanced POD, DMD, Koopman analysis, global linear stability analysis, and resolvent analysis were presented. The modal structures extracted by these techniques could shed light on different aspects of the flowfield with some shared similarities. For the data-based methods, the POD analysis captures the most energetic modes and the balanced POD analysis captures the most controllable and observable modes, providing balancing and adjoint modes. DMD (and Koopman analysis) extracts the dynamic modes along with their growth rates and frequencies from the flowfield data. Although POD and DMD methods can use flowfield data from numerical simulation and experimental measurements, the balanced POD requires adjoint flow data, which makes it applicable to numerical studies only. Operator-based analysis techniques of the Koopman analysis, global linear stability analysis, and resolvent analysis were also presented. Koopman analysis rigorously connected DMD to nonlinear dynamical systems. The global linear stability and resolvent analyses can be used to examine growth or decay characteristics of perturbations with respect to a given base or mean flow. These two approaches require discretized operators from the Navier–Stokes equations to perform the modal stability analysis.

The examples shown in this paper focused mostly on fluid flow analysis. However, modes obtained from these techniques can also be used to develop reduced-order models, which can capture the dynamics of the flow with significantly lower computational costs. Such models can be useful in closed-loop flow control [54,59,84,85,90,93,217] as well as aerodynamic design [218,219]. A large number of references were provided throughout the paper so that readers can seek additional insights as needed. It is hoped that this document can serve as a stepping stone for the readers to become familiar with various modal analysis techniques, analyze a variety of complex flow physics problems, and further advance the developments of these modal analysis techniques.

P. GiviAssociate Editor

* Video recordings of this session have been made available by AIAA on https://www.youtube.com/user/AIAATV/playlists [retrieved 10 October 2017].

† If the square matrix $A$ is real, its eigenvalues are either real or come in complex conjugate pairs.

‡ In such a case, $A$ is called diagonalizable or nondefective. If $A$ is defective, we have $A = V J V^{- 1}$ , with $J$ being the canonical Jordan form [31,33].

§ The discretization of the incompressible Navier–Stokes equations can yield a singular $B$ because the continuity equation does not have a time derivative term. Also see Sec. VII.

¶ Unitary matrices $U$ and $V$ satisfy $U^{*} = U^{- 1}$ and $V^{*} = V^{- 1}$ , with $*$ denoting the conjugate transpose.

** For additional details, see the documentation available on http://www.mathworks.com [retrieved 10 February 2017].

†† Library available online at http://www.netlib.org/lapack [retrieved 10 February 2017].

‡‡ Library available online at http://www.netlib.org/scalapack/ [retrieved 10 February 2017].

§§ Library available online at http://www.caam.rice.edu/software/ARPACK [retrieved 10 February 2017].

¶¶ Because we reserve the symbol $x$ to denote the data vector, we use $ξ$ to represent the spatial coordinates in this paper.

*** Precisely speaking, the covariance matrix is defined as $R \equiv X X^{T} / m$ or $X X^{T} / (m - 1)$ . For clarity of presentation, we drop the factor $1 / m$ and note that it is lumped into the eigenvalue $λ_{k}$ .

††† For the sake of discussion, we consider here that the flowfield data to be placed on a uniform grid such that scaling due to the size of the cell volume does not need to be taken into account. In general, the cell volume for each data point needs to be included in the formulation to represent this inner product (volume integral). Consequently, the covariance matrix [Eq. (24)] should be written as $R \equiv X X^{T} W$ , where $W$ holds the spatial weights. The matrix $X^{T} X$ that later appears in Eq. (30) for the method of snapshot would similarly be replaced by $X^{T} W X$ .

‡‡‡ The matrices $ϕ$ and $Ψ$ are orthonormal, i.e., $ϕ^{T} ϕ = ϕ ϕ^{T} = I$ and $Ψ^{T} Ψ = Ψ Ψ^{T} = I$ .

§§§ Code available at https://github.com/cwrowley/dmdtools [retrieved 10 February 2017].

¶¶¶ Code available at https://github.com/belson17/modred [retrieved 10 February 2017].

**** Code available at http://www.ece.umn.edu/users/mihailo/software/dmdsp/ [retrieved 10 February 2017].

†††† For example, note that $\sin (x) \circ x^{2} = \sin (x^{2})$ .

Appendix: Data Matrix

We present an example of formatting the flowfield data into a matrix form $X$ in preparation to perform the modal decomposition analysis. For simplicity, we assume that the grid is structured with uniform spacing and the data are sampled with a constant time step. Here, we take data from a two-dimensional velocity field $q = (u, v)$ and construct $X$ . In general, the data to be examined can be from an unstructured or nonuniform grid, if scaled appropriately, or processed to reside on a uniform grid. The formulation discussed here can be easily extended to scalar variables or three-dimensional vector fields.

Consider the velocity field $q$ made available from the simulation or experiment:

q_{i j} = (u_{i j}, v_{i j}) = (u (ξ_{i}, η_{j}), v (ξ_{i}, η_{j}))

(A1)

where the spatial coordinates

(ξ_{i}, η_{j})

with

i = 1, \dots, n_{ξ}

and

j = 1, \dots, n_{η}

are for a collocated setup, as illustrated in Fig. A1. To construct the data matrix

X

, we need to organize the data in the form of a column vector at each time instance. An example of stacking the snapshot data to form a column vector

x (t)

at time

t = t_{k}

can be written as

((A2)) Note that the preceding stacking and unstacking operations are for a two-dimensional vector, which yield a column vector size of

n = 2 n_{x} n_{y}

. The unstacking operation is useful for visualizing the decomposition modes once decomposition is performed on the flowfield data. For scalar data, we only need to consider a particular component. For three-dimensional vector data, we can append the third component to the preceding formulation. Once the column data vectors are available at the desired time levels, we can collect them into a data matrix:

X = [\begin{matrix} x (t_{1}) & x (t_{2}) & \dots & x (t_{m}) \end{matrix}] \in R^{n \times m}

(A3)

In general, this data matrix is a tall and skinny matrix for fluid flows (

n ≫ m

). With this data matrix at hand, we are ready to perform modal decompositions through the POD, balanced POD, and DMD.

Fig. A1
Example of discrete velocity data on a uniform collocated grid.

Acknowledgments

This paper was one of the major outcomes from the AIAA Discussion Group (Fluid Dynamics Technical Committee) entitled “Modal Decomposition of Aerodynamics Flows” organized by Kunihiko Taira and Douglas R. Smith [U.S. Air Force Office of Scientific Research (AFOSR) Program Officer, Unsteady Aerodynamics and Turbulent Flows]. D. R. Smith has been instrumental in creating, organizing, and supporting this discussion group at every stage. He, in fact, has been one of the contributors to the present overview paper from conception to final editing. We gratefully acknowledge the contributions made by D. R. Smith to this document and for his continued support of much of this work. The authors also thank the fruitful discussions with the members of the discussion group and greatly acknowledge the generous support from the following agencies: U.S. Air Force Office of Scientific Research, Army Research Office, Defense Advanced Research Projects Agency, U.S. Department of Energy, Deutsche Forschungsgemeinschaft, the European Office of Aerospace Research and Development, the National Science Foundation, and the Office of Naval Research. The authors acknowledge the anonymous referees for providing insightful comments on the paper. K. Taira is grateful to his research group members for providing extensive feedback throughout the preparation of this paper.

References

[1] Holmes P., Lumley J. L., Berkooz G. and Rowley C. W., Turbulence, Coherent Structures, Dynamical Systems and Symmetry, 2nd ed., Cambridge Univ. Press, New York, 2012. Crossref Google Scholar
[2] Rowley C. W. and Dawson S. T. M., “Model Reduction for Flow Analysis and Control,” Annual Review of Fluid Mechanics, Vol. 49, No. 1, 2017, pp. 387–417. doi:https://doi.org/10.1146/annurev-fluid-010816-060042 Crossref Google Scholar
[3] Kutz J. N., Data-Driven Modeling & Scientific Computation: Methods for Complex Systems & Big Data, Oxford Univ. Press, New York, 2013. Google Scholar
[4] Kutz J. N., Brunton S. L., Brunton B. W. and Proctor J. L., Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, SIAM, Philadelphia, PA, 2016. Crossref Google Scholar
[5] Theofilis V., “Global Linear Instability,” Annual Review of Fluid Mechanics, Vol. 43, No. 1, 2011, pp. 319–352. doi:https://doi.org/10.1146/annurev-fluid-122109-160705 Crossref Google Scholar
[6] Schmid P. J. and Henningson D. S., Stability and Transition in Shear Flows, Springer, New York, 2001. Crossref Google Scholar
[7] Taira K., “Proper Orthogonal Decomposition in Fluid Flow Analysis: 1. Introduction,” Journal of Japan Society of Fluid Mechanics (Nagare), Vol. 30, No. 2, 2011, pp. 115–123 (in Japanese). Google Scholar
[8] Taira K., “Proper Orthogonal Decomposition in Fluid Flow Analysis: 2. Applications,” Journal of Japan Society of Fluid Mechanics (Nagare), Vol. 30, No. 3, 2011, pp. 263–271 (in Japanese). Google Scholar
[9] Samimy M., Breuer K. S., Leal L. G. and Steen P. H. (eds.), A Gallery of Fluid Motion, Cambridge Univ. Press, New York, 2003. Google Scholar
[10] Brown G. L. and Roshko A., “On Density Effects and Large Structures in Turbulent Mixing Layers,” Journal of Fluid Mechanics, Vol. 64, No. 4, 1974, pp. 775–816. doi:https://doi.org/10.1017/S002211207400190X Crossref Google Scholar
[11] Strouhal V., “On One Particular Way of Tone Generation,” Annalen der Physik und Chemie (Leipzig), Series 3, Vol. 241, No. 10, 1878, pp. 216–251. doi:https://doi.org/10.1002/(ISSN)1521-3889 Crossref Google Scholar
[12] Rayleigh L., “Acoustical Observations,” Philosophical Magazine and Journal of Science, Vol. 7, No. 42, 1879, pp. 149–162. doi:https://doi.org/10.1080/14786447908639584 Crossref Google Scholar
[13] Bénard H., “Formation de Centres de Giration à L’Arrière d’un Obstacle en Mouvemènt,” Comptes Rendus de l’Académie des Sciences, Vol. 147, 1908, pp. 839–842. Google Scholar
[14] von Kármán T., “Über den Mechanismus des Widerstandes, den ein Bewegter Körper in Einer Flüsseigkeit Erfährt,” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, Vol. 1911, 1911, pp. 509–517. Google Scholar
[15] Taneda S., “Experimental Investigation of the Wakes Behind Cylinders and Plates at Low Reynolds Numbers,” Journal of the Physical Society of Japan, Vol. 11, No. 3, 1956, pp. 302–307. doi:https://doi.org/10.1143/JPSJ.11.302 Crossref Google Scholar
[16] Coutanceau M. and Bouard R., “Experimental Determination of the Main Features of the Viscous Flow in the Wake of a Circular Cylinder in Uniform Translation. Part 1. Steady Flow,” Journal of Fluid Mechanics, Vol. 79, No. 2, 1977, pp. 231–256. doi:https://doi.org/10.1017/S0022112077000135 Crossref Google Scholar
[17] Canuto D. and Taira K., “Two-Dimensional Compressible Viscous Flow Around a Circular Cylinder,” Journal of Fluid Mechanics, Vol. 785, Dec. 2015, pp. 349–371. doi:https://doi.org/10.1017/jfm.2015.635 Crossref Google Scholar
[18] von Helmholtz H., “On Discontinuous Movements of Fluids,” Philosophical Magazine, Vol. 36, No. 244, 1868, pp. 337–346. Google Scholar
[19] Kelvin L., “Hydrokinetic Solutions and Observations,” Philosophical Magazine, Vol. 42, No. 281, 1871, pp. 362–377. Google Scholar
[20] Rosenhead L., “The Formation of Vortices from a Surface of Discontinuity,” Proceedings of the Royal Society of London, Series A: Mathematical, Physical and Engineering Sciences, Vol. 134, No. 823, 1931, pp. 170–192. doi:https://doi.org/10.1098/rspa.1931.0189 Crossref Google Scholar
[21] Winant C. D. and Browand F. K., “Vortex Pairing: The Mechanism of Turbulent Mixing-Layer Growth at Moderate Reynolds Number,” Journal of Fluid Mechanics, Vol. 63, No. 2, 1974, pp. 237–255. doi:https://doi.org/10.1017/S0022112074001121 Crossref Google Scholar
[22] Zaman K. B. M. Q. and Hussain A. K. M. F., “Vortex Pairing in a Circular Jet Under Controlled Excitation. Part 1. General Jet Response,” Journal of Fluid Mechanics, Vol. 101, No. 3, 1980, pp. 449–491. doi:https://doi.org/10.1017/S0022112080001760 Crossref Google Scholar
[23] Melander M. V., Zabusky N. J. and Mcwilliams J. C., “Symmetric Vortex Merger in Two Dimensions: Causes and Conditions,” Journal of Fluid Mechanics, Vol. 195, 1988, pp. 303–340. doi:https://doi.org/10.1017/S0022112088002435 Crossref Google Scholar
[24] Bonnet J. P., Delville J., Glauser M. N., Antonia R. A., Bisset D. K., Cole D. R., Fiedler H. E., Garem J. H., Hilberg D., Jeong J., Kevlahan N. K. R., Ukeiley L. S. and Vincendeaux E., “Collaborative Testing of Eddy Structure Identification Methods in Free Turbulent Shear Flows,” Experiments in Fluids, Vol. 25, No. 3, 1998, pp. 197–225. doi:https://doi.org/10.1007/s003480050224 Crossref Google Scholar
[25] Taira K. and Colonius T., “Three-Dimensional Flows Around Low-Aspect-Ratio Flat-Plate Wings at Low Reynolds Numbers,” Journal of Fluid Mechanics, Vol. 623, 2009, pp. 187–207. doi:https://doi.org/10.1017/S0022112008005314 Crossref Google Scholar
[26] Colonius T. and Taira K., “A Fast Immersed Boundary Method Using a Nullspace Approach and Multi-Domain Far-Field Boundary Conditions,” Computer Methods in Applied Mechanics and Engineering, Vol. 197, Nos. 25–28, 2008, pp. 2131–2146. doi:https://doi.org/10.1016/j.cma.2007.08.014 Crossref Google Scholar
[27] Lumley J. L., “The Structure of Inhomogeneous Turbulent Flows,” Proceedings of the International Colloquium on the Fine Scale Structure of the Atmosphere and Its Influence on Radio Wave Propagation, edited by Yaglam A. M. and Tatarsky V. I., Doklady Akademii Nauk SSSR, Nauka, Moscow, 1967. Google Scholar
[28] Sirovich L., “Turbulence and the Dynamics of Coherent Structures, Parts I–III,” Quarterly of Applied Mathematics, Vol. 45, No. 3, 1987, pp. 561–571. doi:https://doi.org/10.1090/qam/1987-45-03 Crossref Google Scholar
[29] Trefethen L. N. and Bau D., Numerical Linear Algebra, SIAM, Philadelphia, PA, 1997. Crossref Google Scholar
[30] Horn R. A. and Johnson C. R., Matrix Analysis, Cambridge Univ. Press, New York, 1985. Crossref Google Scholar
[31] Golub G. H. and Loan C. F. V., Matrix Computations, 3rd ed., Johns Hopkins Univ. Press, Baltimore, MD, 1996. Google Scholar
[32] Saad Y., Numerical Methods for Large Eigenvalue Problems, 2nd ed., SIAM, Philadelphia, PA, 2011. Crossref Google Scholar
[33] Trefethen L. N. and Embree M., Spectra and Pseudospectra, Princeton Univ. Press, Princeton, NJ, 2005. Crossref Google Scholar
[34] Drazin P. G. and Reid W. H., Hydrodynamic Stability, Cambridge Univ. Press, New York, 1981. Google Scholar
[35] Trefethen L. N., Trefethen A. E., Reddy S. C. and Driscoll T. A., “Hydrodynamic Stability Without Eigenvalues,” Science, Vol. 261, No. 5121, 1993, pp. 578–584. doi:https://doi.org/10.1126/science.261.5121.578 Crossref Google Scholar
[36] Schmid P. J., “Nonmodal Stability Theory,” Annual Review of Fluid Mechanics, Vol. 39, No. 1, 2007, pp. 129–162. doi:https://doi.org/10.1146/annurev.fluid.38.050304.092139 Crossref Google Scholar
[37] Eckart C. and Young G., “The Approximation of One Matrix by Another of Lower Rank,” Psychometrika, Vol. 1, No. 3, 1936, pp. 211–218. doi:https://doi.org/10.1007/BF02288367 Crossref Google Scholar
[38] Anderson E., Bai Z., Bischof C., Blackford L. S., Demmel J., Dongarra J., Croz J. D., Greenbaum A., Hammarling S., McKenney A. and Sorensen D., LAPACK Users’ Guide, SIAM, Philadelphia, PA, 1999. Crossref Google Scholar
[39] Blackford L. S., Choi J., Cleary A., D’Azevedo E., Demmel J., Dhillon I., Dongarra J., Hammarling S., Henry G., Petitet A., Stanley K., Walker D. and Whaley R. C., ScaLAPACK Users’ Guide, SIAM, Philadelphia, PA, 1997. Crossref Google Scholar
[40] Lehoucq R. B., Sorensen D. C. and Yang C., ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, Philadelphia, PA, 1998. Crossref Google Scholar
[41] Trefethen L. N., “Computation of Pseudospectra,” Acta Numerica, Vol. 8, Jan. 1999, pp. 247–295. doi:https://doi.org/10.1017/S0962492900002932 Crossref Google Scholar
[42] Higham D. J. and Trefethen L. N., “Stiffness of ODEs,” BIT Numerical Mathematics, Vol. 33, No. 2, 1993, pp. 285–303. doi:https://doi.org/10.1007/BF01989751 Crossref Google Scholar
[43] Karhunen K., “Zur Spektraltheorie Stochasticher Prozesse,” Annales Academiae Scientiarum Fennicae, Vol. 37, 1946, pp. 1–7. Google Scholar
[44] Loève M., Probability Theory, D. Van Nostrand, Princeton, NJ, 1955. Google Scholar
[45] Berkooz G., Holmes P. and Lumley J. L., “The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows,” Annual Review of Fluid Mechanics, Vol. 25, No. 1, 1993, pp. 539–575. doi:https://doi.org/10.1146/annurev.fl.25.010193.002543 Crossref Google Scholar
[46] Cordier L. and Bergmann M., “Proper Orthogonal Decomposition: An Overview,” Lecture Series 2002-04, 2003-03 and 2008-01 on Post-Processing of Experimental and Numerical Data, Von Karman Institute for Fluid Dynamics 2008, VKI, 2008, p. 46. Google Scholar
[47] Herzog S., “The Large Scale Structure in the Near-Wall Region of Turbulent Pipe Flow,” Ph.D. Thesis, Cornell Univ., Ithaca, NY, 1986. Google Scholar
[48] Aubry N., Holmes P., Lumley J. L. and Stone E., “The Dynamics of Coherent Structures in the Wall Region of a Turbulent Boundary Layer,” Journal of Fluid Mechanics, Vol. 192, July 1988, pp. 115–173. doi:https://doi.org/10.1017/S0022112088001818 Crossref Google Scholar
[49] Moin P. and Kim J., “Numerical Investigation of Turbulent Channel Flow,” Journal of Fluid Mechanics, Vol. 118, May 1982, pp. 341–377. doi:https://doi.org/10.1017/S0022112082001116 Crossref Google Scholar
[50] Glauser M. N. and George W. K., “An Orthogonal Decomposition of the Axisymmetric Jet Mixing Layer Utilizing Cross-Wire Velocity Measurements,” Symposium on Turbulent Shear Flows, Toulouse, 1987. Google Scholar
[51] Ukeiley L., Glauser M. and Wick D., “Downstream Evolution of Proper Orthogonal Decomposition Eigenfunctions in a Lobed Mixer,” AIAA Journal, Vol. 31, No. 8, 1993, pp. 1392–1397. doi:https://doi.org/10.2514/3.11787 Link Google Scholar
[52] Holmes P., Lumley J. L., Berkooz G. and Rowley C., Turbulence, Coherent Structures, Dynamical Systems and Symmetry, 2nd ed., Cambridge Univ. Press, Cambridge, U.K., 2012. Google Scholar
[53] Algazi V. and Sakrison D., “On the Optimality of the Karhunen-Loeve Expansion,” IEEE Transactions on Information Theory, Vol. 15, No. 2, 1969, pp. 319–321. Crossref Google Scholar
[54] Noack B. R., Morzynski M. and Tadmor G. (eds.), Reduced-Order Modelling for Flow Control, Springer, New York, 2011. Crossref Google Scholar
[55] Ukeiley L., Cordier L., Manceau R., Delville J., Glauser M. and Bonnet J. P., “Examination of Large-Scale Structures in a Turbulent Plane Mixing Layer. Part 2. Dynamical Systems Model,” Journal of Fluid Mechanics, Vol. 441, Aug. 2001, pp. 67–108. doi:https://doi.org/10.1017/S0022112001004803 Crossref Google Scholar
[56] Noack B. R., Afanasiev K., Morzynski M., Tadmor G. and Thiele F., “A Hierarchy of Low-Dimensional Models for the Transient and Post-Transient Cylinder Wake,” Journal of Fluid Mechanics, Vol. 497, Dec. 2003, pp. 335–363. doi:https://doi.org/10.1017/S0022112003006694 Crossref Google Scholar
[57] Noack B. R., Papas P. and Monkewitz P. A., “The Need for a Pressure-Term Representation in Empirical Galerkin Models of Incompressible Shear Flows,” Journal of Fluid Mechanics, Vol. 523, Jan. 2005, pp. 339–365. doi:https://doi.org/10.1017/S0022112004002149 Crossref Google Scholar
[58] Rowley C. W., Colonius T. and Murray R. M., “Model Reduction for Compressible Flows Using POD and Galerkin Projection,” Physica D, Vol. 189, Nos. 1–2, 2004, pp. 115–129. doi:https://doi.org/10.1016/j.physd.2003.03.001 Crossref Google Scholar
[59] Nair A. G., Taira K. and Brunton S. L., “Oscillator Network Based Control of Unsteady Fluid Flows,” 20th IFAC World Congress, Paper 4192, Toulouse, July 2017. Google Scholar
[60] Glavaški S., Marsden J. E. and Murray R. M., “Model Reduction, Centering, and the Karhunen–Loeve Expansion,” IEEE Conference on Decision and Control, IEEE Publ., Piscataway, NJ, 1998. doi:https://doi.org/10.1109/CDC.1998.758639 Google Scholar
[61] Noack B. R., “From Snapshots to Model Expansions—Bridging Low Residuals and Pure Frequencies,” Journal of Fluid Mechanics, Vol. 802, Sept. 2016, pp. 1–4. doi:https://doi.org/10.1017/jfm.2016.416 Crossref Google Scholar
[62] Lumley J. L., Stochastic Tools in Turbulence, Academic Press, New York, 1970. Google Scholar
[63] George W. K., “Insight into the Dynamics of Coherent Structures from a Proper Orthogonal Decomposition,” Near Wall Turbulence, edited by Kline S. and Afghan N., Hemisphere, New York, 1988. Google Scholar
[64] Citriniti J. H. and George W. K., “Reconstruction of the Global Velocity Field in the Axisymmetric Mixing Layer Utilizing the Proper Orthogonal Decomposition,” Journal of Fluid Mechanics, Vol. 418, Sept. 2000, pp. 137–166. doi:https://doi.org/10.1017/S0022112000001087 Crossref Google Scholar
[65] Picard C. and Delville J., “Pressure Velocity Coupling in a Subsonic Round Jet,” International Journal of Heat and Fluid Flow, Vol. 21, No. 3, 2000, pp. 359–364. doi:https://doi.org/10.1016/S0142-727X(00)00021-7 Crossref Google Scholar
[66] Suzuki T. and Colonius T., “Instability Waves in a Subsonic Round Jet Detected Using a Near-Field Phased Microphone Array,” Journal of Fluid Mechanics, Vol. 565, No. 1, 2006, pp. 197–226. doi:https://doi.org/10.1017/S0022112006001613 Crossref Google Scholar
[67] Gudmundsson K. and Colonius T., “Instability Wave Models for the Near-Field Fluctuations of Turbulent Jets,” Journal of Fluid Mechanics, Vol. 689, Dec. 2011, pp. 97–128. doi:https://doi.org/10.1017/jfm.2011.401 Crossref Google Scholar
[68] Schmidt O., Towne A., Colonius T., Cavalieri A., Jordan P. and Brès G., “Wavepackets and Trapped Acoustic Modes in a Mach 0.9 Turbulent Jet: A Global Stability Analysis,” Journal of Fluid Mechanics, Vol. 825, Aug. 2017, pp. 1153–1181. doi:https://doi.org/10.1017/jfm.2017.407 Crossref Google Scholar
[69] Towne A., Schmidt O. T. and Colonius T., “Spectral Proper Orthogonal Decomposition and Its Relationship to Dynamic Mode Decomposition and Resolvent Analysis,” Journal of Fluid Mechanics (in review). Google Scholar
[70] Welch P., “The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging over Short, Modified Periodograms,” IEEE Transactions on Audio and Electroacoustics, Vol. 15, No. 2, 1967, pp. 70–73. doi:https://doi.org/10.1109/TAU.1967.1161901 Crossref Google Scholar
[71] Sieber M., Paschereit C. O. and Oberleithner K., “Spectral Proper Orthogonal Decomposition,” Journal of Fluid Mechanics, Vol. 792, April 2016, pp. 798–828. doi:https://doi.org/10.1017/jfm.2016.103 Crossref Google Scholar
[72] Munday P. M. and Taira K., “Wall-Normal Vorticity Injection in Separation Control of NACA 0012 Airfoil,” AIAA Paper 2014-2685, 2014. Google Scholar
[73] Kajishima T. and Taira K., Computational Fluid Dynamics: Incompressible Turbulent Flows, Springer, New York, 2017. Crossref Google Scholar
[74] Munday P. M. and Taira K., “Quantifying Wall-Normal and Angular Momentum Injections in Airfoil Separation Control,” AIAA Journal (in review). Google Scholar
[75] Murray N., Sällström E. and Ukeiley L., “Properties of Subsonic Open Cavity Flow Fields,” Physics of Fluids, Vol. 21, No. 9, 2009, Paper 095103. doi:https://doi.org/10.1063/1.3210772 Crossref Google Scholar
[76] Rowley C. W., “Model Reduction for Fluids Using Balanced Proper Orthogonal Decomposition,” International Journal of Bifurcation and Chaos, Vol. 15, No. 3, 2005, pp. 997–1013. doi:https://doi.org/10.1142/S0218127405012429 Crossref Google Scholar
[77] Camphouse R. C., Myatt J. H., Schmit R. F., Glauser M. N., Ausseur J. M., Andino M. Y. and Wallace R. D., “A Snapshot Decomposition Method for Reduced Order Modeling and Boundary Feedback Control,” AIAA Paper 2008-4195, 2008. Link Google Scholar
[78] Jørgensen B. H., Sørensen J. N. and Brøns M., “Low-Dimensional Modeling of a Driven Cavity Flow with Two Free Parameters,” Theoretical and Computational Fluid Dynamics, Vol. 16, No. 4, 2003, pp. 299–317. doi:https://doi.org/10.1007/s00162-002-0082-9 Crossref Google Scholar
[79] del Sastre P. G. and Bermejo R., “The POD Technique for Computing Bifurcation Diagrams: A Comparison Among Different Models in Fluids,” Numerical Mathematics and Advanced Applications, edited by de Castro A. B., Gómez D., Quintela P. and Salgado P., Springer, Berlin, 2006, pp. 880–888. Crossref Google Scholar
[80] Gordeyev S. and Thomas F. O., “A Temporal Proper Decomposition (TPOD) for Closed-Loop Flow Control,” Experiments in Fluids, Vol. 54, No. 3, 2013, Paper 1477. doi:https://doi.org/10.1007/s00348-013-1477-7 Crossref Google Scholar
[81] Gordeyev S., De Lucca N., Jumper E., Hird K., Juliano T. J., Gregory J. W., Thordahl J. and Wittich D. J., “Comparison of Unsteady Pressure Fields on Turrets with Different Surface Features Using Pressure Sensitive Paint,” Experiments in Fluids, Vol. 55, No. 1, 2014, Paper 1661. doi:https://doi.org/10.1007/s00348-013-1661-9 Crossref Google Scholar
[82] Towne A., Colonius T., Jordan P., Cavalieri A. V. and Brès G. A., “Stochastic and Nonlinear Forcing of Wavepackets in a Mach 0.9 Jet,” 21st AIAA/CEAS Aeroacoustics Conference, AIAA Paper 2015-2217, 2015. Link Google Scholar
[83] Semeraro O., Jaunet V., Jordan P., Cavalieri A. V. G. and Lesshafft L., “Stochastic and Harmonic Optimal Forcing in Subsonic Jets,” 22nd AIAA/CEAS Aeroacoustics Conference, AIAA Paper 2016-2935, 2016. Link Google Scholar
[84] Samimy M., Debiasi M., Caraballo E., Serrani A., Yuan X., Little J. and Myatt J. H., “Feedback Control of Subsonic Cavity Flows Using Reduced-Order Models,” Journal of Fluid Mechanics, Vol. 579, May 2007, pp. 315–346. doi:https://doi.org/10.1017/S0022112007005204 Crossref Google Scholar
[85] Kaiser E., Noack B. R., Cordier L., Spohn A., Segond M., Abel M., Daviller G., Osth J., Krajnović S. and Niven R. K., “Cluster-Based Reduced-Order Modelling of a Mixing Layer,” Journal of Fluid Mechanics, Vol. 754, Sept. 2014, pp. 365–414. doi:https://doi.org/10.1017/jfm.2014.355 Crossref Google Scholar
[86] Moore B. C., “Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction,” IEEE Transactions on Automatic Control, Vol. 26, No. 1, 1981, pp. 17–32. doi:https://doi.org/10.1109/TAC.1981.1102568 Crossref Google Scholar
[87] Juang J.-N. and Pappa R. S., “An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction,” Journal of Guidance, Control, and Dynamics, Vol. 8, No. 5, 1985, pp. 620–627. doi:https://doi.org/10.2514/3.20031 Link Google Scholar
[88] Ma Z., Ahuja S. and Rowley C. W., “Reduced Order Models for Control of Fluids Using the Eigensystem Realization Algorithm,” Theoretical and Computational Fluid Dynamics, Vol. 25, No. 1, 2009, pp. 233–247. Google Scholar
[89] Willcox K. and Peraire J., “Balanced Model Reduction via the Proper Orthogonal Decomposition,” AIAA Journal, Vol. 40, No. 11, 2002, pp. 2323–2330. doi:https://doi.org/10.2514/2.1570 Link Google Scholar
[90] Ilak M. and Rowley C. W., “Modeling of Transitional Channel Flow Using Balanced Proper Orthogonal Decomposition,” Physics of Fluids, Vol. 20, No. 3, 2008, Paper 034103. doi:https://doi.org/10.1063/1.2840197 Crossref Google Scholar
[91] Lall S., Marsden J. E. and Glavaški S., “A Subspace Approach to Balanced Truncation for Model Reduction of Nonlinear Control Systems,” International Journal of Robust and Nonlinear Control, Vol. 12, No. 6, 2002, pp. 519–535. doi:https://doi.org/10.1002/(ISSN)1099-1239 Crossref Google Scholar
[92] Ilak M., Bagheri S., Brandt L., Rowley C. W. and Henningson D. S., “Model Reduction of the Nonlinear Complex Ginzburg-Landau Equation,” SIAM Journal on Applied Dynamical Systems, Vol. 9, No. 4, 2010, pp. 1284–1302. doi:https://doi.org/10.1137/100787350 Crossref Google Scholar
[93] Ahuja S. and Rowley C. W., “Feedback Control of Unstable Steady States of Flow Past a Flat Plate Using Reduced-Order Estimators,” Journal of Fluid Mechanics, Vol. 645, Feb. 2010, pp. 447–478. doi:https://doi.org/10.1017/S0022112009992655 Crossref Google Scholar
[94] Barbagallo A., Sipp D. and Schmid P. J., “Closed-Loop Control of an Open Cavity Flow Using Reduced-Order Models,” Journal of Fluid Mechanics, Vol. 641, Dec. 2009, pp. 1–50. doi:https://doi.org/10.1017/S0022112009991418 Crossref Google Scholar
[95] Bagheri S., Brandt L. and Henningson D. S., “Input-Output Analysis, Model Reduction and Control of the Flat-Plate Boundary Layer,” Journal of Fluid Mechanics, Vol. 620, Feb. 2009, pp. 263–298. doi:https://doi.org/10.1017/S0022112008004394 Crossref Google Scholar
[96] Semeraro O., Bagheri S., Brandt L. and Henningson D. S., “Feedback Control of Three-Dimensional Optimal Disturbances Using Reduced-Order Models,” Journal of Fluid Mechanics, Vol. 677, June 2011, pp. 63–102. doi:https://doi.org/10.1017/S0022112011000620 Crossref Google Scholar
[97] Semeraro O., Bagheri S., Brandt L. and Henningson D. S., “Transition Delay in a Boundary Layer Flow Using Active Control,” Journal of Fluid Mechanics, Vol. 731, Sept. 2013, pp. 288–311. doi:https://doi.org/10.1017/jfm.2013.299 Crossref Google Scholar
[98] Flinois T. L. B., Morgans A. S. and Schmid P. J., “Projection-Free Approximate Balanced Truncation of Large Unstable Systems,” Physical Review E, Vol. 92, No. 2, Aug. 2015, Paper 023012. doi:https://doi.org/10.1103/PhysRevE.92.023012 Crossref Google Scholar
[99] Dergham G., Sipp D., Robinet J.-C. and Barbagallo A., “Model Reduction for Fluids Using Frequential Snapshots,” Physics of Fluids, Vol. 23, No. 6, 2011, Paper 064101. doi:https://doi.org/10.1063/1.3590732 Crossref Google Scholar
[100] Tu J. H. and Rowley C. W., “An Improved Algorithm for Balanced POD Through an Analytic Treatment of Impulse Response Tails,” Journal of Computational Physics, Vol. 231, No. 16, 2012, pp. 5317–5333. doi:https://doi.org/10.1016/j.jcp.2012.04.023 Crossref Google Scholar
[101] Yu D. and Chakravorty S., “A Randomized Proper Orthogonal Decomposition Technique,” 2015 American Control Conference (ACC), IEEE Publ., Piscataway, NJ, 2015, pp. 1137–1142. Crossref Google Scholar
[102] Halko N., Martinsson P.-G. and Tropp J. A., “Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions,” SIAM Review, Vol. 53, No. 2, 2011, pp. 217–288. doi:https://doi.org/10.1137/090771806 Crossref Google Scholar
[103] Schmid P. J. and Sesterhenn J., “Dynamic Mode Decomposition of Numerical and Experimental Data,” 61st Annual Meeting of the APS Division of Fluid Dynamics, American Physical Soc., College Park, MD, 2008. Google Scholar
[104] Schmid P. J., “Dynamic Mode Decomposition of Numerical and Experimental Data,” Journal of Fluid Mechanics, Vol. 656, Aug. 2010, pp. 5–28. doi:https://doi.org/10.1017/S0022112010001217 Crossref Google Scholar
[105] Rowley C. W., Mezić I., Bagheri S., Schlatter P. and Henningson D. S., “Spectral Analysis of Nonlinear Flows,” Journal of Fluid Mechanics, Vol. 641, No. 1, 2009, pp. 115–127. doi:https://doi.org/10.1017/S0022112009992059 Crossref Google Scholar
[106] Mezić I., “Analysis of Fluid Flows via Spectral Properties of the Koopman Operator,” Annual Review of Fluid Mechanics, Vol. 45, No. 1, 2013, pp. 357–378. doi:https://doi.org/10.1146/annurev-fluid-011212-140652 Crossref Google Scholar
[107] Tu J. H., Rowley C. W., Luchtenburg D. M., Brunton S. L. and Kutz J. N., “On Dynamic Mode Decomposition: Theory and Applications,” Journal of Computational Dynamics, Vol. 1, No. 2, 2014, pp. 391–421. doi:https://doi.org/10.3934/jcd Crossref Google Scholar
[108] Chen K. K., Tu J. H. and Rowley C. W., “Variants of Dynamic Mode Decomposition: Boundary Condition, Koopman, and Fourier Analyses,” Journal of Nonlinear Science, Vol. 22, No. 6, 2012, pp. 887–915. doi:https://doi.org/10.1007/s00332-012-9130-9 Crossref Google Scholar
[109] Mann J. and Kutz J. N., “Dynamic Mode Decomposition for Financial Trading Strategies,” Quantitative Finance, Vol. 16, No. 11, 2016, pp. 1643–1655. Crossref Google Scholar
[110] Grosek J. and Kutz J. N., “Dynamic Mode Decomposition for Real-Time Background/Foreground Separation in Video,” 2014. Google Scholar
[111] Kutz J. N., Grosek J. and Brunton S. L., “Dynamic Mode Decomposition for Robust PCA with Applications to Foreground/Background Subtraction in Video Streams and Multi-Resolution Analysis,” CRC Handbook on Robust Low-Rank and Sparse Matrix Decomposition: Applications in Image and Video Processing, edited by Bouwmans T., CRC Press, Boca Raton, FL, 2015. Google Scholar
[112] Erichson N. B. and Donovan C., “Randomized Low-Rank Dynamic Mode Decomposition for Motion Detection,” Computer Vision and Image Understanding, Vol. 146, May 2016, pp. 40–50. doi:https://doi.org/10.1016/j.cviu.2016.02.005 Crossref Google Scholar
[113] Proctor J. L. and Eckhoff P. A., “Discovering Dynamic Patterns from Infectious Disease Data Using Dynamic Mode Decomposition,” International Health, Vol. 7, No. 2, 2015, pp. 139–145. doi:https://doi.org/10.1093/inthealth/ihv009 Crossref Google Scholar
[114] Berger E., Sastuba M., Vogt D., Jung B. and Amor H. B., “Estimation of Perturbations in Robotic Behavior Using Dynamic Mode Decomposition,” Journal of Advanced Robotics, Vol. 29, No. 5, 2015, pp. 331–343. doi:https://doi.org/10.1080/01691864.2014.981292 Crossref Google Scholar
[115] Brunton B. W., Johnson L. A., Ojemann J. G. and Kutz J. N., “Extracting Spatial–Temporal Coherent Patterns in Large-Scale Neural Recordings Using Dynamic Mode Decomposition,” Journal of Neuroscience Methods, Vol. 258, Jan. 2016, pp. 1–15. doi:https://doi.org/10.1016/j.jneumeth.2015.10.010 Crossref Google Scholar
[116] Hemati M. S., Williams M. O. and Rowley C. W., “Dynamic Mode Decomposition for Large and Streaming Datasets,” Physics of Fluids, Vol. 26, No. 11, 2014, Paper 111701. doi:https://doi.org/10.1063/1.4901016 Crossref Google Scholar
[117] Belson B. A., Tu J. H. and Rowley C. W., “A Parallelized Model Reduction Library,” ACM Transactions on Mathematical Software, Vol. 40, No. 4, June 2014, Paper 30. Crossref Google Scholar
[118] Wynn A., Pearson D. S., Ganapathisubramani B. and Goulart P. J., “Optimal Mode Decomposition for Unsteady Flows,” Journal of Fluid Mechanics, Vol. 733, Oct. 2013, pp. 473–503. doi:https://doi.org/10.1017/jfm.2013.426 Crossref Google Scholar
[119] Jovanović M. R., Schmid P. J. and Nichols J. W., “Sparsity-Promoting Dynamic Mode Decomposition,” Physics of Fluids, Vol. 26, No. 2, 2014, Paper 024103. Crossref Google Scholar
[120] Tu J. H., Rowley C. W., Kutz J. N. and Shang J. K., “Spectral Analysis of Fluid Flows Using Sub-Nyquist-Rate PIV Data,” Experiments in Fluids, Vol. 55, No. 9, 2014, pp. 1–13. doi:https://doi.org/10.1007/s00348-014-1805-6 Crossref Google Scholar
[121] Brunton S. L., Proctor J. L., Tu J. H. and Kutz J. N., “Compressed Sensing and Dynamic Mode Decomposition,” Journal of Computational Dynamics, Vol. 2, No. 2, 2015, pp. 165–191. Google Scholar
[122] Gueniat F., Mathelin L. and Pastur L., “A Dynamic Mode Decomposition Approach for Large and Arbitrarily Sampled Systems,” Physics of Fluids, Vol. 27, No. 2, 2015, Paper 025113. doi:https://doi.org/10.1063/1.4908073 Crossref Google Scholar
[123] Erichson N. B., Brunton S. L. and Kutz J. N., “Compressed Dynamic Mode Decomposition for Real-Time Object Detection,” Journal of Real-Time Image Processing, Nov. 2016. doi:https://doi.org/10.1007/s11554-016-0655-2. Crossref Google Scholar
[124] Proctor J. L., Brunton S. L. and Kutz J. N., “Dynamic Mode Decomposition with Control,” SIAM Journal on Applied Dynamical Systems, Vol. 15, No. 1, 2016, pp. 142–161. doi:https://doi.org/10.1137/15M1013857 Crossref Google Scholar
[125] Williams M. O., Kevrekidis I. G. and Rowley C. W., “A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition,” Journal of Nonlinear Science, Vol. 25, No. 6, 2015, pp. 1307–1346. Crossref Google Scholar
[126] Williams M. O., Rowley C. W., Mezić I. and Kevrekidis I. G., “Data Fusion via Intrinsic Dynamic Variables: An Application of Data-Driven Koopman Spectral Analysis,” EPL (Europhysics Letters), Vol. 109, No. 4, 2015, Paper 40007. doi:https://doi.org/10.1209/0295-5075/109/40007 Crossref Google Scholar
[127] Duke D., Soria J. and Honnery D., “An Error Analysis of the Dynamic Mode Decomposition,” Experiments in Fluids, Vol. 52, No. 2, 2012, pp. 529–542. doi:https://doi.org/10.1007/s00348-011-1235-7 Crossref Google Scholar
[128] Dawson S. T. M., Hemati M. S., Williams M. O. and Rowley C. W., “Characterizing and Correcting for the Effect of Sensor Noise in the Dynamic Mode Decomposition,” Experiments in Fluids, Vol. 57, No. 3, 2016, pp. 1–19. doi:https://doi.org/10.1007/s00348-016-2127-7 Crossref Google Scholar
[129] Bagheri S., “Effects of Weak Noise on Oscillating Flows: Linking Quality Factor, Floquet Modes, and Koopman Spectrum,” Physics of Fluids, Vol. 26, No. 9, 2014, Paper 094104. doi:https://doi.org/10.1063/1.4895898 Crossref Google Scholar
[130] Hemati M. S., Rowley C. W., Deem E. A. and Cattafesta L. N., “De-Biasing the Dynamic Mode Decomposition for Applied Koopman Spectral Analysis of Noisy Datasets,” Theoretical and Computational Fluid Dynamics, Vol. 31, No. 4, Aug. 2017, pp. 349–368. doi:https://doi.org/10.1007/s00162-017-0432-2. Crossref Google Scholar
[131] Noack B. R., Stankiewicz W., Morzyński M. and Schmid P. J., “Recursive Dynamic Mode Decomposition of a Transient and Post-Transient Cylinder Wake Flows,” Journal of Fluid Mechanics, Vol. 809, Dec. 2016, pp. 843–872. doi:https://doi.org/10.1017/jfm.2016.678 Crossref Google Scholar
[132] Kutz J. N., Fu X. and Brunton S. L., “Multi-Resolution Dynamic Mode Decomposition,” SIAM Journal on Applied Dynamical Systems, Vol. 15, No. 2, 2016, pp. 713–735. doi:https://doi.org/10.1137/15M1023543 Crossref Google Scholar
[133] Brunton S. L., Brunton B. W., Proctor J. L., Kaiser E. and Kutz J. N., “Chaos as an Intermittently Forced Linear System,” Nature Communications, Vol. 8, May 2017, Paper 19. Crossref Google Scholar
[134] Tu J. H., “Dynamic Mode Decomposition: Theory and Applications,” Ph.D. Thesis, Princeton Univ., Princeton, NJ, 2013. Google Scholar
[135] Holmes P. and Guckenheimer J., Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Vol. 42, Applied Mathematical Sciences, Springer–Verlag, Berlin, 1983. Google Scholar
[136] Koopman B. O., “Hamiltonian Systems and Transformation in Hilbert Space,” Proceedings of the National Academy of Sciences, Vol. 17, No. 5, 1931, pp. 315–318. doi:https://doi.org/10.1073/pnas.17.5.315 Crossref Google Scholar
[137] Mezić I., “Spectral Properties of Dynamical Systems, Model Reduction and Decompositions,” Nonlinear Dynamics, Vol. 41, Nos. 1–3, 2005, pp. 309–325. Crossref Google Scholar
[138] Budišić M., Mohr R. and Mezić I., “Applied Koopmanism,” Chaos: An Interdisciplinary Journal of Nonlinear Science, Vol. 22, No. 4, 2012, Paper 047510. Crossref Google Scholar
[139] Kevrekidis I. G., Gear C. W., Hyman J. M., Kevrekidis P. G., Runborg O. and Theodoropoulos C., “Equation-Free, Coarse-Grained Multiscale Computation: Enabling Microscopic Simulators to Perform System-Level Analysis,” Communications in Mathematical Science, Vol. 1, No. 4, 2003, pp. 715–762. doi:https://doi.org/10.4310/CMS.2003.v1.n4.a5 Crossref Google Scholar
[140] Giannakis D., “Data-Driven Spectral Decomposition and Forecasting of Ergodic Dynamical Systems,” 2015. Google Scholar
[141] Kaiser E., Kutz J. N. and Brunton S. L., “Data-Driven Discovery of Koopman Eigenfunctions for Control,” Automatica, 2017 (in review). Google Scholar
[142] Brunton S. L., Brunton B. W., Proctor J. L. and Kutz J. N., “Koopman Observable Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control,” PLoS ONE, Vol. 11, No. 2, 2016, Paper e0150171. doi:https://doi.org/10.1371/journal.pone.0150171 Crossref Google Scholar
[143] Bagheri S., “Koopman-Mode Decomposition of the Cylinder Wake,” Journal of Fluid Mechanics, Vol. 726, July 2013, pp. 596–623. doi:https://doi.org/10.1017/jfm.2013.249 Crossref Google Scholar
[144] Gaspard P., Nicolis G., Provata A. and Tasaki S., “Spectral Signature of the Pitchfork Bifurcation: Liouville Equation Approach,” Physical Review E, Vol. 51, No. 1, 1995, pp. 74–94. doi:https://doi.org/10.1103/PhysRevE.51.74 Crossref Google Scholar
[145] Susuki Y. and Mezic I., “Nonlinear Koopman Modes and Coherency Identification of Coupled Swing Dynamics,” IEEE Transactions on Power Systems, Vol. 26, No. 4, 2011, pp. 1894–1904. doi:https://doi.org/10.1109/TPWRS.2010.2103369 Crossref Google Scholar
[146] Susuki Y. and Mezic I., “Nonlinear Koopman Modes and a Precursor to Power System Swing Instabilities,” IEEE Transactions on Power Systems, Vol. 27, No. 3, 2012, pp. 1182–1191. doi:https://doi.org/10.1109/TPWRS.2012.2183625 Crossref Google Scholar
[147] Brunton S. L., Proctor J. L. and Kutz J. N., “Discovering Governing Equations from Data by Sparse Identification of Nonlinear Dynamical Systems,” Proceedings of the National Academy of Sciences, Vol. 113, No. 15, 2016, pp. 3932–3937. doi:https://doi.org/10.1073/pnas.1517384113 Crossref Google Scholar
[148] Arbabi H. and Mezić I., “Ergodic Theory, Dynamic Mode Decomposition and Computation of Spectral Properties of the Koopman Operator,” 2016. Google Scholar
[149] Juniper M. P., Hanifi A. and Theofilis V., “Modal Stability Theory,” Applied Mechanics Review, Vol. 66, No. 2, 2014, Paper 024804. doi:https://doi.org/10.1115/1.4026604 Google Scholar
[150] Theofilis V., “Advances in Global Linear Instability of Nonparallel and Three-Dimensional Flows,” Progress in Aerospace Sciences, Vol. 39, No. 4, 2003, pp. 249–315. doi:https://doi.org/10.1016/S0376-0421(02)00030-1 Crossref Google Scholar
[151] Åkervik E., Ehrenstein U., Gallaire F. and Henningson D. S., “Global Two-Dimensional Stability Measures of the Flat Plate Boundary-Layer Flow,” European Journal of Mechanics-B/Fluids, Vol. 27, No. 5, 2008, pp. 501–513. doi:https://doi.org/10.1016/j.euromechflu.2007.09.004 Crossref Google Scholar
[152] Theofilis V. and Colonius T., “An Algorithm for the Recovery of 2- and 3-D BiGlobal Instabilities of Compressible Flow over 2-D Open Cavities,” AIAA Paper 2003-4143, 2003. Link Google Scholar
[153] Zhang W. and Samtaney R., “BiGlobal Linear Stability Analysis on Low-Re Flow Past an Airfoil at High Angle of Attack,” Physics of Fluids, Vol. 28, No. 4, 2016, Paper 044105. Crossref Google Scholar
[154] Sun Y., Taira K., Cattafesta L. N. and Ukeiley L. S., “Spanwise Effects on Instabilities of Compressible Flow over a Long Rectangular Cavity,” Theoretical and Computational Fluid Dynamics, Nov. 2016, pp. 1–11. doi:https://doi.org/10.1007/s00162-016-0412-y Google Scholar
[155] Sun Y., Taira K., Cattafesta L. N. and Ukeiley L. S., “Biglobal Instabilities of Spanwise-Periodic Compressible Open-Cavity Flows,” Journal of Fluid Mechanics, Vol. 826, Sept. 2017, pp. 270–301. doi:https://doi.org/10.1017/jfm.2017.416 Crossref Google Scholar
[156] Rodríguez D. and Theofilis V., “Massively Parallel Solution of the BiGlobal Eigenvalue Problem Using Dense Linear Algebra,” AIAA Journal, Vol. 47, No. 10, 2009, pp. 2449–2459. doi:https://doi.org/10.2514/1.42714 Link Google Scholar
[157] Tezuka A. and Suzuki K., “Three-Dimensional Global Linear Stability Analysis of Flow Around a Spheroid,” AIAA Journal, Vol. 44, No. 8, 2006, pp. 1697–1708. doi:https://doi.org/10.2514/1.16632 Link Google Scholar
[158] Bagheri S., Schlatter P., Schmid P. J. and Henningson D. S., “Global Stability of a Jet in Crossflow,” Journal of Fluid Mechanics, Vol. 624, April 2009, pp. 33–44. doi:https://doi.org/10.1017/S0022112009006053 Crossref Google Scholar
[159] Liu Q., Gómez F. and Theofilis V., “Linear Instability Analysis of Low-Re Incompressible Flow over a Long Rectangular Finite-Span Open Cavity,” Journal of Fluid Mechanics, Vol. 799, July 2016, Paper R2. doi:https://doi.org/10.1017/jfm.2016.391 Crossref Google Scholar
[160] Keller H. B., “Approximation Methods for Nonlinear Problems with Application to Two-Point Boundary Value Problems,” Mathematics of Computation, Vol. 29, No. 130, 1975, pp. 464–464. doi:https://doi.org/10.1090/S0025-5718-1975-0371058-7 Crossref Google Scholar
[161] Bres G. A. and Colonius T., “Three-Dimensional Instabilities in Compressible Flow over Open Cavities,” Journal of Fluid Mechanics, Vol. 599, March 2008, pp. 309–339. doi:https://doi.org/10.1017/S0022112007009925 Crossref Google Scholar
[162] Paredes P., Hermanns M., Chainche S. L. and Theofilis V., “Order $10^{4}$ Speedup in Global Linear Instability Analysis Using Matrix Formation,” Computer Methods in Applied Mechanics and Engineering, Vol. 253, Jan. 2013, pp. 287–304. doi:https://doi.org/10.1016/j.cma.2012.09.014 Crossref Google Scholar
[163] Tatsumi T. and Yoshimura T., “Stability of the Laminar Flow in a Rectangular Duct,” Journal of Fluid Mechanics, Vol. 212, March 1990, pp. 437–449. doi:https://doi.org/10.1017/S002211209000204X Crossref Google Scholar
[164] Paredes P. and Theofilis V., “Centerline Instabilities on the Hypersonic International Flight Research Experimentation HIFiRE-5 Elliptic Cone Model,” Journal of Fluids and Structures, Vol. 53, Feb. 2015, pp. 36–49. doi:https://doi.org/10.1016/j.jfluidstructs.2014.11.002 Crossref Google Scholar
[165] Åkervik E., Brandt L., Henningson D. S., Hœpffner J., Marxen O. and Schlatter P., “Steady Solutions of the Navier-Stokes Equations by Selective Frequency Damping,” Physics of Fluids, Vol. 18, No. 6, 2006, Paper 068102. doi:https://doi.org/10.1063/1.2211705 Crossref Google Scholar
[166] Tuckerman L. and Barkley D., “Bifurcation Analysis for Timesteppers,” Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems, Vol. 19, edited by Doedel E. and Tuckerman L. S., Springer, New York, 2000, pp. 453–466. Crossref Google Scholar
[167] Kelley C. T., Kevrekidis G. I. and Qiao L., “Newton–Krylov Solvers for Timesteppers,” 2004. Google Scholar
[168] Barkley D. and Henderson R. D., “Three-Dimensional Floquet Stability Analysis of the Wake of a Circular Cylinder,” Journal of Fluid Mechanics, Vol. 322, Sept. 1996, pp. 215–241. doi:https://doi.org/10.1017/S0022112096002777 Crossref Google Scholar
[169] He W., Gioria R. S., Pérez J. M. and Theofilis V., “Linear Instability of Low Reynolds Number Massively Separated Flow Around Three NACA Airfoils,” Journal of Fluid Mechanics, Vol. 811, Jan. 2017, pp. 701–741. doi:https://doi.org/10.1017/jfm.2016.778 Crossref Google Scholar
[170] Edstrand A. M., Davis T. B., Schmid P. J., Taira K. and Cattafesta L. N., “On the Mechanism of Trailing Vortex Wandering,” Journal of Fluid Mechanics, Vol. 801, Aug. 2016, Paper R1. doi:https://doi.org/10.1017/jfm.2016.440 Crossref Google Scholar
[171] Sipp D. and Lebedev A., “Global Stability of Base and Mean Flows: A General Approach and Its Applications to Cylinder and Open Cavity Flows,” Journal of Fluid Mechanics, Vol. 593, Dec. 2007, pp. 333–358. doi:https://doi.org/10.1017/S0022112007008907 Crossref Google Scholar
[172] Luchini P., “Reynolds-Number-Independent Instability of the Boundary Layer over a Flat Surface: Optimal Perturbations,” Journal of Fluid Mechanics, Vol. 404, Feb. 2000, pp. 289–309. doi:https://doi.org/10.1017/S0022112099007259 Crossref Google Scholar
[173] Abdessemed N., Sharma A. S., Sherwin S. J. and Theofilis V., “Transient Growth Analysis of the Flow Past a Circular,” Physics of Fluids, Vol. 21, No. 4, 2009, Paper 044103. doi:https://doi.org/10.1063/1.3112738 Crossref Google Scholar
[174] Abdessemed N., Sherwin S. J. and Theofilis V., “Linear Instability Analysis of Low Pressure Turbine Flows,” Journal of Fluid Mechanics, Vol. 628, June 2009, pp. 57–83. doi:https://doi.org/10.1017/S0022112009006272 Crossref Google Scholar
[175] Sharma A., Abdessemed N., Sherwin S. J. and Theofilis V., “Transient Growth Mechanisms of Low Reynolds Number Flow over a Low-Pressure Turbine Blade,” Theoretical and Computational Fluid Dynamics, Vol. 25, Nos. 1–4, 2011, pp. 19–30. doi:https://doi.org/10.1007/s00162-010-0183-9 Crossref Google Scholar
[176] Gaster M., Kit E. and Wygnanski I., “Large-Scale Structures in a Forced Turbulent Mixing Layer,” Journal of Fluid Mechanics, Vol. 150, Jan. 1985, pp. 23–39. doi:https://doi.org/10.1017/S0022112085000027 Crossref Google Scholar
[177] Paredes P., Terhaar S., Oberleithner K., Theofilis V. and Paschereit C. O., “Global and Local Hydrodynamic Stability Analysis as a Tool for Combustor Dynamics Modeling,” Journal of Engineering for Gas Turbines and Power, Vol. 138, No. 2, 2016, Paper 021504. Crossref Google Scholar
[178] Theofilis V., Hein S. and Dallmann U., “On the Origins of Unsteadiness and Three-Dimensionality in a Laminar Separation Bubble,” Philosophical Transactions of the Royal Society of London, Series A: Mathematical, Physical, and Engineering Sciences, Vol. 358, No. 1777, 2000, pp. 3229–3246. doi:https://doi.org/10.1098/rsta.2000.0706 Crossref Google Scholar
[179] Rodríguez D. and Theofilis V., “On the Birth of Stall Cells on Airfoils,” Theoretical and Computational Fluid Dynamics, Vol. 25, Nos. 1–4, 2011, pp. 105–117. doi:https://doi.org/10.1007/s00162-010-0193-7 Crossref Google Scholar
[180] Gaitonde D. V., Canupp P. W. and Holden M. S., “Heat Transfer Predictions in a Laminar Hypersonic Viscous/Inviscid Interaction,” Journal of Thermophysics and Heat Transfer, Vol. 16, No. 4, 2002, pp. 481–489. doi:https://doi.org/10.2514/2.6714 Link Google Scholar
[181] Tumuklu O., Li Z. and Levin D. A., “Particle Ellipsoidal Statistical Bhatnagar–Gross–Krook Approach for Simulation of Hypersonic Shocks,” AIAA Journal, Vol. 54, No. 12, 2016, pp. 3701–3716. doi:https://doi.org/10.2514/1.J054837 Link Google Scholar
[182] Tumuklu O., Levin D. A. and Theofilis V., “Assessment of Degree of Steadiness in Boundary Layer Shock-Interaction Flows,” AIAA Paper 2017-1614, 2017. Google Scholar
[183] Tumuklu O., Pérez J., Levin D. and Theofilis V., “On Linear Stability Analyses of Hypersonic Laminar Separated Flows in a DSMC Framework. Part I: Base Flow Computations in a Double Cone and a ‘Tick’ Model,” 57th IACAS Conference, Tel Aviv/Haifa, Israel, March 2017, Paper ThL1T3.1. Google Scholar
[184] Tumuklu O., Pérez J., Theofilis V. and Levin D., “On Linear Stability Analyses of Hypersonic Laminar Separated Flows in a DSMC Framework. Part II: Residuals Algorithm and the Least Damped Global Modes,” 57th IACAS Conference, Tel Aviv/Haifa, Israel, March 2017, Paper ThL1T3.5. Google Scholar
[185] Paredes P., Gosse R., Theofilis V. and Kimmel R., “Linear Modal Instabilities of Hypersonic Flow over an Elliptic Cone,” Journal of Fluid Mechanics, Vol. 804, Oct. 2016, pp. 442–466. doi:https://doi.org/10.1017/jfm.2016.536 Crossref Google Scholar
[186] Herbert T., “Parabolized Stability Equations,” Annual Review of Fluid Mechanics, Vol. 29, No. 1, 1997, pp. 245–283. doi:https://doi.org/10.1146/annurev.fluid.29.1.245 Crossref Google Scholar
[187] Broadhurst M. S., Theofilis V. and Sherwin S. J., “Spectral Element Stability Analysis of Vortical Flows,” IUTAM Symposium on Laminar-Turbulent Transition, edited by Govindarajan R., Vol. 78, Fluid Mechanics and Its Applications, Springer, Dordrecht, 2006, pp. 153–158. Crossref Google Scholar
[188] Broadhurst M. S. and Sherwin S. J., “The Parabolised Stability Equations for 3D-Flows: Implementation and Numerical Stability,” Applied Numerical Mathematics, Vol. 58, No. 7, 2008, pp. 1017–1029. doi:https://doi.org/10.1016/j.apnum.2007.04.016 Crossref Google Scholar
[189] Paredes P., Theofilis V., Rodrguez D. and Tendero J. A., “The PSE-3D Instability Analysis Methodology for Flows Depending Strongly on Two and Weakly on the Third Spatial Dimension,” AIAA Paper 2011-3752, 2011. Link Google Scholar
[190] De Tullio N., Paredes P., Sandham N. and Theofilis V., “Roughness-Induced Instability and Breakdown to Turbulence in a Supersonic Boundary-Layer,” Journal of Fluid Mechanics, Vol. 735, Nov. 2013, pp. 613–646. doi:https://doi.org/10.1017/jfm.2013.520 Crossref Google Scholar
[191] Martín J. and Paredes P., “Three-Dimensional Instability Analysis of Boundary Layers Perturbed by Streamwise Vortices,” Theoretical and Computational Fluid Dynamics, Aug. 2016, pp. 1–13. doi:https://doi.org/10.1007/s00162-016-0403-z Google Scholar
[192] Galionis I. and Hall P., “Spatial Stability of the Incompressible Corner Flow,” Theoretical and Computational Fluid Dynamics, Vol. 19, No. 2, 2005, pp. 77–113. doi:https://doi.org/10.1007/s00162-004-0153-1 Crossref Google Scholar
[193] Galionis I. and Hall P., “Stability of the Flow in a Slowly Diverging Rectangular Duct,” Journal of Fluid Mechanics, Vol. 555, May 2006, pp. 43–58. doi:https://doi.org/10.1017/S0022112006008664 Crossref Google Scholar
[194] Luchini P. and Bottaro A., “Adjoint Equations in Stability Analysis,” Annual Review of Fluid Mechanics, Vol. 46, No. 1, 2014, pp. 493–517. doi:https://doi.org/10.1146/annurev-fluid-010313-141253 Crossref Google Scholar
[195] McKeon B. J., Sharma A. S. and Jacobi I., “Experimental Manipulation of Wall Turbulence: A Systems Approach,” Physics of Fluids, Vol. 25, No. 3, 2013, Paper 031301. Crossref Google Scholar
[196] Jeun J., Nichols J. W. and Jovanović M. R., “Input-Output Analysis of High-Speed Axisymmetric Isothermal Jet Noise,” Physics of Fluids, Vol. 28, No. 4, 2016, Paper 047101. doi:https://doi.org/10.1063/1.4946886 Crossref Google Scholar
[197] Garnaut X., Lesshafft L., Schmid P. and Huerre P., “The Preferred Mode of Incompressible Jets: Linear Frequency Response Analysis,” Journal of Fluid Mechanics, Vol. 716, Feb. 2013, pp. 189–202. doi:https://doi.org/10.1017/jfm.2012.540 Crossref Google Scholar
[198] Luhar M., Sharma A. S. and McKeon B. J., Repository for all Matlab/Python Code for the Navier-Stokes Resolvent Analysis [online database], 2013, https://github.com/mluhar/resolvent [retrieved 10 Feb. 2017]. Google Scholar
[199] McKeon B. J. and Sharma A. S., “A Critical Layer Model for Turbulent Pipe Flow,” Journal of Fluid Mechanics, Vol. 658, Sept. 2010, pp. 336–382. doi:https://doi.org/10.1017/S002211201000176X Crossref Google Scholar
[200] Farrell B. and Ioannou J., “Stochastic Forcing of the Linearized Navier-Stokes Equations,” Physics of Fluids, Vol. 5, No. 11, 1993, pp. 2600–2609. doi:https://doi.org/10.1063/1.858894 Crossref Google Scholar
[201] Jovanović M. R. and Bamieh B., “Componentwise Energy Amplification in Channel Flows,” Journal of Fluid Mechanics, Vol. 534, July 2005, pp. 145–183. doi:https://doi.org/10.1017/S0022112005004295 Crossref Google Scholar
[202] Towne A., Colonius T. and Schmidt O., “Empirical Resolvent Mode Decomposition,” Bulletin of the American Physical Society, No. G17.00006, 2015. Google Scholar
[203] Beneddine S., Yegavian R., Sipp D. and Leclaire B., “Unsteady Flow Dynamics Reconstruction from Mean Flow and Point Sensors: An Experimental Study,” Journal of Fluid Mechanics, Vol. 824, Aug. 2017, pp. 174–201. doi:https://doi.org/10.1017/jfm.2017.333 Crossref Google Scholar
[204] Beneddine S., Sipp D., Arnault A., Dandois J. and Lesshafft L., “Conditions for Validity of Mean Flow Stability Analysis and Application to the Determination of Coherent Structures in a Turbulent Backward Facing Step Flow,” Journal of Fluid Mechanics, Vol. 798, July 2016, pp. 485–504. doi:https://doi.org/10.1017/jfm.2016.331 Crossref Google Scholar
[205] Gómez F., Blackburn H. M., Rudman M., Sharma A. S. and McKeon B. J., “A Reduced-Order Model of Three-Dimensional Unsteady Flow in a Cavity Based on the Resolvent Operator,” Journal of Fluid Mechanics, Vol. 798, July 2016, Paper R2. doi:https://doi.org/10.1017/jfm.2016.339 Crossref Google Scholar
[206] Lu L. and Papadakis G., “An Iterative Method for the Computation of the Response of Linearised Navier-Stokes Equations to Harmonic Forcing and Application to Forced Cylinder Wakes,” International Journal for Numerical Methods in Fluids, Vol. 74, No. 11, 2014, pp. 794–817. doi:https://doi.org/10.1002/fld.v74.11 Crossref Google Scholar
[207] Luhar M., Sharma A. S. and McKeon B. J., “A Framework for Studying the Effect of Compliant Surfaces on Wall Turbulence,” Journal of Fluid Mechanics, Vol. 768, April 2015, pp. 415–441. doi:https://doi.org/10.1017/jfm.2015.85 Crossref Google Scholar
[208] Jacobi I. and McKeon B. J., “Dynamic Roughness-Perturbation of a Turbulent Boundary Layer,” Journal of Fluid Mechanics, Vol. 688, Dec. 2011, pp. 258–296. doi:https://doi.org/10.1017/jfm.2011.375 Crossref Google Scholar
[209] Luhar M., Sharma A. S. and McKeon B. J., “Opposition Control within the Resolvent Analysis Framework,” Journal of Fluid Mechanics, Vol. 749, June 2014, pp. 597–626. doi:https://doi.org/10.1017/jfm.2014.209 Crossref Google Scholar
[210] Sharma A. S. and McKeon B. J., “On Coherent Structure in Wall Turbulence,” Journal of Fluid Mechanics, Vol. 728, Aug. 2013, pp. 196–238. doi:https://doi.org/10.1017/jfm.2013.286 Crossref Google Scholar
[211] Towne A., Colonius T., Jordan P., Cavalieri A. V. and Bres G. A., “Stochastic and Nonlinear Forcing of Wavepackets in a Mach 0.9 Jet,” AIAA Paper 2015-2217, 2015. Link Google Scholar
[212] Moarref R., Sharma A. S., Tropp J. A. and McKeon B. J., “Model-Based Scaling and Prediction of the Streamwise Energy Intensity in High-Reynolds Number Turbulent Channels,” Journal of Fluid Mechanics, Vol. 734, Nov. 2013, pp. 275–316. doi:https://doi.org/10.1017/jfm.2013.457 Crossref Google Scholar
[213] Moarref R., Jovanovic M. R., Sharma A. S., Tropp J. A. and McKeon B. J., “A Low-Order Decomposition of Turbulent Channel Flow via Resolvent Analysis and Convex Optimization,” Physics of Fluids, Vol. 26, No. 5, 2014, Paper 051701. Crossref Google Scholar
[214] Sharma A. S., Moarref R., McKeon B. J., Park J. S., Graham M. and Willis A. P., “Low-Dimensional Representations of Exact Coherent States of the Navier-Stokes Equations from the Resolvent Model of Wall Turbulence,” Physical Review E, Vol. 93, 2016, Paper 021102(R). Crossref Google Scholar
[215] Sharma A. S., Moarref R. and McKeon B. J., “Scaling and Interaction of Self-Similar Modes in Models of High Reynolds Number Wall Turbulence,” Philosophical Transactions of the Royal Society of London, Series A: Mathematical and Physical Sciences, Vol. 375, No. 2089, March 2017, Paper 20160089. doi:https://doi.org/10.1098/rsta.2016.0089 Crossref Google Scholar
[216] Sharma A. S., Mezić I. and McKeon B. J., “On the Correspondence Between Koopman Mode Decomposition, Resolvent Mode Decomposition, and Invariant Solutions of the Navier-Stokes Equations,” Physical Review Fluids, Vol. 1, No. 3, 2016, Paper 032403(R). doi:https://doi.org/10.1103/PhysRevFluids.1.032402 Crossref Google Scholar
[217] Brunton S. L. and Noack B. R., “Closed-Loop Turbulence Control: Progress and Challenges,” Applied Mechanics Reviews, Vol. 67, No. 5, 2015, Paper 050801. doi:https://doi.org/10.1115/1.4031175 Crossref Google Scholar
[218] Bui-Thanh T., Damodaran M. and Wilcox K., “Aerodynamic Data Reconstruction and Inverse Design Using Proper Orthogonal Decomposition,” AIAA Journal, Vol. 42, No. 8, 2004, pp. 1505–1516. doi:https://doi.org/10.2514/1.2159 Link Google Scholar
[219] LeGresley P. A., “Application of Proper Orthogonal Decomposition (POD) to Design Decomposition Methods,” Ph.D. Thesis, Stanford Univ., Stanford, CA, 2005. Google Scholar

Tables

**Table 1 Summary of the modal decomposition/analysis techniques for fluid flows presented in the present papera**
Techniques	Sections	Inputs	General descriptions
Data-based
POD	III	Data (L or NL flow; C and E)	Determines the optimal set of modes to represent data based on $L_{2}$ norm (energy)
Balanced POD	IV	Data (L forward and L adjoint flow; C)	Gives balancing and adjoint modes based on input–output relation (balanced truncation)
DMD	V	Data (L or NL flow; C and E)	Captures dynamic modes with associated growth rates and frequencies; linear approximation to nonlinear dynamics
Operator-based
Koopman analysis	VI	Theoretical (also see DMD)	Transforms nonlinear dynamics into linear representation but with an infinite-dimensional operator; Koopman modes are approximated by DMD modes
Global linear stability analysis	VII	L NS operators and base flow (C)	Finds linear stability modes about a base flow (i.e., steady state); assumes small perturbations about base flow
Resolvent analysis	VIII	L NS operators and base flow (C)	Provides forcing and response modes based on input–output analysis with respect to a base flow (including time-averaged mean flow); can be applied to turbulent flow

aL (linear), NL (nonlinear), C (computational), E (experimental), and NS (Navier–Stokes).

**Table 2 MATLAB eig and svd examples for eigenvalue and singular value decompositions**
Decomposition	MATLAB code
Eigenvalue decomposition	[V,Lambda] = eig(A)	(AV = VLambda)	Eq. (2)
Generalized eigenvalue decomposition	[V,Lambda] = eig(A,B)	(AV = BV*Lambda)	Eq. (7)
Singular value decomposition	[U,Sigma,V] = svd(A)	(A = USigmaV’)	Eq. (11)
Reduced (economy-sized) SVD	[U,Sigma,V] = svd(A,’econ’)	(A = USigmaV’)	Eq. (12)