4.1. Numerical base-flow solutions for different APG magnitudes

Figure 2 shows the computed laminar base-flow solutions for a range of suction–blowing amplitudes. As the amplitude $v_0$ is gradually increased, the boundary layer becomes more prone to separation. At $v_0=0.125$ we reach incipient separation and at $v_0=0.15$ we obtain the first LSB solution with $u_{\it{rev}}=1.06\, \%$ (see figure 2 b). Figure 2(a) displays four boundary layers whose linear disturbances will be studied in § 5. Case A is the ZPG boundary layer reference, cases B and C are boundary layers under APG and case D shows a separated shear layer hosting a LSB with peak reversed flow $u_{\it{rev}}=2\, \%$ . The boundary layer displacement thickness, $\delta ^{*}$ , and the streamlines show the influence of the APG starting at $x\approx 0.6\times 10^5$ , where the boundary layer deviates from the ZPG reference. Additional analysis, including validation and comparisons with the broader LSB literature, in Appendix A indicates that our model LSB is representative of a ‘short’ bubble whose transition is typically driven by convective disturbances.

Figure 2. (a) Contours of streamwise velocity and streamlines from four base-flows computed at $v_0=0$ (A, ZPG boundary layer), $v_0=0.05$ (B, low APG boundary layer), $v_0=0.1$ (C, medium APG boundary layer) and $v_0=0.16$ (D, LSB with peak reversed flow $u_{\it{rev}}=2\, \%$ ). Displacement thickness from base-flow solution (black-dashed) and from Blasius solution (black-solid) superimposed. The LSB dividing streamline (yellow-solid) is also shown in D. (b) Change of $u_{\it{rev}}$ with varying $v_0$ .

4.2. Global stability

From our GSA calculations (see Appendix B), we observe that the LSB becomes globally unstable due to the emergence of an unstable 3-D, stationary (zero frequency) eigenmode for peak recirculation values beyond 6.60 %, in agreement with the existing literature reporting values in the range 7 %–8 % (Duck et al. Reference Duck, Ruban, Theofilis, Hein and Dallmann2000; Rodríguez et al. Reference Rodríguez, Gennaro and Juniper2013, Reference Rodriguez, Gennaro and Souza2021). We add an indicative neutral stability margin in figure 2 at 7 % to mark the threshold between the stable and unstable regimes. Further investigation on this mode (see Appendix B) reveals the features of a modal centrifugal mechanism, occurring at much lower $u_{\it{rev}}$ than absolutely unstable 2-D K-H waves, which were detected only above 15 % in flat plate (Rodríguez et al. Reference Rodríguez, Gennaro and Juniper2013) and bump (Ehrenstein & Gallaire Reference Ehrenstein and Gallaire2008) geometries. For any LSB base-flow computed within the $v_{0}$ range considered, we have not observed any globally unstable 2-D K-H waves. This is in accordance with literature since we have generated bubbles up to $u_{\it{rev}} \approx 12\, \%$ .

In § 5 we study the optimal linear instability mechanisms of the four base-flows (A to D) using the resolvent approach in § 2.2. Note that all of the four cases investigated further are globally stable, noise-amplifier-type flows.

5.1. Parametric study in the $(\beta ,\omega )$ domain

We perform a parametric study with RA by scanning a combination of frequencies $\omega = [0,50 ]\times 10^{-5}$ and spanwise wavenumbers $\beta = [0,250 ]\times 10^{-5}$ on cases A to D. In other words, we solve the optimisation problem (2.8) to extract the optimal linear forcing/response mode at each $(\beta ,\omega )$ pair to unravel frequencies of high receptivity of the boundary layer to external disturbances and the physics of the mechanisms involved. Maps of the linear gain are shown in figure 3 for cases A to D in the parametric space $(\beta ,\omega )$ . Two main classes of disturbances are identified at selective frequencies and spanwise wavenumbers:

1. (i) 3-D ( $\beta \neq 0$ ) travelling-wave modes at frequency $\omega$ (diamond markers in cases A and B and square markers in cases C and D);

2. (ii) 3-D ( $\beta \neq 0$ ) steady ( $\omega =0$ ) modes (squares in cases A and B and diamonds in cases C and D),

which are listed in table 2. We examine the physics of these disturbances in the remainder of this section.

Figure 3. Contours of ${\mathcal{G}}^2(\beta ,\omega )$ . Square ( $\square$ ) and diamond ( $\lozenge$ ) markers refer to the most and second most amplified $(\beta ,\omega )$ pair, respectively.

Table 2. Summary of linear amplification mechanisms for cases A to D. The symbols have the same meaning as in figure 3.

5.2. From Tollmien–Schlichting to K-H-dominated boundary layer instability

Depending on the magnitude of the APG, the receptivity of the boundary layer to external disturbances changes both quantitatively and qualitatively, meaning that the pressure gradient affects both the energy/amplification and the characteristic physics of the underlying mechanisms. For the former, we refer to the gain in figure 3 and for the latter we plot either the streamwise evolution of the square root of the perturbation kinetic energy, $\sqrt {E(x)}$ , in figure 4, where $E$ is

(5.1) \begin{equation} E(x) = \int _{0}^{\infty } ( |\hat {u}|^2 + |\hat {v}|^2 + |\hat {w}|^2 ) \textrm{d}y \end{equation}

and a planar view of the spatial forcing/response modes in figure 5.

Figure 4. Square root of the kinetic energy associated with the linear response modes marked in figure 3. (a) Unsteady ( $\omega \neq 0$ ) and (b) steady ( $\omega =0$ ) modes. For case D, grey-solid lineis separation; grey-dashed line is reattachment; shaded grey area is the separation region. Case D, panel (a), is down-scaled by one order of magnitude for clearer visualisation. The inset shows a close-up of cases A and B.

Figure 5. Optimal unsteady resolvent modes showing the change of instability from T-S dominated to K-H dominated with increasing APG. Contours of (a) ${Re} \{ \hat {f}_{x} \}$ and (b) ${Re} \{ \hat {u} \}$ . Black-dashed line, displacement thickness; grey-solid lines, base-flow streamwise velocity profiles. Red-solid line, inflection line; red shaded area, critical layer.

From the linear gain in figure 3, we observe that ZPG (A) and low APG (B) attached boundary layers amplify Tollmien–Schlichting (T-S) waves according to the Orr mechanism (Schmid & Henningson Reference Schmid and Henningson1992, Reference Schmid and Henningson2001; Rigas et al. Reference Rigas, Sipp and Colonius2021), which is more energetic for oblique rather than planar waves. We find $(\beta ,\omega )=(30,10)\times 10^{-5}$ for ZPG and $(\beta ,\omega )=(50,17.5)\times 10^{-5}$ for low APG (marked with diamonds in figure 3 a,b) to be the frequency-spanwise wavenumber pairs where this mechanism achieves maximal amplification. While the gain of the oblique T-S wave mechanism is not affected dramatically by the APG, the shape of the spatial growth rate is, as shown in figure 4(a). The inset indicates that in the ZPG boundary layer the energy of T-S waves increases monotonically with increasing $x$ ( $={Re}_x$ in our non-dimensionalisation), while in the presence of low APG there is a peak of energy around ${Re}_x \approx 1.2{-}1.3\times 10^5$ . In figure 5 (a i,b i) we show that the optimal forcing waves located in the first part of the pressure gradient region and within the critical layer excite T-S waves downstream. The spatial arrangement of the response, in agreement with the energy distribution in figure 4(a), indicates that the amplification of T-S waves correlates with the pressure gradient distribution. This is because the pressure gradient enhances the viscous shear in the boundary layer, which is ultimately responsible for the Orr mechanism (Schmid & Henningson Reference Schmid and Henningson2001). Despite all the effects of the pressure gradient on T-S waves discussed thus far, the low APG boundary layer is still reminiscent of the same convective instability mechanisms of the ZPG boundary layer.

Stronger APG levels (case C) lead to (i) boundary layer thickening and (ii) inflectional velocity profiles. The latter is determinant for important changes in the dominant instability mechanism of the boundary layer. The energy of the most amplified mode found at $(\beta ,\omega )=(30,22.5)\times 10^{-5}$ (square marker in figure 3 c, case C) increases by one order of magnitude with respect to the T-S mode of case B but its distribution in $x$ is qualitatively similar – see figure 4(a). To detect any differences between these modes we compare their spatial structures in figure 5. Figure 5(a ii,b ii) show that the optimal forcing waves cluster in a narrower region of the domain more upstream compared with the classic T-S forcing mode and decay upstream of the boundary layer thickening, moreover, the response mode exhibits roller-type structures downstream of boundary layer thickening closely aligned with the inflection line. We ascribe these features to the inviscid K-H instability mechanism, which emerges when the inflection point becomes unstable for sufficiently large APG (Sandham Reference Sandham2008; Diwan & Ramesh Reference Diwan and Ramesh2009; Marxen & Henningson Reference Marxen and Henningson2011; Zhang et al. Reference Zhang, Rowley, Wu, Meneveau and Mittal2019).

When the APG leads to flow separation (case D), the gain of the most amplified K-H mode at $(\beta ,\omega )=(20,25)\times 10^{-5}$ increases by two additional orders of magnitude due to the formation of a separated shear layer which is known to be more unstable than the attached counterpart (Yang Reference Yang2019). The optimal forcing waves live upstream of the separation bubble and feed the response amplification occurring in the aft part of the separated shear layer (figure 5 a iii,b iii), which is a critical site of shear layer instability (Sandham Reference Sandham2008; Jones et al. Reference Jones, Sandberg and Sandham2010; Simoni et al. Reference Simoni, Ubaldi and Zunino2014). The peak response energy is reached shortly downstream of the reattachment point (figure 4 a), where both velocity and pressure gradients are still significant. In case D it is even more evident than case C that the K-H rollers form two streets, one aligned with the inflection line and the other being very close to the wall. This is not observed in the shape of the T-S waves (case B). Because the most amplified K-H mode is 3-D, even though 2-D K-H waves exist in the gain spectrum of figure 3(d) (case D), the angle of the waves relative to the spanwise direction can be computed as $\tan ^{-1}(\beta /\alpha _{x})$ , where $\alpha _{x}$ is the streamwise wavenumber. We find that it is approximately $8^{\circ }$ , meaning that these waves are nearly 2-D. These topological features agree with the broad literature on K-H instability of separated shear layers (Rist & Maucher Reference Rist and Maucher2002; Diwan & Ramesh Reference Diwan and Ramesh2009; Simoni et al. Reference Simoni, Ubaldi and Zunino2013, Reference Simoni, Ubaldi and Zunino2014; Michelis et al. Reference Michelis, Kotsonis and Yarusevych2018a ,Reference Michelis, Kotsonis and Yarusevych b ; Zhang et al. Reference Zhang, Rowley, Wu, Meneveau and Mittal2019; Hosseinverdi & Fasel Reference Hosseinverdi and Fasel2019, Reference Hosseinverdi and Fasel2020).

In summary, the APG, upon modification of the boundary layer velocity profile, enhances the spatial growth of convective disturbances and, if sufficiently strong, activates the K-H mechanism to the point the shear layer evolution is driven by this instability instead of classic T-S waves typical of ZPG and low APG cases. We refer the reader to Appendix C for further discussion.

5.3. Effect of separation on the classic lift-up mechanism

Another class of convective-type instabilities of boundary layers is steady $(\omega =0)$ and 3-D $(\beta \neq 0)$ . In ZPG (case A) and low APG (case B) boundary layers, there is a region in the gain spectrum extending over a large $\beta$ range where the classic lift-up mechanism stands out as the most amplified mode (figure 3 a,b; cases A and B). This non-modal mechanism is well-documented in the literature (Andersson et al. Reference Andersson, Brandt, Bottaro and Henningson2001; Jacobs & Durbin Reference Jacobs and Durbin2001; Schmid & Henningson Reference Schmid and Henningson2001; Symon et al. Reference Symon, Rosenberg, Dawson and McKeon2018; Rigas et al. Reference Rigas, Sipp and Colonius2021) and involves the transient growth of streamwise-oriented $u'$ streaks optimally forced by near-wall streamwise vortices (see figures 28 and 29 in Appendix D). The physical process is characterised by the exchange of mean streamwise momentum across the boundary layer thickness due to the action of streamwise vorticity which pumps high-momentum fluid particles living in the outer shear layer deep into the boundary layer and, conversely, lifts low-momentum particles from the near-wall up to the outer layer. Notably, all the perturbation kinetic energy of the streaky mode plotted in figure 4(b) comes from the $u'$ component, since $v'$ and $w'$ are typically nil – refer to figure 29. The maximum gain of the lift-up mode occurs at $\beta =100\times 10^{-5}$ and drops more steeply for lower $\beta$ than for higher $\beta$ , indicating this instability is strongly 3-D.

The impact of the pressure gradient on the classic lift-up mechanism is now discussed. As long as the boundary layer remains attached, the process is continuous, meaning streamwise vorticity continually forces the $u'$ response whose spatial growth is almost monotonic – and is perfectly monotonic in ZPG boundary layer (figures 28 and 29). The local changes in spatial growth of cases B and C are due to the influence of the APG (figure 4 b). Interestingly, when the separation bubble forms (case D), the spatially continuous action of the streamwise vortices is ‘broken’ into two parts, one upstream of separation that creates local $u'$ disturbances over the bubble, and the other in the reattached shear layer downstream of the bubble leading to elongated streaks (see figure 6). The $u'$ structures in the separation zone imply a local production of perturbation kinetic energy, as seen in figure 4(b), where we introduce the terminology ‘double’ lift-up to indicate the two-stage action of this process due to the presence of the bubble. We notice qualitative similarity to Casacuberta et al. (Reference Casacuberta, Hickel, Westerbeek and Kotsonis2022, Reference Casacuberta, Hickel and Kotsonis2024) as concerns the phase shift of the streaks over and after the separation, which the authors observed before and after the forward-facing step. The mechanism described therein and termed ‘reverse lift-up’ is different, though, since it involves the stabilisation of the streaks not present in our LSB set-up.

In brief, separation affects the spatial organisation of the structures typical of lift-up but the physical mechanism remains the same as the ZPG boundary layer’s.

Figure 6. ‘Double’ lift-up mechanism at $\beta =150\times 10^{-5}$ (diamond marker in figure 3, D). (a) Amplitudes of the streamwise component of the curl of forcing $({\boldsymbol{\nabla}} \times {\boldsymbol{f}}')_{x}$ (black-solid line with $\times$ ) and of the streamwise velocity of the response $u'$ (blue-solid line with $+$ ). (b) Isosurfaces of $({\boldsymbol{\nabla}} \times {\boldsymbol{f}}')_{x}$ (coloured with yellow/black for positive/negative) and $u'$ (coloured with red/blue for positive/negative). The LSB is shown by the $u=0$ grey-coloured isosurface.

5.4. Modal centrifugal instability mechanism in LSB

A mechanism that is absent in attached boundary layers is the modal centrifugal instability related to the separation bubble. As reported in § 4.2 from the calculations presented in Appendix B, this mechanism becomes globally unstable for sufficiently strong recirculating flow inside the bubble ( $u_{\it{rev}}\approx 7{-}8\, \%$ ), which is ultimately dictated by the magnitude of the APG. As shown in figure 7(a), the temporal amplification factor of the least stable eigenvalue (which is stationary, $\textrm{Im} \{ \lambda \}=0$ ) computed over a range of $\beta$ for five LSB base-flows with $2\, \% \leqslant u_{\it{rev}} \leqslant 6.60\, \%$ is always negative, confirming that there is indeed no globally unstable mode in bubbles with $u_{\it{rev}}\lt 7{-}8\, \%$ . However, for $\beta$ values one to two orders of magnitude lower than those typical of streaks, the least stable eigenvalue is less damped. The RA corroborates the GSA data. The gain in figure 7(b) shows two regions of amplification, one with a peak at $\beta =10\times 10^{-5}$ that is very sensitive to the size of the bubble, and another at $\beta =150\times 10^{-5}$ related to lift-up which has been already discussed in § 5.3. Remarkably, GSA fails to capture the lift-up mode because it is non-modal, while the agreement between RA and GSA on the presence of a modal mechanism for low $\beta$ is examined further in figure 7(d) for the $u_{\it{rev}}=2\, \%$ bubble (case D). The similarity between the optimal response from RA and the eigenmode from GSA in figure 7(d) is striking. Either modes display the typical shape of the streamwise vorticity component $(\omega _{x}^{\prime})$ of the centrifugal instability over the separation zone, which is especially intense at the separation and reattachment points of the shear layer.

We highlight the important differences between the structure of the centrifugal and lift-up modes, since these may not be obvious. As far the forcing mode shape is concerned, we refer to figure 7(c) and figure 28(g,h), and note that:

1. (i) in the centrifugal mechanism both $({\boldsymbol{\nabla}} \times {\boldsymbol{f}}')_{x}$ and $({\boldsymbol{\nabla}} \times {\boldsymbol{f}}')_{z}$ are present, while for lift-up it is only $({\boldsymbol{\nabla}} \times {\boldsymbol{f}}')_{x}$ ;

2. (ii) the optimal forcing of the centrifugal mode is particularly strong near the separation point, which is not observed in lift-up.

In terms of the response mode shape (figures 7 d and 29 g,h):

1. (i) the wall-normal vorticity disturbance $(\omega _{y}^{\prime})$ is weak in the centrifugal instability, which is not the case of the lift-up, where $\omega _{y}^{\prime}\sim \omega _{z}^{\prime}$ ;

2. (ii) most importantly, $\omega^{\prime}_{x}$ is virtually negligible in the lift-up mode, but is an essential feature of the centrifugal instability.

Overall, the centrifugal instability originates from the bubble, while lift-up does not. Centrifugal forcing optimally acts on the separation point and the corresponding response is maximal at reattachment. Therefore, this mechanism essentially occurs over one bubble length, as opposed to lift-up which manifests throughout the boundary layer development.

Figure 7. (a) Temporal amplification of the least stable eigenmode from GSA and (b) squared gain from RA at $\omega =0$ for several $\beta$ values and for five bubbles with increasing $u_{\it{rev}}$ . (c) Componentwise amplitudes of curl of forcing computed from RA and (d) response vorticity computed from both RA and GSA at $\beta =10\times 10^{-5}$ , marked on panels (a) and (b) with red-dashed line, and $u_{\it{rev}}=2\, \%$ (case D). The eigenmode is scaled such that its total energy $\int E(x) \: \textrm{d}x$ matches the resolvent response mode’s.

In summary, linear analysis via both GSA and RA revealed that under APG conditions the K-H instability dominates over T-S waves and, when separation occurs, a new instability of centrifugal nature appears at low $\beta$ . The lift-up mechanism leading to streaks is present for any pressure gradient magnitude. We now examine in § 6 the nonlinear, multimodal interactions of the above disturbances to understand the sequence of mechanisms leading to turbulence.

6.1. Quasilinear-in- $\omega$ and HBNS-in- $\beta$ systems

We explore HBNS systems of variable architecture by varying the number of harmonics in the Fourier expansion (2.11) in order to identify and converge the leading modes that carry most of the energy of the disturbance field. To assess convergence of the solution, we track the behaviour of the total drag increase on the plate, $\Delta C_D$ in (2.14), resulting from the mean-flow modification.

In figure 8 we show the $\Delta C_D$ parameter obtained in the range of forcing amplitudes $A=[0.01,10]\times 10^{-7}$ . In figure 8(a) we compare systems with one ( $N=1$ ) against two ( $N=2$ ) time harmonics. We also vary the number of spanwise harmonics from a minimum of two ( $M=2$ ) to a maximum of six ( $M=6$ ). While the number of spanwise harmonics visibly affects the drag, systems with $N=2$ harmonics do not show significant difference to $N=1$ . This implies that the contribution of harmonics with frequency $2\omega$ (and, consequently, also higher-order harmonics $3\omega ,4\omega ,5\omega ,\ldots$ ) to the mean-flow modification is negligible.

Figure 8. (a) Comparison of $\Delta C_D$ for systems with $N=1$ and $N=2$ ; (b) $\Delta C_D$ of quasilinear-in- $\omega$ and HBNS-in- $\beta$ systems.

The overall effect of varying $M$ (at constant $N=1$ ) on $\Delta C_D$ is seen in figure 8(b), where $M$ ranges from 2 to 10. We define these systems as quasilinear-in- $\omega$ and HBNS-in- $\beta$ since they contain harmonics with at most $1\omega$ and a finite number of $m\beta$ . Notably, convergence of $\Delta C_D$ is observed if $M=8$ at least, meaning we resolve the most energetic disturbances and underlying nonlinear interactions within this range. We therefore show that the cascade of energy from low-order (large-scale) to high-order (small-scale) disturbances is dominant in space ( $\beta$ -modes) rather than in time ( $\omega$ -modes). For this reason we perform our analysis on quasilinear-in- $\omega$ ( $N=1$ ) systems for the remainder of the manuscript.

6.2. Nonlinear mechanisms at the early stages of transition

We analyse the onset of nonlinearity in the shear layer dynamics at weakly nonlinear forcing amplitude $A=1\times 10^{-7}$ where the base-flow modification is small but non-zero. This specific study is relevant to the early stages of transition, where the dynamics of the system can no longer be approximated as linear but, at the same time, only few harmonics (or scales) of the flow are active. The HBNS $_{2,1}$ system is suitable for this purpose since, by construction, it can calculate the first nonlinear triadic interaction where $(\pm 1\beta ,\pm 1\omega )$ waves produce the 3-D steady $(\pm 2\beta ,0\omega )$ mode. Actually, the fundamental oblique mode may also produce planar $(0\beta ,\pm 2\omega )$ and oblique $(\pm 2\beta ,\pm 2\omega )$ superharmonics but these are shown to have little impact on $\Delta C_D$ in figure 8(a). Furthermore, a low-order system such as HBNS $_{2,1}$ is sufficiently accurate at this forcing amplitude since all systems converge to the same $\Delta C_D$ . Imposing symmetry about $z=0$ (i.e. $\hat {\boldsymbol{w}}_{-m,n}=[\hat {u}, \hat {v}, -\hat {w}]_{m,n}$ and $\,\hat {\!\boldsymbol{f}}_{-m,n}=[\hat {f}_x, \hat {f}_y, -\hat {f}_z]_{m,n}$ ), the number of coupled equations in this system amounts to $(M+1)\times (N+1)=6$ , including the equation for the modification of the mean-flow by the disturbances.

The optimal nonlinear forcing/response solution in figure 9 shows that the shear layer initially amplifies the K-H mechanism over the separation region which then leads to the $(\pm 2\beta ,0\omega )$ mode through nonlinearity. We now explain this scenario in more detail.

Figure 9. Optimal nonlinear mechanisms computed from the HBNS $_{2,1}$ system by forcing at $(\beta ,\omega )=(20,25)\times 10^{-5}$ and $A=1\times 10^{-7}$ . Componentwise amplitudes of (a) external forcing $(\pm 1\beta ,\pm 1\omega )$ and (b) curl, (c) velocity and (d) vorticity of the K-H response $(\pm 1\beta ,\pm 1\omega )$ , (e) forcing and (f) curl of forcing from nonlinear triadic interaction $(\pm 1\beta ,-1\omega ) + (\pm 1\beta ,+1\omega ) = (\pm 2\beta ,0\omega )$ , (g) velocity and (h) vorticity of the centrifugal/lift-up response $(\pm 2\beta ,0\omega )$ . Time- and spanwise-averaged separation (grey-solid) and reattachment (grey-dashed) locations superimposed. Lines of best fit, $a\exp {(bx)}$ , for $|u'|$ ( $a=3.58\times 10^{-11}, b=2.01\times 10^{-4}, R^2 = 0.9999$ ) and $|\omega _{z}^{\prime}|$ ( $a=8.85\times 10^{-13}, b=1.82\times 10^{-4}, R^2 = 0.9977$ ) are plotted in red-solid.

Forcing optimally originates only from oblique waves (with planar wave and steady disturbance forcing harmonics being virtually zero, therefore not displayed) with a strong $({\boldsymbol{\nabla}} \times \boldsymbol{f}')_{z}$ component located upstream of the separation point – see figure 9(a,b). Their spatial structure, plotted in figure 10(a), is characterised by rollers of finite span tilted against the mean shear. The structure of the optimised forcing is visually similar to the linear resolvent forcing mode of figure 5(a iii). This result indicates that the 3-D K-H instability is the key to trigger mean-flow modification, which is in the cost functional (2.13), therefore justifying the observed similarity. This result also confirms the effectiveness of oblique waves in generating efficient mechanisms for transition to turbulence of both attached and separated shear flows (Schmid & Henningson Reference Schmid and Henningson1992; Marxen et al. Reference Marxen, Lang and Rist2013; Rigas et al. Reference Rigas, Sipp and Colonius2021; Dwivedi, Sidharth & Jovanović Reference Dwivedi, Sidharth and Jovanović2022). The K-H response grows in the separated shear layer with a spatial amplification rate that is well approximated by an exponential curve – see figure 9(c,d), indicating the dynamics of the flow is essentially governed by a linear mechanism over the entire front portion of the separation bubble up to the location of the bubble’s apex, in agreement with Diwan & Ramesh (Reference Diwan and Ramesh2009). Because these waves are spanwise-periodic, they form a chequerboard pattern of rollers with finite spanwise length, as shown in figure 10(a).

Figure 10. The 3-D structure of the modes in figure 9. Isosurfaces of (a) $({\boldsymbol{\nabla}} \times \boldsymbol{f}')_{z}$ (yellow/black for positive/negative) and $\omega _{z}^{\prime}$ (red/blue for positive/negative) from the K-H mechanism, (b) $({\boldsymbol{\nabla}} \times \boldsymbol{f}')_{x}$ (magenta/cyan for positive/negative) from the nonlinear internal forcing and (c) $\omega^{\prime}_{x}$ (magenta/cyan for positive/negative) and $u'$ (red/blue for positive/negative) from the centrifugal/lift-up mechanisms. The instantaneous separation bubble is plotted with the $u=0$ isosurface.

Near the mean reattachment location the K-H waves reach finite amplitude and saturation effects occur. Concurrently, the nonlinear interaction of these waves governed by the $\unicode{x1D649}(\boldsymbol{w},\boldsymbol{w})$ operator in (2.2) produces an internal forcing term (equivalent to internal stresses) with a strong $({\boldsymbol{\nabla}} \times \boldsymbol{f}')_{x}$ component reaching its peak energy shortly downstream of reattachment – figure 9(e,f). This term is computed in the equation for the $(2\beta ,0\omega )$ harmonic ( $m=2$ , $n=0$ ) from the HBNS $_{2,1}$ system,

(6.1) \begin{equation} [ \unicode{x1D647}_{2} + \unicode{x1D649}_{\,0}^{2} ( \hat {\boldsymbol{w}}_{0,0}, \boldsymbol{\cdot} ) ] \hat {\boldsymbol{w}}_{2,0} + \unicode{x1D649}_{1}^{1} ( \hat {\boldsymbol{w}}_{1,-1}, \hat {\boldsymbol{w}}_{1,1} ) = 0, \end{equation}

where $\unicode{x1D649}_{1}^{1} ( \hat {\boldsymbol{w}}_{1,-1}, \hat {\boldsymbol{w}}_{1,1} )$ governs the nonlinear interaction of $(1\beta ,-1\omega )$ with $(1\beta ,1\omega )$ waves through the convective term of the N-S. From now onwards, $\,\hat {\!\boldsymbol{f}}_{\text{K-H}} = [ \hat {f}_{\text{K-H},x}, \hat {f}_{\text{K-H},y}, \hat {f}_{\text{K-H},z} ]^{\top } = -\unicode{x1D64B}^{\top } \unicode{x1D649}_{1}^{1} ( \hat {\boldsymbol{w}}_{1,-1}, \hat {\boldsymbol{w}}_{1,1} )$ will be referred to as ‘K-H forcing’ because it originates from the nonlinear interaction of two K-H waves (for the explicit derivation see Appendix F). $\unicode{x1D64B}^{\top }$ is the restriction operator mapping the full velocity-pressure state $[u,v,w,p]^{\top }$ to the velocity state $[u,v,w]^{\top }$ .

The $(\pm 2\beta ,0\omega )$ response in figure 9(g,h) displays two characteristic regions. By comparison with the linear centrifugal and lift-up modes in figures 7(d) to 29(g,h), we notice that region I displays the features of centrifugal instability, which are that $\omega^{\prime}_{x}$ peaks near the reattachment point and that it is of the same order of magnitude of $\omega _{z}^{\prime}$ , and region II is the site of lift-up where both $\omega _{y}^{\prime}$ and $\omega _{z}^{\prime}$ are amplified but the $\omega^{\prime}_{x}$ response is weak. The structures of $\omega^{\prime}_{x}$ in figure 10(c) further highlight the presence of two mechanisms in the $(\pm 2\beta ,0\omega )$ mode. The first set of streamwise vortices localised around the mean reattachment point is related to the centrifugal mode pertaining to the separation zone (Sansica, Sandham & Hu Reference Sansica, Sandham and Hu2016; Hildebrand et al. Reference Hildebrand, Dwivedi, Nichols, Jovanović and Candler2018; Mauriello, Larchevêque & Dupont Reference Mauriello, Larchevêque and Dupont2022). At reattachment, the flow is strongly curved due to the displacement effect of the bubble, which makes it a susceptive site for the development of centrifugal instabilities of Görtler type (Görtler Reference Görtler1941; Saric Reference Saric1994; Li & Malik Reference Li and Malik1995). A detailed study evaluating sufficient conditions for centrifugal instability is provided in Appendix E, where we essentially demonstrate that streamline curvature at separation and reattachment supports a locally unstable centrifugal mode. On the other hand, the second set of streamwise-elongated vortices and $u'$ streaks are ascribed to lift-up (Jacobs & Durbin Reference Jacobs and Durbin2001; Balamurugan & Mandal Reference Balamurugan and Mandal2017; Rigas et al. Reference Rigas, Sipp and Colonius2021) taking place in the redeveloping boundary layer.

In summary:

1. (i) oblique wave-like disturbances optimally located upstream of separation excite the inviscid, inflectional K-H instability at $(1\beta ,1\omega )$ through a linear mechanism that takes advantage of the shear present in the separated shear layer;

2. (ii) at saturation amplitude, the K-H waves interact nonlinearly to produce a source of internal forcing at $(2\beta ,0\omega )$ that features mainly streamwise vortices;

3. (iii) streamwise vorticity manifests in the $(2\beta ,0\omega )$ response near the reattachment point, indicating the presence of a centrifugal instability;

4. (iv) the presence of $\omega^{\prime}_{x}$ due to the centrifugal response translates into streamwise vortices that trigger lift-up and, therefore, streaks in the redeveloping boundary layer.

While there are commonalities with the oblique wave-to-streaks mechanism found by Schmid & Henningson (Reference Schmid and Henningson1992) and Rigas et al. (Reference Rigas, Sipp and Colonius2021) for attached boundary layers, the presence of separation makes this scenario more complex due to the centrifugal instability.

Finally, it is important to note that the K-H induced forcing is different to the optimal forcing of either the centrifugal mode in figure 7(c) and lift-up in figure 28(g,h). In the former, we have observed a strong streamwise vorticity footprint at the separation point, while in the latter a continuous action of streamwise vorticity in the streamwise direction. Despite these differences, the K-H forcing must have some projection on the optimal forcing mode responsible for the centrifugal instability such that this mechanism can be excited. This will be examined in detail in § 6.3.

6.3. Resolvent-based reconstruction of the weakly nonlinear centrifugal instability mechanism

Broadly following the approach of Dwivedi et al. (Reference Dwivedi, Sidharth and Jovanović2022) who derived the nonlinear $(\pm 2\beta ,0\omega )$ forcing from weakly nonlinear analysis, we investigate the capability of the resolvent basis to approximate the nonlinear forcing/response mechanism at $(\pm 2\beta ,0\omega )$ for low amplitude forcing and low reversed flow. In doing so, we aim to answer the following questions: ‘Is this mechanism low-rank’ and ‘how similar/different is the K-H forcing to the optimal forcing of the $(\pm 2\beta ,0\omega )$ mode’? Our procedure, although sharing many commonalities with Dwivedi et al. (Reference Dwivedi, Sidharth and Jovanović2022), employs the mean-flow (not the laminar base-flow) to extract the resolvent operator. This is to account for any changes in the base-flow topology which are readily calculated by HBNS.

Following the linear resolvent framework from (2.7) to (2.10), the $(2\beta ,0\omega )$ response may be written as

(6.2) \begin{equation} \hat {\boldsymbol{w}}_{2,0} = \unicode{x1D643} \left ( 2\beta ,0\omega \right ) \unicode{x1D648} \unicode{x1D64B} \,\hat {\!\boldsymbol{f}}_{2,0}, \end{equation}

where $\unicode{x1D643} ( 2\beta ,0\omega ) = [ -\unicode{x1D63C}_{2} (\hat {\boldsymbol{w}}_{0,0}) ]^{-1}$ is the resolvent operator extracted from the time- and spanwise-averaged flow at amplitude $A=1\times 10^{-7}$ and $\,\hat {\!\boldsymbol{f}}_{2,0} = \,\hat {\!\boldsymbol{f}}_{\text{K-H}}$ is the K-H forcing acting on the $(2\beta ,0\omega )$ harmonic. Considering the set of orthonormal forcing $ [ \,\hat {\!\boldsymbol{f}}_1,\,\hat {\!\boldsymbol{f}}_2,\ldots ,\,\hat {\!\boldsymbol{f}}_n ]$ and response $ [ \hat {\boldsymbol{u}}_1,\hat {\boldsymbol{u}}_2,\ldots ,\hat {\boldsymbol{u}}_n ]$ modes, where $\hat {\boldsymbol{u}}_{i}$ contains only the three components of velocity, and the associated energy gains $ ( \sigma _1,\sigma _2,\ldots ,\sigma _n )$ (analogous to ${\mathcal{G}}$ ), and further letting

(6.3) \begin{equation} \alpha _{i} = \frac { \langle \,\hat {\!\boldsymbol{f}}_{i},\,\hat {\!\boldsymbol{f}}_{\text{K-H}} \rangle _{\Omega _{\,\hat {\!\boldsymbol{f}}}}}{\sqrt { \langle \,\hat {\!\boldsymbol{f}}_{\text{K-H}},\,\hat {\!\boldsymbol{f}}_{\text{K-H}} \rangle _{\Omega _{\,\hat {\!\boldsymbol{f}}}}}} \end{equation}

be the scalar projection coefficient of the $i$ th forcing vector on the K-H forcing, the nonlinear $(2\beta ,0\omega )$ response can be written as

(6.4) \begin{equation} \hat {\boldsymbol{w}}_{2,0} = \sum _{i=1}^{n\rightarrow \infty } \hat {\boldsymbol{u}}_{i} \sigma _{i} \alpha _{i} \sqrt {\langle \,\hat {\!\boldsymbol{f}}_{\text{K-H}},\,\hat {\!\boldsymbol{f}}_{\text{K-H}} \rangle _{\Omega _{\,\hat {\!\boldsymbol{f}}}}}, \end{equation}

and, similarly, the K-H forcing as

(6.5) \begin{equation} \,\hat {\!\boldsymbol{f}}_{\text{K-H}} = \sum _{i=1}^{n\rightarrow \infty } \,\hat {\!\boldsymbol{f}}_{i} \alpha _{i} \sqrt { \langle \,\hat {\!\boldsymbol{f}}_{\text{K-H}},\,\hat {\!\boldsymbol{f}}_{\text{K-H}} \rangle _{\Omega _{\,\hat {\!\boldsymbol{f}}}}}. \end{equation}

Approximations of the nonlinear K-H forcing $\,\tilde {\!\boldsymbol{f}}_{\text{K-H}}$ and response $\tilde {\boldsymbol{w}}_{2,0}$ can be reconstructed with a finite set of modes ( $n=N_{m}\lt \infty$ ).

A suitable metric to assess the accuracy of the approximation is the relative error,

(6.6) \begin{equation} \varepsilon _{\,\hat {\!\boldsymbol{f}}} = \frac {\sqrt { \langle \,\hat {\!\boldsymbol{f}}_{\text{K-H}}-\,\tilde {\!\boldsymbol{f}}_{\text{K-H}},\,\hat {\!\boldsymbol{f}}_{\text{K-H}}-\,\tilde {\!\boldsymbol{f}}_{\text{K-H}} \rangle _{\Omega _{\,\hat {\!\boldsymbol{f}}}}}}{\sqrt { \langle \,\hat {\!\boldsymbol{f}}_{\text{K-H}},\,\hat {\!\boldsymbol{f}}_{\text{K-H}} \rangle _{\Omega _{\,\hat {\!\boldsymbol{f}}}}}}, \end{equation}

for the K-H forcing, and

(6.7) \begin{equation} \varepsilon _{\hat {\boldsymbol{w}}} = \frac {\sqrt {\left \langle \hat {\boldsymbol{w}}_{2,0}-\tilde {\boldsymbol{w}}_{2,0},\hat {\boldsymbol{w}}_{2,0}-\tilde {\boldsymbol{w}}_{2,0} \right \rangle _{\Omega _{\hat {\boldsymbol{w}}}}}}{\sqrt {\left \langle \hat {\boldsymbol{w}}_{2,0},\hat {\boldsymbol{w}}_{2,0} \right \rangle _{\Omega _{\hat {\boldsymbol{w}}}}}}, \end{equation}

for the nonlinear $(2\beta ,0\omega )$ response. It can also be shown (see Appendix G for the derivation) that the relative error of the response approximation is upper bounded according to

(6.8) \begin{equation} \varepsilon _{\hat {\boldsymbol{w}}} \leqslant \frac {\sigma _{N_{m}+1}}{\overline {\sigma }_{N_{m}}} \frac {\varepsilon _{\,\hat {\!\boldsymbol{f}}}}{\sqrt {1-\varepsilon _{\,\hat {\!\boldsymbol{f}}}^{2}}} \leqslant \frac {\varepsilon _{\,\hat {\!\boldsymbol{f}}}}{\sqrt {1-\varepsilon _{\,\hat {\!\boldsymbol{f}}}^{2}}}, \end{equation}

where

(6.9) \begin{equation} \overline {\sigma }_{N_{m}}^{2} = \frac {\sum _{i=1}^{N_{m}} \sigma _{i}^2 \left | \alpha _{i} \right |^{2}}{\sum _{i=1}^{N_{m}} \left | \alpha _{i} \right |^{2}}, \end{equation}

is the weighted average of the $\sigma _{i}$ distribution in the range $[1,N_{m}]$ and $\sigma _{N_{m}+1}$ is the gain of the first truncated mode sitting outside of the prescribed range. The ratio $\sigma _{N_{m}+1}/\overline {\sigma }_{N_{m}}\leqslant 1$ indicates the distance of the first truncated mode from the weighted average and its value decreases with $N_{m}$ . Effectively, the smaller the ratio, the more energy of the forcing/response mechanism is captured by the resolvent basis expansion.

The domain of the linear forcing is capped at $y=4\times 10^3$ to avoid spurious structures appearing in the far field. We also consider two cases, one where the forcing is not restricted in $x$ and another where it is calculated within $[x_{S}+0.7L_{b},x_{out}]$ , i.e. from 70 % of the bubble length to the outlet. This is an indicative region coinciding with incipient shear layer reattachment and the nonlinear interaction of K-H waves to mimic the true nonlinear mechanism. We outline the reconstruction performance of both test cases in figure 11. As the behaviour of the squared energy gain in the top row suggests, the linearised operator in (6.2) is not low-rank, implying the underlying mechanism cannot be approximated accurately with a low-order linear model. The restricted forcing case achieves better projections overall, especially within the leading three modes – see figure 11(b). While the behaviour of $\alpha$ is locally erratic, we observe a decay on average, meaning the projection degrades for high-order modes. Overall, the relative error of the forcing approximation $\varepsilon _{\,\hat {\!\boldsymbol{f}}}$ is large – see figure 11(c). The restricted forcing case is clearly better than the unrestricted one, however, even with 20 modes, the error settles around 80 % and 90 %, respectively. These results highlight the limitations of linear RA in modelling the internal stresses produced by the nonlinear dynamics of this LSB flow even at low amplitude. In figure 12 we compare the shape of the true and approximated K-H forcing with streamwise spatial restriction for the streamwise component of both the forcing itself and the curl of the forcing. The ground truth (figure 12 a,e) is characterised by structures emerging in the site of reattachment, which is where saturation of K-H waves occurs – refer to figure 9(c–f). The approximation with one mode (figure 12 b,f) falls short in that both the shape and the amplitude are extremely inaccurate, leading to approximately 95 % error. The reconstruction improves with 20 modes (figure 12 c,g) but the error of the $x$ -component of the curl is still large, which is attributed to the $y$ and $z$ components of the forcing vector – see figure 31 in Appendix G for more details.

Figure 11. (a) squared resolvent gain, (b) modulus of the projection coefficient, (c) relative error of the K-H forcing approximation and (d) of the nonlinear response approximation with error bound $ {\sigma _{N_{m}+1}}/{\overline {\sigma }_{N_{m}}} {\varepsilon _{\,\hat {\!\boldsymbol{f}}}}/{\sqrt {1-\varepsilon _{\,\hat {\!\boldsymbol{f}}}^{2}}}$ plotted with a dashed line.

Figure 12. Resolvent-based reconstruction of the weakly nonlinear ( $A=1\times 10^{-7}$ ) K-H forcing $\,\hat {\!\boldsymbol{f}}_{\text{K-H}}$ from the restricted forcing case: (a,e) true nonlinear forcing, reconstruction with (b,f) $N_{m}=1$ mode, (c,g) $N_{m}=20$ modes and (d,h) amplitudes. Streamwise components of the forcing (a–d) and curl of forcing (e–h) shown.

Despite the poor reconstruction of the K-H forcing, the approximation of the nonlinear $(2\beta ,0\omega )$ response is surprisingly good. With only four modes in the unrestricted forcing case, we achieve an error of approximately 25 %, while with only two modes in the restricted forcing case, we reach an error of 15 %. With seven modes the error drops to approximately 15 % and 10 % in the unrestricted and restricted cases, respectively – see figure 11(d). The $u'$ structure of the response is qualitatively captured by only the first resolvent mode, but with 20 modes the improvement is remarkable since we observe virtually no mismatch in the amplitude (figure 13 a–d). The streamwise vorticity of the response computed as $\tilde {\omega }_{x,2,0} = \partial _y \tilde {w}_{2,0} - i\beta \tilde {v}_{2,0}$ is less accurate compared with the $u'$ reconstruction (figure 13 e–h). While one mode is insufficient to recover the vortical structures even qualitatively, 20 modes improve the accuracy but the error is still above 20 % (refer to figure 31 c).

Figure 13. Same as figure 12 but for the nonlinear response $\hat {\boldsymbol{w}}_{2,0}$ .

The results show that many modes are required to reconstruct the full streamwise vorticity triggered by centrifugal instability. We can therefore conclude that centrifugal instability is not due to a single resolvent mode, but to multiple modes, which all exhibit streamwise vorticity.

6.4. Role of centrifugal instability for laminar–turbulent transition

At high forcing amplitude ( $A=10\times 10^{-7}$ ), the centrifugal instability arising in region I of figure 9 becomes more unstable and grows farther upstream inside the separation zone (see figure 14), a feature observed in the global eigenmode representative of modal centrifugal instability of the bubble – figures 7 to 25.

Figure 14. Amplitude of $\omega^{\prime}_{x}$ of the nonlinear $(\pm 2\beta ,0\omega )$ response at low and high forcing amplitudes.

In the remainder of the manuscript, we describe the role of this instability for breakdown of the main K-H roller shed by the shear layer, which paves the way to turbulence. We highlight that at high amplitude weakly nonlinear approximations fail and other methods, such as DNS, become necessary. Our HBNS methodology offers a self-contained tool to also explore high amplitude disturbance mechanisms of the late transition stages. We recall from figure 8 that larger systems are required for these regimes, hence, we employ 10 harmonics in space ( $M=10$ ) and one in time ( $N=1$ ) for the analysis hereafter.

6.5. Spanwise vortex distortion and breakdown

In order to describe the sequence of events leading to breakdown of the shear layer, we reconstruct the spatiotemporal flow obtained from HBNS $_{10,1}$ system at $A=10\times 10^{-7}$ . Spanwise-periodic oscillations of wavelength $\lambda _z=2\pi /\beta$ (approximately equal to the mean separation length, $L_b$ ) appear in the separation bubble and are ascribed to the K-H mode which is activated in the fore part of the bubble (figure 15 a). Mutually, a roller of finite span ( $=\lambda _z/2$ ) detaches from the shear layer at $x\approx 1.1\times 10^5$ (figure 15 b), indicative of the roll-up process (Marxen et al. Reference Marxen, Lang and Rist2013; Michelis et al. Reference Michelis, Kotsonis and Yarusevych2018a,b ; Hosseinverdi & Fasel Reference Hosseinverdi and Fasel2020). Spanwise vortices are shed downstream at fixed temporal frequency equal to the K-H instability frequency, confirming that the main vortex shedding dynamics is driven by the K-H mechanism and hence explaining the success of linear stability analysis in recovering its characteristic frequency (Simoni et al. Reference Simoni, Ubaldi and Zunino2013; Yarusevych & Kotsonis Reference Yarusevych and Kotsonis2017; Ziade et al. Reference Ziade, Feero, Lavoie and Sullivan2018; Michelis et al. Reference Michelis, Kotsonis and Yarusevych2018a , Reference Michelis, Yarusevych and Kotsonisb ).

Figure 15. Reconstruction of the instantaneous flow at $A=10\times 10^{-7}$ from HBNS $_{10,1}$ . Isosurfaces of (a,c) $u=0$ and (b,d) Q-criterion (isovalue, $4\times 10^{-11}$ ) coloured with $u$ within the spanwise domain marked by the red-solid lines. (a,b) Time- and spanwise-averaged flow + $(\pm 1\beta ,\pm 1\omega )$ harmonic, (c,d) time- and spanwise-averaged flow + $(\pm 1\beta ,\pm 1\omega )$ + $(\pm 2\beta ,0\omega )$ + $(\pm 3\beta ,\pm 1\omega )$ harmonics.

Spanwise K-H rollers are prone to destabilise due to the action of the centrifugal instability discussed earlier. More specifically, the net effect of the interaction between $(\pm 1\beta ,\pm 1\omega )$ and $(\pm 2\beta ,0\omega )$ modes, which gives $(\pm 3\beta ,\pm 1\omega )$ through a suitable triad, is the distortion of the spanwise vortex, as seen in figure 15(c,d). The K-H rollers developing within the shear layer are progressively distorted into structures that are reminiscent of $\Lambda$ -vortices (figure 15 d). The uplift of (negative) spanwise vorticity at $x\approx 1.3-1.4\times 10^5$ , which, in turn, promotes mixing between the boundary layer and the outer flow, is a direct consequence of this mechanism. Therefore, it is observed that the $(\pm 1\beta ,\pm 1\omega ) + (\pm 2\beta ,0\omega ) = (\pm 3\beta ,\pm 1\omega )$ triadic interaction is key for the initiation of the main K-H vortex breakdown which paves the way to turbulence. This mechanism is further illustrated in figure 16. Here we provide an overview of the complete sequence of events that characterises the breakdown of the K-H vortex and discuss the underlying physical mechanism. In figure 16(a) we show the 3-D structure of the transitional shear layer and in figure 16(b) four $y$ – $z$ sectional views extracted at $x=[1.1,1.2,1.32,1.5]\times 10^5$ . At the first station, the K-H roller is essentially 2-D. At the second station, the shape of the roller is visibly distorted since spanwise oscillations appear. At the third station, the deformation is enhanced as the roller is no longer aligned in the $z$ -direction. The central region is uplifted and stretched away from the sides along the streamwise direction. At the fourth station we observe a $\Lambda$ -shaped structure. The mechanism by which the roller deforms is driven by the centrifugal instability, appearing in the figure as a pair of counter-rotating streamwise vortices whose distance from the centres is $\lambda _z/4$ , i.e. approximately one quarter of the mean separation length, $L_b$ . These vortices act through the braid regions (located between two consecutive rollers) and induce upward flow that interacts with the spanwise roller, causing its distortion and disintegration (Marxen et al. Reference Marxen, Lang and Rist2013; Kumar & Sarkar Reference Kumar and Sarkar2023). This process shares some similarities with secondary flow structures, called ‘rib’ vortices, which are responsible for the breakdown of spanwise rollers in mixing layers and wakes (Lopez & Bulbeck Reference Lopez and Bulbeck1993; Sun et al. Reference Sun, Schrijer, Scarano and van Oudheusden2014; Yang Reference Yang2019). While we do not observe vortex loops revolving around the rollers, the mechanism whereby streamwise vorticity disintegrates the main vortex is present in our results, the origin of streamwise vortices being a centrifugal-type instability of convective nature discussed earlier.

Figure 16. (a) Isosurfaces of Q-criterion (isovalue, $4\times 10^{-11}$ ) coloured with $u$ from the full flow field and of (magenta/cyan) positive/negative $\omega^{\prime}_{x}$ from the $(\pm 2\beta ,0\omega )$ mode. (b) The $y$ – $z$ planar views extracted at four $x$ -stations $[1.1,1.2,1.32,1.5]\times 10^5$ , labelled with numbers $(1)$ to $(4)$ .

The complete breakdown of $\Lambda$ -structures establishes transition to turbulence of the shear layer. Notably, when the large scales of the flow breakdown to turbulence energy spreads to infinitely many (smaller) scales all the way down to the Kolmogorov microscales (Pope Reference Pope2000; George Reference George2013). Capturing these structures is outside the scope of this work since the proposed methodology aims at characterising the transitional stages of the flow given a finite set of harmonics retained in the model. This means we are able to resolve only part of the energy spectrum, the remainder being the truncation error associated with the expansion (2.11). Regardless of this limitation, the HBNS $_{10,1}$ system is able to capture the final stage of the $\Lambda$ -vortex breakdown – see figure 16. Additional insight into the main sites of turbulence production is offered by plotting the production of turbulent kinetic energy in figure 17. Two regions of the flow are identified, one is along the inflection line near incipient shear layer reattachment, which is where shear layer roll-up and K-H roller shedding take place. The correct estimation of production at this station is paramount for the prediction of the bubble topology and, consequently, the generation of turbulence in the redeveloping boundary layer downstream of reattachment, which constitute two major problems in RANS turbulence models applied to separated transitional shear layers (Bernardos et al. Reference Bernardos, Richez, Gleize and Gerolymos2019a ,Reference Bernardos, Richez, Gleize and Gerolymos b ). Along the vortex shedding line ( $y\approx 2\delta ^{*}$ ), but substantially downstream of reattachment, there are local patches of production, which match the locations of $\Lambda$ -vortex formation ( $x\approx 1.4{-}1.5\times 10^5$ ) and breakdown ( $x\approx 1.6{-}1.7\times 10^5$ ). Much closer to the wall ( $y\approx \delta ^{*}/5$ ), the thin layer of turbulence production is related to near-wall turbulence.

Figure 17. Contours of overall production of turbulent kinetic energy ${\mathcal{P}} = -\overline {u_{i}^{\prime }u_{j}^{\prime }} {\partial \langle u_i \rangle _{z,t}}/{\partial x_j}$ . The dividing streamline of the time- and spanwise-averaged separation bubble and inflection line are plotted in yellow and red, respectively.

6.6. Integral quantities

Although the accurate estimation of the boundary layer transition point is an extremely arduous task, it is often sought to inform turbulence models used in steady simulations for which it is clearly impossible to predict transition otherwise (Bernardos et al. Reference Bernardos, Richez, Gleize and Gerolymos2019a ,Reference Bernardos, Richez, Gleize and Gerolymos b ; Dellacasagrande et al. Reference Dellacasagrande, Lengani, Simoni and Ubaldi2024). Being equipped with spatiotemporal information of the flow, the transition point is yet not fixed in time/space and a single location cannot be identified. Nevertheless, we employ boundary layer integral analysis on the time- and spanwise-averaged flow to possibly draw a connection between the dynamic events of transition described thus far and the statistical characteristics of the boundary layer. In figure 18 we plot the skin friction coefficient $C_f$ of the time- and spanwise-averaged flow for five HBNS system architectures: HBNS $_{2,1}$ , HBNS $_{4,1}$ , HBNS $_{6,1}$ , HBNS $_{8,1}$ and HBNS $_{10,1}$ . The reference laminar (Blasius) and turbulent (White Reference White1991) curves for ZPG boundary layer are also superimposed together with the baseline $C_f$ from the laminar base-flow.

Figure 18. Skin friction coefficient evaluated for different quasilinear-in- $\omega$ systems ( $N=1,M\gt 1$ ) and three forcing amplitudes $A=[2,6,10]\times 10^{-7}$ . The laminar and turbulent (from White (Reference White1991)) ZPG boundary layer profiles are plotted with black dashed and dashed–dotted lines.

Firstly, the turbulent curve is reached only by systems with $M\geqslant 8$ , in agreement with the analysis in § 6.1. When the truncation error is large, the turbulent energy cascade cannot take place and the skin friction levels remain close to the baseline laminar curve (systems HBNS $_{2,1}$ and HBNS $_{4,1}$ ). The HBNS $_{6,1}$ system contains the minimum number of harmonics necessary to reach the turbulent curve at $A=10\times 10^5$ , but trajectories of $C_f$ at high forcing amplitudes converge only for $M \geqslant 8$ . In this case, the critical amplitude is found to be $A=6\times 10^{-7}$ . As the amplitude is progressively increased, the coordinate where the skin friction crosses the turbulent curve moves towards the mean reattachment point. A similar trend is also observed for increasing system’s order, suggesting that in the limit of zero truncation error, transition and reattachment should collapse to the same point, giving rise to the turbulent reattachment scenario observed in Hosseinverdi & Fasel (Reference Hosseinverdi and Fasel2019, Reference Hosseinverdi and Fasel2020).

From the skin friction results we notice some robustness of the HBNS solutions in the prediction of $C_f$ , at least in the transitional regime. Beyond this stage, i.e. where $C_f$ has crossed the turbulent curve, the turbulent reattached boundary layer developing downstream of the bubble requires an intractably large number of harmonics in both $z$ and $t$ to resolve the full spectrum of turbulence down to the dissipative scales. This is where the current state-of-the-art HBNS approach loses robustness and accuracy. An irregular, wavy behaviour of $C_f$ is observed in figure 19(a), which also displays the statistical state of the boundary layer downstream of the transitional region. Despite the oscillations, the skin friction remains around the turbulent level.

Figure 19. (a) Skin friction coefficient, (b) shape factor and (c) boundary layer mean velocity profiles at $A=10\times 10^{-7}$ . Mean reattachment and the streamwise coordinate at which the skin friction reaches the turbulent curve are plotted with grey-dashed and blue-solid lines, respectively. Profiles extracted at (blue) $x=1.2\times 10^5$ , (red) $x=1.35\times 10^5$ and (green) $x=1.5\times 10^5$ .

Another parameter we consider is the boundary layer shape factor $H=\delta ^{*}/\theta$ plotted in figure 19(b). The shape factor gradually departs from the laminar value ( $H=2.6$ for ZPG boundary layer, (Pope Reference Pope2000)) due to the influence of the pressure gradient. High values of $H$ are reached in the separated zone ( $H_{max} \approx 5.5$ ), later followed by a sharp drop due to the coexisting effects of boundary layer reattachment and onset of transition. In the transition region the shape factor crosses the laminar reference level and marches towards the turbulent range ( $1.3 \leqslant H \leqslant 1.5$ for ZPG boundary layer, (Pope Reference Pope2000)). Downstream of the transition region the shape factor settles around $H=1.6$ , which is, together with the $C_f$ curve, another indication of a turbulent-like boundary layer in the postreattachment region.

Finally, we show the boundary layer velocity profile in inner units extracted at three $x$ -stations marked with blue ( $x=1.2\times 10^5$ ), red ( $x=1.35\times 10^5$ ) and green ( $x=1.5\times 10^5$ ) arrows in figure 19(c). The profiles are superimposed on the viscous sublayer ( $y^{+}\lt 5$ ) and log-law ( $30\lt y^{+}\lt 300$ ) curves. A turbulent boundary layer is known to follow the universal log-law profile in the log region (Pope Reference Pope2000), hence we look at whether and how accurately the extracted profiles follow this trend. The velocity profile at $x=1.2\times 10^5$ is extracted from the early-transitional region, and, in agreement with $C_f$ and $H$ , it is not turbulent yet – in fact, it is still very close to the laminar state since the $u^{+} = y^{+}$ law is followed accurately in the viscous sublayer, but the log-layer slope is missing. The boundary layer is definitely not turbulent at this stage, also because, from a spatiotemporal point of view, the process of K-H roller distortion and breakdown is still at an early stage (refer to figures 15 and 16). At $x=1.35\times 10^5$ , the profile bends towards the log-law curve, portending the breakdown event, yet, the slope is not perfectly recovered. At this station, the K-H vortex has deformed sufficiently to form a $\Lambda$ . Ultimately, the log-layer slope is captured more accurately at $x=1.5\times 10^5$ , suggesting a more pronounced turbulent character of the boundary layer. Contextually, the $\Lambda$ -vortex is on the verge of breaking into small-scale turbulence.

It should be noted that a fully developed turbulent profile may not be achievable within our computational domain since the boundary layer relaxes to the turbulent equilibrium several separation bubble lengths downstream of the reattachment point, estimated to be around six by Alam & Sandham (Reference Alam and Sandham2000). Hosseinverdi & Fasel (Reference Hosseinverdi and Fasel2019) hypothesise this slow relaxation may be due to the persistence of coherent vortex structures reminiscent of the large-scale K-H rollers in the boundary layer downstream of reattachment. The final breakup of these large scales into small vortical structures generated near the wall yields profiles in better agreement with the fully developed turbulent boundary layer characteristics.