2.1 Flow geometry and governing equations

The case study selected is the 2D incompressible laminar flow around the backward-facing step. The configuration is sketched in figure 1. It consists of an inlet channel of height $h$ and length $L_{i}=5h$ followed by a step of height $h$ entailing a sudden expansion of height $H$ , where $H=2h$ , hence the expansion ratio is ${\it\Gamma}=h/H=0.5$ as depicted in figure 1. The Reynolds number is defined as $Re=hU_{\infty }/{\it\nu}$ , where $h$ is the inlet channel height, $U_{\infty }$ the centreline (maximum) velocity of the plane Poiseuille inlet boundary condition and ${\it\nu}$ the kinematic viscosity. The non-dimensionalized frequency defined by the Strouhal number is given by $St=fh/(U_{\infty })={\it\omega}/(2{\rm\pi})$ . Throughout the article, all the variables are non-dimensional using length $h$ , density ${\it\rho}_{\infty }$ and velocity $U_{\infty }$ as reference scales.

Figure 1.

Sketch of the flow configuration superimposed onto the velocity field in the $x$ direction of the base flow for the backward-facing step at $Re=500$ .

The flow is governed by the forced 2D incompressible NSEs,

(2.1) $$\begin{eqnarray}\partial _{t}\boldsymbol{u}+\mathscr{N}(\boldsymbol{u})=\boldsymbol{f},\end{eqnarray}$$

where

(2.2) $$\begin{eqnarray}\mathscr{N}(\boldsymbol{u})\equiv (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}+\boldsymbol{{\rm\nabla}}p-Re^{-1}{\rm\Delta}\boldsymbol{u}\end{eqnarray}$$

corresponds to the advective, pressure gradient and diffusive terms. In the entire study, the pressure field $p$ is such that velocity fields are divergence-free $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$ following the incompressibility condition. The term $\boldsymbol{f}$ represents a harmonic body forcing of the form $\boldsymbol{f}(\boldsymbol{x},t)=\boldsymbol{f}_{1}(\boldsymbol{x})\text{e}^{\text{i}{\it\omega}t}+\text{c.c.}$ with a given fixed frequency $St={\it\omega}/(2{\rm\pi})$ .

2.2 Linear transfer function

The steady solution of the NSEs,

(2.3) $$\begin{eqnarray}\mathscr{N}(\boldsymbol{U}_{B})=0,\end{eqnarray}$$

is called the base flow $\boldsymbol{U}_{B}$ and is linearly stable for the 2D backward-facing step with expansion ratio ${\it\Gamma}=0.5$ at the chosen $Re=500$ , with its first unstable mode appearing in 3D for $Re_{cr}\sim 748$ (Barkley et al. Reference Barkley, Gomes and Henderson2002; Lanzerstorfer & Kuhlmann Reference Lanzerstorfer and Kuhlmann2012). The classical approach is to study the linear response to harmonic forcing around this stable base flow. Figure 1 shows the outline of the base flow of the backward-facing step at $Re=500$ .

The exact nonlinear response can be approximated by the linear response, assuming a small amplitude of the harmonic forcing,

(2.4) $$\begin{eqnarray}[\partial _{t}\boldsymbol{u}_{1B}^{\prime }+\mathscr{L}_{\boldsymbol{ U}_{B}}(\boldsymbol{u}_{1B}^{\prime })]=\boldsymbol{f},\end{eqnarray}$$

where higher-order nonlinear terms are neglected as a first approximation, and the general operator $\mathscr{L}_{\boldsymbol{U}}(\boldsymbol{u}^{\prime })$ is the corresponding operator for the NSEs linearized around any $\boldsymbol{U}$ , i.e.

(2.5) $$\begin{eqnarray}\mathscr{L}_{\boldsymbol{U}}(\boldsymbol{u}^{\prime })\equiv (\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}^{\prime }+(\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{U}+\boldsymbol{{\rm\nabla}}p^{\prime }-Re^{-1}{\rm\Delta}\boldsymbol{u}^{\prime }.\end{eqnarray}$$

Furthermore, since the forcing is harmonic, the corresponding response $\boldsymbol{u}^{\prime }(\boldsymbol{x},t)$ will also be harmonic, $\boldsymbol{u}^{\prime }(\boldsymbol{x},t)\simeq \boldsymbol{u}_{1B}^{\prime }(\boldsymbol{x},t)=\boldsymbol{u}_{1B}(\boldsymbol{x})\text{e}^{\text{i}{\it\omega}t}+\text{c.c.}$ , and oscillate purely at the forcing frequency, due to the linearity of the operator (2.4).

The linear equation (2.4) can be rewritten formally as $\boldsymbol{u}_{1B}=\mathscr{R}({\it\omega})\boldsymbol{f}_{1}$ , where $\mathscr{R({\it\omega})}=(\text{i}{\it\omega}+\mathscr{L}_{\boldsymbol{U}})^{-1}$ is the resolvent operator. The amplitude of the response in general can be measured as the square root of the kinetic energy of the purely time-dependent fluctuation averaged over one period $T$ ,

(2.6) $$\begin{eqnarray}R=\sqrt{\frac{1}{T}\int _{0}^{T}\!\int _{{\it\Omega}}\boldsymbol{u}^{\prime 2}\,\text{d}{\it\Omega}}\,\text{d}t.\end{eqnarray}$$

Consequently, the amplitude of the linear response $\boldsymbol{u}_{1B}^{\prime }$ corresponds to the $L^{2}$ norm

(2.7) $$\begin{eqnarray}R=\sqrt{\frac{1}{T}\int _{0}^{T}\!\int _{{\it\Omega}}\boldsymbol{u}_{1B}^{\prime 2}\,\text{d}{\it\Omega}\,\text{d}t}=\sqrt{2}\Vert \boldsymbol{u}_{1B}\Vert ,\end{eqnarray}$$

as it is sinusoidal in time, defining the $L^{2}$ norm as determined by the Hermitian inner product $(\boldsymbol{a}\mid \boldsymbol{b})=\int _{{\it\Omega}}\overline{\boldsymbol{a}}\boldsymbol{\cdot }\boldsymbol{b}\,\text{d}{\it\Omega}$ , for complex fields in the domain ${\it\Omega}$ . In the same fashion, the forcing amplitude is calculated as the $L^{2}$ norm of the harmonic forcing,

(2.8) $$\begin{eqnarray}A=\sqrt{\frac{1}{T}\int _{0}^{T}\!\int _{{\it\Omega}}\boldsymbol{f}^{2}\,\text{d}{\it\Omega}\,\text{d}t}=\sqrt{2}\Vert \,\boldsymbol{f}_{1}\Vert .\end{eqnarray}$$

A natural way of measuring the amplification is the gain,

(2.9) $$\begin{eqnarray}G({\it\omega})=R/A,\end{eqnarray}$$

which is defined as the ratio between the amplitude of the response to the amplitude of the input. Subsequently, for the linear case of harmonic forcing, it reads

(2.10) $$\begin{eqnarray}G({\it\omega})=\Vert \boldsymbol{u}_{1B}\Vert /\Vert \,\boldsymbol{f}_{1}\Vert .\end{eqnarray}$$

In particular, it is relevant to determine the largest energy amplification possible at a given frequency, or in other words the optimal gain $G_{opt}({\it\omega})$ that is associated to the optimal forcing $\boldsymbol{f}_{opt}$ ,

(2.11) $$\begin{eqnarray}G_{opt}({\it\omega})=\max _{\boldsymbol{f}_{1}}\frac{\Vert \boldsymbol{u}_{1B}\Vert }{\Vert \,\boldsymbol{f}_{1}\Vert }=\frac{\Vert \mathscr{R}({\it\omega})\boldsymbol{f}_{opt}\Vert }{\Vert \,\boldsymbol{f}_{opt}\Vert }.\end{eqnarray}$$

Introducing the adjoint of the resolvent operator, the gain can be rewritten as a Rayleigh quotient of the resolvent operator and the forcing,

(2.12) $$\begin{eqnarray}G^{2}({\it\omega})=\frac{(\mathscr{R}\boldsymbol{f}_{1}\mid \mathscr{R}\boldsymbol{f}_{1})}{(\boldsymbol{f}_{1}\mid \boldsymbol{f}_{1})}=\frac{(\mathscr{R}^{H}\mathscr{R}\boldsymbol{f}_{1}\mid \boldsymbol{f}_{1})}{(\boldsymbol{f}_{1}\mid \boldsymbol{f}_{1})}.\end{eqnarray}$$

Subsequently, the optimal gain and forcing correspond to the leading eigenvalue ${\it\lambda}_{1}=G_{1}^{2}$ and eigenvector $\hat{\boldsymbol{f}}_{1}$ of the symmetric eigenvalue problem $\mathscr{R}^{H}\mathscr{R}\hat{\boldsymbol{f}}_{k}={\it\lambda}_{k}\hat{\boldsymbol{f}}_{k}$ . The procedure undertaken to obtain the optimal gain is described in more detail in Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013) and Boujo & Gallaire (Reference Boujo and Gallaire2015). Thus, we obtain for each forcing frequency an optimal forcing spatial distribution with its associated gain. A comparison of the optimal gains at different frequencies enables us to obtain the overall optimal gain $G_{opt}$ with its corresponding forcing $\boldsymbol{f}_{opt}$ . Notice that this optimization is achieved only for the linear equation, and not for the full nonlinear gain of the DNS.

2.3 Numerical methods

The linear and nonlinear NSEs are solved by applying the finite element method representing the spatial discretization of the flow fields $(u_{x},u_{y},p)$ by Taylor–Hood (P2, P2, P1) elements. The software FreeFEM++ is used to generate the domain ${\it\Omega}$ triangulation and to build all the required operators. The linear operators are solved by a sparse solver implemented directly in FreeFEM++, while the singular value decomposition is solved in Matlab. The nonlinear systems as the steady-state solutions are solved by the iterative Newton–Raphson method. The time integration of the DNS of the NSEs is calculated by applying a time scheme based on the characteristics–Galerkin method as described in Benitez & Bermudez (Reference Benitez and Bermudez2011).

The computational domain under consideration is defined in Cartesian coordinates with the origin at the step edge. The nominal domain spans from $x=-5h$ to $x=50h$ with an inlet $y=(0,h)$ and an outlet $y=(-h,h)$ as depicted in figure 1, being the same size as in Boujo & Gallaire (Reference Boujo and Gallaire2015). A 2D plane Poiseuille profile is imposed as the inlet Dirichlet boundary condition $x_{i}=-h$ . No-slip conditions are imposed at the upper $y=h$ and lower $y=-h$ boundary and an outflow boundary condition is imposed at the outlet $x_{o}=50h$ .

The size of the domain was chosen after thorough validation. The length of the domain is selected such that the base flow outlet profile returns to the fully developed Poiseuille profile, presenting a relative error between the base flow and the parabolic profile of less than 1 % in $L^{2}$ and $L_{\infty }$ norm. The influence of the entrance length on the recirculation bubble is summarized in table 1 and compared to Barkley et al. (Reference Barkley, Gomes and Henderson2002), to ensure that the chosen length is large enough and the solutions obtained are general and independent of further increase in domain size. In addition, the length of the recirculation bubble shows an excellent match for the whole range of Reynolds number compared to Barkley et al. (Reference Barkley, Gomes and Henderson2002) and Blackburn et al. (Reference Blackburn, Barkley and Sherwin2008).

Table 1.

Comparison of the recirculation bubble position for different domain geometries, i.e. variation of inlet length $L_{i}$ . The selected domain and mesh is marked in bold.

The nominal domain and mesh M1 is of 33 814 cells and 155 691 degrees of freedom (DoF), giving a gain of 7453 represented in bold in table 2, which compares quantitatively well to the optimal $G_{opt}=7480$ computed by Marquet & Sipp (Reference Marquet and Sipp2010) ensuring the quality of the methodology and mesh used. Furthermore, mesh independence is verified by doubling the number of cells by a uniform refinement, obtaining M6 with 338 247 DoF and 73 982 cells, which provides a gain of 7456, thus a relative variation of less than 0.05 %.

Table 2.

Comparison of the linear optimal gain at $St=0.075$ for different domain geometries. The selected domain and mesh is marked in bold.

First- and second-order methods (Benitez & Bermudez Reference Benitez and Bermudez2011) are compared for the time integration scheme of the DNS in table 3. The second-order method is chosen, as it presents much more robust results in terms of gain when compared to the first-order method in spite of a small increase in the computational time. The DNS time step is chosen by achieving time-step-independent results in terms of nonlinear gain as summarized in table 3. The selected time step ${\it\delta}t=0.02$ is deemed to be small enough as a trade-off between accuracy and computational time since it is located in the time-step-independent region.

Table 3.

Nonlinear saturated gain $G$ of the DNS for different time steps ${\it\delta}t$ . The nonlinear solution is forced by the optimal harmonic forcing of amplitude $A=4\times 10^{-4}$ for $Re=500$ . The selected time scheme is second order and the time step is ${\it\delta}t=0.02$ marked in bold.

5.1 Model solution

The coupled equations (5.1) of the self-consistent model must be solved iteratively for a given target amplitude of the harmonic forcing $A_{f}$ . There are two options, both starting the iterative process from the steady flow base flow $\boldsymbol{U}_{B}$ with its corresponding linear response $\boldsymbol{u}_{1B}^{\prime }$  (A 2).

5.1.2 Pseudo-time relaxation

The steady mean flow equation is modified to a pseudo-time-dependent equation (pseudo-time ${\it\tau}$ ), which represents the variation of the instantaneous mean flow from base flow $\boldsymbol{U}_{B}$ to the saturated mean flow $\boldsymbol{U}$ . The amplitude of the forcing is a smooth ramp function of time $A({\it\tau})$ that goes from 0 to the target forcing amplitude $A_{f}$ , staying at $A_{f}$ until the variations in the flow fields are negligible. At each pseudo-time step $k$ of the instantaneous mean flow $\boldsymbol{U}^{k}$ the linear forcing equation is solved and the new linearly forced response $\boldsymbol{u}_{1}^{k}$ is updated into the Reynolds stress for the next time step. The time integration is performed by a semi-implicit backward Euler method, summarized as

(5.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\boldsymbol{U}^{k+1}-\boldsymbol{U}^{k}}{{\it\delta}{\it\tau}}+\mathscr{L}_{\boldsymbol{U}^{k}}(\boldsymbol{U}^{k+1})-(\boldsymbol{U}^{k}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{U}^{k}=-2\,\text{Re}((\bar{\boldsymbol{u}}_{1}^{k}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}_{1}^{k}) & \displaystyle\end{eqnarray}$$

(5.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle [\text{i}{\it\omega}+\mathscr{L}_{\boldsymbol{U}^{k}}]\boldsymbol{u}_{1}^{k}=A({\it\tau})\boldsymbol{f}_{1,1}, & \displaystyle\end{eqnarray}$$

where $k+1$ represents the solution for the next step. We recall that all the velocity fields are divergence-free to satisfy continuity. Notice how the standard nonlinear equation for $k+1$ is approximated by its linear counterpart and a nonlinear term of the previous time step $k$ is subtracted to retrieve exactly the nonlinear mean flow equation (4.1a ) when the steady mean flow at the saturation is reached $\boldsymbol{U}_{k}=\boldsymbol{U}_{k+1}$ , hence $\mathscr{L}_{\boldsymbol{U}^{k}}(\boldsymbol{U}^{k+1})-(\boldsymbol{U}^{k}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{U}^{k}=\mathscr{N}(\boldsymbol{U}^{k})$ .

This methodology is similar to the amplitude stepping previously described, but reducing the amplitude step between iterations to a very small value ${\it\delta}A$ and adding an extra dynamical term $(\boldsymbol{U}^{k+1}-\boldsymbol{U}^{k})/{\it\delta}{\it\tau}$ . The advantages of this method are two-fold. First, the nonlinear equation of the mean flow is replaced by a linear system, requiring, however, a reduction in the amplitude step. Second, the mean flow is not locked to a steady solution for each amplitude $A_{k}$ , since the role of the previous relaxation factor ${\it\gamma}$ is substituted by a pseudo-time step that allows for physical advection of disturbances and thus stabilizing the convergence until a purely steady solution is attained. Independently of the target forcing amplitude $A_{f}$ , a minimum number of iterations is required, thus making the method not suitable for small amplitudes. The pseudo-time method should not be confused with a DNS; instead it is part of possible iterative methods that could be used to obtain a steady solution to the present nonlinear system (5.1).

5.2 Results: gain saturation and mean flow distortion

The self-consistent model predicts accurately the response and gain of the exact DNS capturing the strong nonlinear effects of the saturation as depicted in figure 5. The coupling between the mean flow and response equation inherent to the semi-linear model overcomes the limitations of the weakly nonlinear theory. The linear response (4.2) around the saturated DNS mean flow $\boldsymbol{U}_{DNS}$ at each forcing amplitude $A$ as already presented in § 4 in figure 4 is also compared to the self-consistent model in figure 5 for larger forcing amplitudes. It overpredicts the saturation (figure 5 b), providing a poorer estimation of the gain and response amplitude when compared to the present semi-linear model. This discrepancy could be in large part because the linear response is calculated a posteriori and decoupled from the mean flow equation, suggesting that the coupling between the mean flow and perturbation equation present in the self-consistent model and also in the exact DNS is relevant. Figure 5 confirms the picture of the proposed self-consistent model in which the Reynolds stress plays the main role in the saturation process while higher harmonic interactions are negligible.

Figure 5.

(a) Gain and (b) response saturation as functions of forcing amplitude $A$ for DNS (triangles) compared to self-consistent results (solid line with circles) and linear response $u_{1DNS}$ around the saturated mean flow from DNS (squares). Forcing with optimal structure at frequency $f=0.075$ and $Re=500$ .

The spatial distribution of the saturated response calculated by the self-consistent model accurately predicts the full nonlinear DNS as presented in figure 6 for a forcing $A=0.01$ . The wavelength and the position of the largest amplification of the self-consistent response match the exact DNS solution for both $u_{x}^{\prime }$ and $u_{y}^{\prime }$ fluctuations. The linear response around the DNS mean flow also resembles the full DNS structure, despite its less accurate approximation of the saturated gain (figure 5 b). It should be noted that these discrepancies are minimal when compared to the structure of the linear response around the base flow as illustrated in figure 3.

Figure 6.

Comparison of the response structure in the $x$ and $y$ directions for the self-consistent model $u_{1}$ (a,d), snapshot of the DNS $u^{\prime }$ (b,e), and linear response around the DNS mean flow $u_{1DNS}$ (c,f). High saturation level with a forcing amplitude $A=0.01$ , frequency $St=0.075$ and $Re=500$ .

As illustrated in figure 3, the Reynolds stress varies strongly with the forcing amplitude; however, thanks to the coupling present in the semi-linear model, the structure of the Reynolds stress divergence of the DNS and self-consistent model compare very closely as depicted in figure 7. In this case, the Reynolds stress divergence calculated from the linear prediction around the mean flow differs substantially from the exact DNS solution. The Reynolds stress forcing is crucial to obtain the accurate mean flow with the correct recirculation region and thus, through a proper coupling, the corresponding response.

Figure 7.

Comparison of the Reynolds stress divergence structure in the $x$ and $y$ directions for the self-consistent model $2\,\text{Re}((\bar{\boldsymbol{u}}_{1}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}_{1})$ (a,d), exact average from the DNS $\langle (\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}^{\prime }\rangle$ (b,e), and construction from the linear response around the DNS mean flow $2\,\text{Re}((\bar{\boldsymbol{u}}_{1DNS}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}_{1DNS})$ (c,f). High saturation level with a forcing amplitude $A=0.01$ , frequency $St=0.075$ and $Re=500$ .

A more quantitative comparison between the self-consistent model and DNS mean flow is summarized in figure 8, where the positions of the recirculation regions are depicted as a function of the forcing amplitude. The self-consistent model approximates precisely the exact recirculation bubbles of the nonlinear DNS. This results from the accurate calculation of the Reynolds stress by the self-consistent model which provides the correct mean flow. At high forcing the saturation is very strong and higher-order nonlinear effects, which are neglected in the semi-linear model, start playing an more important role, explaining the slight divergence between the self-consistent model results and the DNS. Nevertheless, the self-consistent model captures even the non-monotonic trend of the recirculation bubble position present in the exact DNS solutions, which would be missed by a linear prediction around the base flow.

Figure 8.

Positions of the bottom (a) and top (b) recirculation bubbles, as a function of the forcing amplitude $A$ for the DNS and self-consistent model for the optimal forcing frequency $St=0.075$ and $Re=500$ .

The frequency dependence of the saturation process of the gain starting from the linear prediction to a forcing amplitude $A=0.001$ is presented in figure 9. The comparison of the self-consistent saturated gain and the DNS for $A=0.001$ in figure 9 shows a very accurate prediction at all frequencies, with only a slight shift in the optimal frequency. More generally, figure 9 demonstrates that the most responsive frequency does not change much when nonlinear effects are included, as pointed out by the preliminary results of Marquet, Sipp & Lesshafft (Reference Marquet, Sipp and Lesshafft2010) using random noise forced DNS.

Figure 9.

Gain saturation as a function of frequency for several forcing amplitudes $A$ , ranging from 0 to $A=\cdot 10^{-3}$ in equal intervals. The DNS results are compared to the self-consistent and linear response around the base flow, from the literature. The forcing applied has its optimal structure at each frequency and $Re=500$ .