The classical unsteady boundary layer: A numerical study

The transition process from a laminar to turbulent flow near a solid wall (e.g., at the suction side of a slender airfoil) is still not fully understood. We focus on the early stage of the flow transition based on boundary layer (BL) theory, that is, high Reynolds number asymptotic expansions, which reduce the Navier–Stokes equations into simplified forms. While steady solutions of the resulting classical BL equations and corresponding regularisations based on viscous–inviscid interaction theory are studied quite extensively, there is only few information about the fully time‐dependent scenario. Hence, we numerically study the classical unsteady BL and its breakdown structure, the Van Dommelen–Shen singularity, on the simplest example we can imagine, the incompressible planar fluid flow through a rectangular channel with suction. The depicted BL builds at the lower channel wall and follows initially a generic Blasius behaviour, which gets temporally modulated through suction by means of a slot located at the upper wall. As usual, the main part of the flow is inviscid and governed by Laplace's equation. Remarkably, the wall‐slip velocity, or alternatively the pressure gradient imposed on the BL, can be calculated in closed form and we have full control of its temporal dependency through the variation of the applied suction strength. This specific flow configuration allows us to split off the initial Blasius solution and promises high‐accuracy solutions of the remaining unsteady BL equation in the scaled stream function formulation. The Chebyshev collocation is used in the wall‐normal direction, where the infinite domain is mapped onto [ − 1, 1]. On the contrary, we discretise the mainstream direction with a finite difference scheme of second‐order accuracy, which proved to be more precise near the build‐up of singularities. In order to trigger a dynamic breakdown of the unsteady BL equations, we increase the suction rate continuously from zero to a level slightly above the critical value (obtained from the steady case). As a consequence, the system cannot relax to a steady solution and instead blows up in finite time at some mainstream position. We try to resolve this process with high resolution and the final goal is to verify the inner structure (the two‐dimensional blow‐up profile) of the Van Dommelen–Shen singularity and its role in flow separation and transition.


INTRODUCTION
We are interested in understanding the incipient process of laminar-turbulent transition as precisely as possible.For this sake, we study it in the framework of high Reynolds number () asymptotics.In contrast to full computations of the Navier-Stokes equations, which yield only the solution to a specific problem without a clear detailed interpretation, the asymptotic theory aims to give a more universally valid description of the phenomenon of flow transition.In the limit of  → ∞, the friction term in the Navier-Stokes equation becomes only relevant at certain flow areas; these are especially the areas near rigid surfaces, denoting a so-called singular perturbation problem.We imagine that the whole process of laminar-turbulent transition can be described through a cascade of such problems [1,2].The first stage in this cascade is described by the classical boundary layer (short BL) theory.Here, as introduced by Prandtl [3], the flow can be split into an inviscid main zone and a thin boundary layer near the wall, which interplay with each other only over different perturbation orders in a strictly hierarchical manner.Combined with an appropriate initial condition, planar BL flow is described in leading order by the BL equation for the stream function Ψ, where  is the mainstream coordinate and ȳ =  √  is the suitably scaled wall-normal coordinate.The appearing pressure gradient   gets imposed from the main zone's wall-slip velocity   (, ) via Bernoulli's law, −  =     +   . ( Under a large enough adverse pressure gradient (2), classical BL theory will actually break down, resulting typically in a blow-up of the wall-normal velocity  = −Ψ  .Usually, such a breakdown can be regularised via the introduction of new -dependent scales and the next singular perturbation problem.While follow-up stages have been extensively studied for a stationary starting point (i.e., the steady BL theory) [2,4], there are only few studies on the fully time-dependent scenario [5][6][7].
For this reason, we numerically study the classical unsteady BL equation ( 1) and its breakdown structure, the Van Dommelen-Shen singularity [6,8].We consider the simplest problem we can imagine, namely the incompressible planar flow through a rectangular channel with a suction slot, as emphasised in Section 2. We will see in Section 3 based on the first simulation results that the suction strength will serve as an ideal control parameter to locally influence the BL flow and to induce the breakdown of classical BL theory.The theoretical description of the breakdown structure as well as its comparison with our numerical results will be considered in Sections 4 and 5, respectively, followed by the conclusions in Section 6.

Channel flow
The flow scenario we are interested in, motivated by [9,10], is visualised in Figure 1.All quantities are appropriately non-dimensionalised using the reference velocity at the channel entrance  ∞ and the distance  to the suction slot.The channel has some height ℎ, and the relevant volume fluxes of the system are denoted by V.The suction slot of size  is positioned at the upper channel wall.Presuming large Reynolds numbers  ∶=  ∞ ∕ ≫ 1, the considered boundary layer (circled 2) builds at the lower channel wall and depends on the suction strength  0 ().
We want to mention that for one suction slot alone the incoming flow at the channel entrance  = 0 would depend on  0 .We avoid this through the introduction of a mirror suction of the same strength in front of the channel inlet (marked in grey).This gives us the perfect toy example to study the almost local influence of suction on the flow.There is only one control parameter, the relative suction rate (), defined as the ratio between the volume flux removed by the suction slot V and the constant volume flux at the channel inlet Vℎ , In the main zone (circled 1), the flow in leading order is inviscid and assumed irrotational, hence governed by Laplace's equation Δ = 0. We find an appropriate solution for the potential  with the method of distributed singularities and images, where we superpose the flow of artificial mirror suctions outside the domain such that  = 0 and  = ℎ represent straight and parallel streamlines.There are now two benefits.Firstly, since the Laplace operator Δ is inherently independent of time , we are able to freely modulate the time dependency of the flow through a temporal modulation of the suction rate ().And, secondly, the wall-slip velocity   =   (,  = 0, ) has a closed form, Hence, we have a pleasant control over the pressure gradient (2), imposed on the BL equation (1).

Numerical treatment
For the computation of the boundary layer equation ( 1), we bring the system to a better suited form, The new intrinsic unknown is the scaled stream function , which depends on a scaled wall-normal coordinate .There is some kind of freedom in the definition of the auxiliary function , we choose where the integral accounts for the growth of the boundary layer to some extent.In transformation (5), we also separated a similarity solution f(), which is known as the Blasius or flat plate solution [2], f already solves the BL equation in the limit of zero suction  = 0 (i.e., for   = 0,   = 1).
After the substitution into (1), we obtain the non-linear differential equation for , with the no-slip boundary, matching and starting conditions given by Still, for non-zero suction, our system is formulated in such a way that the Blasius solution holds near the channel entrance  → 0 + .In addition, at mainstream positions sufficiently downstream from the suction slot, the flow will only passively depend on the suction.Hence, we restrict our computational area in the mainstream direction to a finite domain  ∈ [  ,   ] and we use a simple finite differences scheme of second-order accuracy with  + 1 grid points.Whereas, in the wall-normal direction, we depict the whole semi-infinite domain  ∈ [0, ∞) through the mapping  =  tan(( + 1)∕4), with typically  = 3.The Chebyshev variable  ∈ [−1, 1] gets discretised by the Gauss-Lobatto points   = − cos(∕),  = 1, … ,  + 1, on which we interpolate Chebyshev polynomials [11].
For the calculation of ( 7) towards finite time blow-up, we will also use a mapping in the -direction just to increase the grid resolution locally near the estimated singularity position   (see Section 5).For the temporal discretisation, we use a finite difference scheme of second-order accuracy with an adaptive time stepping [12] in order to decrease the step size as we approach the singularity.Additionally, to reduce computation costs, we linearise the BL equation ( 7) with respect to the previous time step.Since the error order of the linearisation equates to the error order of the temporal scheme, we do not expect an overall change in accuracy.
It is noteworthy that, in order to use a linear solver, we have to bring the resulting system of equations to a matrix-vector form.Therefore, we lift the almost tri-diagonal finite difference matrix   and the dense Chebyshev differentiation matrix   into a combined vector space via the tensor product, for example, The successful computed solution will then be interpreted by means of the usual BL characteristics.A typical measure of the thickness of the BL is the displacement thickness Another characteristic function is the wall shear stress, with the familiar Blasius value f (0) ≈ 0.4696.Both  * and   depend on the scaled stream function .In (9), we also see that the stream function Ψ diverges linearly in the far field  → ∞ [13] due to the matching with the main zone.This linear growth is already captured by the Blasius solution, f( → ∞) ∼  −  1 with  1 ≈ 1.2168 [14] such that  is indeed regular.

TEMPORAL SETTLING
We show the first computation results and compare steady and unsteady solutions in a short overview.We begin with a stationary scenario, where the right-hand side of the underlying equation ( 7) is set to zero.We use a constant suction rate , (3), and we repeat our simulation several times with a successively higher value of , until we reach a critical value   , associated with the condition of so-called marginal separation [4,15].
In Figure 2, we visualise the suction-induced pressure gradient (2) and the resulting BL characteristics ( 9) and ( 10) for a value of  slightly smaller than   .We can clearly see that the displacement thickness  * and the wall shear stress   reach an extrema caused by suction, positioned slightly downstream of the suction slot.At the critical suction, the wall shear minimum reaches the zero value in the form of a kink, indicating the breakdown of the steady BL theory.For more details and regulations around   consider, for example, [2,4,15].
Things change if we consider the unsteady scenario.We introduce a time dependency in the system through a temporal modulation of the suction strength.We increase  slowly from zero (Blasius case) to some fixed value   .When   <   , then the system will relax ultimately, after a long time into the stationary solution.
But when   >   , then there exists no stationary solution the system can approach to, and we reach a finite-time blowup of the stream function .In contrast to the steady-state case, it is indeed possible to reach negative values of the wall shear stress [16], indicating flow reversal areas, which we discuss in more detail in the next sections.And our final goal is to examine this breakdown structure of the classical unsteady BL.
Naturally, there are many other interesting options for the time dependency of our control parameter  besides the monotonic ramp.For this sake, we want to mention that the evolution towards the blow-up is rather slow.For example, we can raise the suction rate  dynamically for a short period above the critical value and decrease it afterwards without obtaining a finite-time breakdown.10), ( 6) and ( 9).Respective values for the Blasius case are drawn in grey.

ASYMPTOTIC THEORY OF FINITE-TIME BREAKDOWN
There exists an asymptotic theory of the breakdown of (1), formulated as slipping upstream separation.In order to understand this process, before diving into the numerical results of Section 5, we shortly revisit the findings of [5,6,8].
Let us say the blow-up happens at some mainstream position   and time   , then the appropriate scales ( x, ỹ) ∼ ( 1), as we approach the blow-up point (  ,   ) are given by In accordance with the behaviour of the displacement thickness, the blow-up region flattens in the -direction but explodes in the ȳ-direction.The scales actually depend on further flow-specific constants besides   and   .A special role plays the small positive slipping velocity  while  0 and  are just scaling constants.The actual blow-up is located at an inner inviscid zone, which is basically defined by the elongations of the arising flow reversal area.Below and above this region in the wall-normal direction, there exist passive shear layers, which are just present for the matching with the no-slip at the solid wall as well as with the outer potential flow.
In correspondence with (11), the mainstream velocity  = Ψ ȳ in the blow-up zone expands as In leading order again the constant slipping velocity  appears, the next higher order function Ũ is associated with the blow-up profile and can be derived from the BL equation ( 1), yielding a balance between unsteady inertial terms independent of   .For the sake of brevity, an appropriate solution restricted to ỹ ∈ [0, 2 ỹ0 ( x)] has the following implicit form: where ỹ0 ( x) and Ũ0 ( x) are the location and value of the characteristic minima curve, defined through Ũ ỹ = 0 and plotted in Figure 5.The solution of ( 13) is symmetric around ỹ0 ( x) and related to the elliptic integral of the first kind,  (, ) [17].
Qualitatively, the blow-up zone ỹ ∈ [0, 2 ỹ0 ( x)] describes the flow reversal area.The maximum negative mainstream flow is located at ỹ0 ( x) and successively decreases as we approach the boundaries in the wall-normal direction.While the mainstream velocity (12) stays regular, its -derivative and the wall-normal velocity  diverge as  −   → 0 − , depending on Ũ = Ψ ỹ , F I G U R E 3 Adaptation of the time step Δ according to the most diverging quantity   .
F I G U R E 4 Spatial resolution in mainstream direction, adjusted to the singularity position   .
This builds the structure of the Van Dommelen-Shen singularity, [8].The wall-normal velocity  is the most diverging quantity, contributing to the BL equation.With respect to our scaled stream function formulation (5),  has the same temporal divergence properties near the blow-up as Ψ.

TRIGGERING FINITE-TIME BREAKDOWN
As briefly discussed in Section 3, we raise the suction strength to a fixed value   = 0.051 slightly greater than the critical value and observe afterwards how the scaled stream function  or rather its -derivative, evolves to infinity.
In the simulation, we proceed with a constant time step Δ until we reach a certain threshold.Thereafter, we reduce the time step according to the temporal evolution of the maximum of the wall-normal velocity  ∝ −  (see Figure 3).We achieve times near the blow-up of   −  ∼ (10 −2 ) and velocity magnitudes of   ∼  (10 4 ), making it plausible that we reached the asymptotic regime (14) to some extent.For this sake, we also needed naturally a high spatial resolution in the mainstream direction, as depicted in Figure 4. We use  = 10 000 finite difference points and a resolution near the supposed mainstream singularity position of nearly Δ ∼ (10 −7 ).In contrast, the resolution in the wall-normal direction is hardly of relevance during the blow-up, such that  = 80 collocation points suffice.For remarks on the general numerical treatment consider Section 2.
From the tracking of the minima of   , we obtain   ≈ 1.0427 and   ≈ 6.213.However, it is not that easy to estimate the other flow-specific constants in the blow-up scalings (11).To find , we consider the expansion (12) for , evaluated at the characteristic minima curve ỹ = ỹ0 ( x).At x = 0, the graph of Ũ0 has a zero crossing in the form of an inflection point.Hence, we track the inflection point of the numerically calculated velocity minimum and find the slipping velocity  ≈ 7 × 10 −3 as the offset.For the remaining scaling constants  0 and , we try a direct comparison between the asymptotically and numerically calculated minima curves (see Figure 5).Since x and ỹ depend strongly on the estimated constants and there are always higher order contributions in expansion ( 12) present, we obtain just order of magnitude estimates for  0 ≈ 2 × 10 −3 and  ≈ 5 × 10 −3 .
Overall, our numerical results and the theory of breakdown are in good agreement.Still, due to an overall lack of accuracy very near to the blow-up caused by the subtle structure of the problem, we cannot rule out that a slightly different The characteristic minima of the mainstream velocity Ũ, (12), value Ũ0 in (A) and location ỹ0 in (B).Comparison between numerics as we approach   (black lines for  ≈ 6.195, 6.201, 6.203, 6.205, 6.206) and the theoretical limit (red).We adjusted  0 and  in the blow-up scalings (11) such, that a convergence is at least imaginable.
form than the asymptotic structure (11), ( 12) is possible.This seems to be one of the reasons why some earlier numerical works, for example, [8,18], are formulated in Lagrangian coordinates, which have other drawbacks.Finally, let us discuss the behaviour of the BL characteristics ( 9) and (10) as we approach the blow-up time   , visualised in the form of a close-up in Figure 6.The displacement thickness  * in Figure 6A builds a very sharp peak, in accordance with the temporal increase of the stream function (14).To get an impression of how sharp the peak really is, we suggest a comparison to the stationary result in Figure 2.This clearly indicates the local breakdown of the BL theory.Since the thickness of the BL locally grows without bound, the wall-normal coordinate ȳ =

√
is no longer of order one.The wall shear stress   in Figure 6B lags a little bit behind but ultimately peaks also at the same mainstream position   .Its negative values clearly show flow reversal, occurring near the lower channel wall of our system.We would like to point out that unsteady flow separation happens slightly later, as it is defined through the BL breakdown, possibly in agreement with the findings of [16].

CONCLUSION
In summary, we studied the breakdown of the classical unsteady boundary layer equations using a well-defined flow scenario.The advantages of the considered channel flow are the free control over its time dependency as well as over the BL evolution towards breakdown, through the suction rate (), (3), and its modulation in time.
The resulting pressure gradient imposed on the BL is given in the analytic form, and we only have to calculate deviations (7) from the flat plate Blasius similarity solution (i.e., from the zero suction scenario  = 0).Our high-resolution Eulerian scheme combines the benefits of finite differences and Chebyshev collocation, which enables us to reach the asymptotic regime of the finite-time blow-up.
While a precise verification of the structure of the Van Dommelen-Shen singularity is hardly possible, due to the presence of small flow-specific constants in (11), our numerical solution strongly supports the characteristic features of the asymptotic structure.These are the temporal pole structure (14) and the overall shape of the generated flow reversal area.
Based on the examined singularity structure, we may continue to the subsequent singular perturbation problem, in our efforts to understand the incipient process of laminar-turbulent transition in BLs.Regarding the blow-up of the next (or first interactive) stage, there are many open questions (see, e.g., [7]), such as the occurrence of instabilities.

F I G U R E 1
Sketch of the channel.The suction  0 causes a velocity perturbation Δ far from the channel inlet/outlet, and a local growth of the boundary layer thickness (circled 2).The fluid in the main area of the channel (circled 1) undergoes a potential flow in leading order.

A
C K N O W L E D G M E N T S This work was supported by the Austrian Science Fund FWF, Grant No. P 31873-N32.O R C I D Markus Kaczvinszki https://orcid.org/0000-0003-2341-6108Stefan Braun https://orcid.org/0000-0002-7145-1103R E F E R E N C E S