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 Reference Schubauer and Skramstad1948; Boiko et al. Reference Boiko, Westin, Klingmann, Kozlov and Alfredsson1994; Schmid & Henningson Reference Schmid and Henningson2001; Schlichting & Gersten Reference Schlichting and Gersten2016). 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 Reference Klebanoff1955; Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967; Morkovin Reference Morkovin1969). Those streaks support further instabilities (Henningson, Lundbladh & Johansson Reference Henningson, Lundbladh and Johansson1993; Kachanov Reference Kachanov1994; Andersson et al. Reference Andersson, Brandt, Bottaro and Henningson2001; Jacobs & Durbin Reference Jacobs and Durbin2001; Brandt, Schlatter & Henningson Reference Brandt, Schlatter and Henningson2004; Schlatter et al. Reference Schlatter, Brandt, Lange and Henningson2008; Vaughan & Zaki Reference Vaughan and Zaki2011; Hack & Zaki Reference Hack and Zaki2014) leading to the nucleation of turbulent spots that invade the flow (Matsubara & Alfredsson Reference Matsubara and Alfredsson2001; Brandt et al. Reference Brandt, Schlatter and Henningson2004; Kreilos et al. Reference Kreilos, Khapko, Schlatter, Duguet, Henningson and Eckhardt2016). 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 Reference Jordison1970), with moderate corrections when taking non-parallelism into account (Berlin Reference Berlin1998). However bypass transition has been reported for finite-amplitude disturbances introduced as early as
$Re_{\unicode[STIX]{x1D6FF}^{\ast }}\approx 300$
(Jacobs & Durbin Reference Jacobs and Durbin2001). 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.
Reference Eckhardt, Schneider, Hof and Westerweel2007). 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 Reference Schneider and Eckhardt2006). 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 Reference Duguet, Schlatter and Henningson2009; Khapko et al.
Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016). Recently, supported by numerical evidence that even arbitrary initial noise allows for a transient detection of the edge regime (Khapko et al.
Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016), 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.
Reference Kreilos, Khapko, Schlatter, Duguet, Henningson and Eckhardt2016). 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 Reference Biau2012; Khapko et al.
Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013, Reference Khapko, Duguet, Kreilos, Schlatter, Eckhardt and Henningson2014, Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016). 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.
Reference Cherubini, De Palma, Robinet and Bottaro2011; Duguet et al.
Reference Duguet, Schlatter, Henningson and Eckhardt2012), 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. (Reference Schlatter, Deusebio, de Lange and Brandt2010). 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 Reference Zammert and Eckhardt2017a ,Reference Zammert and Eckhardt 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.
Reference Chevalier, Schlatter, Lundbladh and Henningson2007). 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012).
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. Reference Duguet, Schlatter and Henningson2009). 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 (Reference Itano and Toh2001) and Skufca, Yorke & Eckhardt (Reference Skufca, Yorke and Eckhardt2006). 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
$$\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. (Reference Cherubini, De Palma, Robinet and Bottaro2011) and Duguet et al. (Reference Duguet, Schlatter, Henningson and Eckhardt2012) 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012) 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. (Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016) 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012). 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012) 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.
Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016) except for the sparser occurrence of bursts. Compared to the spanwise-symmetric analysed in Duguet et al. (Reference Duguet, Schlatter, Henningson and Eckhardt2012), 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 Reference Jeong and Hussain1995) 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 Reference Toh and Itano2003; Biau Reference Biau2012; Khapko et al.
Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013, Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016), 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
$$\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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012). 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 Reference Zammert and Eckhardt2017b ). 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.
Reference Brandt, Schlatter and Henningson2004). 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 Reference Berlin1998). 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 Reference Schmid and Henningson2001). 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 Reference Cossu and Brandt2004) 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 Reference Gaster1975). 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012). 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:
$$\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}$$
$$\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}$$
$$\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. (Reference Duguet, Monokrousos, Brandt and Henningson2013)).
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. Reference Duguet, Schlatter, Henningson and Eckhardt2012; Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013, Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016), 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. Reference Xu, Mughal, Gowree, Atkin and Sherwin2017; Canton et al. Reference Canton, Rinaldi, Örlü and Schlatter2019) 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012) 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 (Reference Zammert and Eckhardt2017b
) 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 Reference Shahinfar and Fransson2011), 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 (Reference Shahinfar and Fransson2011) 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 Reference Shahinfar and Fransson2011). 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 Reference Narasimha1994). Future efforts should be made to better characterise the complex transition to turbulence in this intermediate regime.
Acknowledgements
We wish to acknowledge the late B. Eckhardt’s contribution to our research over the years, in particular the insightful discussions regarding the edge states in spatially developing flows. Thanks also for Professor Eckhardt’s hospitality during numerous visits to Marburg, the last one occurring in January 2019 where valuable input on the current paper was given by him. Financial support by the Swedish Research Council (VR) grant no. 2016-03541 is gratefully acknowledged. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC and NSC.
Appendix A
In order to implement the moving box it is necessary to do a Galilean change of reference and to update the base flow enforcing correct inflow and outflow conditions while preserving streamwise periodicity. The update of the base flow can be explained as follows. We first consider the incompressible Navier–Stokes equations
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}v_{i}+v_{j}\unicode[STIX]{x2202}_{j}v_{i}=-\unicode[STIX]{x2202}_{i}p+Re_{\unicode[STIX]{x1D6FF}_{0}^{\ast }}^{-1}\unicode[STIX]{x2202}_{j}^{2}v_{i}+F_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{j}v_{j}=0. & \displaystyle\end{eqnarray}$$
The perturbation velocity field
$\boldsymbol{u}$
is deduced from the full velocity field
$\boldsymbol{v}$
by subtracting the analytical Blasius profile
$\boldsymbol{v}^{B}=(u_{x}^{B},v_{x}^{B},0)$
, itself an excellent approximation to the steady base flow for
$x\gg 1$
(Schlichting & Gersten Reference Schlichting and Gersten2016). The field
$\boldsymbol{F}=(F_{x},F_{y},F_{z})$
is a volume force imposed in the fringe region in order to ensure the
$x$
-periodicity in spatial simulations. The form of the forcing is
$$\begin{eqnarray}\boldsymbol{F}=\unicode[STIX]{x1D6FE}(x)(\boldsymbol{{\mathcal{U}}}(t)-\boldsymbol{v}),\end{eqnarray}$$
where
$\boldsymbol{v}$
is the instantaneous flow field and
$\unicode[STIX]{x1D6FE}(x)$
is a non-negative fringe function. The function
$\unicode[STIX]{x1D6FE}(x)$
is non-zero in the fringe region only. The streamwise component of
$\boldsymbol{{\mathcal{U}}}(t)$
is computed as (see Chevalier et al.
Reference Chevalier, Schlatter, Lundbladh and Henningson2007)
$$\begin{eqnarray}{\mathcal{U}}(x,y,z,t)=U(x,y,z,t)+[U(x+x_{L},y,z,t)-U(x,y,z,t)]S\left(\frac{x-x_{blend}}{\unicode[STIX]{x1D6E5}_{blend}}\right),\end{eqnarray}$$
where
$S(x_{blend},\unicode[STIX]{x1D6E5}_{blend})$
is a step function, and where
$x_{blend}$
and
$\unicode[STIX]{x1D6E5}_{blend}$
are chosen so that the flow connects smoothly the outflow to the inflow. The normal component of
$\boldsymbol{{\mathcal{U}}}(t)$
is computed from continuity. The novelty in our approach is that
$U(x,y,z,t)$
, the solution to the boundary layer equations, is allowed to change with time, thereby modifying the base flow on which the fringe acts;
$U(x,y,z,t)=v_{x}^{B}(x,y,x_{0}(t))$
is updated at time
$t$
, such that the inlet of the box
$x_{0}$
for which
$\unicode[STIX]{x1D6FF}_{0}^{\ast }=1$
at
$Re_{\unicode[STIX]{x1D6FF}^{\ast }}=300$
reads
$x_{0}\leftarrow x_{0}+c_{box}t_{m}$
, where
$t_{m}$
is the total time during which the box has moved with velocity
$c_{box}$
.
Open access
