Welcome back. This week, I have prepared a topic that will hopefully be both interesting and exciting. We will cover a computational technique called *resolvent analysis* that allows us to understand how to best control a fluid flow so that it acts in a desirable way. The specific problem can be outlined as follows: For the last eight decades, most aircraft have flown at high Reynolds numbers, typically on the order of 10^{6} or more. In these flow settings, boundary layers are thin and we can model most of the flow around the wing as inviscid. These inviscid approximations have allowed aerodynamicists to develop relatively elegant and well-understood descriptions of the lift, drag, and moments about an airplane. However, in the 2020s, there is an increasing need to understand lower Reynolds number flows, which are less understood than these high Reynolds number flows. The reason for this is that new small, lightweight aircraft have smaller length scales than larger planes, and that new high-flying aircraft cruise in regions of the atmosphere that are less dense. It is important to remember at this point that the Reynolds number is defined by

Re = ρ u L / µ,

where ρ is the fluid density, u is the speed of the aircraft, L is a characteristic length scale, and µ is the dynamic viscosity of the fluid. The Reynolds numbers for these new aircraft are on the order of 10^{4} to 10^{5}, one or two orders of magnitude below Re for faster, conventional aircraft.

The issue is that these moderate Reynolds number flows (10^{4} ≤ Re ≤ 10^{5}) tend to start out laminar, but then separate from the upper, *suction* surface of wings. In general, flow separation is undesirable because it leads to significant increases in drag and a loss of lift. However, if an airfoil is in such a separated flow, it can increase its angle of attack to some critical angle, α_{crit}, where the flow over its upper surface transitions to turbulence and reattaches before the trailing edge. We can note at this point that α_{crit} is typically around 8º and is a function of Re and wing aspect ratio, AR [1]. When the flow reattaches, there is a sharp increase in the lift coefficient, C_{L}, drag coefficient, C_{D}, and lift-to-drag ratio, L/D. Visually-speaking, the flow starts out attached at the leading edge, then separates away from the wing for a while, and then becomes turbulent and reattaches. This “bubble” of air beneath th + e separated flow is called the “laminar separation bubble,” and is an important topic of study in many fluids communities.

As an aerospace engineer designing aircraft to operate at lower (or technically speaking, *moderate*) Reynolds numbers, it would be advantageous to know how to encourage laminar separated flow to reattach on a wing. Accomplishing this is done by *flow control*, the introduction of external forcing to make the flow move in a more favorable way. One way to control the flow over a wing is to apply a very thin band of metal at the leading edge that will heat up in patterns that are periodic in time and in space. The temporal frequency of this forcing can be denoted by the variable ω, and the spatial frequency of this forcing (in the spanwise direction) can be denoted by the variable k_{z}, which is also known as the *wavenumber*. The goal of resolvent analysis is to help us determine what forcing configuration (ie. what ω-k_{z} combination) will allow us to most effectively make the flow reattach before it passes the trailing edge.

The following description of resolvent analysis is taken from the wonderful 2019 paper by Professor Chi-An Yeh of UCLA [2]. To begin our discussion, we can start with the Navier-Stokes equations, which describe the movement of all fluid flows:

∂ρ/∂t + ∇ • (ρ **u**) = 0 (Eq. 1)

ρ ∂**u**/∂t + ρ (**u** • ∇) **u** = – ∇p + µ ∆**u** + ρ **g **+** f**** ^{ +}**. (Eq. 2)

The first equation, “∂ρ/∂t + ∇ • (ρ **u**) = 0,” tells us that the divergence of the velocity vector field is always equal to zero. This is the same as saying that mass is conserved, since zero divergence means that flow cannot appear out of nowhere or disappear into nothingness. The second, more fearsome equation, is really Newton’s Second Law in disguise. The terms on the right-hand-side are forces due to pressure, viscosity, gravity, and external flow control mechanisms, respectively. The terms on the left-hand-side are analogous to the product of mass (ρ is mass per unit volume), and acceleration. Indeed, the total time derivative of the velocity vector can be written as

D**u**/Dt = ∂**u**/∂t + (**u** • ∇) **u**

through a direct application of the multivariable chain rule. The Navier-Stokes equations are a system of nonlinear partial differential equations. Specifically, the nonlinearity appears in the term “(**u** • ∇) **u.**” This nonlinearities makes understanding and solving the Navier-Stokes equations extremely difficult, if not impossible. However, the ingenious strategy of resolvent analysis is to focus on the linear dynamics from the Navier-Stokes equations and model the nonlinear terms as internal forcing.

Let us define a state variable, **q**, which will represent everything we wish to know about our flow at a given position in space and time: the density, velocity in the x, y, and z-directions, and temperature.

**q** = [ρ, u, v, w, T]^{T}.

We can then compactly-rewrite Eq. 2 to read as

∂**q**/∂t = N(**q**) + **f **** ^{+}**, (Eq. 3)

where N is the (nonlinear) Navier-Stokes operator that acts on the state variable **q**. The next step is to perform a Reynolds decomposition and separate the state variable **q **into a mean component, **q**_{avg }, and a fluctuating component, **q**_{fluc}:

**q** = **q**_{avg} + **q**_{fluc}.

While variable **q**_{avg} is associated with the fundamental, average, *base* flow, the variable **q**_{fluc} is sensitive, and wiggles wildly throughout time. We fundamentally want our external forcing to influence the mean flow, and we can analyze how this can happen by focusing on the fluctuating component.

Let’s substitute our Reynolds decomposition of the state variable back into the Eq. 3 to obtain

∂**q**_{avg}/∂t + ∂**q**_{fluc}/∂t = L_{ }_{q}_{,avg} ( **q**_{fluc}) + N(**q**_{avg}) + h(**q**^{ n}_{fluc}) + **f**** ^{ +}**,

where L_{ }_{q,avg} is an operator that acts on **q**_{fluc} and captures all of the terms that are linear with respect to **q**_{fluc}, and “h(**q**^{ n}_{fluc})” represents the nonlinear higher-order terms in **q**_{fluc}. Since **q**_{avg} is a time-averaged value of the state variable, we set set ∂**q**_{avg}/∂t = 0. For the sake of compactness, we can also let

*U*** **=** **N(**q**_{avg}) + h(**q**^{ n}_{fluc}) + **f**** ^{ +}**.

Then, we can arrive at the following homogeneous linear differential equation for **q**_{fluc}:

∂**q**_{fluc}/∂t = L_{ }_{q}_{,avg} ( **q**_{fluc}) + *U*** **. (Eq. 4)

It is important to emphasize at this point that Eq. 4 represents the linearized Navier-Stokes equation for the fluctuating flow about the base flow defined by **q**_{avg}. As discussed in a previous blog post, Fourier transforms are a natural way to pass from time domains to frequency domains. Since we are interested in optimizing the temporal frequency, ω, and the spanwise spatial frequency, k_{z}, of our active flow control forcing, it makes sense for us to Fourier transform Eq. 4 with respect to time, t, and the spanwise spatial coordinate, z. Let **U**(k_{z}, ω) be the Fourier-transformed version of ** U**, and let q(k

_{z}, ω) be the Fourier-transformed version of

**q**

_{fluc}. Our equation of interest in Fourier space is then

-i ω q = L_{ }_{q}_{,avg} (q) + **U**. (Eq. 5)

The general solution for this equation in the variable q will be the sum of a homogeneous solution and a particular solution. Finding the homogenous solution amounts to solving the following equation, which is precisely an eigenvalue problem:

-i ω q = L_{ }_{q}_{,avg} (q). (Eq. 6)

Therefore, q is an eigenvector (otherwise known as an *eigenmode*) of L_{ }_{q}_{,avg}, and -i ω is the eigenvalue associated to q. The set of eigenvalues of L_{ }_{q}_{,avg} is called the *spectrum *of L_{ }_{q}_{,avg}, and this information is useful to know. However, because our flow is shear-dominated, L_{ }_{q}_{,avg} must be a non-normal operator, meaning that we need to know the pseudospectrum, in addition to the spectrum of L_{ }_{q}_{,avg} in order to have a more complete understanding of the dynamics at hand. (As a reminder, the pseudospectrum of an operator is the set of all eigenvalues of that operator, together with the set of all complex numbers that are “close to being eigenvalues”, in a sense that we will have to explore later on). To analyze the pseudospectrum of L_{ }_{q}_{,avg}, we can examine the particular solution to Eq. 5, which can be written as

q = [-i ω – L_{ }_{q}_{,avg}] ^{-1} **U**.** **(Eq. 7)

To ease notation, let’s let

**H** = [-i ω – L_{ }_{q}_{,avg}] ^{-1}, (Eq. 8)

and call **H **the *resolvent operator. *It turns out that analyzing **H** can give us information about the pseudospectrum of L_{ }_{q}_{,avg}.

Now suppose that we have discretized our flow using a Large Eddy Simulation (LES), and have created a matrix representation of L_{ }_{q}_{,avg}, which we will denote by **L**_{q}_{, avg}, and a matrix representation of **H**, which we will denote by **H. **Let’s also define the energy norm between two state vectors via the inner product

< **q**** _{1}**,

**q**

_{2}> =

**q**

_{1}

^{*}**W q**

**,**

_{2}where **q**_{1}^{*}** **denotes the conjugate transpose of **q**** _{1}** and

**W**is a positive-definite matrix that captures information about the energy of the flow and a numerical spatial integration around the wing. Since

**W**is positive-definite, there is a unique matrix

**W**

^{1/2}such that

**W**^{1/2 }**W**^{1/2} = **W**.

Now define

**H**** _{W}** =

**W**

^{1/2}

**H**

**W**

^{ –}^{1/2}.

Then, the pseudospectrum of L_{ }_{q}_{,avg} can be obtained via a singular value decomposition of **H**** _{W}** [3]. (Recall that a singular value decomposition essentially separates a matrix into its different rotating and stretching components).

**H**_{W}** **= **Q**_{W}** Σ U**^{*}** _{W}**. (Eq. 9)

We can interpret the right singular vectors in **U**** _{W}** as forcing modes and the left singular vectors in

**Q**

_{W}**as response modes. Furthermore, the singular vectors in the diagonal matrix**

**Σ**can be interpreted as amplification factors that scale the forcing modes. It turns out that usually the leading singular value, σ

_{1}is an order of magnitude larger than all of the other singular vectors. Hence, σ

_{1}can effectively quantify the gain of our system for a particular ω-k

_{z}pair. To determine the forcing configuration that best affects the flow, we can graph σ

_{1}for different values of ω and k

_{z}and see which frequency-wavenumber combination gives rise to the largest gain.

While there is a lot more to discuss, we will end our discussion at this (hopefully) climactic stopping point. In a future blog post, I look forward to talking more about fluid dynamics and delving deeper into all of the math that is at play here. Until then, take care.

**References:**

- https://www.researchgate.net/publication/354621720_Flow_transitions_on_a_cambered_airfoil_at_moderate_Reynolds_number
- https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/resolventanalysisbased-design-of-airfoil-separation-control/46E822C33CC66853ED4A8A0CFD51CCBB
- https://www.scirp.org/(S(i43dyn45teexjx455qlt3d2q))/reference/ReferencesPapers.aspx?ReferenceID=1161811

## Leave a Reply