# A second-order resolvent formulation for the analysis of turbulent flow structures

Quentin Chevalier, Lutz Lesshafft, André V G Cavalieri

► To cite this version:

Quentin Chevalier, Lutz Lesshafft, André V G Cavalieri. A second-order resolvent formulation for the analysis of turbulent flow structures. Comptes Rendus. Mécanique, 2023, 351, pp.355 - 371. 10.5802/crmeca.193. hal-04224163

HAL Id: hal-04224163

https://hal.science/hal-04224163v1

Submitted on 2 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Comptes Rendus
## Mécanique

Quentin Chevalier, Lesshafft Lutz and André V. G. Cavalieri

A second-order resolvent formulation for the analysis of turbulent flow structures

Volume 351 (2023), p. 355-371

Published online: 26 September 2023

https://doi.org/10.5802/crmeca.193

<figure><img src="image_2.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org
e-ISSN : 1873-7234
Spontaneous Articles / *Articles spontanés*

# A second-order resolvent formulation for the analysis of turbulent flow structures

Quentin Chevalier*, a, Lesshafft Lutz a and André V. G. Cavalieri b

<sup>a</sup> Laboratoire d'Hydrodynamique, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

<sup>b</sup> Instituto Tecnológico de Aeronáutica, 12228-900 São José dos Campos, Brazil

*E-mails*: quentin.chevalier@polytechnique.edu (Q. Chevalier), lutz.lesshafft@polytechnique.edu (L. Lesshafft), andre@ita.br (A.V.G. Cavalieri)

**Abstract.** An attempt to improve the accuracy of resolvent-based predictions by including velocity correlations in the linear model is developed here. Closure assumptions for unresolved nonlinearities are thus pushed back to a higher order. Turbulent channel flow is considered as a test case: response and forcing modes obtained from singular value decomposition of the new resolvent model are compared to Spectral Proper Orthogonal Decomposition (SPOD) modes extracted from a Direct Numerical Simulation (DNS) database. The performance of the approach is also measured against previous resolvent-based models. The new model does not yield significant global improvement, but does improve predictions in some regions. Further work on the method should target the linear modeling of the velocity-pressure gradient correlation tensor.

**Keywords.** Resolvent analysis, Turbulence, Reduced-order modeling, Theoretical fluid mechanics, Linear analysis.

**Funding.** The first author is supported by the *Direction Générale de l'Armement*. This work has also received funding from the CleanSky2 Joint Undertaking (JU) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 785303.

Manuscript received 27 June 2022, revised 14 February 2023 and 5 April 2023, accepted 25 April 2023.

## 1. Introduction

Incompressible turbulent flows are very diverse and cover a broad range of applications. However, they are also notoriously difficult to model. In this work the choice was made to solve the Navier-Stokes equations (NSE) in frequency space instead of approaching time-dependent solutions directly. Here, a novel spectral reduced-order approach to model incompressible turbulent flows will be presented. This method is derived from resolvent modeling with an input from Reynolds-stress models. The resolvent approach proposes a linear model of a fully turbulent flow at a given frequency. Interpreting nonlinearities in the NSE as a forcing input into the linearized

<footnote>* Corresponding author.</footnote>
fluctuation equations, the most amplified response structures are taken as a prediction for energetic turbulent flow structures, for instance in [1, 2]. This approach is most often applied in frequency space and comes with all the advantages and drawbacks of a linear model, namely inability to account for all terms, simplicity of implementation, as well as computational efficiency. Reynolds-stress models originate from a theoretical attempt to improve the quality of Reynolds averaged Navier-Stokes models by adding second-order terms in the relevant equations as in [3, chapter 11].

Resolvent analysis has already demonstrated its ability to accurately predict large-scale structures in turbulent channel flow at a given frequency; [4] already exposed excellent agreement between Direct Numerical Simulation (DNS) data and a combination of twelve resolvent modes with weights derived from an optimisation problem computed at every wavenumber in a spectral framework. Said agreement persisted throughout the DNS energy spectrum and across the channel without any empirical modeling. Quantifying complexity in the same study clearly shows the main advantage of reduced-order modeling - indeed, the computations involved with the resolvent model are more than three orders of magnitude cheaper than the relevant DNS, and they nonetheless produce high quality estimations of the fluctuations. Despite these encouraging results, this approach provides no information on nonlinear terms in the original equations, which can also represent an engineering objective.

Indeed, other authors have given thought to this idea of a low-order modelling of turbulent flows that would be more cost-efficient than DNS - good enough for engineering purposes, and orders of magnitude less expensive. [5] consider a stratified channel flow and prove that on the limited range of scales considered, a resolvent model is able to reproduce an energy budget similar to that of a DNS.

As part of the effort to improve resolvent model performance for wall-bounded flows, [6] apply a decomposition of the velocity field in solenoidal and non-solenoidal parts in the Orr-Sommerfield-Squire (OSS) formulation. Exact Coherent States (ECS) for simple Couette and Poiseuille flows were used to validate the model's performance. An astute method to explicitly compute interaction coefficients as part of an optimisation problem led to a marked increase in performance. A single resolvent mode was demonstrated to accurately represent some ECS solutions. Finally, a characterization of nonlinear terms was performed using the DNS of a channel flow, but stopped at the conclusion that its spectral signature is non trivial. The most important drawback of this approach is that the OSS formulation makes the corresponding nonlinear terms much harder to interpret and model, thus limiting the possibility for further model improvement.

[7] computed Spectral Proper Orthogonal Decomposition (SPOD) and resolvent modes for a wide range of wavenumbers from a DNS database. That study avoids the use of eddy viscosity and presents large discrepancies between SPOD and resolvent modes. Possible explanations include different physical phenomena which dominate various flow regions, such as the lift-up effect which seems to produce favourable conditions for model and data agreement. [7, Figures 3 and 4] display the correlation between SPOD and resolvent modes which clearly illustrate the appeal of a new reduced-order method that would be more general than the current state of the art.

[8] present good agreement between a DNS power spectral density (PSD) and a modeled one, using eddy viscosity and a resolvent formalism in a channel flow, as well as an engineered forcing that corresponds better to true flow statistics than white-noise.

More recently, [9] detailed how to obtain a "coloured" forcing term more representative of flow physics. This approach results in a tool with predictive value, though its generality is somewhat lessened by the use of an eddy viscosity model.
[10] present a resolvent-based estimation of turbulent channel flow using wall fluctuations as an input. This work also highlights the importance of “coloured” forcing in the NSE, as white-noise forcing leads to large fluctuations in the logarithmic layer.

Another recent study by [11] stresses how the addition of eddy viscosity makes a significant difference in the resolvent model's ability to properly capture flow energy. Mainly concerned with the latter quantity, [11] outline the balance between production, dissipation, and nonlinear transfer of turbulent fluctuation kinetic energy.

Indeed, the ability of eddy viscosity to improve resolvent models and enable them to reproduce fine structures close to the wall present in DNS is desirable in channel flows, as stressed by [12]. Although the addition of eddy viscosity makes the resolvent model more robust across wave speeds, the conclusion of [12] was that a scale dependant eddy viscosity should be considered.

A major drawback of resolvent methods listed above is their inability to predict nonlinear terms representative of the real flow. Indeed, the analysis of such a linear operator provides response modes of physical quantities and corresponding optimal forcing terms that represent flow nonlinearities. These forcing terms are rarely representative of actual physical flow quantities. In an effort to address this limitation, and in order to obtain the “coloured forcing” directly in the new approach, a usual resolvent model for the velocities was coupled with equations for velocity correlations such as to obtain a “second-order resolvent model” that would accurately account for both dominant order response and nonlinear terms in the NSE, called forcing terms.

This idea can to some degree be found in the stochastic structure stability theory of [13]. Indeed, this work details how to obtain a linear model for a turbulent mean flow and the associated Reynolds stresses where both evolve in time. The authors acknowledge that for large Reynolds numbers the method becomes highly intractable, preferring instead a reduced order model which uses a few realisations of the fluctuations forced with random white noise to approximate non-linear terms in the baseflow equations as the system becomes large.

Another similar approach goes by the name of generalised quasilinear approximation. This formalism explicitly introduces a range of frequencies that are allowed to interact with one another. [14] explain how this is a generalisation of the quasilinear approach where a single frequency is considered independently of all others. Indeed, in order to close the model, some interactions between waves are cut off from the generalised model. This approach does not always fare well in describing the whole system, but almost always constitutes an improvement compared to the classical quasilinear approximation.

This paper is organised as follows: Section 2 presents the methodology of usual resolvent analysis as well as the additional equations specific to the new approach, Section 3 gives more details about parameter choices and the test dataset, Section 4 presents the two linear systems and links them together, Section 5 exhibits results across the wavenumber space and details a specific set of parameters. Finally, Section 6 summarises the study.

## 2. Methodology

This section aims at laying out necessary formalism and equations. Starting from the NSE, and a given baseflow, analytical tools to obtain insights into the coherent structures of the flow will be built. The first goal is to obtain a linear system $u' = \underline{R}f$ which transforms a given forcing $f$ representing non-linear terms of the relevant equations into a fluctuation velocity field $u'$ through a resolvent operator $\underline{R}$. The standard approach of defining the linear resolvent operator, as widely applied in the recent literature, is presented in Section 2.1. The proposition for an extended linear resolvent operator is introduced in Section 2.2. Subsequently, Section 2.3 will
explain how a linear system can be exploited to obtain its most amplified fluctuations and optimal forcing.

## 2.1. First-order equations

Writing the incompressible NSE for total velocity $\underline{u} = \underline{U} + \underline{u}'$ and pressure $p = P + p'$ with ensemble averages $\underline{\underline{u}} = \underline{U}$ and $\underline{\underline{p}} = P$ and the convention $[\nabla U]_{ij} = \partial U_i / \partial x_j$ gives

$$D\underline{u}' + \underline{\underline{\nabla}}\underline{U} \underline{u}' + \underline{\underline{\nabla}} \cdot (\underline{\underline{\tau}} - \underline{\underline{\tau}}) + \underline{\underline{\nabla}}p' - \frac{1}{Re}\underline{\underline{\Delta}}\underline{u}' = 0, \quad (1)$$

with $\underline{\underline{\tau}} = \underline{u}'\underline{u}'^T$, and the total derivative is defined as $D\underline{u}' = \partial_t \underline{u}' + \underline{\underline{\nabla}}\underline{u}'\underline{U}$. Since the term $\underline{\underline{\nabla}} \cdot (\underline{\underline{\tau}} - \underline{\underline{\tau}})$ is nonlinear in the fluctuation quantities $q_{\underline{u}} = [\underline{u}'^T p']^T$, it may be interpreted as forcing as in [4]. The linear model system for flow fluctuations is written

$$\left\{ \begin{array}{l} \underline{\underline{\nabla}} \cdot \underline{u}' = 0 \\ D\underline{u}' + \underline{\underline{\nabla}}\underline{U} \underline{u}' + \underline{\underline{\nabla}}p' - \frac{1}{Re}\underline{\underline{\Delta}}\underline{u}' = -\underline{\underline{\nabla}} \cdot (\underline{\underline{\tau}} - \underline{\underline{\tau}}) \quad \Leftrightarrow \quad \underline{\underline{L}}u\underline{q}_{\underline{u}} = \begin{bmatrix} 0 \\ \underline{\underline{f}}_{\underline{u}} \end{bmatrix}. \end{array} \right. \quad (2)$$

Note that the first line of the linear system relates to the incompressibility condition, which is not forced. Resolvent analysis thus circumvents the closure problem by labelling unclosed terms as forcing. Conceptually, one goes from trying to solve for $q_{\underline{u}} = [\underline{u}'^T p]^T$ to looking for coherent structures favoured by the linear operator $\underline{\underline{L}}u$.

## 2.2. Second-order equations

The key proposition of this paper is to model the $\underline{\underline{\tau}}$ tensor, in an attempt to increase the accuracy of the resolvent model. From equation (1) one obtains

$$\left( D - \frac{1}{Re} \underline{\underline{\Delta}} \right) \underline{\underline{\tau}} + \underline{\underline{\tau}} \underline{\underline{\nabla}} \underline{\underline{U}}^T + \underline{\underline{\nabla}} \underline{\underline{U}} \underline{\underline{\tau}} \\ = - \underline{u}' \underline{\underline{\nabla}} p'^T - \underline{\underline{\nabla}} p' \underline{u}'^T - \frac{2}{Re} (\underline{\underline{\nabla}} \underline{u}')^T \underline{\underline{\nabla}} \underline{u}' - \underline{u}' [\underline{\underline{\nabla}} \cdot (\underline{\underline{\tau}} - \underline{\underline{\tau}})]^T - \underline{\underline{\nabla}} \cdot (\underline{\underline{\tau}} - \underline{\underline{\tau}}) \underline{u}'^T. \quad (3)$$

The right hand side of equation (3) will be considered as forcing terms, denoted $\underline{\underline{f}}$, their common trait being their nonlinearity with respect to $\underline{u}', p'$ and $\underline{\underline{\tau}}$.

Replacing the symmetric $\underline{\underline{f}}$ and $\underline{\underline{\tau}}$ tensors with vectors $\underline{\underline{f}}^v$ and $\underline{\underline{\tau}}^v$, which each contain six independent tensor components, equation (3) can be cast as a linear system

$$\underline{\underline{L}}\underline{\underline{\tau}}\underline{\underline{\tau}}^v = \underline{\underline{f}}^v. \quad (4)$$

## 2.3. Resolvent analysis

Let us define a generic extensor $\underline{\underline{B}}$, which extends a vector $\underline{\underline{f}}$ with zeros

$$\underline{\underline{B}}\underline{\underline{f}} = [\underline{\underline{0}}^T \underline{\underline{f}}^T]^T. \quad (5)$$

Given a state vector $\underline{\underline{q}}$ that includes the fluctuation velocities $\underline{\underline{u}}'$ (for example $q_{\underline{u}} = [\underline{u}'^T p']^T$), let a generic extractor (or observation matrix) $\underline{\underline{H}}$, be the operator that extracts the velocities from the entire state

$$\underline{\underline{H}}\underline{\underline{q}} = \underline{\underline{u}}'. \quad (6)$$

Assume a linear operator $\underline{\underline{L}}$ that links the state vector $\underline{\underline{q}}$ (which again includes the fluctuation velocities $\underline{\underline{u}}'$) and the forcing vector $\underline{\underline{f}}$ through an extensor $\underline{\underline{B}}$ (as is the case in equation (2)) as follows

$$\underline{\underline{L}}\underline{\underline{q}} = \underline{\underline{B}}\underline{\underline{f}}. \quad (7)$$
Using these general matrices $\underline{B}$, $\underline{H}$ and $\underline{L}$, which will be detailed later for specific systems, let the resolvent matrix be $\underline{R} = \underline{H} \underline{L}^{-1} \underline{B}$, so that one finally has $\underline{u}' = \underline{R} \underline{f}$.

The Singular Value Decomposition (SVD) of the $\underline{R}$ matrix, $\underline{R} = \sum_i \sigma^{(i)} \underline{\psi}^{(i)} \underline{\phi}^{(i)H}$ yields $\sigma^{(i)}$ the gain for mode $i$ of the $\underline{R}$ matrix, ordered in decreasing value.

The normalised vectors $\underline{\phi}^{(i)}$ (right singular vectors) thus computed are called in the following the forcing modes, and $\underline{\psi}^{(i)}$ (left singular vectors) will be referred to as the response modes of the resolvent system. The positive real-valued $\sigma^{(i)}$ represents a gain between the norm of the forcing and that of the response. The first of these modes, $\underline{\phi}^{(1)}$ and $\underline{\psi}^{(1)}$, are optimal in the sense that they represent the forcing/response pair associated with the highest gain.

This entire process is not new in the case of system (2) and was performed by [7, 9], amongst others.

### 3. Test case: channel flow

#### 3.1. Flow configuration

The specific case of a fully developed channel flow between two infinite plates was studied, separated by a distance $2h$. The presence of viscosity imposes $\underline{u}' = \underline{0}$ at both walls. This naturally translates into $\underline{\tau} = \underline{0}$ at the wall. Even though there is physical motivation to impose additional boundary conditions for $\underline{\tau}$, the order of the system prevents the addition of further constraints.

For a wide range of parameters, forcing and response structures are concentrated near the channel walls, with very weak coupling between the top and the bottom halves of the channel. This weak coupling leads to poorly conditioned numerical systems when the full channel height is resolved. To address this problem, following [7], only the half-channel height was resolved, and symmetry conditions at the centre line imposed as follows.

With $x$ streamwise, $y$ wall-normal and $z$ spanwise directions, one obtains $\underline{U}(x, y, z) = U(y)\underline{e}_x$ with $U$ a function made strictly even across the channel. If $u'_z$ is an even (respectively odd) function, and thus symmetric (resp. anti-symmetric) across the channel centre line, then it follows from equation (9) that $p'$ is an even (resp. odd) function. In turn, this gives odd (resp. even) $u'_y$, and even (resp. odd) $u'_x$.

The parity of the six components of $\underline{\tau}$ is constant - $u'_x u'_y$ and $u'_y u'_z$ are always odd, $u'_x u'_z$ is always even, regardless of the parity of $u'_z$. For the sake of conciseness, the analysis was constrained to modes with even $u'_z$, yet in the course of the study this choice had little influence on the results. Symmetry was enforced a posteriori, after DNS calculations and singular value decomposition, as in previous work by [7]. In the following, only a half-channel is studied.

#### 3.2. Fourier transforms

Consider Fourier transforms in space and time for all flow quantities

$$q(x, y, z, t) = \hat{q}(y) \exp (ik_x x + ik_z z - i\omega t) \quad (8)$$

Note that all terms depending solely on averaged flow quantities disappear through Fourier transform at non-zero wavenumber. This is crucial because the forcing of equation (2) becomes $\underline{f}_u = -\underline{\nabla} \cdot \underline{\tau}$ and the link with system (4) is now apparent.

Writing $\delta_i^x$ the Kronecker delta, $dU = dU/dy$ baseflow shear and $F = -i\omega + ik_x U - (\partial_y^2 - k_x^2 - k_z^2)/Re$ gives the channel flow case new directional formulations of equations (1) and (3),

$$F \hat{u}'_i + dU \delta_i^x \hat{u}'_y + ik_x \hat{\tau}_{ix} + \partial_y \hat{\tau}_{iy} + ik_z \hat{\tau}_{iz} + \hat{\partial}_i p' = 0, \quad (9)$$

$$F \hat{\tau}_{ij} + dU (\delta_i^x \hat{\tau}_{ij} + \delta_j^y \hat{\tau}_{iy}) = \hat{f}_{ij}. \quad (10)$$
### 3.3. Wall-normal discretisation

Discretisation of the full channel in the wall-normal direction was performed using $n = 129$ Chebyshev–Gauss–Lobatto nodes, with $y = 0$ being the centre of the channel,

$$ \forall k \in \mathbb{N}, k \le n, y_k = h \cos \frac{k\pi}{n}, \quad (11) $$

to allow for more points closer to the channel wall. The operator $\partial_y$ is implemented to spectral order using Chebyshev polynomials following [15]. A test case for convergence relative to $n$ is detailed in Section 5.2.

In the following, the notation $d$ will be used to denote operators discretised using this spacing.

Consider an integration matrix $\underline{W}_d$ so that $\underline{\hat{u}}_d^H \underline{W}_d \underline{\hat{u}}_d \approx \int_{-h}^h \|\underline{\hat{u}}\|^2(y)dy$ in the chosen discretisation. Modes computed using the process detailed in Section 2.3 satisfy

$$ \underline{\phi}_d^{(i)H} \underline{\phi}_d^{(j)} = \underline{\psi}_d^{(i)H} \underline{\psi}_d^{(j)} = \delta_j^i, \quad (12) $$

so they are orthonormal for an equal spacing discretisation only. This can be corrected as follows.

Using the Cholesky decomposition $\underline{W}_d = \underline{C}_d^T \underline{C}_d$, one may transform the resolvent as

$$ \underline{\tilde{R}}_d = \underline{C}_d \underline{R}_d \underline{C}_d^{-1}. \quad (13) $$

Then one can then take the Singular Value Decomposition (SVD) of the $\underline{\tilde{R}}_d$ matrix, $\underline{\tilde{R}}_d = \sigma^{(i)} \underline{\tilde{\psi}}_d^{(i)} \underline{\tilde{\phi}}_d^{(i)H}$, and obtain the true modes of the system $\underline{\phi}_d^{(i)} = \underline{C}_d^{-1} \underline{\tilde{\phi}}_d^{(i)}$ and $\underline{\psi}_d^{(i)} = \underline{C}_d^{-1} \underline{\tilde{\psi}}_d^{(i)}$. These new modes are indeed orthonormal with respect to the chosen spacing,

$$ \underline{\phi}_d^{(i)H} \underline{W}_d \underline{\phi}_d^{(j)} = \underline{\psi}_d^{(i)H} \underline{W}_d \underline{\psi}_d^{(j)} = \delta_j^i. $$

### 3.4. DNS dataset

The analysis is based on spectral DNS data of a channel flow. The data was used to extract the time-averaged mean flow and to compute the leading SPOD modes following [16], against which the linear model results will be compared. This dataset has been validated in previous studies such as [7, 9].

Two Reynolds numbers are defined as

$$ Re = \frac{U_{bulk}h}{\nu} \quad \text{and} \quad Re_{\tau} = u_{\tau} Re, \quad (14) $$

where $U_{bulk}$ is the mean velocity, $h$ is the channel half-width, $\nu$ is the kinematic viscosity, and $u_{\tau} = \sqrt{\max(dU)/Re}$ is the friction velocity. Inner scalings are further defined as $x^+ = x \cdot u_{\tau}/\nu$ and $u^+ = u/u_{\tau}$. Two DNS datasets are used, with ($Re = 2800$, $Re_{\tau} = 180$) and ($Re = 10000$, $Re_{\tau} = 550$).

The time-wise Fourier transform necessary for SPOD computation, was performed using the Welch method with a second-order sinusoidal window of size $N = 256$, with 75% overlap. The Welch transform process also featured a force correction to counter the spurious influence of the windowing function as in [17]. More information on the chosen DNS parameters can be found in Table 1, which presents the size of the box in streamwise and spanwise directions $L_x$ and $L_z$, the grid steps in these directions in inner scaling $\Delta x^+$ and $\Delta z^+$, as well as the smallest and largest steps in the normal direction $\Delta y_{\min}^+$ and $\Delta y_{\max}^+$.

<table><caption>**Table 1.** Characteristic parameters of the DNS database.</caption><thead><tr><th>Re<sub>τ</sub></th><th>Re</th><th>L<sub>x</sub></th><th>L<sub>z</sub></th><th>Δx<sup>+</sup></th><th>Δz<sup>+</sup></th><th>Δy<sub>min</sub><sup>+</sup></th><th>Δy<sub>max</sub><sup>+</sup></th></tr></thead><tbody><tr><td>179</td><td>2800</td><td>4π</td><td>2π</td><td>11.78</td><td>5.89</td><td>5.42 · 10<sup>-2</sup></td><td>4.42</td></tr></tbody></table>
## 4. Final linear systems

In the previous two sections a variety of transformations and hypotheses that introduced much in the way of notation was detailed. For the rest of this study, the hat notation for all Fourier coefficients (see Subsection 3.2) will be removed, as well as the prime symbols on fluctuations (see Subsection 2.1), the *d* (see Subsection 3.3) and the *v* (see Subsection 2.2) superscript notations in the interest of clarity.

In this section, the final form of the linear systems of interest will be detailed in Sections 4.1 and 4.2. Section 4.3 explains the procedure for the comparison between the linear model results and the leading SPOD modes of the DNS data.

### 4.1. First-order system

Equation 9 can be directly discretised to yield a $4n \times 4n$ discrete system, where $[dU]_{ij} = \delta_i^j dU_i$ and $[F]_{ij} = \delta_i^j F_i$

$$ \left[ \begin{matrix} ik_x \underline{I}_n \underline{\partial}_y & ik_z \underline{I}_n \underline{0}_n \\ \underline{F} & \underline{dU} \\ \underline{0}_n & \underline{0}_n \\ \underline{0}_n & \underline{0}_n \end{matrix} \right] \left[ \begin{matrix} \underline{u}_x \\ \underline{u}_y \\ \underline{u}_z \\ \underline{p} \end{matrix} \right] = \underline{L}_u \underline{q}_u = \left[ \begin{matrix} \underline{0}_n \\ \underline{f}_u \end{matrix} \right]. \quad (15) $$

The extensor $\underline{B}_u$ and extractor $\underline{H}_u$ operators therefore are non-square identity matrices,

$$ \underline{B}_u = \left[ \begin{matrix} \underline{0}_n \\ \underline{I}_{3n} \end{matrix} \right] \qquad \underline{H}_u = \left[ \begin{matrix} \underline{I}_{3n} & \underline{0}_n \end{matrix} \right]. \quad (16) $$

So far, this formulation is not novel. It has been done in [2, 7, 9, 10]. The process detailed in Sections 2.3 and 3.3 yields a resolvent matrix $\underline{R}_{1\text{rst}}$ order for the so-called “first-order system”, as well as associated response modes $\underline{\psi}_{1\text{rst}}^{(i)}$ and forcing modes $\underline{\phi}_{1\text{rst}}^{(i)}$.

The currently most successful resolvent-based modelling of energetic coherent structures in channel flow uses a first-order formulation (Section 2.1) with an eddy viscosity model. For comparison, the Cess model employed in reference studies [9, 12] is reproduced here

$$ \frac{v_t}{v} = \frac{1}{2} \sqrt{1 + \frac{\kappa^2 Re_{\tau}^2}{9} (1 - y^2)^2 (1 + 2y^2)^2 \left[ 1 - \exp\left(-Re_{\tau} \frac{1 - |y|}{A}\right) \right]^2} - \frac{1}{2}, \quad (17) $$

using the friction Reynolds number as defined in equation (14). When applied, this model replaces the molecular diffusion term in equation (9), $T_{\text{diff},v} = \Delta \underline{u}/Re$, by a more complete term $T_{\text{diff},v_t} = \underline{\nabla} \cdot [(1 + v_t/v)(\underline{\nabla} \underline{u} + \underline{\nabla} \underline{u}^T)]/Re$. In order to verify that the above model accurately represents the flow dynamics, a comparison between the mean flow computed using the eddy viscosity and the one obtained by DNS was performed. The computation of the mean flow from the eddy viscosity was done as in [18, section 2]. The velocity profiles are compared in Figure 1. Differences between the two curves exist but remain small.

This is the resolvent system used by [7, 9], with an eddy viscosity representing at least part of the unknown nonlinear terms. Applying once more the process in Sections 2.3 and 3.3 yields a resolvent matrix $\underline{R}_{\text{eddy viscosity}}$, as well as associated response modes $\underline{\psi}^{(i)}$ and forcing modes $\underline{\phi}^{(i)}$.

### 4.2. Second-order system

The original proposition in the method is to include equation (10) inside a unified resolvent operator that includes velocity, pressure, and $\underline{\tau}$ in the state vector. Recall that the motivation is to
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Mean velocity profiles from DNS and computed as a steady state with the Cess eddy viscosity model (equation (17)).</figcaption></figure>

improve the capability of the resolvent operator to correctly predict flow physics both in terms of perturbations and their first moment.

Equation 10 translates to a $6n \times 6n$ linear system with the vector notations $\underline{\tau}$ and $\underline{f}$ (see Section 2.2) written as

$$ \underline{L}_{\tau} \underline{\tau} = \underline{f}_{\tau}. \qquad (18) $$

As before, $\underline{f}_{\tau}$ represents the forcing of the second order equations. Together, equations (15) and 18 form a general system formulation that is $10n \times 10n$ in size,

$$ \begin{bmatrix} \underline{L}_u & \underline{A} \\ \underline{0}_{4n} & \underline{L}_{\tau} \end{bmatrix} \begin{bmatrix} q_u \\ \underline{\tau} \end{bmatrix} = \begin{bmatrix} \underline{0}_{4n} \\ \underline{f}_{\tau} \end{bmatrix} \Leftrightarrow \underline{L} \underline{q} = \underline{f}. \qquad (19) $$

The equation above makes use of the coupling matrix $\underline{A}$ of size $4n \times 6n$ defined as

$$ \underline{A} \underline{q} = \begin{bmatrix} \underline{0}_n \\ \underline{\nabla} \cdot \underline{\tau} \end{bmatrix}. \qquad (20) $$

The extensor $\underline{B}$ and extractor $\underline{H}$ may be explicitly derived

$$ \underline{B} = \begin{bmatrix} \underline{0}_{4n} \\ \underline{I}_{6n} \end{bmatrix}, \qquad \underline{H}_u = \begin{bmatrix} \underline{I}_{3n} & \underline{0}_{7n} \end{bmatrix}. \qquad (21) $$

The process detailed in Sections 2.3 and 3.3 yields a new resolvent matrix $\underline{R}_{2nd}$ order, with its associated response modes $\underline{\psi}^{(i)}$ and forcing modes $\underline{\phi}^{(i)}$. Note that the latter correspond to the forcing terms of equation (10) instead of those of equation 9.

### 4.3. Comparison of model versus data

The process detailed above was implemented in a Matlab code (version R2019b) using the svds (subset singular value decomposition) function library to obtain resolvent modes.

The most amplified response and forcing modes obtained by the resolvent model was compared with the most energetic SPOD modes obtained from DNS data. Said SPOD modes are computed using the method of snapshots detailed in [16]. This method allows us to obtain SPOD velocity modes, denoted $\underline{\xi}^{(i)}$. The same approach was used to compute SPOD nonlinear modes, obtained directly from the nonlinear term of equation (1) in the DNS.
SPOD modes represent the orthonormal basis that captures the most of the flow energy for any given rank as established in [2]. They are obtained as eigenvectors of a Cross-Spectral Density (CSD) matrix computed at a given frequency. SPOD may be performed for both velocities and nonlinear terms of the NSE. While the gain-based ranking of the modes gives an indication of each mode's prevalence in the flow, there is no causal association between SPOD nonlinear and forcing modes obtained in this manner. In other words, $\underline{\zeta}^{(i)}$ is not related to $\underline{\zeta}^{(i)}$, and both modes are computed independently.

The matrix system detailed in equation (19) is exact but the flow response $q$ can only be determined if the forcing $f_{\underline{\mu}}$ is known. Because of the large gain $\sigma^{(1)}$ and large gain separation $\sigma^{(1)}/\sigma^{(2)}$ typically observed in a resolvent operator, it is expected that the first left-singular vector $\underline{\psi}^{(1)}$ approximates $\underline{\zeta}^{(1)}$, the actual flow fluctuations obtained from SPOD. This has been observed in previous studies such as [7, 9], and the coherent structures in a flow are of engineering and scientific interest by themselves, especially if they are provided by a numerically cheap linear analysis process.

However, as pointed out for instance by [19, section III.c] and by [20], there is no reason why the nonlinear terms in the actual flow $\underline{\zeta}^{(i)}$ should be approximated by the right-singular vectors $\underline{\phi}^{(i)}_{\underline{1}^{\text{rst}}\text{order}}$ of the resolvent matrix $\underline{R}_{\underline{1}^{\text{rst}}\text{order}}$. Indeed $\underline{\zeta}^{(1)}$ represent the best low-rank approximation of flow nonlinear terms, whereas $\underline{\phi}^{(1)}_{\underline{1}^{\text{rst}}\text{order}}$ only represents the optimal input to obtain $\underline{\psi}^{(1)}_{\underline{1}^{\text{rst}}\text{order}}$ with respect to $\underline{R}_{\underline{1}^{\text{rst}}\text{order}}$.

Thus, there is no real physical insight to explain why $\underline{\zeta}^{(1)}$ and $\underline{\phi}^{(1)}_{\underline{1}^{\text{rst}}\text{order}}$ would have a similar structure, yet it is accurate to say that only the part that is strongly amplified by the linear operator is expected to have an important influence on flow behaviour. By contrast, $\underline{\nabla} \cdot \underline{R}_{\underline{\tau}} \underline{\phi}^{(1)}_{\underline{2}^{\text{nd}}\text{order}}$ includes additional physics that might yield a better approximation quality of $\underline{\zeta}^{(1)}$.

# 5. Results

## 5.1. Parameter space exploration

For a first comparison with results from [7, 9], the agreement between SPOD data and resolvent predictions was evaluated by varying $\underline{\lambda}^{+} = [\underline{\lambda}_{x}^{+} \ \underline{\lambda}_{z}^{+}]^{T}$. The + superscript denotes inner scaling as defined in Section 3.4. This agreement is quantified in terms of the projection coefficient $\underline{\beta} = |\underline{\psi}^{(1)} \underline{W} \underline{\zeta}^{(1)H}|$, i.e. between first resolvent mode and first SPOD velocity mode. As stated in Section 4.3, the dominant velocity mode obtained through SPOD from DNS data is compared with a left-singular vector of the resolvent operator.

Only the $Re_{\tau} = 180$ ($Re = 2800$) case will be discussed here, but the $Re_{\tau} = 550$ ($Re = 10000$) case was also explored with similar conclusions. Isocontours of $\underline{\beta}$ for the range of $\underline{\lambda}^{+}$ available at the fixed frequency $f^{+} = 10^{-2}$ are shown in Figure 2a for the new model and Figure 2b for the eddy viscosity enhanced first model.

This graph indicates no clear incentive to use the “second-order” resolvent model over the first-order eddy viscosity resolvent model throughout the $\underline{\lambda}^{+}$ plane at this frequency. Even though the new model performs slightly better over a restricted area in the low $\underline{\lambda}_{x}^{+}$, high $\underline{\lambda}_{z}^{+}$ regime, its performance drops elsewhere. The value $f^{+} = 10^{-2}$ was chosen as representative of the near-wall cycle, for consistency with [7].

Results obtained with different frequencies ($10^{-4}$, $10^{-3}$, $10^{-1}$, 1) show a similar pattern. The isocontours are also shown at $f^{+} = 10^{-4}$ in Figures 2c and 2d, where one can see that the eddy viscosity outperforms the second-degree resolvent model throughout the domain overall.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Projection coefficient of first SPOD against resolvent response mode for a range of wavenumbers $\beta = |\psi^{(1)} \underline{W} \xi^{(1)H}|$ at $Re = 2800$ for two different frequencies.</figcaption></figure>

## 5.2. Near wall structures

In the following, a case study for ($\lambda_x^+ = 1130$, $\lambda_z^+ = 113$, $c^+ = \omega^+ \lambda_x^+ / (2\pi) = 12$) is presented. This choice of parameters corresponds to the “near-wall” case discussed in [9, paragraph 3.1]. Indeed, these wavenumbers represent a particularly energetic point in frequency space for $y^+ \approx 15$, so relatively close to the wall.

To verify that results are independent of the number of points $n$, an additional set of calculations at $n = 257$ was performed for the near wall regime. This calculation with roughly double the $y$ resolution as the one presented below differed by less than 2% from the following results in all aspects. Therefore, it is considered proven in the rest of this study that the results are converged with respect to $n$.

Figure 3 shows the absolute value of the first SPOD velocity and resolvent response modes for a variety of models along the three flow directions. It features the “first-order model” from Section 4.1 with and without eddy viscosity, as well as the “second-order model” detailed in Section 4.2. Likewise, Figure 4 presents the first SPOD nonlinear and resolvent forcing modes for all directions. As detailed in Sections 2 and 4.3, the modes compared here are the first modes
with respect to gain obtained from eigenvalue decomposition of a CSD matrix, and from the SVD of a resolvent matrix for a specific model. There is a noteworthy exception to this, “second-order forcing terms” computed as normalised $\bar{\nabla} \cdot \bar{R}_t \bar{\phi}_{2\text{nd order}}^{(1)}$. Thus, all modes presented below verify $\underline{x}^H \underline{W} \underline{x} = 1$. As discussed in Subsection 3.1, only half a channel was considered supposing the velocity along $z$ to be even for all modes.

All in all, Figure 3 yields overall positive results. In the streamwise direction (Figure 3a), where the fluctuations are most intense, the “second-order” resolvent model succeeds in capturing most of the SPOD mode. The model captures very well close wall behaviour in this direction with no need for empirical eddy viscosity or wall functions. However, the new resolvent model does tend to underestimate the fluctuation amplitudes in the other two directions (Figures 3b and 3c), performing rather poorly in the spanwise direction, where it smoothens out the flow oscillations to an extreme degree. As in [9], the introduction of an eddy viscosity model greatly improves the performance of the resolvent model (consider the difference between the eddy viscosity and the first-order curves).

Figure 4 shows that the nonlinear terms of equation (1) are better represented overall along the three directions using the second-order model than the first-order one. This result was to be expected as there is no reason to expect that $\underline{\phi}_{\text{1st order}}^{(1)} \approx \bar{\zeta}^{(1)}$ (see Section 4.3). Notice though that the forcing components in the wall-normal and spanwise directions are strongly underpredicted by the new resolvent model. The new model allows again for an accurate restitution of actual flow statistics close to the wall in the streamwise direction. The capability to predict, to some degree, the properties of the nonlinear terms may be useful to construct estimators, as proposed in [10].

### 5.3. Large scale structures

In the following, a case study for ($\lambda_x = 4.19$, $\lambda_z = 1.26$, $c^+ = 16$) is explored. As before, this choice of parameters and the name of this test case correspond to that of [9, paragraph 3.1]. The denomination “large scale” is inferred from the energetically significant dynamics associated with these wavenumbers, which exhibit structures that span the entire channel. The scaling is identical to the previous section, and only the absolute value of the modes is presented. Figure 5 presents the first SPOD velocity and resolvent response modes, whereas Figure 6 presents the first SPOD nonlinear and resolvent forcing modes.

Overall, this case study yields more mixed results. In two of the three directions, the model performs worse than the basic first-order resolvent model or the eddy viscosity system. The new model consistently underestimates real statistics for low $y^+$, and does not improve predictions compared to the simpler first-order model, even in the dominant direction. Worse, the second-order model smoothens out variations of velocity amplitude in the spanwise direction. Here again, the addition of an eddy viscosity yields significant improvements in the predictions from the first-order resolvent model.

Interestingly, the model overestimates the forcing in the $x$ direction in Figure 6 and underestimates the rest. In comparison, an eddy viscosity model yields somewhat relevant predictions in the spanwise direction but otherwise fails to capture the dynamics of the nonlinear terms in the flow. Again, this is to be expected, as optimal forcing of the resolvent matrix cannot be expected to correspond to actual flow statistics.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison between the first SPOD mode of the DNS data $\xi_x^{(1)}$ and predictions from various resolvent models $\psi^{(1)}$. Velocity components are shown for $(Re_\tau, \lambda_x^+, \lambda_z^+, c^+) = (180, 1130, 113, 12)$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison between the first SPOD mode of the DNS data $\zeta^{(1)}$ and prediction from the second-order model $\nabla \cdot \underline{R}_T \phi_{2^{\text{nd}} \text{order}}^{(1)}$. Optimal forcing obtained from first-order resolvents are also shown for reference. Nonlinear and forcing components are shown for $(\mathrm{Re}_T, \lambda_x^+, \lambda_z^+, c^+) = (180, 1130, 113, 12)$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between the first SPOD mode of the DNS data $\xi_y^{(1)}$ and predictions from various resolvent models $\psi^{(1)}$. Velocity components are shown for $(Re_\tau, \lambda_x, \lambda_z, c^+) = (180, 4.19, 1.26, 16)$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison between the first SPOD mode of $\nabla \cdot \underline{\tau}$ in the DNS data $\zeta^{(1)}$ and prediction from second-order resolvent analysis $\nabla \cdot \underline{R_{\tau}}\phi^{(1)}$. Optimal forcing obtained from first-order resolvents are also shown for reference. Nonlinear and forcing components are shown for ($Re_{\tau}$, $\lambda_x$, $\lambda_z$, $c^+$) = (180, 4.19, 1.26, 16).</figcaption></figure>
# 6. Conclusion

The resolvent operator, which relates forcing input of a linear system to the associated response, has been formulated for the Navier-Stokes equations such that the forcing input includes third-order turbulent correlations, which drive a velocity response indirectly via the correlation tensor $\underline{\tau}$. In contrast, classical resolvent analysis, classified as “first-order” in this paper, defines forcing input as a substitute for second-order correlations, which drive the velocity response directly.

This study tested whether the shift of forcing to higher order improves the linear response predictions in the case of turbulent plane channel flow. The relevant quantities to compare are the leading SPOD modes extracted from turbulent DNS data, which represent energetic coherent turbulence structures, and the leading SVD modes from the linear resolvent model. In the present case of a flow with streamwise and spanwise invariant statistics, such modes are obtained for a range of wavenumbers for a logarithmic range of frequencies. Reasoning was done on a half-channel and presented results for an even $u_z$.

Classical first-order resolvent analysis generally yields meaningful predictions of the leading SPOD modes, but with much room for improvement in quantitative accuracy. The limiting assumption of this approach is that the nonlinear terms in the fluctuation equations can be replaced by white noise. The currently best available strategy is to use this assumption in combination with the inclusion of an eddy viscosity in the linear operator. Much recent discussion in the community revolves around the question of how to improve the predictions further by replacing the white noise with “coloured” noise input. The present approach instead aims to produce the “colour” of $\underline{\nabla} \cdot \underline{\tau}$ by computing the $\underline{\tau}$ in a spectral setting as a response to white noise forcing at a higher order.

The success of this strategy, as far as it has been possible to establish, is not as clear as might have been hoped. Streamwise velocity fluctuation amplitudes are very well predicted in a “near-wall” setting in the sense of [7, 9], and the streamwise component of $\underline{\nabla} \cdot \underline{\tau}$ obtained from the model also compares well to the corresponding SPOD mode from the DNS. On the other hand, wall-normal as spanwise velocity and $\underline{\nabla} \cdot \underline{\tau}$ amplitudes are strongly underpredicted. In the “large-scale” setting in the sense of [7, 9], our second-order model clearly gives less accurate results than the standard first-order resolvent, especially when eddy viscosity is included in the latter. This observation stands even if an imperfect model is used for the eddy viscosity, namely one that produces a baseflow that differs slightly from the DNS one. When the model performance is measured in terms of projection with SPOD modes throughout the parameter space, the second-order model cannot compete with the first-order model with eddy viscosity.

The roadblock to accurate second-order resolvent modelling is a typical closure problem, probably caused by the velocity-pressure-gradient correlation detailed in Section 2.1. Efforts to cast this tensor in the form of a linear function of $u$, $p$ or $\underline{\tau}$ so that it could be incorporated into the resolvent operator, have not been fruitful. Yet the present study represents a first step in a new direction for improved resolvent analysis applied to turbulent flow, from where new paths are hoped to emerge.

# Conflicts of interest

The authors declare no competing financial interest.

# Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.
Acknowledgments

Results reflect only the authors' view and the JU or DGA are not responsible for any use that
may be made of the information it contains. The authors convey their thanks to the anonymous
reviewer who significantly contributed to improving our work.

References

[1] Y. Hwang, C. Cossu, “Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow”, *J. Fluid Mech.* **664** (2010), p. 51-73.

[2] A. Towne, O. T. Schmidt, T. Colonius, “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis”, *J. Fluid Mech.* **847** (2018), p. 821-867.

[3] S. B. Pope, *Turbulent Flows*, Cambridge University Press, 2000.

[4] R. Moarref, M. R. Jovanović, J. A. Tropp, A. S. Sharma, B. J. McKeon, “A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization”, *Phys. Fluids* **26** (2014), no. 5, article no. 051701.

[5] M. A. Ahmed, H. J. Bae, A. F. Thompson, B. J. McKeon, “Resolvent analysis of stratification effects on wall-bounded shear flows”, *Phys. Rev. Fluids* **6** (2021), no. 8, article no. 084804.

[6] K. Rosenberg, B. J. McKeon, “Computing exact coherent states in channels starting from the laminar profile: A resolvent-based approach”, *Phys. Rev. E* **100** (2019), no. 2, article no. 021101.

[7] L. I. Abreu, A. V. G. Cavalieri, P. Schlatter, R. Vinuesa, D. S. Henningson, “Resolvent modelling of near-wall coherent structures in turbulent channel flow”, *Int. J. Heat Fluid Flow* **85** (2020), article no. 108662.

[8] P. Morra, O. Semeraro, D. S. Henningson, C. Cossu, “On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows”, *J. Fluid Mech.* **867** (2019), p. 969-984.

[9] P. Morra, P. A. S. Nogueira, A. V. G. Cavalieri, D. S. Henningson, “The colour of forcing statistics in resolvent analyses of turbulent channel flows”, *J. Fluid Mech.* **907** (2021), article no. a24.

[10] F. R. Amaral, A. V. G. Cavalieri, E. Martini, P. Jordan, A. Towne, “Resolvent-based estimation of turbulent channel flow using wall measurements”, *J. Fluid Mech.* **927** (2021), article no. a17.

[11] S. Symon, S. J. Illingworth, I. Marusic, “Energy transfer in turbulent channel flows and implications for resolvent modelling”, *J. Fluid Mech.* **911** (2021), article no. a3.

[12] S. Symon, A. Madhusudanan, S. J. Illingworth, I. Marusic, “On the use of eddy viscosity in resolvent analysis of turbulent channel flow”, preprint, arXiv:2205.11216, 2022.

[13] B. F. Farrell, D. F. Gayme, P. J. Ioannou, “A statistical state dynamics approach to wall turbulence”, *Philos. Trans. R. Soc. Lond., Ser. A* **375** (2017), no. 2089, article no. 20160081.

[14] J. B. Marston, G. P. Chini, S. M. Tobias, “Generalized Quasilinear Approximation: Application to Zonal Jets”, *Phys. Rev. Lett.* **116** (2016), no. 21, article no. 214501.

[15] S. C. Reddy, J. A. C. Weideman, “The Accuracy of the Chebyshev Differencing Method for Analytic Functions”, *SIAM J. Numer. Anal.* **42** (2005), no. 5, p. 2176-2187.

[16] O. T. Schmidt, T. Colonius, “Guide to Spectral Proper Orthogonal Decomposition”, *AIAA J.* **58** (2020), no. 3, p. 1023-1033.

[17] E. Martini, A. V. G. Cavalieri, P. Jordan, L. Lesshafft, “Accurate Frequency Domain Identification of ODEs with Arbitrary Signals”, preprint, arXiv:1907.04787, 2020.

[18] W. C. Reynolds, W. G. Tiederman, “Stability of turbulent channel flow, with application to Malkus's theory”, *J. Fluid Mech.* **27** (1967), no. 2, p. 253-272.

[19] A. Towne, T. Colonius, P. Jordan, A. V. G. Cavalieri, G. A. Brès, “Stochastic and nonlinear forcing of wavepackets in a Mach 0.9 jet”, in *21st AIAA/CEAS Aeroacoustics Conference*, American Institute of Aeronautics and Astronautics, Inc., 2015.

[20] U. Karban, E. Martini, A. V. G. Cavalieri, L. Lesshafft, P. Jordan, “Self-similar mechanisms in wall turbulence studied using resolvent analysis”, *J. Fluid Mech.* **939** (2022), article no. A36.