2.1. Governing equations

The non-dimensional Navier–Stokes and continuity equations for an incompressible, constant viscosity flow are

(2.1) \begin{align}& \frac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u} + \boldsymbol{\nabla } p = \frac {1}{Re _{\delta _0^*}} \nabla ^2 \boldsymbol{u} + \boldsymbol{f}, \end{align}

(2.2) \begin{align}&\qquad\qquad\qquad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} = 0, \end{align}

where $\boldsymbol{u}=[u,v,w]^T$ are velocity components in Cartesian $x$ , $y$ and $z$ directions, respectively, $p$ is the pressure, $\boldsymbol{f}$ is the forcing term and $Re _{\delta _0^*}=(U_\infty \delta _0^*)/\nu$ is the Reynolds number based on the reference length, $\delta _0^*$ , which is taken as the displacement thickness at $x=x_0$ ; $U_\infty$ is the free-stream velocity and $\nu$ is the kinematic viscosity.

For each definition of the wall waviness, two different sets of simulations are performed using the two computational domains shown in figure 1. A preliminary simulation considers the larger domain (denoted by subscript 1) and solves the two-dimensional steady-state form of the Navier–Stokes equations. Uniform velocity ( $u=1$ ) at inlet $\Omega _{i,1}$ , free-stream conditions ( $u=1, \partial v / \partial y = 0)$ on the top boundary $\Omega _{top,1}$ , no-slip conditions on the wall $\Omega _{wall}$ , symmetry boundary condition upstream of the flat plate $\Omega _{sym}$ and a standard outflow condition ( $p = 0, \partial \boldsymbol{u} / \partial x = 0$ ) at the outlet $(\Omega _{o,1})$ are used as boundary conditions for these simulations. The height and downstream extent of the larger domain are checked to be enough to ensure that the pressure on the top and outlet boundaries of the inner domain, i.e. $p_{top}$ and $p_{outlet}$ , do not change by further increasing the size of the larger domain. The solutions on the larger domain are computed using the commercial software COMSOL using a $P_3 - P_2$ finite-element formulation, i.e. third-order elements for the velocity field and second-order elements for the pressure.

The rest of the simulations are performed in a reduced computational domain (enclosed by red lines in figure 1), hereafter referred to as the DNS domain, using the spectral element method solver Nek5000 (Deville et al. Reference Deville, Fischer and Mund2002; Fischer et al. Reference Fischer, Lottes and Kerkemeier2008). The $P_N - P_{N-2}$ formulation in Nek5000 with seventh-order elements for velocity and fifth-order elements for pressure is used. In figure 1 the height and length of the DNS domain are denoted by $H$ and $L$ , respectively. The boundary conditions for the DNS domain are summarised in table 1. For three-dimensional simulations, periodic boundary conditions are imposed on the (spanwise) lateral boundaries. For each DNS case (both two- and three-dimensional simulations), the prescribed profiles for $[u,v]_{inlet}$ , $[u,p]_{top}$ and $p_{outlet}$ (table 1) are obtained from the simulation on the larger domain. Note that, the solutions obtained with COMSOL on the larger domain are only used to provide boundary conditions for DNS simulations. This choice of boundary conditions is based on two reasons: (i) due to the presence of the waviness, the velocity at the inlet might differ from the Blasius solution; and (ii) by setting $p_{top}$ and $p_{outlet}$ , the two-dimensional base flows become independent of the height and downstream extent of the DNS domain (Shahriari & Hanifi Reference Shahriari and Hanifi2016). The effect of the waviness on the inlet profiles (for case E – see § 3.2 for the definition of different cases) is shown in figure 2. Finally, a fringe region is added at the outlet of the DNS domain to avoid reflections from the outlet into the domain.

Table 1. Boundary conditions imposed on the DNS domain.

Direct numerical simulations using the smaller domain are performed in this work with different objectives. One is the computation of steady-state solutions of the Navier–Stokes equations to be subsequently subject to linear stability analyses. If the nonlinear system of equations is globally unstable, no steady-state solution can be found using time-marching techniques. Instead, the steady-state solution of the unstable system can be found using techniques such as selective frequency damping (SFD) (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006) or Boostconv (Citro et al. Reference Citro, Luchini, Giannetti and Auteri2017). The Boostconv method aims to enhance the convergence of iterative solvers by modifying the part of the eigenspectrum dictated by the unstable or slowly decaying modes. This is achieved by adjusting the residual at each step to ensure faster reduction of the residual in the subsequent iteration. This method is used in this work to find the spanwise-homogeneous two-dimensional base flows, denoted by $\boldsymbol{u}$ $_0=[u_0,v_0,0]^T$ . SFD is used in § 4.2.1 to find the spanwise-modulated fully three-dimensional base flows, denoted by $\boldsymbol {u}_{3D}=[u_{3D},v_{3D},w_{3D}]^T$ . SFD adds a forcing term to the momentum equations (using term $f$ in (2.1)) that, in practice, acts as a low-pass filter efficiently damping the unsteadiness present in the flow.

Figure 2. Effect of waviness on the inlet boundary condition. The DNS inlet profile is for case E (see tables 2 and 3).

2.2. Linearised equations

Considering a steady solution of the governing equations $\bar {\boldsymbol{q}} = [\bar {\boldsymbol{u}},\bar {p}]^T$ as base flow, the linearised Navier–Stokes and continuity equations around $\bar {\boldsymbol{q}}$ can be written as

(2.3) \begin{align}& \frac {\partial \boldsymbol{u}'}{\partial t} + (\boldsymbol{u}' \boldsymbol{\cdot }\boldsymbol{\nabla })\bar {\boldsymbol{u}} + ( \bar {\boldsymbol{u}} \boldsymbol{\cdot } \boldsymbol{\nabla })\boldsymbol{u}' +\boldsymbol{\nabla } p' = \frac {1}{Re _{\delta _0^*}} \nabla ^2 \boldsymbol{u}', \end{align}

(2.4) \begin{align}&\qquad\qquad\qquad\qquad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}' = 0, \end{align}

where $(')$ denotes perturbation quantities. When the linear solver in Nek5000 is used, to ensure consistency with the boundary conditions satisfied by $\bar {\boldsymbol{q}}$ , homogeneous Dirichlet boundary conditions are applied for $\boldsymbol{u}$ $'$ on $\Omega _{i,2}$ and $\Omega _{wall}$ , while on $\Omega _{top,2}$ and $\Omega _{o,2}$ , the respective boundary conditions are

(2.5) \begin{equation} u'=w'=0, \quad \frac {1}{Re _{\delta _0^*}}\frac {\partial v'}{\partial y} - p' = 0 \end{equation}

and

(2.6) \begin{equation} \frac {1}{Re _{\delta _0^*}}\frac {\partial u'}{\partial x} - p' = 0, \quad \frac {\partial v'}{\partial x} = \frac {\partial w'}{\partial x} = 0. \end{equation}

2.3. Modal stability analysis

Considering modal analysis, the perturbations are written in the form of

(2.7) \begin{equation} \boldsymbol{q}'(x,y,z,t) = \hat {\boldsymbol{q}}(x,y,z)e^{-i{\Lambda t}} + c.c., \end{equation}

where $\Lambda =\omega + i\sigma$ and $c.c.$ denotes the complex conjugate. By substituting (2.7) into (2.3) and (2.4), a generalized eigenvalue problem of the form

(2.8) \begin{equation} \unicode{x1D647} \hat {\boldsymbol{q}}=\Lambda \unicode{x1D64D} \hat {\boldsymbol{q}} \end{equation}

is found, where

(2.9) \begin{equation} \unicode{x1D64D} = \begin{pmatrix} \unicode{x1D644} & 0\\ 0 & 0 \\ \end{pmatrix}, \quad \unicode{x1D647} = \begin{pmatrix} -\boldsymbol{\nabla }\bar {\boldsymbol{u}} - \bar {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla } + \frac {1}{Re _{\delta _0^*}} \nabla ^2 & -\boldsymbol{\nabla } \\ \boldsymbol{\nabla }\boldsymbol{\cdot }& 0 \\ \end{pmatrix}. \end{equation}

The eigenvalues $\Lambda = \omega + i\sigma$ determine the asymptotic behaviour of the system to small amplitude disturbances, where $\sigma$ and $\omega$ are the exponential growth rate and frequency of perturbations, respectively. If $\sigma \gt 0$ , the system is globally unstable and any small amplitude perturbation will grow exponentially over time.

To obtain the two-dimensional global modes of the two-dimensional base flows $\bar {\boldsymbol{q}}$ , a matrix-free approach (Bagheri et al. Reference Bagheri, Åkervik, Brandt and Henningson2009) is employed that uses the implicitly restarted Arnoldi method (Lehoucq et al. Reference Lehoucq, Sorensen and Yang1998) implemented in Nek5000. In this case, the two-dimensional base flows ( $\boldsymbol{u_0}$ ) are obtained with Nek5000 (on the DNS domain) using Boostconv. The base flows will be introduced in § 3.2, while the stability analysis results corresponding to this part are presented in § 4.1.

Figure 3. Mapping of the physical domain to the computational domain used in the matrix-forming eigenvalue analysis. The computational domain corresponds to the shaded area of the physical domain.

To obtain the three-dimensional global modes of two-dimensional base flows, a matrix-forming approach is used. The base flows here are obtained by interpolating the base flows calculated using Nek5000 on a smaller subdomain (shaded area in figure 3). In this case, instead of solving the three-dimensional eigenvalue problem, perturbations can be assumed to be periodic in the spanwise direction, with spanwise wavenumber $\beta$ . In this case, the ansatz (2.7) is rewritten as

(2.10) \begin{equation} \boldsymbol{q}'(x,y,z,t) = \hat {\boldsymbol{q}}(x,y)e^{i(\beta z - {\Lambda t})} + c.c. \end{equation}

By substituting the ansatz (2.10) into (2.3) and (2.4), a generalized eigenvalue problem similar to (2.8) is obtained. The corresponding terms for linear operators $\unicode{x1D647}$ and $\unicode{x1D64D}$ of the resulting eigenvalue problem can be found in Rodríguez & Theofilis (Reference Rodríguez and Theofilis2009). Fourth-order central finite differences are used to discretise the derivatives of the operators. Moreover, in order to account for the curvature of the geometry as well as non-orthogonal grids, the analytical mapping

(2.11) \begin{equation} \xi = x, \quad \eta = {H_0}\dfrac {y(x)-y_w(x)}{y_{top}(x) - y_w(x)}, \end{equation}

is used to map the physical domain $(x,y)$ into the computational one ${(\xi, \eta })$ as shown in figure 3. Homogeneous Dirichlet boundary conditions are imposed for perturbations at the inlet, top boundary of the computation domain and at the wall, while homogeneous Neumann boundary conditions are imposed at the outlet of the domain. In (2.11), $H_0$ is the height of the subdomain, while $y_w(x)$ and $y_{top}(x)$ are the coordinates of the wall and top boundary in the physical domain, respectively. Note that the top boundary of the physical subdomain chosen for the eigenvalue problem (grey area in figure 3) is formed by simply moving the shape of the wall $H_0$ units in the $y$ direction. This particular choice of the top boundary of subdomain for eigenvalue analysis guarantees an orthogonal grid in the computational domain, and both the physical and computational domains will have the same grid clustering near the wall without any further transformations.

Finally, the resulting two-dimensional eigenvalue problem is solved in parallel with sparse linear-algebra techniques using the shift-and-invert Krylov–Schur algorithm implemented in SLEPc (Hernandez et al. Reference Hernandez, Roman and Vidal2005). Also, the MUMPS package (Amestoy et al. Reference Amestoy, Duff, L’Excellent and Koster2000) is used for performing the lower–upper decomposition of the sparse matrices. The results corresponding to this section are presented in § 4.2. This approach is validated with Nek5000 as shown in Appendix A.2.

3.1. Onset of flow separation

The wall waviness presents two consecutive changes in curvature that introduce favourable pressure gradient and APG regions, and cause flow separation if the intensity of the APG is large enough. Different works proposed boundary layer separation criteria based on dimensionless parameters that consider the thickness of the boundary layer, the kinematic viscosity and a measure of the streamwise velocity gradient at the boundary layer edge (Thwaites Reference Thwaites1949; Gaster Reference Gaster1966); thus, a priori knowledge of velocity profiles is required. Ultimately, these quantities are a consequence of the geometrical parameters of the waviness ( $h,\lambda$ ) and the Reynolds number based on the boundary layer thickness. This suggests an alternative criterion for the prediction of the boundary layer separation for the present configuration, based on the dimensionless numbers $h/\delta$ and $Re _{\delta _{x_c}^*}={(U_\infty {\delta _{x_c}^*})}/{\nu }$ . Here $\delta$ and $\delta ^*$ are the boundary layer and displacement thickness of the corresponding Blasius solution at $x = x_c$ , respectively. To find such a dependence, the wavelength $\lambda$ for which the wall shear stress vanishes has been determined for simulations with different Reynolds numbers and waviness heights. These are shown in the left panel of figure 4. Next, all the curves are scaled with $(Re _{\delta _{x_c}^*})^{0.5}$ , obtaining an excellent collapse as shown in the right panel of figure 4. Then, a power-law curve is fitted to the data to obtain an equation to predict the waviness parameters that causes flow separation:

(3.1) \begin{equation} \frac {\delta }{\lambda }(Re _{\delta _{x_c}^*})^{0.5}=0.1028\left (\frac {h}{\delta }\right )^{-1.2892}+0.1611. \end{equation}

For any combination of parameters that lies above this curve, flow separation is expected to happen.

Figure 4. Scaling of waviness parameters for incipient separation.

Table 2. Geometrical parameters defining the cases considered in this work.

3.2. Base flow cases

Equation (3.1) allows constructing different geometries resulting in LSBs with different (not known a priori) length, thickness and reverse flow intensity. For all the geometries considered in the rest of this work, the Reynolds number based on local Blasius displacement thickness at the centre of the waviness is $Re _{\delta _{x_c}^*}=1720.8$ , which corresponds to $Re _{{x_c}}={(U_\infty x_c)}/{\nu }=1 \times 10^6$ . For all cases, the location of the inlet $(x=x_0)$ is chosen such that the Reynolds number based on the local Blasius displacement thickness is $Re _{\delta _{x_0}^*}=1332.925$ . The value of $\delta _{x_0}^*$ is taken as the reference length to make all coordinates non-dimensional. Finally, the origin of the dimensionless streamwise coordinate is displaced to the inlet of the DNS computational domain, so that $x=0$ is the inlet coordinate and $x_c \approx 300$ . For reference, the dimensionless Blasius boundary layer thickness at $x_c$ is $\delta _{x_c}=3.7512$ . Table 2 summarises the geometrical parameters for different geometries considered; these cases are also shown on the separation diagram in figure 5. Separated flow is present for all cases except case B, which is at the edge of separation. Note that, all two-dimensional base flows are obtained with a DNS domain size of $H=700$ and $L=1400$ .

Figure 5. Geometrical parameters defining the cases considered in this work. The dashed line shows the separation curve for $Re _{\delta _{x_c}^*}=1720.8$ .

Figure 6. Base flows for cases A, C and E. The colour map corresponds to streamwise velocity $u$ . Here ${\rm d}C_p/{\rm d}s$ is shown by a dashed black line and the values are indicated by the right vertical axis. White and red dashed lines show reverse flow $y_r$ and dividing streamlines $y_d$ , respectively. The figure is not drawn to scale.

Figure 6 shows the streamwise velocity $u$ for cases A, C and E. The figure is not drawn to scale in order to enhance the visibility of the recirculation bubble. The derivative of the pressure coefficient $C_p$ with respect to the wall arc length $s$ , $dC_p/ds$ , is shown by the dashed black line. The red dashed line shows the dividing streamline $y_d$ , defined as the $y$ coordinate below which the streamwise mass flux is zero, i.e. $\int _{y_{wall}}^{y_d}u(x,y){\rm d}y=0$ . The white dashed line corresponds to the coordinate $y_r$ that delimits the reverse flow region ( $u\lt 0$ ). Table 3 summarises different characteristics of the separation bubble for all cases.

Table 3. Maximum reverse flow $(u_{rev})$ , $x$ coordinate of separation $(x_s)$ and reattachment $(x_r)$ points, length of separation bubble $(L_s = x_r - x_s)$ , momentum thickness at separation point $(\theta _s)$ , maximum vertical distance between the wall and dividing streamline $(h_1)$ , maximum vertical distance between the wall and separation point $(h_2)$ , Reynolds number based on momentum thickness at the separation point ( $Re _{\theta, s}$ ) and Reynolds number based on the length of the separation bubble $(Re _{L})$ . For all cases, $Re _{\delta _{x_c}^*}=1720.8$ .

Figure 6, illustrates how the change of wall curvature, especially near the reattachment, plays a significant role in the shape and length of the separation bubble. The flow experiences an acceleration in the first part of the waviness, followed by an APG region that results in flow separation. A pressure plateau is evident within the recirculating region, except for the proximity of the reattachment region. For cases A and E, near the rear end of the bubble the flow experiences a sudden deceleration and the dividing streamline bends towards the wall. The reattachment occurs either at the end point of the wall waviness or downstream, in the flat plate. After reattachment, the flow accelerates until a zero pressure gradient region is recovered over the flat plate. Conversely, in case C where the reattachment is upstream of the end of the waviness, the flow experiences an acceleration after the reattachment.

Regarding the maximum value of reverse flow, both $h$ and $\lambda$ play a role. This can be seen by comparing maximum reverse flow magnitude for cases G, E and H, where $h/\lambda$ is the same; however, the reverse flows are different. The reverse flow increases by increasing the height of waviness while keeping constant $\lambda$ . At a constant height of the waviness, increasing $h/\lambda$ (i.e. decreasing $\lambda$ ), increases the reverse flow magnitude. The effect of height of the waviness is found to be more pronounced: $h/\lambda$ for case A is almost $16\,\%$ higher than that for case H, but case H has a higher maximum reverse flow than case A.

The wall curvature also affects the streamlines of the flow, especially near the rear part of the waviness where the flow undergoes a sudden deceleration. This is shown in figure 7 that depicts the streamlines for cases G, A, E and H. Among the four cases shown in this figure, the curvature of the streamlines in case G does not change sign inside the recirculation bubble. The same happens for cases C, D and F (not shown here). On the other hand, the curvature of streamlines in cases A, E, H and I (not shown here) changes sign in the vicinity of the reattachment point. This change in the sign of streamline curvature is found to directly impact the three-dimensional global centrifugal instabilities, as will be discussed in § 4.2. Note that in the base flow for case H a second recirculation region inside the main one is formed, as shown in figure 7(e). The same streamline topology has been reported for different two-dimensional steady solutions of the Navier–Stokes equations recovering fully laminar recirculation bubbles (e.g. Alizard et al. Reference Alizard, Cherubini and Robinet2009; Balzer & Fasel Reference Balzer and Fasel2016). This seems to be an effect of the rapid increase of pressure gradient near the reattachment point. In what follows, the spanwise-homogeneous (two-dimensional) base flows presented in this section will be referred to by $\boldsymbol{u}_0=[u_0,v_0,0]^T$ .

Figure 7. Flow streamlines near reattachment point for cases (a) G, (b) A, (c) E and (d) H. The figure is not drawn to scale. (e) Velocity vectors within the marked dashed rectangular region in case H.

4.1. Two-dimensional instability: vortex shedding

In this section the results for global stability analysis for two-dimensional perturbations (introduced in § 2.3) are presented. The numerical details and domain sensitivity analysis are given in Appendix A.1. Figure 8 shows the eigenspectra for all the cases. Note that the eigenspectrum being symmetric about $\omega = 0$ , only the side of it with positive frequencies is shown. The eigenspectra show two different families of eigenmodes, characterised as low-frequency ( $\omega \lt 0.06$ ) and high-frequency ( $\omega \gt 0.06$ ) regions, respectively. For case B, which is the case without flow separation, only the low-frequency part of the eigenspectrum appears. Among all cases studied, only cases E, F and H are globally unstable with respect to two-dimensional instabilities; a self-sustained oscillator-type instability is expected to happen for these cases in nonlinear simulations without external forcing. The peak reverse flow for these unstable cases is $u_{rev}\approx 10-12\,\%$ , which is lower than the threshold of $u_{rev}\approx 16 \,\%$ (Rist & Maucher Reference Rist and Maucher2002; Alam & Sandham Reference Alam and Sandham2000), and the inflection point does not cross the dividing streamlines of these bubbles (as proposed by Avanci et al. Reference Avanci, Rodríguez and Alves2019), suggesting that the global instability does not originate from the local absolute inflectional instability, but by fully non-parallel effects.

Figure 8. Eigenspectra of two-dimensional global eigenmodes ( $\beta = 0$ ) of the two-dimensional base flows.

Figure 9. Real part of the streamwise velocity component of global eigenmodes for case E, corresponding to the four eigenvalues marked with the letters $a$ , $b$ , $c$ and $d$ in figure 8.

The growth rate and frequency of high-frequency modes show a higher sensitivity to the geometry of the waviness. This is in agreement with the results of Ehrenstein & Gallaire (Reference Ehrenstein and Gallaire2008) for a separated flow over a hump. Figure 9 shows the real part of the streamwise velocity component of global modes for case E, corresponding to the four eigenvalues marked with the letters $a$ , $b$ , $c$ and $d$ in figure 8. Figure 9 confirms that that low-frequency and high-frequency modes belong to different families. Low-frequency modes, corresponding to convectively unstable waves, resemble Tollmien–Schlichting waves downstream of the wall waviness. The corresponding eigenfunctions are also similar to eigenfunctions of a family of low-frequency modes found downstream of the reattachment point of an LSB caused by a smooth hump, as shown by Ehrenstein & Gallaire (Reference Ehrenstein and Gallaire2008). The modes in the high-frequency family exhibit a notable spatial amplification in the separated shear layer and inside the recirculation region, upstream of the reattachment point. The spatial shape of all the high-frequency modes (either stable or unstable) shares the same features; however, their spatial extension downstream of the reattachment point decreases from eigenmode to eigenmode as their frequency increases (modes $c$ and $d$ ). The eigenmodes shown in figure 9 are similar to the two-dimensional eigenmodes found in the LSBs caused by a smooth cavity (Åkervik et al. Reference Åkervik, Hœpffner, Ehrenstein and Henningson2007) and smooth hump (Ehrenstein & Gallaire Reference Ehrenstein and Gallaire2008).

4.2. Three-dimensional instability: base flow three-dimensionalisation

Since the base flows $\boldsymbol{u}_0$ are homogeneous in the $z$ direction, their three-dimensional instability is studied using global stability analysis of perturbations of the form (2.10). Details of the numerical set-up and cross-validations using Nek5000 are given in Appendix A.2.

For an LSB on a flat plate, Theofilis et al. (Reference Theofilis, Hein and Dallmann2000) reported the existence of a single steady unstable centrifugal three-dimensional eigenmode for a specific spanwise wavenumber. Recently, Rodríguez et al. (Reference Rodríguez, Gennaro and Souza2021a ) showed that this eigenmode is responsible for spanwise modulation of an initially two-dimensional base flow and can potentially lead to transition through a self-excited secondary instability mechanism. Here, on a curved surface, we show that for a specific $\beta$ , more than one zero-frequency unstable centrifugal mode can coexist. The effect of spanwise wavenumber on the growth rate of the leading unstable eigenmode(s) is shown, for all cases, in figure 10. Note that all the eigenmodes shown in this figure have zero frequency. In the case of having two unstable eigenmodes for a specific $\beta$ , the lower-growth-rate modes are shown by the same symbol as the higher-growth-rate modes, but are connected with a dashed line. No more than two unstable eigenmodes have been found in any of the cases.

Figure 10. Dependence of the maximum growth rate of steady ( $\omega =0$ ) eigenmodes on the spanwise wavenumber for different cases. In the case of coexistence of two eigenmodes at one spanwise wavenumber, the growth rate of the second mode is shown by a dashed line.

Figure 11. (a) Stability diagram in $h/\delta$ – $\delta / \lambda$ parameter space. Plus and circles show the stability characteristics of two- and three-dimensional instabilities, respectively (blue: stable, red: unstable). The black dashed line shows the separation curve for $Re _{\delta _{x_c}^*}=1720.8$ (see figure 5). (b) Stability diagram of three-dimensional instabilities in the $Re ^{1/2}_{hh}$ – $h / \lambda$ parameter space. The blue circle shows the stable case, black circles show the cases with one unstable eigenvalue and red circles show the cases where two unstable eigenvalues are present for a particular value of $\beta$ .

Figure 11 summarises the stable and unstable cases in two different parameter spaces. Figure 11(a) shows the stability characteristics of cases in $h/\delta$ – $\delta / \lambda$ parameter space, where the plus and circle symbols are for two- and three-dimensional instabilities, respectively. Blue and red colours also show stable and unstable cases, respectively. From the figure it is clear that, at a fixed height of the waviness, by decreasing the wavelength of the waviness, i.e. increasing $\delta / \lambda$ , three-dimensional instabilities become unstable prior to two-dimensional instabilities. Figure 11 (b) shows the stability characteristics of three-dimensional instabilities in the $Re ^{1/2}_{hh}$ – $h / \lambda$ parameter space, where $Re _{hh}=u(h)h/\nu$ is the local Reynolds number at the $x$ location of the maximum height of the waviness $h$ and $u(h)$ is the value of the Blasius velocity profile at height $h$ (Von Doenhoff & Braslow Reference Von, Albert and Braslow1961; Bucci et al. Reference Bucci, Cherubini, Loiseau and Robinet2021). Although analysing more cases is required to properly identify the neutral curves at this parameter space, from the figure it is clear that at a fixed $Re _{hh}$ , by increasing the aspect ratio $h / \lambda$ of the waviness, first the flows become unstable with respect to one unstable global mode (black circles) and then with respect to two global eigenmdoes (red circles) for a fixed value of spanwise wavenumber.

The existence of multiple steady centrifugal eigenmodes is found to be directly related to the number of sign changes in the curvature of the recirculating streamlines of the base flow. As explained in § 3.2 and shown in figure 7, cases A, E, H and I present more than one change of curvature sign due to the presence of the curved wall. For these four cases, as figure 10 shows, the coexistence of two unstable centrifugal modes is possible in a finite range of $\beta$ . Note that the unstable range of $\beta$ for the higher-growth-rate eigenmode is almost twice that for the lower-growth-rate mode.

Figure 12 shows the real part of the streamwise $u$ (left panels) and $v$ (right panels) velocity components for the two unstable eigenmodes of case A with $\beta =0.7$ , and the single unstable eigenmode for case G with $\beta =0.4$ ; the top and middle panels show the eigenfunctions for the higher-growth-rate (M1) and lower-growth-rate mode (M2) for case A, and the bottom panels show the same eigenfunctions for the unstable mode (M1) in case G. Flow streamlines in the vicinity of the change of streamline curvature are also plotted as thin grey lines in the figure. The higher-growth-rate mode is very localised in the region where the curvature of streamlines changes. For a given $\beta$ , the growth rate of this mode, which will be referred to as the localised mode, is one order of magnitude higher than that of the other mode, which will be referred to as the spread mode. It should be noted that, at a fixed $h/\delta$ when both the localised and spread modes coexist, the growth rate of both modes increases with decreasing wavelength of the waviness, i.e. increasing $\delta / \lambda$ , as can be seen from figure 10 (cases A and I). The spatial structure of the spread mode extends over the whole recirculation bubble and shares the same features as the three-dimensional global mode found in LSBs formed on a flat plate (Theofilis et al. Reference Theofilis, Hein and Dallmann2000; Cherubini et al. Reference Cherubini, Robinet, De Palma and Alizard2010; Rodríguez & Theofilis Reference Rodríguez and Theofilis2010), smooth humps (Gallaire et al. Reference Gallaire, Marquillie and Ehrenstein2007) and S-shaped ducts (Marquet et al. Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009). When only one unstable mode is found (e.g. case G), the eigenfunction identifies it as a spread mode. There are two important differences between the properties of the localised and the spread eigenmodes that are consistently observed in their respective eigenfunctions for all the cases. First, unlike the spread mode, the localised mode nearly vanishes along the dividing streamline (the red line in figure 12), except for a small region centred above the dominant recirculation centre. Second, the $u$ component of the localised mode changes sign inside the bubble and near the reattachment point. Conversely, the $u$ component of the spread mode (even when it coexists with the localised mode) changes sign only downstream of the reattachment point; the flow disturbance induced by this mode produces a spanwise modulation of the entire recirculation region (Rodríguez & Theofilis Reference Rodríguez and Theofilis2010) that is found here to be more relevant to the subsequent nonlinear evolution than that induced by the localised mode, as is discussed in the following sections.

Figure 12. Real part of the $u$ component (left panels) and $v$ component (right panels) of global modes for case A with $\beta =0.7$ and case G with $\beta =0.4$ . The most unstable mode for case A is denoted by M1 and the second unstable mode is denoted by M2. The red solid line shows the dividing streamline for each case. The flow streamline in the rear part of the bubble is shown by grey lines. Red and blue colours show positive and negative values, respectively. The figure is not drawn to scale.

4.2.1. Nonlinear saturation of the three-dimensional instability

To investigate the nonlinear response and evolution of the system due to the growth of global unstable modes until nonlinear saturation, a set of DNS is performed. The simulations are initialised from the corresponding two-dimensional base flow for each case. Note that only cases in which all two-dimensional eigenmodes were stable are considered. Table 4 shows the cases considered. Periodicity is imposed on the spanwise boundaries. The spanwise extent of the computational domain ( $L_z=2\pi / \beta$ ) in the simulations is chosen in each case to impose a fundamental wavenumber $\beta$ of particular relevance based on figure 10. For case G, the spanwise size of the domain ( $L_z$ ) is chosen such that the corresponding spanwise wavenumber of domain ( $\beta =2 \pi / L_z$ ) matches the spanwise wavenumber of the most unstable mode for this case, i.e. $\beta =0.4$ , as shown in figure 10. For case A, two spanwise sizes are considered, $L_z=2 \pi /1.2$ and $L_z=2 \pi /0.7$ , where 0.7 is the spanwise wavenumber of the most unstable mode for case A, and 1.2 is the spanwise wavenumber where only one localised mode is unstable and the spread mode is stable (see figure 10). In the following, the value of $\beta$ indicates the spanwise size of the domain, i.e. $L_z=2\pi /\beta$ , in each case. Convergence studies are reported in Appendix A.3.

Table 4. Cases considered for nonlinear simulations. Here $L_z=2 \pi /\beta$ denotes the spanwise size of the domain.The forth column describes the number and type of unstable eigenmodes at the fundamental spanwise wavenumber $\beta$ ; SP and LM denote the spread and localised modes, respectively. The other columns show the maximum reverse velocity for initial two-dimensional flow, $u_0$ ; three-dimensional saturated flow, $u_{3D}$ ; spanwise averaged saturated flow, $\bar {u}_{2D}$ ; and the peak negative base flow distortion at nonlinear saturation, $u_{bd}(t_{sat})$ .

Figure 13. Time evolution of the peak of the absolute value of streamwise perturbation ( $\|u'\|_{\infty }$ ) (a) and integrated kinetic energy of perturbations (b) for case A with $\beta =0.7$ (solid black line), case A with $\beta =1.2$ (dashed black line) and case G with $\beta =0.4$ (solid blue line). The blue dashed line correspond to case G after applying SFD. The inset in panel (a) shows the oscillations that happen for case G before using SFD. The red dashed lines show the growth rate predictions based on LST. Note that $t_0$ is the transient time that is different for each case. Note that the simulations are run for a much longer time than the part shown in this figure.

Figure 13(a) shows the temporal evolution of the absolute value of the peak streamwise perturbation component $\|u'\|_{\infty }$ for each case, where ${u}'(x,y,z,t) = {u}_{3D}(x,y,z,t) - {u}_0(x,y)$ is the difference between the instantaneous velocity and the two-dimensional base flow. All the simulations present an initial transient that is not shown here. In what follows, the time $t' = t-t_0$ is considered, where $t_0$ is an arbitrary time after the initial transient but before important dynamics appears. The evolution of disturbance kinetic energy $E(t) = \sqrt {(\boldsymbol{u}' \cdot \boldsymbol{u}')}/2$ is shown in figure 13(b). Note that the energy is integrated in the whole domain. The instantaneous growth rate of the disturbances can be computed from the DNS data using the disturbance energy as

(4.1) \begin{equation} \sigma _{DNS} = \frac {1}{E(t')}\frac {{\rm d}E(t')}{{\rm d}t'}. \end{equation}

In figure 13 the red dashed lines show the growth of the most unstable mode for each case obtained from linear theory (figure 10). In the following, an overview of the nonlinear saturation process for each of these three cases will be given.

1. (i) Case G, $\beta =0.4$ : this case aims to investigate the dynamics when only a single spread mode is present. The perturbation experiences an exponential growth for a long time (at $t'\approx 3300$ , $\sigma _{DNS}\approx 0.001200$ and $\sigma _{LST}\approx 0.001204$ ) followed by a slower growth due to nonlinear interactions of higher harmonics of the spanwise wavelength starting from $t'=(t-t_0) \approx 5000$ , as is clear from the reduced slope of the blue line in figure 13(b). At $t' \approx 9000$ , oscillations starts to appear in the shear layer close to reattachment point, as shown in the inset of figure 13(a). These oscillations are an indicator of a self-excited secondary instability that could lead to transition (Rodríguez et al. Reference Rodríguez, Gennaro and Souza2021a ). Note that, figure 13(a) shows the time evolution of the peak of the absolute value of $u'$ in the whole domain, i.e. $\|u'\|_{\infty }$ . Thus, the slope of amplitude growths in figure13(a) should not be considered as the growth rate of perturbations, since the point corresponding to peak perturbation might vary over time. However, this figure is useful to identify self-excited oscillations in the domain when the kinetic energy of the domain is not a good metric (since the oscillations do not appear in such an integral metric), and the location of oscillations is not known a priori. The subsequent transition process will be investigated in § 6. Here, to obtain the saturated base flow for this case for further analysis, when oscillations start to appear in the flow, SFD is activated. The saturated base flow is assumed to be converged when no significant change up to the seventh digit was found in $\| u' \|_\infty$ . The blue dashed line in figures 13(a) and 13((b) shows the evolution of $\|u'\|_{\infty }$ and kinetic energy after activating the SFD.

2. (ii) Case A, $\beta =1.2$ : this case aims to investigate the dynamics when only a single localised mode is present. The black dashed line in figure 13 shows that, after the initial transient, perturbations grow exponentially in the domain until nonlinear saturation gives rise to a steady flow. As opposed to the previous case, new fluctuations do not appear after the saturation. The growth rate of the exponential growth phase is approximated at $t' \approx 550$ . For this case, $\sigma _{DNS} \approx 0.003473$ which is close to that obtained using linear stability analysis, $\sigma _{LST}=0.003488$ .

3. (iii) Case A, $\beta =0.7$ : in this case, the localised mode and its first harmonic ( $\beta =1.4$ ) coexist with an unstable spread mode with $\beta = 0.7$ . The same scenario as in the previous case is recovered. However, the growth rate of perturbations is larger than that for $\beta =1.2$ , as figure 10 also suggests. Surprisingly, for this case, the saturated state is achieved for later times compared with $\beta =1.2$ . The reason for this is that the rapid initial exponential growth shown in figure 13 (at $t'_0 \approx 300$ , $\sigma _{DNS}=0.0.005320$ and $\sigma _{LST}=0.005439$ ) is due to the localised mode. However, the spread mode that is also unstable keeps growing with a growth rate much lower than that of the localised mode. Nevertheless, a steady saturated state is obtained that again does not present signs of subsequent oscillations.

Table 4 summarises the maximum reverse flow for spanwise-homogeneous base flow ( $u_{0,rev}$ ), spanwise averaged saturated base flow ( $\bar {u}_{2D,rev}$ ) and three-dimensional saturated base flow ( $u_{3D,rev}$ ), for the three cases considered here. The difference between $u_{0,rev}$ and $\bar {u}_{2D,rev}$ is an indicator of the magnitude of the mean flow distortion due to spanwise modulation of the initially two-dimensional base flow. In all three cases, $u_{3D,rev}$ is larger and $\bar {u}_{2D,rev}$ is smaller than $u_{0,rev}$ . However, for case A with $\beta =1.2$ , where only the localised unstable mode is present, the bubble experiences a weak mean flow distortion as the difference between peak reverse flows is very small. Case G, despite the comparatively lower growth rate of its eigenmode, presents the largest peak mean flow distortion.

4.2.2. Base flow distortion induced by the three-dimensional instability

Figure 14. Contours of the base flow distortion $u_{bd}$ for case G ( $\beta =0.4$ ) at four different times. The streamwise planes are plotted at the spanwise location with a peak negative base flow distortion value. The red solid and dashed lines show the location of the dividing streamline for the initial two-dimensional base flow and instantaneous flow, respectively. The dashed black line shows the location of the streamwise inflection point. The thin grey dashed line shows the location of the reverse flow streamline. The colour bars are symmetric and the negative limit at each panel shows the peak negative base flow distortion value at the corresponding time. The figure is not drawn to scale.

Figure 15. Spatial evolution of amplitudes of individual spanwise harmonics modes ( $A_{\hat {u}}$ ) for case G at four time instances shown in figure 14. The separation and reattachment points are marked by vertical lines with labels $x_s$ and $x_r$ , respectively. The amplitude of the fundamental mode ( $1\beta$ ) and first five harmonics ( $2\beta -6\beta$ ) are shown in the figure.

The nonlinear simulations discussed in the previous section recover steady saturated flows, either because of the absence of spontaneously growing oscillations or because SFD is used to damp the oscillations. In consequence, the simulated disturbance field $\boldsymbol{u}'(x,y,z,t)$ corresponds to the nonlinear distortion of the two-dimensional base flow $\boldsymbol{u}_0$ induced by the three-dimensional instability. To remark this and distinguish it from the disturbance fields appearing in the simulations in later sections, the base flow distortion will be referred to as $\boldsymbol{u}_{bd}(x,y,z,t)= \boldsymbol{u}_{3D}(x,y,z,t) - \boldsymbol{u}_0(x,y)$ . Note that the base flow distortion evolves in time until the final saturated state is reached at time $t'=t'_{sat}$ , which is different for each case. Figures 14, 16 and 17 show the evolution of the base flow distortion for the cases in table 4. The streamwise velocity field is shown at the spanwise location where the base flow distortion at $t'=t'_{sat}$ has its peak negative value. In the plots, the dividing streamline for $u_0$ , the section of the dividing streamsurface for $u'$ and time $t'$ , the in-plane streamwise inflection points ( $\partial ^2 u / \partial y^2=0)$ and the section of the $u'= 0$ surface are also shown. Table 4 summarises the peak negative base flow distortion at the saturated state for the three cases considered here.

1. (i) Case G, $\beta =0.4$ (figure 14). At $t'=3306$ , the perturbations are still experiencing exponential growth (figure 13 b) and the base flow distortion is very similar to the global mode as expected. At $t'=4746$ , the nonlinear interactions start to become visible as the change in slope of the blue curve in figure 13(b) shows and, as a result, the bubble starts to deform. At $t'=6186$ , nonlinear interactions have further distorted the bubble and moved the location of the peak negative $u_{bd}$ away from the wall. Finally, at saturation ( $t'_{sat} \approx 55\,000$ ) the bubble is distorted completely and nonlinear interactions around the reattachment point introduce flow deformations that extend downstream of the bubble. At saturation, the line of external inflection points intersects the dividing streamline at $x\approx 334$ . Avanci et al. (Reference Avanci, Rodríguez and Alves2019) showed that, under the parallel flow assumption, this is a necessary condition for absolute instability of Kelvin–Helmholtz instability and self-excited oscillations. Note that downstream of the reattachment point, based on the shape of the global mode, the base flow distortion should be positive in the plane where the base flow distortion inside the bubble is negative. However, from the figure it can be seen that the bubble has a strong distortion in the vicinity of the reattachment point, and negative base flow distortion appears even further downstream. Thus, it can be conjectured that the strong deformation of the bubble, and also the negative base flow distortion downstream of the reattachment point, is mainly due to nonlinearities resulting from the growth of the spread mode.

   Further evidence is shown in figure 15, which illustrates the nonlinear growth of harmonics of the wavenumber corresponding to the unstable eigenmode. First, the spanwise Fourier decomposition of the streamwise perturbation velocity is obtained, i.e.

   (4.2) \begin{equation} u'=\sum _{n=0}^{N}\hat {u} \exp {(in \beta z)+\textrm{c.c.},} \end{equation}
   where $\beta$ is the fundamental spanwise wavenumber of the simulation domain. The amplitude of the individual spanwise harmonic mode is then defined at each streamwise coordinate $x$ as
   (4.3) \begin{equation} A_{\hat {u}}(x)=\max \limits _{y} |\hat {u}(x,y)|. \end{equation}
   Figure 15 shows the amplitudes of the fundamental ( $\beta$ ) and first five harmonics ( $2\beta -6\beta$ ) for the same time instances as in figure 14. While the amplitudes of the harmonics at the earlier times are nearly negligible, they grow over time showing a progressive transfer of energy from lower to higher harmonics and their eventual saturation at comparable amplitudes. Interestingly, the harmonics present their peak amplitudes in the vicinity of the reattachment point, explaining the localised changes during saturation observed in figure 14.

2. (ii) Case A, $\beta =1.2$ (figure 16). At $t'=500$ , the perturbations are still in their linear regime. Note that, as figure 12 also shows, the components of the global eigenmodes, unlike for the spread mode, vanish along the shear layer, except in a small region in the vicinity of the location where the global mode is maximum. For this reason, even at saturation, the bubble does not experience noticeable modulation and expansion except for a very weak distortion around $x \approx 360$ where the localised mode is present, and also upstream due to nonlinear interactions. This also explains the weak mean flow distortion shown in table 4 for this case. An important feature of the saturated base flow when only the localised mode is present happens in the vicinity of the reattachment point, where the $u$ component of the global eigenmode changes sign (figure 12). This is consistent with the base flow distortion observed at all times. As a result, in the vicinity of reattachment for this case, nonlinear interactions are in favour of increasing the streamwise velocity perturbations and cancelling the negative base flow distortion. Thus, negative-velocity base flow distortion does not happen outside of the base flow’s recirculating region and the bubble does not expand. Unlike for case G, the dividing streamline is practically unaltered.

   

   Figure 16. Contours of base flow distortion $u_{bd}$ for case A ( $\beta =1.2$ ) at four different times. The annotations are the same as for figure 14.

   

   Figure 17. Contours of base flow distortion $u_{bd}$ for case A ( $\beta =0.7$ ) at four different time instances. The annotations are the same as for figure 14.

3. (iii) Case A, $\beta =0.7$ (figure 17). In this case, two unstable eigenmodes coexist. Initially (e.g. $t'=300$ ), only the growth of the localised mode can be observed. No noticeable change in the dividing streamline is visible for this time. At $t'=600$ , nonlinear interactions start to become relevant as the dividing streamline is deformed in the vicinity of the localised mode. Until this time, no sign of growth of the spread mode is visible. However, after $t' \approx 2000$ , the spread mode becomes visible as shown by the base flow distortion at $t'=2640$ ; the base flow distortion resembles $u_{bd}(t'=6186)$ for case G. The interaction of the two unstable modes shows itself near the reattachment point. It was shown for case G that the nonlinearities following the growth of the spread mode distorts the bubble, especially near the reattachment point, pushing the peak reversed flow region away from the wall. On the other hand, in the vicinity of the reattachment point, nonlinearities following the localised mode intensify the positive base flow distortion and prevent the bubble from getting distorted near the reattachment point. Here, the interaction of these two modes is such that the bubble deforms slightly at the reattachment point, but still the distortion is much weaker than in case G. One explanation can be that the growth rate of the localised mode is almost one order of magnitude higher than the growth of the spread mode and, in consequence, the nonlinearities due to rapid growth of the localised mode prevent the slower-growing spread mode from distorting the flow effectively near the reattachment point. In summary, the base flow distortion in this case results from the combined effects of spread and localised modes: the rapid growth of the localised mode suppresses bubble deformation near the reattachment point, while the slower growth of the spread mode causes the bubble to expand upstream of the reattachment point.

Figure 18 shows the base flow distortion in the $z$ – $y$ cross-stream plane at the streamwise location of the peak negative base flow distortion for three cases, together with the sections of the dividing streamsurface, the reverse flow surface and the local surface of inflection points, following the same convention as in figures 14, 16 and 17. The spanwise modulation of the bubbles is clear from this plot, showing that case G experiences the most intense spanwise modulation, while case A with $\beta =1.2$ is only slightly modulated. Interestingly, due to strong nonlinear interactions for case G, the reverse flow region becomes very localised in the spanwise direction, which explains the large value of negative base flow distortion for this case compared with other cases, as given in the last column of table 4.

Figure 18. Contours of the base flow distortion in the $z$ – $y$ plane for $(a)$ case G ( $\beta =0.4$ ) at $x\approx 378$ , $(b)$ case A ( $\beta =1.2$ ) at $x\approx 363$ and ( $c$ ) case A ( $\beta =0.7$ ) at $x\approx 364$ . The annotations are the same as for figure 14. The thin grey solid line shows the location of the reverse flow streamline for the initially two-dimensional base flow.

5.1. Secondary global stability analysis

Figure 19. Secondary eigenvalue spectrum corresponding to the saturated flow for case G with $\beta =0.4$ . The inset shows the spectrum close to the origin.

Figure 20. Real part of the $u$ component of the most unstable eigenmode corresponding to the saturated $\boldsymbol{u}_{3D}$ for case G with $\beta =0.4$ . The yellow isosurface corresponds to $u_{bd}=-0.25$ . The red and blue colours show the positive and negative values, respectively.

Figure 21. Same as figure 20 for the second most unstable eigenmode corresponding to the saturated $\boldsymbol{u}_{3D}$ for case G with $\beta =0.4$ .

A three-dimensional global stability analysis of the saturated base flow $\boldsymbol{u}_{3D}$ for case G, which is obtained using SFD, is performed using the implicitly restarted Arnoldi method implemented in Nek5000. A Krylov subspace size of 100 is used to ensure convergence of at least the first 20 eigenvalues. The eigenspectra, shown in figure 19, contains two unstable modes with ( $\omega, \sigma )=(0.1479,0.0065)$ and ( $\omega, \sigma )=(0.1555,0.0198)$ . In the proximity of $\omega =0$ (the zoomed inset in the figure), a family of stable low-frequency modes appears. The shape and role of these modes in the dynamics is explained in § 6.1. The spatial structure of the $u$ component of the most unstable and second unstable modes are shown in figures 20 and 21, respectively. In these figures, the yellow isosurface corresponds to $u_{bd}=-0.25$ , i.e. $86 \,\%$ of the peak negative base flow distortion $u_{bd,min}$ . Both unstable modes are varicose (symmetric) on the spanwise direction and are rooted around the peak negative base flow distortion. Note that both the peak reverse flow and base flow distortion $u_{bd,min}$ occur at the same spanwise plane; however, the location of $u_{bd,min}$ , $(x,y)_{bd}=(378.4,0.97)$ , is located slightly downstream and above the location of the peak negative reverse flow ( $(x,y)_{rev}=(369.8,0.62)$ ). The second unstable mode has a slightly longer wavelength compared with the most unstable one, and it has a larger downstream spatial extension. Using a weakly non-parallel Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) analysis, Rodríguez et al. (Reference Rodríguez, Gennaro and Souza2021a ) approximated the properties of the most unstable eigenmode responsible for self-excited secondary instability in a pressure-gradient-induced LSB over a flat plate. The spatial structure of the eigenmode found in that work looks qualitatively similar to both eigenmodes found here using global stability analysis. Furthermore, the varicose eigenmodes here look similar to the eigenmodes found in the wake region of a cylindrical isolated roughness element with a diameter-to-height aspect ratio larger than one (Loiseau et al. Reference Loiseau, Robinet, Cherubini and Leriche2014), and also to those recovered for a cuboid isolated roughness element when the height of the roughness is large compared with the boundary layer displacement thickness (Ma & Mahesh Reference Ma and Mahesh2022). It is anticipated here that these eigenmodes will give rise to self-excited oscillations in the rear part of the bubble, in the vicinity of $u_{bd,min}$ , and will cause a rapid transition to turbulence. The transition process is investigated in detail in § 6.1.

5.2. Impulse response analysis

To investigate the linear response of the saturated flow $\boldsymbol{u}_{3D}$ to perturbations, a linear simulation is performed in which a spatially localised initial disturbance is introduced at $x\approx 370$ , close to the location of the peak reversed flow. The spatial structure of the initial disturbance is described in Appendix B.

The solid black line in figure 22 shows the evolution of the disturbance kinetic energy $E$ and the red dashed line shows the growth rate of the most unstable eigenmode calculated using global stability analysis confirming that, after an initial transient growth, the linear disturbance is governed by the most unstable eigenmode. Due to the large gap between the growth rates of the two unstable modes, the most unstable mode becomes dominant very quickly. As a result, a beating phenomena at the frequency difference between those of the two unstable modes, as that recovered by, for example, Ehrenstein & Gallaire (Reference Ehrenstein and Gallaire2008), does not happen here. Instead, an exponential disturbance growth in which the disturbance field (figure 20) and the oscillation frequency closely follow the prediction of the global stability analysis is recovered and it continues until nonlinearities set in.

Figure 22. The time evolution of the perturbation energy (solid black line) as a result of an impulse perturbation field imposed at $t=0$ , considering the saturated flow $\boldsymbol{u}_{3D}$ for case G with $\beta =0.4$ as the base flow. The red dashed line shows the growth rate of the most unstable eigenmode in figure 19.

6.1. Case G, $\beta =0.4$

Figure 23. Spatio-temporal diagram of spanwise perturbation velocity for case G with $\beta =0.4$ , along a curve parallel to the wavy wall that passes through the location of $u_{bd,min}$ . The dashed black lines show the location of the separation and reattachment point, and the red dashed line shows the location of peak negative base flow distortion. Note that the perturbations are calculated based on the saturated flow $\boldsymbol{u}_{3D}$ .

Figure 24. Top panel: the time evolution of the perturbation kinetic energy (calculated with respect to the saturated $\boldsymbol{u}_{3D}$ ) for case G with $\beta =0.4$ . The inset zooms in on region ‘S’ corresponding to $11\,950\lt t'\lt 12\,150$ . Bottom panel: maximum $y$ coordinate of the dividing streamline in the streamwise plane with peak negative base flow distortion. The horizontal dashed line indicates the maximum $y$ coordinate of the dividing streamline in the saturated state. The shaded areas indicate the P, S and T temporal phases for the second turbulent burst.

The simulation shown in figure 13 is continued without applying SFD after the appearance of self-excited oscillations at $t'\approx 9000$ . Note that the oscillations start to appear before the primary instability fully saturates and the peak negative base flow distortion ( $u'_{bd,min}\approx 0.27$ ) is slightly lower than its value at the fully saturated state that is used to perform the secondary stability analysis in § 5, namely $u'_{bd,min}\approx 0.29$ . Figure 23 shows the $x{-}t$ diagram for spanwise perturbation velocities for a line ( $x \in [275,485]$ ) parallel and close to the wavy wall in the streamwise plane that passes through the location of the peak negative base flow distortion corresponding to the saturated state. The perturbation velocity here is calculated using the saturated flow $\boldsymbol{u}_{3D}$ as reference, not the two-dimensional initial base flow $\boldsymbol{u}_0$ . The spatio-temporal diagram shows two specific features of this flow. First, after the appearance of oscillations at $t' \approx 9000$ , breakdown to turbulence happens, starting around $x \approx 378$ . However, the flow undergoes a relaminarisation after that. This turbulent burst-relaminarisation cycle repeats itself subsequently, with an angular repetition frequency $\omega \approx 0.0048-0.0057$ , estimated from figure 23. Second, after each turbulent burst, spanwise perturbations travel upstream within the bubble. In the following, each of these observations will be explained in detail.

6.1.1. Turbulent burst-relaminarisation cycle

The top panel in figure 24 shows the kinetic energy of the perturbations integrated in the domain for $280\lt x\lt 490$ (normalised by the energy at $t'=4000$ ). The bottom panel shows the maximum $y$ coordinate of the dividing streamline of the bubble in the spanwise plane with peak negative base flow distortion, which is an indicator of the instantaneous height of the bubble. Starting from $t'=4000$ , the height of the bubble increases slowly until $t'\approx 4900$ , when the nonlinear interactions become relevant and induce a strong spanwise modulation of the bubble, as explained in § 4.2.2. Then, the height of the bubble keeps increasing, which results in a strong peak negative base flow distortion and triggers a self-excited secondary instability. Soon after the appearance of oscillations (figure 13), for $9300\lt t'\lt 9600$ , the secondary instabilities grow exponentially in the domain that cause the rapid increase of perturbation kinetic energy. Once the flow oscillations arising from the secondary instability reach nonlinear amplitudes, for $9600\lt t'\lt 10\,000$ , a rapid transition to turbulence happens (which is also observable in figure 23). During this lapse of time, the vortical structures formed inside the bubble pull the fluid out of the bubble and cause the height of the bubble to decrease rapidly. Eventually, at $t' \approx 10\,000$ , the vortical structures leave the bubble and are advected downstream until they exit the computational domain. This is also visible in figure 23 where, at $t'\approx 10\,000$ , the leading edge of the spanwise perturbations is leaving the domain, causing a rapid decrease in perturbation energy. As a result, the bubble relaminarises and shrinks until it reaches an intermediate state with a minimum height ( $t'=10\,200$ ). Following from this initial turbulent burst, a cycle of spanwise deformation, secondary instability and turbulent burst is repeated in a quasi-periodic manner for the rest of the simulation.

The shaded areas in figure 24 show three time lapses named P, S and T, which denote the main phases of the subsequent turbulent bursts. The regions are associated with primary (P) and secondary (S) instabilities, and the transition (T) region. After the first turbulent burst leaves the region adjacent to the wall waviness, the bubble reaches an intermediate state, at $t'\approx 10\,200$ , between its initial two-dimensional state and the saturated state where no secondary instability is active anymore; however, the primary instability is still globally unstable and drives the flow towards its saturated state again. This is clear from region P in figure 24, where the kinetic energy growth rate changes and the bubble height starts to increase until it reaches its maximum again. After that, in phase S, the negative base flow distortion is strong enough to trigger the secondary instabilities and the ensuing turbulent burst reoccurs in phase T.

Recently, in an experiment of a pressure-gradient-induced LSB on a flat plate (Aniffa & Mandal Reference Aniffa and Mandal2023), a similar phenomenon was reported: a low-frequency oscillation of the shear layer and vertical motion of the inflection points that is accompanied by an intermittent vortex shedding where the bubble alternates between a non-vortex-shedding state and a shedding state. Here, it is also observed at $t'\approx 8400$ that the external streamwise inflection point crosses below the dividing streamsurface of the recirculation bubble. However, no sign of oscillations was observed until $t'\approx 9000$ , when the height of the bubble gets close to its maximum and the peak negative base flow distortion becomes strong enough to trigger the secondary instabilities. This is in agreement with the necessary condition for absolute instability of inflectional instability proposed by Avanci et al. (Reference Avanci, Rodríguez and Alves2019) and its three-dimensional analog seen in Rodríguez et al. (Reference Rodríguez, Gennaro and Souza2021a ).

Figure 25. Instantaneous $\lambda _2$ visualizations (coloured by streamwise velocity) of the vortical structures in the flow for case G with $\beta = 0.4$ . The red (positive) and blue (negative) isosurfaces show the spatial structure of the second unstable global mode as shown in figure 21. The energy growth shown in the inset of figure 24 corresponds to the time interval shown in this figure. Note that different values of $\lambda _2$ are used in the different plots for better visualization.

Figure 25 shows the vortical structures using $\lambda _2$ visualizations (Jeong & Hussain Reference Jeong and Hussain1995) coloured by streamwise velocity for six different time instants during $11\,954\lt t'\lt 12\,110$ , corresponding to phase S before the third peak of disturbance energy in figure 24. The $u$ component of the second unstable secondary instability (same as figure 21) is also shown. The inset in the top panel of figure 24 shows the growth of the disturbance energy for the same time span. Figures 24 and 25 show that, after the appearance of oscillations, a train of hairpin vortices starts to develop as a nonlinear consequence of the global eigenmode. The train of hairpin vortices here is similar to that reported downstream of a cylindrical roughness element by Loiseau et al. (Reference Loiseau, Robinet, Cherubini and Leriche2014), where the authors concluded that a varicose unstable global mode contributes to nonlinear generation of the hairpin vortices. From $t'\approx 12\,014$ the disturbance energy increases exponentially until saturation, which occurs at $t' \approx 12\,110$ due to nonlinear interactions. The growth rate of the secondary instabilities estimated by using the disturbance energy from DNS for $12\,014\lt t'\lt 12\,074$ (the slope of the red dashed line in the inset of figure 24) is $\sigma _{DNS}\approx 0.0161$ , which is slightly lower than the value $\sigma = 0.0198$ found using global stability analysis. The growth of secondary instabilities is also calculated for all other S phases, which are found to be in the range $\sigma _{DNS} = 0.014-0.016$ . The early presence of hairpin vortices, even at $t'=11\,954$ , shows the importance of nonlinear effects as soon as self-excited secondary instabilities are triggered. This can contribute to the lower growth rates approximated from DNS compared with those calculated from global stability analysis. Moreover, oscillations appear at a peak negative base flow distortion that is slightly lower than its value at the fully saturated state found with SFD. Also, as the bottom panel in figure 24 shows, the state of the bubble at the end of each different P phase is different: for each S phase, secondary instabilities are triggered for different values of peak negative base flow distortion, which might explain the range of growth rates found from DNS.

To assure the quality of the results, simulations have been performed for a finer and bigger domain, which the corresponding results are given in Appendix A.3.

6.1.2. Upstream propagation of perturbations

In figure 23 it was shown that after the occurrence of a turbulent burst, spanwise perturbations propagate upstream within the bubble. The fastest speed of upstream propagation of perturbations is estimated from the $x{-}t$ diagram to be in the range [−0.051, −0.048]; however, as they get closer to the separation point, their propagation speed decreases. Figure 26 shows the spanwise velocity perturbations (with respect to the saturated flow $\boldsymbol{u}_{3D}$ ) in the spanwise plane with peak negative base flow distortion ( $u_{bd,max}$ ) while figure 27 shows the same flow quantity for a plane parallel to the wall that contains $u_{bd,max}$ . The panels in the first column show the instantaneous perturbation field for three time instances before (top panel) and during (middle panel) the first burst, and slightly after the leading edge of perturbations has left the bubble (bottom panel). The panels in the right column show the same thing but for the second turbulent burst ( $10\,625\lt t'\lt 11\,550$ ).

Figure 26. Spanwise perturbation velocity for six different times plotted at the spanwise plane with peak negative base flow distortion. The red dashed line shows the location of instantaneous dividing streamline. The colour bar is saturated to visualise the flow structures. The figure is not drawn to scale.

Figure 27. Spanwise perturbation velocity for six different times plotted for a plane parallel to the wall that passes through the peak negative base flow distortion. The colour bar is saturated to visualise the flow structures. The figure is not drawn to scale.

The first panel, at $t'=9401$ , shows the evolution of self-excited secondary instabilities while still they are growing exponentially (figure 24). At this time instance the bubble height is close to its maximum. The shape of the perturbations is the same as the global eigenmode obtained by the secondary instability analysis. Note that since the $u$ component is global mode being symmetric (figure 20), the spanwise component must be antisymmetric, which is visible in figure 27 at $t'=9401$ . Shortly after the secondary instability saturates and nonlinear interactions become more intense, the disturbance energy reaches a saturated state and a turbulent burst is formed in the aft part of the recirculation bubble, as can be seen in the second panel at $t'=9665$ . Close to the reattachment point ( $x\approx 370$ ), spanwise perturbations start to create a boomerang-shape structure, which becomes more evident as the leading edge of the perturbation leaves the bubble and the height of the bubble is close to its minimum. The boomerang-like structures then start to propagate upstream within the bubble (third panel at $t'=10\,025$ ). Note that after the boomerang-shape structures are formed in the flow, the symmetry of the flow structure is broken as can be seen in figure 27 at $t'=10\,025$ . Afterwards, the bubble relaminarises and gradually expands in the wall-normal direction (phase P in figure 24), while the spanwise perturbations that emerged from the vicinity of reattachment in the previous turbulent burst keep propagating upstream in the bubble. Eventually, self-excited oscillations reappear as can be seen in the first panel of the right column at $t'=10\,625$ . The formation of the boomerang-like spanwise perturbations is visible from the second panel of the right column at $t'=10\,949$ . Since the propagation speed of the perturbations decreases as they get closer to the separation point, the new perturbations originated in the second turbulent burst reach and interact with those from the first burst and create more complex structures, as can be seen from the third panel at $t'=11\,549$ . The same cycle is repeated for later times as shown by figure 23.

The upstream-propagating boomerang-like structures resemble a family of stable, low-frequency eigenmodes recovered by the global stability analyses. Figure 28( $a$ – $c$ ) show the spanwise velocity component of three global eigenmodes of the saturated flow $\boldsymbol{u}_{3D}$ , marked by Ms1, Ms2 and Ms3 in the part of the spectrum shown in panel ( $e$ ). Panel ( $d$ ) shows the three-dimensional shape of the Ms3 mode. The structures are very similar to the shape of the perturbations shown in figures 26 and 27, although the streamwise wavelength observed in the DNS is not the same as the wavelength of the individual eigenmodes shown in figure 28.

Figure 28. ( $a$ – $c$ ) Spanwise velocity component (real part) of global eigenmodes labelled by Ms1, Ms2 and Ms3 in panel ( $e$ ), at the spanwise plane with peak negative base flow distortion. ( $d$ ) Three-dimensional view of the spanwise velocity field of the Ms3 mode. The red and blue colours show positive and negative values, respectively. ( $e$ ) Eigenspectrum corresponding to global stability analysis of the saturated flow $\boldsymbol{u}_{3D}$ (same as the inset in figure 19). Panels ( $a$ - $c$ ) are not drawn to scale.

The same family of eigenmodes is found in the global stability analysis of the two-dimensional base flow $\boldsymbol{u}_0$ (figure 36). However, the shape of the eigenmodes in the spanwise direction is modified on account of the absence of the spanwise distortion that defines the saturated flow. An adjoint sensitivity analysis (analogous to that by Giannetti & Luchini Reference Giannetti and Luchini2007) has been done for the same family of eigenmodes corresponding to the two-dimensional base flow; the region slightly upstream of the reattachment point, located in the vicinity of the peak negative base flow distortion was found to be the most sensitive region. It should be noted that the adjoint sensitivity analysis is not done for the saturated state due to its very high computational cost. However, since the low-frequency eigenmodes for both two-dimensional and saturated flows share many similarities, it is conjectured that the sensitivity characteristics of them are also similar.

These observations suggest that the upstream-propagating spanwise perturbations are associated with the low-frequency stable eigenmodes, which get excited through the interactions with the self-excited secondary instabilities upstream of the reattachment point. For further analysis, a three-dimensional resolvent analysis could be done, but this is out of the scope of the present study. After their nonlinear excitation, these eigenmodes can describe a non-modal behaviour or interact nonlinearly, explaining why the wavelength of the perturbations from DNS is not constant and does not match the individual wavelength of the global stability eigenmodes.

6.2. Case E, $\beta =0.6$

The global stability analyses presented in § 4 show that the two-dimensional base flow for case E sustains both two-dimensional oscillatory and three-dimensional steady unstable eigenmodes. Nonlinear DNS simulations with $\beta =0.6$ $(L_z=2\pi /0.6)$ are performed to investigate the interactions between the unstable modes and their subsequent evolution. Details of the DNS domain size and grid are given in Appendix A.3.

Figure 29. Maximum $y$ coordinate of dividing streamline in the streamwise plane with peak negative base flow distortion (left panel). Spatio-temporal diagram for the $u$ component of perturbation velocity for case E, $\beta =0.6$ , for a line parallel to and with distance $\Delta y \approx 1.45$ from the wall at the same plane (right panel). Perturbations are calculated based on the two-dimensional base flow.

Figure 29 shows the $x$ – $t$ diagram for streamwise perturbation velocity (calculated with respect to the initial two-dimensional base flow $\boldsymbol{u}_0$ ) at the spanwise plane with peak negative base flow distortion, for a plane parallel to the wavy wall and separated by a distance $\Delta y \approx 1.45$ from it. An initial transient has been eliminated from the figure. Self-excited vortex shedding, as predicted by global stability analysis (§ 4), develops spontaneously and is identified by the dark blue spots in the figure.

The growth of three-dimensional instability can be inferred from $t' \approx 2000$ onwards, when the negative base flow distortion increases upstream of the reattachment point. As shown in § 4.2, nonlinear interactions originated by the growth of the spread mode give rise to a spanwise modulation of the bubble. Thus, the total streamwise velocity decreases in the spanwise locations where the base flow distortion inside the bubble is negative, which simultaneously increases the $y$ coordinate of the dividing streamline of the bubble. As a result of the nonlinear distortion of the underlying base flow, the self-excited oscillations originating from the initially two-dimensional unstable eigenmodes experience a spanwise shear and different advection speeds at different spanwise locations. This leads to the formation of $\Lambda$ vortices with a spanwise wavelength equal to the fundamental spanwise wavelength of the dominant three-dimensional global mode, where the legs of the vortex lie at the spanwise planes with peak negative perturbation inside the bubble. The $\Lambda$ vortices turn into hairpin vortices as they travel downstream and their breakdown triggers a rapid transition to turbulence. The formation of these vortices is shown in figure 30.

Figure 30. Instantaneous $\lambda _2$ visualizations of the vortical structures in the flow (green isosurfaces) at $t'=2570$ . The blue and red colour show the negative ( $u'=-0.05$ ) and positive ( $u'=0.03$ ) isosurfaces of streamwise perturbation velocity inside the bubble. The horizontal plane shows the streamwise perturbation velocity for a plane at $y=0.2$ . The figure is shown only for a reduced part of the domain ( $300\lt x\lt 485$ ). For $\lambda _2$ visualizations, for $x\lt 435$ , isosurfaces of $\lambda _2=-0.0006$ and for $x\gt 435$ , isosurfaces corresponding to $\lambda _2=-0.005$ are plotted. Two spanwise periods of the simulation domain are shown.

6.2.1. Unsteady flow characteristics

To analyse the complex three-dimensional vortex shedding and ensuing laminar–turbulent transition, spectral proper orthogonal decomposition (SPOD) is used (Schmidt & Colonius Reference Schmidt and Colonius2020). The cross-spectral density employed is based on the kinetic energy. Here, the data from $t'\approx 2400$ to the end of the simulation (figure 29), which consist of 11 500 three-dimensional snapshots, are divided into blocks of 1150 snapshots with $50 \,\%$ overlap (19 blocks in total). This leads to a frequency-domain discretization with $\Delta \omega =0.0091$ and Nyquist frequency $\omega _{max} = 5.23$ . Figure 31 shows the SPOD energy of the leading mode (the black solid line), while the red dashed line in the figure shows the sum of the energy of all SPOD modes for each frequency.

Figure 31. Leading SPOD mode energy in the domain. The red dashed line represents the sum of the energy of all SPOD modes for each frequency.

The SPOD spectrum presents several clear peaks. Vertical solid lines denote specific low and high frequencies, while dashed lines represent higher harmonics of each frequency indicated by the same colour. Starting from high-frequency peaks, the black solid lines specify two frequencies named $\omega _{\textit{KH}1,0}=0.1278$ and $\omega _{\textit{KH}2,0}=0.1661$ that are very close to the frequency of two unstable two-dimensional global modes for case E ( $\omega _{\textit{KH}1,LST}=0.1280$ and $\omega _{\textit{KH}2,LST}=0.1650$ ), as shown in figure 8. Note the higher peak for $\omega _{\textit{KH}1,0}$ than for $\omega _{\textit{KH}2,0}$ that is in agreement with the higher growth rate of $\omega _{\textit{KH}1,LST}$ compared with $\omega _{\textit{KH}2,LST}$ , as figure 8 shows. The black dashed line specifies the first harmonic $(\omega _{\textit{KH}1,1})$ of $\omega _{\textit{KH}1,0}$ . The frequency peak at $\omega =0.0894$ is also very close to the frequency of a globally stable mode with $\omega =0.0882$ in the eigenmode spectrum shown in figure 8. Figure 32 shows the real part of the leading SPOD mode at $\omega =0.1278$ for two spanwise periods and for $x\lt 485$ . This mode consists of spanwise-modulated rollers appearing just downstream of the reattachment point and extending further downstream. Its spatial shape resembles the hairpin vortices from figure 30, where the heads of the hairpin vortices are again visible at $z\approx \pm 5$ and $z\approx 15$ . Moreover, the SPOD mode exhibits other sets of hairpin-like structures that appear farther downstream (from $x\approx 470$ ) and present a different spanwise phase compared with the fundamental vortices: i.e. their heads are in the spanwise plane with peak negative base flow distortion inside the bubble ( $z\approx 0$ ). These secondary hairpin vortices are also observed in DNS and can be seen in the magnified view shown in figure 30. Note that the SPOD mode corresponds to $\omega _{\textit{KH}2,0}=0.1661$ (not shown here) is very similar to the one for $\omega _{\textit{KH}1,0}=0.1278$ .

Figure 32. Real part of the $u$ component of the leading SPOD mode for $\omega =0.1278$ . Red and blue colours correspond to ${u}'=\pm 0.0015$ and $u'=-0.0015$ , respectively. Two spanwise periods of the simulation are shown for visualization.

The low-frequency oscillations, on the other hand, cannot be directly related to two-dimensional self-excited mechanisms as there is no low-frequency unstable mode in the global spectrum for this case. In the absence of external forcing, as is the case here, they must be due either to interactions between high-frequency oscillations or to the nonlinearities resulting from the growth of three-dimensional instabilities and base flow distortion. The blue solid line specifies a low-frequency peak with $\omega =0.0383$ . Interestingly, this frequency is equal to the difference between the frequency of $\omega _{\textit{KH}1,0}$ and $\omega _{\textit{KH}2,0}$ , i.e. the beating frequency due to the interaction of $\omega _{\textit{KH}1,0}$ and $\omega _{\textit{KH}2,0}$ . The first harmonic of this frequency is also specified by the blue dashed line in the figure. The high energy content in the very low-frequency range ( $\omega \lt 0.03$ ) suggests that there might still be low-frequency dynamics in the SPOD spectrum that requires much longer time series to be resolved. The low-frequency, high energy content in the separation bubbles is attributed to the so-called breathing or flapping of the bubble (Zaman et al. Reference Zaman, Mckinzie and Rumsey1989; Sigurdson Reference Sigurdson1995), i.e. a quasi-periodic vertical movement of the separated shear layer from the wall. This very low-frequency vertical movement of the bubble is shown in the left panel of figure 29. Furthermore, when the height of the bubble increases, a strong vortex is shed from the bubble that results in the formation of a $\Lambda$ vortex (with head located at $z\approx 0$ ) just downstream of the reattachment point. This corresponds to the dark blue spots downstream of the reattachment point coincident with the increase of the bubble height (marked by horizontal solid lines in figure 29). Figure 33 shows $\lambda _2$ visualization of the vortical structures at $t'=6000$ . The hairpin vortices downstream of the reattachment point are highlighted in the magnified view of the figure.

Figure 33. Instantaneous $\lambda _2$ visualizations of the vortical structures in the flow (green isosurfaces) at $t'=6000$ . The blue and red colours show the negative ( $u'=-0.05$ ) and positive ( $u'=0.055$ ) isosurfaces of streamwise perturbation velocity inside the bubble. The horizontal plane shows the streamwise perturbation velocity for a plane at $y=0.2$ . The figure is shown only for a part of the domain ( $300\lt x\lt 485$ ). For $\lambda _2$ visualizations, for $x\lt 435$ , isosurfaces corresponding to $\lambda _2=-0.002$ are plotted and, for $x\gt 435$ , isosurfaces corresponding to $\lambda _2=-0.08$ are plotted. Note that two spanwise periods of the simulation are shown for visualization.