Experimental and numerical investigations of transition in a pressure-gradient-induced laminar separation bubble

Abstract

A pressure-gradient-induced laminar separation bubble (LSB) was examined using wind-tunnel experiments, direct numerical simulations (DNS) and linear local/global stability analysis. The LSB was experimentally generated on a flat plate using the favourable-to-adverse pressure gradient imposed by an inverted modified NACA $64_3-618$ airfoil. Direct numerical simulation was performed using boundary conditions extracted from a steady precursor simulation of the entire flow field. Despite good agreement in the upstream boundary layer with the experiment, DNS exhibited an approximately 25 % longer mean separation bubble, attributed to an earlier onset of transition due to the free-stream turbulence (FST) in the experiment. Introducing a very low level of isotropic FST in the DNS, similar to that measured in the experiment, caused earlier transition, decreased the mean bubble length and led to a remarkably good agreement between the DNS and experiments. Differences were observed for the dominant frequencies in the experiment and DNS, but both were within the band of most amplified frequencies predicted by LST. Proper orthogonal decomposition confirmed that dominant coherent structures from DNS and experiments are related to the inviscid Kelvin–Helmholtz instability and have similar characteristics despite slight differences in frequency. Local and global stability and dynamic mode decomposition analysis corroborated the convective nature of the dominant two-dimensional (2-D) instability and showed that the LSB is globally unstable to a range of 3-D wavenumbers, in agreement with 3-D structures observed in experiments. Results confirmed the strong impact of very low FST levels on the LSB and indicate a close agreement of the time-averaged and instability characteristics between the experiments and DNS.

1. Introduction

When subjected to a strong adverse pressure gradient (APG), the laminar boundary layer developing in a low Reynolds number, low free-stream turbulence environment can no longer follow the contour of the surface and separates. Downstream of the location of laminar separation, hydrodynamic instability of the inviscid shear layer causes rapid formation of intense spanwise vortical structures (i.e. Rist & Maucher 1994; Diwan & Ramesh 2009; Postl, Balzer & Fasel 2011; Balzer & Fasel 2016). Subsequent three-dimensional deformation of these structures (i.e. Postl et al. 2011; Marxen, Lang & Rist 2013; Michelis, Yarusevych & Kotsonis 2018) leads to breakdown to turbulence, further enhancing mixing and leading to turbulent reattachment of the flow. In a time-averaged sense, this causes a closed recirculation region, called a laminar separation bubble (LSB, Gaster 1967; Horton 1968). The transition process is susceptible to a variety of parameters, such as, disturbances in the approaching boundary layer and free-stream turbulence (FST). Relevant technical applications include, but are not limited to turbo-machinery, laminar airfoils used in uncrewed aerial vehicles and wind turbines. The displacement effects associated with the region of recirculating fluid can drastically reduce the aerodynamic efficiency in practical applications and result in highly unsteady flow fields that can increase noise emissions and structural vibrations. A better understanding of the underlying mechanisms responsible for the transition processes in LSBs is required for the development of effective and efficient flow control strategies for delaying transition in the future.

Laminar separation bubbles have been studied extensively dating back to the seminal work of Gaster (1967) and Horton (1968). Early studies were predominantly experimental and focused on the topology and bursting behaviour of LSBs, in particular with respect to variations in Reynolds number, FST and pressure gradients. In the past decades, experimental efforts shifted towards laminar to turbulent transition within the bubble (Häggmark et al. 2001; Diwan & Ramesh 2009; Radi & Fasel 2010; Simoni et al. 2017; Yarusevych & Kotsonis 2017; Michelis et al. 2018). The experimental studies above as well as numerical investigations (Rist & Maucher 1994; Postl et al. 2011; Marxen, Lang & Rist 2012; Balzer & Fasel 2016) show initial growth of Kelvin–Helmholtz (K–H) instabilities in the separated shear layer, progression of amplitude growth towards nonlinear interaction and subsequent abrupt transition to turbulence (Alam & Sandham 2000; Marxen et al. 2013; Hosseinverdi & Fasel 2018, 2019). In particular Marxen et al. (2009) and Hosseinverdi & Fasel (2018, 2019) show a highly unstable shear layer sensitive to small disturbances within the incoming flow, i.e. in the upstream boundary layer and in the FST. They found amplified instability waves rapidly reach large (nonlinear) amplitudes in the separated region. Time-resolved data from simulations and experiments show periodic shedding of spanwise coherent (two-dimensional) vortical structures (Postl et al. 2011; Balzer & Fasel 2016; Simoni et al. 2017; Istvan & Yarusevych 2018; Hosseinverdi & Fasel 2019). Spanwise modulation of these structures is a result of secondary instability mechanisms (Balzer & Fasel 2016; Michelis et al. 2018; Hosseinverdi & Fasel 2019). The subsequent breakdown into small-scale turbulence was observed upstream of the mean reattachment location (figure 1).

Figure 1. Schematic of a pressure-gradient-induced LSB on a flat plate.

Instability mechanisms are generally classified as either ‘convective’ or ‘absolute/global’ Huerre & Monkewitz (1990). Several attempts have been made to determine a criterion for the onset of global instabilities in LSBs based on the work related to free shear layers following Huerre & Monkewitz (1985). The suggested threshold for the onset of global instability was on the maximum local reverse flow $u_{rev}$ magnitude in the separated region of the LSB, $u^*_{rev} = u_{rev}/U_\infty$ , normalized by the free-stream velocity, $U_\infty$ . Based on two-dimensional (2-D) shear-layer profiles in the presence of a wall, values are in the range $u^*_{rev} \approx 17{-}25\, \%$ (Allen & Riley 1995; Hammond & Redekopp 1998; Diwan & Ramesh 2009), in agreement with direct numerical simulation (DNS) results by Alam & Sandham (2000) and Fasel & Postl (2006). However, unsteadiness and 3-D breakdown were observed in simulations of 3-D LSBs for reverse flow as low as $u^*_{rev} \approx 7\, \%$ (Spalart & Strelets 2000; Hosseinverdi & Fasel 2013; Rodríguez et al. 2013; Rodríguez et al. 2021), causing transition in the absence of external disturbances and casting doubt on the validity of the suggested threshold. Rodríguez et al. (2021) found that flows with maximum reverse flow velocities above $u^*_{rev} \approx 7\, \%$ can lead to 3-D absolute instability of the K–H waves and thus cause a self-sustained transition in the LSB. While these results can explain the lower threshold for the absolute instability, it remains an unresolved topic at this time, especially regarding the influence of external disturbances due to FST.

In contrast to linear stability calculations and numerical simulations, wind-tunnel experiments are never void of FST. Even at low levels, FST has a significant impact on the transition process in both experiments (Klebanoff 1971; Simoni et al. 2017; Jaroslawski et al. 2023) and simulations (Jacobs & Durbin 2001; Balzer & Fasel 2016; Istvan & Yarusevych 2018; Hosseinverdi & Fasel 2019). For low levels of FST ( $Tu \lesssim 0.5\, \%$ ), the disturbance development due to the K–H-type instability in the shear layer remains largely two-dimensional. Increased levels of FST cause the formation of streamwise elongated streaks inside the boundary layer as a consequence of the so-called Klebanoff modes (K-modes, Klebanoff, Tidstrom & Sargent 1962; Klebanoff 1971; Kendall 1985, 1990; Meitz 1996; Fasel 2002; Marxen et al. 2013; Hosseinverdi & Fasel 2017a ). The K-modes cause significant distortion in the streamwise velocity component in the spanwise and wall-normal directions. Contrary to Tollmien–Schlichting waves, K-modes are characterized by their low frequency, low spatial growth rates and spanwise spacing of a few boundary-layer thicknesses. Linear stability investigation and comparison with DNS discussed in Hosseinverdi & Fasel (2019) indicate that slow (algebraic) growth of the K-mode is followed by a strong exponential amplification in the region of APG. The growth rates observed for the 2-D instability waves, attributed to the K–H instability in the shear layer yield significantly larger initial growth rates than the K-modes. The combined effect of these two instability mechanisms, K–H and K-mode, is that in the presence of FST transition is accelerated, leading to a smaller mean separated region. Increased FST amplifies this behaviour and eventually leads to bypass transition (Jacobs & Durbin 2001; Jaroslawski et al. 2023).

In this paper, experimental and numerical (DNS) investigations of LSBs on a flat plate are presented. The LSBs were generated by an inverted wing, and results are supported by linear stability analysis (local and global). The focus of the investigation is on (i) ‘natural’ transition in a low FST environment in both the wind tunnel and DNS and (ii) extensive, physics-based comparisons of the time-averaged flow and stability characteristics in the LSB between the experiment, stability analysis and DNS with and without levels of FST comparable to those observed in the wind-tunnel experiments. Detailed comparisons of this type (experiments, DNS and stability analysis) are uncommon, especially when considering the effects of FST. The current work explores the unforced conditions and the effects of FST, while providing a physics-based comparison between stability analysis, DNS and experiment for both the time-average and unsteady flow. Local and global stability analyses address the instabilities in the flow relevant to the ’natural’ transition. Direct numerical simulation and experiment aim to confirm that the same instability mechanisms are present.

2. Windtunnel experiments

2.1. Experimental set-up and instrumentation

The experiments were conducted in the Arizona Low Speed Wind Tunnel (ALSWT) situated in the Department of Aerospace and Mechanical Engineering at the University of Arizona. The closed-loop wind tunnel has a test section of 0.9 m $\times$ 1.2 m $\times$ 3.65 m (height x width x length, 3 ft $\times$ 4 ft $\times$ 12 ft). A Pitot tube mounted 0.4 m downstream of the test section entry at the tunnel sidewall extends into the free stream to acquire total and static pressures to measure the flow speed. The maximum flow speed in the test section is $U_{\infty } \approx 80\,\mathrm {m\,s^-{^1}}$ (262.5 ft s^–1). Uniformity of the mean flow over the test section is at or better than $\pm 0.5\,\%$ . The turbulence intensity is less than Tu = 0.05 % in the range of 1 Hz to 10 kHz and considerably lower (e.g. Tu = 0.035 %) for the flow speed of the experiments discussed in this paper (7 m s^–1, see § 2.2). A more extensive investigation of the free-stream properties can be found in Borgmann et al. (2020). The temperature inside the tunnel is regulated by a heat exchanger with a chilled water supply. Throughout the experiments, temperature was held within the range of $\pm$ 0.44 $^{\circ }$ C (1 $^{\circ }$ F) of 22.2 $^{\circ }$ C (72 $^{\circ }$ F).

A LSB is generated on a flat plate by an inverted wing: a NACA $64_3-618$ airfoil with a chord length of $c = {8}\,\mathrm {in.}$ (203.2 mm) shown in figure 2. The favourable-to-adverse pressure gradient along the suction side of the inverted wing leads to the formation of a LSB for a wide range of Reynolds numbers. Projection of the pressure gradient onto a flat plate neglects the influence of surface curvature to the transition process, analogous to early work by Gaster (1967). In addition to the flexibility in the mounting location and incidence angle of the inverted wing, the relatively small streamwise extent and high aspect ratio provide a quasi-2-D LSB across a large spanwise extent (more than four times the time-averaged separation length). In addition, the low total blockage in the wind tunnel provides a swift pressure recovery downstream of the LSB. This contrasts with the large extent of APG generally observed downstream of reattachment when using a contoured ceiling (Watmuff 1999; Diwan & Ramesh 2009). Varying the distance and the incidence angle of the wing relative to the flat plate allows for a wide range of possible free-stream pressure distributions acting on the flat-plate boundary layer. For the investigation discussed here the inverted wing is placed at a distance of 3.2 in. (81.3 mm) from the plate surface; 17.5 in. (0.44 m) downstream of the leading edge; at an incidence angle of ${2}^{\circ }$ (relative to the inverted airfoil, see drawing in figure 2). The wing is equipped with boundary-layer suction along $50\,\%{-}60\,\%$ -chord and $85\,\%{-}90\,\%$ -chord, on the suction side and a zigzag turbulator tape to trip the flow close to the leading edge at the pressure side to avoid separation on either surface of the inverted wing. Careful distribution of suction holes in a PVC pipe inside the plenum of the hollow downstream half of the inverted wing ensure even suction across the span. Suction is provided by LAMB AMETEK 117500-12 3-STAGE vacuum motor, connected to both ends of the inverted wing.

Figure 2. Flat-plate model and inverted airfoil used for LSB transition experiments: (a) isometric view showing coordinate system (red), pressure tap locations (orange) and suction direction (blue); (b) cross-section of airfoil showing PVC pipe outline, plenum in the hollow wing and holes along the suction side; planform view shows suction distribution along the airfoil suction side where red lines highlight the active sections; (c) side view, highlighting the incidence angle between the inverted wing and flat plate (blue text); Coloured planes in (a) and (c) mark the PIV planes in the x-y (green) and x-z planes (orange); The origin of the coordinate system shown in (a) and (c) places $z = 0$ at the centre of the span and $y = 0$ at the flat-plate surface.

The flat-plate model shown in figure 2 is made from an aluminium honeycomb structure, which is covered on both sides by a 1 mm aluminium skin. The width of the plate is 47.75 in. (1.21 m), spanning nearly the entire width of the test section. Gaps between the sides of the plate and the walls of the test section were sealed with a soft expanding foam. The plate is 78 in. (1.98 m) long with a total thickness of 0.5 in. (13 mm) and a flatness of $\pm {0.008}\,\rm {in\, ft^-{^1}}$ ( $\pm {0.665}\,\rm {mm\, m^-{^1}}$ ). Streamlined support structures elevate the plate to 8.25 in. (210 mm) above the wind-tunnel floor to avoid wall effects. The interchangeable leading edge is configured as a 1 : 20 super-ellipse, see Lin, Reed & Saric (1992). Its total length of 9.5 in. (241 mm) and smooth transition to the plate top surface combined with the adjustable flap at the trailing edge of the plate ensure consistent and repeatable inflow conditions, as the flap compensates for blockage and allows for positioning of the stagnation point along the leading edge. The flow along the bottom side of the plate is tripped just downstream of the leading edge to reduce unsteadiness in the flow underneath the plate. The origin of the coordinate system used throughout the investigation is located at the centre of the leading edge, with the x axis in the streamwise, y axis in the vertical and z axis in the spanwise direction. The plate surface is considered zero in the vertical (y) direction. All results in this work are presented in dimensionless form where the length scales, $x$ , $y$ and $z$ , are scaled with a reference length of $L_\infty ^\ast = {0.0254}\,\mathrm {m}$ ( ${1}\,\mathrm {inch}$ ), in accordance with the earlier experiments by Gaster (1967) and velocities are scaled with the free-stream velocity $U_\infty ^\ast$ , respectively

Pressure taps are located along two parallel lines in the streamwise direction, 1.5 in. (38 mm) offset to either sides of the centreline. The streamwise spacing is 12.5 mm (0.5 in.) in the vicinity of the LSB and 25 mm (1 in.) elsewhere. The spacing is gradually refined towards the leading edge with the first tap at 15 mm (0.6 in.) and an initial spacing of 5 mm 0.2 in.). Scanivalve Corp. ZOC33 pressure scanners (full scale range [FD] and accuracy of 2.49 k Pa (10 in. $\text {H}_2\text {O}$ ) $\pm 0.15\, \%$ FS) in combination with an ERAD Remote A/D module are used to record surface pressure. Pressure data are sampled at 625 Hz and averaged over 96 s.

Time-resolved velocity measurements were recorded using constant temperature anemometry (CTA). A Dantec Dynamics StreamlinePro frame was operated with a 55P11 miniature wire probe and a 55H21 support. The hot-wire probe was attached to a two axis traverse along the wind-tunnel ceiling providing variable positioning of the sensor along the x-y plane. The CTA probe was calibrated using a StreamLine Pro Automatic Calibrator and referenced to a temperature probe positioned above the CTA in the free stream on the traverse.

The CTA data in the free stream and along the LSB were recorded at 40 kHz for 20 s and analogue filtered between 1 Hz and 10 kHz. Spectral content was calculated using Welch’s method, with a resolution of 0.61 Hz averaged over 21 windows. Free-stream turbulence intensity was calculated as the root-mean-square of the velocity fluctuations scaled by the local free-stream velocity (2.1) over the same frequency range (1 Hz–10 kHz).

(2.1) \begin{equation} Tu = \frac {\sqrt {|u'|^2}}{U_\infty } .\end{equation}

In addition to hot-wire data, a LaVision GmbH particle image velocimetry (PIV) system was used to obtain spatially resolved velocity data in x-y planes of the wake of the inverted wing and in the LSB (PIV plane shown in figure 2). A Quantel Evergreen HP dual-head Nd:YAG laser illuminates submicron di-ethyl-hexyl-sebacat (DEHS) seed particles. In the wake of the inverted wing, 300 image pairs were acquired at 10 Hz with a 5-megapixel 16-bit Imager sCMOS camera, incorporating a 2x teleconverter lens, 50 mm lenses and narrow band-pass optical filters. For each preprocessed and time-filtered image pair, subregions were cross-correlated using decreasing window sizes ( $64^2$ – $24^2$ pixel) and multipass processed with $50\, \%$ overlap followed by a median filter and smoothing function. To increase the resolution in the LSB, two cameras with the same optics set-up were placed in parallel. The data in the x-y plane along the LSB consist of 1500 image pairs, recorded at 14 Hz. Images were preprocessed and time filtered, and subregions were cross-correlated with window sizes decreasing from $64^2$ – $16^2$ pixel with $50\, \%$ overlap. The final images were stitched and smoothed by a single pass with a median filter (size 3 × 3). The error based on a 95 % confidence interval for the velocity magnitude relative to the free-stream velocity (7 m s^–1) is ${\lt}3\,\%$ in the wake of the inverted wing and ${\lt}2\,\%$ in the separation bubble. The larger number of images recorded for the time-averaged data in the separation bubble is motivated by convergence of second-order statistical quantities (e.g. Reynolds stress, turbulent kinetic energy, etc.). Assuming a normal distribution of all statistical quantities the relative error in the turbulent kinetic energy based on a 95 % confidence interval is ${\lt}5.5\,\%$ in the LSB.

Stereoscopic PIV was used in the x-z plane of the LSB. Two 5-megapixel 16-bit Imager sCMOS cameras incorporating scheimpflug mounts, 85 mm lenses and narrow band-pass optical filters are used to record 800 image pairs at 14 Hz. Using the same processing routine described for the x-y plane data, image pairs are processed separately, using a final window size of $32^2$ pixel. The resulting velocity fields are further postprocessed to remove vectors with correlation peak ratios of less than 1.5 and correlation coefficients below 0.1. After applying the stereo calibration, a 3 $\times$ 3 Gaussian smoothing algorithm and polynomial filters were employed for the final data. The time-averaged results have less than 3.5 % statistical uncertainty (based on free stream) using a 95 % confidence interval. The spatial resolution is 1.6 mm (0.064 in.).

2.2. Free-stream conditions

The free-stream flow in the test section with and without the model (airfoil and flat plate) installed was measured upstream of the location of the inverted wing (when installed), at the mid point (vertically) between the flat plate and tunnel ceiling in both cases ( $x = 9$ , $y = 13$ , $z = 0$ ). The tunnel operating conditions were identical with respect to fan frequency and fan blade angle, resulting in a slight increase in free-stream velocity ( ${\approx}4\, \%$ ) when the model was installed due to the blockage effect of the airfoil and flat plate. This generates a slightly elevated noise floor in the PSD with the model installed. Prior to mounting the model in the wind-tunnel test section, assessment of the free-stream of the empty test section showed low levels of turbulence intensity (Tu $\leqslant 0.034\,\%$ , in the range of 1 Hz to 10 kHz). Data were collected at a Reynolds number of $Re_C = 90\,000$ , based on the chord of the inverted wing ( $c = {203.2}\,\mathrm {mm}$ ). Major contributions to the turbulence intensity are observed in the low frequency range ( $\lt {5}\,\mathrm {Hz}$ ), resulting in Tu $\leqslant 0.015\,\%$ if the signal is band pass filtered between 5 Hz to 10 kHz). Therefore, the PSD of the streamwise velocity is shown for reference in figure 3. With the flat-plate model and inverted wing installed, a minimal increase to Tu $\leqslant 0.035\, \%$ is observed in the range of 1 Hz to 10 kHz (or Tu $\leqslant 0.017\, \%$ in the range of 5 Hz to 10 kHz), while spectral content remains nearly unchanged. Peaks in the spectrum at 36 Hz and 72 Hz can be attributed to the fan blade passing frequency and its harmonic.

Figure 3. Power spectral density (PSD) of the free-stream velocity in the empty test section and at test conditions for the LSB with the model installed, $Re_C = 90\,000$ .

2.3. Flow around the inverted wing

At low chord based Reynolds number ( $Re_C = 90\,000$ ) of the inverted wing and low incidence angles ( ${2}^{\circ }$ ) between the airfoil and the flat plate (figure 2), laminar separation was observed at the suction (bottom) and pressure (top) side of the NACA $64_3{-}618$ airfoil. Velocity magnitude contours in figure 4 show the necessity of flow control to avoid separation leading to unsteadiness in the free stream. Without control the flow separates on the wing just downstream of the maximum thickness at $x = 21.45$ ( $35\, \% c$ ) on the suction side facing the flat plate, figure 4(a). The boundary layer does not reattach to the wing and causes a significant wake. In addition, a separation bubble forms on the pressure side between $x = 22.05$ and $x = 24.8$ (figures 4 a and 4 c). The same was observed in precursor computational fluid dynamics (CFD) calculations shown in figures 4 b and 4 d. Steady 2-D boundary-layer suction through the mechanism described in § 2.1 was successful in preventing separation on the suction side. The method is different from the jet entrainment scheme at the trailing edge of the inverted wing used in the seminal work by Gaster (1967, 2006). Successful application of quasi-continuous boundary-layer suction (plenum suction behind a perforated section of the wing) in a similar set-up in a water tunnel experiment was previously shown by Radi & Fasel (2010) and Jagadeesh & Fasel (2013). Figures 4(c) and 4(d) demonstrate the effectiveness of suction to suppress boundary-layer separation as the flow remained attached along the entire chord.

Figure 4. Velocity magnitude $|U|$ along the inverted wing in the experiment (a, c, e) and from the base flow calculations in ANSYS Fluent withthe $\gamma -Re_\theta -SST$ transition model (b) and (d); (a) experiment with no boundary-layer suction and no trip on pressure side; (b) ANSYS fluent without suction; (c) experiment with suction, but no trip on pressure side; (d) ANSYS Fluent calculation with suction; (e) experiment with suction and pressure side tripped; (f) PSDs from measurement in the wake of the pressure side with and without trip ( $x = 25.5$ ; $y = 4.75$ ).

The small separation bubble at the pressure side had negligible influence on the time-averaged flow field between the wing and flat plate. Spectral analysis of CTA measurements taken above the trailing edge downstream of the separated region reveals two frequency peaks between 200 and 300 Hz (figure 4 f). To avoid interference between these periodic disturbances emanating from the pressure side of the inverted wing and the LSB on the flat plate in the experiment, the boundary layer on the pressure side was tripped just downstream of the leading edge. Then the turbulent boundary layer no longer separated, the frequency peaks in the power spectra were suppressed (figure 4 f) and the separation on the pressure side was no longer present in the time-averaged flow field along the inverted wing (figure 4 e). With boundary-layer tripping (pressure side) and suction (suction side), the flow remained attached along the entire chord on both sides of the wing (figure 4 e).

2.4. Laminar separation bubble on the flat plate

The pressure-gradient-induced LSB at the flat plate was recorded for the flow field shown in figure 4(e), for which the application of flow control to fully attach the flow along the entire inverted wing. Pressure measurements on the flat plate in figure 5 confirm the effects of suction at the wing to the surface pressure at the flat plate. In figure 5 the downstream $C_p$ development corresponds to figures 4(a) and 4(e), without and with flow control (suction and trip). As expected, the separation without flow control at the wing (figure 4 a) result in a weak favourable-to-adverse pressure gradient in the time-averaged sense. With suction, a significant plateau of constant $C_p$ between $x = 22.5$ and $x = 26$ indicates the presence of a LSB. The effect of boundary-layer tripping at the pressure side of the inverted wing was found to be negligible in the time-averaged pressure measurements (not shown here). The current experimental data compare very well with the pressure gradients in the experiments by Gaster (1967) (Series I, Case 6 in figure 5) for a similar free-stream velocity 21.8 ft s^–1 (6.64 m s^–1)) compare well with the current experiment. To compare the pressure gradients of the current work with results from Gaster (1967), the boundary layer along the flat plate was tripped, thus resulting in a quasi-inviscid pressure development. The resulting turbulent boundary layer stays attached and surface pressure measurements resemble the inviscid case (red and grey line in figure 5). In both cases (laminar and turbulent boundary layer along the flat plate), the NACA $64_3-618$ causes a more gradual favourable pressure gradient compared with the inverted wing used by Gaster (1967). However, the APG and separation location match remarkably well, with a slightly larger LSB in the current study. Results from Gaster (1967) (Series I, Case 6 in figure 5) were scaled in the streamwise direction to compensate for the smaller chord length of $c_{Gaster} = {139.7}\,\mathrm {mm}$ (5.5 in., chord-based Reynolds number of $63\,000$ ). The boundary-layer momentum thickness at separation was lower in Gaster’s experiments, $Re_{s,\theta }=218$ (Series I, Case 6 – 21.8 ft s^–1 (6.64 m s^–1)).

In figure 5(b) spanwise pressure measurements at the flat plate are provided for select downstream locations, for the most unsteady regions close to separation and reattachment of the LSB. The results in figure 5(b) indicate small variations support the assumption of essentially 2-D conditions in a time-averaged sense.

Figure 5. Flat-plate surface $C_p$ development in the experiment. (a) Chordwise $C_p$ with and without boundary-layer suction at the inverted wing, pressure side tripped in both cases. Comparison with Series I, Case 6 (Gaster 1967), for quasi-inviscid (tripped flat-plate boundary layer). (b) Spanwise $C_p$ at select chordwise locations for suction and tripping at the inverted wing.

An overview of time-averaged quantities for the LSB is provided in figure 6. Normalized velocities ( $u/U_{\infty }$ and $v/U_{\infty }$ ) are shown next to the root-mean-square (r.m.s.) of the respective fluctuating quantities. In the region of APG, the incoming laminar boundary layer separates at $x=21.8$ . Mean reattachment was found at $x=27.6$ resulting in a bubble length of $L=5.8$ and a maximum bubble height of $H=0.35$ at $x = 26$ . Separation and reattachment are determined based on the dividing streamline extracted from the time-averaged flow field (figure 1). Velocity vectors along the wall ( $y\lt 0.04$ ) were omitted to avoid the impact of laser reflections from PIV. Vector fields were interpolated to ensure no-slip conditions at the wall. The dividing streamline is identified as the streamline between the time-resolved separation and reattachment location, closing the separated region (Kurelek et al. 2020). The fluctuating quantities in figure 6 were obtained from the r.m.s. of the $u'$ and $v'$ -velocity fluctuations. Low levels of streamwise velocity fluctuations near the separation location and along the dividing streamline upstream of maximum thickness agree with results in Zaman, Mckinzie & Rumsey (1989); Michelis & Kotsonis (2016) and appear to be related to low frequency ’flapping’, which was observed in separation bubbles at similar Reynolds numbers. Significant levels of streamwise and wall-normal fluctuations near the maximum bubble height and downstream of it are related to formation of strong periodic vortex shedding, as a result of the inviscid instability in the shear layer developing between the free stream and the region of reverse flow. Based on the bubble topology, the two parameter bursting criterion proposed by Gaster (1967) was calculated as $\overline {P} = -0.206$ , based on the momentum thickness at separation $\theta = 0.024$ (0.613 mm), the dynamic viscosity $\nu = 1.658 \times 10^{-5}$ , the velocity decrease across the separation ( $\triangle U = {1.335}\,\mathrm {m\,s^-{^1}}$ ) and the separation length ( $L_s$ ). The resulting value is indicative of a short LSB and is almost identical to the LSB in Gaster (1967) Series I, Case 6 as shown in figure 5 ( $\overline {P}_{Gaster,I,6} = -0.199$ ) at similar free-stream conditions.

Figure 6. Time-averaged velocity components (a) $U/U_{\infty }$ and (b) $V/U_{\infty }$ , r.m.s. of velocity fluctuations (c)  $u'_{rms}/U_{\infty }$ and (d) $v'_{rms}/U_{\infty }$ across the separation bubble derived from PIV. Streamline overlay for reference. Dividing streamline (dashed line), separation ( $\blacktriangleleft$ ) and reattachment location ( $\blacktriangle$ ).

3. Numerical/theoretical investigations

3.1. Set-up for stability investigations and DNS

The numerical investigations, consisting of stability analysis and DNS, were carried out to complement the experimental work in order to extract and better explain the relevant physics. The set-up was guided by the wind-tunnel experiments, as discussed in § 2. The computational approach involves two steps (see figure 7 a): in step $I$ , a ‘precursor calculation’ of the entire flow field based on the set-up of the wind-tunnel experiments including the wind-tunnel walls, inverted wing and the top of the flat plate was performed using ANSYS Fluent. The ‘precursor calculation’ uses a transitional $\gamma -Re_{\theta }-SST$ (shear stress transport) model and a slip wall along the suction side of the inverted wing. In step $II$ , high resolution DNS of the flow along the flat plate downstream of the leading edge was performed. Inflow and free-stream boundary conditions for the DNS are extracted from the flow field computed in the precursor calculations.

Figure 7. Schematic of the computational set-up for precursor calculations, step $I$ , according to ALSWT experiments (a). Also shown are contours of streamwise velocity. The area indicated by dashed lines corresponds to the computational domain for DNS, step $II$ (b).

The computational domain for the presented stability analysis and DNS was defined as $10\leqslant x\leqslant 36$ (figure 7 b), the domain height is $l_y=2.8$ and the domain width in the spanwise direction is $l_z=2$ . At the outflow boundary, $x=x_o$ , all second derivatives in the streamwise direction are set to zero. In addition, a buffer domain in region $x_b\lt x\lt x_o$ smoothly dampens out the fluctuations generated inside the domain and prevents upstream contamination (figure 7 b). No-slip and no-penetration conditions are enforced on the surface of the flat plate.

For all DNS results presented in this paper, the same computational grid with $n_x\times n_y\times n_z=2601\times 200\times 220\approx 115$ million grid points was used. The grid spacing is uniform in the streamwise direction, and the grid points are exponentially displaced from the wall in the wall-normal direction to improve the resolution near the wall. In the spanwise direction, 139 Fourier modes (220 collocation points) were employed for the domain width of $l_z=2$ . The spanwise domain width was selected based on the maximum bubble height (Jones, Sandberg & Sandham 2008; Balzer & Fasel 2016; Hosseinverdi & Fasel 2019). It was found that, to resolve the largest spanwise structures, a domain width of at least four times the maximum bubble height is necessary. The ratio of domain width to the maximum bubble height was found to be approximately six. The grid resolution in wall units within the separation bubble and in the redeveloping turbulent boundary layer downstream of the reattachment point was then $\triangle x^+\leqslant 5.8$ , $\triangle z^+\leqslant 5.2$ and $\triangle y^+_w\leqslant 0.5$ .

Direct numerical simulation of transition requires numerical methods with low numerical dispersion and dissipation errors. An extensively validated high-order accurate Navier–Stokes solver developed in the University of Arizona CFD Laboratory was employed for the present numerical simulations (Meitz & Fasel 2000; Hosseinverdi, Balzer & Fasel 2012; Hosseinverdi & Fasel 2017b ). In this code, the 3-D unsteady incompressible Navier–Stokes equations are solved in vorticityvelocity formulation

(3.1) \begin{align} \frac {\partial \boldsymbol {\omega }}{\partial t} &=\boldsymbol {\nabla }\times (\boldsymbol {u}\times \boldsymbol {\omega })+\frac {1}{Re}\triangle \boldsymbol {\omega }, \end{align}

(3.2) \begin{align} \triangle \boldsymbol {u}&=\boldsymbol {\nabla }\times \boldsymbol {\omega }, \end{align}

where the vorticity vector, $\boldsymbol {\omega }$ , is defined as the negative curl of the velocity vector $\boldsymbol {\omega }=-\boldsymbol {\nabla }\times \boldsymbol {u}$ . In (3.1), the global Reynolds number is defined as $Re=U_\infty ^\ast L_\infty ^\ast /\nu ^\ast$ , where $\nu ^\ast$ is the kinematic viscosity. The asterisk denotes dimensional quantities. Here, $U_\infty ^\ast$ and $L_\infty ^\ast$ are reference velocity and length scales, respectively.

The governing equations are solved in a 3-D Cartesian coordinate system, where the streamwise, wall-normal and spanwise directions are denoted by $x$ , $y$ and $z$ , respectively (figure 7), and the corresponding velocity and vorticity components, ( $u$ , $v$ , $w$ ) and ( $\omega _x$ , $\omega _y$ , $\omega _z$ ). The governing equations are integrated in time using an explicit fourth-order RungeKutta scheme. All derivatives in streamwise and wall-normal directions are approximated with fourth-order accurate compact finite differences. Note that the finite-difference approximations for the derivatives with respect to $y$ are constructed for a non-equidistant grid instead of using a coordinate transformation. The flow field is assumed to be periodic in the spanwise direction, therefore, it can be expanded in Fourier cosine and sine series with a pseudo-spectral treatment of the nonlinear terms. The velocity Poisson equations are solved by a combination of a fourth-order standard compact difference in the wall-normal direction and Fourier sine transform in the streamwise direction. For details of the numerical method see Meitz & Fasel (2000), Hosseinverdi et al. (2012)and Hosseinverdi & Fasel (2017b ).

3.2. Direct numerical simulations without free-stream turbulence

Contrary to the wind-tunnel environment, without the addition of FST and other artificial disturbances in the computational domain, the low numerical noise of the high-fidelity DNS leads to quasi-zero FST conditions in the calculations, providing ideal conditions to investigate the natural transition in the LSB. Time-averaged velocity profiles from the DNS without FST in figure 8(a) show the streamwise component in the boundary layer along the flat plate. The results are in good agreement with those extracted from PIV data in the experiment (figure 8 a), suggesting near identical boundary-layer development upstream of laminar separation between experiment and numerical calculations in a time-averaged sense, independent of the difference in FST. As expected, the APG decelerates the flow in the boundary layer downstream of the apex of the inverted wing ( $x = 20.42$ ), resulting in an inflection point in the boundary layer in the vicinity of laminar separation. Experimental data outside the symbols in figure 8(a) were interpolated, enforcing no-slip conditions at the wall.

Figure 8. (a) Normalized boundary-layer profiles between the onset of APG and laminar separation, extracted from PIV and DNS. (b) Instantaneous flow visualizations for 3-D DNS for $Tu=0\, \%$ . Shown are contours of spanwise vorticity (averaged in spanwise direction) together with the mean dividing streamline identified by black dashed lines.

Instantaneous contours of spanwise vorticity (averaged in the $z$ -direction) obtained from 3-D DNS without FST are presented in figure 8(b). The location and size of the bubbles were obtained from the time- and spanwise-averaged data to extract the dividing streamlines in figure 8(b). The inflection of the velocity profile in the shear layer between free-stream and recirculating flow causes growth of disturbances upstream of the maximum bubble height due to the K–H instability. Subsequent ‘roll-up’ and shedding of vortical structures were observed near the location of the maximum bubble height. The shear-layer roll-up produces large wall-normal momentum which facilitates entrainment and aids in the reattachment process. Small vortical structures that are generated near the wall due to the breakdown of these large-scale structures further contribute to the reattachment process. It should be noted that transition to turbulence appears to be self-sustained for the DNS even without FST, i.e. no external forcing or disturbances were required to initiate or sustain the transition process. This observation agrees with the bounded range of globally unstable spanwise wavenumbers identified in § 3.3. This is in contrast to zero-pressure-gradient boundary layers and implies that the development of 3-D disturbances is a direct result of an absolute/ global instability mechanism.

Comparison of the mean flow via streamlines and time-averaged surface pressure coefficient shows the importance of FST for the topology of the LSB (9). The development of $C_p$ calculated using the transitional $\gamma -Re_{\theta }-SST$ model in the Reynolds-averaged Navier–Stokes (RANS) calculations which is known to struggle with LSBs in low Reynolds number flows, predicts – as expected – a shorter bubble than the experiments. The high-fidelity DNS is in good agreement with the experiment for $x\leqslant 25.5$ . However, streamlines and surface pressure indicate that reattachment occurs farther downstream compared with the experiments. The presence of minute ‘numerical noise’ in the DNS (with no added FST) reduces the amplitudes of the disturbances in the upstream boundary layer. The dominant 2-D mode in the shear layer requires extended growth downstream before reaching critical amplitudes (compared with the case with FST), necessary for the interaction with the global mode in the LSB to cause transition to turbulence. This leads to a later transition along the shear layer and results in significant deviation from the wind-tunnel experiments for $x \geqslant 25.5$ , despite the low levels of FST observed in the wind-tunnel experiments (§ 2.2). The influence of FST in the transition process, especially within the LSB, is well documented (Jacobs & Durbin 2001; Balzer & Fasel 2016; Hosseinverdi & Fasel 2019; Simoni et al. 2017; Istvan & Yarusevych 2018) and results in this section corroborate that meaningful comparison of numerical results with experimental efforts will have to account for the disturbances present in the respective wind-tunnel or flight environment (§§ 3.4 and 4).

Figure 9. Comparison of mean results between experiments and DNS. Plotted in (a) are streamlines obtained from experiments and DNS without FST. (b) Comparison of flat-plate surface pressure coefficient, $C_p$ , between experiment, DNS without FST as well as RANS calculations.

3.3. Linear local and global stability analyses

In addition to DNS, local and global linear stability analyses were carried out to investigate the initial stages of the transition process, and to help interpret observations from the experiments and DNS (§ 4) such as for example the self-sustained transition in the DNS without FST as discussed in the previous section. First, linear stability theory (LST) analysis was performed based on the mean flow field (time- and spanwise averaged) obtained from DNS without FST (§ 3.2) in order to map out the stability behaviour with respect to 2-D disturbance waves. The stability diagram from the solution to the OrrSommerfeld equation for the mean streamwise velocity profiles (see figure 8 a) is presented in figure 10, which shows contours of the spatial growth rates, $\alpha _i$ , in the $St$ -x plane. Since the spatial stability problem was considered ( $\alpha$ -complex), the local velocity profiles are either convectively stable ( $\alpha _i\gt 0$ ) or unstable ( $\alpha _i\lt 0$ ). The solid black contour line corresponds to the neutral curve ( $\alpha _i=0$ ). LST results reveal that the LSB is unstable with respect to disturbances in a broad range of frequencies.

Figure 10. Linear stability theory result for the time- and spanwise-averaged $u$ -velocity of the DNS without FST. Shown are the contours of the amplification rate, $\alpha _i$ , together with contour lines of $N$ -factor for 2-D distributed waves in the x- $St$ plane. The solid contour line represents the neutral curve. The dominant shedding frequency observed in the experiment ( $St_{EXP}$ ) and DNS ( $St_{DNS}$ ) are indicated by horizontal dotted lines. Here, $x_{apg}$ and $x_s$ correspond to the $x$ -location of the onset of APG and separation point, respectively.

For the linear global stability analysis (LGSA), the linearized Navier–Stokes equations (LNSE) with a mean flow obtained from DNS were directly solved as an initial value problem

(3.3a) \begin{align} \frac {\partial \boldsymbol {\omega }^{\prime }}{\partial t}&={\nabla }\times \left (\boldsymbol {u}^{\prime }\times \boldsymbol {\Omega } + \boldsymbol {U}\times \boldsymbol {\omega }^{\prime }\right )+\frac {1}{Re}\triangle \boldsymbol {\omega }^{\prime }, \end{align}

(3.3b) \begin{align} \triangle \boldsymbol {u}^{\prime }&={\nabla }\times \boldsymbol {\omega }^{\prime }, \end{align}

where $\boldsymbol {U}(x,y)=(U,\,V,\,0)^T$ and $\boldsymbol {\Omega }(x,y)=(0,\,0,\,\Omega _z)^T$ correspond to the the steady baseflow, and $\boldsymbol {u}^\prime (\boldsymbol {x},t)$ and $\boldsymbol {\omega }^\prime (\boldsymbol {x},t)$ denote the disturbance flow. This framework allows us to capture temporally decaying/growing global modes in the context of convective/absolute instability mechanisms. Note that this approach fully takes into account all the non-parallel effects with respect to the baseflow and the form of the disturbances.

Since the base flow is only a function of $x$ and $y$ and homogeneous in the third direction, 3-D (oblique) disturbance waves are further decomposed in the $z$ -direction as

(3.4a) \begin{align} \left [u^\prime ,v^\prime , \omega _z^\prime \right ]^{T}(\boldsymbol {x},t)&=\left [ \hat {u}, \hat {v}, \hat {\omega }_{z}\right ]^T (x,y,t)\cos (\gamma z)\,, \end{align}

(3.4b) \begin{align} \left [w^\prime ,\omega _x^\prime , \omega _y^\prime \right ]^{T}(\boldsymbol {x},t)&=\left [ \hat {w}, \hat {\omega }_{x}, \hat {\omega }_{y}\right ]^T (x,y,t)\sin (\gamma z)\,, \end{align}

where $\gamma$ is the spanwise wavenumber. In the LGSA, disturbances with a specified $\gamma$ were perturbed through a narrow blowing and suction slot on the flat-plate surface by pulsing the wall-normal velocity of the following form:

(3.5) \begin{equation} \hat {v}_k(|x-x_f|\leqslant b/2,y=0,t)= S(x) \exp \left ( \left [\frac {t-t_p}{\delta t_p}\right ]^2 \right ), \end{equation}

where $x_f$ and $b$ correspond to the centre and streamwise width of the slot, respectively, and are set to $x_f=19$ and $b=0.2$ . The shape function, $S(x)$ , is a polynomial which is zero outside the forcing slot such that smooth derivatives at boundaries of the suction/blowing slot are obtained (Hosseinverdi & Fasel 2019). In (3.5), the temporal function is realized by a Gaussian distribution, where $t_p=10 \triangle t$ and $\delta t_p=0.25\triangle t$ , for which the corresponding frequency spectrum is broadband with approximately the same amplitude levels for a frequency range up to $St = f L_\infty ^\ast / U_\infty = 100$ .

For a specified $\gamma$ , the temporal growth rate of the least stable mode (or the most unstable one in case of an absolute instability) will determine the long-time behaviour of the flow with respect to 2-D/3-D disturbances. Figure 11(a) shows the temporal evolution of the maximum disturbance kinetic energy, $\hat {q}_{k}$ , for a few selected values of $\gamma$ for the mean flow obtained from the DNS with zero FST, where $\hat {q}_{k}$ is computed as $\hat {q}^2_k= (\hat {u}^2 + \hat {v}^2 + \hat {w}^2 )$ . The linear response of the mean flow to the wave packet disturbances can be classified into two regions: (i) a short-time behaviour for $0\lt t\lt 30$ , which is associated with the convective instability characteristics of the base flow (cf. LST) and (ii) a long-time behaviour which indicates the global stability of the mean bubble, i.e. being globally stable or unstable.

Figure 11. (a) Time evolution of maximum disturbance $\hat {q}_k$ in the logarithmic scale for DNS without FST. (b) Temporal growth rate as a function of spanwise wavenumber.

In figure 11(a), the decay/growth of the maximum $\hat {q}_k$ for 2-D/3-D modes indicate that the mean LSB obtained from the DNS without FST is globally stable/unstable with respect to the 2-D/3-D disturbances. Analysis for 2-D and 3-D disturbances in the range of $0\leqslant \gamma \leqslant 10\pi$ for the base flow from the DNS and the temporal growth rate associated with each $\gamma$ is then computed based on $\tilde {\sigma }=\log (\hat {q}_{k,max})/\triangle t$ . The distribution of the temporal growth rate as a function of $\gamma$ is displayed in figure 11(b). It can be observed that the LSB is asymptotically stable with regard to $\gamma =0$ (2-D mode) and becomes globally unstable for a range of $\gamma$ , i.e. $\pi \leqslant \gamma \leqslant 9\pi$ , where the maximum $\tilde {\sigma }$ is attained for $\gamma =3\pi$ . The presence of this global instability explains the self-sustaining transition to turbulence in the LSB, consistent with findings by Rodríguez et al. (2021).

Subsequently, using the time-dependent data available from the LNSE, dynamic mode decomposition (DMD) was employed to extract the relevant unstable modes. When a discrete dataset of the flow field is taken at $N$ snapshots separated by a constant sampling time interval $\triangle t$ , DMD allows for an expansion of the velocity data in the form

(3.6) \begin{equation} \boldsymbol {u}(\boldsymbol {x},t)=\sum _{j=1}^{N} \boldsymbol {\psi }_j(\boldsymbol {x}) \exp {\left ([\sigma _j +i\omega _j] t\right )} , \end{equation}

where $\boldsymbol {\psi }_j$ are the spatial DMD modes and $\sigma _j$ and $\omega _j$ represent the growth rates and oscillatory radial frequencies, respectively.

For example in figure 12(a) the DMD eigenvalue spectrum in theSt- $\sigma$ plane is displayed for $\gamma =3\pi$ (i.e. most unstable $\gamma$ found, see figure 11(b) based on the mean flow from the DNS without FST. The DMD eigenvalues were found to be in complex conjugate pairs, therefore, the spectrum is symmetric about $St=0$ . Note that the data in the asymptotic state ( $45\lt t\lt 120$ ) were used for the DMD. The first three most unstable modes are highlighted as solid (red) symbols. It is worth noting that the largest growth rate obtained from the DMD ( $\sigma \approx 0.2$ ) is very close to $\tilde {\sigma }$ computed based on $\hat {q}_{k,max}$ for $\gamma =3\pi$ (see figure 11 b). The spatial eigenfunctions of the most unstable modes are plotted in figure 12(b), which indicate the trapped modes inside the separated region of the mean bubble.

Figure 12. Dynamic mode decomposition analysis of the LNSE results for $\gamma =3\pi$ based on the mean flow obtained from the DNS with zero FST: (a) eigenvalue spectrum and (b) eigenfunctions of the most unstable modes. Shown are contour lines of $\Re (v)$ -velocity.

A similar analysis was performed for 2-D disturbances ( $\gamma =0$ ) using the mean flow obtained from the DNS with zero FST. The resulting DMD eigenvalue spectrum is displayed in figure 13. Consistent with the observation in figure 11(b), for this case the mean separation bubble is asymptotically stable ( $\sigma \lt 0$ ). Of special interest is the frequency associated with the least stable mode ( $St=0.67-0.68$ ), which is highlighted in figure 13. The relevance of this will be re-adressed later in § 4.

Figure 13. Dynamic mode decomposition eigenvalue spectrum based on the LNSE data for $\gamma =0$ . Mean flow obtained from the DNS without FST. (b) The DMD eigenmode corresponding to $St=0.67$ as obtained based on LNSE for $\gamma =0$ ; results based on the mean flow field from the DNS without FST.

3.4. Direct numerical simulations with free-stream turbulence

For meaningful comparisons of numerical simulations and wind-tunnel experiments the influence of FST on the LSB has to be addressed. In the DNS results presented and discussed in this paper, low levels of FST were introduced at the inflow boundary by superimposing the velocity and vorticity components of the steady base flow with velocity and vorticity fluctuations. The methodology adopted here to generate free-stream disturbances at the inflow boundary is based on a superposition of eigenmodes from the continuous spectra of the OrrSommerfeld (OS) and homogeneous Squire (SQ) operators (Grosch & Salwen 1978; Jacobs & Durbin 2001; Brandt, Schlatter & Henningson 2004; Hosseinverdi et al. 2012; Hosseinverdi & Fasel 2018, 2019)

(3.7) \begin{equation} \boldsymbol {u}^\prime (x_i, y, z, t)=\sum _\omega \sum _{\gamma } A(|\boldsymbol {k}|)\,\boldsymbol {\phi }(y; \omega , \beta , \gamma )\,e^{i(\gamma z-\omega t)}\,. \end{equation}

Here, $\omega$ is the angular disturbance frequency, $\beta$ and $\gamma$ are the wall-normal and spanwise wavenumbers, respectively, with $|\boldsymbol {k}|= (\omega ^2+\gamma ^2+\beta ^2 )^{0.5}$ . The dispersion relation $\alpha =\omega /U_\infty$ was used in (3.7) to express the streamwise wavenumber ( $\alpha$ ) in terms of the angular frequency. The inflow disturbance vorticity field is calculated from the disturbance velocity field as $\boldsymbol {\omega }^\prime =-\nabla \times \boldsymbol {u}^\prime$ and the eigenfunction $\boldsymbol {\phi }$ is a normalized, randomized and weighted superposition of OS and SQ eigenmodes. The normalization ensures that the energy of each disturbance mode is unity. The coefficients $A(|\boldsymbol {k}|)$ are used to determine the contribution of the eigenfunctions to the total turbulent kinetic energy. For the simulations presented in this work, the amplitudes of the individual inlet perturbations are assigned such that a von Kármán energy spectrum is obtained. A detailed description of the implementation and validation results are provided in Hosseinverdi et al. (2012) and Hosseinverdi & Fasel (2018, 2019). Several levels of FST have been considered, where the best agreement with the wind-tunnel experiments was obtained for $Tu=0.02\, \%$ , which is comparable to the measured levels of FST in the experiments (§ 2.2). Addition of small levels of FST in the DNS did not change the velocity profiles upstream of the LSB (figure 8 a), thus maintaining the good agreement for the time-averaged flow field upstream of the LSB.

When isotropic FST with an intensity of $Tu=0.02\, \%$ is introduced at the inflow boundary, owing to the large growth rates associated with the inviscid instability mechanism, the vortical disturbances due to the FST are rapidly amplified by several orders of magnitude (Balzer & Fasel 2016; Hosseinverdi & Fasel 2019). This leads to an upstream movement of the transition onset location, which in turn reduces the mean bubble size (both in the $x$ - and $y$ -directions) compared with the zero FST case. Analogous to figure 8(b), instantaneous contours of spanwise vorticity (averaged in the $z$ -direction) from the 3-D DNS with FST are shown in figure 14 together with the time- and spanwise-averaged dividing streamline.

Figure 14. Instantaneous flow visualizations for $Tu=0.02\, \%$ . Shown are contours of spanwise vorticity (averaged in spanwise direction) together with the mean dividing streamline identified by black dashed lines.

As for the case without FST (see figure 11), the LSB obtained from the DNS with FST was found to be globally stable/unstable with respect to the 2-D/3-D disturbances, respectively. Both LSBs (with and without FST) are asymptotically stable with regard to the 2-D mode and become globally unstable for a bounded range of spanwise wavenumbers ( $\gamma$ ), with an increase in maximum temporal growth rates in the DNS with FST compared with the case without FST. Dynamic mode decomposition results for the case with FST were found to be very similar to those presented in § 3.3 for the case without FST and therefore not repeated here. The DMD eigenvalues associated with the least stable mode ( $St=0.67{-}0.68$ ) agree well with those for the case without FST (see figure 13) which confirms that the least stable DMD mode in the presence of FST still corresponds to the dominant shear-layer mode captured in the DNS.

The structure, location and frequency of the most dominant mode in the stability calculations and DNS provide strong evidence of an inviscid shear-layer instability mechanism based on the inflectional velocity profiles. To confirm this conjecture, LST analysis based on the mean flow obtained from the DNS with FST was performed. The downstream evolution of a 2-D $u^{\prime }$ -velocity disturbance wave with $St=0.67$ obtained from the DNS with FST is compared with the LST results in figure 15(a). Based on the stability diagram in figure 10 the location corresponding to the critical Reynolds number for $St=0.67$ is located between the onset of the APG and the separation location ( $x_{cr}\approx 20.7$ ). While the onset of the exponential amplification of the dominant 2-D mode in the DNS is observed further downstream, results for $x\geqslant 24$ inside the separated region are in excellent agreement with the LST.

In figure 15(b) the wall-normal Fourier amplitude distributions of the $u^{\prime }$ -velocity between DNS and LST are compared at several streamwise locations within the region of exponential growth. The amplitude distributions from both LST and DNS are normalized by their respective maximum values within the boundary layer. The amplitude distributions from the DNS and LST agree remarkably well within this region, supporting the fact that the same flow features are captured. The maximum of the $u^{\prime }$ -velocity occurs close to the displacement thickness ( $y/\delta _1=1$ ) of the local boundary-layer profile, which in turn is very close to the local inflection point in the velocity profile. This further supports that the exponential growth of this most dominant mode is due to the inviscid shear-layer instability.

Figure 15. Comparison of disturbance evolution and eigenfunction distributions between LST and DNS with FST intensity of $Tu=0.02\, \%$ . (a) Downstream development of the Fourier amplitude (maximum inside the boundary layer) for the $u^{\prime }$ -velocity of the 2-D disturbances with $St=0.67$ ( $f={185}\,\mathrm {Hz}$ ). (b) Wall-normal distribution of $u^{\prime }$ -velocity Fourier amplitude taken at the several downstream locations. Disturbance amplitudes are scaled by their respective maximum value and $y$ -coordinate is normalized with the local displacement thickness.

4. Comparison between experiments and simulations

Results for the LSBs from the experiment and numerical simulations are compared for the time-averaged flow field and the dynamics in the separated shear layer. In figure 16 the surface pressure development along the flat plate from the experiment is shown next to numerical results from the ‘precursor calculation’ and the DNS with and without FST. Surface pressure for all cases is in excellent agreement upstream of $x = 25$ . The RANS model, known to struggle with separated flows, underpredicts the size of the LSB. Without FST the DNS result deviate from the wind-tunnel experiments for $x \geqslant 25.5$ , predicting a significantly longer separated region (see § 3.2) As expected, the curve for the DNS with Tu = 0.02 % (§ 3.4), indicates a smaller size for the mean bubble than the case for zero FST. With the addition of FST, the pressure development shows excellent agreement between the DNS and the experimental results (§ 2.4).

Figure 16. Comparison of flat-plate surface pressure coefficient, $C_p$ , between experiment, DNS without and with FST as well as RANS calculation.

More detailed comparisons of the mean flow via streamlines is presented in figure 17. Streamlines of the time-averaged LSB obtained in the DNS with FST compare remarkably well with the experiment, showing a near identical bubble topology in both cases. This similarity is further corroborated by the integral boundary-layer development in figure 18. The streamwise development of the displacement and the momentum thickness Reynolds number, $Re_{\delta _1}$ and $Re_\theta$ , is compared for the experiment and the DNS with FST.

Figure 17. Comparison of time-averaged streamlines for the experiments and DNS with $Tu = 0.02\, \%$ .

Figure 18. Comparison of time-averaged displacement and momentum thickness Reynolds number, derived from semi-empiric similarity solver upstream of separation (black), extracted from PIV in the experiment (blue, symbols) and DNS with $Tu = 0.02\, \%$ (red). Shaded area corresponds to mean separated region.

The boundary-layer parameters from the experiment were obtained from $u$ -velocity profiles extracted from PIV, enforcing the no-slip condition at the wall. Upstream of the PIV measurements, a semi-empirical solver suggested by Cebeci & Bradshaw (2012) is used to calculate the boundary layer based on the wind-tunnel conditions and measurements of the surface pressure of the flat plate in the experiment (figure 5). Downstream of $x = 18.5$ and upstream of the time-averaged separation location, boundary-layer data extracted from PIV and predicted values are in good agreement. Furthermore, predictions based on surface pressure and results from PIV agree well with the DNS (with FST) as can be observed in figure 18. Downstream of $x=21$ , the rapid growth of the $Re_{\delta _1}$ corresponds to the onset of the boundary-layer separation. In contrast, $Re_\theta$ varies slowly downstream of separation and across the separated flow region of the bubble. Both parameters, $Re_{\delta _1}$ and $Re_{\theta }$ , show excellent agreement between experiment and DNS with FST in the entire mean separated region. The Reynolds numbers based on displacement and momentum thickness at separation are $Re_{\delta _1,s}=909$ and $Re_{\theta ,s}=270$ in the DNS and $Re_{\delta _1,s}=919$ and $Re_{\theta ,s}=314$ in the experiment. The small discrepancy in $Re_{\theta ,s}$ is attributed to the small boundary-layer thickness at separation, which is challenging to resolve quantitatively with PIV.

The velocity profiles for the streamwise velocity component ( $u$ ) for various downstream locations are provided in figure 19. Upstream of $x=25.5$ , the velocity profiles essentially collapse for the experiments and the DNS with and without FST. The time-averaged velocity profiles upstream of the LSB were unaffected when low levels of FST were added. Earlier onset of transition in the experiment and DNS with FST lead to significant deviation of the velocity profiles for $x \geqslant 25.5$ (and streamlines, figure 9 a). All time-averaged quantities, streamlines (figure 17), boundary-layer parameters (figure 18) and velocity profiles (figure 19) from experiments are in close agreement with the DNS when FST is added, at integral levels comparable to the FST measured in the experiment (§ 2.2).

Figure 19. Wall-normal distributions of $u$ -velocity of mean results between experiments and DNS at several $x$ locations. Here, $u$ is normalized with the local boundary-layer edge velocity.

Further relevant comparison between experiment and DNS can be made by scrutinizing the dynamics in the transitional flow in the LSB. Large-scale energetic structures can be identified using proper orthogonal decomposition (POD) Lumley (1967). Towards this end, the snapshot method introduced by Sirovich (1987) was implemented following Weiss (2019). Decomposition of the time ( $\textit {t}$ ) and space ( $\boldsymbol {x}$ ) dependent flow field yields

(4.1) \begin{align} u(\boldsymbol {x},t_n) = \sum ^{N-1}_{i=0} a_i(t_n) \Phi _i(\boldsymbol {x}), \end{align}

where $\Phi _i$ represents the spatial POD eigenfunctions and $a_i$ the POD time coefficients.

The resulting eigenfunctions (modes) are orthogonal and sorted by their contribution to the total kinetic energy of the flow field: the first mode representing the largest fraction of the total kinetic energy. Small mode numbers in the POD will show the most dominant structures in the flow field, assuming the existing coherent structures in the flow are of high kinetic energy. For the POD analysis presented and discussed in this paper, 1000 snapshots of the $v$ -velocity components from the DNS, and 1500 vector fields from the PIV measurements were used. As previously discussed in § 2.4, analysis of the $v$ -velocity component in the x-y plane in the experiment highlights the 2-D modes that develop in the shear layer (Lengani et al. 2014; Simoni et al. 2017). It should be noted that while the POD analysis was carried out for the 2-D field in the experiments, the POD modes from the DNS are based on the spanwise-averaged results.

The POD energy spectrum is shown in figure 20. The left vertical axis represents the fraction of energy, $E_i$ , and the right vertical axis indicates the cumulative energy sum, $S_i$ , defined as

(4.2) \begin{equation} E_i=\left(\lambda _i/\sum _{j=1}^{N-1} \lambda _j\right)\times 100, \,\, S_i=\left(\sum _{k=1}^i \lambda _k/\sum _{j=1}^{N-1} \lambda _j\right)\times 100. \end{equation}

Figure 20. The POD energy spectrum obtained from experiments and DNS with $Tu=0.02\,\%$ . Plotted are fraction of energy (left vertical axis) and cumulative energy sum (right vertical axis).

In the above equation, $\lambda _i$ is the magnitude of the eigenvalue of the POD mode, which corresponds to twice the kinetic energy content of the respective POD mode, and $N$ is the total number of POD modes. The POD eigenvalue spectrum in figure 20 reveals that most of the energy is contained in the lower mode numbers for both the experiments and DNS, indicating the existence of significant coherent structures in the flow. In particular, there are several mode pairs with almost identical magnitudes. These modes are linked to downstream travelling structures. Each pair representing the same structure with a phase shift of approximately $1/4$ of the respective wavelength. Several mode pairs are also observed in the experiments (1 and 2, 3 and 4, 5 and 6) with only mode 1 and 2 apparent as a pair in the energy spectrum of the DNS. Differences between the amplitude levels of the dominant POD modes and number of mode pairs between the experiments and DNS are likely a result of the 3-D field used in the POD analysis of the DNS, whereas only PIV measurements in the 2-D x-y plane were available in the experiment.

Comparison of the eigenfunctions of the most energetic POD modes obtained from the experiments and DNS with $Tu=0.02\, \%$ are presented in figure 21. For better visual comparison, for each mode pair, the mode is chosen such that the phase difference between the observed structures is minimal. For reference, the dividing streamline is also shown in to indicate figure 21 the bubble size and shape. In figure 21, the entire mode pair 1 and 2 is presented and compared for the experiment and DNS with FST in order to also allow observation of the phase shift and the convective nature of the structures. For brevity, only one mode is shown for each of the subsequent mode pairs. Modes 1 and 2 originate well upstream of the saturation of the primary instability (figure 15) and indicate strong coherent structures in the separated shear layer. This mode is directly related to the shear-layer mode (K–H instability), showing remarkable agreement between the experiments and DNS in the LSB region. Mode 9 (experiment) and mode 8 (DNS) capture similar size structures of weaker energy downstream of reattachment. Modes 3 and 5 in the experiment as well as modes 4 and 5 in the DNS capture larger structures downstream of reattachment. The onset of these larger-scale structures coincides with the maximum of fluctuating quantities and develops downstream of estimated transition locations around reattachment.

Figure 21. Dominant POD eigenfunctions from experiments (left) and DNS with $Tu=0.02\, \%$ (right). Plotted are contours of $v$ -velocity component.

In figure 22, stereoscopic PIV measurements are presented for a plane parallel to the wall at a constant wall distance (close to the maximum bubble height, $y = 0.35$ ), next to corresponding results from the DNS with FST. The first four POD modes based on the $v$ -velocity component show the spanwise extent of the dominant flow structures. The first two modes capture the shedding vortex structures in the LSB shear layer. The phase shift between mode 1 and 2 confirms the convective nature of the observed structures. Modes 3 and 4 show spanwise variations of the flow, resembling a checkerboard-like pattern. Modes 3 and 4 also show a phase shift indicating the convective nature of these structures in the streamwise direction with a spanwise wavelength of approximately half of the domain width ( $\triangle z = 1$ ). Around mean reattachment, several larger-scale structures are apparent which may be related to those found in the wake of the LSB in modes 3 and 5 in the x-y plane (figure 21). Comparison of the $v$ -component to respective POD modes extracted from DNS simulations with FST ( $Tu = 0.02\, \%$ ) show good agreement in the spanwise coherent shear-layer shedding in modes 1 and 2. Modes 3 and 4 show downstream travelling structures, similar in streamwise and spanwise extent, resembling checkerboard-like patterns downstream of maximum bubble height (spanwise shifted with respect to the experiment). Spanwise modulation of the 2-D rollers in the POD, resembles the peak valley formation and spanwise vortex splitting observed in Michelis et al. (2018).

Figure 22. Dominant POD eigenfunctions from experiments (left) and DNS with $Tu=0.02\,\%$ (right). The location of the wall-parallel PIV plane relative to the LSB streamlines in the x-y plane is shown as the dash-dot line in the top figures, at $y = 0.35$ . Contours of dominant $v$ -eigenfunctions derived from the $v$ -velocity field are shown in (a) and (b) and contours of dominant $u$ -eigenfunctions based on the $u$ -velocity in (c) and (d), respectively.

The first two POD modes derived from the $u$ -velocity component are shown in figure 22. Mode 1 is a result of the measurement planes proximity to the maximum gradient in the $u$ -velocity profile at maximum bubble height. Small low frequency temporal variations in bubble extent are dominant in this plane. More interesting is the observation of streamwise streaks in mode 2 in the experiment and DNS (and mode 1 at and downstream of reattachment). The onset of these structures inside the bubble and downstream of maximum bubble height coincide with the aforementioned spanwise modulation of the spanwise coherent structures in the early stages of transition. The wavelength of the observed structures is in good agreement with the wavenumber associated with the the largest temporal growth rates figure 11, possibly showing the presence of such a global instability in the experiment, however, additional work is needed to corroborate this.

Figure 23. Power spectral density of the $u$ -velocity at various streamwise locations inside the separated region extracted along $y=0.55$ , obtained in (a) experiments and (b) DNS with $Tu=0.02\,\%$ . Included is the result obtained from the LST at $x=26$ . Vertical black dashed-dotted lines correspond to the dominant shedding frequency. Additional light blue vertical dashed lines in (b) indicate relevant frequencies based on the distinct wavenumbers chosen for free-stream generation in the DNS.

Fluctuating quantities in the shear layer near the maximum bubble height were measured using CTA in the experiment. Results presented in figure 23(a) show a broad peak centred at 250 Hz ( $St = f L_\infty ^\ast / U_\infty = 0.88$ ). The probe is placed in the high velocity region of the shear layer, $y=0.55$ , above the inflection point ( $y\approx 0.37$ ) to avoid obstruction of the measurements due to reverse flow. Along the separated shear layer the dominant frequency represents the natural shedding of the characteristic spanwise vortical ‘rollers’. Following a growth in amplitude towards the maximum bubble height, subsequent turbulent breakdown leads to approximately turbulent spectra at the location of maximum fluctuations (see figure 6), indicated by $f^{-5/3}$ slope in figure 23(a) (grey dashed line). Corresponding results from LST capture a similar broad peak of most unstable frequencies with the maximum located at a slightly lower frequency. Spectra of the $u$ -velocity from the DNS with $Tu=0.02\, \%$ were extracted at the same locations. The PSD in figure 23(b) exhibits multiple distinct peaks with its maximum at 185 Hz ( $St=0.671$ ) slightly below the LST peak but well within the envelope. Additional secondary peaks were found at higher frequencies filling up the broad band of amplified frequencies predicted by LST. The dominant frequency is very close to the most unstable mode observed in the DMD (see figure 13). The corresponding DMD eigenfunction, in figure 13(b), associated with $St=0.67$ as obtained from the LNSE and DNS without FST is in good agreement with the dominant POD modes from experiment and DNS with FST (figure 21) and confirms that the least stable DMD mode indeed corresponds to the dominant shear-layer instability captured in the DNS and experiment.

The presence of distinct peaks in the DNS in contrast to the continuous spectra in the experiment is a result of the generation of FST based on a finite number of eigenmodes of the OS and SQ operators (dashed vertical lines in figure 23 b). Harmonic content of the most dominant modes is found in both experiment and DNS. In addition, at locations downstream of maximum bubble height, subtle subharmonic content is present in the PSD derived from the experiment. Differences in the dominant frequency between the experiment and DNS are likely a result of mismatches in the FST spectra. The addition of isotropic turbulence in the DNS (§ 3.4), results in an ideal von Kármán energy spectrum along a finite number of modes. The spectra of the free-stream measurements in the experiments are, unsurprisingly, different from the theoretical spectrum used in the DNS (figure 3). The maximum in the broad band of most unstable frequencies estimated by LST is $St = 0.78$ , which is closer to the maximum found in the experiment. An increase in frequency resolution of the FST in the DNS would be expected to fill in the broad peak predicted by LST, similar to observations in Grille Guerra et al. (2022), where broad band random disturbance input in the boundary layer upstream of separation (bypassing the receptivity process) resulted in the same broad band frequency peak as in the corresponding experiments. The LST analysis predicts that in the base flow of the LSB in the DNS with FST, disturbances in the range of $0.5\leqslant St\leqslant 1$ exhibit very similar stability characteristics (figure 10). Due to the similarity between the DNS with FST and the experiments in the time-averaged flow the same characteristics can be expected in the experiment. Therefore, minute differences in the disturbance characteristics in the free stream and subsequent exponential amplification in the shear layer can lead to the differences in the observed frequency between the experiment and DNS.

5. Conclusions

Stability and transition in the presence of a LSB was investigated in a combined approach using experiments DNS and local and global stability analysis. A LSB on a flat plate was generated by an inverted wing (with a NACA $64_3-618$ airfoil) placed in proximity to a flat plate. In contrast to recent studies using contoured wind-tunnel ceilings (Diwan & Ramesh 2009; Michelis et al. 2018; Dellacasagrande et al. 2020), the inverted wing in the present work was motivated by Gaster (1967). However, contrary to the jet entrainment scheme at the trailing edge of the airfoil used by Gaster (1967), boundary-layer suction was employed on the suction side of the wing. Despite using a different airfoil and chord length, the same APG along the flat plate as in Gaster (1967), Series I, Case 6, for similar free-stream velocity could be obtained. Streamlines of the time-averaged flow field around the inverted wing in the experiment agree with the steady RANS simulations, providing boundary conditions for high-fidelity DNS.

Comparison of DNS results with and without FST confirm the influence of FST to the topology of the time-averaged LSB and in particular on the transition physics. Low levels of FST ( $Tu = 0.02\, \%$ ) in the DNS, comparable to the measured FST levels in the wind tunnel (see § 2.2), lead to earlier onset of transition and resulted in a time-averaged flow field that matched the experiments remarkably well. Comparison of the streamwise velocity fluctuations in the separated shear layer reveals a difference in the unsteady frequency content ( $St = 0.88$ , ${250}\,\mathrm {Hz}$ in experiments vs $St = 0.671$ , ${185}\,\mathrm {Hz}$ in DNS with FST). The LST analysis of the time-averaged flow field from the DNS with FST also indicated that both frequencies are within the range of the most unstable disturbances. The LST predicted similar amplification rates for a broad spectrum of most unstable frequencies between $0.5 \lt St \lt 1$ ( ${150}\,\mathrm {Hz}$ and ${300}\,\mathrm {Hz}$ ). The observed differences in the dominant frequency are likely a result of the modelling of the FST in the DNS. Despite the difference in frequency, POD analysis identified similar dominant coherent structures in DNS with FST and experiment. The most energetic downstream travelling structures were found in the shear layer in the region of exponential disturbance growth, confirming the relevance of the inviscid shear-layer instability. This was further corroborated by LST analysis which showed excellent agreement with the DNS with FST. Global stability analysis shed light on the K–H-type instability and explained the self-sustained transition mechanism for the LSB without FST. While both LSBs from DNS (with and without FST) are asymptotically stable, with respect to the 2-D mode, they become globally unstable for a range of 3-D wavenumbers. Analysis of the 2-D modes show the least stable mode coincides with the dominant shedding frequency in the DNS, and that the respective DMD mode compared well with the respective POD modes and the K–H mode in the shear layer of the LSB in the DNS and experiment.

The experimental and numerical investigation in this study were based on the same geometry and free-stream conditions. However, it should be noted that the measurements and numerical calculation results of the respective flow fields are entirely independent from one other, i.e. the wind-tunnel experiments did not require input from the numerical calculations and vice versa. Therefore the very good agreement between experiment, simulations and theory is indeed remarkable for both the time-averaged flow field and the time-dependent dynamical characteristics of the LSB. The detailed understanding of the transition physics gained from this investigation can serve as a basis for future studies of active flow control.