2.1. Numerical set-up

Each simulation was set according to Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2020), based on the work of Brandt et al. (Reference Brandt, Schlatter and Henningson2004), and consists of a $231 \times 121 \times 108 (x \times y \times z)$ Cartesian grid constructed with Chebyshev nodes in the $y$ direction, perpendicular to the wall, and homogeneously spaced nodes in the streamwise and spanwise directions, $x$ and $z$. The boundary layer is started with a finite thickness. All variables are non-dimensionalised by the reference length $\delta ^*_0$, the boundary-layer displacement thickness at the intake, and a time scale $t = \delta ^*_0/U_\infty$, where $U_\infty$ is the free-stream velocity. The numerical domain is a box of size $x \in [0,1000]$, $y \in [0,60]$ and $z \in [-25,25]$. Both the $x$ and $z$ directions are periodic and decomposed in Fourier modes, while the $y$ direction uses a Chebyshev polynomial basis. Periodicity in the streamwise direction is achieved by the introduction of a fringe region contained in the range $x \in [910,1000]$.

At the intake, $Re_{\delta _0^*}= U_\infty \delta ^*_0/\nu = 300$, with $\nu$ being the kinematic viscosity of the fluid. At this Reynolds number, the boundary layer is linearly stable and, thus, TS waves are not expected to be significant over the relatively short extent of the domain. Instead, streamwise elongated (streaky) structures are observed in the simulations at the highest FST levels investigated, as shown in figure 1.

Figure 1. Snapshot of a simulation with $Tu=3.5\,\%$ incoming FST level, where streaky structures can be identified. Slice at $y = 0.8$, inside the boundary layer.

The no-slip condition

(2.1)\begin{equation} \boldsymbol{u^\prime}(x,0,z,t) = 0, \end{equation}

is imposed on the wall and the Neumann condition

(2.2)\begin{equation} \frac{\partial}{\partial y} \boldsymbol{u^\prime}(x,60,z,t) = 0, \end{equation}

is applied on the upper boundary, with $\boldsymbol {u^\prime }(x,y,z,t)$ representing velocity fluctuations with respect to the two-dimensional Blasius base flow, $\boldsymbol {U}_{BL}(x,y)$ (figure 2). For all performed simulations, the physical domain ends before the development of turbulent spots in the boundary layer, i.e. before the transition to the turbulent regime. The use of a short spatial domain reduces the computational cost of the present study, which involves detailed post-processing of several numerical simulations. Moreover, by restricting the domain to the initial development of streaks we can focus on the receptivity stage, before the actual transition to turbulence that would occur at larger values of $x$.

Figure 2. Diagram of the boundary-layer set-up, showing boundary conditions. The $x$ and $z$ directions are periodic and $\delta _0$ is the initial boundary-layer thickness. Here, $\boldsymbol {U_{BL}}(x,y)$ is the Blasius base flow and $\boldsymbol {u^\prime }(x,y,z)$ are velocity fluctuations. Legend: (dotted line) forcing modes from the continuous branch of the OSS operator.

Concerning the time evolution, linear terms of the Navier–Stokes (NS) equations are implicitly marched with a second-order Crank–Nicolson scheme, while an explicit third-order, four-stage, Runge–Kutta scheme is applied over nonlinear terms. For each simulation, we compute a total of 2000 snapshots, taken in time steps of $\Delta t=10$, of fully developed, statistically stationary, flow.

2.2. Fringe region forcing

Some assumptions are made to synthesise valid inflow conditions at $x=0$ and circumvent the need to compute a turbulent field upstream of the flat plate or the flow around a leading edge. Isotropic and homogeneous FST is introduced in the simulations by forcing several modes in the continuous branch of the linearised OSS operator within the fringe region, as illustrated in figure 2.

The FST generation procedure is referred to in Schlatter (Reference Schlatter2001) and Brandt et al. (Reference Brandt, Schlatter and Henningson2004), based on the methods presented in Grosch & Salwen (Reference Grosch and Salwen1978) and Jacobs & Durbin (Reference Jacobs and Durbin2001). Considering the linearised NS (LNS) momentum equations in perturbation form around a base flow and nonlinear fluctuation terms gathered into the function $f(\boldsymbol {u^\prime })$, we force a desired velocity vector $\boldsymbol {\zeta }(x,y,z,t)$ inside the fringe following the formulation

(2.3)\begin{equation} \left. \begin{array}{c@{}} \dfrac{\partial \boldsymbol{u^\prime}}{\partial t} = LNS (\boldsymbol{u^\prime},\boldsymbol{U_{BL}}) + {f(\boldsymbol{u^\prime})} + \sigma(x)(\boldsymbol{\zeta}-\boldsymbol{u^\prime}),\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u^\prime} = 0, \end{array}\right\} \end{equation}

where $\sigma (x)$ is a gain function, which is positive inside the fringe region and null everywhere else (figure 3). The term $\sigma (x)(\boldsymbol {\zeta }-\boldsymbol {u^\prime })$ is thus responsible for smoothly changing the fluctuation field entering the left side of the fringe region towards the desired reference forcing vector $\boldsymbol {\zeta }$ introduced in the fringe.

Figure 3. Fringe gain function, the same as that described in Chevalier et al. (Reference Chevalier, Lundbladh and Henningson2007). The maximum gain inside the fringe is set to 0.8.

Isotropic homogeneous turbulence can be represented as a sum of Fourier modes with random amplitude (Rogallo Reference Rogallo1981). In the boundary-layer case, however, this approach is not capable of modelling the presence of the wall, as the $y$ direction is non-homogeneous. For this application, a basis composed of modes in the continuous spectrum of the OSS operator is assumed to be a reasonable choice to satisfy all the necessary boundary conditions. These modes tend to Fourier modes far from the wall and decay to zero near it, generating a perturbation field mainly localised outside the boundary layer, as shown in Appendix A. On the other hand, modes of the discrete spectrum are only significant inside the boundary layer and decay exponentially farther from the wall, not being suitable in this application (Grosch & Salwen Reference Grosch and Salwen1978).

By computing eigenfunctions in the continuous branch of the OSS spectrum, $\boldsymbol {u^\prime _{OSS}}$, normalised to unit energy, we can write the expansion for an arbitrary perturbation vector

(2.4)\begin{equation} \boldsymbol{\zeta}(x,y,z,t) = {\rm Re}\left\{\sum_{\omega} \sum_{\gamma} \sum_{\beta} \varPhi(\omega,\gamma,\beta) \boldsymbol{u^\prime_{OSS}}(\omega,y,\beta) \, {\rm e}^{{\rm i}({\rm Re}\{\alpha(\omega,\gamma,\beta)\} x + \beta z - \omega t )}\right\}, \end{equation}

where $\omega$, $\gamma$ and $\beta$ are real parameters and $\alpha (\omega,\gamma,\beta )$ is the complex eigenvalue of $\boldsymbol {u^\prime _{OSS}}$ computed via spatial stability (Jacobs & Durbin Reference Jacobs and Durbin2001). Here, $\alpha,\gamma,\beta$ are respectively the wavenumbers in the $x,y,z$ directions and $\omega$ the frequency. The factor $\varPhi$ is the energy scaling applied to match the von Kármán spectrum, discussed in the following paragraphs. Only the real part of $\alpha$ is taken inside the exponent to maintain the forcing fluctuation at a fixed magnitude throughout the fringe zone's streamwise extension, ignoring in practice the effects of viscous attenuation (Brandt et al. Reference Brandt, Schlatter and Henningson2004).

We consider wavenumbers $\kappa = \sqrt {\textrm {Re}\{\alpha \}^2+\gamma ^2+\beta ^2}$ equally spaced within the range limited by the numerical resolution of the simulations, $\kappa \in [\kappa _l,\kappa _u]$. In general, $\kappa _l$ is a function of the domain size, while $\kappa _u$ is bounded by the resolution of the grid. For simplification, we replace $\omega = \alpha U_\infty$, considering that modes of the continuous spectrum have phase speed equal to $U_\infty$, to define a tridimensional space of parameters $(\omega,\gamma,\beta )$ for which a given value $\kappa$ is represented by a spherical shell (Brandt et al. Reference Brandt, Schlatter and Henningson2004). We select $N_s$ shells, within which we include $N_\kappa$ combinations of the $(\omega,\gamma,\beta )$ parameters of constant $\kappa$, filling the surface with equally spaced points (Schlatter Reference Schlatter2001). The value $\gamma =0$ is avoided. A random rotation is applied to each shell at every time step to further improve isotropy. In this work, we adopt the values $\kappa _l=0.23$, $\kappa _u=3.0$, $N_s = 20$ and $N_\kappa = 10$, in a total of $N_s N_\kappa = 200$ eigenfunctions, the same as in Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2020).

Once the suitable modes are chosen, the energy scale needs to be applied. Considering the von Kármán spectrum for isotropic homogeneous turbulence and following the three-dimensional spectrum construction in Tennekes & Lumley (Reference Tennekes and Lumley1972), we have the formula for turbulent energy as a function of wavenumber

(2.5)\begin{equation} E(\kappa) = \frac{2}{3} \frac{a (\kappa L)^4}{(b + (\kappa L)^2 )^{{17}/{6}}} L \cdot Tu^2, \quad L = \frac{1.8}{\kappa_{max}}, \end{equation}

where $a = 1.606$, $b = 1.35$ and $Tu$ is the turbulence intensity level defined as

(2.6)\begin{equation} Tu = \sqrt{\frac{(u'^2_{rms}+v'^2_{rms}+w'^2_{rms})}{3}}. \end{equation}

In this equation, the integral length scale $L=7.5 \delta _0^*$ is set to the same value considered in Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2020), yielding a wavenumber of maximum energy, $\kappa _{max}$, near the minimum allowed value of $\kappa _l$. According to the results shown in Brandt et al. (Reference Brandt, Schlatter and Henningson2004), the increase of the turbulence integral length reduces the turbulence intensity decay at the free stream and promotes transition in positions further upstream. Therefore, this choice of integral length scale consists of a worst case scenario, which allows a shorter domain size in the streamwise direction.

Concerning the energy scaling, it is demonstrated in Schlatter (Reference Schlatter2001) that the factor $\varPhi$ in (2.4) can be then expressed as

(2.7)\begin{equation} \varPhi(\kappa) = \sqrt{\frac{E(\kappa) \Delta \kappa}{N_s}}, \end{equation}

where $\Delta \kappa$ is the difference between consecutive values of $\kappa$.

Finally, the amplitudes of OSS modes in the continuous branch of the spectrum must be addressed at the top boundary of the domain. To prevent numerical instabilities, we multiply the eigenfunctions by a smooth step function $S(\kern0.7pt y)$ (Brandt et al. Reference Brandt, Schlatter and Henningson2004) to dampen forcing perturbations above the position $y_d = 0.8 y_{max}$.

A more detailed discussion concerning the properties of the inflow perturbations generated using OSS modes in the continuous branch is presented in Appendix A.

3.1. Input–output formulation

To apply the resolvent analysis framework over NS equations, we separate the velocity field into a two-dimensional, time-invariant, laminar solution (Jovanović & Bamieh Reference Jovanović and Bamieh2005) or ensemble average flow (McKeon & Sharma Reference McKeon and Sharma2010)

(3.1)\begin{equation} \boldsymbol{U} = [U(x,y),V(x,y),0]^T, \end{equation}

and fluctuations

(3.2)\begin{equation} \boldsymbol{u^\prime} = [u^\prime(x,y,z,t),v^\prime(x,y,z,t),w^\prime(x,y,z,t)]^T, \end{equation}

in order to write the linearised equations around $\boldsymbol {U}$, as described in (2.3). Using tensor formulation, the system can be written as

(3.3)\begin{equation} \left. \begin{array}{c@{}} \displaystyle \dfrac{\partial u^\prime_i}{\partial t} + U_j \dfrac{\partial u^\prime_i}{\partial x_j} + u^\prime_j \dfrac{\partial U_i}{\partial x_j} ={-} \dfrac{\partial p^\prime}{\partial x_i} + \dfrac{1}{Re} \dfrac{\partial^2 u^\prime_i}{\partial x_j \partial x_j} + f_i + \sigma(\zeta_i - u^\prime_i), \\ \displaystyle \dfrac{\partial u^\prime_j}{\partial x_j} = 0, \end{array}\right\} \end{equation}

with nonlinear terms grouped in $f_i = - u^\prime _j ({\partial u^\prime _i}/{\partial x_j})$, considered in the resolvent framework as a forcing that drives the linear dynamics. The term $\zeta _i$ is the forcing vector defined in (2.3), which guarantees that inlet conditions in the model statistically match those observed in the boundary-layer simulations. The function $\sigma (x)$ is the same as presented in figure 3.

Next, we apply the normal mode ansatz $\boldsymbol {u^\prime } = \hat {\boldsymbol {u}}(x,y) \, \textrm {e}^{\textrm {i} (\beta z - \omega t)}$ over the velocity, pressure and forcing fields to expand (3.3) as

(3.4) \begin{align} \left. \begin{array}{c@{}} \displaystyle -{\rm i} \omega \hat{u} + U \dfrac{\partial \hat{u}}{\partial x} + V \dfrac{\partial \hat{u}}{\partial y} + \hat{u} \dfrac{\partial U}{\partial x} + \hat{v} \dfrac{\partial U}{\partial y} + \dfrac{\partial \hat{p}}{\partial x} - \dfrac{1}{Re}\left(\dfrac{\partial^2}{\partial x^2} -\beta^2 + \dfrac{\partial^2}{\partial y^2} \right) \hat{u} = \hat{f}_x + \sigma (\hat{\zeta}_x-\hat{u}), \\ \displaystyle -{\rm i} \omega \hat{v} + U \dfrac{\partial \hat{v}}{\partial x} + V \dfrac{\partial \hat{v}}{\partial y} + \hat{u} \dfrac{\partial V}{\partial x} + \hat{v} \dfrac{\partial V}{\partial y} + \dfrac{\partial \hat{p}}{\partial y} - \dfrac{1}{Re}\left(\dfrac{\partial^2}{\partial x^2} -\beta^2 + \dfrac{\partial^2}{\partial y^2} \right) \hat{v} = \hat{f}_y + \sigma (\hat{\zeta}_y-\hat{v}), \\ \displaystyle -{\rm i} \omega \hat{w} + U \dfrac{\partial \hat{w}}{\partial x} + V \dfrac{\partial \hat{w}}{\partial y} + {\rm i} \beta \hat{p} + \dfrac{1}{Re} \left(\dfrac{\partial^2}{\partial x^2} -\beta^2 + \dfrac{\partial^2}{\partial y^2} \right) \hat{w} = \hat{f}_z + \sigma (\hat{\zeta}_z-\hat{w}), \\ \displaystyle \dfrac{\partial \hat{u}}{\partial x} + \dfrac{\partial \hat{v}}{\partial y} + {\rm i} \beta \hat{w} = 0, \end{array} \right\} \end{align}

where the nonlinear term $\widehat {f_i}$ can be written as a convolution of the Fourier transform of velocity components

(3.5)\begin{equation} \widehat{f_i}(\beta,\omega) ={-} \hat{u}_j \ast \frac{\partial \hat{u}_i}{\partial x_j} ={-} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \hat{u}_j(\beta_0,\omega_0) \frac{\partial }{\partial x_j} \hat{u}_i(\beta-\beta_0,\omega-\omega_0) \, {\rm d} \beta_0 \, {\rm d} \omega_0. \end{equation}

This implies that $\widehat {f_i}$ are the only terms responsible for energy transfers between different wavenumbers ($\beta,\beta _0,\beta -\beta _0$) and frequencies ($\omega,\omega _0,\omega -\omega _0$), in triads related to the turbulent energy cascade (Cheung & Zaki Reference Cheung and Zaki2014; Moffatt Reference Moffatt2014).

In practice, the fringe perturbation vector in Fourier space, $\hat {\boldsymbol {\zeta }}$, is approximated by the velocity fluctuation field computed from the simulations, denoted as $\hat {\boldsymbol {u}}_{\boldsymbol r}$, which is substituted in (3.4) for all three spatial components. Next, this equation is discretised reproducing the same grid of the LES and equivalent boundary conditions to write the system in state-space form

(3.6)\begin{equation} \left. \begin{array}{c@{}} (\boldsymbol{\varOmega}+{\boldsymbol{\mathsf{L}}}) \hat{\boldsymbol{q}} = {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat{\boldsymbol{u}}_{\boldsymbol r} + {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{\boldsymbol{f}}, \\ \hat{\boldsymbol{y}} = {\boldsymbol{\mathsf{H}}} \hat{\boldsymbol{q}}, \end{array} \right\}, \end{equation}

and obtain

(3.7)\begin{equation} {\boldsymbol{\mathsf{R}}} = {\boldsymbol{\mathsf{H}}}(\boldsymbol{\varOmega}+{\boldsymbol{\mathsf{L}}})^{{-}1} \implies \hat{\boldsymbol{y}} = {\boldsymbol{\mathsf{R}}} ({\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat{\boldsymbol{u}}_{\boldsymbol r} + {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{\boldsymbol{f}}), \end{equation}

where ${\boldsymbol{\mathsf{R}}}$ is the resolvent operator and $\hat {\boldsymbol {q}} = [\hat {u},\hat {v},\hat {w},\hat {p}]^T$, $\hat {\boldsymbol {y}} = [\hat {u},\hat {v},\hat {w}]^T$, $\hat {\boldsymbol {u}}_{\boldsymbol r} = [\hat {u}_r,\hat {v}_r,\hat {w}_r]^T$, $\hat{\boldsymbol{f}} = [\,\hat{f}_x,\hat{f}_y,\hat{f}_z]^T$ are vectors composed of row-wise stacked components. Operators $\boldsymbol {\varOmega }$, ${\boldsymbol{\mathsf{L}}} {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}}$, ${\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}}$, ${\boldsymbol{\mathsf{H}}}$ and the boundary conditions are defined in Appendix B. The operator ${\boldsymbol{\mathsf{H}}}$ simply removes $\hat {p}$ from the output while the operator ${\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}}$ restricts the application of the respective input to the region displayed in figure 4. It should be noted that the inclusion of the pressure, $p$, in the state $\hat {\boldsymbol {q}}$ removes the need to explicitly project the nonlinear forcing, $\hat{\boldsymbol{f}}$, into a solenoidal space since the incompressible LNS system will redirect any non-solenoidal component in the nonlinear forcing to the pressure field, as described in Rosenberg & McKeon (Reference Rosenberg and McKeon2019). This formulation allows for the separation of contributions from external forcing, $\hat {\boldsymbol {y}}_{\boldsymbol {L}}$, and nonlinear forcing, $\hat {\boldsymbol {y}}_{\boldsymbol N}$, as

(3.8)\begin{equation} \left. \begin{array}{c@{}} \hat{\boldsymbol{y}}_{\boldsymbol{L}} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat{\boldsymbol{u}}_{\boldsymbol r}, \\ \hat{\boldsymbol{y}}_{\boldsymbol N} = {\boldsymbol{\mathsf{R}}} {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \kern0.06em \hat{\boldsymbol{f}}, \end{array} \right\} \end{equation}

while still considering a single resolvent operator, such that $\hat {\boldsymbol {y}} = \hat {\boldsymbol {y}}_{\boldsymbol {L}} + \hat {\boldsymbol {y}}_{\boldsymbol {N}} \approx \hat {\boldsymbol {u}}_{\boldsymbol r}$. Even though the full response, $\hat {\boldsymbol {y}}$, is a superposition of linear and nonlinear components, it is not the case that $\hat {\boldsymbol {y}}_{\boldsymbol {L}}$ and $\hat {\boldsymbol {y}}_{\boldsymbol {N}}$ evolve in a dynamically independent way, since $\hat{\boldsymbol{f}}$ is a function of the field fluctuations and needs to be computed beforehand from NS simulations in the context of the resolvent framework. The component ${\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat {\boldsymbol {u}}_{\boldsymbol r}$ (linear input) accounts for the external flow perturbations coming through the domain upstream boundary and acts only in the fringe region, within a given pair $(\beta,\omega )$. On the other hand, ${\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{\boldsymbol{f}}$ (nonlinear input) acts everywhere and accounts for the energy transfers between different wavenumbers and frequencies, due to the convolutional nature $\hat{\boldsymbol{f}}$, as described by (3.5).

Figure 4. Diagram of the geometric distribution of input terms. While the nonlinear term acts everywhere, the linear term is only present inside the fringe region. Legend: (light grey) ${\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}}\hat{\boldsymbol{f}}$; (grey) ${\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat {\boldsymbol {u}}_{\boldsymbol r} + {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}}\hat{\boldsymbol{f}}$.

3.2. Spectral estimation

Both $\hat {\boldsymbol {u}}_{\boldsymbol r}$ and $\hat{\boldsymbol{f}}$ are computed directly from velocity fluctuations $\boldsymbol {u^\prime _r}$ from the simulation. Given the velocity fluctuation field $\boldsymbol {u^\prime }(x,y,z,t)$ at each snapshot, we compute nonlinear terms $\boldsymbol {f}(x,y,z,t) = -(\boldsymbol {u^\prime _r} \boldsymbol {\cdot } \boldsymbol {\nabla } ) \boldsymbol {u^\prime _r}$. Next, we apply the fast Fourier transform (FFT) in the periodic direction, $z$, to obtain $\bar {\boldsymbol {u}}_{\boldsymbol r}(x,y,\beta,t)$ and $\bar {\boldsymbol {f}}(x,y,\beta,t)$. These are organised in data matrices

(3.9a,b) \begin{equation} \bar{{\boldsymbol{\mathsf{U}}}}_{{\boldsymbol{\mathsf{r}}}} = \left[ \begin{array}{@{}cccc@{}} \mid & \mid & & \mid \\ \bar{\boldsymbol{u}}_{\boldsymbol{r}}^{\boldsymbol{(1)}} & \bar{\boldsymbol{u}}_{\boldsymbol{r}}^{\boldsymbol{(2)}} & \cdots & \bar{\boldsymbol{u}}_{\boldsymbol{r}}^{\boldsymbol{(N_t)}} \\ \mid & \mid & & \mid \end{array}\right], \quad \bar{{\boldsymbol{\mathsf{F}}}} = \left[ \begin{array}{@{}cccc@{}} \mid & \mid & & \mid \\ \bar{\boldsymbol{f}}^{\boldsymbol{(1)}} & \bar{\boldsymbol{f}}^{\boldsymbol{(2)}} & \cdots & \bar{\boldsymbol{f}}^{\boldsymbol{(N_t)}} \\ \mid & \mid & & \mid \end{array}\right], \end{equation}

each containing $N_t$ time-ordered snapshot column vectors. The spectral estimation in frequency is performed using the Welch method (Welch Reference Welch1967) via the algorithm presented in Towne et al. (Reference Towne, Schmidt and Colonius2018). This procedure returns the quantities $\hat {\boldsymbol {u}}_{\boldsymbol r}(x,y,\beta,\omega )$ and $\hat{\boldsymbol{f}}(x,y,\beta,\omega )$, which are assembled in the final spectral data matrices

(3.10a,b) \begin{equation} \hat{{\boldsymbol{\mathsf{U}}}}_{{\boldsymbol{\mathsf{r}}}} = \left[ \begin{array}{@{}cccc@{}} \mid & \mid & & \mid \\ \hat{\boldsymbol{u}}_{\boldsymbol r}^{\boldsymbol{(1)}} & \hat{\boldsymbol{u}}_{\boldsymbol r}^{\boldsymbol{(2)}} & \cdots & \hat{\boldsymbol{u}}_{\boldsymbol r}^{\boldsymbol{(N_b)}} \\ \mid & \mid & & \mid \end{array}\right], \quad \hat{{\boldsymbol{\mathsf{F}}}} = \left[ \begin{array}{@{}cccc@{}} \mid & \mid & & \mid \\ \hat{\boldsymbol{f}}^{\boldsymbol{(1)}} & \hat{\boldsymbol{f}}^{\boldsymbol{(2)}} & \cdots & \hat{\boldsymbol{f}}^{\boldsymbol{(N_b)}} \\ \mid & \mid & & \mid \end{array}\right], \end{equation}

for each pair wavenumber and frequency $(\beta,\omega )$, containing $N_b$ columns which correspond to the number of blocks used in the windowing procedure.

3.3. Spectral correction due to windowing

The presence of the windowing function in the spectral estimation adds new terms to the response of the LNS equations written in (3.3), as pointed out by Martini et al. (Reference Martini, Cavalieri, Jordan and Lesshafft2020a). Considering the operators defined in Appendix B and the matrices in (3.9a,b), we write equation (3.6) in the time domain as

(3.11)\begin{equation} \left. \begin{array}{c@{}} {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{y}}}} \dfrac{\partial \boldsymbol{y}}{\partial t} + {\boldsymbol{\mathsf{L}}} \boldsymbol{q} = {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \bar{\boldsymbol{u}}_{\boldsymbol{r}} + {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{f}}}} \bar{\boldsymbol{f}}, \\ \boldsymbol{y} = {\boldsymbol{\mathsf{H}}} \boldsymbol{q}, \end{array} \right\} , \quad {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{y}}} = {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{f}}}}. \end{equation}

Applying the Welch method for spectral estimation implies that each data block is multiplied by a windowing function $w(t)$ so that (3.11) becomes

(3.12)\begin{equation} \left. \begin{array}{c@{}} w {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{y}}}} \dfrac{\partial \boldsymbol{y}}{\partial t} + w {\boldsymbol{\mathsf{L}}} \boldsymbol{q} = w {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \bar{\boldsymbol{u}}_{\boldsymbol r} + w {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{f}}}} \bar{\boldsymbol{f}} \\ w \boldsymbol{y} = w {\boldsymbol{\mathsf{H}}} \boldsymbol{q} \end{array} \right\}. \end{equation}

The windowing function $w(t)$ commutes with all time-invariant operators, for instance,

(3.13)\begin{equation} w {\boldsymbol{\mathsf{L}}} \boldsymbol{q} = {\boldsymbol{\mathsf{L}}} (w \boldsymbol{q}), \end{equation}

but not with the time derivative, which obeys the identity

(3.14)\begin{equation} w \frac{\partial \boldsymbol{y}}{\partial t} = \frac{\partial}{\partial t} (w \boldsymbol{y}) - \frac{{\rm d} w}{{\rm d} t} \boldsymbol{y}. \end{equation}

These relations imply that (3.12) can be rewritten in the form

(3.15)\begin{equation} \left. \begin{array}{c@{}} \displaystyle {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{y}}}} \dfrac{\partial}{\partial t} (w \boldsymbol{y}) + {\boldsymbol{\mathsf{L}}} (w \boldsymbol{q}) = {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} (w \bar{\boldsymbol{u}}_{\boldsymbol r}) + {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} (w \bar{\boldsymbol{f}}) + {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{y}}}} \left(\dfrac{{\rm d} w }{{\rm d} t} \boldsymbol{y}\right),\\ w \boldsymbol{y} ={\boldsymbol{\mathsf{H}}} (w \boldsymbol{q}), \end{array} \right\} \end{equation}

and transformed into frequency space, as

(3.16)\begin{equation} \left. \begin{array}{c@{}} (\boldsymbol{\varOmega}+{\boldsymbol{\mathsf{L}}}) \hat{\boldsymbol{q}} = {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat{\boldsymbol{u}}_{\boldsymbol r} + {\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{\boldsymbol{f}} + \hat{\boldsymbol{q}}_{\boldsymbol c}, \\ \hat{\boldsymbol{y}} = {\boldsymbol{\mathsf{H}}} \hat{\boldsymbol{q}}, \end{array} \right\} \end{equation}

which contains a windowing correction term

(3.17)\begin{equation} \hat{\boldsymbol{q}}_{\boldsymbol c} = {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{y}}}} \mathcal{F} \left\{\frac{{\rm d} w }{{\rm d} t} \boldsymbol{y}\right\}, \end{equation}

where $\mathcal {F}$ denotes the Fourier transform in time. In practice, $\hat {\boldsymbol {q}}_{\boldsymbol c}$ is constructed using available simulation data

(3.18)\begin{equation} \hat{\boldsymbol{q}}_{\boldsymbol c} \equiv {\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{y}}}} \mathcal{F} \left\{\frac{{\rm d} w }{{\rm d} t} \bar{\boldsymbol{u}}_{\boldsymbol r}\right\}, \end{equation}

with $\bar {\boldsymbol {u}}_{\boldsymbol r}$ representing the column vectors of $\bar {\boldsymbol {U}}_{\boldsymbol r}$, in (3.9a,b). The term ${\textrm {d} w }/{\textrm {d} t}$ is computed directly from the analytical formula of the windowing function used in the Welch method.

Physically, $\hat {\boldsymbol {q}}_{\boldsymbol c}$ is related to transients that are inevitably introduced when the signal is windowed (i.e. inputs necessary to match initial and final conditions of each data block) and implies that windowed spectral estimations create a mismatch between inputs and outputs, even in the case of perfectly converged statistics. Even though the windowing procedure cannot be avoided when dealing with large datasets, due to computer memory constraints, the magnitude of the correction term $\hat {\boldsymbol {q}}_{\boldsymbol c}$ can be reduced by increasing the size of the data block, which tends to proportionally decrease the value of $d w /d t$ since longer blocks imply wider windows with smaller derivatives.

3.4. Response reconstruction from inputs

From (3.7) and the spectral data matrices in (3.10a,b), we can compute the reconstructed response in Fourier space

(3.19)\begin{equation} \hat{{\boldsymbol{\mathsf{Y}}}}= \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{L}}} + \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{N}}} + \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{C}}} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{{\boldsymbol{\mathsf{u}}}} \hat{{\boldsymbol{\mathsf{U}}}}_{{\boldsymbol{\mathsf{r}}}} + {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{{\boldsymbol{\mathsf{F}}}} + {\boldsymbol{\mathsf{R}}} \hat{{\boldsymbol{\mathsf{Q}}}}_{\boldsymbol{\mathsf{c}}}, \end{equation}

where $\hat {{\boldsymbol{\mathsf{Q}}}}_{\boldsymbol{\mathsf{C}}}$ is the correction due to windowing, discussed in § 3.3, in order to obtain $\hat {{\boldsymbol{\mathsf{Y}}}} \approx \hat {{\boldsymbol{\mathsf{U}}}}_{{\boldsymbol{\mathsf{r}}}}$ by construction. In other words, the sum of all inputs with the proper correction of the distortions generated by the windowing procedure leads, in principle, to the recovery of the simulated velocity fluctuation fields, by means of the resolvent operator. This allows us to calculate separate contributions of linear mechanisms resulting from the upstream fluctuations, related to $\hat {{\boldsymbol{\mathsf{U}}}}_{{\boldsymbol{\mathsf{r}}}}$, and nonlinear receptivity due to triadic interactions, related to $\hat {{\boldsymbol{\mathsf{F}}}}$.

The cross-spectral density (CSD) matrix of $\hat {{\boldsymbol{\mathsf{Y}}}}$, can be estimated from the ensemble as

(3.20)\begin{equation} \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}} = \frac{1}{N_b} \hat{{\boldsymbol{\mathsf{Y}}}}\hat{{\boldsymbol{\mathsf{Y}}}}^H = \frac{1}{N_b} \hat{{\boldsymbol{\mathsf{Y}}}}(\hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{L}}} + \hat{{\boldsymbol{\mathsf{Y}}}}_{{\boldsymbol{\mathsf{N}}}} + \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{C}}})^H, \end{equation}

with the superscript $\{{\cdot }\}^H$ representing the conjugate transpose, and can be rewritten as

(3.21)$$\begin{gather} \left. \begin{array}{c@{}} \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{L}}}} = \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{L}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{L}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{NL}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{L}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{C}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{L}}}} \\ \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{{\boldsymbol{\mathsf{NL}}}}} = \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{L}}} {\boldsymbol{\mathsf{Y}}}_{{\boldsymbol{\mathsf{NL}}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{NL}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{NL}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{C}}} {\boldsymbol{\mathsf{Y}}}_{{\boldsymbol{\mathsf{NL}}}}}\\ \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{{\boldsymbol{\mathsf{C}}}}} = \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{L}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{C}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{{\boldsymbol{\mathsf{NL}}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{C}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{C}}} {\boldsymbol{\mathsf{Y}}}_{\boldsymbol{\mathsf{C}}}} \end{array} \right\}, \end{gather}$$

(3.22)$$\begin{gather}\hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}} = \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{L}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{N}}}} + \hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{C}}}}. \end{gather}$$

Each one of the three factors in (3.22) computes the coherence between the respective response component and the reconstructed signal. Even though these are not independent quantities, since factors contain cross-products between components, this formulation constitutes a budget measure of how each component contributes to the spectrum of the reconstructed signal.

In practice, the CSD matrix, $\hat {{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}}$, is never fully assembled due to its huge size. Since we are interested in the kinetic energies at each pair $(\beta,\omega )$, we only effectively compute the power spectral density (PSD), defined as the diagonal of the CSD matrix. Considering that the PSD is always positive and real, we obtain the relations

(3.23)\begin{gather} \boldsymbol{P_U} = {\rm Re}\{ \text{diag} (\hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{U}}}_{\boldsymbol{\mathsf{r}}} {\boldsymbol{\mathsf{U}}}_{\boldsymbol{\mathsf{r}}}})\},\end{gather}

(3.24)\begin{gather} \boldsymbol{P_Y}= {\rm Re}\{ \text{diag} (\hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{L}}}}) \} + {\rm Re}\{ \text{diag} (\hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{N}}}}) \} + {\rm Re}\{ \text{diag} (\hat{{\boldsymbol{\mathsf{C}}}}_{{\boldsymbol{\mathsf{YY}}}_{\boldsymbol{\mathsf{C}}}}) \}\nonumber\\ = \boldsymbol{\varPi_L} + \boldsymbol{\varPi_N} + \boldsymbol{\varPi_C}, \end{gather}

where $\boldsymbol {P_U} \approx \boldsymbol {P_Y}$ by construction. The term $\boldsymbol {P_U}$, computed directly from the velocity fluctuation fields of the simulation, data matrix $\hat {{\boldsymbol{\mathsf{U}}}}_{\boldsymbol{\mathsf{r}}}$ in (3.10a,b), is called statistical PSD. Then, $\boldsymbol {P_Y}$, computed through the sum of components of the input–output model is called reconstructed PSD. Because of the cross-products, $\boldsymbol {\varPi }$ components are not PSDs and can assume either positive or negative values, which are interpreted, respectively, as inflows or outflows of energy at a given pair $(\beta,\omega )$, i.e. energy exchanges between linear, nonlinear and correction components.

The equivalence between $\boldsymbol {P_U}$ and $\boldsymbol {P_Y}$ is verified numerically by the reconstruction coefficient defined as

(3.25)\begin{equation} \gamma(\beta,\omega) = \frac{\boldsymbol{P_U}^T \boldsymbol{P_Y}}{\boldsymbol{P_U}^T \boldsymbol{P_U}}. \end{equation}

Within this metric, a coefficient $\gamma \approx 1$ indicates that the reconstruction $\boldsymbol {P_Y}$ has the correct magnitude and shape, implying that the input–output model is accurate. Thus, linear and nonlinear components, $\boldsymbol {\varPi _L}$ and $\boldsymbol {\varPi _N}$ respectively, are representative in the system's response, assuming they are individually more significant than the windowing correction term, $\boldsymbol {\varPi _C}$. We may thus assess, using simulation data and resolvent analysis, the relative contribution of linear and nonlinear mechanisms in disturbance growth.

To reduce the quantity of data presented, only $\boldsymbol {\varPi _L}$ and $\boldsymbol {\varPi _N}$ components of $\boldsymbol {P_Y}$ will be displayed in corresponding results. Proof that conditions exposed in the previous paragraph are met is given by presenting the associated coefficient $\gamma$ and the magnitude of the correction component, defined as $\max |\boldsymbol {\varPi _C}|$, for each spatial direction. A more complete comparison between statistical and reconstructed PSDs for selected pairs $(\beta,\omega )$ is exposed in Appendix C.

3.5. Resolvent-based extended spectral POD

The resolvent-based extended spectral POD (RESPOD) presented in Karban et al. (Reference Karban, Martini, Cavalieri, Lesshafft and Jordan2022) is a form of extended POD (Borée Reference Borée2003) which exploits the dynamical properties of spectral POD (Towne et al. Reference Towne, Schmidt and Colonius2018) to statistically correlate inputs and outputs of a linear system in frequency space. The method can be viewed as a procedure to obtain forcing modes, ranked by their effect on the most energetic flow structures. These can be employed, for instance, in turbulence control models, as in Chevalier et al. (Reference Chevalier, Hæpffner, Bewley and Henningson2006).

Given input and output spectral data matrices, respectively $\hat {{\boldsymbol{\mathsf{F}}}}$ and $\hat {{\boldsymbol{\mathsf{U}}}}$, related linearly in the resolvent framework by

(3.26)\begin{equation} \hat{{\boldsymbol{\mathsf{U}}}} = \mathcal{R} \hat{{\boldsymbol{\mathsf{F}}}}, \end{equation}

we define an augmented state

(3.27)\begin{equation} \hat{{\boldsymbol{\mathsf{Q}}}} = \left[ \begin{array}{@{}c@{}} \hat{{\boldsymbol{\mathsf{U}}}} \\ \hat{{\boldsymbol{\mathsf{F}}}} \end{array} \right], \end{equation}

over which we apply the spectral POD method using the snapshot algorithm (Sirovich Reference Sirovich1987). By computing the weighted CSD matrix $\hat {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\mathsf{Q}}}$ in the row space of $\hat {{\boldsymbol{\mathsf{Q}}}}$, we have

(3.28)\begin{equation} \hat{{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\mathsf{Q}}} = \frac{1}{N_b} \hat{{\boldsymbol{\mathsf{Q}}}}^H \left[ \begin{matrix} {\boldsymbol{\mathsf{W}}} & 0 \\ 0 & 0 \end{matrix} \right] \hat{{\boldsymbol{\mathsf{Q}}}} = \frac{1}{N_b} \hat{{\boldsymbol{\mathsf{U}}}}^H {\boldsymbol{\mathsf{W}}} \hat{{\boldsymbol{\mathsf{U}}}}, \end{equation}

where matrix ${\boldsymbol{\mathsf{W}}}$ represents the grid quadrature weights. Next, we arrive at the eigenproblem

(3.29)\begin{equation} \hat{{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\mathsf{Q}}} \boldsymbol{\varTheta} = \boldsymbol{\varTheta} \boldsymbol{\varLambda}, \end{equation}

where $\boldsymbol {\varTheta }$ and $\boldsymbol {\varLambda }$ are, respectively, spectral POD expansion coefficients and energies of $\hat {{\boldsymbol{\mathsf{U}}}}$. The respective eigenvectors in the column space of $\hat {{\boldsymbol{\mathsf{Q}}}}$ are then given by

(3.30)\begin{equation} \tilde{\boldsymbol{\varPsi}} = \left[ \begin{matrix} \boldsymbol{\varPsi} \\ \boldsymbol{\varPhi} \end{matrix} \right] = \frac{1}{\sqrt{N_b}} \hat{{\boldsymbol{\mathsf{Q}}}} \boldsymbol{\varTheta} \boldsymbol{\varLambda}^{{-}1/2}, \end{equation}

showing that the augmented eigenvector $\tilde {\boldsymbol {\varPsi }}$ is composed of spectral POD modes $\boldsymbol {\varPsi }$ and forcing modes $\boldsymbol {\varPhi }$, which are both directly computed from the expansion coefficients $\boldsymbol {\varTheta }$ and energies/eigenvalues $\boldsymbol {\varLambda }$. Finally, by substituting equations (3.26) and (3.27) into (3.30), we get the RESPOD relation

(3.31)\begin{equation} \boldsymbol{\varPsi} = \mathcal{R} \boldsymbol{\varPhi}, \end{equation}

which shows that response and forcing modes are related by the resolvent operator.

If $\mathcal {R}$ is non-singular, $\boldsymbol {\varPhi }$ is simply the application of the inverse operator $\mathcal {R}^{-1}$ over $\boldsymbol {\varPsi }$. However, if $\mathcal {R}$ is singular, $\boldsymbol {\varPhi }$ can be shown to contain both minimal-norm forcing components, the same computed by resolvent-based estimation from Towne, Lozano-Durán & Yang (Reference Towne, Lozano-Durán and Yang2020) and Martini et al. (Reference Martini, Cavalieri, Jordan, Towne and Lesshafft2020b), and dynamically unobservable components, which are correlated to the minimal-norm forcing in the subspace spanned by the input signal (Karban et al. Reference Karban, Martini, Cavalieri, Lesshafft and Jordan2022). One thus obtains forcing modes $\boldsymbol {\varPhi }$, taken from data, which drive the observed response spectral POD modes $\boldsymbol {\varPsi }$. This yields a ranked modal decomposition of nonlinear forcing data statistics, which is a particularly useful modelling tool in cases where the response results predominantly from the nonlinear dynamics.

3.6. Spectral parameters

Spectral estimation via the Welch method is performed using blocks of $N_{FFT}=192$ realisations, a value defined via the cross-correlation procedure detailed in Blanco et al. (Reference Blanco, Martini, Sasaki and Cavalieri2022). In all the following analyses, we employ a windowing function

(3.32)\begin{equation} w(t) = \sin\left(\frac{{\rm \pi} t}{T}\right), \quad t \in [0,T], \end{equation}

and overlap of $O_{FFT} = 3/4$ between consecutive blocks, based on the guidance given in the work of Antoni & Schoukens (Reference Antoni and Schoukens2009).

8.1. Coherent structures generated by a linear mechanism

First, we focus on the pair $(\beta,\omega ) = (1.131,0.000)$ at $Tu=3.5\,\%$ and compute spectral POD modes of the data matrix $\hat {{\boldsymbol{\mathsf{Y}}}}_L$, defined in (3.19). The resulting leading mode, displayed in figure 14, features elongated streaky structures, with alternating regions of positive and negative streamwise velocity inside the boundary layer, in between counter-rotating vortices bringing high-speed flow towards the boundary layer and ejecting low-speed flow from it, in a clear instance of the lift-up effect.

Figure 14. First spectral POD mode for $(\beta,\omega ) = (1.131,0.000)$ at $Tu=3.5\,\%$, scaled by the respective eigenvalue. (a) Velocity profile of leading spectral POD mode at inlet and peak amplitude positions. (b) Real part of the leading spectral POD mode; (top) cross-section at $x = 80$; (bottom) slice at $y = 1$, inside the boundary layer.

Spatial transient growth is clear in figure 15, which displays the maximum magnitudes of each velocity component at each streamwise position. While both spanwise and vertical components only decay, the streamwise component grows before exponentially decaying. Since these structures are spatially stable and only active in upstream positions, they cannot trigger the transition to turbulence in the present simulations, even though they might contribute to it through nonlinear energy transfers.

Figure 15. Maximum velocity amplitudes of the scaled spectral POD mode for $(\beta,\omega )=(1.131,0.000)$ at $Tu=3.5\,\%$. Transient growth of streaks generated by the linear mechanism (streamwise velocity component $u$), with streamwise vortices (spanwise and vertical components $v$ and $w$) that spatially decay.

This spatial stability can be explained through transient growth theory. When the spanwise wavenumber, $\beta = 1.131$, at zero frequency, is introduced in the formulation presented in Andersson et al. (Reference Andersson, Berggren and Henningson1999), we conclude that optimal perturbations reach a maximum amplification before decaying in the streamwise direction. In a more complete analysis, we confirm that this is the case for all spanwise wavenumbers present in the FST spectrum at near-zero frequencies, as seen in figure 16(a). The largest linear amplification is found to be at the parameter corresponding to the structures previously shown in figure 13, which are, nevertheless, still less energetic than the ones presented in figure 11.

Figure 16. Amplification of incoming perturbations. Amplitude is defined as the square root of the kinetic energy at a given position $x$. (a) Optimal growth according to Andersson et al. (Reference Andersson, Berggren and Henningson1999) at spanwise wavenumbers corresponding incoming FST perturbations, computed for $\omega = 0$; (b) Amplifications measured for $Tu=3.5\,\%$. For $\beta =0.377$, not present in the OSS spectrum, the maximum optimal growth amplification is around 25, while the maximum measured amplification reaches 120.

It is worth noting that, since we can only introduce weak perturbations inside the boundary layer at the intake, which are not optimal in generating streaks through the linear amplification mechanism, the actual observed amplifications, shown in figure 16(b), are weaker than those predicted by the optimal growth theory for cases displaying linear receptivity, but larger for cases where nonlinear interactions are important, such as $\beta = 0.503$ (figure 13) and $\beta =0.377$ (figure 11, wavenumber not present in the OSS spectrum and weak in the FST energy spectrum at the intake, as seen in Appendix A).

8.2. Nonlinear coherent structures

Next, we focus on the pair $(\beta,\omega ) = (0.377,-0.003)$ at $Tu=3.5\,\%$ and follow the procedure described in § 3.5, to define an augmented state $\hat {{\boldsymbol{\mathsf{Q}}}}$, composed only by the nonlinear contributions of the model, which were identified to be the most important in this case. According to the notation of (3.19), we have

(8.1)\begin{equation} \hat{{\boldsymbol{\mathsf{Q}}}} = \left[ \begin{array}{@{}c@{}} \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{N}}} \\ \hat{{\boldsymbol{\mathsf{F}}}} \end{array} \right], \end{equation}

and, thus, spectral POD and forcing modes are linked by the relation $\varPsi = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \varPhi$.

The leading spectral POD mode, shown in figure 17, has the same overall shape found in streaks generated through linear receptivity: elongated structures, alternating streamwise velocity and counter-rotating vortices. However, significant velocity amplitudes are only present near the wall, below the position $y=15$. Besides, the characteristic wavenumber $\beta$ is smaller, such that the spacing between alternating regions is larger. As the frequency of this mode is not zero, the streaks appear inclined due to the perceived phase velocity; a similar mode is obtained for negative wavenumber, with mirrored inclination. These structures are the most energetic in cases where $Tu \geq 2\,\%$.

Figure 17. First spectral POD mode for $(\beta,\omega ) = (0.377,-0.003)$ at $Tu=3.5\,\%$, scaled by the respective eigenvalue. (a) Velocity profile of leading spectral POD mode at inlet and position $x \approx 700$. (b) Real part of the leading spectral POD mode; (top) cross-section at $x = 700$; (bottom) slice at $y = 1$, inside the boundary layer.

The streamwise evolution of streaks generated by the nonlinear mechanism is not as steep as the one generated by linear growth. Streamwise velocity amplitudes are lower and quasi-streamwise vortices are weaker than those found at the intake for linear streaks. Contrary to their linear counterpart, the amplification is sustained over the whole length of the domain. Streamwise perturbations fit an algebraic growth pattern, while quasi-streamwise vortices scale proportionally to $\sqrt {Re_x}$ as seen in figure 18. In practice, the nonlinear interactions promote the necessary conditions to counteract the dampening effect of viscosity via a continuous forcing originating from the FST outside the boundary layer.

Figure 18. Maximum velocity amplitudes. Algebraic growth of streaks generated by the nonlinear mechanism (streamwise velocity component $u$) with vortices (spanwise and vertical components, $v$ and $w$) growing proportional to $\sqrt {Re_x}$.

As previously noted by Sasaki et al. (Reference Sasaki, Morra, Cavalieri, Hanifi and Henningson2020), the velocity profile of streaks generated by the nonlinear mechanism closely matches the optimal response of the linear amplification theory, computed for $(\beta,\omega ) = (0.377,0.000)$ according to the procedure described in Andersson et al. (Reference Andersson, Berggren and Henningson1999). Especially good agreement is achieved around $Re^*=600$, based on the boundary-layer displacement thickness (see figure 19). The fact that these streaks tend to conform to the same overall shape predicted by an optimal linear mechanism and, in turn, match the experimental profiles in Westin et al. (Reference Westin, Boiko, Klingmann, Kozlov and Alfredsson1994) and Matsubara & Alfredsson (Reference Matsubara and Alfredsson2001), corroborates the conjecture concerning the existence of a strong dynamical attractor capable of ‘bringing near to itself the velocity profile under most initial conditions’, as mentioned in Luchini (Reference Luchini2000). In practice, it indicates the impossibility of asserting the linear or nonlinear nature of streaky perturbations based solely on measurements of the corresponding velocity profiles at a given position.

Figure 19. Comparison between optimal response profile, computed for $(\beta,\omega ) = (0.377,0.000)$, and leading spectral POD mode. A vertical line is drawn at $\delta _{99}$.

In a last analysis, the shape of nonlinear interactions can be analysed by looking at the first forcing mode, shown in figure 20(b). From the RESPOD formulation, spectral POD and forcing modes are phase synchronised, such that by superposing the spectral POD mode shown in figure 17, we observe that the nonlinear forcing acts by feeding streamwise vortices just outside the edge of the boundary layer, while directly weakening the streaks inside of it.

Figure 20. Forcing mode for $(\beta,\omega ) = (0.377,-0.003)$ at $Tu=3.5\,\%$, scaled by respective eigenvalue. (a) Forcing profile of leading mode at the inlet and at $x \approx 700$. (b) Real part of the leading forcing mode; (top) cross-section at $x \approx 700$; (bottom) slice at $y = 3$, inside the boundary layer.

This effect can be better described by first decomposing the nonlinear forcing data into their components in each spatial direction and then reapplying the RESPOD analysis. Following the notation adopted in (8.1), we construct

(8.2a–c) \begin{equation} \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{1}}} = \left[ \begin{array}{@{}c@{}} \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{x}}}\\ \emptyset\\ \emptyset \end{array}\right], \quad \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{2}}} = \left[ \begin{array}{@{}c@{}} \emptyset\\ \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{y}}}\\ \emptyset \end{array}\right], \quad \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{3}}} = \left[ \begin{array}{@{}c@{}} \emptyset\\ \emptyset\\ \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{z}}} \end{array}\right], \end{equation}

such that

(8.3)\begin{equation} \left\{ \begin{array}{@{}c} \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{1}}} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{1}}} \\ \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{2}}} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{2}}} \\ \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{3}}} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{3}}} \end{array} \right., \quad \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{N}}} = \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{1}}} + \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{2}}} + \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{3}}}, \end{equation}

in order to obtain the component-wise augmented state $\hat {{\boldsymbol{\mathsf{Q}}}}_{{\boldsymbol{\mathsf{a}}}}$, for which (3.30) gives the component-wise response and forcing modes,

(8.4)\begin{equation} \hat{{\boldsymbol{\mathsf{Q}}}}_{{\boldsymbol{\mathsf{a}}}} = \left[ \begin{array}{@{}c@{}} \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{1}}}\\ \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{2}}}\\ \hat{{\boldsymbol{\mathsf{Y}}}}_{\boldsymbol{\mathsf{3}}}\\ \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{1}}}\\ \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{2}}}\\ \hat{{\boldsymbol{\mathsf{F}}}}_{\boldsymbol{\mathsf{3}}} \end{array} \right] \xrightarrow{(3.30)} \tilde{\boldsymbol{\varPsi}} = \left[ \begin{array}{@{}c@{}} \boldsymbol{\varPsi_1} \\ \boldsymbol{\varPsi_2} \\ \boldsymbol{\varPsi_3} \\ \boldsymbol{\varPhi_1} \\ \boldsymbol{\varPhi_2} \\ \boldsymbol{\varPhi_3} \end{array} \right] \implies \left\{ \begin{array}{@{}c@{}} \boldsymbol{\varPsi_1} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \boldsymbol{\varPhi_1} \\ \boldsymbol{\varPsi_2} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \boldsymbol{\varPhi_2} \\ \boldsymbol{\varPsi_3} = {\boldsymbol{\mathsf{R}}}{\boldsymbol{\mathsf{B}}}_{\boldsymbol{\mathsf{f}}} \boldsymbol{\varPhi_3} \end{array} \right. .\end{equation}

In (8.4), the vectors $\boldsymbol {\varPhi _1}$, $\boldsymbol {\varPhi _2}$ and $\boldsymbol {\varPhi _3}$, are respectively the separated components in $x$, $y$ and $z$ of the forcing mode presented in figure 20. Therefore, the vectors $\boldsymbol {\varPsi _1}$, $\boldsymbol {\varPsi _2}$ and $\boldsymbol {\varPsi _3}$, displayed in figure 21, are the corresponding phase-synchronised responses to each forcing component. Indeed, the results imply that the action of $\hat{f}_x$ generates streamwise structures acting in opposition of phase concerning the streaks that are mainly generated by $\hat{f}_y+\hat{f}_z$. This streamwise dampening effect is not enough to counteract streak growth, even though $\hat{f}_x$ steadily grows, reaching larger amplitudes than the other two components. This observation is supported by results from optimal growth theory, which indicate linear amplification from vortical, $\hat{f}_y$ and $\hat{f}_z$, forcing is stronger than from pure streamwise, $\hat{f}_x$, forcing.

Figure 21. Component-wise RESPOD analysis. (a) Maximum amplitudes of the forcing components along the streamwise direction. (b) Phase-synchronised responses due to the nonlinear forcing at $x=700$. Horizontal line indicates $\delta _{99}$ thickness. (Top) Response generated by $\hat{f}_x$. (Bottom) Response generated by the composition of $\hat{f}_y+\hat{f}_z$.

Besides, the feature of opposing effects acting in a non-optimal manner to form the most energetic structures in the flow has already been observed in previous works (Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021; Nogueira et al. Reference Nogueira, Morra, Martini, Cavalieri and Henningson2021), where it was found that forcing fields computed from the nonlinear terms of the NS equations tend to project poorly into the optimal input mode computed using resolvent analysis. In the present context, the results indicate that the streaks lose energy, through the $\hat{f}_x$ component, to other wavenumbers, in what could potentially be an initial stage of streak instability and breakdown (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995).

8.3. Nonlinear receptivity mechanism

Given the spanwise length of the domain, $L_z$, and $\beta _0 = 2 {\rm \pi}/ L_z$, streaks generated by the nonlinear mechanism appear at approximately $(\beta,\omega ) = (3 \beta _0,0)$. This seems to suggest, at least for the flow case considered here, a different receptivity mechanism than the classical oblique wave set-up described in § 1, since the wavenumber in question cannot be reached by triadic interaction of the type $(\pm \beta,\omega ) \rightarrow (2\beta,0)$ (Berlin et al. Reference Berlin, Wiegel and Henningson1999; Brandt et al. Reference Brandt, Henningson and Ponziani2002). This could, however, imply an interaction between oblique waves of different wavenumbers, such as $[(\beta,\omega ),(-2\beta,\omega )] \rightarrow (3\beta,0)$.

A meaningful analysis of this mechanism would require a decomposition of the nonlinear convection term into its triadic components in both spanwise wavenumber and frequency, as described in (3.5). The identification of a set of triads linking a nonlinear pair $(\beta,\omega )$ localised inside the boundary layer to two linear pairs predominantly present in the free stream would constitute a useful data-driven approach to identify nonlinear receptivity mechanisms induced by FST in a statistically stationary set-up. Moreover, the ranked nonlinear forcing modes could be employed to characterise perturbation–perturbation interactions neglected in restricted nonlinear models (Farrell et al. Reference Farrell, Ioannou, Jiménez, Constantinou, Lozano-Durán and Nikolaidis2016).

This is not accomplished in the present work for two main reasons: (i) the databases were set up to resolve mainly the low-frequency dynamics, only a small part of the full frequency spectrum of the incoming FST perturbations, as shown in Appendix A; (ii) the spectral decomposition of less energetic pairs $(\beta,\omega )$ inevitably encounters significant windowing correction components $\boldsymbol {\varPi _{C}}$, violating the rule established in § 3.4. Arguably, a triadic analysis could be performed with a properly time-resolved database.