Directed percolation in aerodynamics: resolving laminar separation bubble on airfoils

In nature, phase transitions prevail amongst inherently different systems, while frequently showing a universal behavior at their critical point. As a fundamental phenomenon of fluid mechanics, recent studies suggested laminar-turbulent transition belonging to the universality class of directed percolation. Beyond, no indication was yet found that directed percolation is encountered in technical relevant fluid mechanics. Here, we present first evidence that the onset of a laminar separation bubble on an airfoil can be well characterized employing the directed percolation model on high fidelity particle image velocimetry data. In an extensive analysis, we show that the obtained critical exponents are robust against parameter fluctuations, namely threshold of turbulence intensity that distinguishes between ambient flow and laminar separation bubble. Our findings indicate a comprehensive significance of percolation models in fluid mechanics beyond fundamental flow phenomena, in particular, it enables the precise determination of the transition point of the laminar separation bubble. This opens a broad variety of new fields of application, ranging from experimental airfoil aerodynamics to computational fluid dynamics.


I. INTRODUCTION
Three decades ago, Pomeau described the flow of a fluid as a collection of oscillators that interact with each other. When observing that [1] "each oscillator if in a turbulent state may either relax spontaneously toward its quiescent state or contaminate its neighbors" he concluded that "this is precisely the definition of the process called directed percolation in statistical physics" and therefore raised the possibility of laminarturbulent transition belonging to the same universality class as directed percolation (DP).
The idea is quite reasonable in the sense that transition from laminar to turbulent flow (hereafter referred to as transition) can be described via the so-called spatio-temporal intermittency [2,3]. While since then several simulations have supported Pomeau's conjecture, only in the last few years it was possible to provide experimental evidence [4], due to novel possibilities of extracting accurate measurements from a turbulent flow with sufficient spatial and temporal resolution. Recent studies approached transition from different angles by means of low-order models [5][6][7], sophisticated simulations [8][9][10], as well as experiments [11][12][13]. They concordantly indicate non-equilibrium phase transition occurring in basic shear flows, i.e. pipe, channel and Couette flows. While investigating the hypothesis of transition belonging to the universality class of directed percolation, they treat transition into turbulence as a fundamental phenomenon. Up to now, directed percolation has not been used in flows with a direct aim of helping to solve specific needs in applied sciences that deal with turbulence by extension.
Here, we argue that the applicability of directed percolation model can be extended from fundamental fluid dynamics to practical aerodynamics relevant for engineering problems and, thus, to a more generally valid concept. More specifically, we provide evidence that directed percolation is capable of characterizing the onset of a laminar separation bubble Photograph of an LSB on airfoil CK220. Laser light sheet illuminates the particle seeded inflow from right-hand side and perpendicular to the airfoil's surface resulting in a shadow at the leading edge. Inflow is coming from left-hand side. LSB can be identified as a region without smoke between the airfoil's surface and the ambient flow. False colors are used for better visualization. (b) Labeled zoom, emphasizing the flow topology at the LSB. LSB occurs just behind the airfoil's thickest cross section, where streamlines (in blue) separate from the surface. Further downstream, when LSB is at its maximum expansion, the laminar shear layer destabilizes resulting in transition to turbulence. High momentum perpendicular to the airfoil's surface contained in the shear layer enables reattachment and the formation of a turbulent boundary layer.
(LSB) on the suction side of an airfoil (see Fig. 1) [14,15]. This is of particular importance since the aerodynamic functionality of an airfoil depends sensitively on such transitions like the LSB. Furthermore, the precise determination of the onset of transition into LSB as well as transition in general are still an open problem in aerodynamics. So far, there are no methods available that reliably capture the nonlinear nature of an LSB. In this way, DP is potentially valid for nonlinear flow phenomena encountered in the class of external flows over curved surfaces, where the flow becomes more and more unstable as it is advected downstream. Based on our findings, we discuss how DP properties can also be of use in applications, ranging from computational fluid dynamics (CFD), through the design process of airfoils, to maintenance of rotor blades in wind turbines.
The paper is structured as follows. At first, the experimental setup is described in Sec. II. This includes an airfoil subjected to laminar inflow and the optical flow visualization method. In Section III, a procedure is explained how to map the measured flow field into binary cells, either laminar or turbulent, which enable for estimating a critical Reynolds number identifying the onset of LSB. Characteristic exponents at this critical value are determined in Sec. IV and discussed subsequently in Sec. V. Finally, Section VI concludes the paper.

II. THE EXPERIMENT
The phenomenon of LSB can be qualitatively visualized by means of a laser illuminated smoke photography (see Fig. 1). The seeded flow is coming from left-hand side, separates just after the airfoil's thickest cross section, destabilizes and undergoes transition resulting in reattachment of the turbulent boundary layer. The recirculation zone between the airfoil's surface and the shear layer is the one usually addressed as LSB. In Figures 1a and 1b, the LSB appears as the smoke free region above the airfoil. While transition, taking place within the shear layer between LSB and free flow, is threedimensional, the LSB's onset happens in a linearly stable laminar flow region whose boundary layer has a thickness small compared to the dimensions of the LSB. Therefore, transition into LSB can be approximated locally as a quasi twodimensional process.
The formation of an LSB is a highly unsteady process which emerges stochastically over a certain region on airfoils. In order to properly investigate this complex and delicate flow topology, it is crucial to use non intrusive methods featuring high spatio-temporal resolution. Therefore, stereoscopic high-speed particle image velocimetry (HSPIV) is used to visualize an LSB on the suction side of a CK220 airfoil in the wind tunnel. A schematic representation of the experiment is given in Fig. 2. The experimental setup consists of two Phantom Miro M320S high-speed cameras and a Litron LDY303 laser. This enables HSPIV measurements with a recording frequency of 2,000 velocity fields per second and a recording length of T ≈ 3 s at reduced resolution of 896 × 792 px 2 . For a measurement region of ∆x × ∆y = 40 × 40 mm 2 , where x and y are in chordwise, respectively spanwise, direction, a spatial resolution of dx = dy < 0.4 mm is obtained with a sufficient accuracy (stereo residue below 0.5 px). This resolution corresponds to dx/c < 2 × 10 −3 at given airfoil dimension of c × s = 220 × 250 mm 2 , where c denotes the airfoil's chord length and s is its span.
As illustrated in Fig. 2, the light sheet is adjusted tangen- tially to the airfoil's surface approximately at the onset of the LSB at a distance dz/c < 3 × 10 −3 or dz = 0.5 mm, respectively. Due to the airfoil's curvature, this distance varies slightly, by less than 5% of the LSB's thickness. The onset of LSB is thus captured properly. At the same time, the measurement region covers the LSB completely.
In the present experiment, the LSB appears at a global Reynolds number of Re chord = 160, 000 and an angle of attack of α = 5 • . Re chord is obtained with respect to the chord length and a free stream velocity of u ∞ = 11 m/s. In this way, the ambient conditions of the experiment are described in general. Moreover, any transition process is strongly dependent on detailed flow conditions, such as free stream properties. To minimize this impact, the experiments are performed in a wind tunnel with inlet 250 × 250 mm 2 , length 2000 mm and a closed test section with a low turbulence intensity of TI < 0.004.

III. FROM PARTICLE IMAGE VELOCIMETRY DATA TO LAMINAR AND TURBULENT STATES
Following the objective of this study, the LSB's onset needs to be derived from experimental data. HSPIV results provide instantaneous velocity fields of the whole LSB and the flow upstream as well as downstream of the LSB. Characteristic differences between those three regions are reflected by velocity magnitude u and the fluctuations u . Both quantities define the turbulence intensity TI. The value of TI is usually FIG. 3. Illustration of measured and computed quantities: The plane of space and control parameter (Rex) is formed by the airfoil's spanwise and chordwise direction, where Rex = Re(x, u∞) is based on the distance x from the leading edge and the inflow velocity u∞. Rex provides a measure for the flow's likelihood to undergo transition from laminar flow into an LSB and, therefore, serves as a control parameter. The presence of LSB is indicated by means of turbulence intensity (TI) exceeding a certain threshold denoted as turbulent phase. The instantaneous turbulence intensity is derived from 10 temporally consecutive velocity fields at one respective position Rex in a plane crossing the laminar separation bubble. Phase transition from laminar flow into LSB happens once turbulent clusters merge significantly. At this critical point, Rec, characteristic turbulent clusters occur, depicted as an evolution in time (vertical plane). Averaging over time at each position Rex, the fraction of turbulent cells ρ is shown next to the plane of TI. Following ρ in streamwise direction, phase transition from laminar boundary layer flow into the LSB and subsequent reattachment are obvious. The first phase transition represents the LSB's onset which is important for our work. derived locally from the whole data set losing its temporal information, whereas directed percolation analysis is based on the time-resolved evolution of laminar and turbulent regions. Therefore, we define the TI at the time t as the TI computed locally at the position (x, y) from 10 consecutive velocity fields in the following manner: TI(x, y, t) = 1 The general mapping procedure into laminar and turbulent states is visualized in Fig. 3. The flow is moving along the airfoil's surface (gray) in the x-direction which is chordwise. In the horizontal color-coded plane, the LSB can be clearly identified as a high level of turbulence intensity. Mainly a velocity magnitude close to zero within the LSB contributes to the sharp rise of TI, while velocity fluctuations increase only slightly. It is important to notice that the term "laminar separation bubble" is misleading at this point. In the context of the present work, the LSB is revealed by a turbulent flow region, whereas it is referred to as laminar in airfoil aerodynamics. This is due to the fact that velocity fluctuations within the LSB are small compared to fluctuations present in a fully turbulent boundary layer. With that in mind, we use a threshold value of the turbulence intensity TI th to distinguish between laminar and turbulent regions, called clusters. Each measurement point of the complete dataset is set to 0 if TI(x, y, t) < TI th (laminar phase) or 1 if TI(x, y, t) ≥ TI th (turbulent phase). The turbulent phase corresponds to the LSB whereas laminar phase identifies ambient flow.
The onset of the LSB is accurately determined by the evo-lution of laminar and turbulent clusters. The distribution of clusters is evaluated at each position x along the spatial and temporal dimensions y and t, as exemplary illustrated close to the LSB's onset by the bicoloured vertical plane in Fig. 3.

IV. UNIVERSAL PHASE TRANSITION INTO LAMINAR SEPARATION BUBBLE
The flow along the airfoil is more likely to separate from the surface (onset of LSB) the further downstream from the leading edge it has been evolved due to an increasing receptivity for perturbations. In analogy to well known transition results from flat plates [16], this is expressed in terms of local Reynolds number Re(x) = Re x = x u ∞ /ν kin , where ν kin denotes the kinematic viscosity. This locally changing Re x represents the state of the boundary layer and must not be confused with the fixed Re chord reflecting the experiment globally. Therefore, Re x is used as the control parameter in the framework of DP for the present experiment. Notice that one HSPIV snapshot is composed of the control parameter Re xaxis and the spatial dimension y of our percolation model. The time dimension t in this percolation model corresponds to the physical time in the experiment.
In accordance with DP, a critical value Re c of the Reynolds number Re x is observed, where phase transition to LSB takes place. In general, this value is characterized by turbulent clusters merging into one infinitely connected cluster. Figure 4a shows three illustrative extracts containing laminar and turbulent clusters in space y and time t. While below/above the critical Reynolds number laminar/turbulent clusters dominate, at Re c the number of laminar and turbulent clusters span the entire system in both directions.
The DP phase transition from the laminar boundary layer into the LSB is characterized by three critical exponents at  Re c . The first exponent, β, describes the critical behavior of the so-called turbulent fraction, as a function of the reduced Reynolds number, ε := (Re x − Re c )/Re c , and a proportionality factor, ρ 0 . The mean fraction is determined from turbulent cells over space y and time t for each value Re x . The other two critical exponents, ν ⊥ and ν , characterize the diverging correlations of cluster sizes at the Re c in space and time. According to the so-called hyperscaling relation, µ = 2 − β/ν [17], the correlation lengths are univocally expressed by the transverse and longitudinal fractal dimensions, µ ⊥ and µ . These are defined as the expo-  and L represent the size of laminar clusters measured in the spatial and temporal directions, y and t respectively. In analogy of the DP analysis presented by [11] and [12], we estimate from our measured data first the critical Reynolds number Re c and subsequently the set of critical exponents (β, µ ⊥ , µ ). The results will be discussed in the following Sec. V. The turbulent fraction over local Reynolds number is shown for three descriptive thresholds TI th ∈ [0.8, 0.93, 1.42] in Fig. 4b. Phase transition occurs within the narrow range of values of Re x where an abrupt increase of the turbulent fraction is observed separating a laminar phase (ρ ∼ 0) from a turbulent phase (ρ ∼ 1). For the whole DP analysis, it is important to determine precisely the critical value of Re c . Re c and β can be obtained simultaneously by a best fit according to Eq. (2). The measurement points used for the best fit are selected in two steps. At first, the derivative of the turbulent fraction with respect to the Reynolds number dρ dRex is computed, as shown in Fig. 5a. Considering the largest value of the derivative together with the derivative values at the two nearest measured Reynolds numbers (red stars in Fig. 5a), a parabola with negative concavity is defined, whose mathematical maximum is taken as an initial estimate of the critical value. While starting the best fit at this estimated location, in step two, only measurement points are taken into account meeting 0.001 < ε < 0.01, in accordance with the experimental spatial resolution. As shown in Fig. 4b, this procedure yields values of Re c = 1.061(9) × 10 5 , β = 0.28 (4) for TI th = 0.8, Re c = 1.063(6) × 10 5 , β = 0.28(0.05) for TI th = 0.93 and Re c = 1.07(1) × 10 5 , β = 0.28(5) for TI th = 1.42. In comparison, the predicted theoretical value of 1 + 1D directed percolation is β DP = 0.276.
In order to evaluate the validity of the determined characteristic values Re c and β, the rescaled turbulent fraction,ρ, defined asρ is shown in Fig. 5b as a function of the reduced Reynolds number using the theoretical value of β DP . In this representation, scaling is readily apparent as a horizontal line equal to unity. Dashed lines also show the rescaled fraction obtained for slight deviations from the exponent, β DP ± 0.05.
To estimate the other two critical exponents describing the fractal dimensions, the distributions of the sizes of the laminar clusters at the critical Reynolds number have to be known. Based on our finite experimental resolution, we have no direct access to these critical distributions. To overcome this problem we take the corresponding distributions at the two Reynolds numbers close to Re c and sum them up in a weighted manner using the normalized inverse of their distances a 1 and a 2 to Re c (cf. inset of Fig. 5a).
While the procedure of estimating the critical size distributions is not standard, it results in robust estimates of spatial critical exponents µ ⊥ . The distributions of cluster sizes in space and time at the critical Reynolds number are shown in Figs. 6a and 6b, respectively. The evaluations are performed with same three threshold values TI th as previously. Each best fit for the spatial fractal dimension (smallest three cluster sizes are disregarded) yields µ ⊥ = 1.75 (9) and compares well with the theoretical value predicted by directed percolation, µ ⊥,DP = 1.748. For the temporal direction, experimental data is inconclusive so that a best fit is not performed. The theoretical exponent, µ ,DP = 1.84 is shown in Fig. 6b as a dashed line.
As a last step in our data analysis we investigate the effect of the choice of the threshold values TI th for this phase transition. Beyond the three threshold values considered up to now, Fig. 7 shows the results for two of the three critical exponents with 95% confidence intervals for 200 values of TI th covering a range of 0.5 < TI th < 2.5. The theoretically predicted values of the critical exponents are shown by horizontal dashed lines. It is evident that the values of β and µ ⊥ exponents are robust against the change of the turbulence intensity threshold for TI th > 0.64. A variation of a factor of 5 in the turbulence intensity threshold implies a variation of less than 2% of the critical Reynolds number. Thus, we conclude that the critical Reynolds number marking the onset of the LSB, depends only slightly on TI th .

V. DISCUSSION
The main aspect of our work is to show evidence that the onset of an LSB on an airfoil is linked to directed percolation. This evidence is shown by estimated critical exponents. It is well known that this estimation is very sensitive to the choice of the critical point, or here, the critical Reynolds number. To show the quality of our estimate, we present compensated plots in Fig. 5b. Another point supporting our choice of Re c is that the cluster size distributions become exponentially shaped (not shown here) when cluster size distributions are computed at Reynolds numbers which differ by less than 1000 from the critical value. For the critical exponents of the turbulent fraction (β) as well as the transverse fractal dimension (µ ⊥ ), Figs. 5b, 6a and 7 show a robust accordance with the predicted values of 1 + 1D directed percolation. This also holds true if the confidence intervals for our error estimation are put into question. The results on longitudinal fractal dimension (µ ) are inconclusive but can be taken as consistent with (or not contradicting) 1 + 1D directed percolation, see Fig. 6b. Particularly cluster distributions in time are known to be more difficult to estimate as these seem to suffer more from finite size effects and experimental uncertainty.
The obtained scaling ranges for our systems are not much more than one decade. Definitely, it would be desirable to investigate finite size scaling. In contrast to numerical investigations, we are limited by our experimental setup. Neither   the spatial dynamic range of HSPIV nor the size of the wind tunnel can easily be changed. Taking also into account that estimation of power laws may lead to biased estimation and pitfalls [18], our experimental results indicate, to the best of our belief, astonishingly consistent results with the 1 + 1D directed percolation. We also see that the quality of our results are comparable with those of other groups, like the work [12] on directed percolation in a channel flow.
Based on these findings, it is now possible to introduce an alternative procedure to determine the onset of the LSB. Under the assumption of 1 + 1D DP holding true, the critical Reynolds number can be determined following the compensated representation shown in Fig. 5b, but now fixing β DP to its theoretical value and varying Re c , For a variety of possible critical Reynolds numbers,ρ is shown as a function of ε and Re c in Fig. 8a. While increasing Re c , the slope,ρ = ∂ρ/∂ε, changes from positive to negative. In analogy to Fig. 5b, the horizontal line inρ represents the case where DP properties are found. Thus, from the zero-crossing ofρ the best estimation Re c is obtained. In comparison to the estimation of the critical Reynolds number by fitting a power law about the maximum derivative of ρ (see Figs. 4b and 5a), the uncertainty in Re c decreases by two orders of magnitude from Re c = 1.063(6) × 10 5 to Re c = 1.06302(2) × 10 5 for TI th = 0.93. This shows that the concept of DP enables to determine the critical Reynolds number for the LSB with very high precision. While such high precision concepts are very rare in fluid mechanical research, the introduced procedure may serve as a new benchmark.

VI. CONCLUSION
This work presents a first experimental evidence that directed percolation, as a more general concept, is also valid for practical relevant aerodynamics, namely the flow over the suction side of an airfoil. Our work has been inspired by recent achievements in fundamental turbulence research that link the onset of turbulence to directed percolation phase transition. In comparison to the flow situations investigated up to now, the flow over an airfoil changes its Reynolds number along with its stability while evolving downstream. In this sense, a new kind of spatially dependent DP is present, for which an adiabatic approximation has been anticipated.
Applying a bond directed percolation model to characterize transition from a laminar boundary layer into an LSB on an airfoil, one obtains values for the critical exponents consistent with those in 1 + 1D directed percolation. The physical implication of this universality class indicates that the boundary layer at the LSB's onset is slender compared to the dimension of the LSB and, thus, flow instabilities cannot spread perpendicular to the surface. As an important aspect for practical applications, with the assumption of an 1 + 1D directed percolation, a new method is introduced to determine the transition point into the LSB with very high precision of better than 10 in Reynolds numbers. This is for fixed TI less than 1 per mill or for all TI less than 1 per cent of the cord length. For our profile, the precision exceeds the optical resolution of 0.4 mm.
Since instabilities like the LSB have essential impact on the performance of airfoils, it is of great importance to know how and where they emerge. From the findings of this work, new directions for future applications and investigation are now open. First, in CFD, several frameworks need a model that delivers the location of transition [19] and, thus, the LSB's onset. Since the directed percolation framework is now shown to be able to retrieve a consistent determination where the LSB onset is located, it can be used to validate current transition models. Additionally, the temporal evolution of the LSB's onset parameterized by DP might be of use for unsteady loworder models [20] which often couple details of the boundary layer flow with CFD approaches used for the ambient flow [21,22].
Second, in a lot of applications, transition to turbulence is a phenomenon that needs to be controlled properly, also when taking place within an LSB. For instance, vortex generators are applied to rotor blades of wind turbines close to the location of an LSB in order to avoid aerodynamic instabilities that cause high fatigue loads. Nowadays, this is done based on efficient engineering models used in the design process of a wind turbine. As revealed by different studies [23,24], these models cannot account for nonlinear flow behavior that inherently governs the emergence of LSBs. Particularly when an LSB emerges very close to the natural onset of laminar-turbulent transition, models fail to correctly predict both phenomena. The more reliably such situations can be identified and characterized, the better vortex generators can be applied to rotor blades, resulting in more efficient wind turbine operation along with reduction of destructive aerodynamic loads.
Third, the wind tunnel experiment here investigated considered an airfoil subjected to a constant ambient inflow. In reality, airfoils of a plane or a wind turbine are subjected to non-stationary inflows. Assuming that the non-stationarity of such flow occurs at a smaller time scale than the time scale needed for the LSB onset to take place as a percolation transition, the concept of directed percolation has the potential to derive a model for the (non-stationary) dynamics of the LSB's onset in time.
All in all, results of the present study indicate a more comprehensive significance of percolation models in fluid mechanics beyond fundamental laminar-turbulent transition phenomena.