1 Introduction
Turbulent flows are encountered in a wide variety of industrial applications, and they often present low-frequency unsteadiness. One of the main goals is the prediction of these unsteady features, such as the characteristic frequencies or the spatial structure of the unsteadiness. A substantial quantity of research has been devoted to the use of stability theory to achieve this goal. The classical linear stability theory around a steady base flow fails to do so, by not accounting for the nonlinearities that have a strong influence on the unsteady behaviour of the flow (Barkley Reference Barkley2006). However, there are many papers (Mattingly & Criminale Reference Mattingly and Criminale1972; Pier Reference Pier2002; Mittal Reference Mittal2008, for instance) showing that those nonlinear effects may be factored in by using the time-averaged field (the mean flow) instead of the base flow. For example, in the cylinder case, a linear stability analysis around the mean flow gives an excellent prediction of the vortex shedding frequency, even far from criticality (Pier Reference Pier2002; Barkley Reference Barkley2006).
The same mean flow stability approach has been successfully used for flows exhibiting a broadband spectrum originating from convective instabilities. For instance, Gudmundsson & Colonius (Reference Gudmundsson and Colonius2011) focused on turbulent round jets: using an array of microphones, they experimentally measured the pressure fluctuations outside the jet shear layer and showed that the pressure amplitude and phase of the data accurately matched predictions from a local stability analysis around the mean flow. The same conclusions were obtained by Oberleithner, Rukes & Soria (Reference Oberleithner, Rukes and Soria2014) in the case of a transitional jet, where the mean flow and unsteady structures were fully characterised by particle image velocimetry measurements: at given frequencies, a local stability analysis around the mean velocity profiles led to the full reconstruction of the spatial structure (in magnitude and phase) of the perturbation field.
These examples are just a part of the large number of studies that demonstrate the capability of a mean flow stability analysis to predict the unsteady features of a transitional or turbulent flow. On the other hand, a formal justification for the well foundedness of such an approach is still wanting, and only a few studies have been carried out to address this question. Among those, Sipp & Lebedev (Reference Sipp and Lebedev2007) focused on self-excited systems which present a strong dominating frequency. They showed that, in the vicinity of the bifurcation threshold, if the mean flow harmonic dominates the second harmonic, a global linear stability analysis around the mean flow yields a marginally stable mode whose frequency matches the natural frequency of the flow. Recently, Turton, Tuckerman & Barkley (Reference Turton, Tuckerman and Barkley2015) more generally demonstrated that if the flow exhibits monochromatic harmonic oscillations, then the linearised operator around the mean flow indeed exhibits a purely imaginary eigenvalue equal to the frequency of the flow. This criterion has even been considered by Mantič-Lugo, Arratia & Gallaire (Reference Mantič-Lugo, Arratia and Gallaire2014, Reference Mantič-Lugo, Arratia and Gallaire2015) to build a self-consistent model of cylinder flow that predicts the frequency of the vortex shedding for Reynolds numbers up to 110. But the conclusions of these studies cannot be extended to the more general case where a flow presents a broadband spectrum. On account of this, Mantič-Lugo (Reference Mantič-Lugo2015) studied the prediction of the nonlinear dynamics of a laminar backward facing step by a self-consistent model that predicts the saturation quite well.
In this article, we aim at providing, in the general case of a flow field that presents either a broadband or a peaked spectrum, mathematical insights to justify the efficiency of a mean flow stability approach to predict the spatio-temporal features of a flow field. Our objective is not to tackle the turbulence closure problem, as in the stochastic structural stability theory (SSST) of Farrell & Ioannou (Reference Farrell and Ioannou2012): the mean flow is here considered as known a priori (for example, from experimental measurements) and we aim at extracting the spatial structures of the underlying fluctuation field (the coherent structures) that are consistent with this mean flow. The theoretical developments are expressed in a general framework, without consideration of the turbulent/laminar or amplifier/oscillator nature of the flow, and are based on the optimal response modes computed from a singular value decomposition (SVD) of the resolvent operator around the mean flow. The use of the SVD in hydrodynamic stability has first been considered as a way to study energy amplification due to non-modal mechanisms in transitional flows (Trefethen et al. Reference Trefethen, Trefethen, Reddy, Driscoll and others1993; Schmid & Henningson Reference Schmid and Henningson2012). Non-modal analysis with turbulent mean flows has been introduced by Butler & Farrell (Reference Butler and Farrell1993), Chernyshenko & Baig (Reference Chernyshenko and Baig2005), Del Álamo & Jimenez (Reference Del Álamo and Jimenez2006), Cossu, Pujals & Depardon (Reference Cossu, Pujals and Depardon2009), Pujals et al. (Reference Pujals, García-Villalba, Cossu and Depardon2009) and Moarref & Jovanović (Reference Moarref and Jovanović2012). McKeon & Sharma (Reference McKeon and Sharma2010) showed that the first optimal response computed from an SVD around the mean flow dominates the velocity field of the full flow field when the first singular value is significantly larger than all others. However, they restricted their analysis to a turbulent pipe configuration, which is invariant (homogeneous) in the streamwise direction. Hence, the mean flow was constant in the streamwise direction and all fluctuating quantities were Fourier transformed in this direction. Therefore, only a local stability analysis was required. The present study is an extension of their work, where we consider more general configurations, in particular open flows which are not invariant in the streamwise direction, such as a backward facing step or jet configurations. In such cases, a global stability analysis is required, in which all operators and variables depend on the streamwise coordinate in a non-specified (general) manner. As a consequence, the convective non-normality (Marquet et al. Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009) becomes a dominant physical mechanism, which is responsible for the downstream location of the optimal response modes and the upstream location of the optimal forcing modes. Finally, we also aim at elucidating the link between the global stability results and those provided by local stability approaches, such as spatial stability or PSE (parabolised stability equations, see Herbert Reference Herbert1997) analyses.
The present article is divided into four sections. In § 2, we highlight the role of the resolvent operator around the mean flow for the determination of the temporal spectral (Fourier) mode of the fluctuation field. Then, we introduce a rank 1 approximation of the resolvent operator to link this spectral mode to the dominant optimal response, and relate this result to local and global stability analyses (§ 3). This relation yields models for the prediction of the frequency spectrum at arbitrary points of a flow (§ 4). The last part is dedicated to an application of these results to a turbulent backward facing step.
2 Resolvent-based equation for the Fourier mode of the fluctuation field
We consider the incompressible, homogeneous Navier–Stokes equations:
$$\begin{eqnarray}\displaystyle \partial _{t}\boldsymbol{u}=-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}-\boldsymbol{{\rm\nabla}}p,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & & \displaystyle\end{eqnarray}$$
with
$\boldsymbol{u}$
,
${\it\nu}$
,
$p$
respectively standing for the velocity, kinematic viscosity and pressure of the fluid. If we consider a Reynolds decomposition of the flow variables, the mean quantities (denoted by an overline) verify:
$$\begin{eqnarray}\displaystyle -\overline{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\overline{\boldsymbol{u}}-\overline{\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }}+{\it\nu}{\rm\nabla}^{2}\overline{\boldsymbol{u}}-\boldsymbol{{\rm\nabla}}\overline{p}=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\overline{\boldsymbol{u}}=0, & & \displaystyle\end{eqnarray}$$
and the fluctuating quantities, denoted by a prime:
$$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}\boldsymbol{u}^{\prime }=-\overline{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\overline{\boldsymbol{u}}+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}^{\prime }-\boldsymbol{{\rm\nabla}}p^{\prime }-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }+\overline{\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0. & \displaystyle\end{eqnarray}$$
$\boldsymbol{f}^{\prime }=\overline{\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }}-\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\prime }$
, and the linear operators
$L$
,
$B$
and
$P$
by: 
$$\begin{eqnarray}\displaystyle L=\left(\begin{array}{@{}cc@{}}-\overline{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}()-()\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\overline{\boldsymbol{u}}+{\it\nu}{\rm\nabla}^{2}() & -\boldsymbol{{\rm\nabla}}()\\ \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }() & 0\end{array}\right),\quad B=\left(\begin{array}{@{}cc@{}}1 & 0\\ 0 & 0\end{array}\right),\quad P=\left(\begin{array}{@{}c@{}}1\\ 0\end{array}\right), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
equations (2.3), (2.4) can be rewritten in the compact form:
$$\begin{eqnarray}\displaystyle B\partial _{t}\left(\begin{array}{@{}c@{}}\boldsymbol{u}^{\prime }\\ p^{\prime }\end{array}\right)=L\left(\begin{array}{@{}c@{}}\boldsymbol{u}^{\prime }\\ p^{\prime }\end{array}\right)+P\boldsymbol{f}^{\prime }. & & \displaystyle\end{eqnarray}$$
The temporal Fourier transform of (2.6) gives, for each frequency
${\it\omega}$
:
$$\begin{eqnarray}\displaystyle \text{i}{\it\omega}B\left(\begin{array}{@{}c@{}}\hat{\boldsymbol{u}}\\ \hat{p}\end{array}\right)=L\left(\begin{array}{@{}c@{}}\hat{\boldsymbol{u}}\\ \hat{p}\end{array}\right)+P\hat{\boldsymbol{f}}, & & \displaystyle\end{eqnarray}$$
with
$\hat{\boldsymbol{u}}$
,
$\hat{p}$
,
$\hat{\boldsymbol{f}}$
as respectively the Fourier transform of
$\boldsymbol{u}^{\prime }$
,
$p^{\prime }$
and
$\boldsymbol{f}^{\prime }$
. Finally, by using the restriction operator
$R=(1~0)$
to eliminate the pressure term, (2.7) can be rearranged as:
$$\begin{eqnarray}\hat{\boldsymbol{u}}=\mathscr{R}({\it\omega})\hat{\boldsymbol{f}},\end{eqnarray}$$
where
$\mathscr{R}({\it\omega})=R(\text{i}{\it\omega}B-L)^{-1}P$
is the resolvent around the mean flow of the system (Sipp et al.
Reference Sipp, Marquet, Meliga and Barbagallo2010). This equation is the cornerstone of the analysis that will follow. It highlights the central role that the linear operator
$\mathscr{R}$
plays in the natural dynamics of a flow. Note that
$\hat{\boldsymbol{f}}$
is not an external forcing, and cannot be seen as such, since it depends on the perturbation field. In this study, we avoid any closure problem by assuming that the mean flow is known a priori, for example by experimental data (similarly to McKeon & Sharma (Reference McKeon and Sharma2010)). In this situation, (2.8) can be seen as a linear system governed by
$\mathscr{R}$
relating the spectral mode
$\hat{\boldsymbol{u}}$
to a specific forcing
$\hat{\boldsymbol{f}}$
.
3 Relation between the spectral mode and the dominant optimal response
This section focuses on the study of the linear operator
$\mathscr{R}$
introduced in the previous section. For a given harmonic forcing
$\tilde{{\bf\phi}}$
and corresponding harmonic response
$\tilde{{\bf\psi}}=\mathscr{R}\tilde{{\bf\phi}}$
, we define the gain function
${\it\mu}^{2}(\tilde{{\bf\phi}})$
as:
$$\begin{eqnarray}{\it\mu}^{2}(\tilde{{\bf\phi}})=\frac{\Vert \tilde{{\bf\psi}}\Vert ^{2}}{\Vert \tilde{{\bf\phi}}\Vert ^{2}}=\frac{\langle \mathscr{R}\tilde{{\bf\phi}},\mathscr{R}\tilde{{\bf\phi}}\rangle }{\langle \tilde{{\bf\phi}},\tilde{{\bf\phi}}\rangle }=\frac{\langle \mathscr{R}^{\dagger }\mathscr{R}\tilde{{\bf\phi}},\tilde{{\bf\phi}}\rangle }{\langle \tilde{{\bf\phi}},\tilde{{\bf\phi}}\rangle },\end{eqnarray}$$
with
$\Vert \cdot \Vert$
a relevant Hermitian norm,
$\langle \cdot ,\cdot \rangle$
the corresponding scalar product and
$\mathscr{R}^{\dagger }$
the adjoint of
$\mathscr{R}$
, defined as the operator
$\mathscr{R}^{\dagger }$
that verifies
$\forall \boldsymbol{a},\boldsymbol{b};~\langle \boldsymbol{a},\mathscr{R}\boldsymbol{b}\rangle =\langle \mathscr{R}^{\dagger }\boldsymbol{a},\boldsymbol{b}\rangle$
. For example, we may consider the energy norm defined as
$\Vert \tilde{{\bf\psi}}\Vert ^{2}=\langle \tilde{{\bf\psi}},\tilde{{\bf\psi}}\rangle =\int (\tilde{{\bf\psi}}\boldsymbol{\cdot }\tilde{{\bf\psi}})\,\text{d}\boldsymbol{x}$
, where
$()\boldsymbol{\cdot }()$
is the local Hermitian inner product.
The set of optimal forcings of unit norm
$({\bf\phi}_{j},j\geqslant 1)$
, that maximise
${\it\mu}^{2}$
, can be obtained by solving
$\mathscr{R}^{\dagger }\mathscr{R}{\bf\phi}_{j}={\it\mu}_{j}^{2}{\bf\phi}_{j}$
. The solution gain values
$({\it\mu}_{j}^{2},j\geqslant 1)$
and the corresponding forcings
$({\bf\phi}_{j},j\geqslant 1)$
are sorted such that
$({\it\mu}_{1}^{2}\geqslant {\it\mu}_{2}^{2}\geqslant \cdots \,)$
, and so the optimal forcing
${\bf\phi}_{1}$
corresponds to the harmonic forcing which verifies
${\it\mu}_{1}^{2}={\it\mu}^{2}({\bf\phi}_{1})=\sup _{\tilde{{\bf\phi}}}{\it\mu}^{2}(\tilde{{\bf\phi}})$
.
Note that the operator
$\mathscr{R}^{\dagger }\mathscr{R}$
is Hermitian, consequently
$\langle {\bf\phi}_{i},{\bf\phi}_{j}\rangle ={\it\delta}_{ij}$
, meaning that
$({\bf\phi}_{j},j\geqslant 1)$
is an orthonormal basis of the forcing space. The projection of
$\hat{\boldsymbol{f}}$
in (2.8) onto this basis then yields:
$$\begin{eqnarray}\hat{\boldsymbol{u}}=\mathscr{R}({\it\omega})\mathop{\sum }_{j\geqslant 1}\langle {\bf\phi}_{j},\hat{\boldsymbol{f}}\rangle {\bf\phi}_{j}.\end{eqnarray}$$
For each optimal forcing
${\bf\phi}_{j}$
, the corresponding optimal response of unit norm
${\bf\psi}_{j}$
is defined by
${\bf\psi}_{j}={\it\mu}_{j}^{-1}\mathscr{R}({\it\omega}){\bf\phi}_{j}$
. Similarly to the set of optimal forcings, the set of optimal responses of unit norm
$({\bf\psi}_{j},j\geqslant 1)$
defines an orthonormal basis of the response space (
$\langle {\bf\psi}_{i},{\bf\psi}_{j}\rangle ={\it\delta}_{ij}$
). Using this basis, (3.2) can be rewritten as:
$$\begin{eqnarray}\hat{\boldsymbol{u}}=\mathop{\sum }_{j\geqslant 1}{\bf\psi}_{j}{\it\mu}_{j}\langle {\bf\phi}_{j},\hat{\boldsymbol{f}}\rangle .\end{eqnarray}$$
This final relation is related to the singular value decomposition (SVD) of the operator
$\mathscr{R}$
, with
$({\it\mu}_{j},j\geqslant 1)$
the so-called singular values of
$\mathscr{R}$
, and
$({\bf\psi}_{j},j\geqslant 1)$
and
$({\bf\phi}_{j},j\geqslant 1)$
respectively the left and right singular vectors of
$\mathscr{R}$
. In the case where, for a given frequency
${\it\omega}$
, the first singular value
${\it\mu}_{1}({\it\omega})$
is significantly larger than the others, a rank 1 approximation of (3.3) yields:
$$\begin{eqnarray}\hat{\boldsymbol{u}}(\boldsymbol{x},{\it\omega})\approx {\bf\psi}_{1}(\boldsymbol{x},{\it\omega}){\it\mu}_{1}({\it\omega})\langle {\bf\phi}_{1}({\it\omega}),\hat{\boldsymbol{f}}({\it\omega})\rangle .\end{eqnarray}$$
The approximation error reads:
$$\begin{eqnarray}\frac{\Vert \hat{\boldsymbol{u}}-{\bf\psi}_{1}{\it\mu}_{1}\langle {\bf\phi}_{1},\hat{\boldsymbol{f}}\rangle \Vert ^{2}}{\Vert \hat{\boldsymbol{u}}\Vert ^{2}}=1-\frac{{\it\mu}_{1}^{2}|\langle {\bf\phi}_{1},\hat{\boldsymbol{f}}\rangle |^{2}}{\displaystyle \mathop{\sum }_{i\geqslant 1}{\it\mu}_{i}^{2}|\langle {\bf\phi}_{i},\hat{\boldsymbol{f}}\rangle |^{2}}.\end{eqnarray}$$
The condition
${\it\mu}_{1}\gg {\it\mu}_{j\geqslant 2}$
is therefore not a strictly sufficient condition for the approximation to hold. A more rigorous one is:
$$\begin{eqnarray}{\it\mu}_{1}^{2}|\langle {\bf\phi}_{1},\hat{\boldsymbol{f}}\rangle |^{2}\gg \mathop{\sum }_{i\geqslant 2}{\it\mu}_{i}^{2}|\langle {\bf\phi}_{i},\hat{\boldsymbol{f}}\rangle |^{2}.\end{eqnarray}$$
However, under the reasonable assumption that
$\hat{\boldsymbol{f}}$
does not display any preferential direction toward one of the suboptimal forcings
${\bf\phi}_{j\geqslant 2}$
, this stricter condition reduces to
${\it\mu}_{1}\gg {\it\mu}_{j\geqslant 2}$
. Note that the accuracy of the approximation is ensured in a global sense (the error is evaluated with a 2-norm), which yields a good approximation in the high-energy regions. In contrast, this might be locally inaccurate, especially in low-energy regions, which are hopefully, in most applications, of lesser interest. Note that building approximations that hold at any point requires
$\infty$
-norm-based approaches, which have been considered in hydrodynamic stability only very recently (Foures, Caulfield & Schmid Reference Foures, Caulfield and Schmid2013).
We have explicitly indicated the variable dependency of each term of (3.4), highlighting that the only term on the right-hand side depending on the spatial position
$\boldsymbol{x}$
is
${\bf\psi}_{1}$
, since
${\it\mu}_{1}$
is a function of the frequency only, and the result of the scalar product
$\langle {\bf\phi}_{1}({\it\omega}),\hat{\boldsymbol{f}}({\it\omega})\rangle$
is space independent. This implies that for the frequencies where the dominant singular value hypothesis holds,
$\hat{\boldsymbol{u}}$
is directly proportional to the dominant optimal response
${\bf\psi}_{1}$
. In other words, this result states that the spatial structure of the unsteadiness at these frequencies may be closely related to the dominant singular mode around the mean flow. From the literature (for instance Dergham, Sipp & Robinet Reference Dergham, Sipp and Robinet2013), the dominant singular value hypothesis seems to be fulfilled when a flow presents one dominant convective instability mechanism (for example, the Kelvin–Helmholtz mechanism in a shear layer). This point is discussed in further detail in § 5.3, and the rest of the present study is dedicated to such flows. This excludes flows displaying a range of frequencies for which the dominant singular value hypothesis is not verified, due, for example, to the existence of several strong instability mechanisms. These cases could be treated by a higher rank approximation over this particular frequency range.
The present results justify why a mean flow stability analysis may accurately predict unsteady features of a flow. The findings of the present study can also relate to the global stability modes (eigenvectors of the linearised Navier–Stokes operator) around a mean flow. Turton et al. (Reference Turton, Tuckerman and Barkley2015) showed that if the flow field is monochromatic, then the mean flow exhibits a marginal global mode at the frequency of the flow. Indeed, for the frequency response in (2.8) to reproduce a monochromatic tone at a frequency
${\it\omega}_{0}$
, there is only one possibility: the resolvent operator should display a very strong response in the very vicinity of this frequency, which requires that an eigenvalue of the Jacobian is close to
$\text{i}{\it\omega}_{0}$
. We are therefore led to the same conclusion as Turton et al. (Reference Turton, Tuckerman and Barkley2015), which appears as a particular case of the more general approach presented here. Note that the reciprocal also holds: if a quasi-marginal eigenvalue exists in the vicinity of some frequency, then the first optimal gain function
${\it\mu}_{1}^{2}$
will become extremely strong in the vicinity of this frequency. Therefore, a clear separation of the singular values near this frequency is expected, which would justify the validity of a mean flow stability analysis.
4 Predictive models for the local frequency spectra
Equation (3.4) is a function of the frequency and of space. Predicting the frequency spectrum of the flow at different points requires the determination of all the terms of the equation, including
$\hat{\boldsymbol{f}}$
. Two models are proposed in this section, which will provide a prediction of the frequency spectrum. Defining
${\it\Lambda}({\it\omega})\equiv {\it\mu}_{1}\langle {\bf\phi}_{1},\hat{\boldsymbol{f}}\rangle$
, (3.4) is re-written:
$$\begin{eqnarray}\hat{\boldsymbol{u}}(x,{\it\omega})\approx {\it\Lambda}({\it\omega}){\bf\psi}_{1}(x,{\it\omega}).\end{eqnarray}$$
The function
${\it\Lambda}$
is directly related to the energy spectral density
$e({\it\omega})$
of the fluctuation field integrated over the whole domain
${\it\Omega}$
. Indeed, Parseval’s theorem yields:
$$\begin{eqnarray}\frac{1}{2}\int _{{\it\Omega}}\left(\int _{0}^{\infty }\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime }\,\text{d}t\right)\,\text{d}\boldsymbol{x}=\frac{1}{4{\rm\pi}}\int _{{\it\Omega}}\left(\int _{0}^{\infty }\hat{\boldsymbol{u}}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\,\text{d}{\it\omega}\right)\,\text{d}\boldsymbol{x}=\int _{0}^{\infty }\frac{\Vert \hat{\boldsymbol{u}}\Vert ^{2}}{4{\rm\pi}}\,\text{d}{\it\omega}.\end{eqnarray}$$
Consequently,
$e({\it\omega})=\Vert \hat{\boldsymbol{u}}({\it\omega})\Vert ^{2}/4{\rm\pi}$
, and since
${\bf\psi}_{1}$
is of unit norm, (4.1) gives:
$$\begin{eqnarray}e({\it\omega})=\frac{1}{4{\rm\pi}}\Vert \hat{\boldsymbol{u}}({\it\omega})\Vert ^{2}\approx \frac{1}{4{\rm\pi}}|{\it\Lambda}({\it\omega})|^{2}.\end{eqnarray}$$
A first simple assumption consists in considering that
$e$
is proportional to the dominant optimal energetic gain function
${\it\mu}_{1}^{2}$
, which can be achieved by posing
${\it\Lambda}={\it\mu}_{1}$
(model 1). Note that from the definition of
${\it\Lambda}$
, this implies that
$\langle {\bf\phi}_{1},\hat{\boldsymbol{f}}\rangle =1$
for all
${\it\omega}$
, which is a very strong assumption since the projection of the turbulent forcing term on
${\bf\phi}_{1}$
must depend on the frequency and is not equal to 1 in general. Hence, we expect that model 1 may provide (at best) only a qualitative shape of the frequency spectrum. Moreover, this model presents one strong limitation that stems from considering
${\it\Lambda}$
as a positive real-valued function (the actual function
${\it\Lambda}$
is complex valued). By doing so, we lose all phase-related information, and can only predict the amplitude of the spectra. The main advantage of this model is that it does not require any knowledge of
$\hat{\boldsymbol{u}}$
. Results in § 5.4 confirm that this model provides good qualitative insights.
A more precise and rigorous model can be built by assuming a partial local knowledge of
$\hat{\boldsymbol{u}}$
(model 2). For example, from (4.1), the function
${\it\Lambda}$
may be defined using one component
$\hat{u} _{i}$
of
$\hat{\boldsymbol{u}}$
, as
${\it\Lambda}({\it\omega})=\hat{u}_{i}(\boldsymbol{x},{\it\omega})/{\it\psi}_{1,i}(\boldsymbol{x},{\it\omega})$
. Considering that
$\hat{u} _{i}({\it\omega})$
is known at a given point
$\boldsymbol{x}_{0}$
yields
${\it\Lambda}({\it\omega})=\hat{u}_{i}(\boldsymbol{x}_{0},{\it\omega})/{\it\psi}_{1,i}(\boldsymbol{x}_{0},{\it\omega})$
. This approach consists in rescaling the response mode for each frequency, such that its amplitude and phase match those of
$\hat{\boldsymbol{u}}$
at a given point. Since the spatial shape of this response mode and
$\hat{\boldsymbol{u}}$
are close, this provides a reconstruction of the spectra over the whole domain. Contrary to model 1 that only predicts a qualitative shape of the spectra, model 2 provides a quantitative prediction of the amplitude, but requires the knowledge of the spectrum at at least one point
$\boldsymbol{x}_{0}$
. The issue of the choice of
$\boldsymbol{x}_{0}$
is addressed in § 5.4, where the two models are compared in the turbulent backward facing step case.
5 Application to a turbulent backward facing case
5.1 Computation of the mean flow
We consider the case of a two-dimensional (2-D) backward facing step. From now on, we consider all the quantities to be non-dimensionalised by the step height
$h$
and the free-stream velocity
$u_{\infty }$
. We performed an unsteady three-dimensional (3-D) simulation by using the compressible flow solver ElsA (Cambier, Heib & Plot Reference Cambier, Heib and Plot2013), at low Mach number
$M=0.09$
(the flow field may therefore be considered as homogeneous and incompressible) and
$Re=57\,500$
. These parameters correspond to those of an experiment that has been carried in ONERA wind tunnels (Gallas et al.
Reference Gallas, Arnault, Monnier, Delva, Dandois and Mialon2015). The solid walls of the step are modelled by adiabatic no-slip conditions. We impose a fluid injection upstream of the step on the left boundary of the computational domain, such that the boundary layer thickness is
${\it\delta}=0.35$
at the step location (the momentum thickness is
${\it\theta}=0.041$
and the displacement thickness is
${\it\delta}^{\ast }=0.056$
). Periodic boundary conditions are used in the spanwise direction. All other boundary conditions are treated by imposing an output pressure. The computational domain dimensions are 38 in length (23 downstream from the step edge, 15 upstream), 4 in width and 10 in height (downstream from the step). The mesh contains approximately 19 million cells (block-structured mesh,
$L\times H\times W$
:
$597\times 185\times 161$
for the downstream block,
$x>0$
, and
$91\times 79\times 161$
for the upstream block,
$x<0$
), and satisfies
$y^{+}\leqslant 1$
at the wall boundaries. We then performed a zonal detached eddy simulation (ZDES) (Deck Reference Deck2005) for which the region upstream of the step edge is described by a steady Reynolds-averaged Navier–Stokes model, and the region downstream by an unsteady large eddy simulation model. A 2-D mean flow was obtained by averaging the results over 310 convection time units, and over the spanwise direction.
To assess the validity of the mean flow, three methods have been used to compute the streamwise reattachment location
$x_{r}$
: (a) the location where the mean streamwise velocity
$U$
is equal to 0 at the first cell away from the wall, (b) the location of the dividing streamline of the mean flow and finally (c) the location where the wall shear stress
$\partial U/\partial y=0$
. These three methods yield
$x_{r}=6.45\pm 0.2\,\%$
, which corresponds to the results from the literature for comparable Reynolds number and boundary layer thickness (Driver, Seegmiller & Marvin Reference Driver, Seegmiller and Marvin1987; Adams & Johnston Reference Adams and Johnston1988, for instance). The results of (a) can be seen in figure 1. A secondary reattachment point is located at
$x=0.88$
, which is consistent with the work of Spazzini et al. (Reference Spazzini, Iuso, Onorato, Zurlo and Di Cicca2001), who found a secondary separation at
$x=1.1$
for a Reynolds number of 16 000. They also showed that this length tends to decrease as the Reynolds number increases. Figure 1 displays a third corner recirculation region near
$(x,y)=(0,0)$
, with a size of 0.053. As a comparison, this value is close to that found by Le, Moin & Kim (Reference Le, Moin and Kim1997) at a lower Reynolds number (0.042 for
$Re=5100$
).
The unsteady features computed from the simulation also present a good agreement with the literature. The spectrum at
$(x=2,y=1.5)$
is shown in figure 2(a), and displays a clear peak at approximately
$St_{{\it\theta}}=0.013$
(
${\it\omega}=2.1$
), where
$St_{{\it\theta}}$
is the Strouhal number based on the free-stream velocity and the momentum thickness at the step location. This is the frequency expected for the Kelvin–Helmholtz instability (see for example Chun & Sung Reference Chun and Sung1996). Further downstream, near to the reattachment, the flow presents another energetic frequency. The associated Strouhal number
$F^{+}$
based on the free-stream velocity and the reattachment length is approximately 0.65 (
${\it\omega}=0.65$
) (see figure 2
b), which is consistent with the values that can be found in other studies (values between 0.5 and 1, see for example Dandois, Garnier & Sagaut Reference Dandois, Garnier and Sagaut2007).
Figure 1. Streamwise velocity of the mean flow, computed by averaging in time and in the spanwise direction the results of the 3-D unsteady simulation. Contours are in equal increments from
$-$
0.2 to 1. Negative contours are dashed, the thickened contours correspond to a zero velocity.
Figure 2. Spectra of the streamwise (spanwise-averaged) velocity computed from the unsteady ZDES simulation at two locations: (a)
$x=2$
,
$y=1.5$
, (b)
$x=7.0$
,
$y=0.1$
.
$St_{{\it\theta}}$
is the Strouhal number based on the free-stream velocity and the momentum thickness at the step location, and
$F^{+}$
the Strouhal number based on the free-stream velocity and the reattachment length.
5.2 Fluctuation field versus first optimal response
Since the configuration is invariant in the spanwise direction, the fluctuation field can be expressed via its spatio-temporal Fourier coefficients, depending on the frequency
${\it\omega}$
and on the spanwise number
${\it\beta}$
. The resolvent operator therefore also depends on
${\it\beta}$
and we expect that the Fourier mode
$\hat{\boldsymbol{u}}({\it\omega},{\it\beta})$
is proportional to
${\bf\psi}_{1}({\it\omega},{\it\beta})$
if
${\it\mu}_{1}({\it\omega},{\it\beta})\gg {\it\mu}_{2}({\it\omega},{\it\beta})$
. In the following, we will focus on 2-D perturbations (
${\it\beta}=0$
), for which the spanwise component
$u_{3}^{\prime }$
is null. The Jacobian matrix for 2-D perturbations around the 2-D mean flow has been explicitly computed by finite differences, following the same procedure as in Mettot, Sipp & Bézard (Reference Mettot, Sipp and Bézard2014). An SVD of the resolvent matrix is performed using the ARPACK library (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998), coupled with the parallel MUMPS solver (Amestoy et al.
Reference Amestoy, Duff, L’Excellent and Koster2001). Figure 3 presents the three highest singular values over a wide range of frequencies. It shows that the first singular value
${\it\mu}_{1}$
is several orders of magnitude larger than all others, as expected from a flow presenting a strong convective instability mechanism (see § 5.3). Consequently, following the results of § 3, at all these frequencies, the spectral mode
$\hat{\boldsymbol{u}}$
is expected to be proportional to the dominant optimal response
${\bf\psi}_{1}$
.
Figure 3. Comparison of the first three energetic gain functions in the case of two-dimensional perturbations (
${\it\beta}=0$
): one can see that the first singular value
${\it\mu}_{1}$
is significantly higher than the others, which ensures the validity of the leading singular value assumption (§ 3).
Checking this assertion requires the computation of
$\hat{\boldsymbol{u}}$
from the simulation. To obtain smooth spectral modes, we divided the time series of spanwise-averaged velocity snapshots
$\boldsymbol{u}^{\prime }(x,y,t,{\it\beta}=0)$
, computed from the simulation, into
$N_{b}$
bins. For the present study,
$N_{b}=30$
, with a 50 % overlap between the bins. Each bin has a frequency resolution of
${\rm\Delta}{\it\omega}={\rm\pi}/10$
. For each bin
$k$
, the associated spectrum
$\hat{\boldsymbol{u}}^{k}(x,y,{\it\omega})$
has been computed by fast Fourier transform. The final spectrum may be obtained by two approaches: (i) the amplitude of the spectrum is computed from the root mean square (r.m.s.) of all the bins, which corresponds to the classical Welch algorithm:
$$\begin{eqnarray}|\hat{\boldsymbol{u}}_{welch}|=\left(\frac{1}{N_{b}}\mathop{\sum }_{k=1}^{N_{b}}|\hat{\boldsymbol{u}}^{k}|^{2}\right)^{1/2};\end{eqnarray}$$
(ii) the amplitude and phase of the spectrum are computed by performing a proper orthogonal decomposition (POD) based filtering of these bins. Method (ii) aims at isolating the most energetic, spatially correlated spectral field. It has been used for instance by Gudmundsson & Colonius (Reference Gudmundsson and Colonius2011), and consists of computing, for each bin
$k$
and for each frequency
${\it\omega}$
, the cross-spectral tensor
$\unicode[STIX]{x1D64F}^{(k,{\it\omega})}$
defined as:
$$\begin{eqnarray}\unicode[STIX]{x1D61B}_{ij}^{(k,{\it\omega})}=\hat{\boldsymbol{u}}^{k}({\it\omega},\boldsymbol{x}_{i})\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{k}({\it\omega},\boldsymbol{x}_{j}),\end{eqnarray}$$
where
$\boldsymbol{x}_{i,j}$
are the discretisation points. We then compute the ensemble-averaged tensor
$\unicode[STIX]{x1D64F}_{0}^{{\it\omega}}$
:
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{0}^{{\it\omega}}=\frac{1}{N_{b}}\mathop{\sum }_{k=1}^{N_{b}}\unicode[STIX]{x1D64F}^{(k,{\it\omega})},\end{eqnarray}$$
and finally solve the eigenvalue problem:
$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{0}^{{\it\omega}}\boldsymbol{y}^{{\it\omega}}={\it\lambda}\boldsymbol{y}^{{\it\omega}}.\end{eqnarray}$$
By considering the unit norm eigenvector
$\boldsymbol{y}_{m}^{{\it\omega}}$
associated with the highest eigenvalue, we obtain the POD-filtered spectral mode
$\hat{\boldsymbol{u}}_{pf}({\it\omega})=A_{{\it\omega}}\boldsymbol{y}_{m}^{{\it\omega}}$
, where the amplitude
$A_{{\it\omega}}$
is defined as:
$$\begin{eqnarray}A_{{\it\omega}}=\frac{1}{N_{b}}\mathop{\sum }_{k=1}^{N_{b}}|\langle \hat{\boldsymbol{u}}^{k}({\it\omega}),\boldsymbol{y}_{m}^{{\it\omega}}\rangle |.\end{eqnarray}$$
Figure 4. Comparison between the amplitude of the normalised streamwise velocity field of (a) the spectral mode computed by POD filtering, (c) the optimal response, (e) the spectral mode computed by r.m.s., for
${\it\omega}=2.1$
and two-dimensional perturbations
${\it\beta}=0$
. The three vertical dashed lines represent the locations where profiles were extracted. Figures (b), (d) and (f) compare the POD-filtered profile (thick continuous line), the optimal response profile (red dashed line) and the spectral r.m.s. profile (thin continuous line) respectively for
$x=2$
,
$x=3.5$
and
$x=5$
.
Figure 4 compares the amplitude of the streamwise component of
$\hat{\boldsymbol{u}}$
computed with method (a) (figure 4
e) and method (b) (figure 4
a), with the optimal response
${\bf\psi}_{1}$
(figure 4
c), for
${\it\omega}=2.1$
(the frequency associated with the Kelvin–Helmholtz mechanism, see the peak in figure 2
a). To ease the comparison, each field has been normalised such that its maximal amplitude is 1. Profiles extracted from these fields at three streamwise locations (
$x=2$
,
$x=3.5$
,
$x=5$
) are respectively compared in figure 4(b,d,f). These figures show that for
${\it\omega}=2.1$
, the agreement between the optimal response and the spectral mode of the flow is significantly improved when POD filtering is used. The same conclusions were obtained for all the tested frequencies of the present study. This result is reminiscent of the work of Gudmundsson & Colonius (Reference Gudmundsson and Colonius2011), who found a strongly improved agreement between the experimental spectra and the local stability results computed for a turbulent jet when these spectra were obtained by POD filtering instead of the classical Welch method.
Method (b) also gives the phase structure of the spectral modes, and for all the frequencies investigated, it displays a strong similarity with the phase of the optimal response. As an example, figure 5 illustrates this statement again for
${\it\omega}=2.1$
. Note that from (4.1),
$\hat{\boldsymbol{u}}$
and
${\bf\psi}_{1}$
are approximately equal up to a complex multiplicative constant, therefore their phase fields are equal up to an additive constant. To allow comparison, the phase of the optimal response has been shifted such that it matches the one of the spectral mode at an arbitrary position, here
$\boldsymbol{x}=(3,1.5)$
. Note that despite a very good overall agreement, a careful comparison of figures 4(a,c) and 5(a,b) shows that the strongest discrepancy between the first optimal response and the spectral mode appears in low-energy regions, where the spectral mode is not smooth (for instance the left bottom corner in figure 5(a)). Increasing the simulation time should yield better agreement, but there is no guarantee of the capability of the rank 1 approximation to finely describe these low-energy regions. This issue is addressed in § 5.4.
Since the spectral mode computed with method (b) compares better with the first optimal response, the remaining discussions of the paper concerning the spectral mode all refer to the mode computed with method (b).
Figure 5. Comparison between the phase of the streamwise velocity field of (a) the spectral mode computed by POD filtering, (b) the optimal response, for
${\it\omega}=2.1$
and two-dimensional perturbations
${\it\beta}=0$
. The figures display ten equally spaced contours ranging from
$-{\rm\pi}$
to
${\rm\pi}$
. The three vertical dashed lines represent the locations where profiles were extracted. (c–e) Compare the unwrapped POD-filtered phase profile (continuous line) and the unwrapped optimal response phase profile (red dashed line) respectively for
$x=2$
,
$x=3.5$
and
$x=5$
.
5.3 Link with local stability analysis
The question of the link between a resolvent analysis and a local stability analysis has been addressed by Sipp & Marquet (Reference Sipp and Marquet2013). They showed that a flow presenting a strong convective instability, such as a backward facing step flow, displays respectively an upstream/downstream location of the optimal forcings/responses. In such a situation, the spatial growth of the first optimal response compares favourably with the growth rate computed from a local spatial stability analysis: by defining the streamwise energy density
$d_{{\bf\psi}1}(x,{\it\omega})$
of
$\boldsymbol{{\bf\psi}}_{1}({\it\omega})$
as
$$\begin{eqnarray}d_{{\bf\psi}1}(x,{\it\omega})=\int _{y_{min}}^{y_{max}}\boldsymbol{{\bf\psi}}_{1}(x,y,{\it\omega})\boldsymbol{\cdot }{\bf\psi}_{1}(x,y,{\it\omega})\,\text{d}y,\end{eqnarray}$$
then the energy growth
${\it\sigma}_{{\bf\psi}1}=(\partial d_{{\bf\psi}1}/\partial x)/2d_{{\bf\psi}1}$
at a given frequency is close to the spatial amplification rate obtained from a local spatial stability analysis. Figure 6(a) compares these two quantities in the backward facing step case for
${\it\omega}=2.1$
, and shows that they display a good agreement in the region where the optimal forcing is weak (the local spatial stability results have been obtained by following the procedure detailed in appendix A). This link between a local stability analysis and the dominant optimal response also justifies the argument that when a convective instability mechanism of a flow is strong, then the dominant optimal forcing
${\bf\phi}_{1}$
will produce a particularly energetic response (see Garnaud et al.
Reference Garnaud, Lesshafft, Schmid and Huerre2013; Sipp & Marquet Reference Sipp and Marquet2013). Therefore, the associated gain
${\it\mu}_{1}^{2}$
will be high. Thus, the configurations presenting a strong convective instability mechanism are good candidates to fulfil the dominant singular value assumption required for the results of § 3. This conjecture is supported by the literature, which shows that such flows indeed present a clear separation of singular values (Dergham et al.
Reference Dergham, Sipp and Robinet2013; Boujo & Gallaire Reference Boujo and Gallaire2015; Sartor et al.
Reference Sartor, Mettot, Bur and Sipp2015).
We will now show that not only the energy growth rates of the spatial most amplified
$k^{+}$
mode and of the dominant optimal response mode compare well, but also the shape of the modes, at least in the locally unstable region (where the local mode is amplified). In this region (
$x<3$
for
${\it\omega}=2.1$
, see figure 6
a), the shape of the most amplified local
$k^{+}$
mode at a given streamwise location is close to the leading optimal response profile at the same location (figure 6
b). However, in the region where the mode is spatially damped, this agreement may strongly deteriorate due to the fact that the spatial structure of the leading optimal response mode may no longer be dominated by one single spatial mode. In the backward facing step flow, this leads to a strong discrepancy between the first optimal response and the local eigenmode in the spatially damped regions of the mode, as illustrated in figure 6(c) for
$x=3.1$
and
${\it\omega}=2.1$
. Note that this failure of the local spatial stability analysis was observed for all the tested frequencies of the present study.
Figure 6. Comparison between dominant optimal response mode (black continuous line) and most amplified
$k^{+}$
mode obtained with local spatial stability analysis (red dashed line) for
${\it\omega}=2.1$
: (a) streamwise energy growth of the optimal response versus spatial growth rate computed from a local spatial stability analysis, (b) amplitude profile of the optimal response mode at
$x=0.75$
versus local spatial mode at the same location (streamwise velocity component), (c) amplitude profile of the optimal response mode at
$x=3.1$
versus local spatial mode at the same location (streamwise velocity component). The local stability profiles are normalised such that their maxima match those of the optimal response profile.
This problem can be overcome by solving the parabolised stability equations (see appendix B), initialised just downstream of the step (amplified region) by the most amplified local
$k^{+}$
mode. This procedure yields an accurate reconstruction of the perturbation field in both the amplified and damped region of the mode. Figure 7(a,b) illustrates the similarity between
${\bf\psi}_{1}$
and the PSE field, for
${\it\omega}=2.1$
(to ease the comparison, both fields have been normalised such that their maximum amplitude is 1, and such that their phase match at an arbitrary point, here
$(x=3,y=1.5)$
). This agreement was observed for all the investigated frequencies. This shows that under the condition that the resolvent displays a dominant singular value (or that the flow displays a strong local spatial growth rate), solving the PSE equations around the mean flow may provide a good prediction of the spatial structure of the unsteadiness. Note that, even if here it gives excellent results, PSE is not guaranteed to work for any configuration since it may damp any behaviour not consistent with a single spatial wavelength.
As a conclusion, the link between the resolvent analysis and the local spatial stability analysis shows that the findings of the present study also hold when based on quantities obtained from a local stability analysis: if a mean flow exhibits a strong convective instability at some frequency (the spatial growth rate is strictly positive in some region of the flow), then the structure of the spectral mode at that frequency should correspond to the structure obtained with a PSE analysis based on the mean flow.
Figure 7. Comparison between the real part of the streamwise velocity of (a) the dominant optimal response and (b) the PSE field, for
${\it\omega}=2.1$
(normalised in phase and amplitude). Each figure displays ten equally spaced contours ranging from
$-$
1 to 1.
5.4 Comparison of the predictive models for the local frequency spectra
The ability of model 1 and model 2 (defined in § 4) to predict the frequency spectrum of the flow field at given points is assessed by comparing the resulting estimated energy spectrum with results obtained from the unsteady simulation data.
The question of the choice of the matching points for model 2 needs to be addressed beforehand (see § 4). Since model 2 consists of predicting the frequency spectra at every point of the domain by matching, at all frequencies, the optimal response mode and the spectral mode at given spatial locations, it may lead to poor predictions if these locations are ill chosen. In § 5.2, we showed that, for a given frequency, the energetic regions were well predicted by the dominant optimal response. But the prediction was not so accurate in low-energy regions. This probably stems from the inability of the rank 1 approximation to capture the low-energy regions of the spectral mode. Hence, at a given frequency, low-energy regions should be avoided for the choice of the matching point. This issue is illustrated in figure 8(a): for an arbitrary frequency
${\it\omega}=5$
, the amplitude profiles at
$x=1$
of
$\hat{\boldsymbol{u}}$
and
${\bf\psi}_{1}$
are rescaled such that their amplitudes match at
$\boldsymbol{x}_{0}=(1,1.4)$
, which is a point displaying a low level of energy at this particular frequency. This results in two quite different profiles. However, when the rescaling is based on a point whose energy is high, for instance
$\boldsymbol{x}_{1}=(1,0.9)$
, the overall agreement is strongly improved (figure 8
b). A robust way to compute the function
${\it\Lambda}$
of model 2 is to consider several points, such that at least one point is located in an energetic area of the flow at all the frequencies of interest, and then solving the overdetermined set of
$n$
equations
$\{{\it\Lambda}({\it\omega})=\hat{u}_{i}(\boldsymbol{x}_{k},{\it\omega})/{\it\psi}_{1,i}(\boldsymbol{x}_{k},{\it\omega}),1\leqslant k\leqslant n\}$
by least squares, where
$n$
is the number of points considered. Figure 8(c) shows the result of least-squares fitting based on the two points
$\boldsymbol{x}_{0}=(1,1.4)$
and
$\boldsymbol{x}_{1}=(1,0.9)$
. The resulting profile is close to what is obtained when
$\boldsymbol{x}_{0}$
is not considered. In brief, the least-squares technique chooses automatically the best matching point at each frequency (several matching points are generally mandatory due to the fact that the energetic regions of the spectral modes are not the same at all frequencies).
Figure 8. Illustration of the issue of the choice of the matching point for the determination of
${\it\Lambda}$
. The three figures display profiles at
$x=1$
extracted from
$\hat{\boldsymbol{u}}$
(continuous line) and the optimal response (red dashed line), for
${\it\omega}=5$
: (a) and (b) respectively compare the resulting profiles when the optimal response is scaled to match the amplitude of the low-energy point
$(x=1,y=1.4)$
and the higher-energy point
$(x=1,y=0.9)$
; (c) shows the result after a scaling that minimizes the square of the error at those two points.
For the present study, we used two points for the determination of
${\it\Lambda}$
, located at
$\boldsymbol{x}_{0}=(4,1.5)$
and
$\boldsymbol{x}_{1}=(7,0.1)$
. These points have been chosen in order to correctly capture both the high-frequency behaviour near the edge and the low-frequency unsteadiness in the reattachment region. Figure 9 presents the results for the prediction of the spectra of the streamwise component
$\hat{u} _{1}$
of the velocity, at five points
$\boldsymbol{x}_{a}=(3,1.5)$
,
$\boldsymbol{x}_{b}=(5,3)$
,
$\boldsymbol{x}_{c}=(6,1.5)$
,
$\boldsymbol{x}_{d}=(8,0.1)$
and
$\boldsymbol{x}_{e}=(9,0.1)$
. Note that since model 1 is unable to predict the amplitude of the spectra, it has been rescaled such that its maximum matches the one predicted by model 2. Model 1 gives a reasonable order of magnitude for the energetic frequencies of the flow, and in particular it is able to qualitatively discriminate if a point displays either a high- or a low-frequency content. On the other hand, model 2 led to a much more accurate prediction of those frequencies. The overall shape of the spectra is well predicted, as well as the amplitude that remarkably agrees with the simulation data for all five tested points.
Figure 9. Comparison of the estimate of the streamwise velocity spectra, from model 1 (dotted line) and model 2 (red dashed line), with the simulation (continuous line), at five points: (a)
$\boldsymbol{x}_{1}=(3,1.5)$
, (b)
$\boldsymbol{x}_{2}=(5,3)$
, (c)
$\boldsymbol{x}_{3}=(6,1.5)$
, (d)
$\boldsymbol{x}_{4}=(8,0.1)$
, (e)
$\boldsymbol{x}_{5}=(9,0.1)$
.
6 Concluding remarks
The resolvent operator around the mean flow appears naturally in the full nonlinear equations that govern the perturbation field around the mean flow. If the mean flow exhibits a strong convective instability, the first singular value of this operator is expected to dominate all the others. Consequently, by assuming that the turbulent forcing does not display any preferential direction toward one of the suboptimal forcings, the resolvent can be approximated by using this singular value and the associated optimal response and forcing modes. This approximation yields a mathematical link between the actual dynamics of the flow (the spectral mode) and the dominant optimal response mode. In the case of weakly non-parallel flows, we have shown that the dominant singular value assumption generally holds at some frequency if a strong convective instability mechanism affects the mean flow at that frequency. In such a case, the spatial structure of the spectral mode should be very close to the spatial structure obtained by a PSE analysis around the mean flow.
This article also focuses on the development of methods to predict the frequency spectrum at every point of a flow. It presents a model (model 2) based on partial information of the flow, in order to reconstruct the complete fluctuation field at any point. In experimental applications, it is often easy to accurately obtain the mean flow as well as an accurate spectrum at a few points in the flow, for example wall measurements such as pressure or skin friction. In the context of compressible flows, the experimental jet set-up from Gudmundsson & Colonius (Reference Gudmundsson and Colonius2011) uses an array of microphones measuring the pressure fluctuations just outside the jet shear layer. The approach that is developed in the present study could give a reconstruction of the local spectra at any point from this data and from a measurement of the mean flow. The extent of the region in which mean flow measurements are required also depends on the method chosen to compute the dominant optimal response modes. This region is smaller when using the PSE, since the computation is restricted to the spatial support of those modes. In contrast, a much larger domain is required in the streamwise direction for a global stability computation in order to adapt to the upstream and downstream boundary conditions. Finally, time-resolved data for frequency spectra are only needed at a few points; the mean flow may be obtained by less sophisticated means, such as 5 hole probes or traditional (non-time-resolved) particle image velocimetry measurements. The methodology presented in this article may be well suited for the exploitation of experimental results.
Appendix A. Local spatial stability analysis
The 2-D incompressible equations of the local stability analysis may be derived by considering a parallel base flow whose velocity is written in Cartesian coordinates
$\boldsymbol{U}=(U(y),0)$
. Considering small fluctuations
$(u^{\prime },v^{\prime },p^{\prime })$
around the base flow, where
$u^{\prime }$
,
$v^{\prime }$
and
$p^{\prime }$
are respectively the fluctuating streamwise velocity, the fluctuating cross-stream velocity and the fluctuating pressure around the base flow, gives the following linearised equations:
$$\begin{eqnarray}\displaystyle & u_{x}^{\prime }+v_{y}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & u_{t}^{\prime }+Uu_{x}^{\prime }+U_{y}v^{\prime }=-p_{x}^{\prime }+{\it\nu}(u_{xx}^{\prime }+u_{yy}^{\prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & v_{t}^{\prime }+Uv_{x}^{\prime }=-p_{y}^{\prime }+{\it\nu}(v_{xx}^{\prime }+v_{yy}^{\prime }), & \displaystyle\end{eqnarray}$$
${\it\nu}$
is the inverse of the Reynolds number. Assuming that every fluctuating quantities is of the form
$X^{\prime }(x,y,t)=\hat{X}(y)\exp ({\it\alpha}x-\text{i}{\it\omega}t)$
yields the final set of equations: 
$$\begin{eqnarray}\displaystyle & {\it\alpha}\hat{u} +D\hat{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (-{\it\nu}D^{2}+U{\it\alpha}-{\it\nu}{\it\alpha}^{2})\hat{u} +DU\hat{v}+{\it\alpha}\hat{p}=\text{i}{\it\omega}\hat{u} , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (-{\it\nu}D^{2}+U{\it\alpha}-{\it\nu}{\it\alpha}^{2})\hat{v}+D\hat{p}=\text{i}{\it\omega}\hat{v}, & \displaystyle\end{eqnarray}$$
$D$
designates the differentiation operator with respect to
$y$
. By setting
$\hat{\boldsymbol{q}}=(\hat{u} ,\hat{v},\hat{p})^{\text{T}}$
and
$A_{1}$
,
$A_{2}$
,
$A_{3}$
,
$A_{4}$
as 
$$\begin{eqnarray}A_{1}=\left(\begin{array}{@{}ccc@{}}0 & 0 & 0\\ -{\it\nu} & 0 & 0\\ 0 & -{\it\nu} & 0\end{array}\right),\quad A_{2}=\left(\begin{array}{@{}ccc@{}}1 & 0 & 0\\ U & 0 & {\it\alpha}\\ 0 & U & 0\end{array}\right),\end{eqnarray}$$
$$\begin{eqnarray}A_{3}=\left(\begin{array}{@{}ccc@{}}0 & D & 0\\ -{\it\nu}D^{2} & DU & 0\\ 0 & -{\it\nu}D^{2} & 0\end{array}\right),\quad A_{4}=\left(\begin{array}{@{}ccc@{}}0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 1 & 0\end{array}\right),\end{eqnarray}$$
equations (A 4); (A 5); (A 6) may be recast as
$$\begin{eqnarray}{\it\alpha}^{2}A_{1}\hat{\boldsymbol{q}}+{\it\alpha}A_{2}\hat{\boldsymbol{q}}+A_{3}\hat{\boldsymbol{q}}=\text{i}{\it\omega}A_{4}\hat{\boldsymbol{q}}.\end{eqnarray}$$
A linear spatial stability analysis consists of computing
${\it\alpha}$
and
$\hat{\boldsymbol{q}}$
for a given fixed real value of
${\it\omega}$
. This can be achieved by considering the augmented eigenvalue problem for
${\it\alpha}$
:
$$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}0 & -A_{3}+\text{i}{\it\omega}A_{4}\\ I & -A_{2}\end{array}\right)\left(\begin{array}{@{}c@{}}X\\ Y\end{array}\right)={\it\alpha}\left(\begin{array}{@{}cc@{}}I & 0\\ 0 & A_{1}\end{array}\right)\left(\begin{array}{@{}c@{}}X\\ Y\end{array}\right),\end{eqnarray}$$
with
$I$
the identity operator. The pair
$({\it\alpha},\hat{\boldsymbol{q}}=Y)$
formed from an eigenvalue and eigenvector of (A 10) is also a solution of (A 9).
For the backward facing step case, the boundary conditions were as follows:
$\hat{u} =\hat{v}=0$
at the wall location, and
$\hat{u} =\hat{v}=\hat{p}=0$
at the other end. For the present study, the equations have been discretised using second-order centred finite differences for the derivatives, then the eigenvalue problem (A 10) has been solved at each streamwise location using the ARPACK library (Lehoucq et al.
Reference Lehoucq, Sorensen and Yang1998), with a shift-and-invert strategy. For each streamwise location, the computation of the eigenvalue was initialised by the value of
${\it\alpha}$
computed just upstream from the current location. The first computation was initialised using the energy growth computed from the SVD (see § 5.3).
Appendix B. PSE
In the following, we use notation similar to that in appendix A. The PSE accounts for the effect of slight base flow spreading (Herbert Reference Herbert1997). Therefore, the base flow is now taken to be of the form
$\boldsymbol{U}=(U(x,y),V(x,y))$
and the perturbation
$(u^{\prime },v^{\prime },p^{\prime })$
is separated into an amplitude function and an exponential function. For instance, the streamwise component is written as:
$$\begin{eqnarray}u^{\prime }=\hat{u} (x,y)\exp \left(\int _{x_{0}}^{x}{\it\alpha}({\it\xi})\,\text{d}{\it\xi}-\text{i}{\it\omega}t\right).\end{eqnarray}$$
Substituting this form into (A 1); (A 2); (A 3) yields the PSE equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{u} _{x}=-{\it\alpha}\hat{u} -D\hat{v}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{v}_{x}=\frac{1}{U}(\text{i}{\it\omega}-{\it\alpha}U-VD-V_{y}+{\it\nu}(D^{2}+{\it\alpha}^{2}))\hat{v}-\frac{D}{U}\hat{p}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{p}_{x}=(\text{i}{\it\omega}-VD-U_{x}+{\it\nu}(D^{2}+{\it\alpha}^{2}))\hat{u} +(UD-U_{y})\hat{v}-{\it\alpha}\hat{p}. & \displaystyle\end{eqnarray}$$
$\text{d}/\text{d}x,V\sim O({\it\nu})$
, and then neglecting all the terms of order
${\sim}{\it\nu}^{2}$
and higher. In the present study, the discretisation of these equations is based on second-order centred finite differences. The previous equations may be recast in the compact matrix form: 
$$\begin{eqnarray}\hat{\boldsymbol{q}}_{x}=\mathscr{L}\hat{\boldsymbol{q}},\end{eqnarray}$$
with
$\hat{\boldsymbol{q}}=(\hat{u} ,\hat{v},\hat{p})^{\text{T}}$
, and
$\mathscr{L}$
the operator corresponding to the linear equations (B 2); (B 3); (B 4). Since both the amplitude functions and the streamwise wavenumber depend on
$x$
, an additional equation is needed to remove all ambiguity. We assume that the amplitude functions verify the auxiliary condition:
$$\begin{eqnarray}\int _{y_{min}}^{y_{max}}(\hat{u} ^{\ast }\hat{u} _{x}+\hat{v}^{\ast }\hat{v}_{x}+\hat{p}^{\ast }\hat{p}_{x})\,\text{d}y=0.\end{eqnarray}$$
The boundary conditions are the same as those mentioned in appendix A. We then solve the equations starting from an upstream location initialised with a solution computed from a local spatial stability analysis, and then marching downstream by solving the first-order approximation:
$$\begin{eqnarray}\frac{\hat{\boldsymbol{q}}(x+{\it\delta}x)-\hat{\boldsymbol{q}}(x)}{{\it\delta}x}=\mathscr{L}\hat{\boldsymbol{q}}(x+{\it\delta}x).\end{eqnarray}$$
The auxiliary condition is ensured at each location by iteratively solving (B 7) then updating
${\it\alpha}$
following the algorithm described for example in Gudmundsson & Colonius (Reference Gudmundsson and Colonius2011). Note that for the present study, we used the stabilising procedure described in Andersson, Henningson & Hanifi (Reference Andersson, Henningson and Hanifi1998), which consists of rewriting (B 5) as:
$$\begin{eqnarray}\hat{\boldsymbol{q}}_{x}=\mathscr{L}\hat{\boldsymbol{q}}+s\mathscr{L}\hat{\boldsymbol{q}}_{x},\end{eqnarray}$$
where
$s$
is a positive real number based on the step
${\it\delta}x$
of the marching method and on the local value of
${\it\alpha}$
(see Andersson et al.
Reference Andersson, Henningson and Hanifi1998).
