Optimal perturbations and transition energy thresholds in boundary layer shear flows

Subcritical transition to turbulence in spatially developing boundary layer flows can be triggered efficiently by finite amplitude perturbations. In this work, we employ adjoint-based optimization to identify optimal initial perturbations in the Blasius boundary layer, culminating in the computation of the critical energy threshold for subcritical transition and the associated fully localized critical optimum in a spatially extended configuration, the so-called minimal seed. By dynamically rescaling the variables with the local boundary layer thickness, we show that the identified edge trajectory approaches the same attracting phase space region as previously reported edge trajectories, and reaches the region faster.

Shear flow transition to turbulence is a complex classical physics problem demonstrating rich spatiotemporal dynamics [1][2][3]. The inherently nonlinear process of subcritical bifurcations and the resulting formation of spatially localized states play a dominant role both in fluid dynamics [4,5] and in other areas of physics [6,7].
The study of subcritical transition in shear flows has seen significant progress via the application of concepts and methods derived from a dynamical systems viewpoint, which naturally arise from the phase space coexistence of two attracting regimes. Such concepts are the edge manifold, which is the invariant set separating trajectories approaching the laminar and turbulent attractors, and edge states, saddle states present on this manifold [8,9]. Another important phase space concept is the minimal seed (or critical optimum), which is the point on the edge manifold that is closest to the laminar attractor in the L 2 (energy) norm [10]. Perturbing the minimal seed infinitesimally along its unstable direction yields the lowest-energy initial perturbation that can trigger subcritical transition to turbulence. This framework has been studied in detail in spatially homogeneous flows, including the identification of the minimal seed for pipe flow [11] and plane Couette flow [12,13].
Previous work has shown that the application of dynamical systems concepts can be fruitful in spatially developing flows too, albeit with increased computational complexity and a more complicated theoretical interpretation [14,15]. The canonical example of spatially developing flows is boundary layers, and specifically the Blasius boundary layer, which approximates the incompressible high-Reynolds laminar flow over a horizontal semi-infinite flat plate [16,17]. No fixed length scale exists and the Reynolds number increases along the streamwise direction x; it is usually defined with the boundary layer's displacement thickness δ * (x) := ∞ 0 (1−U (x, y)/U ∞ )dy, where U denotes the streamwise velocity field, U ∞ the freestream velocity, and y the wall-normal direction. The Blasius solution is unstable to perturbations of sufficiently high amplitude, but becomes linearly unstable too at R δ * 520 [18,19]. Thus, depending on the local Reynolds number, transition can be triggered both subcritically and supercritically. A key difference between the two transition routes is their characteristic timescales, with the subcritical route being faster [15]. Exploiting this timescale separation, we focus on Blasius transition for low-to-moderate R δ * values, where the subcritical route is dominant.
Previous work has attempted to identify optimal initial perturbations for the Blasius boundary layer [20]. Although these authors do indeed find nonlinear optima in a narrow domain, they employ a different definition for the minimal seed. More specifically, they identify a "basic building block" [20] at the core of the found optima for initial energies above a "nonlinearity threshold" [20]. This, however, is not consistent with the above definition. Additionally, no critical threshold for transition was identified in [20].
In this work, we identify for the first time the critical energy threshold and fully localized critical optimum for subcritical transition in a spatially developing boundary layer flow. Moreover, we investigate the edge manifold's dynamics right from its point closest to the laminar solution, and demonstrate that trajectories starting from different regions of the manifold do indeed approach the same attracting region in phase space.
Within the subcritical transition timescales, we use adjoint-based optimization to identify optimal initial perturbations and the critical energy threshold E c for transition. We employ a spatially extended domain, allowing for full localization of the perturbations. We time-march the Blasius solution until a converged state is reached and use it as the baseflow; we consider the evolution of perturbations about this baseflow. For the optimization we employ a well established variational framework [21]. The optimizing functional is the perturbation kinetic energy gain at time horizon T where u(t, x) denotes the perturbation velocity field and · the L 2 norm. The governing equations are the Navier-Stokes and continuity equations. To update the initial field u 0 we employ the gradient-rotation method [22], which ensures that E 0 remains constant without an explicit energy constraint. To relax the high memory demand during the adjoint integration, we incorporate a memory checkpointing routine [23]. The use of state-of-the-art methods is critical in such a computationally demanding undertaking: compared to previous studies [20], we use a 2.5-times longer and 10-times wider domain, as well as a 5-times longer optimization time horizon. Additionally, on average 154 optimization iterations are performed for each reported optimum. Numerical convergence is checked by the residual r n := (δ u0 L n ) ⊥ 2 / δ u0 L n 2 , where L is the problem's Lagrangian [22]. The parameters explicitly appearing in the algorithm are the initial energy E 0 , the time horizon T , and the inflow Reynolds number R in . Due to the spatial variation, the streamwise domain length L x implicitly affects the optimization as well. For fixed R in , the values of T and L x control the range of local Reynolds numbers the flow evolution will span. Therefore, they determine the dynamical behavior present within the optimization horizon.
The streamwise, normal and spanwise directions are denoted by x, y and z, respectively. The presented results are obtained for R in = 275, T = 400 and domain dimensions (L x , L y , L z ) = (500, 60, 100). An inflow fringe region of length l x = 25 is used, with l x being part of L x . The length and velocity scales used for nondimensionalization are the inflow boundary layer displacement thickness δ * 0 (upstream of the fringe) and U ∞ ; the used timescale is thus δ * 0 /U ∞ . The reported R in value is attained right downstream of the fringe. The simulations are done using the open source spectral element code Nek5000 [24]. In terms of boundary conditions, at the inflow we set the velocity to its base value (zero perturbation), while a no-slip condition is imposed at the bottom wall. A stress-free (Neumann) condition is imposed at the outflow. At the top boundary: the streamwise and spanwise velocities are set to their base values (zero perturbation), and a stress-free condition is used for the normal velocity. Periodicity is imposed in the spanwise direction. The initial perturbation field consists of a divergence-free pair of tilted streamwise vortices, so that no symmetries are introduced into the resulting solutions [25]. We do not require transition to be triggered within the time horizon T . To establish whether an optimum triggers transition or not, we perform long-time runs with time horizon t f = 1500, using the volume-averaged streamwise vorticity ω x := (1/V |ω x | 2 dv) 1/2 as an observable, where dv = dxdydz and V is the domain's volume. ω x (t) characterizes the streamwise vortices' strength, and is equal to zero for the baseflow. Transitioning optima present an ever-increasing ω x (t) evolution; non-transitioning ones transiently grow and then decay while slowly relaminarizing.
To calculate E c we follow a bisection process. Starting from two E 0 levels, one leading to subcritical transition, E 0,T , (d) t = 2500; the perturbation is considerably longer, while maintaining its streaky structure; the λ2 contours help identify its active core [14], typical of edge trajectory perturbations. and one to relaminarization, E 0,L , we choose their average value, E 0,i = 1/2(E 0,T +E 0,L ), and perform an optimization for this energy level. If the calculated optimum triggers transition, then E 0,T ← E 0,i ; otherwise E 0,L ← E 0,i . We then set E 0,i+1 in the same fashion and continue the iterative process [26]. Having approached E c with an accuracy of O(10 −3 ), we employ an edge tracking algorithm [8,15,27] to refine it and ensure that the calculated critical threshold and optimum do indeed lie on the edge manifold. The edge tracking is performed using the in-house pseudospectral code SIMSON [28]. Fig. 1a presents the optimization results, with each point corresponding to the identified optimal initial perturbation. Increasing E 0 leads to increased gain, although with different increasing rates for energies below and above E c . For E 0 < E c the gain increases slowly with E 0 , with the rate getting larger as E c is approached. For E 0 > E c the gain increases rapidly with E 0 , signifying that a qualitatively different range of dynamics is now within reach. Crossing the energy threshold E c means that the optima can now access the turbulent dynamics, leading to the attainment of highly energetic endflows via a much more effective extraction of energy from the baseflow.
We identify the critical energy threshold separating the optima undergoing subcritical transition and those that do not, converging to the value E c = 3.644 · 10 −2 ± 2 · 10 −5 . E c is obtained via an edge tracking computation, resulting in the closest edge-bracketing trajectories presenting a separation of 2% at t = 824.6 [29]. In phase space terms, E c characterizes the lowest-energy iso-hypersurface that intersects the edge manifold; the critical optimum is the initial perturbation corresponding to the intersection point. The critical optimum identified in the present configuration is the point on the edge manifold that is closest to the laminar solution at R δ * (x m ) = 425.1, where x m corresponds to the center of mass of the critical optimum's spanwise velocity. Due to the spatial development, E c and the critical optimum depend on the local streamwise coordinate x, since x sets the local Reynolds number R δ * (x). Fig. 1b shows the edge tracking results, with the green full line being the identified edge trajectory; the vorticity curves for the optima for E 0 = 3.5 · 10 −2 , 3.7 · 10 −2 are also shown. Despite their small separation in G T , the optima quickly diverge from one another and from the edge trajectory in their vorticity timeseries. Moreover, their early-time vorticity peaks are reached at different times, a result of their different asymptotic behavior and spatial positioning. Since for E 0 < E c triggering subcritical transition is not possible, the corresponding optimum tries to delay its inevitable viscous decay to achieve the maximal G T possible. To do so, it also positions itself further downstream, where the local R δ * (x) values are higher, leading to more favorable growth conditions. On the other hand, transition can be triggered for E 0 > E c , thus the respective optimal perturbation tries to trigger it as early as possible, in order to take advantage of the turbulent dynamics and attain a higher G T . Due to its higher energy content and its ability to trigger transition, it can also survive relatively lower R δ * (x) values; thus, it positions itself further upstream, leaving bigger part of the domain for its energy growth and spatial spreading. Because of the spatial variation, the streamwise position the localized optima will assume within the domain is also an output of the optimization process, and is not a priori set by the algorithm. It is worthwhile comparing Fig. 1a to the energy-gain diagram found for plane Couette flow (Fig. 3a in [13]). First, the gain discontinuity at E c found in [13] does not appear in the present study. This occurs because we do not require transition to be triggered within the optimization time horizon, so no turbulent endflow needs to be attained. Thus, optima on either side of the edge manifold do not diverge considerably from it, and hence are not so strongly separated in terms of gain at t = T . Second, G T keeps increasing for E 0 > E c . In [13], the final energy E T saturates for E 0 ≥ E c , leading to the decrease of gain for higher E 0 . On the contrary, in our case E T does not saturate, because no turbulent domain-filling endflow is attained; the higher the E 0 , the higher the resulting E T . The gain keeps increasing for E 0 E c because E T grows "more quickly" than E 0 as the latter is increased.
The critical optimum is fully localized, non-symmetric, and consists of alternating areas of positive and negative streamwise velocity tilted against the shear (Fig. 2a). Investigating its evolution, we identify the occurrence of the Orr [30] and liftup [31] energy growth mechanisms. The Orr is a linear mechanism responsible for the early energy growth experienced by the perturbation, t ∈ [0, 60], and can be characterized by the tilting of the structure about the spanwise axis. Fig. 2a-b demonstrate that, although initially leaning backwards, the perturbation has been rotated and is leaning forward at t = 60. The Orr is followed by the liftup mechanism, which acts via streamwise vortices displacing areas of high and low streamwise-momentum fluid in the wall-normal direction. This creates distinct areas of low-and high-speed fluid, leading to the formation of streaks (spanwise variation of streamwise velocity) and overall elongation of the advected structure [32]; the high-and low-speed streaks making up the perturbation can be observed in Fig. 2c-d. Both the tilting of the perturbation caused by the Orr mechanism for t ∈ [0, 60] and the elongation caused by the liftup mechanism for t ∈ [70, 200] can be seen in the short-time simulation movies available in the supplement [33]. The above characteristics of the Orr and liftup mechanisms have been described in greater detail in [34] for plane Couette flow. For longer times, we observe an approximately self-similar evolution, with the perturbation maintaining its streaky structure while growing considerably in size due to the boundary layer's spatial growth, and resembling the edge trajectory identified in [15] (Fig. 2d).
The rescaled edge trajectory is shown in Fig. 3, along with part of the edge trajectory found in [15] using a nonoptimized initial perturbation. Both trajectories approach the same relative attractor in phase space, despite their different initial conditions and early-time behavior. It has been shown in [15] that the corresponding edge trajectory stays in the identified attracting region for times up to t 8000. Moreover, both trajectories form similar loops in the same phase-space region and mirror one another in the last part of their dynamics. These results indicate that memory of the trajectories' initial conditions has indeed been lost [14]. The formation of loops was also identified in [14,15], and stems from a self-sustained physical process alternating between the generation of vorticity and the creation of low-speed streaks [14]. The stages of this process can also be identified by the phases of acceleration and deceleration the advected perturbation goes through, captured in the long-time simulation movie available online. The timestamps provided in Fig. 3a show that the critical optimum needs about half the time to reach the attracting region, compared to the non-optimized initial perturbation; it also requires an approximately 1700-times lower initial energy level, demonstrating its effectiveness in accessing the edge manifold's dynamics.
To summarize, we have used adjoint-based optimization together with edge tracking to calculate optimal initial perturbations for subcritical transition in a spatially developing boundary layer flow, and have identified the critical energy threshold E c = 3.644 · 10 −2 ± 2 · 10 −5 corresponding to the critical optimum at R δ * (x m ) = 425.1. We have observed the main mechanisms responsible for its early-time energy growth, also present in parallel flows [34]. Moreover, the differences in the positioning of the localized optima above and below E c have been explained in terms of the optimization process. Finally, by dynamically rescaling the variables, we have shown that two very different initial conditions on the edge manifold approach the same attracting phase space region already identified for longer times [15]. The critical optimum reaches this region more effectively than non-optimized initial conditions in terms of both initial energy and time needed. Future work includes calculating E c for additional local Reynolds numbers R δ * (x) to reveal the scaling relationship of the two, as has already been done for parallel flow configurations [34].