Numerical simulation of flow characteristics behind the aerodynamic performances on an airfoil with leading edge protuberances

ABSTRACT This article presents a numerical investigation of the effects of leading-edge protuberances on airfoil stall and post-stall performance. An improved delayed detached eddy simulation (IDDES) method was adopted. As a result, to clarify the effects of ‘bi-periodic’ phenomenon around stall region, it was found that the flow separation at troughs was the main inducement of aerodynamic lift degradation within pre-stall regime and the flow pattern where vortices diverged was predominant. It was also found that the variations in flow patterns led to the gentle stall process. Furthermore, to study the statistical characteristics of unsteady vortex shedding, corresponding spectrum characteristics were also analyzed from another perspective, suggesting that the vortex shedding frequency was higher where vortices converged. Eventually, the improved performances of tubercled airfoil within post-stall regime could be attributed to the strong streamwise vortices generated by the leading-edge protuberances. Deploying the methods of vortex dynamics, the generation and evolution of the streamwise vortices were depicted. It turned out that the primary and secondary vortices were induced by spanwise pressure gradient at airfoil surface; meanwhile, vortex stretching played a key role in primary vortex evolution, which initially enhanced the strength of vortices corresponding to the acceleration of streamwise velocity.


Introduction
The research interests in the civil and military utilities at low Reynolds (Re) number with an order of magnitude of 10 5 , e.g. small wind turbine, unmanned aerial vehicles (UAV) and micro air vehicles (MAVs), have been greatly motivated in the last decade (Mueller & DeLaurier, 2003;Wright & Wood, 2004). Nevertheless, the engineering structures will be subject to poor aerodynamic performances caused by flow separation or stall on the airfoils, leading to shortened fatigue lives and even serious accidents at low Reynolds number. Therefore, it is important to find an effective flow control method.
Recently, a new kind of passive flow control technique, i.e. 'leading-edge protuberance', has received many interests, which was inspired by the work of marine biologists who studied the morphology of humpback whales' pectoral flippers (Bolzon, Kelso, & Arjomandi, 2015;Fish & Battle, 1995). It was reported that the airfoil with leadingedge protuberances could provide more aerodynamic lift than the smooth one within the post-stall region and the stall process was rather gentle, which brought the benefit of longer fatigue lives (Custodio, 2007;Drop-kin, Custodio, Henoch, & Johari, 2012;Hansen, Kelso, & Dally, 2011;Zhang, Wang, & Xu, 2014).
Deploying this bionic means, many investigations have been conducted concerning with the improved aerodynamic performance at low Reynolds number (Johari, Henoch, Custodio, & Levshin, 2007;Lin, Lam, Zou, & Liu, 2013;Pedro & Kobayashi, 2008;Skillen, Revell, Pinelli, Piomelli, & Favier, 2015). However, the underlying mechanisms are still not clearly understood. In fact, it is widely accepted that the enhanced momentum exchange accompanied by the vortices induced by tubercles play a key role in the improvement of aerodynamic performances (Hansen et al., 2011;Zhang et al., 2014), and the energization process is more obvious around the peak sections of tubercles. Nevertheless, van Nierop, Alben, and Brenner (2008) claimed that it was not possible for the tubercles to act in such a way, since the wavelength and amplitude were much larger than the boundary layer thickness. To clarify this issue, Hansen et al. (2011) proposed the effective height of tubercles to quantitatively evaluate the effect of tubercles. Furthermore, Zhang et al. (2014) found that the ratio of effective height of tubercles to boundary layer thickness should lie in the interval between 0.1 and 0.5, similar to micro vortex generators. Besides, other hypotheses of tubercle mechanism were also proposed. For example, Custodio (2007) postulated that the counter-rotating streamwise vortices induced by protuberances provided additional vortex lift due to the low pressure in core of these vortices, similar to that observed on a delta wing; meanwhile, analogous conclusion was also drawn by Dropkin et al. (2012) in which the improved aerodynamic lift was attributed to the low-pressure pockets resulting from the leading-edge troughs.
Recently, 'bi-periodic' phenomenon had attracted much attention and was thought to be another possible physics to effectively affect airfoil stall performance. In fact, accompanied by the non-equilibrium of velocity and static pressure, the path lines and vortex structures converged and diverged in adjacent trough sections of tubercles, creating a bi-periodic pattern around stall region (Dropkin et al., 2012;Zhao, Zhang, & Xu, 2015). Although several experimental studies have been conducted and some valuable findings have been obtained concerning with 'bi-periodic' phenomenon (Johari et al., 2007), more information could be achieved through numerical studies since details of flow field could be depicted (Camara & Sousa, 2013;Dropkin et al., 2012;Malipeddi, Mahmoudnejad, & Hoffmann, 2012;Rostamzadeh, Hansen, Kelso, & Dally, 2014;Zhao et al., 2015). First, Reynolds-averaged Navier-Stokes (RANS) methods have been proved to have the capacity of capturing 'bi-periodic' phenomenon. Consequently, Dropkin et al. (2012) found that 'bi-periodic' phenomenon dominated the flow field around stall region and the flow pattern reverted to periodic status within post-stall region. Meanwhile, accompanied by the occurrence of 'bi-periodic' phenomenon, Rostamzadeh et al. (2014) found that the attached flow resided behind the peaks or alternate troughs actually improved the aerodynamic performance right after airfoil stall. However, the details of unsteady turbulent flows were missed due to the time-averaged property of RANS methods and the massive flow separations could not be accurately depicted (Shur, Spalart, Strelets, & Travin, 2008). To overcome this, the detached-eddy simulation (DES) method, a combination of RANS and LES methods, has been implemented. As a result, Camara and Sousa (2013) found that 'bi-periodic' phenomenon could only be observed at high angle of attack (18°) around stall region for specific configuration, consistent with the experimental findings. Simultaneously, Malipeddi et al. (2012) proposed a hypothesis of the inducement of 'bi-periodic' phenomenon based on the observation of vortex structures, suggesting that the interactions between vortices of neighboring valleys triggered 'bi-periodic' phenomenon.
Although a lot of work has been conducted on the performances of tubercled airfoil, the influences of 'biperiodic' phenomenon were still unclear, which would have great impact on stall performance; meanwhile, the generation and evolution mechanisms of predominant streamwise vortices within post-stall region were also unknown, which would improve the overall performances.
In present study, an improved delayed detached eddy simulation (IDDES) method was utilized to simulate the flow field of NACA 63 4 -021 airfoil with and without protuberances. Moreover, the effects of 'bi-periodic' phenomenon on airfoil stall performance and relative physics were discussed in detail. As a result, it was found that the lift degradation within pre-stall regime could be attributed to the flow separation at the leadingedge of troughs and the influence of the flow pattern where vortices diverged was predominant. Accordingly, the complicated variations in the flow patterns led to little change in total lift coefficient around stall region. Furthermore, the unsteady characteristics within time and spectrum domains were also analyzed, suggesting that the vortex shedding frequency was higher where vortices converged. Eventually, the improved performances of tubercled airfoil within post-stall regime could be attributed to the inducements of streamwise vortices and thus the energization process around peak sections. Deploying the methods of vortex dynamics, the generation and evolution mechanisms of the predominant streamwise vortices were depicted. It turned out that the primary and secondary vortices were induced by spanwise pressure gradient at airfoil surface. Meanwhile, vortex stretching played a key role in primary vortex development, which initially enhanced the strength of vortices corresponding to the acceleration of streamwise velocity.

Turbulence modeling
As mentioned in Section I, unsteady detached vortex motion was inevitable in the flow field of tubercled airfoil; RANS method was thus not applicable. Considering feasibility and accuracy, an IDDES method was adopted in present study, which has recently become much favored in the study of the unsteady and geometry-dependent separated flows. In practice, IDDES method based on kω SST turbulence model was adopted (Shur et al., 2008), which can be constructed through modifying the dissipation term of the turbulent kinetic energy equation (k-equation). After introducing a length scale, L hybrid , the turbulent kinetic energy and specific dissipative rate equations can be given in tensor form as (1) Here, S ij = (1/2)(∂u i /∂x j ) + ∂u j /∂x i is the strain tensor and τ ij = 2μ t S ij − (2/3)μ t (∂u i /∂x i )δ ij is the stress tensor. F 1 is the blending function (Shur et al., 2008). The IDDES length-scale can be implemented as (2) r dt (= μ t /ρκ 2 d 2 · max{[ ij (∂u i /∂x j ) 2 ] 1/2 , 10 −10 }) is borrowed from Spalart-Allmaras turbulence model, where d is the distance to wall surface and κ = 0.41. f B = min{2 exp(−9α 2 ), 1.0}, α = 0.25 − d/h max , and h max is set equal to the maximum local grid spacing. The elevating-function f e (= max{(f e1 − 1), 0}f e2 ) is aimed at preventing the reduction of the RANS Reynolds stresses, the details of which is not presented here for brevity (Shur et al., 2008).
Actually, above IDDES method includes two branches, delayed detached-eddy-simulation method (DDES) and wall-modeled large-eddy-simulation method (WMLES). Their coupling ensures a favorable response of the combined model as DDES or WMLES depending on the inflow (or initial) conditions used in the simulation. Here, DDES method is actually a combination of RANS and LES method.

Numerical methods and verifications
In practice, a symmetric total-variation-diminishing (STVD) scheme was adopted to discretize the flux in present study, which can be written as where the inviscid flux is discretized by the sixth-order symmetric scheme. The numerical dissipation is the original Roe's dissipation with fifth-order WENO interpolation. In detail, A inv is the matrix of Roe average, q R and q L are obtained through the fifth-order WENO interpolation. Then, this STVD scheme can be shortened as S6WENO5 and Xiao, Jian, Luo, Huang, and Fu (2013) discussed the influence of blending function φ on DES results. It turned out that S6WENO5 schemes with constant-φ (0.12) and adaptive-φ had similar performances, but S6WENO5 scheme with constant-φ was empirical, and small constant-φ could not generally suppress the numerical oscillations. Consequently, adaptiveφ was adopted in present study, which was closely related to local flow field. where Consequently, the net effect of above formulas is that φ approaches zero in the separation region and it is close to unity in the irrotational region and near the wall. Therefore, the robustness of present numerical method could be guaranteed without the impairment of turbulent structures.
In addition, a modified fully implicit lower-upper symmetric Gauss-Seidel (LU-SGS) with Newton-like subiteration in pseudotime was taken as the timemarching method. To overcome the stiffness of timedependent equations and accelerate the convergence, preconditioning algorithm was adopted here (Weiss, Maruszewski, & Smith, 1999). The approach was in parallel algorithm using domain decomposition and messagepassing-interface strategies for the platform on PC clusters, and the computations were all based on in-house solver UNITs (Xiao et al., 2013).
The flow fields over the smooth and tubercled airfoil with mean chord length c = 100 mm and span length s = 100 mm were investigated. The wavy leading edge had an amplitude of A = 12 mm and a wavelength of λ = 25 mm (Figure 1a), leading to a smooth stall process among various geometric shapes (Johari et al., 2007). Free-stream velocity (U ∞ ) and Reynolds number based on chord length (Re c ) were 30 m/s and 2 × 10 5 , respectively. Figure 1(b) demonstrates the sketch of computational domain and Figure 1(c) illustrates a slice of the computational mesh. Wall surfaces were regarded as noslip and adiabatic for viscous simulation and the lateral sides were regarded as periodic conditions in the spanwise direction. Moreover, the one-dimensional Riemann characteristic analysis was employed to construct a non-reflection boundary condition at far-field boundary, which was placed 20c upstream and downstream of the airfoil, satisfying the requirement of computational domain under similar circumstances (Rostamzadeh et al., 2014;Yoon, Hung, Jung, & Kim, 2011). To properly resolve the wall boundary layers, mesh nodes were clustered near wing surfaces using a geometric expansion and a grid independence test has been conducted for tubercled airfoil. The amounts and distribution of grids are given in Table 1.
As a result, the numerical work was validated using previous experimental results (Johari et al., 2007;Zhang et al., 2014). To this end, six angles of attack (α) were selected, i.e. 6°, 12°, 18°, 21°, 23°and 45°. In addition, the time step was set to 1.667 × 10 −5 s, and the solution proceeded until 0.3s. To obtain the time-mean flow properties and prevent any effects of the initial flow conditions, the first 1500 time steps were removed from the results.
The variations of lift and drag coefficients with α are displayed in Figure 2. Here, lift coefficient was defined as C L = 2L/(ρU 2 ∞ sc), and drag coefficient was defined as C D = 2D/(ρU 2 ∞ sc). L and D were aerodynamic lift and drag; meanwhile, s and c stood for the span and chord of airfoil. As a result, stall process could be identified for smooth airfoil with an abrupt change of C L and C D around 21°, which was coincident with experimental results. Meanwhile, for tubercled airfoil, it could be observed that the results of moderate and fine meshes were in good agreements with experimental results. In fact, the error was constrained within 5% (Table 1). Accordingly, moderate mesh with y + ≈ 1, x + ≈ 20, z + ≈ 40 and growth ratio ≈ 1.12 was selected, satisfying the delicate requirement recommended (Spalart, 2001). Furthermore, the distribution of y + at peak and trough is shown in Figure 2(c) when α = 6 • and 45 • , verifying the quality of present mesh system.   Moreover, the accuracy of present method was also validated through the quantitative comparison of velocity profiles with experimental data (Zhang et al., 2014) at 18°, as shown in Figure 3. In general, the numerical results were in acceptable agreement with experimental results, although the velocity gradients of numerical results were smaller due to numerical dissipation. Besides, the corresponding results of 45°were also available, which are not presented here for brevity.
Eventually, time-averaged results were qualitatively compared with corresponding RANS results (Cai, Zuo, Liu, & Wu, 2015) for tubercled airfoil. The 3D vortex structures identified by Q criterion (Q = 1/2( ij ij − S ij S ij ) 0.5 ) and colored by streawise vorticity ( 23 )are given in Figure 4. Analogous results were achieved, although some distinctions could be observed as the wavelength of tubercles was larger in the work of Cai et al. Actually, 'bi-periodic' phenomenon could be identified in  both results; in other words, the vortex structures converged and diverged around adjacent troughs, creating a bi-periodic pattern around stall region. Meanwhile, the streamwise vortices resided in the vicinity of tubercles in both cases.

Descriptions of 'bi-periodic' phenomenon
Using the aforementioned numerical methods, the flow patterns behind the 'bi-periodic' phenomenon were studied. Figure 5 gives the corresponding instantaneous flow structures using Q criterion (Q = 1/2( ij ij − S ij S ij ) 0.5 ), from which it could be observed that the processes of vortex evolution were rather complicated.
At α = 6 • , the flow pattern was periodic from one trough to the next ( Figure 5a); nevertheless, at α = 18 • , the flow pattern was no longer periodic and 'bi-periodic' phenomenon emerged, as shown in Figure 5(b). The fact that vortices dispersed at one trough (trough-A) as soon as it left the leading edge indicated that there was flow interaction between neighboring troughs, consistent with the findings of Malipeddi et al. (2012). At the neighboring trough (trough-B), vortices were confined to the valley of current trough and above flow pattern repeatedly appeared in spanwise direction. For convenience, the trough where vortices dispersed was named trough-A, and the neighboring trough was named trough-B (Figure 5b).
In fact, the patterns of vortex evolution at α = 21 • and 23 • were similar to that of 18 • , as shown in Figure 5(c and d), and the non-equilibrium between neighboring troughs became more obvious. Evidently, it could be inferred that the occurrence of 'bi-periodic' phenomenon was closely related to α, consistent with previous findings (Camara & Sousa, 2013;Dropkin et al., 2012).
To understand the formation of the 'bi-periodic' phenomenon, we carried out some preliminary studies as displayed in Figure 6 (Zhao et al., 2015). Here, M denotes the Mach number and non-dimensional pressure p(= p/ρ ∞ U 2 ∞ ). As we can see, vortex structures was periodical at initial time instant, e.g. t = 0.0067s, and no difference in the flow patterns between neighboring troughs existed, as shown in Figure 6(a). Then the vortices near trough-A section crossed the neighboring peak and forced the fluid to squeeze into the region near trough-B (Figure 6b) at later time instant, e.g. t = 0.023s. This would lead to an increase of velocity and a decrease of static pressure at trough-B section. To have an intuitive understanding, the slices of velocity and pressure field in x-z plane at y = 0.11c were provided in Figure 6(c and d), respectively. The flow with higher momentum entered the area behind trough-B (bounded by red dash line), while lower static pressure was also observed (bounded by blue dash lines). The non-equilibrium of velocity and pressure would be gradually exacerbated until biperiodic pattern was reached. A more rigorous study based on stability analysis concerned with small disturbance equation might further reveal the mechanism of 'bi-periodic' phenomenon formation.

Effects of 'bi-periodic' phenomenon on stall performance
Interests were more motivated to investigate the flow feature after the bi-periodic pattern was stabilized. Clearly, 'bi-periodic' phenomenon remarkably changed the flow pattern above the airfoil suction side and this would undoubtedly affect the overall aerodynamic performance, especially around the stall region. For this purpose, the corresponding flow field was first analyzed in detail. The flow field characteristics of the tubercled airfoil would lead to the corresponding changes of C L and C D within stall region in Figure 2. To prove this, Figure 7 illustrates the distributions of pressure coefficients (C p = (p − p ∞ )/(ρ ∞ U 2 ∞ )) along x-direction. In general, compared with the results of smooth airfoil, aerodynamic suction around the leading-edge of troughs decreased due to the flow separation within pre-stall region (Figure 7a-c). Furthermore, 'bi-periodic' phenomenon aggravated the lift degradation within pre-stall regime, which removed the aerodynamic suction at the leading-edge of trough-A, especially at α = 18 • and 21 • (Figure 7b and c). Meanwhile, the pressure drag was thus enlarged due to above mechanisms. On the other hand, the advantages of turbercled airfoil emerged right after airfoil stall at α = 23° (Figure 6d), since the aerodynamic suction at the leading-edge of smooth airfoil vanished and the aerodynamic lift dramatically reduced. Together with above results, the flow patterns within pre-stall region (α = 18 • ) were analyzed using 3D timeaveraged static pressure coefficient and streamlines, as shown in Figure 8(a). Flow separations and attachments were indicated by blue and red dashes, respectively. As a result, it could be observed that flow separation was more severe at the leading edge of trough-A than that of trough-B due to 'bi-periodic' phenomenon; therefore, aerodynamic suction almost vanished around trough-A, consistent with the findings in Figure 7(b). Moreover, after flow reattachment, flow became separated again in a short distance around trough-A section; while the location of corresponding flow separation lagged behind around trough-B. Evidently, the massive flow separation at troughs probably imposed negative effects on aerodynamic forces within pre-stall region and it tended to be more remarkable around trough-A due to 'bi-periodic' phenomenon.
Simultaneously, right after smooth airfoil stall, corresponding results of α = 23 • are demonstrated in Figure 8(b). It turned out that flow separation was more severe compared with α = 18 • . Although flow separation dominated the flow field around trough-A section, flow attachment covered a large portion of the suction side around peak sections, which should be attributed to the inducement of tubercles. Therefore, the improved aerodynamic performances of tubercled airfoil right after airfoil stall could be interpreted as flow became completely separated for smooth airfoil (Figure 7d).
To further interpret the flow attachment around peak sections, a streamwise slice (x = 0.2c) of pseudostreamlines colored by Mach number, in Figure 9, demonstrated that the induced velocity of the dominant primary vortices, leading to attached flow near peak sections. As a result, the enhanced momentum transfer mechanism, brought by the presence of streamwise vorticity, was believed to account for higher lift coefficients achieved by the tubercled wings right after stall. Meanwhile, 'bi-periodic' phenomenon could also be detected as the vorticity strength was higher near trough-B.
In addition, as shown in Figure 2, lift coefficient of tubercled airfoil leveled off around stall region, which was a favorite feature for fatigue lives. It could be interpreted by local flow patterns around troughs and peaks. In detail, as flow attachment was remarkable around peak sections, local aerodynamic lift increased with angle of attack. On the contrary, flow separation dominated the flow field around troughs, therefore the degradation of local aerodynamic lift deteriorated along with the increase of angle of attack. In general, these variations in the flow patterns led to little change in total lift coefficient around stall region.
To understand this, Figure 10 gives the C l -distributions along z-direction, and it could be inferred that the lift coefficients of α = 18 • , 21 • and 23 • were approximately equal to 0.85 in consideration of the contributions of these spanwise sections (i.e. Peak, Middle, Trough-A and Trough-B). Here, 'M' indicates the locations of middle sections between troughs and peaks. Actually, the variations of lift coefficients at peaks and troughs supported the aforementioned hypotheses, as lift coefficient of peak section increased with angle of attack and an opposite tendency emerged at trough section.
Meanwhile, the influence of 'bi-periodic' phenomenon certainly extended to the middle sections, as the lift coefficients of the middle sections adjacent to one peak were different. In fact, the lift coefficient of the middle section adjacent to trough-B was close to that of peak; while the lift coefficient of the middle section near trough-A was relatively lower due to the influence of the massive flow separation near trough-A section, which shared the same variation tendency with trough-A and enhanced the stability of stall process.  To further evaluate the effects of different spanwise sections, it could be observed that the lift deficits of troughs and the middle sections adjacent to trough-A should be responsible for the lift degradation within prestall regime and the negative influences of trough-A were more obvious (α = 18 • and 21 • ). On the other hand, contributions of peaks and middle sections evidently improved the overall aerodynamic lift after airfoil stall (α = 23 • ). Actually, the contributions of the middle sections adjacent to trough-B were comparable to those of peaks since the extreme value of negative static pressure resided near the leading edge of trough-B (Figure 8b). In addition, above conclusions were consistent with the descriptions of flow patterns in Figures 7 and 8.
Furthermore, to analyze the fluctuation characteristics, distributions of non-dimensional resolved turbulent kinetic energy (k resolved = (u 2 + v 2 + w 2 )/2) at different spanwise sections of smooth and tubercled airfoils are shown in Figure 11 (α = 23 • ). Here, u , v and w represented the components of fluctuating velocity. As a result, it was found that the turbulent fluctuations were relatively weak for smooth airfoil as shown in Figure 11(a). At peak section, the turbulent fluctuations were weak near the leading edge and flow became disordered around x = 0.2c due to the impact of neighboring troughs (Figure 11b). Meanwhile, 'bi-periodic' phenomenon was reflected in Figure 11(c and d). In particular, the influence range of turbulent flow at trough-A was greater; while the stress of turbulent fluctuations was stronger near the leading-edge of trough-B.
Based on above analyses, the turbulent characteristics within time and spectrum domains were analyzed and the locations of samples were selected near the leading edge of troughs, as shown in Figure 12(a). Actually, taking trough-A for example, there were two samples named sampleA1 and sampleA2, respectively. SampleA1 and sampleB1 share the same xand ycoordinates, i.e. x 1 and y 1 . Here, x 1 = 0.18c, y 1 = 0.232c, x 2 = 0.135c and y 2 = 0.075c, with origin at the smooth airfoil leading edge (Figure 12b and c). In fact, sampleA1 and sam-pleB2 resided near the paths of shear layer and vortex   shedding according to the streamlines; while sampleA2 and sampleB1 were selected for comparison.
As a result, the histories of non-dimensional pressure (p/(ρ ∞ U 2 ∞ )) fluctuations exhibited great unsteadiness as shown in Figure 13. Comparing the results of sam-pleA2 and B2, it could be observed that high-frequency oscillations existed atsampleB2 and the amplitude was larger (Figure 13b), leading to the high-resolved turbulent kinetic energy near the leading edge of trough-B (Figure 11d). Furthermore, the distinctions between sampleA1 and sampleB1 were more obvious. In particular, turbulent fluctuations at sampleA1 had the largest amplitude among the samples since it resided near the path of vortex shedding (Figure 13a); on the contrary, turbulent fluctuations almost vanished at sampleB1 since it was beyond the turbulent flow region.
On the other hand, to study the turbulent characteristics within spectrum domain and capture the primary frequency of vortex shedding, the power spectral density (PSD) information of above pressure fluctuations is demonstrated in Figure 14. Consistent with the findings in Figure 13(a), a primary frequency of vortex shedding could be observed around 600 Hz at sam-pleA1 (Figure 14a). Corresponding results at sampleB2 was also demonstrated in Figure 14(d), and a primary frequency could be observed around 4000 Hz, indicating a vortex shedding process with higher frequency compared with trough-A. Meanwhile, no obvious primary  frequency could be observed at sampleA2 and sampleB1 (Figure 14b and c).

Analyses of performances within post-stall region
Likewise, the flow pattern within post-stall regime was analyzed using the results in Figures 15-19. Figure 15 displayed the instantaneous flow structures over tubercled airfoils using Q criterion at α = 45 • . Obviously, the separation over tubercles was not apparent, suggesting the improvement in the airfoil performances; meanwhile, single periodic vortex structures appeared between neighboring troughs. To further identify the dominant flow structure within post-stall region, the time-averaged non-dimensional vorticity and its 3D decompositions over tubercled airfoil were shown in Figure 16. Evidently, it could be observed that the streamwise vortices dominated the flow field above tubercle surface, which probably delayed the flow separation over tubercles; meanwhile, the influences of longitudinal (ω y ) and spanwise vortices (ω z ) were not notable. Therefore, the generation and evolution of streamwise vortices (ω x ) were of great importance within post-stall region, which would impact the overall performances. In detail, ω x comprised three parts, i.e. primary vortex, secondary vortex and shear layer (Figure 16b), analogous with a delta wing.
To further identify the effect of the dominant streamwise vortices, a streamwise slice (x = 0.2c) of pseudostreamlines are demonstrated in Figure 17, indicating that the induced velocity of the strong streamwise vortices actually injected flow with high momentum toward airfoil surface. As a result, boundary layers would have the capacity to resist strong adverse pressure gradient and flow separation would be delayed.
Moreover, deploying the vortex dynamic methods, the generation mechanisms of streamwise vortices could be depicted. To this end, boundary vorticity flux (BVF, σ = −ν( n · ∇) ω) was computed, which could theoretically denote the generation and viscous diffusion flux of vorticity (Wu, Ma, & Zhou, 2006). Actually, BVF was consist of the pressure component ( σ p = − n × ∇p) and the viscous component ( σ τ = 1/ρ n × ∇) × (μ ω)), among which σ p played a major role under current circumstances (Wu et al., 2006). Here, n is the normal vector of airfoil surface and p is the non-dimensional static pressure. Thus, the streamwise component of σ p (i.e., σ px = (∂p/∂z)n y − (∂p/∂y)n z ) could depict the generation of streamwise vorticity at airfoil surface. Clearly, the contour of σ px was illustrated to indicate the generation of shear layer, primary and secondary vortices (Figure 17a). Furthermore, σ px could be divided into two parts, i.e. σ px1 = (∂p/∂z)n y and σ px2 = (−∂p/∂y)n z , which represented the influences of spanwise and longitudinal pressure gradients on vortex generation. As a result, it could be observed that the primary and secondary vortices were induced by spanwise pressure gradient (Figure 17b), while the formation of shear layer was attributed to the combined effect of spanwise and longitudinal pressure gradients (Figure 17b and c).
Actually, the former one accounted for the effect of streamwise vortex stretching on ω x and the latter one described the production of ω x through vorticity skewing. Therefore, P1 + P2 could represent the overall effect of vorticity convection on ω x . It was found that the primary vortices were enhanced before x = 0.07c since P1 + P2 and ω x were of the same sign there ( Figure 18a); thereafter, the sign of P1 + P2 changed, leading to an impairment of primary vortices. In detail, vorticity stretching term played the major role    in primary vortex evolution because the magnitudes of P1 + P2 and P1 were close near primary vortices (Figure 18a and b); meanwhile, the contributions of vortex stretching in shear layer and secondary vortex evolution were almost offset by vorticity turning mechanism ( Figure 18c). Furthermore, the vortex stretching mechanism was closely related to the streamwise velocity gradient. As an acceleration of streamwise velocity took place near the primary vortices before x = 0.07c, ∂U/∂x was positive there. As a result, P1 = ω x (∂U/∂x) and ω x were of the same sign, resulting in the enhancement of primary vortices. Afterwards, a process of primary vortex impairment was triggered due to the elimination of positive velocity gradient.
As a result of above mechanisms, the distributions of non-dimensional static pressure and streamlines within post-stall region were demonstrated in Figure 19. For tubercled airfoil, the flow attachment was maintained until x = 0.5c around peak sections due to the energization process of primary streamwise vortices mentioned above, and the location of flow separation was consistent with the experimental results (Zhang et al., 2014). Meanwhile, flow separation took place at the leading edge of trough sections because of the adverse pressure gradient. Afterwards, flow reattached around x = 0.157c due to the impact of streamwise vortices and it became detached again beyond x = 0.4c. Simultaneously, the flow separations and attachments corresponding to the secondary vortices could also be detected.

Conclusions
A numerical investigation of the stall performances of NACA 63 4 -021 airfoil with protuberances and its smooth counterpart was presented at low Re c (2 × 10 5 ). An IDDES method was utilized to simulate the unsteady flow field. The associated flow field and inherent physics were analyzed, leading to the following conclusions: (1) The accuracy of present IDDES method was validated through quantitative comparison with experimental results, including aerodynamic coefficients, velocity profiles and time-averaged vortex structures. It was confirmed that the tubercled airfoil could provide more aerodynamic lift than the smooth one within the post-stall region and the stall process was rather gentle.
(2) Compared with the results of smooth airfoil, aerodynamic suction decreased due to the massive flow separation at troughs within pre-stall regime. Furthermore, 'bi-periodic' phenomenon aggravated the lift degradation, which removed the aerodynamic suction around the leading-edge where vortices diverged. On the other hand, the improved aerodynamic performances of tubercled airfoil after airfoil stall could be attributed to the induced velocity of the dominant primary vortices and thus the flow attachment around peak sections. (3) The gentle stall process could be interpreted by local flow conditions. As flow separation was mostly attached around peak sections, local aerodynamic lift increased with angle of attack. On the contrary, flow separation dominated the flow field around trough sections, therefore the local lift degradation deteriorated along with the increase of angle of attack. Meanwhile, 'bi-periodic' phenomenon further enhanced the stability of stall process with sophisticated flow patterns. (4) The distribution of resolved turbulent kinetic energy and the turbulent characteristics of time and spectrum domains were also analyzed. It was found that the vortex shedding where vortices converged had a higher frequency than that where vortices diverged; meanwhile, no obvious primary frequency could be observed away from the paths of vortex shedding. (5) Within post-stall region, the original detached poststall flow over the smooth airfoil was effectively impaired by the dominant strong streamwise vortices generated by the leading-edge protuberances, leading to the greatly enhanced airfoil performances. Moreover, the generation and evolution of predominant streamwise vortices were depicted with vortex dynamics results. In fact, the primary and secondary vortices were induced by spanwise pressure gradient at airfoil surface, while the formation of shear layer was attributed to the combined effect of spanwise and longitudinal pressure gradients. Meanwhile, vortex stretching played a key role in primary vortex evolution, which initially enhanced the strength of vortices corresponding to the acceleration of streamwise velocity. Afterwards, a process of primary vortex impairment was triggered due to the elimination of positive velocity gradient.
The study also suggested that a more rigorous study based on stability analysis concerned with small disturbance equation might further reveal the mechanism of 'bi-periodic' phenomenon formation and IDDES method might be improved considering flow transition. Eventually, the influences of wavelength and amplitude, together with the influences of airfoil shape will be comprehensively discussed in the future.