A tale of two airfoils: resolvent-based modelling of an oscillator vs. an amplifier from an experimental mean

TL;DR

This study uses resolvent analysis combined with experimental data to model flow dynamics around a NACA 0018 airfoil at different angles of attack, revealing distinct oscillator and amplifier behaviors.

Contribution

It introduces a resolvent-based modeling approach that incorporates experimental mean flows and nonlinear forcing to distinguish flow mechanisms at different angles of attack.

Findings

A0 flow behaves as an oscillator with intrinsic instability.

A10 flow acts as an amplifier with unstructured forcing.

Resolvent modes align with SPOD modes when low-rank, and parasitic modes improve higher harmonic modeling.

Abstract

The flows around a NACA 0018 airfoil at a Reynolds number of 10250 and angles of attack of alpha = 0 (A0) and alpha = 10 (A10) are modelled using resolvent analysis and limited experimental measurements obtained from particle image velocimetry. The experimental mean velocity profiles are data-assimilated so that they are solutions of the incompressible Reynolds-averaged Navier-Stokes equations forced by Reynolds stress terms which are derived from experimental data. Spectral proper orthogonal decompositions (SPOD) of the velocity fluctuations and nonlinear forcing find low-rank behaviour at the shedding frequency and its higher harmonics for the A0 case. In the A10 case, low-rank behaviour is observed for the velocity fluctuations in two bands of frequencies. Resolvent analysis of the data-assimilated means identifies low-rank behaviour only in the vicinity of the shedding frequency for A0 and none of its harmonics. The resolvent operator for the A10 case, on the other hand, identifies two linear mechanisms whose frequencies are a close match with those identified by SPOD. It is also shown that the second linear mechanism, corresponding to the Kelvin-Helmholtz instability in the shear layer, cannot be identified just by considering the time-averaged experimental measurements as a mean flow due to the fact that experimental data are missing near the leading edge. The A0 case is classified as an oscillator where the flow is organised around an intrinsic instability while the A10 case behaves like an amplifier whose forcing is unstructured. For both cases, resolvent modes resemble those from SPOD when the operator is low-rank. To model the higher harmonics where this is not the case, we add parasitic resolvent modes, as opposed to classical resolvent modes which are the most amplified, by approximating the nonlinear forcing from limited triadic interactions of known resolvent modes.

Tables

Table 1. Table 1: Properties of the cross-spectral density tensors and resolvent operators for the A0 and A10 cases.

$\omega$-range	${𝑺}_{\widehat{𝒖}\widehat{𝒖}}\left(\omega \right)$ (SPOD)	$\mathscr{H}\left(\omega \right)$ (Resolvent)	${𝑺}_{\widehat{𝒇}\widehat{𝒇}}\left(\omega \right)$	modelling type
A0 case (oscillator)
1. near ${\omega }_s$	low-rank	low-rank	low-rank	first resolvent mode
2. near harmonics of ${\omega }_s$	low-rank	high-rank	low-rank	parasitic mode
3. all other $\omega$	high-rank	high-rank	high-rank	not modelled
A10 case (amplifier)
4. wake frequencies	low-rank	low-rank	high-rank	first resolvent mode
5. shear layer frequencies	low-rank	low-rank	high-rank	first resolvent mode
6. all other $\omega$	high-rank	high-rank	high-rank	not modelled

Table 2. Table 2: Error for the A0 case as a function of point selection.

1 Point	Error	2 Points	Error
1	344%	1 + 2	165%
2	988%	1 + 3	36%
3	152%	2 + 3	32%

Table 3. Table 3: Error for the A10 case as a function of point selection.

Point 4	Point 5	Points 4 $+$ 5
$\gg$ 100%	956%	47%

Equations

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+R{e}^{-1}{\nabla }^2\mathbf{u},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\overline{\mathbf{u}}\cdot \nabla \overline{\mathbf{u}}+\nabla \overline{p}-R{e}^{-1}{\nabla }^2\overline{\mathbf{u}}=\mathbf{f}=-\nabla \cdot \mathbf{R},$$

$$\nabla \cdot \overline{\mathbf{u}}=0.$$

$$\mathcal{L}({\overline{\mathbf{u}}}^{†},{\overline{p}}^{†},\overline{\mathbf{u}},\overline{p},\mathbf{f})=\mathcal{E}(\overline{\mathbf{u}})-⟨{\overline{\mathbf{u}}}^{†},\overline{\mathbf{u}}\cdot \nabla \overline{\mathbf{u}}+\nabla \overline{p}-R{e}^{-1}{\nabla }^2\overline{\mathbf{u}}-\mathbf{f}⟩-⟨{\overline{p}}^{†},\nabla \cdot \overline{\mathbf{u}}⟩,$$

$$\mathcal{E}(\overline{\mathbf{u}})=\frac{1}{2}\Vert \overline{\mathbf{m}}-\mathcal{M}(\overline{\mathbf{u}}){\Vert }_M^2.$$

$$\boldsymbol{Q}=\left[\begin{array}[]{cccc}\boldsymbol{q}_{1}&\boldsymbol{q}_{2}&\cdots&\boldsymbol{q}_{n}\end{array}\right]\in\mathbb{R}^{m\times n},$$

$$\boldsymbol{Q}^{(p)}=\left[\begin{array}[]{cccc}\boldsymbol{q}_{1}^{(p)}&\boldsymbol{q}_{2}^{(p)}&\cdots&\boldsymbol{q}^{(p)}_{n_{b}}\end{array}\right]\in\mathbb{R}^{m\times n_{b}}.$$

$$\hat{\boldsymbol{Q}}^{(p)}=\left[\begin{array}[]{cccc}\hat{\boldsymbol{q}}_{\omega_{1}}^{(p)}&\hat{\boldsymbol{q}}_{\omega_{2}}^{(p)}&\cdots&\hat{\boldsymbol{q}}_{\omega_{n_{b}}}^{(p)}\end{array}\right]\in\mathbb{C}^{m\times n_{b}}.$$

$$\hat{\boldsymbol{Q}}(\omega)=\left[\begin{array}[]{cccc}\hat{\boldsymbol{q}}_{\omega}^{(1)}&\hat{\boldsymbol{q}}_{\omega}^{(2)}&\cdots&\hat{\boldsymbol{q}}_{\omega}^{(b)}\end{array}\right]\in\mathbb{C}^{m\times b},$$

$$\hat{\mathbf{S}}(\omega )=\hat{\mathbf{Q}}(\omega ){\hat{\mathbf{Q}}}^{\ast }(\omega )\in {\mathbb{R}}^{m\times m}.$$

$$\hat{\mathbf{S}}(\omega )\hat{\mathbf{V}}(\omega )=\hat{\mathbf{V}}(\omega )\mathbf{\Lambda }(\omega ).$$

$${\hat{\mathbf{Q}}}^{\ast }(\omega )\hat{\mathbf{Q}}(\omega )\hat{\mathbf{\Theta }}(\omega )=\hat{\mathbf{\Theta }}(\omega )\mathbf{\Lambda }(\omega ),$$

$$\hat{\mathbf{V}}(\omega )=\hat{\mathbf{Q}}(\omega )\hat{\mathbf{\Theta }}(\omega )\mathbf{\Lambda }(\omega {)}^{-1/2}.$$

$${\partial }_t{\mathbf{u}}^{\prime }+\overline{\mathbf{u}}\cdot \nabla {\mathbf{u}}^{\prime }+{\mathbf{u}}^{\prime }\cdot \nabla \overline{\mathbf{u}}+\nabla p^{\prime }-R{e}^{-1}{\nabla }^2{\mathbf{u}}^{\prime }=-{\mathbf{u}}^{\prime }\cdot \nabla {\mathbf{u}}^{\prime }+\overline{{\mathbf{u}}^{\prime }\cdot \nabla {\mathbf{u}}^{\prime }}={\mathbf{f}}^{\prime },$$

$$\nabla \cdot {\mathbf{u}}^{\prime }=0.$$

$$\hat{\mathbf{u}}={\mathbf{C}}^T(i\omega \mathbf{B}-\mathbf{L}{)}^{-1}\mathbf{C}\hat{\mathbf{f}}=\mathcal{H}(\omega )\hat{\mathbf{f}}.$$

$$\boldsymbol{C}=\left(\begin{array}[]{c}1\\ 0\end{array}\right),~{}~{}~{}\boldsymbol{B}=\left(\begin{array}[]{cc}1&0\\ 0&0\end{array}\right),~{}~{}~{}\boldsymbol{L}=\left(\begin{array}[]{cc}-\overline{\boldsymbol{u}}\cdot\nabla()-()\cdot\nabla\overline{\boldsymbol{u}}+Re^{-1}\nabla^{2}()&-\nabla()\\ \nabla\cdot()&0\end{array}\right).$$

$$\mathcal{H}(\omega )=\sum _{j=1}^{\infty }\hat{\mathbf{\Psi }}(\omega )\mathbf{\Sigma }(\omega ){\hat{\mathbf{\Phi }}}^{\ast }(\omega ).$$

$${\mathcal{H}}^{\ast }(\omega )\mathcal{H}(\omega ){\hat{\mathbf{\varphi }}}_i(\omega )={\sigma }_i^2(\omega ){\hat{\mathbf{\varphi }}}_i(\omega ).$$

$$\mathcal{H}(\omega )\approx {\hat{\mathbf{\psi }}}_1(\omega ){\sigma }_1(\omega ){\hat{\mathbf{\varphi }}}_1(\omega ).$$

$$\hat{\mathbf{u}}(\omega )\approx {\hat{\mathbf{\psi }}}_1(\omega ){\sigma }_1(\omega ){\hat{\mathbf{\varphi }}}_1^{\ast }(\omega )\hat{\mathbf{f}}(\omega )={\hat{\mathbf{\psi }}}_1(\omega ){\chi }_1(\omega ),$$

$$\hat{\mathbf{u}}(\omega )\approx \mathcal{H}(\omega ){\hat{\mathbf{f}}}_a(\omega ),$$

$$\hat{\mathbf{f}}({\omega }_c)=\sum _{\begin{array}{c}{\omega }_a+{\omega }_b={\omega }_c\\ {\omega }_a,{\omega }_b\ne 0\end{array}}-\nabla \cdot (\hat{\mathbf{u}}({\omega }_a){\hat{\mathbf{u}}}^T({\omega }_b)).$$

$${\chi }_1(\omega )={\hat{\mathbf{\psi }}}_1(\omega ,{\mathbf{x}}_0)/\check{u}(\omega ,{\mathbf{x}}_0).$$

$${\mathbf{u}}^{\prime }(t)=\sum _{j=1}^{N_{\omega }}{\hat{\mathbf{\psi }}}_1({\omega }_j){\chi }_1({\omega }_j)e^{i({\omega }_{j}t)}.$$

$${\hat{\mathbf{S}}}_{\hat{\mathbf{u}}\hat{\mathbf{u}}}(\omega )=\mathbb{E}\{\hat{\mathbf{u}}(\mathbf{x},\omega ){\hat{\mathbf{u}}}^{\ast }(\mathbf{x},\omega )\}.$$

$${\hat{\mathbf{S}}}_{\hat{\mathbf{u}}\hat{\mathbf{u}}}(\omega )=\mathbb{E}\{\mathcal{H}(\omega )\hat{\mathbf{f}}(\mathbf{x},\omega ){\hat{\mathbf{f}}}^{\ast }(\mathbf{x},\omega ){\mathcal{H}}^{\ast }(\omega )\}.$$

$${\hat{\mathbf{S}}}_{\hat{\mathbf{u}}\hat{\mathbf{u}}}(\omega )=\mathcal{H}(\omega ){\hat{\mathbf{S}}}_{\hat{\mathbf{f}}\hat{\mathbf{f}}}(\omega ){\mathcal{H}}^{\ast }(\omega ),$$

$$\text{rank}({\hat{\mathbf{S}}}_{\hat{\mathbf{u}}\hat{\mathbf{u}}}(\omega ))\le \text{min}\left[\text{rank}(\mathcal{H}(\omega )),\text{rank}({\hat{\mathbf{S}}}_{\hat{f}\hat{f}}(\omega ))\right].$$

$${\hat{\mathbf{f}}}_a(2{\omega }_s)\sim {\hat{\mathbf{\psi }}}_1({\omega }_s)\cdot \nabla {\hat{\mathbf{\psi }}}_1({\omega }_s).$$

$$\hat{\mathbf{f}}(\omega =3.70)=-\sum _{\begin{array}{c}{\omega }_a+{\omega }_b=3.70\\ {\omega }_a,{\omega }_b\ne 0\end{array}}{\chi }_1({\omega }_a){\chi }_1({\omega }_b)\nabla \cdot ({\hat{\mathbf{\psi }}}_1({\omega }_a){\hat{\mathbf{\psi }}}_1^T({\omega }_b)).$$

Full text

1 Introduction

Recent studies, e.g. Gómez et al. (2016a); Beneddine et al. (2017); Symon (2018); Symon & McKeon (2018); He et al. (2019), have shown that it is possible to construct a reduced-order model of a flow using the time-averaged (mean) flow field and a small set of time-resolved velocity data at isolated measurement points. For low Reynolds number bluff body wakes, this could be as few as a single point (see Gómez et al., 2016b; Symon et al., 2019). In a similar vein, Beneddine et al. (2016) and Thomareis & Papadakis (2018) have been able to estimate the spectra at various points of the flow using only the aforementioned measurements. The success of the method relies on the low-rank nature of the (linear) resolvent operator, which one obtains after linearising the Navier-Stokes equations around the (turbulent) mean flow, and treating the nonlinear terms as a source of intrinsic forcing (McKeon & Sharma, 2010). The exact form of the forcing can often be neglected due to a dominant amplification mechanism at particular spatial wavenumbers and frequencies. Linear analysis alone does not reveal the amplitude and phase of the most amplified disturbances yet they can be calibrated using limited unsteady data (Gómez et al., 2016a; Beneddine et al., 2016).

One of the major advantages of this approach is that time-resolved data at all locations in the flow field is not required. From the governing equations, it is possible to obtain a prediction of the temporal Fourier mode in the actual fluctuation field if the dominant singular value of the resolvent operator is much larger than the others. The nonlinear forcing field, furthermore, must not be preferentially biased towards the forcing mode corresponding to lower singular values and this was demonstrated by Beneddine et al. (2016) for turbulent flow over a backward-facing step. The goal of Dynamic Mode Decomposition (DMD) (Schmid, 2010; Rowley et al., 2009) is also to compute modes which oscillate at a distinct frequency, suggesting that resolvent analysis is attempting to approximate DMD modes from the equations instead of time-resolved data. A striking similarity between the most amplified resolvent mode and the DMD mode was noted by Gómez et al. (2014) in turbulent pipe flow and the relationship between resolvent analysis and DMD was formalised in Sharma et al. (2016).

Connections between resolvent analysis and other data-driven approaches have been discussed by Towne et al. (2018). In particular, they examined modes extracted by Spectral Proper Orthogonal Decomposition (SPOD), introduced by Lumley (1970), which produces a basis for the fluctuation field at each temporal frequency and optimally represents the space-time flow statistics. Towne et al. (2018) determined that SPOD modes and resolvent modes are equivalent if the nonlinear forcing field is white in space and time. In practice, this might not be the case since the resolvent framework is a closed loop (McKeon et al., 2013; McKeon, 2017) where the nonlinear interactions between the velocity fluctuations constitute the forcing to the resolvent operator. Despite this fact, resolvent analysis is quite often able to predict the structure of energetic velocity fluctuations. Moreover, if the singular values of the resolvent operator are well separated, it often suggests a similar separation in the SPOD eigenspectra as demonstrated by Schmidt et al. (2018) for a turbulent jet.

Resolvent analysis may be used for predicting structures even if the singular values are not well separated. Recent work by Rosenberg et al. (2019) proved that considering a small subset of triadic interactions of ‘correct’ resolvent modes (i.e. modes computed at frequencies where the singular values are well separated) can approximate the nonlinear forcing, a process reminiscent of weakly nonlinear analysis (Sipp & Lebedev, 2007) without being limited to the vicinity of a critical Reynolds number. When these interactions are fed back through the resolvent operator, they produce the correct structure. These modes were labelled parasitic in the sense that they fed off nonlinear interactions of other modes rather than emerge naturally from a linear amplification mechanism embedded in the operator.

The majority of the aforementioned studies have utilised a mean obtained from simulations although a growing number have considered experimental means (Sasaki et al., 2017; Beneddine et al., 2017; Symon & McKeon, 2018; He et al., 2019). The objective in this article is to obtain a reduced-order representation using a partially known mean obtained from experiments. We use data-assimilation to recover the mean profile on a larger domain and enforce no-slip boundary conditions on the airfoil surface. Data-assimilation is a process whereby experimental measurements are merged with computational fluid dynamics (CFD) (see Hayase, 2015) and has been a growing field in fluid mechanics. A list of recent developments can be found in da Silva (2019). The particular framework used here was developed by Foures et al. (2014) and later applied to flow around an idealised airfoil by Symon et al. (2017).

In this article, we consider the flow around a NACA 0018 airfoil at a chord-based Reynolds number of $Re=10250$ at two angles of attack: $\alpha=0^{\circ}$ (A0 case) and $\alpha=10^{\circ}$ (A10 case). The flows are obtained experimentally from particle image velocimetry (PIV). As will be discussed later in the article, the A0 and A10 cases behave like an oscillator and amplifier, respectively. The former has intrinsic dynamics which are insensitive to background noise while the latter filters and amplifies upstream noise in the downstream direction (Huerre & Rossi, 1998; Sipp et al., 2010). Furthermore, the resolvent norm, or largest singular value of the resolvent operator, tends to have sharp peaks isolated at resonant frequencies in the cases of oscillators, e.g. Symon et al. (2018). Amplifier flows, on the other hand, may have separation between the first and second singular values over a broad ranges of frequencies, e.g. Dergham et al. (2013); Beneddine et al. (2016). Of particular interest to us is the nature of the nonlinear forcing for each type of flow. If the nonlinear forcing is structured, then it can be deduced from a limited number of nonlinear interactions between coherent structures. This was the case in the oscillator flows studied by Rosenberg et al. (2019), where the nonlinear forcing at energetic frequencies consisted of coherent structure at the oscillation frequency and its harmonics. In the case of a turbulent jet (Towne et al., 2015), which is more representative of an amplifier flow, the nonlinear forcing field resembles incoherent turbulent fluctuations.

The principal objectives of this article can be summarized as follows: (1) does data-assimilation of the mean velocity profile make a measurable difference on predicting the velocity fluctuations? (2) what does SPOD reveal when applied to the nonlinear forcing? (3) can additional structure be extracted from resolvent analysis if the singular values are not separated but the nonlinear forcing is approximated? (4) is there a fundamental difference between oscillator and amplifier type flows? (5) how far can we go in modelling with just a few mean flow measurements and some unsteady point measurements? These questions are addressed in the rest of the paper which is organised as follows. In §2, we describe the PIV experiments and in §3 we discuss the mathematical background of the analysis tools used to construct the resolvent-based model. In §4, the mean velocity profiles are data-assimilated onto a larger, more resolved mesh. In §5, we apply SPOD to the velocity fluctuations and nonlinear forcing of both airfoils to highlight the low-rank behaviour in each flow. In §6, the data-assimilated mean profiles are used as an input to resolvent analysis and the results are compared to those of SPOD. We also try to establish new links between resolvent analysis and SPOD which depend on whether the flow is an oscillator or an amplifier. In §7, we use resolvent modes and, where appropriate, parasitic modes to construct a model of the coherent fluctuations in the two airfoil flows before concluding in §8.

2 Experimental setup

The flow around a NACA 0018 airfoil at a Reynolds number of $Re=10250$ and two angles of attack, $\alpha=0^{\circ}$ (A0) and $\alpha=10^{\circ}$ (A10), is measured using PIV. The experiments are performed in a free-surface water facility with a test section measuring 1.6 m in length, 0.46 m in width and 0.5 m in height. The NACA 0018 airfoil is mounted vertically in the tunnel, so that its span is parallel to the test section height, and has a chord length of 10 cm. Its spanwise length is 48 cm, resulting in an aspect ratio of $AR=4.8$. The laser sheet is provided by a YLF dual cavity solid-state laser and is centred at a height of 220 mm. Dunne (2016) measured the out-of-plane velocity component to be less than 4% of the free-stream velocity, allowing the mean flow at this height to be considered two-dimensional.

The PIV setup consists of two Phantom Miro 320 cameras with 50 mm focal length Nikon lenses and 1:1.8 aperture. They have an overlap of 18% in the streamwise direction and sample the flow at a frequency of 125 Hz. The camera resolution is 1920 $\times$ 1200 pixels and they are calibrated at 8.2 px/mm. The seeding particles are hollow glass spheres (reference 110P8 with an average diameter of 11.7 $\mu$m and a specific gravity of 1.1) and the seeding density is approximately 0.1 particles per square pixel. To reduce the impact of surface reflections, the image intensity is calibrated using white-image subtraction and background-image subtraction, which are taken before each run and averaged over 100 snapshots.

LaVision’s software package DaVis computes the velocity vectors where a standard cross-correlation technique via Fast Fourier Transformation is applied to each sequential image. The window-size is reduced from 32 $\times$ 32 px2 to 16 $\times$ 16 px2 over three passes. Once the data are post-processed, missing vectors are interpolated using an average of all non-zero neighbourhood vectors and a median filter (Westerweel & Scarano, 2005) is used for outlier detection. A sample of the post-processed data for both cases is displayed in figure 1. The mean profile, obtained after time-averaging two runs of 3,499 snapshots, is subtracted so that the coherent structures can be more readily identified. As seen in figure 1, velocity vectors below the airfoil cannot be obtained since the laser sheet is obstructed in this region. For the A0 case, the mean velocity profile is symmetric so the mean profile can be reflected over the centreline (the transverse velocity changes sign) to obtain the mean below the centreline and in the shadow region. This is not possible for the A10 case, so the mean profiles for $\alpha=\pm 10^{\circ}$, which are obtained separately during different runs, are stitched together to obtain data around the entire airfoil.

3 Analysis tools

3.1 Data-assimilation of the mean velocity profiles

The experimental profiles are data-assimilated using the framework devised by Foures et al. (2014) so that they can later be used as an input to resolvent analysis. Only the equations are reviewed here as its application to experimental data is discussed by Symon et al. (2017). The flow is governed by the incompressible Navier-Stokes equations, which are non-dimensionalised by the free-stream velocity $U_{\infty}$ and the airfoil chord $c$:

[TABLE]

where $\boldsymbol{u}$ and $p$ are the velocity and pressure, respectively. Substituting $\boldsymbol{u}=\overline{\boldsymbol{u}}+\boldsymbol{u}^{\prime}$, where an overbar denotes a mean and a prime denotes a fluctuation, into (1) and averaging in time yields the incompressible Reynolds-averaged Navier-Stokes (RANS) equations:

[TABLE]

The right-hand side of (2a), or $\boldsymbol{f}$, represents the divergence of the Reynolds stress tensor $\boldsymbol{R}$.

The data-assimilation algorithm treats $\boldsymbol{f}$ as an unknown and leverages partial knowledge of $\overline{\boldsymbol{u}}$ to reconstruct not only $\boldsymbol{f}$ but also $\overline{\boldsymbol{u}}$ on a larger domain which obeys prescribed boundary conditions. A constrained optimisation problem is formulated as follows

[TABLE]

to minimise the difference between a measurement of $\overline{\boldsymbol{u}}$, which satisfies (2), and the experimental mean profile $\overline{\boldsymbol{u}}_{exp}$ measurement. In (3), $(\cdot)^{\dagger}$ signifies an adjoint quantity and $\left<\boldsymbol{a}\cdot\boldsymbol{b}\right>=\int_{\Omega}\boldsymbol{a}\cdot\boldsymbol{b}d\Omega$. The scalar $\mathcal{E}$ quantifies the discrepancy and is given by

[TABLE]

The operator $\mathcal{M}$ projects the numerical data to the subspace spanned by the known PIV measurements $\overline{\boldsymbol{m}}$.

The variations of $\mathcal{L}$ with respect to each dependent variable are set to zero to obtain the equations which are solved iteratively until $\mathcal{E}$ is minimised. Since $\mathcal{E}\neq\mathcal{E}(\overline{p})$, only the solenoidal component of $\boldsymbol{f}$ is recovered (see Symon et al., 2019, for details). To remove the irrotational component and facilitate comparisons to the experimental Reynolds stresses, only $\nabla\times\boldsymbol{f}$ is reported in §4.

The direct and adjoint equations are solved in FreeFem++ (Hecht, 2012) on the computational domain $\Omega\in-5\leq x\leq 15\cup-5\leq y\leq 5$ with the leading edge of the airfoil being centred at the origin. The equations are spatially discretised using quadratic basis functions for the velocity and linear basis functions for the pressure resulting in approximately 360,000 and 1,000,000 degrees of freedom for the A0 and A10 cases, respectively. The reader is referred to Foures et al. (2014); Symon et al. (2017) for more details about the algorithm and its implementation.

3.2 Spectral proper orthogonal decomposition

The velocity fluctuations obtained from the processed PIV datasets are analysed using two modal decomposition techniques. The first, introduced by Lumley (1967, 1970), is spectral proper orthogonal decomposition (SPOD). From the raw, experimental data, SPOD calculates energy-ranked modes which are orthogonal to one another and oscillate at a single frequency. The SPOD modes are computed using the procedure described in Towne et al. (2018) and the code developed by Schmidt & Towne (2019). A brief summary of the steps is presented below.

The PIV data are rearranged into the matrix

[TABLE]

where $m$ represents the number of states, i.e. the streamwise and transverse velocity of all PIV vectors, and $n$ is the total number of snapshots. Thus, the overall size of $\boldsymbol{Q}$ is 49,104 total states by 6,800 snapshots. Following Welch’s method (Welch, 1967), $\boldsymbol{Q}$ is split into $p$ smaller, overlapping segments, or blocks, containing $n_{b}$ snapshots where $n_{b}<n$

[TABLE]

Each row of $\boldsymbol{Q}^{(p)}$ is then Fourier-transformed in time using the discrete Fourier transform (DFT), yielding an ensemble of Fourier realizations

[TABLE]

The columns of $\hat{\boldsymbol{Q}}^{(p)}$ represent a Fourier mode instead of a snapshot in time as they did in $\boldsymbol{Q}^{(p)}$. The Fourier modes for a specific frequency $\omega$ can be arranged into a new data matrix

[TABLE]

where $b$ is the number of blocks.

The cross-spectral density matrix for a specific temporal frequency $\hat{\boldsymbol{S}}(\omega)$ is

[TABLE]

The SPOD eigenvectors $\hat{\boldsymbol{V}}(\omega)$ and eigenvalues $\boldsymbol{\Lambda}(\omega)$ can be obtained via an eigenvalue decomposition of the cross-spectral density matrix

[TABLE]

Since the number of blocks is much smaller than the number of states in the problem, i.e. $b\ll m$, it is computationally less demanding to consider the following eigenvalue problem instead

[TABLE]

as $\hat{\boldsymbol{Q}}^{*}(\omega)\hat{\boldsymbol{Q}}(\omega)\in\mathbb{R}^{b\times b}$, and compute the eigenvector of the full system using

[TABLE]

It should be noted that the choice of $\omega$ is limited by $n$ and the sampling frequency which, for both experiments, is 125Hz. As will be seen in the next section, this is a major disadvantage of SPOD, which requires a lot of data to obtain converged modes and high frequency resolution, in comparison to resolvent analysis.

The PIV snapshots are divided into blocks containing $n=1700$ snapshots with 50% overlap, resulting in $b=7$ realisations of the flow. The resulting frequency resolution is $\Delta\omega=0.46$. SPOD is performed using a rectangular windowing function (nearly identical results are obtained using a Hamming window). The choice of $n$ and $b$ is a compromise between suitable resolution of the lower frequencies, where the energy is most concentrated, and obtaining a sufficient number of realisations to identify low-rank behaviour. We find that the trends are unaffected as long as we consider at least 3 realisations of the flow.

3.3 Resolvent analysis

Once the mean profile has been data-assimilated, it can be used as an input to resolvent analysis, which has been shown to provide an efficient basis for velocity fluctuations by McKeon & Sharma (2010) and many others. Subtracting (2) from (1) results in the equations governing the fluctuations

[TABLE]

The left-hand side of (13a) contains the linear terms and the right-hand side is lumped into a single forcing term $\boldsymbol{f}^{\prime}$ to account for the nonlinear terms. Substituting $\boldsymbol{f}^{\prime}=\hat{\boldsymbol{f}}e^{i\omega t}$, where $\omega$ is a real-valued frequency, implies $\boldsymbol{u}^{\prime}=\hat{\boldsymbol{u}}e^{i\omega t}$ and leads to the following input-output form of the linearised Navier-Stokes equations

[TABLE]

The matrix $\mathcal{H}(\omega)$ is the resolvent operator and the linear operators $\boldsymbol{C}$, $\boldsymbol{B}$, and $\boldsymbol{L}$ are defined as

[TABLE]

$\boldsymbol{L}$ is the linear Navier-Stokes (LNS) operator, $\boldsymbol{B}$ restricts the nonlinear forcing to the velocity subspace, and $\boldsymbol{C}$ can be adjusted to restrict the spatial domain of the input and output vectors.

The singular value decomposition (SVD) of the resolvent operator yields basis functions for the nonlinear term and velocity fluctuations at each temporal frequency $\omega$

[TABLE]

The matrices $\hat{\boldsymbol{\Phi}}(\omega)$ and $\hat{\boldsymbol{\Psi}}(\omega)$ contain the optimal forcing and response modes, respectively, which are ranked by their kinetic energy gains contained in the diagonal matrix $\boldsymbol{\Sigma}(\omega)$. Following Sipp & Marquet (2013), the singular values $\sigma_{i}(\omega)$ are computed by generating the relevant operators in FreeFem++ and reformulating (16) into the following eigenvalue problem:

[TABLE]

The largest eigenvalues of the Hermitian operator $\mathcal{H}^{*}(\omega)\mathcal{H}(\omega)$ are computed using the ARPACK library (Lehoucq & Sorensen, 1996) and the MUMPS parallel solver (Amestoy et al., 2001).

3.4 Reduced-order modelling of the fluctuation fields

In this section, we discuss how an efficient basis for the fluctuations can be educed from resolvent analysis by approximating the linear operator $\mathcal{H}(\omega)$, if it is low-rank, or the nonlinear forcing $\hat{\boldsymbol{f}}(\omega)$. We also explain how time-resolved probe points are used to calibrate the amplitude and phase of the resolvent modes.

3.4.1 Approximation of the (linear) resolvent operator

The resolvent operator’s singular values are said to be well separated if $\sigma_{1}(\omega)\gg\sigma_{2}(\omega)$, in which case the following approximation holds

[TABLE]

(18) can be interpreted as an approximation of the (linear) resolvent operator and signifies that the first response mode $\hat{\boldsymbol{\psi}}_{1}(\omega)$ is sufficient to represent the velocity fluctuations at a frequency $\omega$. In other words, it is possible to write

[TABLE]

where the procedure for computing the complex weight $\chi_{1}(\omega)$ is described in §3.4.3. It should be noted that, in general, the Fourier mode $\hat{\boldsymbol{u}}(\omega)$ is not well defined unless the flow is truly periodic like cylinder flow. Hence, we refer to $\hat{\boldsymbol{u}}(\omega)$ and $\hat{\boldsymbol{f}}(\omega)$ as an arbitrary sample related to a specific bin.

3.4.2 Approximation of the (nonlinear) forcing

When there is an energetic frequency in the flow where the singular values are not well separated, the resolvent modes are no longer an efficient basis to represent the velocity fluctuations. An alternative approach is to approximate the nonlinear forcing $\hat{\boldsymbol{f}}(\omega)$ instead of the linear operator itself (Rosenberg et al., 2019). (19) is modified into the following form

[TABLE]

where $\hat{\boldsymbol{f}}_{a}(\omega)$ is the approximated nonlinear forcing that yields the parasitic mode. As stated in McKeon et al. (2013), the full expression for $\hat{\boldsymbol{f}}(\omega)$ is a convolution over all triadically consistent Fourier modes

[TABLE]

In this article, we will explain the conditions under which the sum in (21) can be severely truncated to a few key triadic interactions and for which types of flows this procedure is applicable.

3.4.3 Weighting the modes

Once a resolvent or parasitic mode for each frequency of interest has been computed using one of the two methods described above, its complex weight $\chi_{1}(\omega)$ is calibrated using a time-resolved probe $\check{u}(t,\boldsymbol{x}_{0})$ in a region of the flow where $\hat{\boldsymbol{\psi}}_{1}(\omega)$ is energetic as demonstrated by Gómez et al. (2016a) and Beneddine et al. (2016). The signal is Fourier-transformed in time so that $\chi_{1}(\omega)$ can be computed in the frequency domain through the following expression

[TABLE]

It should be noted that while $\hat{\boldsymbol{\psi}}_{1}$ has been used in (22), the equation is still applicable to parasitic modes once their shapes have been computed using (20). After the complex amplitude of each mode has been determined, the linear superposition of all modes can be inverse Fourier-transformed to obtain a reduced-order model of the velocity fluctuations in the time domain, i.e.

[TABLE]

3.5 Summary

The resolvent-based modelling procedure is summarised in figure 2, which contains the relevant steps starting with the acquisition of the experimental measurements (top row). The mean profile is data-assimilated onto a FreeFem mesh with the correct boundary conditions and the time-resolved probe measurements are Fourier-transformed in time (row 2). The structure associated with frequency $\omega$ is obtained from resolvent analysis of the data-assimilated mean profile in one of two ways: (1) approximating the operator via the singular value decomposition if there is a large separation of singular values ($\sigma_{1}\gg\sigma_{2}$) or (2) approximating the nonlinear forcing if the frequency is energetic but there is no separation of singular values. These two methods are demarcated by solid and dotted lines, respectively, in row 3 of figure 2. The resolvent modes are weighted from the Fourier-transformed probe points and inverse Fourier-transformed to obtain a model of the velocity fluctuations in time (bottom row).

4 Assimilated flow fields

In this section, the mean profiles obtained from the time-averaged PIV are data-assimilated using the method introduced by Foures et al. (2014). Since data-assimilation attempts to reconstruct the optimal forcing $\boldsymbol{f}$ in the RANS equations, the quantity $\nabla\times\boldsymbol{f}$, in addition to the assimilated mean fields, are compared to their experimental counterparts.

4.1 A0 Case

As mentioned in §2, the data below the airfoil are obstructed by the airfoil’s shadow. To overcome this problem, the data above the centreline are reflected across, keeping the streamwise velocity constant and reversing the sign of the transverse velocity. Both the streamwise and transverse component of the data-assimilated fields are compared to their experimental counterparts in figure 3. The fields are nearly indistinguishable although the data-assimilated contours are slightly more smooth than those in the experiment.

The curl of the forcing $\nabla\times\boldsymbol{f}$, rather than $\boldsymbol{f}$ itself, is shown in figure 4 to eliminate the irrotational component which cannot be captured. The agreement between the experimental field, obtained after computing two gradients of the velocity fluctuations, and the data-assimilated field is good despite some minor discrepancies. One notable aspect about figure 4 is that $\nabla\times\boldsymbol{f}$ for this flow closely resembles that of cylinder flow (see figures 4 and 5.9 of Foures et al. (2014) and Symon et al. (2018), respectively).

4.2 A10 Case

Data-assimilation of the A10 case is slightly more difficult since the flow is not symmetric with respect to the centreline. Two separate experiments are conducted at both $\alpha=-10^{\circ}$ and $\alpha=10^{\circ}$. In the first experiment, the data below the airfoil (suction side), are obstructed by the airfoil shadow so only the pressure side can be measured. In the second experiment, the pressure side is obstructed by the shadow, so only the suction side is measured. The final mean velocity field can be obtained by stitching together the two mean profiles such that the data around the entire airfoil are known. The stitching process, if imperfect, may result in non-physical measurements although the data-assimilation algorithm is able to compensate for this.

Data-assimilation is performed using the stitched profile, i.e. from measurements at $\alpha=-10^{\circ}$ and $\alpha=10^{\circ}$, and the resulting fields are compared to their experimental counterparts in figure 5. In addition to comparing the two components of the mean velocity field, the mean spanwise vorticity $\overline{\omega}_{z}$ is plotted in the bottom row of figure 5. Streamlines of the mean velocity field indicate the size and location of the recirculation bubble. The streamwise and transverse velocities match very well, particularly with respect to the reverse flow region. The separated flow region is accurately reproduced in the data-assimilated fields.

There are two regions of the flow field which are visibly improved by data-assimilation. The first is near the leading edge of the airfoil, where the flow separates and vortices are formed due to the rollup of the shear layer. This improvement is particularly visible in the vorticity plots since the negative contours of the data-assimilated field are continuous unlike in the experiment where there is a gap at $x\approx 0.2$. As noted by Raffel et al. (2018), there are difficulties associated with PIV for laminar flow around airfoils. The effects of strong centrifugal forces around the airfoil leading edge and strong shear result in the outward movement of tracer particles in a direction perpendicular to the curved streamlines. This loss of seeding compromises measurements in the region close to the solid body.

The second area of improvement is at the trailing edge where there are imperfections associated with the stitching together of the mean profiles. Contours of negative mean streamwise velocity protrude below the trailing edge which corresponds to non-physical behaviour. Data-assimilation is able to correct this problem and eliminate inaccuracies associated with the stitching process. These changes are more apparent by examining the difference between the assimilated and experimental mean velocities, i.e. $\Delta\overline{\boldsymbol{u}}=\overline{\boldsymbol{u}}-\overline{\boldsymbol{u}}_{exp}$, which is illustrated in figure 6 for both velocity components. The two regions of disagreement, notably, coincide with the areas of improvement discussed earlier.

It should be emphasised that the benefit to data-assimilating fields is not only filling in missing data, but also in obtaining gradient quantities with less noise. This is particularly noticeable when comparing $\nabla\times\boldsymbol{f}$ in figure 7. Even though the experimental data do not have much noise to begin with, the data-assimilated field is more smooth overall.

Data-assimilation leads to improvement for both the A0 and A10 cases although the changes for A0 are fairly unremarkable. For the A10 case, on the other hand, there are features missing from the experimental profile that data-assimilation fills in, particularly at the shear layer originating near the leading edge. As will be seen in §6, this has major implications for the identification of linear amplification mechanisms in the flow.

5 SPOD results and the difference between an oscillator and an amplifier

In this section, SPOD is applied to the velocity fluctuations obtained from the raw PIV data. We also apply SPOD to the nonlinear term $\boldsymbol{u}^{\prime}\cdot\nabla\boldsymbol{u}^{\prime}$ to investigate the most energetic structures associated with the nonlinear forcing field. Towne et al. (2015) and Towne (2016) applied SPOD to the nonlinear forcing in a turbulent jet and found that the educed structures were primarily incoherent turbulent fluctuations and not nonlinear interactions between coherent structures. We will compare the results of SPOD applied to the velocity fluctuations and nonlinear forcing for both airfoils and highlight fundamental differences between oscillator and amplifier flows. These will have consequences on when we expect SPOD and resolvent analysis to coincide in §6 as well as our ability to model the flows in §7.

5.1 SPOD of the velocity fluctuations

The SPOD eigenspectra for A0 are plotted in figure 8 as a function of $\omega$. The largest peak in $\lambda_{1}$ occurs at $\omega=12.0$, which corresponds to the shedding frequency of the airfoil. Red patches highlight frequencies at which there is significant separation between the first and second eigenvalues, the former of which exhibits sharp peaks at the shedding frequency and its harmonics.

The $v$-component of the first two SPOD modes is shown in figure 9 for the first three harmonics. In all cases, the first eigenvalue is at least one order of magnitude larger than the second eigenvalue. The first SPOD mode is significantly more structured than the second mode for all three harmonics in figure 9, an observation which is consistent with the fact that nearly all the energy at these frequencies is concentrated in the first mode. For frequencies where $\lambda_{1}\approx\lambda_{2}$ (not shown in the interest of brevity), the mode shapes are less structured for all eigenvalues.

Let us now consider the A10 case. The PIV data are subdivided into the same number of blocks containing the same number of snapshots as in the A0 case to allow for direct comparisons between both angles of attack. The SPOD eigenspectra are plotted in figure 8. It is immediately apparent that the SPOD eigenvalues are far less peaked than they were for the A0 case. Instead of observing a fundamental frequency and its harmonics, there are bands of frequencies where low-rank behaviour can be observed, a feature reminiscent of an amplifier flow. The first region, denoted by the red patch in figure 8, is situated around $4<\omega<8$ while the second one resides roughly between $10.5<\omega<14.5$.

The modes for A10 are plotted in figure 10 for the highest eigenvalues in the low- and high-frequency bands. Similar to the A0 case, the leading SPOD mode is highly structured and contains about one order of magnitude more energy than the second mode for $\omega=6.01$ and $\omega=7.85$. Most of the energy in these modes, furthermore, resides in the wake region, or $x>1$, although there is some structure in the shear layer. The second SPOD mode, on the other hand, is more energetic in the shear layer although it is far less structured. The final frequency considered is $\omega=12.47$, which is roughly in the middle of the second band of frequencies identified earlier. The first mode is primarily in the shear layer while the second mode resembles a mix of both the shear layer and wake dynamics. Since SPOD takes into consideration the kinetic energy of the entire domain, the SPOD eigenvalues for frequencies corresponding to the shear layer assume lower values. They are dynamically important, nonetheless, since there is non-negligible separation between the first and second eigenvalue. They are also essential to model the fluctuations above the airfoil.

5.2 SPOD of the nonlinear forcing

We now utilise SPOD to characterise the structure of the nonlinear forcing, which is computed from the PIV snapshots since the data are sufficiently well-resolved in space and time. The SPOD eigenspectra for both cases (A0 & A10) are plotted together in figure 11 to underscore the contrast between the two cases. The eigenspectra for A0 are still peaked in the sense that the energy is high at the vortex shedding frequency and its harmonics. There is also a considerable gap between the first and second eigenvalues at these frequencies, similar to the trends observed for the velocity fluctuations. The spectra for the A10 case, on the other hand, are relatively flat and have no discernible peaks. The behaviour of the nonlinear forcing bears no resemblance to those of the velocity fluctuations as seen in figure 10.

The differences are even more evident when considering the mode shapes. For clarity, we will plot modes at the same frequencies as those in figures 9 and 10 since they are dynamically important. The modes for A0 are plotted in figure 12. It is clear that the nonlinear forcing for the first two harmonics is very structured, as the contours are almost perfectly symmetric (or anti-symmetric in the case of the second harmonic) with respect to the centreline. The structure corresponding to the third harmonic looks slightly more disordered, since it is not particularly energetic, but it can be remarked that the mode shape is approximately symmetric with respect to the centreline and has a well-defined wavelength which is approximately constant, i.e. not changing as a function of $x$.

The mode shapes for the A10 case are plotted in figure 12. The results are reminiscent of those from the jet flow study of Towne et al. (2015) in that there is a distinct lack of spatial coherence. While there may be some notion of a wavelength associated with these structures, it would be far more difficult to quantify in comparison to the results in figure 12. There is also no noticeable difference between the spatial support of the lower frequency modes in figure 13 and that of the higher frequency mode even though the choice of frequency had an impact on the spatial location and distribution of the velocity fluctuations.

5.3 Discussion and classification of the two cases

SPOD of the velocity fluctuations reveals low-rank behaviour for both the A0 and A10 cases although there are some major differences. The SPOD eigenspectra of A0 are highly peaked at the vortex shedding frequency and its harmonics. The A0 case, thus, can be considered sparse in that there are few energetic frequencies in the flow. The SPOD eigenspectra for A10 is more broadband, with separation between the first and second eigenvalues occurring at two distinct ranges of frequencies making the flow less sparse, i.e. more frequencies are required to model the fluctuating flow field. Within the lower range, the SPOD modes are concentrated in the wake whereas within the higher range, they are concentrated in the shear layer. Based on these trends alone, it seems that the A0 case resembles an oscillator-type flow while the A10 case behaves more like an amplifier since, in the case of an oscillator, the fluctuations are characterised by a well-defined frequency which can be readily identified in figure 9. There is no clear peak in figure 10.

Application of SPOD to the nonlinear forcing confirms this appraisal. In the case of A0, there is low-rank behaviour at all harmonics and the mode shapes are highly structured as they contain a well-defined wavelength and symmetry. This suggests they can be represented by a limited number of nonlinear interactions whereas in the A10 case, there is no apparent low-rank behaviour and the forcing field contains no coherent structure. These results are consistent with those of Rosenberg et al. (2019), who show that the forcing for the cylinder wake is highly structured at the shedding frequency and its harmonics, as well as Towne et al. (2015), who illustrate that the nonlinear forcing is comprised of incoherent turbulent motions for a turbulent jet. It is also in tune with the idea that the input to an amplifier is less important, as the linear dynamics have a preferential bias towards one particular structure.

In theory, this makes modelling of amplifier flows easier as the linear dynamics will do the work for us while nonlinearity needs to be taken into account for oscillator flows. Consequently, amplifier flows are more amenable to stochastic forcing of an approximated linear operator. Alternatively, the modelling of oscillator flows may benefit from approximating the nonlinear forcing due to the emergence of parasitic modes. These ideas are tested in §7.

6 Resolvent analysis of the data-assimilated profiles

For each angle of attack, two mean profiles are used as an input to resolvent analysis. The first is from the raw experimental mean, which is interpolated onto the same FreeFem mesh used in the data-assimilation algorithm. Due to the sparsity of data near the surface of the airfoil, the interpolated mean profile is not likely to satisfy the no-slip boundary condition, so this is manually enforced by setting the velocity at all points along the airfoil surface to zero. The points outside the PIV domain, where velocity vectors are not available, are set to uniform flow, i.e. $\overline{\boldsymbol{u}}=(1,~{}0)$. This may seem like an oversimplification, but it has no discernible impact on the shape of the resolvent modes in the region where the experimental mean is known. The second profile is the data-assimilated mean, whose results are compared to those of the interpolated experimental mean in order to reinforce the benefits of data-assimilation before performing resolvent analysis. This section is divided into three parts: the results for the A0 and A10 cases are discussed in §6.1 and §6.2, respectively. In §6.3, the results from resolvent analysis are compared to those from SPOD and we remark on the different behaviours observed for an oscillator versus an amplifier flow.

6.1 A0 case: a single linear mechanism

The first three singular values of the resolvent operator are plotted in figure 14 using both the interpolated experimental and data-assimilated means. The trends for both cases are quite similar to cylinder flow (Symon et al., 2018) in that there is a distinct range of frequencies where the first singular value is an order of magnitude higher than the second singular value. Data-assimilating the mean profile has a negligible impact on the absolute and relative magnitudes of the singular values with respect to one another. There is also no impact on the mode shapes (see Symon, 2018).

6.2 A10 case: multiple linear mechanisms

Unlike the A0 case, the singular values, illustrated in figure 15, of the resolvent operator for the A10 case depend strongly on the mean profile. The interpolated mean has a single peak at $\omega=5.9$ whereas the data-assimilated mean has two peaks at $\omega=7.3$ and $\omega=15.4$. The frequency range over which the resolvent is low-rank is also much wider for the data-assimilated mean.

The identification of a secondary peak can be attributed to filling in data in the shear layer near the leading edge. As seen in figure 16, which contains the most amplified resolvent mode at two frequencies corresponding to the peaks in figure 15 as well as one intermediate frequency ($\omega=12.0$), two separate linear mechanisms are identified. The first corresponds to shedding- or wake-type modes which occur at lower temporal frequencies (top two rows of figure 16). The second corresponds to a Kelvin-Helmholtz instability mechanism in the shear layer which occurs at higher temporal frequencies (bottom row of figure 16.) The absence of PIV vectors near the leading edge results in resolvent analysis failing to capture the second linear mechanism when the interpolated experimental mean is analysed instead of the data-assimilated mean. The mode shapes for intermediate frequencies have a mix of both linear mechanisms as seen for $\omega=12$ in figure 16. These results are consistent with Thomareis & Papadakis (2018) and Yeh & Taira (2019), who computed resolvent modes for a NACA 0012 airfoil at similar Reynolds numbers and angles of attack from time-averaged DNS and large-eddy simulation (LES) data, respectively. The latter study also identified two branches in the eigenspectrum of the LNS operator which were grouped into wake or shear layer modes.

6.3 Comparison to SPOD

Even though data-assimilation adds no new insight for the A0 case, there are several noteworthy differences between the results of resolvent analysis and SPOD for this flow. The most significant one is that resolvent analysis identifies a single peak at the shedding frequency and low-rank activity is limited to nearby frequencies only. There is no low-rank behaviour at higher harmonics whereas in SPOD, there is an order of magnitude separation between the first and second eigenvalue of the first three harmonics.

These observations can be explained by accounting for the nature of the nonlinear fluctuations. As in Towne et al. (2018), let us rewrite the cross-spectral density at an arbitrary frequency $\omega$ in terms of an expectation operator $\mathbb{E}(\cdot)$:

[TABLE]

Noting that $\hat{\boldsymbol{u}}(\omega)=\mathcal{H}(\omega)\hat{\boldsymbol{f}}(\omega)$, (24) can be rewritten as

[TABLE]

Instead of expressing $\mathcal{H}(\omega)$ in terms of its singular value decomposition, the resolvent operators are moved outside the expectation operator since they are deterministic to yield the following equation

[TABLE]

where $\hat{\boldsymbol{S}}_{\hat{\boldsymbol{f}}\hat{\boldsymbol{f}}}(\omega)=\mathbb{E}\{\hat{\boldsymbol{f}}(\boldsymbol{x},\omega)\hat{\boldsymbol{f}}^{*}(\boldsymbol{x},\omega)\}$ is the cross-spectral density of the nonlinear forcing as in Towne et al. (2019). $\hat{\boldsymbol{S}}_{\hat{\boldsymbol{u}}\hat{\boldsymbol{u}}}(\omega)$ may be low-rank in two scenarios.

The first is when the resolvent operator is low-rank, which occurs when there is a linear amplification mechanism. This is observed at the shedding frequency of the A0 case, frequencies corresponding to wake modes in the A10 case, and frequencies of the shear layer modes in the A10 case as summarised in rows 1, 4, and 5 in table 1. Both SPOD and resolvent analysis identify low-rank behaviour, similar to the case of the turbulent jet in Schmidt et al. (2018) where the Kelvin-Helmholtz instability accounts for the low-rank nature of the cross-spectral density tensor.

The second source of low-rank behavior can originate from the cross-spectral density of the nonlinear forcing, or $\hat{\boldsymbol{S}}_{\hat{\boldsymbol{f}}\hat{\boldsymbol{f}}}(\omega)$, being low-rank. Even if $\mathcal{H}(\omega)$ is nearly full-rank, as it is for higher harmonics of the shedding frequency, the left-hand side of (26) will be low-rank as long as at least one of the terms on the right-hand side is low-rank, i.e.

[TABLE]

For the A0 case, the nonlinear forcing is very structured at the shedding frequency and its harmonics; consequently, SPOD of the velocity fluctuations reveals large spectral gaps at these frequencies as seen in rows 1 and 2 of table 1. We theorise that the nonlinear forcing can be estimated by limited triadic interactions of resolvent modes to educe the structure at these higher harmonics, a hypothesis which is tested in the next section.

The SPOD eigenspectra of the nonlinear forcings in figures 12 and 13 verify that the cross-spectral density is low-rank at the shedding frequency and its harmonics in the A0 case and at no frequencies for the A10 case. The trends reinforce the notion that A0 behaves like an oscillator while A10 behaves like an amplifier. It also suggests that resolvent modes and SPOD modes will deviate substantially for the higher harmonics in the A0 case since the resolvent is not low-rank and the forcing is anything but white noise in space and time. This will be seen in §7.1 for the second harmonic of the shedding frequency. The modes for the A10 case, on the other hand, can be expected to match at all frequencies where the the resolvent is low-rank since these coincide with the low-rank frequencies identified by SPOD.

We compare the dominant SPOD mode to the optimal resolvent mode in figure 17 to reinforce the notion that the first resolvent response mode is representative of the true fluctuations in the flow when the resolvent is low-rank. One difference, which is particularly notable for the A0 shedding mode, is that the amplitude of the resolvent mode increases as a function of $x$. By comparing the final red contour in the top row of figure 17, we note that the amplitude of the SPOD mode is quite weak at $x=2.8$, whereas it is very strong for the resolvent mode. A similar phenomenon can be observed for the A10 wake mode. The shear layer mode, on the other hand, is a very close match. Despite these discrepancies, the wavelengths of the modes are in complete agreement, signifying that the resolvent modes correctly predict the convection velocities of the structures in the flow.

7 Resolvent-based modelling of the velocity fluctuations

We now move on to determine how much of the full fluctuating field can be deduced starting from only a limited mean profile and data at a few probe points. We search for an efficient basis for a reduced-order model of the energetic coherent structures using the data-assimilated mean and very limited time-resolved probes (i.e. knowledge of the velocity at one or two points). This section is divided into three parts: in §7.1 we attempt to approximate the structure of the nonlinear forcing using limited triadic interactions at frequencies where the linear dynamics are not low-rank for the A0 case. The reduced-order models with and without these parasitic modes are compared. In §7.2, we consider the A10 case and examine whether just one or multiple probes are needed to account for the presence of two linear mechanisms in the flow. Finally, we discuss our results in the context of experimental limitations encountered in §7.3.

7.1 Parasitic modes and modelling of the A0 case

7.1.1 Approximation of the nonlinear forcing

To improve the basis for frequencies outside the range where the resolvent operator is low-rank, the nonlinear forcing is approximated using very limited triadic interactions. Let us consider two examples. The simplest is the second harmonic of the shedding frequency $\omega_{s}$, whose nonlinear forcing is approximated in exactly the same manner as that outlined in Rosenberg et al. (2019), i.e.

[TABLE]

Many other modes may contribute to the right-hand side of (28) although this approximation is adequate since the singular values associated with $\hat{\boldsymbol{\psi}}_{1}(\omega_{s})$ are very large. Very low frequencies can also be recovered by postulating that they result from the beating of two highly amplified frequencies. The nonlinear forcing for the frequency $\omega=3.70$, for example, is

[TABLE]

Since the singular value associated with the shedding frequency at $\omega_{s}=12.0$ is very large, it can be reasoned that the magnitude of $\chi_{1}$ for this frequency will outweigh all other terms in (29). To obtain the correct frequency, therefore, it only needs to interact with the frequency $\omega=-8.30$, where the resolvent operator is still approximately rank-1, to produce a mode at the desired frequency $\omega=+3.70$. (29), therefore, reduces to

[TABLE]

where the $\chi_{1}$’s have been dropped since $\chi_{1}(\omega=12.0)\chi_{1}(\omega=-8.3)$ appears before each term in (30) and, at this point, we are only interested in improving the basis for this particular frequency. If we were to include more triads in either (28) or (30), we would need to know the phase of the modes since $\chi_{1}$ appears before each mode in (29). We argue here that the high amplifications associated with the shedding mode allow us to ignore the triadic interactions of suboptimal modes and frequencies which are less amplified.

The parasitic modes, or ‘forced’ resolvent modes (i.e. those obtained after right-multiplying the resolvent operator by the approximated nonlinear forcing), are compared to the optimal response mode $\hat{\boldsymbol{\psi}}_{1}$ from the singular value decomposition and the first SPOD mode in figure 18.

Despite how simple the modelling for the nonlinear forcing is, the structure at all frequencies considered is a major improvement over the predicted structure from the first resolvent response mode. In particular, the parasitic mode obeys the correct symmetries while the optimal response mode resembles a stretched or compressed version of the vortex shedding mode.

7.1.2 Modelling preliminaries

Several candidate points to calibrate the complex weights of the resolvent modes are displayed in figure 19. They are referred to as points 1, 2, and 3, and are ordered by their streamwise position, i.e. point 1 is the furthest upstream.

The power spectrum of P2 (P1 and P3 are very similar) for the transverse velocity component, plotted in figure 19, contains a large peak at $\omega=12.0$ and a secondary one at $\omega=24.0$. These are the first and second harmonics of the shedding frequency, respectively. There is another discernible peak at $\omega=10.7$ resulting from jitter in the shedding frequency. This is also present in the SPOD results for the leading eigenvalue, although it is less noticeable. The jitter does not influence the resolvent norm.

The unweighted resolvent modes are first normalised using the standard energy inner product based on the two velocity components. Once the resolvent modes are weighted using the calibration point, they should reflect the true magnitude of the fluctuations. To make a fair comparison between the model and snapshots obtained from PIV, the experimental data are filtered to exclude temporal frequencies outside the range used in the model. The filtering involves Fourier-transforming the raw PIV data in time, truncating Fourier modes at unwanted frequencies, and inverse Fourier-transforming back into the time domain. Another possible way to compare the resolvent model to the experiment is to project the raw data onto the first SPOD mode at all frequencies of interest. Generally these would be where SPOD is approximately rank-1. Since the frequency resolution of the SPOD modes is only $\Delta\omega\approx 0.5$, we choose to filter based on the first procedure.

7.1.3 Results and quantification of error

We begin by comparing the performance of each candidate point. The instantaneous error is computed between the model and PIV as follows

[TABLE]

where $v$ denotes the transverse velocity component and $N$ is the number of PIV vectors. $I_{v}(t)$ is integrated in time to obtain a global error metric

[TABLE]

$E_{v}$ is normalised by the total kinetic energy contained in the filtered PIV and the results for each point are summarised in table 2. Even though the performance of each point is good at some time instants, the overall performance for all three is relatively poor although P3 is the best. If two points are considered, the performance is nearly an order of magnitude better. Furthermore, points which are situated further downstream in the flow yield better results, which is a consequence of the resolvent modes failing to capture the correct streamwise decay of the fluctuations. The combination of points P2 and P3 yields the lowest amount of error, so the following discussion considers this choice only.

The values computed in table 2 factored in all resolvent modes computed, i.e. both classical and parasitic. The first modelling attempt, however, only uses modes where the resolvent operator is nearly rank-1, i.e. $6.73<\omega<15.7$, henceforth referred to as rank-1 modes. These 41 modes are computed using the first singular response mode of the resolvent operator. Harmonics of the shedding frequency where SPOD reveals low-rank behaviour are excluded for now. A representative snapshot is presented in figure 20 and the flow resembles the traditional von Kármán vortex street behind a cylinder. The raw PIV data are contained in the top row while the filtered PIV are compared to the resolvent model using rank-1 modes only in the second row.

The velocity fluctuations are symmetric with respect to the centreline since the most amplified structure predicted by resolvent analysis resembles the von Kármán vortex street at all rank-1 frequencies. Similar to the cylinder, there is only one linear mechanism, so modes at suboptimal frequencies end up being a stretched or compressed (in the streamwise direction) version of the shedding mode (see Symon, 2018; Symon et al., 2019). The true velocity fluctuations are not symmetric with respect to the centreline if all frequencies are considered. The difference between the filtered PIV and resolvent model is plotted in figure 21. The largest discrepancies occur immediately behind the trailing edge ($1.0<x<1.5$), where the fluctuations are under predicted by the resolvent model. By calibrating at points P2 and P3, where the resolvent modes over predict the fluctuations, the resulting model has a harder time reconstructing the region further upstream.

When parasitic modes are included, the model is no longer constrained to be symmetric since the second harmonic and low frequency modes are antisymmetric with respect to the centreline. The flow is now modelled using 91 modes ranging from frequencies as low as $\omega=2.72$ to frequencies as high as $\omega=39$, in order to encompass the third harmonic. The resolvent-model is animated in time and compared to the filtered PIV in movies 1 and 2, respectively. Not all frequencies within this range are included since many of them are not energetic as indicated by the probe points. The model is compared to the filtered PIV in the bottom row of figure 20.

Because the modes outside the rank-1 range are not particularly energetic, the filtered PIV in the bottom row resembles the raw PIV snapshot in the top row fairly well (it contains 80% of the kinetic energy). The addition of parasitic modes does not dramatically improve the model, but it does capture more detailed features of the flow field. Further downstream, the negative velocity perturbations tilt upstream while the positive velocity perturbations tilt downstream and this is consistent with the PIV data. The discrepancy between the improved model and the filtered PIV is illustrated in figure 21. Small improvements can be seen in the region beyond $x=1.5$. For this particular snapshot, the addition of parasitic modes reduces the error from 22% to 18%.

7.2 Two linear mechanisms and modelling of the A10 case

7.2.1 Selection of probe points

The challenge associated with the A0 case was to incorporate parasitic modes using nonlinear interactions. For the A10 case, the difficulty stems from being able to capture the two linear mechanisms present in the flow. Two candidate points are proposed in figure 22. P4 is located in the shear layer and P5 is in the airfoil wake.

The power spectrum for P4 contains two broad peaks at frequencies similar to those identified by resolvent analysis and SPOD. This strengthens the argument that the resolvent tends to be low-rank at energetic frequencies as remarked by Moarref et al. (2013) for wall-bounded turbulent flows. The power spectrum for P5, on the other hand, contains one broad peak centred around the wake frequencies. It seems reasonable, therefore, to suggest that calibrating the modes using this point alone would not successfully model the fluctuations in the shear layer.

The performance using a single point or both points is summarised in table 3 using 66 modes in the range $3.14<\omega<17.7$. The reconstruction does very badly when only P4 is used as it overpredicts the fluctuations in the wake substantially. P5 is an improvement but it underpredicts the fluctuations in the shear layer. The model performs best when both points are used to calibrate the resolvent modes.

7.2.2 Modelling of the A10 velocity fluctuations

For the sake of brevity, only results for the use of both points are summarised in figure 23 as they help explain why a single point is not sufficient to obtain a good model of the flow. The raw PIV data are visualised in the top row. It is fairly evident that the flow is more complicated for this angle of attack as the presence of both linear mechanisms can be observed in relatively distinct regions of the flow. There is also an area immediately behind the trailing edge where they overlap leading to more complex flow structures.

Similar to the A0 case, the PIV data are frequency-filtered. We begin by considering 31 frequencies where roughly wake modes only are present, i.e. $3.15<\omega<9.65$, as seen in row 2 of figure 23. Even though there is some activity in the shear layer, it is small in comparison to the unfiltered case. The resolvent model is plotted alongside it. While the agreement between the wavelengths of the structures is excellent, the model overpredicts the intensity of the fluctuations as far upstream as $x=1.5$. One potential method for combating this problem is to consider a discounted resolvent, i.e.

[TABLE]

where $s$ is a complex number with nonzero real part. Yeh & Taira (2019) applied the discounted resolvent to flow around a NACA 0012 airfoil since they encountered unstable modes. Since we do not find unstable eigenmodes near the wake frequencies, we opt to not apply this procedure.

The PIV data are next filtered to retain 31 frequencies corresponding to the shear layer dynamics, i.e. $11.0<\omega<17.7$, and are compared to the resolvent model in row 3 of figure 23. The model is very good, both in terms of the wavelengths of the structures predicted and their amplitudes. The signature of these modes at P5 is almost negligible, however, so the model misses them altogether when only P5 is considered. This explains why the use of two points fares significantly better than just one.

Finally, the model is compared to the filtered PIV for $3.14<\omega<17.7$ (all 66 frequencies) in the bottom row of figure 23. Consistent with earlier results, the agreement is best in the shear layer while the performance of the model degrades further downstream. The model does capture some of the complex behaviour immediately behind the airfoil and the inclination of the fluctuations is correct overall. Despite poor performance further downstream, potential reasons for which are discussed in the next subsection, the resolvent model does quite well at reproducing the coherent structures of the flow. While it is beyond the scope of this study, the very good agreement in the shear layer for the A10 case suggests that the fluctuating forces could be predicted by the resolvent modes. Furthermore, these could be tested by measuring the fluctuating forces on the airfoil with a force sensor.

7.2.3 Dynamics and convection velocities of the structures

Before discussing experimental limitations of the procedure, it is worth analysing the convection velocities of the structures in the resolvent model. To aid the reader, both the filtered PIV fluctuations ($3.14<\omega<17.7$) and resolvent model are animated and uploaded as supplementary movies 3 and 4, respectively, to visualise the evolution of the flow in time. The movies contain 1000 snapshots and show that there is a coupling between the shear layer and wake modes in the region of the flow where the recirculation bubble is situated. Despite jitter in the flow’s behaviour, the vortices which form in the shear layer meet those which emerge from the trailing edge and, in some cases, they merge. The convection velocity of the fluctuations, consequently, is relatively constant as a function of $x$ although the perturbations begin to slow down further downstream. The reduced-order model from resolvent analysis is able to capture these key features despite poor predictions beyond $x=2$.

7.3 Experimental limitations

Before we conclude, it is worth discussing how the results could be improved and why this is difficult to achieve due to the limitations of the current experimental data sets. Based on the SPOD results in this study and previous literature, it is likely that A0 and A10 have separate solutions.

The main problem associated with the A0 case is that the resolvent modes do not correctly model the streamwise decay of the fluctuations, a phenomenon also remarked upon by Rosenberg et al. (2019) for cylinder flow. In order to combat this problem, Rosenberg et al. (2019) approximated the nonlinear forcing for the shedding frequency using the triadic interaction between the shedding mode, predicted by the singular value decomposition, and the parasitic mode at the second harmonic, i.e.

[TABLE]

We attempted to do the same for the A0 case, since we have a good prediction of the fluctuation at the shedding frequency from the optimal resolvent mode and the second harmonic by approximating its nonlinear forcing. SPOD of the nonlinear forcing, furthermore, strongly suggests that it is necessary to account for this feedback since there is separation between the first two SPOD eigenvalues and the first mode is very structured. The approximated nonlinear forcing for the airfoil at the shedding frequency, however, did not improve the mode shape.

A possible reason for why this failed might be related to the enormous amplification at the shedding frequency due to the proximity of a stable eigenvalue to the neutral axis. Similar to cylinder flow, it is necessary to obtain a very accurate estimate of the nonlinear forcing in order to obtain the correct mode since the resolvent operator is so biased in favour of the first response mode. Lesshafft et al. (2018) also remarked that spurious numerical noise can be sufficient to trigger the optimal mode and thus overwhelm suboptimal modes. The data-assimilated mean profile is not completely free of noise so the resolvent modes and their triadic interactions are not sufficiently accurate to alter the mode shape predicted by the resolvent at the shedding frequency. It would be worth investigating low Reynolds number DNS data, which are uncontaminated by noise, to understand the exact source of the deviations and how to alleviate these problems when using experimental data.

The main region of disagreement for the A10 case is in the far wake, where the model overestimates the velocity fluctuations. Since the nonlinear forcing is not structured at this angle of attack, it suggests the root cause to be three-dimensional effects which the model does not incorporate. Three-dimensionality results from the formation of stall cells, the number of which is given by $n_{s}=AR/2.28$ (Weihs & Katz, 1983), resulting in $n_{s}\approx 2$ for this setup. We observed earlier that the mean spanwise velocity component and its spanwise gradient are approximately zero. The mean flow, consequently, is nearly two-dimensional, hence the success of the data-assimilation (particularly with respect to the captured Reynolds stress gradients). In order to consider spanwise wavenumbers, however, we would need to obtain the spanwise-averaged mean profile (similar to Beneddine et al., 2016), which is inaccessible from planar PIV. By considering spanwise wavenumbers other than $k_{z}=0$, which is more suitable for the region of the flow dominated by the Kelvin-Helmholtz instability, it may be possible to obtain resolvent modes which more accurately depict the far wake dynamics. The investigation of these proposed solutions is an area of current study.

8 Conclusions

The flows around an airfoil at two angles of attack, $\alpha=0^{\circ}$ and $\alpha=10^{\circ}$, have been modelled from resolvent analysis using the procedure summarised in figure 2. The experimental measurements which were necessary included a crude mean profile obtained by PIV and no more than two probe measurements. The experimental data were first analysed using SPOD, which suggested that the angle of attack resulted in significant changes in the flow. SPOD of the velocity fluctuations indicated low-rank behaviour at the shedding frequency and its harmonics for the A0 case. The first eigenvalue of the cross-spectral density tensor was significantly higher than the second eigenvalue, resulting in large but narrow peaks. The A10 case, on the other hand, did not contain well-defined peaks although there was separation between the first two eigenvalues around two ranges of frequencies. The mode shapes within these ranges resembled wake and shear layer modes, respectively.

We also used SPOD to analyse the nonlinear fluctuations, which were directly computed from the PIV data. In the A0 case, the cross-spectral density of the nonlinear forcing exhibited similar trends to the velocity fluctuations in that a spectral gap existed at the shedding frequency and its harmonics. No such similarity existed for the A10 case as there was no low-rank behaviour at any frequency. The mode shapes observed for the A10 case, consequently, were highly unstructured, suggesting they did not emerge from a small number of modal interactions. As expected, the mode shapes for the A0 case were very structured for the shedding frequency and its harmonics, suggesting that the flow behaves like an oscillator with intrinsic dynamics which are insensitive to background noise. It followed that the A10 case behaved more like an amplifier.

In order to confirm this interpretation, the flows were analysed from a resolvent point of view. Before using the mean profile as an input, they were data-assimilated to decrease noise, fill in missing data, and correct artificial errors which resulted from stitching together multiple mean profiles (in the A10 case). The effect of data-assimilation was far more evident for the A10 case as it managed to identify two linear mechanisms which peaked at frequencies very similar to those identified by SPOD. When the mean profile was interpolated onto the CFD mesh, it was unable to identify the shear layer peak due to missing data near the separation point. Both the data-assimilated and interpolated means only identified a single linear mechanism for the A0 case at the shedding frequency, confirming our earlier interpretation that the A0 case behaved like an oscillator which revolved around the shedding mode while the A10 case behaved like an amplifier.

By expanding the cross-spectral density tensor into the product of the resolvent operator, its Hermitian transpose, and the cross-spectral density tensor of the nonlinear forcing (similar to Towne et al., 2018), we noted that low-rank behaviour may emerge when either the resolvent is low-rank (a linear amplification mechanism) or the cross-spectral density tensor of the nonlinear fluctuations is low-rank (structured nonlinear forcing due to modal interactions). The instances of low-rank behaviour due to either cause were summarised in table 1 for both airfoils. Notably in the A0 case, we educed the structure of the fluctuations at higher harmonics, where the resolvent operator is not low-rank, by considering limited triadic interactions of highly amplified resolvent modes.

When it came to resolvent-based modelling of the A0 flow, we no longer had to rely exclusively on the resolvent operator being low-rank to predict fluctuations. Indeed when the model was limited to rank-1 modes, we only recovered fluctuations which are symmetric with respect to the centreline since all rank-1 modes resemble the main instability mechanism. By adding parasitic modes, we were able to recover asymmetric fluctuations and capture finer details of the velocity fluctuations. Parasitic modes were not featured in the reconstruction of the A10 case. Instead, we determined that two probe points were necessary to capture both linear mechanisms in the flow and prevent overestimation of the fluctuations in the far wake region.

The combination of data-assimilation and resolvent analysis have shown their potential to fill in data from incomplete measurements and it is for this reason that they are a topic of ongoing study. The governing equations were exploited to obtain a mean velocity field which obeyed prescribed boundary conditions, removed experimental noise, and increased resolution. The performance of interpolation is clearly inferior as seen in the A10 case. Moreover, resolvent analysis produces a good basis which can be calibrated from limited time-resolved measurements.

It can also be concluded that the combination of SPOD and resolvent analysis can shed light on the nature of the nonlinear forcing in complex flows. In particular, they highlight when the nonlinear forcing is structured and the computation of parasitic modes in this study from experimental data suggest that this process can be leveraged for reduced-order modelling in the future. In cases where the nonlinear forcing is less structured, e.g. the A10 case, it could be possible to devise new strategies for predicting their statistics, e.g. Towne et al. (2019), or modelling their effect on the large-scale structures, e.g. Illingworth et al. (2018). It may even suffice to rely on linear mechanisms only and treat the forcing as stochastic. The results from the paper also confirm the potential of resolvent analysis to predict acoustic fields, e.g. Schmidt et al. (2018), or the forces acting on a solid body, e.g. Gómez et al. (2016b). The A10 model, in particular, is quite promising in the vicinity of the airfoil so it would be possible to verify whether the predicted pressure field is valid by comparing the forces to those measured experimentally using a force sensor.

9 Acknowledgments

The support of the Army Research Office (ARO) through grant number W911NF-17-1-0306 and the Office of Naval Research (ONR) through grant number N00014-17-1-3022 is gratefully acknowledged.