Edge tracking in spatially developing boundary layer flows

Abstract

Recent progress in understanding subcritical transition to turbulence is based on the concept of the edge, the manifold separating the basins of attraction of the laminar and the turbulent state. Originally developed in numerical studies of parallel shear flows with a linearly stable base flow, this concept is adapted here to the case of a spatially developing Blasius boundary layer. Longer time horizons fundamentally change the nature of the problem due to the loss of stability of the base flow due to Tollmien–Schlichting (TS) waves. We demonstrate, using a moving box technique, that efficient long-time tracking of edge trajectories is possible for the parameter range relevant to bypass transition, even if the asymptotic state itself remains out of reach. The flow along the edge trajectory features streak switching observed for the first time in the Blasius boundary layer. At long enough times, TS waves co-exist with the coherent structure characteristic of edge trajectories. In this situation we suggest a reinterpretation of the edge as a manifold dividing the state space between the two main types of boundary layer transition, i.e. bypass transition and classical transition.

1 Introduction

Understanding the onset of transition in boundary layer flows has always been an important challenge in aerodynamics because of the high drag associated with turbulent flows. It is crucially dependent on the level of ambient turbulence present in the free stream, both for aeronautic applications and in wind tunnels. An archetype for the investigation of boundary layer transition is the incompressible Blasius boundary layer. The corresponding laminar base flow is two-dimensional and develops spatially with the distance downstream. Two main transition scenarios dominate this flow case depending on the incoming level of perturbations. For weak incoming levels corresponding to flight mode in a calm atmosphere, the classical scenario predicts linear instability of the base flow to two-dimensional, spanwise-invariant Tollmien–Schlichting (TS) waves, followed by their further destabilisation (Schubauer & Skramstad 1948; Boiko et al. 1994; Schmid & Henningson 2001; Schlichting & Gersten 2016). For higher turbulence levels (typically for perturbation velocities exceeding 2 % of the free-stream velocity) akin to those found in turbomachines, the so-called bypass route starts at earlier streamwise locations. It features streamwise velocity defects (‘streaks’) sustained by streamwise vorticity via the lift-up mechanism, and growing rapidly downstream (Klebanoff 1955; Kline et al. 1967; Morkovin 1969). Those streaks support further instabilities (Henningson, Lundbladh & Johansson 1993; Kachanov 1994; Andersson et al. 2001; Jacobs & Durbin 2001; Brandt, Schlatter & Henningson 2004; Schlatter et al. 2008; Vaughan & Zaki 2011; Hack & Zaki 2014) leading to the nucleation of turbulent spots that invade the flow (Matsubara & Alfredsson 2001; Brandt et al. 2004; Kreilos et al. 2016). For sufficiently perturbed inflow, no TS wave has been convincingly reported in the bypass picture.

Linear stability of the Blasius boundary layer profile can be carried out relatively easily by freezing the base flow and parametrising it by a Reynolds number $Re_{\unicode[STIX]{x1D6FF}^{\ast }}=U_{\infty }\unicode[STIX]{x1D6FF}^{\ast }/\unicode[STIX]{x1D708}$ , where $\unicode[STIX]{x1D6FF}^{\ast }$ the displacement thickness for the laminar base flow, $U_{\infty }$ is the free-stream velocity and $\unicode[STIX]{x1D708}$ the kinematic viscosity of the fluid. This Reynolds number is understood here as a parameter whereas for genuine spatially developing flows it increases with the distance downstream and should rather be interpreted as a streamwise coordinate. The onset of TS waves is reported from linear stability analysis at $Re_{\unicode[STIX]{x1D6FF}^{\ast }}\approx 520$ (Jordison 1970), with moderate corrections when taking non-parallelism into account (Berlin 1998). However bypass transition has been reported for finite-amplitude disturbances introduced as early as $Re_{\unicode[STIX]{x1D6FF}^{\ast }}\approx 300$ (Jacobs & Durbin 2001). This suggests that bypass transition falls a priori into the category of subcritical flows such as pipe flow, channel flow, plane Couette flow or the asymptotic suction boundary layer (ASBL), for which turbulence can exist despite the linear stability of the base flow (Eckhardt et al. 2007). Borrowing concepts and toolboxes from subcritical transition hence appears as promising. In particular, a useful and recent nonlinear concept developed in wall-bounded flows is the edge manifold, the state space boundary that separates the respective basins of attraction of the laminar and the turbulent state (assumed here an attractor). This manifold $\unicode[STIX]{x1D6F4}$ is formally of codimension one, and is the stable manifold of a simpler low-dimensional invariant regime called the edge state. In principle all trajectories starting on $\unicode[STIX]{x1D6F4}$ reach asymptotically this edge state, while any small perturbation to such an edge trajectory is on the verge of both relaminarising and transitioning (Schneider & Eckhardt 2006). This property makes edge states, together with the set of all edge trajectories, the alternative nonlinear base flow of interest to predict bypass transition and therefore to prevent or to control it. Besides, edge states and their approximations display robust localisation in physical space as soon as the numerical domain allows for it (Duguet, Schlatter & Henningson 2009; Khapko et al. 2016). Recently, supported by numerical evidence that even arbitrary initial noise allows for a transient detection of the edge regime (Khapko et al. 2016), the concept of the edge state and its instability have been used to model the nucleation process of turbulent spots in the presence of strong free-stream turbulence (Kreilos et al. 2016). It is thus essential for the understanding of spot nucleation to gather additional knowledge about edge states in spatially developing boundary layer flows.

Former attempts at identifying edge states in boundary layer flows fall into two categories. The first category relies on different parallel approximations of the Blasius profile (Biau 2012; Khapko et al. 2013, 2014, 2016). Robust slow dynamics emerges in all cases: low- and high-speed streaks burst quasi-periodically before they switch location. The second category is free from the parallel assumption and gathers edge computations in longer domains (Cherubini et al. 2011; Duguet et al. 2012), although with some limitations, being the domain size or spanwise symmetry. The coherent structure emerging in all cases consists of an elongated localised pair of streaks, whose three spatial dimensions grow self-similarly with the variable $\unicode[STIX]{x1D6FF}^{\ast }$ . Being computationally much more demanding, this approach shows severe limitations in time. The longest bisection to date, originally performed with imposed spanwise symmetry, suggests that the symmetric edge trajectory also features quasi-cyclic regeneration of the low-speed streak in terms of rescaled time variables. Importantly, none of these investigations has reported the presence of TS waves, though such waves are expected to grow exponentially in time provided $Re_{\unicode[STIX]{x1D6FF}}^{\ast }$ can get large enough.

The present study aims at revisiting edge computations in the Blasius boundary layer, by extending their time horizon sufficiently far that the nonlinear interaction between TS waves and the streaks manifests itself. An earlier study of the interaction between streaks and TS waves by superimposing both disturbances can be found in Schlatter et al. (2010). In our case, the interaction is revealed by combining the classical bisection method with a numerical technique tailored to deal with moving localised disturbances. As we shall see, this calls for a generalisation of the concept of edge state to cases where the governing base flow instability is supercritical rather than subcritical. This situation is not specific to the Blasius profile and is expected in other shear flows such as channel flow, which also loses its stability at a finite value of the Reynolds number (Zammert & Eckhardt 2017a , b ).

The structure of the paper is at follows. Section 2 is devoted to the numerical aspects of this work, whereas the main results are given in § 3. Eventually, all the results and their implications are discussed in § 4.

2 Computational methodology

2.1 Flow set-up

The Blasius boundary layer is the incompressible boundary layer flow over a flat plate with zero pressure gradient. Let $x,y,z$ denote respectively the streamwise, wall-normal and spanwise coordinates, $x$ being measured from the inlet of the domain, itself located at a distance $x_{0}$ from the leading edge of the plate. Let $\boldsymbol{v}=(v_{x},v_{y},v_{z})$ be the corresponding velocity field. Let $\boldsymbol{u}=\boldsymbol{v}-\boldsymbol{v}^{B}$ be the velocity perturbation to the Blasius solution $\boldsymbol{v}^{B}$ , itself a self-similar two-dimensional solution of the incompressible boundary layer equations for the rescaled variable $\unicode[STIX]{x1D702}=y(U_{\infty }x)/\unicode[STIX]{x1D708}$ . A local Reynolds number can be defined as $Re_{\unicode[STIX]{x1D6FF}^{\ast }}=U_{\infty }\unicode[STIX]{x1D6FF}^{\ast }(x)/\unicode[STIX]{x1D708}$ , where $\unicode[STIX]{x1D6FF}^{\ast }=\int _{0}^{\infty }(1-v_{x}^{B})/U_{\infty }\,\text{d}y$ is the thickness of the undisturbed Blasius flow.

Direct numerical simulation of the incompressible Navier–Stokes equations is performed using SIMSON (Chevalier et al. 2007). The equations are solved in the velocity–vorticity formulation using a pseudo-spectral method, and are advanced in time using an explicit low storage fourth-order Runge–Kutta method for the nonlinear term and a second-order implicit Crank–Nicolson method for the linear terms. The velocity field $\boldsymbol{v}(x,y,z,t)$ is expanded in $N_{x}$ and $N_{z}$ Fourier modes in the $x$ and $z$ directions, respectively. The wall-normal expansion is based on $N_{y}$ Chebyshev modes. To satisfy the periodicity required by the Fourier expansion in the streamwise direction despite the spatial development, a fringe region is considered at the downstream end of the computational domain. Within the fringe an artificial volume force damps all velocity disturbances and re-establishes the correct inflow profile (see appendix A for details). The boundary layer thickness at the inlet of the domain $\unicode[STIX]{x1D6FF}_{0}^{\ast }$ is chosen so that $Re_{\unicode[STIX]{x1D6FF}_{0}^{\ast }}=300$ . Velocities and lengths are made non-dimensional by using respectively the quantities $U_{\infty }$ and $\unicode[STIX]{x1D6FF}_{0}^{\ast }$ . In these units the computational domain is a rectangular box $V$ of size $(L_{x},L_{y},L_{z})=(6000,60,100)$ . The spectral resolution is $(N_{x},N_{y},N_{z})=(6144,201,256)$ , excluding the additional modes in $x$ and $z$ used for dealiasing by the $3/2$ -rule. A no-slip/no-penetration boundary condition is enforced at the wall ( $y=0$ ) and a Neumann condition is used in the free stream ( $y=L_{y}$ ) for the three velocity components. The local resolution is comparable to the one used in Duguet et al. (2012).

2.2 The moving box technique

Numerical simulation of localised coherent structures usually requires computational domains at least one order of magnitude larger than their typical size (Duguet et al. 2009). When such coherent structures are allowed to move within the domain, length requirements become even more severe, while most of the flow domain contains only laminar, largely undisturbed flow.

Figure 1. Sketch of the flow geometry including the moving box.

In order to circumvent this issue, we suggest a moving box technique, not used previously in the case of spatially developing flows, in which the computational domain moves in the streamwise direction with piecewise constant velocity $c_{box}$ (see figure 1). The method then consists of Galilean changes of the reference frame of the type $v_{x}\leftarrow v_{x}-c_{box}$ , where $c_{box}$ is the instantaneous speed of the moving box (see appendix A for details). In the present simulations, the box velocity takes alternatively the values $c_{box}=0$ or $0.8$ depending on the time elapsed. Different histories of $c_{box}$ yield exactly the same results as long as all travelling coherent structures remain within the computational domain. The present computation of the edge trajectory using the moving box technique was validated against the case of a non-moving domain of size $(L_{x},L_{y},L_{z})=(12\,000,60,100)$ , simulated with a spectral resolution of $(N_{x},N_{y},N_{z})=(12\,288,201,256)$ .

2.3 Edge tracking algorithm

The dynamics on the edge is tracked iteratively using the bisection technique originally introduced by Itano & Toh (2001) and Skufca, Yorke & Eckhardt (2006). A single scalar time-dependent observable $a(t)$ is selected, together with two threshold values $a_{L}$ and $a_{T}$ allowing to discard unambiguously trajectories evolving towards the laminar or the turbulent regime, respectively. The choice of the observable is here

(2.1) $$\begin{eqnarray}a=\left(\frac{1}{\text{vol}(V)}\int _{V}|\unicode[STIX]{x1D714}_{x}|^{2}\,\text{d}v\right)^{1/2},\end{eqnarray}$$

a normalised integral measure of the streamwise vorticity $\unicode[STIX]{x1D714}_{x}$ in the computational domain of volume $\text{vol}(V)$ , where $\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D714}_{x},\unicode[STIX]{x1D714}_{y},\unicode[STIX]{x1D714}_{z})$ is the vorticity field. This observable is not affected by the choice of the reference frame. It is a measure of the amplitude of streamwise vortices, known to play a crucial role in the streak formation. It is zero for the laminar Blasius profile. Moreover the observable is also zero when the aforementioned profile is perturbed by TS waves in their linear stage (with no spanwise dependence and zero spanwise velocity). For a trajectory starting from a given initial perturbation velocity field $\boldsymbol{u}_{0}$ , localised in space, the observable $a(t)$ is monitored until its value crosses one of the two thresholds, $a_{L}$ from above, or $a_{T}$ from below. In the present study we used $a_{L}=8.74\times 10^{-5}$ and $a_{T}=2.68\times 10^{-3}$ , although different values of $a_{L}$ have been tested without significant influence. Depending on which threshold is crossed, the field $\boldsymbol{u}_{0}$ is re-labelled respectively $\boldsymbol{u}_{L}$ and $\boldsymbol{u}_{T}$ . From the history of the bisection and the knowledge of the last fields $\boldsymbol{u}_{L}$ and $\boldsymbol{u}_{T}$ , a new update $\boldsymbol{u}_{0}\leftarrow \frac{1}{2}(\boldsymbol{u}_{L}+\boldsymbol{u}_{T})$ is constructed (for simplicity one can start with $\boldsymbol{u}_{L}=\mathbf{0}$ ). This is a convergent iterative process that brackets the edge manifold. The iterative bisection process yields a sequence of positive numbers $\unicode[STIX]{x1D706}^{(n)}$ , $n=1,2,\ldots$ and an associated sequence of initial conditions $\boldsymbol{u}^{(n)}$ . For $n$ large enough $\unicode[STIX]{x1D706}^{(n)}\approx \unicode[STIX]{x1D706}^{\ast }$ , such that in the $n\rightarrow \infty$ limit $\boldsymbol{u}_{\unicode[STIX]{x1D706}^{\ast }}=\unicode[STIX]{x1D706}^{\ast }\boldsymbol{u}_{0}^{(0)}$ lies exactly on the edge manifold $\unicode[STIX]{x1D6F4}$ , and $\boldsymbol{u}_{\unicode[STIX]{x1D706}^{(n)}}$ is an approximation of such a state. In practice, when machine precision is met regarding the accuracy of $\unicode[STIX]{x1D706}^{\ast }$ , the whole bisection process is restarted from another state further along the edge trajectory. This bisection technique is similar to that used in previous subcritical systems, except that crossing the bound $a_{L}$ is here not necessarily interpreted as relaminarisation. In Cherubini et al. (2011) and Duguet et al. (2012) the time horizon in the bisection was not sufficiently long for the linear instability of the base flow to manifest itself, therefore crossing the threshold $a_{L}$ from above necessarily corresponded to an effective relaminarisation. We demonstrate here that this picture needs be revised when bisecting is no longer limited in time.

3 Results

3.1 The edge for moderate times

The initial condition $\boldsymbol{u}_{0}$ consists of a pair of counter-rotating vortices located close to the wall at $x=50$ , $z=0$ , and not aligned with the main flow direction; unlike in Duguet et al. (2012) the evolution of the flow is here free from any discrete symmetry constraint. For $0\leqslant t\leqslant 4700$ (in units of $\unicode[STIX]{x1D6FF}_{0}^{\ast }/U_{\infty }$ ) it is possible to effectively bracket the edge using the algorithm described above and the current definition of the observable $a(t)$ . Time series of $a(t)$ during the bisection process are shown in figure 2. The dynamics along the edge trajectory is globally comparable to that in Khapko et al. (2016) though there are also clear differences. The velocity field along the edge trajectory is characterised by long streamwise streaks and streamwise vortices, both characteristics of a self-sustained process described in Duguet et al. (2012). It remains spatially localised in three dimensions at all times. Three-dimensional visualisations of three specific snapshots are shown in figure 3, where a typical regeneration cycle of the edge dynamics is displayed. The whole localised structure is unsteady; similarly to Duguet et al. (2012) it is possible to identify a robust streaky core upstream, together with recurrent secondary structures that detach from the main core, travel further downstream and dissipate. The global structure and the temporal dynamics along this edge trajectory is qualitatively closer to that computed in other flows, notably the localised edge state obtained in the ASBL (Khapko et al. 2016) except for the sparser occurrence of bursts. Compared to the spanwise-symmetric analysed in Duguet et al. (2012), the streak organisation is different. The robust pair of symmetric low-speed streaks flanking a high-speed streak is lost in favour of an unsteady dynamics, where the sandwiched streak is alternatively the high-speed and the low-speed one. As this process does not occur simultaneously at all locations, it is possible to visualise different stages of the streak switching cycle at the same time by comparing for instance the initial and final frames in figure 4. Visualisation of the $\unicode[STIX]{x1D706}_{2}$ criterion (Jeong & Hussain 1995) in figure 3 shows the presence of slightly asymmetric vortical structures above the regions where the streaks pinch together. Our observations suggest that the strongest of these approximately hairpin-shaped vortices, because of the wall-normal flows they induce, act as precursors of the switch events such as the one in figure 4. Streak switching appears recurrently in the literature as characteristic of unsteady edge states (Toh & Itano 2003; Biau 2012; Khapko et al. 2013, 2016), and could now be confirmed in spatially developing boundary layer flows as well.

Figure 2. Observable $a(t)$ versus $t$ during application of the bisection algorithm over moderate time horizons. The thick line (blue online) represents the edge trajectory found by bisection, the two dotted lines are the observable bounds $a=a_{L}=8.74\times 10^{-5}$ and $a_{T}=2.68\times 10^{-3}$ (see text).

Figure 3. Three-dimensional perspective view from above of the edge trajectory at times $t=1700,2750$ and $3550$ (from top to bottom), showing isosurfaces of streamwise perturbation velocity with respect to spanwise mean with values $0.06$ and $-0.08$ (red and blue, respectively), together with vortical structures $\unicode[STIX]{x1D706}_{2}=-1.5\times 10^{-5}$ (green). Flow from left to right. The black lines are separated by a distance of 200 in units of $\unicode[STIX]{x1D6FF}_{0}^{\ast }$ . The last snapshot includes the cross-sections presented in figure 4.

Figure 4. Illustration of a streak switching event along the edge trajectory at $t=3550$ ; $(y,z)$ cross-sections of $u_{x}$ . From left to right $x=2140,2420$ and $2980$ .

The spatially developing nature of the boundary layer suggests the use of a ‘local’ boundary layer thickness as a characteristic length scale. To define a length scale rigorously it is necessary to refer to an unambiguous streamwise location. Exploiting the spatial localisation of the coherent structure under study, we define its streamwise position $x_{G}$ in analogy with the centre of mass by

(3.1) $$\begin{eqnarray}x_{G}(t)=\frac{\displaystyle \iint x|u_{x}(x,y_{p},z,t)|\,\text{d}x\,\text{d}z}{\displaystyle \iint |u_{x}(x,y_{p},z,t)|\,\text{d}x\,\text{d}z},\end{eqnarray}$$

where $y_{p}$ is the wall-normal coordinate of a given plane known to intersect the coherent structure at all times (in practice $y_{p}=1.5\unicode[STIX]{x1D6FF}_{0}^{\ast }$ ). The weighting function in the integrals in (3.1) is the absolute value of the streamwise perturbation velocity $u_{x}$ , $x_{G}(t)$ is an estimation of the displacement of the main coherent structure with time from its initial location. The local displacement thickness $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FF}^{\ast }(x_{G})$ , evaluated at this time-dependent location, provides the scaling factor $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{0}^{\ast }$ for the boundary layer growth.

Figure 5. The $(z,t)$ space–time diagram of $u_{x}(x,y,z,t)$ along the edge trajectory, in a frame moving with the centre of mass located at $x=x_{G}(t)$ and $y=y_{p}$ , of the perturbation streamwise velocity. Both $z$ and $t$ are rescaled by the local boundary layer thickness $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{0}^{\ast }$ , causing the apparent narrowing of the domain.

Monitoring $x_{G}(t)$ and its temporal derivative ${\dot{x}}_{G}(t)$ is a convenient way to define the advection speed of the structure as a whole. After a relatively short initial transient of $t\approx 500$ , ${\dot{x}}_{G}(t)$ reaches an almost constant value of 0.625 in units of $U_{\infty }$ , comparable to the values found in Duguet et al. (2012). Velocities scale with $U_{\infty }$ and are unaffected by the rescaling, while vorticities and the domain dimensions are now time dependent. The rescaled space–time diagram for the streamwise velocity perturbation $u_{x}(x=x_{G}(t),y=y_{p},z,t)$ (evaluated at the centre of mass of the coherent structure) is shown in figure 5. The approximately constant width of the streak core in time is a confirmation that the $\unicode[STIX]{x1D6FF}^{\ast }(x_{G})/\unicode[STIX]{x1D6FF}_{0}^{\ast }$ -scaling is relevant at least for the spatial coordinates. Similarly, the same scaling for the time coordinate suggests that the apparent period for the streak switching becomes inherently larger as the boundary layer grows. This diverging time scale, together with the diverging length of the streaks, is the main quantitative obstacle to significantly longer edge tracking. In practice it is possible to observe only two full streak switching cycles, during the time intervals $t\approx$ 600–1300 and 1900–3300, i.e. in the rescaled variables $t\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF}\approx$ 300–420 and 630–770. Furthermore the whole state grows in length, from ${\approx}400\unicode[STIX]{x1D6FF}_{0}^{\ast }$ at $t=1750$ to ${\approx}1200\unicode[STIX]{x1D6FF}_{0}^{\ast }$ at $t=4700$ , as shown in figure 3.

3.2 Coexistence of Tollmien–Schlichting waves and streaks for long times

3.2.1 Physical space

The edge trajectory is by construction linearly unstable. So is the laminar base flow, except that the instabilities manifest themselves over different time scales, with the base flow instability being typically slower (Zammert & Eckhardt 2017b ). If the time horizon in the bisection is large enough there is a competition between bypass transition (expected from the linear instability of the edge) and the classical transition scenario (featuring Tollmien–Schlichting waves, due to the linear instability of the base flow). We report now the unforced simultaneous occurrence of these two transition mechanisms.

Figure 6. The $(z,t)$ space–time diagrams of $u_{x}(x,y,z,t)$ evaluated at $x=x_{G}(t)$ and $y=y_{p}$ , for two initially nearby trajectories bracketing the edge trajectory. Both $z$ and $t$ are rescaled by the local boundary layer thickness $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{0}^{\ast }$ . (a) Bypass transition route. (b) Streak decay followed by classical transition based on the growth of the Tollmien–Schlichting waves.

Figure 6 shows that a relative difference $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D706}^{\ast }$ of magnitude $O(10^{-4})$ in the choice of the bisection parameter $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}^{\ast }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}$ leads to different transition scenarios in the long-time horizon. For $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}>0$ destabilisation of the edge trajectory leads to breakdown of the streaks akin to that reported in bypass transition (Brandt et al. 2004). The breakdown of the streak is local in $z$ , as visible in figure 6(a). On the other hand, for $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D706}<0$ the streaky structure is not strong enough to self-sustain, instead it starts to decay viscously. During the streak decay the linear instability of the flow manifests itself in the form of TS waves emerging upstream of the streak core, well visible in figure 6(b). The apparent wavelength as well as the propagation speed of these waves makes it clear that they correspond to TS waves. Regarding wavelengths, these waves have $k_{z}=0$ and $k_{x}<0.35$ at all times, which matches the linear stability curve from the literature (cf. Berlin 1998). The speed of propagation of the waves can be measured from space–time diagrams such as figure 7. It is approximatively $0.32\pm 0.005$ , matching quantitatively the speed of the TS waves obtained from linear stability analysis of the frozen Blasius profile at their onset (Schmid & Henningson 2001). A more detailed study of the TS waves is shown in figure 8(b). The figure shows the comparison between the nonlinear TS waves in our simulation and the results from linear stability analysis, for a streamwise wavenumber $k_{x}\approx 2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{TS}=0.037$ , where $\unicode[STIX]{x1D706}_{TS}=170\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF}_{x,TS}^{\ast }$ . The wavelength $\unicode[STIX]{x1D706}_{TS}$ is obtained from the intercrest distance 170 scaled with the local boundary thickness $\unicode[STIX]{x1D6FF}_{x,TS}^{\ast }$ , itself an evaluation of $\unicode[STIX]{x1D6FF}^{\ast }$ at the position $x=x_{TS}$ . The position $x_{TS}$ is obtained by tracking in time and space the specific crest passing through $x=5186$ for time $t=8500$ (indicated in figure 7 by a cross). Moreover figure 8 (a) shows the position of $x_{G}(t)$ from the trajectory bracketing the edge from below. The position of the edge is well fitted by the straight line $x_{G}(t)=x_{G0}+v_{Gx}t$ , where $x_{G0}=70$ and $v_{Gx}=0.625$ . It is clear that the TS waves travel at a velocity approximately half that of the streak core. This slower propagation speed makes the associated TS wavepacket detach progressively from the decaying streak core, while as a result the position of the centre of mass drifts upstream of the streak core as seen in figure 6(b). The exact way how the TS wavepacket is triggered is not fully understood, preliminary examination points towards disturbances growing from the wake of the streaks. The waves grow upstream of the coherent structure, and they do not interact with the streak core. This is consistent with previous observations that finite-amplitude streaks act as additional damping of TS-wave growth (Cossu & Brandt 2004) and that streaks and TS waves are not seen to overlap spatially. The TS wavepacket grows in size as it travels downstream, with the number of individual rollers growing too (Gaster 1975). At later times the waves inside the wavepacket undergo a secondary instability with spanwise wavenumber $k_{z}\neq 0$ (of fundamental Klebanoff type), followed by a rapid local breakdown into turbulence.

Figure 7. Space–time diagram of $u_{x}(x,y,z,t)$ for $z=-10$ and $y=y_{p}$ for the trajectory away from the edge manifold, as displayed in figure 6(b). The streaks decay while the TS wavepacket grows in amplitude, forming a turbulent spot. The red cross indicates the initial position of the crest tracked in figure 8.

Figure 8. Streamwise propagation velocities. (a) Position of the centre of mass of the coherent structure $x_{G}(t)$ for a trajectory below the edge versus time. The vertical line indicates the time limitation of edge tracking (see text). Circles: data, solid line (blue online): linear fit. (b) Phase velocity of TS waves, numerical simulation versus linear instability analysis (see text). The error bars represent the accuracy in estimating the convection speed from figure 7.

3.2.2 State space

Further complementary information about the global temporal dynamics can be obtained from phase portraits. We use here the $L^{2}$ -norms of the three integral observables $\unicode[STIX]{x1D714}_{x}$ , $\unicode[STIX]{x1D714}_{y}$ and $u_{z}$ . In order to take into account the evolution of the boundary layer thickness downstream and the variable size of the domain when expressed in units of $\unicode[STIX]{x1D6FF}(t)$ rather than $\unicode[STIX]{x1D6FF}_{0}^{\ast }$ , the root mean square (r.m.s.) vorticities need to be rescaled using the correcting factor $(\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF})^{1/2}$ , and the r.m.s. velocities by $(\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF})^{(3/2)}$ , as explained in Duguet et al. (2012). We hence consider the three phase portraits observables $\unicode[STIX]{x1D6FA}_{x}$ , $\unicode[STIX]{x1D6FA}_{y}$ and $W$ , defined by (3.2)–(3.4) and parametrised by time only:

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{x}=(\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF})^{1/2}\left(\frac{1}{\text{vol}(V)}\int _{V}|\unicode[STIX]{x1D714}_{x}|^{2}\,\text{d}v\right)^{1/2}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{y}=(\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF})^{1/2}\left(\frac{1}{\text{vol}(V)}\int _{V}|\unicode[STIX]{x1D714}_{y}|^{2}\,\text{d}v\right)^{1/2}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle W=(\unicode[STIX]{x1D6FF}_{0}^{\ast }/\unicode[STIX]{x1D6FF})^{3/2}\left(\frac{1}{\text{vol}(V)}\int _{V}|u_{z}|^{2}\,\text{d}v\right)^{1/2}. & \displaystyle\end{eqnarray}$$

The reduced phase portrait is displayed in figure 9 (a). The initial condition associated with $\unicode[STIX]{x1D706}^{\ast }$ corresponds to the point $(\unicode[STIX]{x1D6FA}_{x},\unicode[STIX]{x1D6FA}_{y},W)=(1.72,6.62,1.80)\times 10^{-4}$ . Trajectories bracketing the edge have an initial condition very close to it, and are indistinguishable in the figure. Past a initial transient of approximatively $500\unicode[STIX]{x1D6FF}_{0}^{\ast }/U_{\infty }$ , all trajectories approach a single recurrent region bounded from above by $W=6\times 10^{-5}$ and characterised by apparent swirls in all displayed variables. At later times the trajectories leave this part of state space and evolve towards higher $W$ values associated with turbulent flow. The recurrence region corresponds to the neighbourhood of the unstable edge state. It is displayed in the zoom in figure 10, where up to three recurrent cycles can be distinguished. Each of the swirls in the ( $\unicode[STIX]{x1D6FA}_{x},\unicode[STIX]{x1D6FA}_{y}$ ) plane corresponds in physical space to a streak switching event and a release of vortical perturbations downstream. The dynamics of the trajectories leaving the edge is interesting, because independently of the route to turbulence all trajectories point towards the same direction, i.e. towards a common turbulent attractor. This last aspect, sets this flow apart from the usual subcritical picture characterised by the coexistence of two different attractors with clearly distinct basins of attraction (see e.g. the introductory sketch in Duguet et al. (2013)).

Figure 9. Three-dimensional phase portrait using global variables $\unicode[STIX]{x1D6FA}_{x}$ , $\unicode[STIX]{x1D6FA}_{y}$ and $W$ . The initial condition associated lies at approximately $(\unicode[STIX]{x1D6FA}_{x},\unicode[STIX]{x1D6FA}_{y},W)=(1.72,6.62,1.80)\times 10^{-4}$ but the first 500 time units are not shown. The blue dot at the $(0,0,0)$ is the laminar state. The red trajectory corresponds to bypass transition, the blue one to classical TS transition, and their points in common define the converged part of the edge trajectory. In grey, the projections on the different two-dimensional planes.

Figure 10. Zoom on the edge region using only $\unicode[STIX]{x1D6FA}_{x}$ and $\unicode[STIX]{x1D6FA}_{y}$ (same data as figure 9). The dotted lines correspond to two other bracketing trajectories. The first 300 time units are omitted. Swirls are indicated using closed dashed lines (orange online).

4 Discussion and conclusions

Large-scale computational edge tracking in a spatially developing Blasius boundary layer has been revisited, considering much longer time horizons than previously. In addition to a better understanding of the short-time dynamics, the results call for a revision of the concept of edge at large times, in the special case where the base flow instability is of a supercritical nature rather than subcritical. On the shorter time scales, a regeneration process is confirmed, whose characteristic time scale increases as the boundary layer thickness grows. Before the edge instability, recurrent visits to a streaky active core are observed as in most previous studies (Duguet et al. 2012; Khapko et al. 2013, 2016), including the slow phenomenon of streak switching, observed here in a spatially developing configuration. Furthermore, concomitant scenarios of bypass and classical transition have been observed on longer time scales, and their possibly simultaneous occurrence blurs the long-time output of the edge tracking algorithm. To our knowledge, whereas instances of such a coexistence have been reported in more complicated geometries (Xu et al. 2017; Canton et al. 2019) this coexistence is reported for the first time that in a simulation of an unforced boundary layer flow.

Figure 11. Observable $a(t)$ versus $t$ during application of the bisection algorithm. The green vertical line marks the maximum time for the edge tracking algorithm using that observable, and the two dotted lines stand for the observable bounds $a=a_{L}$ and $a_{T}$ (see text).

A moving box technique allows for a more efficient usage of the computational domain. This moving box technique is a simple generalisation of the Galilean changes of reference frequently used in simulations of parallel flows. It proves efficient for the tracking of all spatially localised coherent structures; it is not limited to edge tracking and can also handle minimal seeds, linear or nonlinear localised wavepackets and incipient turbulent spots. The longest individual simulations requested more than 12 000 time units for one of the two bounds $a_{L}$ or $a_{T}$ to be reached, this is at least four to six times more than in Duguet et al. (2012) due to the use of the moving box technique. If both the streaky state and the TS waves, which travel at markedly different speeds, are contained in the same computational domain (the worst-case scenario), the improvement in tracking time is 40 % compared to the fixed case, and much more otherwise. In total, edge tracking has been achieved here over a time horizon three times longer than in previous investigations of the Blasius boundary layer. The temporal limit for the edge tracking is met when the relative difference in the main observable between the closest edge-bracketing trajectories is 2 % or less, and we have no further information on the nature (bypass or classical) of the trajectories within these 2 %. This limit corresponds here $t\approx 4700$ and the corresponding distance, expressed in non-dimensional units of $Re_{x}=U_{\infty }x/\unicode[STIX]{x1D708}$ , is $Re_{x}^{max}\approx 9\times 10^{5}$ .

From the perspective of dynamical system theories, the reinterpretation of the edge manifold $\unicode[STIX]{x1D6F4}$ as a boundary between two transition scenarios leads to interesting theoretical questions. One of the objectives of the description of the transition in terms of dynamical systems is to map out the state space of the system together with its invariant sets, at least in the context of an initial-value problem. Note that the Reynolds number is here not a parameter. The only fixed point of the system is the laminar state, which happens to be an unstable one, in the sense that there always exist perturbations likely to grow exponentially in time (without any parameter threshold to exceed). The turbulent attractor may not be within reach or even well defined, but the edge state as relative attractor is well defined in the asymptotic limit $t\rightarrow +\infty$ . In particular it is an invariant set for the temporal dynamics. It is also a self-similar object in terms of spatial dimensions and temporal dynamics. However the edge tracking technique, which by construction only identifies finite-time trajectories belonging to the stable manifold of the edge state, finds here a natural time limitation. This algorithm is based on bounds for a given observable $a$ , originally designed to discard the basins of attraction of the other attracting sets. Here there is only one attracting turbulent set, therefore $\unicode[STIX]{x1D6F4}$ does not separate two distinct basins of attraction like in bistable systems. This is the main difficulty in interpreting phase portraits such as those in figure 9. In particular, no matter how small or how large the bound $a_{L}$ , provided the computational box is long enough there will always be another trajectory with $a<a_{L}$ also experiencing transition, possibly featuring destabilisation of TS waves. A similar situation was described by Zammert & Eckhardt (2017b ) in the context of plane Poiseuille flow in a small periodic domain, where the type of transition was determined from the knowledge of pre-defined transition times. In the present study, a suitable choice of the observable can extend the time over which the bisection stays valid. For instance, the use of the $\unicode[STIX]{x1D714}_{x}$ component to define the observable $a$ in (2.1) is a convenient way to ignore Tollmien–Schlichting waves in their initial stage of their growth, because for TS waves it is non-zero only after they have undergone a secondary instability. Other one-dimensional observables can be envisioned. However, without extensive monitoring of the properties of the trajectories generated during the bisection, a one-dimensional observable by itself does not contain enough information on the route to turbulence (bypass or classical) that is followed. Multi-dimensional observables, containing e.g. information about amplitude or growth of the TS waves, are a possible alternative but they do not make bisection generically possible, since bisection is essentially a one-dimensional search process. For investigations where the asymptotic nature of the edge state matters, alternative methods are now welcome, specifically local methods that do not require the knowledge of how trajectories behave far away in state space from the edge manifold under study.

Figure 12. Value of $Tu$ (%) versus  $Re_{x}$ for experimental bypass transition data (Shahinfar & Fransson 2011), showing intermittency of 10 % and 90 % (thick lines) and 50 % (dashed line). The arrows mark the $Re_{x}$ -limitation of the bisection algorithm. For larger times and $Re_{x}$ , mixed transition is expected (grey area). The theoretical onset for TS waves lies at $Re_{x}=9.1\times 10^{4}$ .

The possibility for proper coexistence of the two routes to turbulence sets a natural time limitation to the current algorithm of edge tracking. This time limit, estimated in § 3 as approximately $4700\unicode[STIX]{x1D6FF}_{0}^{\ast }/U_{\infty }$ , is highlighted in figure 11, which shows the time series of the observable $a(t)$ during the bisection process. It corresponds to a limitation of the maximum distance downstream along which edge tracking can be achieved, estimated here as $Re_{x}^{max}\approx 9\times 10^{5}$ . This is close to 10 times the value of $Re_{x}=9.1\times 10^{4}$ corresponding to the onset of TS waves according to linear stability analysis, and much larger than the value of $Re_{x}=3\times 10^{4}$ corresponding to the inlet defined by $Re_{\unicode[STIX]{x1D6FF}_{0}^{\ast }}=300$ . We now demonstrate that if the asymptotic edge state remains out of reach, the computed edge trajectories are relevant in practice for most instances of bypass transition under the influence of free-stream turbulence. We use for this purpose the experimental database by Shahinfar & Fransson (2011) from the KTH wind tunnel with parametrisable free-stream turbulence, measured by $Tu=u_{rms}/U_{\infty }$ and expressed in %. The $Re_{x}$ -intervals where intermittency (and hence nucleation of turbulent spots) occurs are reported in figure 12 as functions of the parameter $Tu$ , where the dependence on the integral length scale of the incoming turbulence is already taken into account.

The lower and upper bounds of these intervals are defined respectively by the values of $Re_{x}$ for which the intermittency factor is 0.1 and 0.9. They both scale like $Re_{x}=O(Tu^{-2})$ (Shahinfar & Fransson 2011). For $Tu\approx 2\,\%$ and above, the corresponding values of $Re_{x}$ in the diagram are always below the value of $Re_{x}^{max}$ . This means that the edge state is known – and already computed – for these cases, and that it can be used as an alternative base flow for a stability analysis (generalised to unsteady flows). The only difficult situation corresponds to $Tu<2\,\%$ , when the bisection technique is not able to track the edge all the way into the zone where spots nucleate. In this mixed area, both routes to turbulence are indistinguishable, at least from the monitoring of one scalar observable only. This mixed area is the parameter area where the turbulent patches observed in practice can be due to the destabilisation of either TS waves or streaky perturbations. It has long been known as a delicate range of parameters and has sometimes been labelled ‘weak bypass regime’, as opposed to the ‘strong bypass regime’ where streak breakdown is the sole cause for transition (Narasimha 1994). Future efforts should be made to better characterise the complex transition to turbulence in this intermediate regime.