The role of parasitic modes in nonlinear closure via the resolvent feedback loop

TL;DR

This paper investigates how parasitic modes influence nonlinear closure in fluid flows by analyzing the feedback loop of the resolvent operator, revealing their role in nonlinear forcing and flow fluctuation structures.

Contribution

It introduces a method to identify the influence of parasitic modes within the resolvent framework, enhancing understanding of nonlinear flow interactions beyond dominant linear mechanisms.

Findings

Parasitic modes contribute significantly to nonlinear forcing.

Resolvent modes align with DMD modes at fundamental frequencies.

Discrepancies occur at higher harmonics due to lack of low-rank structure.

Abstract

We use the feedback loop of McKeon \& Sharma (2010), where the nonlinear term in the Navier-Stokes equations is treated as an intrinsic forcing of the linear resolvent operator, to educe the structure of fluctuations in the range of scales (wavenumbers) where linear mechanisms are not active. In this region, the absence of dominant linear mechanisms is reflected in the lack of low-rank characteristics of the resolvent and in the disagreement between the structure of resolvent modes and actual flow features. To demonstrate the procedure, we choose low Reynolds number cylinder flow and the Couette equilibrium solution EQ1, which are representative of very low-rank flows dominated by one linear mechanism. The former is evolving in time, allowing us to compare resolvent modes with Dynamic Mode Decomposition (DMD) modes at the first and second harmonics of the shedding frequency. There is a match between the modes at the first harmonic but not at the second harmonic where there is no separation of the resolvent operator's singular values. We compute the self-interaction of the resolvent mode at the shedding frequency and illustrate its similarity to the nonlinear forcing of the second harmonic. When it is run through the resolvent operator, the `forced' resolvent mode shows better agreement with the DMD mode. A similar phenomenon is observed for the fundamental streamwise wavenumber of the EQ1 solution and its second harmonic. The importance of parasitic modes, labeled as such since they are driven by the amplified frequencies, is their contribution to the nonlinear forcing of the main amplification mechanisms as shown for the shedding mode, which has subtle discrepancies with its DMD counterpart.

Equations

$${\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 }},\nabla \cdot {\mathbf{u}}^{\prime }=0.$$

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

$$\mathcal{H}(\mathbf{k})=\sum _p{\hat{\mathbf{\psi }}}_p(\mathbf{k}){\sigma }_p(\mathbf{k}){\hat{\mathbf{\varphi }}}_p^{\ast }(\omega ),$$

$$\hat{\mathbf{u}}(\mathbf{k})\approx {\hat{\mathbf{\psi }}}_1(\mathbf{k}){\sigma }_1(\mathbf{k}){\chi }_1(\mathbf{k})={\hat{\mathbf{\psi }}}_1(\mathbf{k}){\sigma }_1(\mathbf{k})⟨{\hat{\mathbf{\varphi }}}_1(\mathbf{k}),\hat{\mathbf{f}}(\mathbf{k})⟩.$$

$$\hat{\mathbf{f}}({\mathbf{k}}_3)=\sum _{{\mathbf{k}}_1+{\mathbf{k}}_2={\mathbf{k}}_3}-\hat{\mathbf{u}}({\mathbf{k}}_1)\cdot \nabla \hat{\mathbf{u}}({\mathbf{k}}_2)$$

$$\hat{\mathbf{f}}({\mathbf{k}}_3)=\sum _{{\mathbf{k}}_1+{\mathbf{k}}_2={\mathbf{k}}_3}\sum _a\sum _b-{\sigma }_a({\mathbf{k}}_1){\sigma }_b({\mathbf{k}}_2){\chi }_a({\mathbf{k}}_1){\chi }_b({\mathbf{k}}_2)\left[{\hat{\mathbf{\psi }}}_a({\mathbf{k}}_1)\cdot \nabla {\hat{\mathbf{\psi }}}_b({\mathbf{k}}_2)\right].$$

$$iq{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}\hat{\boldsymbol{u}}(q\boldsymbol{k}_{h})=\boldsymbol{L}\hat{\boldsymbol{u}}(q\boldsymbol{k}_{h})+\mathcal{N}(q\boldsymbol{k}_{h}),$$

$$\mathcal{N}(q{\mathbf{k}}_h)=\sum _{r\ne 0,q}\mathcal{N}(\hat{\mathbf{u}}(r{\mathbf{k}}_h),\hat{\mathbf{u}}((q-r){\mathbf{k}}_h).$$

$$\underbrace{2i{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}\hat{\boldsymbol{u}}(2\boldsymbol{k}_{h})}_{\epsilon^{2}}=\underbrace{\boldsymbol{L}\hat{\boldsymbol{u}}(2\boldsymbol{k}_{h})}_{\epsilon^{2}}+\underbrace{\mathcal{N}(\hat{\boldsymbol{u}}(\boldsymbol{k}_{h}),\hat{\boldsymbol{u}}(\boldsymbol{k}_{h}))}_{\epsilon^{2}}+\underbrace{\mathcal{N}(\hat{\boldsymbol{u}}(3\boldsymbol{k}_{h}),\hat{\boldsymbol{u}}(-\boldsymbol{k}_{h}))}_{\epsilon^{4}}+\underbrace{\mathcal{N}(\hat{\boldsymbol{u}}(-\boldsymbol{k}_{h}),\hat{\boldsymbol{u}}(3\boldsymbol{k}_{h}))}_{\epsilon^{4}}+\cdots.$$

$$\hat{\mathbf{u}}(2{\mathbf{k}}_h)\approx \mathcal{H}(2{\mathbf{k}}_h)\mathcal{N}(\hat{\mathbf{u}}({\mathbf{k}}_h),\hat{\mathbf{u}}({\mathbf{k}}_h))\approx \hat{\mathbf{u}}(2{\mathbf{k}}_h)\approx -{\sigma }_1^2{\chi }_1^2\mathcal{H}(2{\mathbf{k}}_h)({\hat{\mathbf{\psi }}}_1({\mathbf{k}}_h)\cdot \nabla {\hat{\mathbf{\psi }}}_1({\mathbf{k}}_h)),$$

$$\underbrace{i{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}\hat{\boldsymbol{u}}(\boldsymbol{k}_{h})}_{\epsilon}=\underbrace{\boldsymbol{L}\hat{\boldsymbol{u}}(\boldsymbol{k}_{h})}_{\epsilon}+\underbrace{\mathcal{N}(\hat{\boldsymbol{u}}(2\boldsymbol{k}_{h}),\hat{\boldsymbol{u}}(-\boldsymbol{k}_{h}))}_{\epsilon^{3}}+\underbrace{\mathcal{N}(\hat{\boldsymbol{u}}(-\boldsymbol{k}_{h}),\hat{\boldsymbol{u}}(2\boldsymbol{k}_{h}))}_{\epsilon^{3}}+\cdots.$$

$$\begin{split}\hat{\boldsymbol{f}}(\boldsymbol{k}_{h})&=\chi_{1}|\chi_{1}|^{2}\left\{\left[\sigma_{1}^{2}\mathcal{H}(2\boldsymbol{k}_{h})(\hat{\boldsymbol{\psi}}_{1}\cdot\nabla\hat{\boldsymbol{\psi}}_{1})\right]\cdot\nabla\left[\sigma_{1}\hat{\boldsymbol{\psi}}_{1}^{*}\right]+\left[\sigma_{1}\hat{\boldsymbol{\psi}}_{1}^{*}\right]\cdot\nabla\left[\sigma_{1}^{2}\mathcal{H}(2{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}})(\hat{\boldsymbol{\psi}}_{1}\cdot\nabla\hat{\boldsymbol{\psi}}_{1})\right]\right\}\\ &=\chi_{1}|\chi_{1}|^{2}\hat{{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{f}}}}^{\prime}(\boldsymbol{k}_{h}),\end{split}$$

$$|{\chi }_1|={\left[\frac{1}{⟨{\hat{\mathbf{\varphi }}}_1({\mathbf{k}}_h),{\hat{\mathbf{f}}}^{\prime }({\mathbf{k}}_h)⟩}\right]}^{1/2}.$$

Full text

The role of parasitic modes in nonlinear closure via the resolvent feedback loop

Kevin Rosenberg, Sean Symon, and Beverley J. McKeon

{addmargin}

[3em]3em

We use the feedback loop of McKeon & Sharma, (2010), where the nonlinear term in the Navier-Stokes equations is treated as an intrinsic forcing of the linear resolvent operator, to educe the structure of fluctuations in the range of scales (wavenumbers) where linear mechanisms are not active. In this region, the absence of dominant linear mechanisms is reflected in the lack of low-rank characteristics of the resolvent and in the disagreement between the structure of resolvent modes and actual flow features. To demonstrate the procedure, we choose low Reynolds number cylinder flow and the Couette equilibrium solution EQ1, which are representative of very low-rank flows dominated by one linear mechanism. The former is evolving in time, allowing us to compare resolvent modes with Dynamic Mode Decomposition (DMD) modes at the first and second harmonics of the shedding frequency. There is a match between the modes at the first harmonic but not at the second harmonic where there is no separation of the resolvent operator’s singular values. We compute the self-interaction of the resolvent mode at the shedding frequency and illustrate its similarity to the nonlinear forcing of the second harmonic. When it is run through the resolvent operator, the ‘forced’ resolvent mode shows better agreement with the DMD mode. A similar phenomenon is observed for the fundamental streamwise wavenumber of the EQ1 solution and its second harmonic. The importance of parasitic modes, labeled as such since they are driven by the amplified frequencies, is their contribution to the nonlinear forcing of the main amplification mechanisms as shown for the shedding mode, which has subtle discrepancies with its DMD counterpart.

1 Introduction

The resolvent framework of McKeon & Sharma, (2010) can be used for low-order modeling (e.g. Rosenberg,, 2018) and reconstruction of flows from limited measurements (e.g. Symon,, 2018; Gómez et al.,, 2016). Unlike data-driven methods such as dynamic mode decomposition DMD, which require snapshots of the fluctuating flow fields, only knowledge of the mean velocity profile is needed to form the resolvent operator which maps nonlinear terms, or forcing, to velocity fluctuations. One of the primary reasons for its success is the tendency for the resolvent operator to be low-rank at energetic frequencies. Under such conditions the optimal mode, which can be obtained via a singular value decomposition (SVD) of the operator, is significantly more amplified than its suboptimal counterparts and hence the flow structure can be approximated without precise knowledge of the nonlinear forcing. Effectively, then, mode structure is obtained by approximating the (action of the) resolvent operator.

Sharma et al., (2016) and Towne et al., (2018) have noted a correspondence between resolvent modes, DMD, and spectral Proper Orthogonal Decomposition (SPOD) modes. In particular, Towne et al., (2018) noted that resolvent modes are equivalent to SPOD modes when the nonlinear forcing is uncorrelated in space and time. When the nonlinear forcing is not white noise, resolvent modes deviate from their data-driven counterparts and this is particularly evident at frequencies where the amplification of the optimal mode is not sufficiently higher than the suboptimal modes.

The objective here is to take advantage of the feedback loop in McKeon et al., (2013) and determine under what conditions mode shapes can be obtained by the action of the full resolvent on nonlinear forcing computed from resolvent modes, i.e. when it is appropriate to approximate the forcing rather than the resolvent, and how well these modes compare from data-driven approximations. In this respect, this study will share similarities to weakly nonlinear analyses (Herbert,, 1983; Sipp & Lebedev,, 2007) in the sense of considering a limited set of finite-amplitude perturbations and their interactions. However, while weakly nonlinear analyses typically consider the neutral stability of a base flow and expansions in the vicinity of a critical Reynolds number, here we will consider the neutral stability of a mean flow and make no restriction on the Reynolds number range. The flows analyzed herein are cylinder flow and the EQ1 equilibrium solution of Couette flow. Both flows exhibit strong amplification at a single frequency/spatial wavenumber which is linked to the neutral stability of the temporal mean profile in the case of the cylinder (Barkley,, 2006) and the streamwise-averaged mean in the case of EQ1 (Hall & Sherwin,, 2010). Mantič-Lugo et al., (2014) leveraged this property to develop a self-consistent model for the cylinder where the Reynolds stresses are approximated by the most amplified mode whose amplitude is adjusted until the flow is neutrally stable. In a similar vein, the choice of EQ1 is motivated by previous work analyzing the emergence of streamwise-constant streaks, the scale interactions which give rise to them, and their connections to the self-sustaining process of turbulent shear flows (Waleffe,, 1997; Hall & Sherwin,, 2010; Farrell & Ioannou,, 2012).

2 Approach

We consider the incompressible Navier-Stokes equations (NSE) and decompose the velocity field into a mean component $\overline{\boldsymbol{u}}$ and fluctuating component $\boldsymbol{u}^{\prime}$. In the case of the cylinder, we use a temporal mean with $\overline{\boldsymbol{u}}=\left[\overline{u}(x,y),\overline{v}(x,y),0\right]^{T}$ and for EQ1 we use a temporal and streawise-averaged mean with $\overline{\boldsymbol{u}}=\left[\overline{u}(y,z),\overline{v}(y,z),\overline{w}(y,z)\right]^{T}$. The resulting equations for the fluctuations are

[TABLE]

Following McKeon & Sharma, (2010), we define $\boldsymbol{f}^{\prime}=-\boldsymbol{u}^{\prime}\cdot\nabla\boldsymbol{u}^{\prime}+\overline{\boldsymbol{u}^{\prime}\cdot\nabla\boldsymbol{u}^{\prime}}$ and treat $\boldsymbol{f}^{\prime}$ as an unknown forcing. Upon Fourier transforming in the appropriate directions, we can recast equation 1 as

[TABLE]

where $\mathcal{H}(\mathbf{k})$ is the resolvent operator for a particular wavenumber/frequency vector $\boldsymbol{k}$, $\boldsymbol{L}$ is the linear Navier-Stokes operator associated with the mean, and $\boldsymbol{B},\boldsymbol{C}$ are operators which define the inputs and outputs; a complete mathematical definition of these operators for these flows is found in Symon, (2018), Rosenberg, (2018). It should be noted for the cylinder flow that $\boldsymbol{k}$ is simply defined by the temporal frequency $\omega$ and in the case of the EQ1 by the streamwise wavenumber $k_{x}$ (and temporal frequency $\omega=0$, as it is an equilibrium solution). Subsequently, we will analyze the resolvent operators $\mathcal{H}(\omega)$ and $\mathcal{H}(k_{x})$ for the two flows respectively (where we have dropped the $\omega$ dependence in the latter). However, for generality and to facilitate comparisons between the flows going forward, we will maintain the notation $\mathcal{H}(\mathbf{k})$.

The resolvent operator relies on the mean velocity profile $\overline{\boldsymbol{u}}$ as an input and acts as a transfer function from nonlinearity $\hat{\boldsymbol{f}}(\mathbf{k})$ to velocity fluctuations $\hat{\boldsymbol{u}}(\mathbf{k})$ in Fourier space. A singular value decomposition (SVD) of the resolvent operator

[TABLE]

leads to optimal sets of orthonormal singular response and forcing modes, or $\hat{\boldsymbol{\psi}}_{p}(\mathbf{k})$ and $\hat{\boldsymbol{\phi}}^{*}_{p}(\mathbf{k})$, respectively. The modes are ranked by the forcing-to-response gain $\sigma_{p}(\mathbf{k})$ under the $L_{2}$ norm. When the first singular value is much larger than the second one, i.e., $\sigma_{1}(\mathbf{k})\gg\sigma_{2}(\mathbf{k})$, the velocity response may be approximated as

[TABLE]

In cases where the singular values are not sufficiently separated and/or the forcing is highly structured, Equation 4 will not accurately reproduce the velocity fluctuations at a given $\mathbf{k}$.

One course of action is to consider suboptimal resolvent modes ($p>1$) and project the resolvent forcing modes onto the nonlinear forcing $\hat{\boldsymbol{f}}(\mathbf{k})$ to determine the complex weights $\chi_{p}$ associated with each resolvent response mode. Here we approximate the structure of the fluctuations when the resolvent operator is not low-rank by considering the feedback loop introduced by McKeon et al., (2013), where the nonlinearity is formed by triadic interactions of the velocity modes, i.e.,

[TABLE]

The associated nonlinear forcing, consequently, is given by triadically consistent resolvent response modes obtained from the previous singular value decompositions

[TABLE]

All terms on the right-hand side of Equation 6 can be determined by approximation of the resolvent with the exception of $\chi_{a}$ and $\chi_{b}$, which may be obtained only under special circumstances, which we investigate here. In such cases, the forcing from Equation 6 can then be used to obtain better approximations of the velocity fluctuations for the non-low rank case.

We compare the mode shapes arising from a rank-one approximation of the resolvent operator with those from approximation of the forcing for two flows. Both have a strong physical mechanism underlying a dominant coherent structure in a narrow band of scales even though energetic activity is not limited to that band.

A direct numerical simulation (DNS) of two-dimensional cylinder flow is performed in FreeFem++ (Hecht,, 2012) using the same mesh, boundary conditions, and discretizations as those in Symon et al., (2018) at $Re=UD/\nu=100$, based on cylinder diameter $D$, inlet velocity $U$, and kinematic viscosity $\nu$. The time-step is $\Delta t=0.02$ and the mean is obtained after time-averaging over 25 vortex shedding cycles once the flow achieves a steady limit cycle. The linear operators needed for resolvent analysis are formed in FreeFem++ around the two-dimensional mean flow, and resolvent modes are computed using the procedure outlined in Sipp & Marquet, (2013). The DNS snapshots are analyzed using DMD (e.g Schmid,, 2010; Rowley et al.,, 2009) to isolate the true structure of the flow at individual frequencies.

The Couette equilibrium solution EQ1 (Nagata,, 1990; Gibson et al.,, 2009) is considered under the same domain and discretization as outlined in Gibson et al., (2009) with Reynolds number of $Re=1000$. The discrete spatial operators needed to form the resolvent utilize Chebyshev and Fouier differentation matrices in the wall-normal and spanwise directions respectively (Weideman & Reddy,, 2000).

3 Approximation of Cylinder Flow & Exact Coherent State

Let us denote $\boldsymbol{k}_{h}$ as the frequency/wavenumber of the first energetic harmonic for these two flows, i.e. the shedding frequeny $\omega_{s}$ and the fundamental streamwise wavenumber $k_{x_{f}}$ respectively. Symon et al., (2018) showed that low-rank behavior for cylinder flow at $Re=100$ is limited to a bandwidth of frequencies in the vicinity of $\boldsymbol{k}_{h}$. The singular values for $\boldsymbol{k}_{h}$ are plotted in Figure 1(a), which shows $\sigma_{1}(\boldsymbol{k}_{h})$ is two orders of magnitude larger than the rest. While progress can be made with this assumption, the forcing is, in fact, structured and influences the streamwise decay of the resolvent modes. This can be seen in Figure 1(b-c), which compares the streamwise component of the resolvent and DMD modes. The overall agreement is good although the mode shapes begin to diverge around $x=5$.

Similar observations can be made for EQ1. Figure 1(j) indicates very low-rank behavior for $\boldsymbol{k}_{h}$ and, just like the cylinder case, the resolvent accurately predicts the mode shape as seen in Figure 1(k-l) where the $u$-component of the first resolvent mode is compared to the true solution.

When the resolvent is not low-rank, such as at $\boldsymbol{k}=2\boldsymbol{k}_{h}$, as seen in Figure 1(d,m), respectively, the leading mode from a SVD is unlikely to predict the structure of the velocity fluctuations. Consequently, there is no agreement between the resolvent and DMD modes, which are plotted in Figure 1(e-f), respectively, for the cylinder. The resolvent mode and true velocity fluctuations for EQ1 in Figure 1(n-o), respectively, also disagree.

In both cases, there is no linear mechanism active at $2\boldsymbol{k}_{h}$. The analysis of Turton et al., (2015) can be extended to argue which scale interactions are needed in Equation 6 to recover $\hat{\boldsymbol{f}}(2\boldsymbol{k}_{h})$. The linearized NSE are rewritten as

[TABLE]

[TABLE]

For $\|\hat{\boldsymbol{u}}(q\boldsymbol{k}_{h})\|\sim\epsilon^{|q|}$, $\epsilon\ll 1$, as it must be for large singular values (Figures 1 (a,j)), and considering $q=2$,

[TABLE]

Retaining terms on the order of $\epsilon^{2}$ results in

[TABLE]

where Equation 4 is used to approximate $\hat{\boldsymbol{u}}(\boldsymbol{k}_{h})$. Terms with the subscript 1 are associated with the wavenumber vector $\boldsymbol{k}_{h}$.

Approximating the forcing $\hat{\boldsymbol{f}}(2\boldsymbol{k}_{h})$ by the self-interaction of $\hat{\boldsymbol{\psi}}_{1}(\boldsymbol{k}_{h})$ leads to a reasonable comparison with the true forcing for both the cylinder (Figures 1(g-h)) and EQ1 (Figures 1(p-q)). The modes excited by these forcings, $\hat{\boldsymbol{u}}_{f}(2\boldsymbol{k}_{h})=\mathcal{H}(2\boldsymbol{k}_{h})\hat{\boldsymbol{f}}(2\boldsymbol{k}_{h})$ show much improved agreement with the true response at $2\boldsymbol{k}_{h}$ (Figures 1(i,r).)

The EQ1 structures in Figures 1(r) and 1(o) are nearly identical. While there is significant improvement in the mode shape for the cylinder, the agreement between Figure 1(i) and the DMD mode in Figure 1(f) is not perfect. The residual discrepancy for the cylinder is due to the fact \hat{\boldsymbol{\psi}}_{1}({{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}})\neq\hat{\boldsymbol{u}}({{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}) and only one triad was considered. Nevertheless, the forced resolvent mode captures the correct symmetries and is much more representative of the fluctuations at 2{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}.

The mode at $2\boldsymbol{k}_{h}$ can be considered ‘parasitic’ in that it is driven by the dominant linear mechanism at $\boldsymbol{k}_{h}$ and emerges from the resolvent operator being forced by the $\boldsymbol{k}_{h}$ velocity modes. It plays a role in the nonlinear forcing of the $\boldsymbol{k}_{h}$ mode. Setting $q=1$ in Equation 7 results in

[TABLE]

It was argued by Turton et al., (2015) that the $\epsilon^{3}$ terms could be neglected resulting in $\boldsymbol{L}$ containing a marginally stable eigenvalue as discovered by Barkley, (2006). Retaining those terms and substituting equations 4 and 10, we can express the forcing as

[TABLE]

where notably \hat{{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{f}}}}^{\prime}(\boldsymbol{k}_{h}) can be computed from the resolvent with only knowledge of the mean. From Equation 4 we can solve for the amplitude as

[TABLE]

From Equation 11, retaining terms at $\mathcal{O}(\epsilon)$ and treating the higher-order terms as an unknown forcing, the resolvent operator can be formulated and used to compute the spatial structure of the mode \hat{\boldsymbol{\psi}}_{1}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\boldsymbol{k}_{h})}}. Consequently, after cascading down the nonlinear interactions, at $\mathcal{O}(\epsilon^{3})$ an equation can be derived for the corresponding amplitude of the mode, in close analogy to a weakly-nonlinear analysis. For EQ1 at $Re=1000$, this approximation of $|\chi_{1}|$ has a relative error with respect to the true value computed directly from the solution of less than $1\%$. As Equation 12 and the subsequent calculation of $|\chi_{1}|$ relies on \hat{\boldsymbol{\psi}}_{1}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\boldsymbol{k}_{h})}} closely approximating the true spatial structure, the agreement in the case of the cylinder is not as good. However if the DMD modes, which are the true velocity fluctuations at {{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}} and 2{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}, are used, they can almost exactly reproduce \hat{\boldsymbol{f}}({{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}) (Symon,, 2018), confirming the dominance of this interaction in the full flow. Moreover, if \mathcal{H}({{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}) is driven by this forcing, it results in a forced resolvent mode which very accurately approximates the DMD mode in Figure 1(c). Thus it is possible to improve predictions from resolvent analysis, even when the operator is low-rank, by approximating the nonlinear forcing.

The results for the circular cylinder and EQ1 can be summarized in Figure 2. A SVD of the resolvent operator at $\boldsymbol{k}_{h}$ yields a highly amplified mode which interacts nonlinearly with its complex conjugate to produce the Reynolds stresses that sustain the mean profile. It also interacts with itself to produce the forcing at $2\boldsymbol{k}_{h}$. When $\mathcal{H}(2\boldsymbol{k}_{h})$ is driven by this forcing, it results in a ‘parasitic’ mode which nonlinearly interacts with the $-\boldsymbol{k}_{h}$ mode to produce the structured nonlinear forcing for $\boldsymbol{k}_{h}$. Thus in the case of the cylinder, where there is a greater mismatch between the resolvent mode and the DMD mode, the original \boldsymbol{\psi}_{1}({{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{k}_{h}}}) can be iteratively improved using the feedback loop in Figure 2.

4 Summary

Resolvent analysis was used to identify the linear mechanisms and dominant nonlinear interactions required to provide a self-consistent description of low-Reynolds number cylinder flow and the Couette equilibrium solution EQ1. It was demonstrated that the leading resolvent mode, obtained by approximation of the resolvent in a region where it is low-rank, is able to correctly capture the spatial structure of the most energetic mode. Considering the relevant scale interactions, i.e. approximating the forcing and leveraging the low rank nature of the energetic mode generates improved approximations of the second harmonic, here termed a parasitic mode because of its link to the interactions of the energetic mode and the lack of low-rank behavior of its resolvent. These ideas can be consolidated to provide a means to determine the amplitude of the dominant energetic mode, and hence effectively close the system, using knowledge of only the mean profile, in cases where the resolvent is not low-rank for all scales. It is hypothesized that this framework may also have relevance for flows with multiple linearly dominant energetic modes. The hope would be to use the resolvent operator to obtain low-dimensional representations of the important linearly-amplified structures, which are hypothesized to correspond to identifiable physical mechanisms, and use the notion of parasitic modes to provide an accurate representation of the weakly energetic but nonlinearly relevant scales.

The support of ONR under grants N00014-17-1-2307 and N00014-17-1-3022 is gratefully acknowledged.