The triple-deck stage of marginal separation.

The method of matched asymptotic expansions is applied to the investigation of transitional separation bubbles. The problem-specific Reynolds number is assumed to be large and acts as the primary perturbation parameter. Four subsequent stages can be identified as playing key roles in the characterization of the incipient laminar-turbulent transition process: due to the action of an adverse pressure gradient, a classical laminar boundary layer is forced to separate marginally (I). Taking into account viscous-inviscid interaction then enables the description of localized, predominantly steady, reverse flow regions (II). However, certain conditions (e.g. imposed perturbations) may lead to a finite-time breakdown of the underlying reduced set of equations. The ensuing consideration of even shorter spatio-temporal scales results in the flow being governed by another triple-deck interaction. This model is capable of both resolving the finite-time singularity and reproducing the spike formation (III) that, as known from experimental observations and direct numerical simulations, sets in prior to vortex shedding at the rear of the bubble. Usually, the triple-deck stage again terminates in the form of a finite-time blow-up. The study of this event gives rise to a noninteracting Euler-Prandtl stage (IV) associated with unsteady separation, where the vortex wind-up and shedding process takes place. The focus of the present paper lies on the triple-deck stage III and is twofold: firstly, a comprehensive numerical investigation based on a Chebyshev collocation method is presented. Secondly, a composite asymptotic model for the regularization of the ill-posed Cauchy problem is developed. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s10665-021-10125-3.

Introduction

The present paper is concerned with the investigation of the early stages of laminar–turbulent transition induced by localized boundary layer separation. Notably in this regard, the analysis will be approached from the perspective of high Reynolds number () asymptotics. The mentioned type of transition may be observed e.g. in the flow past the leading-edge suction side of an airfoil. Here, the angle of attack is kept almost fixed at a relatively small value (or may change in a quasi-steady fashion) and the main, predominantly inviscid, part of the flow remains steady and follows the streamlined body contour. This in various practical applications quite common setting is distinct from the more complex ‘dynamical stall’ flow phenomenon that manifests itself on an airfoil during rapid, transient motion in which the angle of attack exceeds the static stall limit value [1, 2]. Whereas the latter flow situation is characterized by the formation of large dynamic stall vortices of size comparable with the chord, associated massive alterations of the aerodynamic loads (e.g. augmented lift) and stall hysteresis effects, vortical structures in the present context are small-scaled, remain confined to the thin shear layer adjacent to the body surface and lift is affected insignificantly only.

There are many possible routes to turbulence. In the present context, one may roughly distinguish between two: In the first, ‘natural’, scenario, initially two-dimensional instabilities, the so-called Tollmien–Schlichting (TS) waves, develop in a fully attached laminar boundary layer due to the presence of external disturbances (e.g. sound waves, surface roughnesses, structure vibrations, freestream turbulence, etc.). The downstream travelling TS waves gradually grow in amplitude until they break down to form three-dimensional vortical structures and, eventually, turbulent spots appear. Since this transition mechanism typically is triggered by small perturbations, at least the initial stage can be described by linear stability theory based on the parallel flow assumption, which is well established in the literature, see e.g. [3]. A rigorous high asymptotic analysis of the TS wave generation process, the so-called receptivity problem, leads to the conventional triple-deck scalings originally discovered by Stewartson [4], Stewartson and Williams [5], Neiland [6] and Messiter [7] in completely different contexts, see Goldstein [8, 9] and Ruban [10]. However, a thorough theoretical understanding of the evolution towards full transition via this route is still lacking.

In contrast to this seemingly ‘slowly’ evolving scenario, the second route, i.e. the one we focus on, starts from a laminar boundary layer, but on the verge of separation. Whereas separation and turbulence are two different things, it turns out that separation very greatly enhances the propensity towards the formation of turbulence. Due to the presence of an internal inflection point in the velocity profile of the separated shear layer, the flow is highly unstable from the very beginning and inherently requires a nonlinear description to fully resolve the underlying physical mechanisms. Above all, this pathway to turbulence bypasses the development of TS waves (therefore may be called ‘bypass transition’) and runs through a cascade of characteristic stages with distinct length and time scales which become progressively shorter [11].

The increased use of high-tenacity materials like carbon-fibre-reinforced composites allows the realization of relatively thin airfoils with the objective to minimize viscous drag. Their lift coefficient is proportional to the angle of attack for . However, they are prone to the occurrence of a laminar separation bubble (LSB) at their leading-edge suction side when exceeds a certain value which still is small. Since LSBs may trigger rapid transition under certain conditions, their development crucially impacts the drag because the transition location then ‘jumps’ from its usual, ‘natural’ position at the mid to rear section to almost the leading edge of the airfoil. A further increase of eventually leads to the entire disintegration of the LSB with the consequence of abrupt and complete loss of lift (stall) and a substantial gain of drag. These circumstances explain the undiminished research interest in this phenomenon over more than eighty-five years, cf. the paper of Jones [12], who was the first to discover the phenomenon of shallow separation bubbles on airfoils experimentally in a wind tunnel facility.

As is well known, in practice high flows are turbulent, and hence one can expect that there will be a natural tendency for structures with small length and short time scales to develop. Scientific observations indicate the formation of well-defined, coherent vortical structures, the so-called -vortices or ‘rollers’, in the rear region of the bubble that are convected downstream, see e.g. [13-15], and Fig. 1. Therefore, it is obvious that there is an inherent unsteady behaviour related to LSBs.

Fig. 1

Colour visualization of a transitional separation bubble at the bottom wall of a laminar water tunnel: successive snapshots during a full -vortex generation/disintegration cycle (side view, flow from left to right, characteristic Reynolds number ). Source U. Rist & M. Lang, Institute of Aerodynamics and Gas Dynamics, University of Stuttgart

The mentioned structures, of size comparable with the laminar boundary layer thickness, convey the early stage of the turbulent part of the flow and persist a certain distance downstream of the mean reattachment point until they disintegrate into small-scale vortices (the ‘cloudy’ pattern at the top right of Fig. 1) characteristic of the developing turbulence. Again, this laminar–turbulent transition scenario generated by a locally detached shear layer is not yet fully understood, but it is commonly accepted that the repetitive bursting of a LSB, i.e. the (largely) periodic vortex shedding at its rear has to play a key role insofar as it precedes and initiates rather than follows full transition, see e.g. [16].

The genesis of a rigorous high asymptotic framework for the description of (steady, planar) laminar separation bubbles dates back to the early 1980s in the course of the vibrant development of triple-deck theory and its application to fundamental flow problems in those days [17-19]. Considerable efforts were undertaken to extend the promising approaches with regard to time dependence [20, 21], and the third spatial dimension (spanwise direction), e.g. [22-25]. One of the great challenges in dealing with high- flows is their susceptibility not only to separation but also to instability. Specifically, the nonlinear evolution equations which arise from the asymptotic analysis typically exhibit stiff behaviour and, even when starting from perfectly smooth initial data, their solutions may eventually develop singularities [26]. In the present context, the progress was hampered by the occurrence of difficulties in the sense that the simplified mathematical models resulting from the perturbation approach turned out to be ‘ill-posed’ [27]. Regularization methods in this regard were not yet known and moreover the computing power available at that time was still insufficient to solve the occasionally delicate nonlinear evolution equations numerically. Of course, high-end computing power is available now, and as far as the ill-posedness is concerned, we believe these difficulties can be overcome as well [28].

From the current asymptotic viewpoint, the bursting process of separation bubbles evolves through a sequence of stages, Fig. 2. Each of them can be posed as a singular perturbation problem to be solved by application of matched asymptotic expansions. The involved spatial and temporal scales relevant in each stage and, consequently, the dominating physical mechanisms appear in a natural way. Most important, the perturbation approach leads to universally valid results in the form of similarity laws independent of the specific flow under consideration, which, nevertheless, require numerical treatment. As shown below in more detail, the successive passage of these distinct stages clearly exhibits features of a developing turbulent boundary layer: the shift of viscous effects to a wall layer of decreasing thickness and the release of vortices into a detached, predominantly inviscid, rotational region of increasing thickness. During these early stages of the transition process the overall flow structure, comprising an outer inviscid flow field and a viscous boundary layer adjacent to the solid wall, is retained and not destroyed by bubble bursting (as for e.g. in the case of massive, break-away separation, not discussed here).

Fig. 2

High Reynolds number asymptotic (layer) structure of incipient laminar–turbulent transition in laminar separation bubbles (schematic, cf. Fig. 1). Classical laminar boundary layer at the verge of separation (black): inviscid outer flow , viscous region . Consecutive interactive stages of separation (red) and spike formation (blue), followed by self-induced vortex wind-up (magenta, shifted to the right for better illustration) accompanied by finite-time blow-up events, here modelled up to the Euler–Prandtl stage. Long-term objective: linkage to the time-mean two-tiered turbulent boundary layer flow description [30, 31], (purple)

Section 2 serves to present essential preliminary work regarding the high asymptotic description of LSBs in the perspective of the present approach and to motivate the starting point of our study. Readers who are familiar with marginal separation theory can skip this chapter and start with page 8 or 9. Section 3 examines and extends the pivotal work of Elliott and Smith [29] as related to a systematic numerical investigation of the spike formation stage and its finite-time breakdown. Moreover, a composite asymptotic model is presented in Sect. 4 for the purpose of resolving the ill-posed Cauchy problem associated with the current triple-deck formulation. The noninteractive Euler–Prandtl stage which is initiated by the finite-time blow-up of the preceding triple-deck stage is outlined in Sect. 5. The main findings of the present work are briefly summarized and discussed in Sect. 6 and more detailed descriptions of the asymptotic layer structure of the triple-deck stage and the applied numerical technique are provided in Appendices A and B, respectively.

From classical to interaction boundary layer theory

Here we briefly outline the current status of how the onset of laminar–turbulent transition associated with localized boundary layer separation can be described by an asymptotic analysis. To this end the most important scalings, similarity laws (respective leading-order evolution equations), and main characteristics of each subsequent stage are summarized and discussed. The fundamental perturbation parameter of the present study is the Reynolds numberwhich is assumed to be large. More precisely, we consider the singular limit problem of the Navier–Stokes equations applied to incipient separation and, for the sake of clarity and simplicity, restrict the treatment to planar, incompressible flows. Extensions of the theory to incorporate (quasi-) three-dimensional flow cases may be found in [22, 25, 32, 33]. The inclusion of compressibility effects is treated in [19, 34, 35] for strictly sub- and/or supersonic flows and in [36] for transonic flow conditions. Here and below the superscript tilde indicates dimensional quantities and , and denote the freestream velocity, the kinematic viscosity of the fluid and a suitable reference length. The computational results are most conveniently expressed in terms of the elementary boundary layer characteristics displacement thickness and wall shear stress ,where , , , and denote the arc length measured from the leading edge, the usual boundary layer wall-normal coordinate, the density, the prescribed inviscid wall velocity and the boundary layer velocity profile in streamwise direction which satisfies the matching condition as .

Prandtl’s steady boundary layer at the verge of separation (stage I)

A comprehensive presentation of the subsequent stages of the early transition process necessarily includes the solution of Prandtl’s classical boundary layer equations as a starting point. To this end, the practically important example of the steady flow past the suction side of an airfoil at various small angles of attack is considered, Fig. 3. As can be seen from Fig. 3a, an increase of causes the wall shear distribution, indicated in blue, to form a local minimum in a region of adverse pressure gradient conditions (at about of the chord of the examined airfoil).

Fig. 3

Steady classical boundary layer flow at the leading edge LE (suction side) of a specific symmetric airfoil with relative thickness (airfoil data can be found in Supplementary online materials) in dependence of the angle of attack : a wall shear stress , b displacement thickness according to (2) downstream of the front stagnation point

The distinct value of the control parameter is characterized by vanishing wall shear stress at a single point , referred to as the limiting case of marginal separation (indicated in red). The corresponding kink in the displacement thickness distribution, Fig. 3b, already indicates a local breakdown of Prandtl’s hierarchical procedure: in this limit, the pressure correction due to the displacement effect results in the asymptotic expansion for the flow ceasing to be uniformly valid. For values of , the numerical computations lead to a severe breakdown at in the form of the Goldstein square-root singularity (magenta curves), where the wall-normal velocity component locally behaves like [37, 38]. As shown by Stewartson [39], viscous–inviscid interaction theory typically is not capable of removing (strong) Goldstein singularities (for a counterexample see e.g. [40]) and later Sychev [41] concluded that a proper flow description must avoid the appearance of these strong singularities in general. For instance, free streamline (Helmholtz–Kirchhoff) theory and the concept of self-induced separation is one possible way out to describe massive flow separation from smooth surfaces in the limit as , see e.g. [42].

Marginal separation (stage II)

Here, however, we are not concerned with the latter phenomenon, but rather focus on the weaker singularity in the limit , of incipient separation: if the local flow description enables the displacement and the thereby induced pressure to strongly interact on the same level of approximation, a self-consistent asymptotic triple-deck framework, commonly referred to as the theory of marginal separation, can be set up. It was originally introduced independently by Ruban [17, 18] and Stewartson, Smith and Kaups [19] and later generalized to include temporal dependence [20, 21], and also flow control devices. The pivotal outcome of this research is the fundamental equationfor the displacement function A(X, T), which is related to the boundary layer characteristics (2) through the local expansionsHere, is the local coordinate in streamwise direction, a measure for the deviation of the control parameter from its limiting value and is the time. The effect of external disturbances is incorporated by the quantities h(X, T) and , representing the shape of a surface mounted obstacle1 [44], and the velocity distribution in a suction slot [45], respectively. All quantities have been nondimensionalized and suitably scaled to represent O(1) quantities, see [43, 46] for details. Furthermore, and are positive constants. In (4), represents the value of the displacement thickness according to classical boundary layer theory in the limit at . The values of the constants , , as well as all other constants of proportionality entering the scalings depend on the specific flow problem under consideration. Representative candidates of flow problems whose theoretical description also can be condensed into a similarity law of the form (3) include channel flows with a suction slot or a displacement body at the upper wall [13, 15, 47], smooth backward facing step flows [3], retarded boundary layer flows of dense gases [48, 49], sub- or supersonic flows past flared cylinders [50] and viscous wall jets that are forced to deflect [51].

Comprehensive investigations of (3) for steady as well as unsteady, forced and unforced flows led to the following important conclusions, see Fig. 4: weak Goldstein singularities, represented by , are regularized by means of the interaction strategy and may result in localized reverse flow zones characterized by , whereas for no (real) solution can be obtained, reflecting Stewartson’s finding [39]. Consequently, there exists an upper bound for admissible steady solutions: . Furthermore, the (steady, unforced flow) solutions, most conveniently depicted in the ‘parameter space’ versus , Fig. 4b, bifurcate three times (at points/curves labelled c, e and g) and therefore are (multiple) nonunique in the range [19, 52]. A bifurcation analysis based on the limit process as leads to a reduction of (3) to a nonlinear differential equation of Fisher–KPP type whose specific characteristics indicate a significant change of the stability behaviour at the bifurcation points [32, 53, 54]: upper/lower branch solutions are stable/unstable in general and this allows for the distinction between short and long separation bubbles (grey toned areas in Fig. 4a represent the stable range), in qualitative agreement with experimental findings. Moreover, the critical value may effectively be altered by means of suitably chosen passive flow control devices h(X) and/or , [43, 53]. In case of time-harmonic forcing, the theory reveals what has been observed numerically as well as experimentally, cf. [55, 56]: the higher the forcing frequency of the control device, the higher the stabilizing effect on the separation bubble, see [57] for details.

Fig. 4

a Steady wall shear distributions according to (3) for unforced flow . Dashed lines: local solutions of classical boundary layer theory , cf. Fig. 3a. b Fundamental curve of marginal separation; dashed line: local solution of classical boundary layer theory (asymptote for ), dotted line: parabola approximation in the vicinity of the bifurcation point

The solution of the initial value problem for equation (3) induced by unsteady forcing and appropriate initial conditions shows the following characteristics: subcritical flows may remain bounded or blow up, i.e. develop singularities, within finite time (dependent on device location, forcing amplitude and frequency), whereas supercritical flows invariably lead to finite-time blow-up(s) [21, 57]. These observations, which essentially rest on the nonuniqueness of the steady two-dimensional base flow and its saddle-node bifurcation at , Fig. 4b, allow the formulation of a precise vortex shedding criterion: subcritical flows require a certain perturbation amplitude to trigger bubble bursting, i.e. the (repeated) release of a vortex, cf. Fig. 1. In sharp contrast, supercritical conditions are always accompanied by self-sustained vortex shedding, even in the absence of any perturbations.

As pointed out by Ryzhov and Smith [27], the Cauchy problem associated with (3) is ill-posed, i.e. the solution (scheme) is prone to short-scale instabilities. However, the study by Braun and Scheichl [28] shows that the application of a composite asymptotic model which includes higher order effects (such as e.g. the streamline curvature) successfully leads to a corresponding regularization.

A distinct feature of bubble bursting and likewise an important ingredient of the laminar–turbulent transition process is the tendency of the flow to develop progressively shorter scales in time and space. In the present context we consequently focus on solutions which lead to finite-time singularities since it is solely these which enable the scale transformations necessary for an asymptotic description. Interestingly, any of the above-mentioned blow-up solutions ultimately approaches a unique, self-similar structure, entirely independently of the choice of initial data, below or above critical flow conditions (), and, if present at all, the type of forcing [58]. The finite-difference scheme used to compute the spatio-temporal evolution, described in [58], turned out to suppress short-scale instabilities without any specific regularization measures. In the final phase of the calculation, the solution develops large amplitudes, which become more and more pronounced in a region of decreasingly spatial extent. Eventually, the numerical results suggest the occurrence of a singularity at a single point , at the blow-up time in accordance with the predicted scalings of the terminal structureas , [21]. Here the blow-up profile is governed by the homogeneous ordinary nonlinear integro-differential equationIts unique nontrivial solution is depicted in Fig. 5, see also [58].

Fig. 5

Various quantities determining the initial condition (11) of the triple-deck stage and the asymptotic behaviour of A and , (34) in the limit as

Spike formation (triple-deck stage III)

As shown by Smith [21], and Elliott and Smith [29], the subsequent stage following the finite-time breakdown of marginal separation characterized by (5) and (6) is a fully nonlinear, unsteady triple-deck interaction. Its essential equation, formulated in terms of the stream function , readssee Appendix A for details. The induced pressure and the displacement function are related to each other by the Hilbert transform, i.e. the familiar interaction law for incompressible flows, according toMoreover, the problem is subject to the no-slip boundary conditions at the solid wall , the matching conditions with the main part of the boundary layeras and , as . Here , and are the local coordinates in streamwise direction, wall-normal direction and the time. Again, all quantities have been suitably scaled to be of O(1), see Appendix A. In the triple-deck stage of the bursting event, the length and time scales are already rather short, and—in contrast to the conditions before—the induced pressure gradient reaches the same order of magnitude as the imposed (adverse) pressure gradient. Moreover, the local displacement effect and the action of the wall shear stress (2) become intensified,cf. (4). Here, and again denote problem specific constants.

Most important, the connection of (7) to the terminal structure of (3) is ensured by specification of the initial (or equivalently matching) condition with the Hilbert pairsand as . Here we have used the scalings , , and . For details, consult the investigation of Braun and Scheichl [28], but please note that in (11) a slightly different formulation of the stream function is used for reasons of the numerical implementation. Substitution of (11) into (7) recovers the blow-up structure of the preceding stage in form of the staggered set of boundary layer equations linearized about the sublayer separation profile ,subject to the adapted no-slip boundary conditions and the matching condition determinable from (9). Here the leading order term represents the nontrivial eigensolution of the singular operator equation which meets all the boundary conditions. Application of the adjoint operator approach and Fredholm’s alternative as described in [28] to the second-order problem subject to the no-slip conditions , at and the interaction law eventually leads to the solvability condition (6) which determines the unique blow-up profile , Figs. 5 and 8a. The presence of the eigenfunctions , with the arbitrarily choosable amplitudes (reflecting the ‘history’ of the flow, or, in other words, which allows the embedding of the local flow description into a global one) shows that the spatio-temporal evolution of is not uniquely determined by the local analysis of the flow but, nevertheless, predictable. A proper start of the triple-deck computations requires the setup of an accurate numerical approximation of the time derivative term in (7) in the limit as . To this end, the particular solution to (12), chosen such that it remains bounded as , and the associated homogeneous part is also determined. The latter of course results from the solvability condition for the third-order problem with , at . Within the scope of the present study, also the higher order terms and were determined (similarly to ), see Figs. 5 and 6a, b.

Fig. 8

a Evolution of the displacement function (blue) and the scaled wall shear stress (red) according to (7)–(31) for ; each minimum location of A is denoted by a full circle, its asymptotic prediction (35) by a dotted line. b Continuation of (a) for . Grey lines depict the results for where the used resolution fails to meet the imposed accuracy requirements, cf. Figs. 10 and 11. c Corresponding interaction pressure (magenta) for of (a) and (b)

Fig. 6

Contour plots of the perturbation stream functions (streamlines) according to (12): a , in increments of (black lines); —minimum: . Grey lines and full circles indicate the computational grid with resolution and according to (16), (22); only every seventh grid line is plotted in wall-normal direction. b , in increments of and (black lines); —minimum: , —maximum:

With these prerequisites at hand we turn to the main objective of the present work, namely the presentation of our numerical approach for solving the initial value problem (7)–(12).

Numerical treatment of the triple-deck stage

The numerical method for solving the problem (7)–(12) comprises several components and is structured as follows.

Spatial discretization and differentiation. Initial attempts to approximate the spatial derivatives with a finite-difference scheme proved to be totally unreliable in terms of accuracy, especially regarding the derivatives of second and third order in the wall-normal direction. A significant increase in accuracy was achieved by the application of a spectral collocation method based on Chebyshev polynomials both in wall-normal and streamwise direction. The use of (pseudo-)spectral methods for the solution of triple-deck problems, especially in the (nonlinear) flow separation regime, originally goes back to the work of Burggraf and Duck [59], see also [33, 60]. We essentially follow the comfortable differentiation matrices approach introduced by Gajjar [61] for the treatment of viscous–inviscid interaction problems. To this end, the (semi-)infinite physical domains are mapped to the Chebyshev domain , see [62, 63]. Specifically, any (unknown) field variable f is interpolated by a polynomial p(s) in Lagrange form, where the basis polynomials are represented by the barycentric formula, see the very instructive article by Berrut and Trefethen [64] for details:Here , and denote the intrinsic unknowns of the problem, the unequally spaced Gauss–Lobatto grid points (nodes) for the interpolation process to be well conditioned and the weights with indices , respectively. Furthermore, we make use of the usual Kronecker delta when and 0 otherwise.

Spatial differentiation of order m then is easily performed via differentiation matrices whose entries can be computed analytically [64], e.g. for the differentiation order To minimize the effects of roundoff errors, we follow [65] and compute the higher order differentiation matrix entries iteratively viaIn the context of numerical schemes involving Chebyshev polynomials the interested reader is referred to a variety of useful references in the open-source software project Chebfun2 [66-68].

The domain mappings are accomplished by means of the nonlinear transformationswhere and B, C denote appropriate origin shifts and positive stretching parameters, respectively; see Fig. 7 and Table 1 for the typical numerical parameter settings. As a consequence, the elements of the differentiation matrices in the physical space, e.g. for the wall-normal direction, are then given by application of the chain rule (Faà di Bruno’s formula)Here the subscripts on s denote derivatives with respect to y, e.g.Spatial integration. In a similar fashion, a matrix–vector notation is also used for integration, albeit the determination of the integration weights (matrix entries) requires numerical evaluation in general. Here, the aim is to apply the above introduced polynomial representation (13) of a field quantity to the Weyl fractional integral operators appearing e.g. in the equation of the blow-up profile (6). To this end, integration by parts is performed in order to ensure integrability, i.e.Here and sufficient decay as is assumed. Using (13) and the convention of vanishing values of f at the boundaries, , we finally obtainwith according to (16). By analogy,For the Hilbert transform one deduceswith the symmetry property . The determination of the matrix entries and is performed in a preprocessing step by means of numerical integration (quadrature) routines of the NAG C Library (Mark 26) suitable for highly oscillatory integrands (d01skc) and which allow for (integrable) singularities at the user-specified break-points (d01slc).

Fig. 7

Computational domain, Gauss–Lobatto grid and allocation of the equations of the triple-deck stage

Table 1

Standard specification of the numerical parameters (unless otherwise stated)

m	n≈m/2	om	B	x^∗	C	y^∗	D	E_1	E_2
100…250	50…125	3n/4	0.65	0.9	0.7	0.	1.	1.	1.

Singular parts, similarity properties and temporal discretization. In order to capture the singular behaviour of as , (9), and its similarity properties as and , we writewhere , , and solve for the (bounded) intrinsic unknowns , and in a fully implicit manner. The thereby used scaling functionshave the properties as , when and as , when . Moreover, the time derivative is replaced by a backward finite-difference approximation of second- or third-order accuracy such that adaptive time stepping on a mapped domain is possible [58],Again, here denotes a stretching parameter, typically of O(1), table 1. However, for the asymptotic regime of large negative times finite differencing, since derived from truncated Taylor series expansions, is inappropriate and hence the time derivative for the first couple of time steps is computed analytically based on the ‘matching-in-time’ condition (11). Using the relations , , etc., cf. (11), one obtains e.g. for as . Here , and the function values of , , etc. at the temporally fixed grid points (cf. Fig. 7) are determined by means of Chebyshev interpolation in each time step. A similar relation is applied for the time derivative of , see [69] for details.

Algebraic decay. The simple application of the above-introduced differentiation and integration matrices in some way obfuscates the properties of the underlying interpolation polynomials. If the Chebyshev polynomials of the first kind and degree are combined with the transformation (16) for the streamwise coordinate, then expanding for large |x| yieldsi.e. algebraic decay behaviour with negative integer powers of x. The unknown field quantities displacement function A, pressure P and perturbation stream function g indeed decay algebraically to zero (or a bounded value) in the far field, but not necessarily with a negative integer power of the respective coordinate. The quantity concerned, say f(x), therefore is written as a product [70],where is the unknown part to be interpolated and a weighting function to be specified such that the decay of as is sufficiently strong with . In practice, however, the decay of f may not be known (exactly) or may change with time. In this case, the specific far-field properties of the weighting function are adjusted by numerical experiments such that satisfactory results are obtained.

Boundary conditions and implementation of equations. After application of (21), the modified no-slip boundary conditions at the solid wall readFrom the far-field condition (9) it follows that as , whereandOn writing and using the mapping (16) for the wall-normal coordinate , one obtains the exact relationfor the interaction pressure P. Here we want to point out that for the first derivative the differentiation matrix defined by (14) on the original Chebyshev interval has to be used. For more details see Appendix B, in which a test example (Blasius boundary layer flow) serves to demonstrate the spectral method on an unbounded domain. Moreover, in (30) the interaction law (8) is used to eliminate the unknown variable P. Since A and g decay to zero as , and to account for the elliptic character of the Hilbert transform (8), we setThe assignment of the individual equations to the grid points , , of the computational domain for the determination of the unknowns displacement function and stream function at each time step can be seen in Fig. 7. The use of has proven advantageous in practice. Obviously, the equation of motion (7)—in its via (21) transformed version—cannot be evaluated at the upper edge , ; instead, the far-field condition (30) is prescribed there. In order to ensure the fulfilment of the boundary condition (27), the implementation of (7) is omitted at a more or less arbitrarily selected vertical position . The entire nonlinear algebraic system of equations is solved in a fully implicit manner with the NAG3 c05qbc routine in each time step. On the computer platform used (Intel® Xeon® CPU E5-2680 v3 @ , up to 16 threads), the computing time for a single time step with a grid resolution of is initially (i.e. for large negative times) about , towards the end of the computation (before the suspected singularity occurs and the scheme no longer converges at all) up to about with a maximal memory usage of . The full computation (typically time steps) for the mentioned resolution takes approximately 1400 hours CPU time (effectively hours on our server).

Reliability and consistency checks of the numerical solution. The following quantitative criteria are used to evaluate and monitor the quality of the solutions found in each time step. An obvious test is to determine the maximum absolute values of the residuals (res) of the (transformed and) discretized versions of equation (7) omitted at the vertical position and of (28) at , , respectively:Typical temporal evolutions of and for different spatial resolutions are depicted in Figs. 10 and 16 . Sufficiently small values of and , say less than , are required for an acceptable solution, but the decisive numerical solution acceptance criterion is derived from the assessment of the decay behaviour of the amplitude spectrum of relevant quantities at the high wave number regime. The amplitude spectrum —i.e. the magnitudes of the Chebyshev expansion coefficients—of the concerning approximated function f is given bywhere for and else, see e.g. [71]. An alternative method to efficiently compute the coefficients is to use the discrete cosine transform (DCT). The selected spatial resolution of the numerical method (maximum polynomial degrees m, n) is considered satisfactory if the amplitudes of the quantities sought, e.g. , , etc., decrease sufficiently with increasing ‘wave numbers’ k; the corresponding quantity is then referred to as ‘fully resolved’. Typical amplitude spectra are displayed in Figs. 11 and 17 : a resolution that at first is sufficient to compute the spatio-temporal solution behaviour in the initial phase ultimately always proves to be insufficient when a blow-up event is emerging.

Fig. 10

Error indicators (red) and (blue) according to (32) for various spatial resolutions : —, —, —. Full symbols indicate the upper limit of reliable solutions. Adaptive time step distribution (green). By means of (36) estimated blow-up point: , (black dashed line)

Fig. 16

Effect of regularization, (61), on the error indicators (red) and (blue) according to (32) for a spatial resolution of : —, —. Full symbols indicate the upper limit of reliable solutions; predicted blow-up time

Fig. 11

Modulus of the amplitude spectrum (33) of and for various spatial resolutions and specific time steps , . For better readability, the lower wavenumber range has been left out for the higher resolutions

Fig. 17

Effect of regularization: modulus of the amplitude spectra (33) of and for fixed grid resolution at various time steps. (coloured), ; (grey), ,

Numerical results

A representative numerical solution of (7)–(31) based on the proven parameter settings listed in Table 1 is depicted in Fig. 8: from (4) it is evident that the blow-up profile is the starting form as for both the displacement function A and the (scaled) wall shear stress , cf. (10),As time passes, these quantities deviate from each other, grow in amplitude and finally run into yet another finite-time blow-up scenario. Please notice that according to (10) the correction to the displacement thickness is proportional to . The thus changed presentation of the results in Fig. 8b renders their interpretation considerably easier. In particular, the close resemblance between the temporal development of the displacement profile with the spike formation in the separating streamline visible in the third picture from above in Fig. 1 becomes apparent. Figure 9 shows the corresponding development of instantaneous contours of the stream function for the terminal time steps.

Fig. 9

Instantaneous contour lines of the stream function according to (7) at , , cf. Fig. 8b, c (see the Supplementary material for a video file). Spatial resolution ; for , the solution no longer meets the accuracy requirements. Isolines: . —location of minimum of at :

It is important to ask whether every prescribed initial condition (11), which is essentially determined by specifying the amplitude values and of the eigenfunctions and , respectively, must ultimately lead to a blow-up event. To put it the other way round: is it possible to choose and such that the evolution of , and governed by (7) and (8) converges to the trivial (steady) solution , without the occurrence of any finite-time singularity? In order to gain some insight into this matter, we examine the development of the minimum of the displacement function and its location over time. We start with large negative values of time when (11) is valid. Setting in (11) we findas . Here the specific numerical values of the involved quantities are , , , , , and , cf. Figs. 5 and 8a. The inspection of (35) shows that all expansion coefficients for up to are negative regardless of the sign of . It, therefore, seems impossible that even large negative values of at can change anything in the pronounced focussing tendency of the displacement function which is dominated by the leading order terms.

Eventually, one can ask whether the inclusion of properly chosen flow control measures (e.g. surface mounted obstacles of height of and/or suction slots with wall-normal velocity component of ) can avoid the otherwise (most likely) inevitable finite-time breakdown of the triple-deck formulation (7). This task has not yet been investigated.

Finite-time blow-up

The currently manageable maximum resolution in space is a limiting factor. It causes the computation to terminate prior to the entry into the asymptotic blow-up regime, described below, due to the occurrence of spurious oscillations (‘wiggles’). They can definitely be ruled out as being a result of short-scale instabilities, their appearance has to be attributed to the limited spatial resolution. The finer the grid is chosen, the later the numerical solution starts to oscillate, in case of short-scale instabilities, the reverse would happen. In particular, inspection of the modulus of the amplitude spectrum , Fig. 11, shows that the individual three different spatial resolutions yield reliable solutions up to . However, up to time it is only possible to find an acceptable solution with the finest grid . To accurately predict the solution behaviour for even the finest resolution is no longer sufficient.

Nevertheless, the present results heavily underpin the similarity structure put forward in [29],as the blow-up point , is approached. The interaction relation between and reads as in (8) and the blow-up profiles and are governed bysubject to as , as . Interestingly, this terminal structure is inviscid in nature and therefore a (passive) viscous sublayer is required such that the no-slip condition at the solid wall can be fulfilled [29]. Introducing the sublayer scalingsone obtains the modified boundary layer equationsubject to at , as [70]. Here the slip velocity and the pressure are imposed and related to each other by Bernoulli’s equation .

A detailed numerical analysis of (37) and (39) based on a method similar to that described above showed that the solution to (37) is not unique, but depends on two parameters, say, and [70]. For each set of the two parameters for which a solution to (37) could be found, the discretized version of (39) proved to be solvable as well, see Figs. 12 and 13 for a typical example. Without doubt, the solution depicted in Fig. 12 resembles features already seen in the final development of the corresponding quantities in Fig. 8b and therefore justifies the scalings (36). From a physical point of view, the terminal form (36)–(39) of the triple-deck stage indicates the emergence of a predominantly inviscid region where the crucial dynamics of the bursting process take place. Viscous effects are then confined to a narrow region adjacent to the solid wall. Nonetheless, a proper overall flow description requires a compatible viscous sublayer solution to exist.

Fig. 12

Triple-deck blow-up profiles, (37), (39): displacement function , pressure , slip velocity and wall shear stress . Two-parameter solution with prescribed values and [70]

Fig. 13

Triple-deck blow-up profiles: contour lines of a according to (37); resolution , isolines , minimum and b according to (39); resolution , isolines . Two-parameter solutions with prescribed values and [70]

Treatment of the ill-posed Cauchy problem (stage III)

Worth mentioning is the fact that some inconsistencies associated with the proposed Cauchy problem (7), (11) were discovered in [29], and their repercussions on the solvability in general necessitate a comprehensive discussion. In that paper, the occurrence of instabilities was put in the context of an ill-posed Cauchy problem. More precisely, it was shown that the instability of the velocity field against short-scale disturbances and, entailed by that, the incorrectness of the problem are a direct consequence of the abnormal dispersion relation governing the linearized problem in the limit of very high wave numbers k.

To reconsider this issue, we start from the assumption that the solution to (7) near a position and for a short time interval might be expanded aswhere , and denote the concerning function values at the sampling point . The additional terms , and supplementing the Taylor series expansions of the background flow account for the emergence of possibly unstable waves with amplitudes , and , respectively, the sizes of which are defined below. For the description of these short-scale instabilities the distinguished limit is chosen in which a balance between pressure gradient, local time derivative and convective terms exists. Therefore the necessary scalingscan be readily obtained from (36). Here, however, the disturbances are assumed to occur prior to the blow-up structure and to be smaller in magnitude, implying the order-of-magnitude relationshipRecalling that as , one then obtains from (7) together with (40) and (41) the linear leading-order equationsubject to the conditions as , as and . The viscous sublayer necessary to bring the slip velocity satisfying at rest is of thickness . However, the latter problem will not be investigated in the following.

Using a Fourier approach in the sense of a stability analysis such that each quantity f is expressed aswhere is real, and solving the ordinary differential equation resulting from (43) then yieldswith understood as the principal value. In the limit as , behaves asApplication of the constraint , i.e. and the far-field condition stated above gives , which, when and are replaced by their unscaled counterparts and , readsThis local dispersion relation for the high-wave-number limit has two complex roots, one is associated with stable and the other with unstable modes. However, as also pointed out in [29], the growth rates of the thus occurring short-scale instabilities obey and are therefore not bounded from above. Consequently, the given Cauchy problem (7), (11) turns out to be ill-posed and has to be regularized. It should be noted that the above analysis is similar to the triple-deck study of the Blasius boundary layer with respect to the so-called upper-branch stability, see [72], insofar as, to leading order, viscosity does not enter the dispersion relation. However, for the case of an attached boundary layer, the leading-order dispersion relation for high wave numbers leads to neutral modes only, whereas here the situation is obviously completely different.

Composite asymptotic model

The essence of the above result is that higher order terms in not included so far must play a distinctive role in the evolution of the separation process, given the bounded spectrum the unsteady Navier–Stokes equations are assumed to have, but which the leading order problem fails to deliver. However, to the authors’ knowledge, there does not exist a rigorous asymptotic theory for the regularization of otherwise ill-posed Cauchy problems. In other words, the filtering process entailed by the temporal and spatial scales the asymptotic analysis is based on may lead to the undesired result that the terms needed for regularizing the Cauchy problem at a certain level of approximation cannot be incorporated into it. They will thus enter the equations of higher order and, by that, form the higher order forcing terms. Consequently, the only way to regularize the original problem is to set up a composite asymptotic model, where these corrections are taken as part of the leading order problem such that they can effectively contribute to stabilizing its homogeneous solution against short-scale disturbances.

This method has been proven very successful in the past, see e.g. [73] and [74] as well as [75] and the references therein, and, in particular, for studying the unsteady marginal separation stage preceding the triple-deck stage discussed here, see [28]. As revealed in these studies, for higher order asymptotic terms to play the role of regularization terms, they must at least generate derivatives with respect to t and x of orders that are higher than those already present in the system. Furthermore, in [28] it was shown that the induced streamwise pressure gradient in the lower deck which replaces the original one in order to form a composite model equation has to include the effects of the streamline curvature in the main part of the boundary layer. Surprisingly, a second correction term, which is due to the unsteady nature of the main deck, was found. Whereas both appear naturally in many triple-deck problems as higher order contributions to the leading-order pressure gradient, see e.g. [72, 76, 77], only the one stemming from the streamline curvature is usually used as regularization term, [73–75, 78]. The following analysis of the triple-deck problem considered here will show that, in principle, the regularization presented in [28] can be taken over almost unchanged.

Inspection of the far-field condition of the lower-deck solution as , see (9), shows that in the main-deck region MD, where , the expansion of the stream function assumes the formPlease notice that in the following, we will use the unscaled versions of the field quantities and length scales, i.e. those before the affine transformations (72) were applied. Substitution of (48) into the streamwise momentum equation yields to leading orderand additionally,The last equation in (50) is interesting insofar in that here, the induced pressure component of leading order enters. Furthermore, the pressure generated by in the upper deck and acting in the lower deck will form the component of in the pressure expansion. This fact suggests a pressure variation in transverse direction to come into play, since the normal gradient of precisely this term of is comparable in size with the leading-order convective term in the transverse momentum equation. Consequently, the extended version of the streamwise pressure gradient in the main deck is given bywhere all terms, except for , are identical to the corresponding components in the lower deck. As already indicated by the above order-of-magnitude argument, the asymptotic expansion of the transverse momentum equation then leads to the relationshiprevealing the intrusion of a normal pressure gradient in the MD. As will be shown in the following, the last equation in (50) together with (52) will induce a pressure response in the lower deck that is proportional to the curvature of the streamlines, i.e. the second derivative of the displacement function with respect to the streamwise coordinate.

After applying the rules for matching with the lower deck and extracting the singular parts of the resulting integral, one can write the solutions to (49) and (50) asPlease note that for the derivation of these expressions, the expansions and have been used (see e.g. [42]). In addition, the equation for the normal pressure gradient (52) results in the expressionHere, denotes the induced pressure of order acting in the lower deck, which remains to be determined later.

As mentioned above, the purpose of the analysis presented here is the identification of those pressure terms that are promising candidates for regularizing the Cauchy problem stated in (7) and (11). Hence, they have to be related directly to the leading order effects represented by . For simplicity, therefore, lower-deck displacement functions of higher order, such as , and , have been omitted in the main-deck solutions.

The upper deck (UD) expansions, with the appropriate coordinate normal to the wall now being , areandHere, the quantities , represent the induced pressure fields, which decay for large distances from the origin. Through substitution into the momentum equations, it is easy to verify that the terms in the stream function and the pressure expansions satisfy the linearized Euler equationswhere is the Laplacian .

In general, the solutions for the pressure fields are sought in the semi-infinite half-plane with the the relations for the wall-normal pressure gradients taken as Neumann boundary conditions at . For the calculation of the streamwise pressure gradient at , a simplified approach is possible by exploiting the fact that the derivatives and , satisfy the Cauchy–Riemann equations such that the complex functions are analytic. As a consequence, the imaginary part of each function at is given by the Hilbert transform (8) of the real part. Relying on the matching principles and taking into account that apart from the MD pressure terms do not depend on , one then obtains from (57) and (53)As expected, the well-known leading order result for is recovered. Moreover, the last equations in (57) and (53), respectively, together with (54) lead toBy invoking the relation , one can infer from (58) and (59) that the fourth component of the lower-deck pressure gradient is given byFurthermore, it should be noted here that as a matter of course, the singular parts of the solutions (53) and (54) as and match seamlessly with the corresponding counterparts in the lower and upper deck, respectively.

As a consequence, when the terms from (58) and (60) are combined, the affine transformations (72) are reapplied and is redefined, the induced pressure gradient replacing in the lower-deck triple-deck problem (7) in order to form a composite model equation can be written asThe constant C has to be calculated from the unperturbed boundary layer separation profile with and the prescribed pressure gradient at . It turns out to be positive in all cases investigated, e.g. for the Falkner–Skan similarity separation profile characterized by the critical Hartree parameter one obtains .

On substituting into (45) and rescaling k and the angular frequency in a very similar way as in (41) but now with respect to such that and , we obtain the modified and perturbation parameter-free dispersion relation for the compound modelwith the last two terms on the left-hand side resulting from the higher order terms in (61). Inspection shows that for , the modes associated with both roots become neutral, see Fig. 14, and thus the growth rates are bounded from above. This, in turn, gives good reason to assume that the nonlinear Cauchy problem based on this corrected model is well posed. Given the fact that all terms in (62) are of the same order of magnitude, one can furthermore conclude that as soon as the temporal and spatial scales become so small that and the problem will be well posed by itself, i.e. without the necessity to include additional regularization terms. This occurs when and , which are precisely the scalings of the Eulerian stage assumed to follow the present one, see Sect. 5.

Fig. 14

The two roots of the dispersion relations (47) and (62). Solid lines: (unstable) and (stable) for the composite model with ; if , origin: . Dashed lines: asymptotes corresponding to the original relation (47)

Numerical evidence of regularization

Detailed investigations of the numerical experiments that led to the results presented in Sect. 3.1 indicate that for the spatio-temporal ranges and resolutions used so far, the small but, nonetheless, always present terms resulting from the discretization based on Chebyshev polynomials are very effective in suppressing short-scale instabilities with unbounded growth rates. However, the computation of solutions that blow up in finite time is not feasible with standard methods, due to the required fine resolution of the domains in the neighbourhood of the singularity. In particular, approximations based on polynomial representations are known not to be well suited in this regard. Consequently, the apparently inherent numerical regularization cannot be taken for granted when the solutions approach their assumed finite-time breakdown and grid refinements are necessary. In order to investigate the predicted regularization effect of the modified interaction law (61) according to the composite model derived above, the numerical scheme already developed for the analysis of the two-dimensional triple-deck problem (7) was extended appropriately. Above all, care has to be taken that the damping of small-wavelength disturbances by incorporating higher order terms is not overdone as this could give rise to nonphysical delays in the development of the blow-up solution or nonphysical results at all. The crucial points to be considered in this respect are, first, to quantify the regularization mechanisms already contributed by the numerical scheme. Secondly, then the size of the higher order corrections as given in (61) have to be adjusted in an adaptive way in order for the problem to remain well posed in each time step.

Numerical experiments with and showed practically no differences to the case and therefore indicate an (unintentional) intrinsic regularization of the numerical scheme used. An additional filtering is observed when is raised above about 0.01...0.05 (depending on the resolution). A concrete sample calculation demonstrating the effect of regularization is depicted in Figs. 15–17. The results for the time steps shown clearly indicate the effective suppression of grid oscillations, i.e. the expected low-pass filter behaviour of the extended interaction law (61). Since the value selected here for the regularization parameter is relatively large, there is certainly cause for concern as to whether the solutions in question can still be classified as physically correct. For obvious reasons, should not be larger than absolutely necessary.

Fig. 15

Effect of regularization, (61), with on the evolution of and at . Grey lines display the corresponding results for at . Spatial resolution

Outlook: Vortex wind up (Euler–Prandtl stage IV)

The theory, outlined in Sects. 2 and 3, so far appears to agree with experimental observations regarding the early stages of the laminar–turbulent transition in separation bubbles (cf. Figs. 1 and 2 ), in the sense that it shows how an initially slender, nearly steady, reverse flow eddy of length scale much greater than the boundary layer thickness changes dynamically, through viscous–inviscid interaction mechanisms, into a thick blunt eddy spanning the whole boundary layer in the noninteracting, incompressible, fully nonlinear Euler–Prandtl stage formulated below.

Following the lines of reasoning put forward in [29], one can show that, in the triple-deck stage, the evolution of the flow towards the finite-time blow-up as is accompanied by a growth of the (originally viscous) sublayer at such a rate that , until its thickness is indistinguishable from the majority of the boundary layer, i.e. , cf. Fig. 2. This breakdown implies a new, predominantly inviscid, structure and occurs when the time interval reduces to . With new (local) scalings then given by , the relevant velocity components as well as the pressure become O(1) quantities and are subjected to a fully nonlinear evolution where . The governing equations of this new stage in terms of the stream function simplify to the Euler equationssubject to the boundary conditions at , as . Here, the classical boundary layer separation profile at the point of marginal separation acts as a rotational ‘convective background’ to the inviscid main flow. Again, a viscous sublayer is necessary to bring the slip at the bottom to rest: introducing the appropriate wall-normal coordinate results in the spatio-temporal evolution of the stream function and the velocities , being determined by Prandtl’s boundary layer equationsubject to at . Furthermore, as , where B denotes the slip velocity which is related to the imposed pressure via Bernoulli’s equation . The initial (matching) conditions of (63), (64) as are fixed byand the similarity variables , , and . Here, the leading order terms are determined by the triple-deck blow-up profiles , , and , which satisfy (37) and (39), respectively, see Figs. 12 and 13 . More precisely, the explicit quotation of in (63) leads to a slightly modified formulation of the blow-up profiles, but without loss of generality.

Current efforts aim, at first, to compute (at least) the next-order corrections in (65) in order for the associated initial value problems to be defined properly, and, secondly, to construct a suitable numerical scheme able to capture the vortex wind-up process governed by (63) and (64). As in the triple-deck stage, the infinite spatial and temporal domains and the initial similarity structure require special treatment, cf. (21). According to the local expansions of the boundary layer characteristics (2) in the Euler–Prandtl stagewe expect O(1) corrections to the original displacement thickness and a substantially enhanced action of the wall shear stress.

Although we started the investigation of transitional separation bubbles with the analysis of formally steady classical boundary layer flows on the verge of separation, rather weak unsteady forcing (blowing) in the first, large-scale, interaction stage of marginal separation eventually led us (via a second, shorter-scale triple-deck interaction) to a scenario that in literature is commonly referred to as unsteady separation, see the reviews [79, 80] and references therein. In sharp contrast to the usual setting underlying the model problems studied in the context of unsteady separation, the outer inviscid flow field is not given in advance in the current formulation (63), (64). The spatio-temporal evolution of the flow is entirely self-induced once kicked by the initial condition (65). Nevertheless, from solving the Cauchy problems (63)–(65) we expect the occurrence of the generic behaviour known from the investigation of unsteady boundary layer flows with prescribed pressure gradient: the abrupt eruption of near-wall fluid (vorticity) out of the viscous boundary layer in form of a finite-time blow-up. The existence of this singular behaviour was first reported by Blasius [81] in his study of an impulsively started circular cylinder. Another example was given by Walker [82] who studied the finite-time singularity in the boundary layer exposed to a rectilinear vortex. Then these results were re-examined numerically and analytically by various authors, including [83-89] and [90].

As pointed out in [79], however, a Rayleigh-type instability may arise in the Euler region well before a singularity develops in the noninteractive boundary layer evolution (64), see also [91]. From a computational point of view, a Lagrangian-type scheme might be advantageous compared to the more conventional Eulerian formulations [87, 92], although more recently the latter approach also proved suitable [93].

Concluding remarks

The present paper essentially is a revision of the fundamental work of Elliott and Smith [29]. A spectral collocation method based on Chebyshev polynomials has been developed for treating the second—shorter—interactive (triple deck) stage of marginally separated flows. Special emphasis was placed on solutions that blow up within finite time. The predicted blow-up structure could be confirmed with a high level of evidence. To this end, the high wave number behaviour of the amplitude spectrum of relevant quantities proved to be particularly suitable for assessing the reliability of the numerical solutions. The numerical scheme implemented turned out to regularize the proven ill-posedness of the underlying Cauchy problem, i.e. the tendency of solutions to unbounded growth of short-scale disturbances is sufficiently suppressed. Additionally, we have shown that a suitable modification of the usual interaction law by means of a so-called composite asymptotic model again contributes to a successful suppression of the predicted short-scale instability.

Although the asymptotic theory of marginal separation and its extension to a triple-deck stage is established for more than three decades, we believe it is a key approach able to unveil the secrets of the laminar–turbulent transition process in separation bubbles. Still, its potential is by far not exploited and important issues still remain open, which keeps this research area in a state of excitement. From our current perspective, these points comprise the following: (i) Treatment of ill-posed Cauchy problems associated with the asymptotic theory of separated flows, (ii) tracing the breakdown cascade further: vortex wind-up and shedding in the Euler–Prandtl stage, (iii) Providing tools for transition detection and control.

In order to keep track of the big picture and, at the same time, to shed light on the achievements of the asymptotic theory, we inspect the -dependence of typical quantities characterizing the status of the boundary layer with respect to laminar–turbulent transition. Such an approach may also serve as source of suggestions and reference points for laboratory experiments and computational fluid dynamics. For example, the theory of marginal separation (stage II) predicts a separation bubble length of , an initial (‘steady’) bubble height of and the slope of the (almost) steady separation streamline to be of , which are facts that still do not seem to be commonly known, cf. the comments in [94]. From the triple-deck model (stage III), a length of and a height of can be estimated for the emerging spike at the rear of the bubble, which is consequently still characterized as a slender phenomenon. In the Euler–Prandtl stage IV, both the length and height of the shedding vortex then reach , comparable to the thickness of the original boundary layer. Table 2 illustrates the successive order-of-magnitude reductions of the relevant time scales, the streamwise length scales, the viscous wall layer thickness and the enhancement of the wall friction velocity during incipient transition. Our admittedly ambitious long-term objective is to fully bridge the gap between the established concepts of both the laminar boundary layer and the time-averaged (fully developed, attached) two-tiered turbulent boundary layer, [30, 31] (and the references therein), via a correct asymptotic description of the whole transition process.

Table 2

Order of magnitude alterations of selected quantities during the early stages of bypass transition triggered by incipient separation in the limit as (, ), cf. Fig. 2

Stage	I	II	III	IV	Final^†
Time scale, t~u~_∞/L~	–	Re^+1/20	Re^-1/7	Re^-1/2	–
Length scale, x~/L~	Re^0	Re^-1/5	Re^-2/7	Re^-1/2	Re^0
Viscous layer thickness, y~/L~	Re^-1/2	Re^-11/20	Re^-4/7	Re^-3/4	Re^-1lnRe
Pressure, (p~-p~_∞)/(ρ~u~∞2)	Re^0	Re^-1/2	Re^-2/7	Re^0	Re^0
Wall friction velocity, u_τ	Re^-1/4	Re^-7/20	Re^-2/7	Re^-1/8	ln^-1Re

I—steady classical boundary layer, II—marginal separation, III—triple deck stage (spike formation), IV—Euler–Prandtl stage (vortex wind-up), Time-averaged, fully developed wall-bounded turbulence, see e.g. [30]

Table 3

Blasius flow characteristics displacement thickness and wall shear stress in dependence of the number of Chebyshev collocation points n with

n	10	20	40	80	Reference value [96]
β_1	1.2	1.2165	1.216776	1.216780622	1.216780621614862
f^″(0)	0.52	0.46957	0.4696001	0.4695999884	0.4695999883610133