Direct numerical simulation of transitional flow past an airfoil with partially covered wavy roughness elements Open Access

Flow over rough surfaces occurs in a wide range of engineering applications, such as piping systems, marine vehicles, and urban atmospheric boundary layers.^1,2 The roughness of these surfaces can significantly affect the boundary layer development and, thus, influence the surface drag, heat, and momentum transfer, etc.³ The pioneering work of quantifying the roughness effect in pipe flows was led by Nikuradse,⁴ and Moody⁵ proposed the Moody diagram that characterized the friction factor as a function of Reynolds number and relative roughness. As the Reynolds number increases, three distinct flow regimes were identified: a laminar flow regime where roughness has a negligible effect on the friction, a transitional regime, and a fully turbulent regime where roughness becomes a significant factor once it exceeds a critical size. In the past decades, many investigations have been devoted to the roughness effects in canonical wall-bounded flows in simple geometries, for both laminar^6,7 and turbulent regimes.^8,9

The flow past bluff bodies is characterized by laminar-turbulent transition and boundary layer separation, and hence, often, to non-equilibrium effects (curvature and pressure gradient). Thus, the roughness effects on the flow behavior would be more complex. In fact, the roughness elements on the bluff body modify the surrounding flow structure, transition scenario, and separation behavior, affecting the aerodynamic performance¹⁰ and yielding significant departures from Nikuradse's initial proposal. Cheng et al.¹¹ performed a wall-resolved large-eddy simulation (LES) of flow past a grooved cylinder up to transcritical Reynolds numbers, with azimuthal sinusoidal-shaped groove height ranging from $D/16$ to $D/640$ (D is the diameter of the cylinder). The authors found that the pressure and skin friction coefficients and the near-wall field are significantly affected by the groove, with locally trapped separation bubbles formed as cavity flows within grooves. Thus, it is found that the drag coefficient can be increased or reduced depending on both roughness height and flow regime. Zhou et al.¹² experimentally investigated rectangular-shaped grooved cylinders in the subcritical regime with a fixed groove height of $0.05D$ and found that the roughness can reduce the drag coefficient and narrow the wake. Different types of roughness elements, such as dimples¹³ and irregular ice¹⁴ partially or entirely placed on the cylinder, have been investigated, which show different effects on the aerodynamic performance case by case.

Unlike the cylinder with constant curvature, the airfoil geometry may be considered a combination of large curvature (near the front stagnation part) and small curvature near the trailing edge portion of the airfoil, which will no doubt result in rich flow physics with roughness elements. The study of airfoils with roughness in different configurations is motivated by the analysis of surface imperfections due to manufacturing (gaps between metal plates, rivets, etc.) or ice accretions, which form with various shapes near the leading edge. The latter is typical when an aircraft enters a region containing supercooled water droplets under different meteorological and flight conditions.¹⁵ Based on representative ice geometries, Bragg et al.¹⁶ qualitatively classified ice accretions into four major categories: roughness ice, horn ice, streamwise ice, and spanwise ridge ice. They also observed that many ice shapes are complex, which do not purely fall into a single category but may have features representative of two or more. Most experimental investigations have been conducted, mainly quantifying aerodynamic performance degradation with ice accretions. The general consensus is that ice roughness, regardless of the ice roughness shape, results in increased drag and reduction in lift.¹⁷ However, there are minimal numerical investigations of the impact of accretions on airfoils (refer to the review paper by Stebbins et al.¹⁸). The main reason is that the ice shapes can vary widely and result in complex three-dimensional surfaces with intricate geometries, which makes the high-quality mesh generation challenging. Consequently, the mesh resolution would also dramatically increase, and the necessity for more computational resources would also increase. Meshing all the detailed roughness features is generally prohibitive for practical computational studies, especially for DNS, which resolves a full range of turbulent length scales. In industrial simulations, most calculations of the airfoil flow are performed by solving the Reynolds-Averaged Navier–Stokes (RANS) equations with turbulence models.^18,19 In those settings, the near-wall grid size is of order $0.01–0.02C$ (C is the chord length), which could not even resolve small-scale geometric characteristics.²⁰ 

Serson et al.²¹ investigated the effect of spanwise waviness elements on the leading edge of the NACA0012 airfoil with DNS and found that the wavy roughness elements may decrease or increase the lift coefficient, depending on the Reynolds number and flow regime (pre-stall or post-stall). Kumar et al.²⁰ studied effects of leading-edge roughness designed to mimic ice accretion with LES, and the ice accretion is modeled with cylindrical (uniform in the spanwise direction) and randomly distributed bump-like roughness elements. They found that the roughness elements can accelerate the transition to turbulence and produce different 3D structure development. Kim et al.²² investigated the spanwise wavy roughness effect on the airfoil–turbulence interaction by solving the compressible Euler equations and found that the wavy roughness can reduce the surface pressure fluctuation. Thomas et al.²³ investigated the development of Tollmien–Schlichting (TS) wave instabilities with the effect of streamwise sinusoidal surface waviness elements and found that long-wave roughness can stabilize TS disturbance. Goodhand et al.²⁴ investigated the impact of streamwise geometric variation near the leading edge and found that even small-scale variation can modify the separation and reattachment locations drastically and, thus, affect the aerodynamic performance of the airfoil, which is in accordance with the Stebbins et al.,¹⁸ i.e., even smallest roughness elements could cause noticeable reductions in aerodynamic capabilities.

Despite many studies on this topic, few have attempted to describe the flow in detail, explaining the mechanisms responsible for the changes in the aerodynamic forces caused by the roughness elements. In the present research, we are focusing on the effect of surface roughness elements on the aerodynamic performance of low Reynolds number airfoil as well as the underlying mechanism. To that end, we perform high-resolution DNS of the NACA0012 airfoil near stall condition, with streamwise wavy roughness elements located near the leading edge. The length of the roughened region and the height of the roughness elements are fixed while varying the wavelength. The Reynolds number is set at $Rec=5×104$⁠, providing a good balance between the tractability of the underlying mechanisms and the computational cost of the highly resolved simulations.

The paper is organized as follows. The computational details, including the flow configuration, numerical scheme, and boundary conditions, are given in Sec. II. The DNS code is validated in Sec. III. The results are presented and analyzed in Sec. IV, including the lift and drag coefficients, pressure and skin friction coefficients, mean flow field, unsteady separation bubble development, and scaling of wake dynamics. Finally, some conclusions are drawn in Sec. V.

The artificial wavy roughness elements are placed near the leading edge on the suction side, as shown in Fig. 1. The tuple $(s0,yn)$ denote the wall-tangential and wall-normal coordinates to the smooth airfoil surface, with the origin located at $x/C=0.037$ and the arc length from the leading edge to the origin is $0.05Lu$⁠, where L_u is the length of the upper airfoil surface. These roughness elements are characterized by the sinusoidal-wave structure as follows:

where h and λ are the roughness elements' amplitude (height) and wavelength, respectively. The roughened region starts from $x/C=0.037$ and ends at $x/C=0.136$⁠, occupying 10% of the upper surface length (i.e., $0.1Lu$⁠). The roughness height is set to $h=2×10−3C$⁠, and the roughness wave number (RWN) is $k=0.1Lu/λ$⁠.

The location of the roughened region follows the experimental work of Fujino et al.²⁵ which is also adopted by Bagade et al.²⁶ This choice is prompted by the fact that the most common surface roughness, such as ice accretion, occurs in this region, which is also susceptible to many important flow phenomena, such as dynamic stall. This small-scale roughness height is larger than Fujino et al.²⁵ but slightly smaller than the cylindrical element (⁠$h=3×10−3C$⁠) as adopted by Kumar et al.²⁰ The angle of attack (AoA) and Reynolds number are fixed, i.e., $AoA=10°$ and $Rec=5×104$⁠. This Reynolds number is smaller than the realistic (high speed) flight (⁠$Rec∼106$⁠); however, it is applicable to the small-scale autonomous aircraft, wind turbines, etc.²⁷ The angle of attack considered here is close to the stall condition, where the aerodynamic performance can be strongly modified by roughness elements. We perform DNS of six configurations, i.e., $k=0, 4, 6, 8, 10, and 12$⁠, to investigate the effects of different roughness geometries, in which k = 0 corresponds to the smooth case.

The incompressible Navier–Stokes equations in the generalized curvilinear coordinates read

where $(ξ1,ξ2,ξ3)=(ξ,η,ζ)$ denote the generalized curvilinear coordinates, and U^m and $Fim$ are given by

The triplet $(x1,x2,x3)=(x,y,z)$ denote the Cartesian coordinates with corresponding velocity components denoted as $(u1,u2,u3)=(u,v,w)$⁠. The symbol p indicates the pressure, ν is the kinematic viscosity, $J−1$ is the Jacobian of the transformation, U^m is the volume flux normal to the surface of constant ξ^m, and G^mn is the mesh skewness tensor. It should be noted that the $ξ3$ direction is congruent with the z direction since the spanwise geometry is uniform in the present research.

The computational domain is discretized in hexahedral elements using a C-type collocated body-fitted mesh. The governing equations given by (2) are spatially discretized using the energy-conservative fourth-order finite difference scheme proposed by Morinishi et al.²⁸ Specifically, the convective term is discretized in the skew-symmetric form to minimize the aliasing error.²⁹ The pressure–velocity coupling is achieved by the semi-implicit fractional step method,^30,31 which follows the predictor-corrector procedure. The solution is advanced in time by the Adams–Bashforth method for the explicit terms and the Crank–Nicolson scheme for the implicit terms. The incompressibility condition is enforced by solving a Poisson-type equation for the pressure, p. Both the modified Helmholtz equations resulting from the implicit treatment of viscous terms in the predictor step and the pressure Poisson equation are solved using the multigrid solver with line-relaxed Gauss–Seidel smoothers. The DNS code has been parallelized using the standard MPI protocol. To achieve optimal load balancing, the meshes are divided into blocks of equal size, and each of them is assigned to a unique processor. The in-house developed DNS code is described in detail by Zhang et al.³² for DNS of airfoil flows.

The numerical setup and domain size are illustrated in Fig. 2. The computational domain size should be large enough to minimize the effects of the artificially imposed boundary conditions. In the present simulations, the domain size (wake length $W=10C$ and domain radius $R=6C$⁠) is identical to the choice of Zhang and Samtaney³³ for DNS of flow past the NACA0018 airfoil at $Rec=5×104$⁠. We note that a sufficiently long spanwise domain size (L_z) is important to resolve the largest turbulent structures in the spanwise direction. Here, we use $Lz=0.8C$⁠, as recommended by Zhang and Samtaney,³³ which is larger than most of the DNS and LES simulations of an airfoil flow with similar Re_c.^34–37 In addition, we impose a uniform flow $(u,v,w)=(U0,0,0), U0=1$ at the inlet, the convective boundary condition $∂u/∂t+UB∂u/∂x=0$ on the outflow plane, where U_B is the bulk mean velocity, and a no-slip wall boundary condition on the airfoil surface. Finally, we use periodic boundary condition in the spanwise direction.

All the simulations are performed with a mesh $Nξ×Nη×Nz=2048×256×256$⁠. The time step is restricted by the advective Courant–Friedrichs–Lewy (CFL) condition. Here, the CFL number is set to $CFL=0.8$⁠. The mesh is clustered around the airfoil surface, and the maximum wall unit $Δη+$ is around 0.9 (based on the local shear velocity and first wall-normal grid size) to ensure that the near-wall region is fully resolved. The maximum wall units in ξ and z directions, i.e., $Δξ+$ and $Δz+$ are less than 10 in the turbulent region. The mesh resolution has been carefully verified in Refs. 33 and 38 and follows the well accepted criterion for DNS of wall-bounded turbulent flows.³⁹ The total simulation time is approximately $60C/U0$⁠, and only the last $45C/U0$ time units are collected for statistical analysis. All the simulations are performed on the Cray XC40 Shaheen II supercomputer at KAUST. The total computational cost is around 15 million core hours.

This section presents several statistically averaged quantities to validate the DNS code. All these quantities are averaged in time and spanwise direction. The latter average procedure is used throughout the paper unless stated otherwise.

Figure 3 shows the distribution of time- and spanwise-averaged pressure coefficient C_p for the smooth airfoil case (k = 0). Here, C_p is computed as

where p_w is the time- and spanwise-averaged wall pressure, and $p∞$ is the reference pressure taken at the outlet boundary. It can be found that the DNS results compare well with the experimental results from Yamaguchi et al.⁴⁰ However, there is some discrepancy near the pressure plateau, which is usually used to identify the transition onset and separation point in experiments.⁴¹ 

Figure 4 shows the mean and root mean square (RMS) velocity profiles along the suction side of the smooth airfoil. The DNS predicted mean velocity profiles (⁠$u¯/U0$⁠) are in good agreement with the experimental results from Ohtake,⁴² and the reverse flow near the leading edge, i.e., $x/C=0.1$ is well captured. For the RMS velocity profiles (⁠$u′¯/U0$⁠), the near-wall gradient of $u′¯/U0$ of the present DNS is slightly smaller than the experimental results within the region $x/C∈[0.4,0.9]$ but compare well in the outer region. Overall, the DNS predicted results generally agree with the experimental results.

In this section, we investigate the roughness effect on the aerodynamic performance of the airfoil, (i.e., the lift and drag coefficients) and the flow field near the body.

The sectional lift and drag coefficients are computed as

where the lift and drag forces are computed by integrating the wall pressure and shear stress over the airfoil surface.

Figure 5 shows the mean and RMS values of the lift and drag coefficients for different roughness wave numbers. In what follows, $ΔCL$ denotes the decrement of the lift coefficient and the superscript indicates the roughness wave number, k. The lift coefficient, C_L, is decreasing from k = 0 to k = 8, with a smooth variation from k = 0 to k = 6 (⁠$ΔCL/CL0=1%$⁠) but a sharp decrease from k = 6 to k = 8 (⁠$ΔCL/CL6=9.8%$⁠). At k = 8, C_L reaches its minimum value, beyond which we observe only small variations. The correlation of C_L with k is similar to the drag crisis phenomenon, as observed in bluff body flows, in which the drag coefficient falls suddenly as the Reynolds number increases beyond the sub-critical regime.^43–45 The drag crisis phenomenon is attributed to either transition or dynamics of separation bubbles.⁴⁵ We will further discuss the behavior of C_L in Secs. IV B and IV C.

The drag coefficient, C_D, shows a monotonic correlation with k, which indicates that a larger roughness wave number makes the airfoil's surface rougher and increases the drag. Furthermore, the variation of C_D with k is much stronger than that of the lift coefficient, C_L. In fact, from k = 0 to k = 6 we observe that $ΔCD/CD0=21.1%$⁠. Finally, a sharp increase from k = 6 to k = 8 is observed (⁠$ΔCD/CD6=85.1%$⁠), similar to the lift coefficient.

The RMS of the lift and drag coefficients, $CL,rms$ and $CD,rms$⁠, show similar trend with k. For $k≤8$⁠, they grow with k with much stronger variations, i.e., $CL,rms6/CL,rms0=4.3$ and $CD,rms6/CD,rms0=11.7$⁠. A sharp increase from k = 6 to k = 8 is also observed for both $CL,rms$ and $CD,rms$⁠, i.e., $CL,rms8/CL,rms6=4.1$ and $CD,rms8/CD,rms6=2.1$⁠. For $k≥8$⁠, the variations in the mean and RMS of lift and drag coefficients are not large.

Figure 6 shows the contribution ratio of the pressure to the lift and drag coefficients, for different roughness wave numbers. We note that the pressure contributes more than $99%$ to the lift coefficient while the contribution ratio of pressure to the drag coefficient is strongly increasing with k, from $82.9%$ at k = 0 to $96%$ at k = 12, indicating that the form drag is dominant for all k cases. Therefore, it is important to examine the distributions of pressure and skin friction coefficients, since the lift and drag coefficients are based on the integration of them around the airfoil. We emphasize that although the contribution of the skin friction coefficient is small (see Fig. 6), the skin friction behavior remains quite important since it is the principal indication of near-wall dynamics including near-wall separation and reattachment.

Based on the above-mentioned analysis, we plot the distributions of the pressure and skin friction coefficients in Fig. 7. The pressure coefficient, C_p, is computed from Eq. (4) and the skin friction coefficient, C_f, is calculated as

where $τw,s$ is the time- and spanwise-averaged streamwise (parallel to the airfoil surface) wall shear stress, computed from the wall-normal differentiation using the third-order accuracy one-sided velocity-derivative method.⁴⁶ 

From the C_p profile shown in Fig. 7(a), we observe that the magnitude of the pressure coefficient decreases with increasing roughness wavenumber, k, of the roughened region. The sudden change in C_p from k = 6 to k = 8 is in accordance with the lift and drag coefficient behaviors as shown in Fig. 5. Downstream the roughened region, the recovery of C_p for $k≤6$ is much faster and stronger than $k≥8$⁠, which produces larger adverse pressure gradients for $k≤6$ except for the zone near the trailing edge, i.e., $x/C>0.9$⁠. For $x/C>0.3$⁠, the C_p-profiles for $k≤6$ and $k≥8$ cases seem collapse on two different curves, respectively. Moreover, the C_p profile from k = 8 is slightly larger than k = 10 and 12. Not surprisingly, these differences in the coefficient are not obvious on the pressure side, except for $x/C>0.9$⁠, which would affect the wake development downstream of the trailing edge (see Sec. IV D for more discussions).

The locations of separation and reattachment can be detected from the plot of the friction coefficient, C_f, as shown in Fig. 7(b). We plot the time- and spanwise-averaged streamlines to further visualize the separation bubbles in Fig. 8. For the smooth airfoil case (i.e., k = 0), a small separation bubble is observed near the leading edge, also visible in Fig. 8(a) and in agreement with observations by Eljack et al.⁴⁷ (compressible DNS with Mach number of 0.4). With the increase in the roughness wave number (⁠$k≤6$⁠), the reattachment location is moving downstream, and the primary separation bubble is growing (see Fig. 8). More secondary separation bubbles are formed around the roughened region with increasing k. From k = 6 to k = 8, the primary separation bubble is suddenly enlarged, and massive separation occurs on the suction side beyond k = 8, covering approximately $99%$ of the suction surface. For k = 8, there are two bubble centers detected from the closed streamlines, beyond which the bubbles are relaxed, and one primary bubble is formed for k = 10 and 12, with the bubble center moving downstream. The separation strength can be measured by the magnitude of C_f inside the separation zone. We can find that although the size of the primary separation bubble is much larger for $k≥8$⁠, the strength of separation is weaker than $k≤6$⁠. In the separation zone, the peak magnitude of C_f for $k≤6$ is decreasing with k, and “oscillating” C_f-profiles are observed in rough airfoil cases [see Fig. 7(b)], which indicate formations of secondary bubbles inside the wavy grooves. The undershoots of local C_f after the roughened region (⁠$x/C>0.24$⁠) are enhanced with increasing k. This behavior is consistent with the observations in the rough-to-smooth transition of turbulent boundary layer.⁴⁸ The memory effect of the upstream roughness condition is minimized for $x/C≥0.5$ with attached boundary layer (⁠$k≤6$⁠). Still, it does not show much variation for the detached boundary layer (⁠$k≥8$⁠) due to massive separation.

The sudden growth of the primary separation bubble from k = 6 to k = 8 may also explain the variations of lift and drag coefficients with roughness wave numbers greater than zero, as shown in Fig. 5. These separation and reattachment behaviors imply strong effects on aerodynamic performance. However, we would like to check if the transition also plays a key role. The contour plot of the time- and spanwise-averaged off diagonal Reynolds stress, $−u′v′¯/U02$ are also shown in Fig. 8. Here, the threshold value of 0.001 is chosen as an indicator for transition onset, which is commonly used in various numerical⁴⁹ and experimental investigations⁵⁰ of flow past bluff bodies. The locations of transition onset do not show obvious differences, occurring at $x/C≈0.09 ± 0.02$⁠, on the top of primary separation bubbles for all k cases (see Fig. 8). This feature indicates that the aerodynamic performance (lift and drag coefficients) is dominated by these large-scale separation bubbles rather than transition with different roughness configures.

Figure 9 shows RMS of wall pressure (⁠$pw,rms$⁠) and wall shear stress (⁠$τw,rms$⁠) along the suction side of the airfoil for different roughness wave numbers. For $k≤6$⁠, both $pw,rms$ and $τw,rms$ are increasing inside the roughened region, and decreasing rapidly to some stable level after the reattachment location. For $k≥8$⁠, the k = 8 case produces slightly larger $pw,rms$ and $τw,rms$ for $x/C≤0.85$⁠. Meanwhile, $τw,rms$ is an order of magnitude smaller than $pw,rms$⁠, indicating that the contribution of wall pressure fluctuations to the RMS of lift and drag coefficients is much larger than wall shear stress fluctuations. The latter is same as the contribution of wall pressure to lift and drag coefficients, as shown in Fig. 6. The mean level of $pw,rms$ on the suction side for $k≥8$ is larger than smaller k cases and produces the sharp increase in $CL,rms$ and $CD,rms$ from k = 6 to 8 as shown in Fig. 5.

In Sec. IV B, the mean skin friction coefficient and streamline from the DNS have been discussed, which provide helpful information for interpreting the effects of the separation and reattachment on the aerodynamic performance (in the time-averaged sense) due to the wavy roughness. However, these mean-flow concepts are insufficient to fully understand the near-wall flow physics, including the unsteady load force and dynamics of separation bubbles that important in aerodynamic flows. In the present section, we discuss the unsteady separation and reattachment, first on the convective motion of the separation bubble, and then examine the proportion of reverse flow duration in the near-wall field.

Figure 10 shows the space-time diagram of the instantaneous spanwise-averaged skin friction coefficient for different roughness wave numbers, which also represents the evolution of the multiple vortical structures along the suction side.³³ The evolution of separation and reattachment, although shown only for a short duration of $5C/U0$⁠, reveals the primary features. It is evident that the location of the first separation point remains virtually unchanged for the entire time sequence, the same as the mean separation point [see Fig. 7(b)]. Meanwhile, unsteady separation and reattachment can be observed near the trailing edge, which are missing in the mean C_f-profiles. The flow reversal is strongest around $x/C=0.15$⁠, 0.16, and 0.2 for k = 0, 4, and 6, respectively, and moves downward the airfoil with increasing k. This is also consistent with the locations of mean negative C_f peak as shown in Fig. 7(b). The local convective motions (around the right boundary of the roughened region) of separation bubbles can be identified by the convective speed V_c, based on the slope of the stripes. This velocity corresponds to the signature of shear layer structures on the wall.⁵¹ It is found that the convective velocities for different roughness wave numbers are $Vc=(0.6±0.06)U0$⁠, which is similar to the convective velocity of the reattaching shear layer in backward-facing step flows (⁠$Vc=0.6U0$⁠).⁵² For the smooth airfoil case (k = 0), $Vc=0.66U0$ is larger than the NACA0012 case with the same Re_c but $AoA=5°$⁠, i.e., $Vc=0.56U0$ as reported by Zhang and Samtaney,³³ which confirms the argument that larger AoA can produce larger V_c (see Ref. 53). Dandois et al.⁵¹ summarized the convective velocities in flows of different configurations (mixing layer, backward-facing step, mixing layer past a cavity, fence, flat-plate leading edge) and found a clear convective velocity, $Vc∈[0.4,0.65]U0$ in these flows. These convective velocities are weakly correlated with the lift coefficients for different roughness wave numbers (see Fig. 5), but the sharp decrease in V_c from k = 6 to 8 corresponds to the sharp variation of C_L and C_D. Meanwhile, the convective velocity is increasing to $Vc=0.62U0$ at k > 8, which is slightly larger than V_c of k = 4 and 6 (⁠$Vc≈0.6U0$⁠). This acceleration may be attributed to the interactions of secondary bubbles (merging and breakup of bubbles with decreased $|Cf|$⁠) in the roughened region. For $k≥8$ cases, although some clear convective velocities can be observed beyond the right boundary of the roughened region, the convection motions of the shear layer structures are complex due to massive separation, which is linked to the convective phenomenon of different frequencies.⁵¹ The interactions of instantaneous separation bubbles are also reported by Gao et al.³¹ for the NACA0018 case (⁠$AoA=5°, Rec=105$⁠), which may also contribute to the wavy roughness effects on the aerodynamic performance. The convective separation bubbles downward the trailing edge will also change the evolution of wake flows, which can explain the roughness effects based on scaling analysis (see Sec. IV D for further discussions).

Figure 11 shows the ratio of reversed flow for all k cases (time of near-wall reversed flow $(Cf(x,z,t)<0)$ over the total simulation time). The flow is fully attached to the surface upstream of the leading edge detachment point, which is the same as the instantaneous spanwise-averaged skin friction coefficient as shown in Fig. 10. However, the unsteady behavior of the separation and reattachment beyond that detachment point is absent in the mean C_f-profiles. The roughened region is alternately occupied by separation and reattachment, which is obvious near the last mean reattachment lines. For $k≤6$⁠, the mean attached zones are nibbled, especially near the trailing edge. For the mean attached zone of k = 6 case, the number of reverse flow events is small compared to the total number of events. However, the gentle nibble grows into massive annexation, and most of the suction side is dominated by separation for $k≥8$⁠. The spanwise variations of the reverse fraction represent the meandering of the wall streamlines, which consist flow structures of different scales and frequencies. All these unsteady motions will contribute to larger RMS of lift and drag coefficients with growing k.

The wake downstream the trailing edge of an airfoil is asymmetric due to the effects of pressure gradients and curvature.^54,55 Based on the momentum conservation law, the drag force acting on an airfoil is related to the wake structure.^56,57 Presently, we try to explore the drag mechanism behind the wavy roughness effects from the perspective of wake evolution, which is strongly affected by different separation behaviors (see Secs. IV B and IV C) on the suction side of the airfoil with partially covered roughness elements.

The first quantity of interest is the mean velocity profile in the wake region, where the classical self-similarity wake profile fails for the asymmetric wake.^55,58,59 Following the scaling approach proposed by Hah and Lakshminarayana,⁵⁹ we define two characteristic length scales l_s, l_p and two characteristic velocity, U_c and U_d, where U_c is the lowest wake velocity in the x direction; $Ud=U0−Uc$ is the maximum velocity defect; l_p and l_s denote the distance from the wake center (location of U_c) to the points where the velocity defect, $u¯−Uc$⁠, is equal to $Ud/2$ on the pressure and suction sides, respectively (see Fig. 12). Hah and Lakshminarayana⁵⁹ proposed the following Gauss profile for the wake defect velocity:

where y_c is the coordinate of the wake center, and Λ is a constant equal to –0.6993.⁵⁹ 

Figure 13 shows the scaled wake profiles for different roughness wave numbers at two typical locations, i.e., $x/C=1$ close to the trailing edge and $x/C=2$ away from the trailing edge. In the very-near-wake region, i.e., $x/C=1$⁠, it is evident that the mean velocity profiles on the suction side follow the scaling equation (7) well but deviate sharply on the pressure side. For $k≥8$⁠, the deviation is even pronounced, which might be caused by the stronger pressure gradient effects at the trailing edge (see Fig. 7). The separation bubbles close to the trailing edge for $k≥8$ (see Fig. 8) have substantial effects on the evolution of the wake and make the wake profiles more asymmetric. For $x/C=2$⁠, the “history effect” from the trailing edge is vanishing, the wake structure is becoming more symmetric, and the wake profiles follow the scaling equation (7) well.

Based on the above-mentioned scaling approach, the maximum velocity defect, $(U0−Uc)/U0$ and wake width, $lw=ls+lp$ for different roughness wave numbers are shown in Fig. 14. It is evident that the wake maximum velocity defect at the wake center is decaying with increasing downstream distances. In the very-near-wake region, i.e., $x/C=1$⁠, the defect is sharp and increasing with k, but this trend disappears for $x/C≥1.1$⁠. For $x/C≥1.2$⁠, the defect is decreasing with k, which means that the velocity recovery at the wake center is faster for smaller roughness wave numbers. The decay rate of velocity defect is increasing with k for $x/C≤1.4$⁠, but almost the same for $x/C≥1.5$⁠. The wake width is increasing with k, and a sharp increase can be observed from k = 6 to k = 8, which is in accordance with the drag coefficients as shown in Fig. 5. The growth rate of the wake width is also increasing with k and much faster for $k≥8$⁠, which confirms that k = 8 is a critical case.

Qualitatively, Pope⁶⁰ proposed that the momentum deficit flow rate, which equals the drag on the bluff body, follows:

It should be noted that although the wake defect is decreasing with k, the variation is much weaker than the wake width, which to some extent, explains the drag mechanism with different roughness wave numbers, as shown in Fig. 5(b).

We perform a direct numerical simulation of flow past the NACA0012 airfoil at 10° angle of attack and Reynolds number of $5×104$ with partially covered wavy roughness elements to investigate the roughness effect on the aerodynamic performance and the underlying mechanism. The simulations are performed with an energy-conservative fourth-order parallel code solving the incompressible Navier–Stokes equations in generalized curvilinear coordinates. The length of the roughened domain and roughness height are fixed and located near the leading edge, whereas the roughness wavenumber k varies from 4 to 12. The smooth airfoil is used as a baseline case.

The roughness elements degrade the aerodynamic performance, i.e., reduce the lift coefficient and increase the drag coefficient. The drag coefficient increases monotonically with k, while the variation of lift coefficient with k is similar to the drag crisis phenomenon as observed in cylinder flows. Both the sharp variations of lift and drag coefficients from k = 6 to k = 8 indicate that k = 8 is a critical case. The latter result is also confirmed in the RMS of the lift and drag coefficients.

To reveal the underlying mechanism, we examine the mean pressure, skin friction coefficient, and the mean flow field. The roughness elements can strongly affect the separation behavior but have little effect on the location of transition onset. For $0≤k≤6$⁠, the main separation bubble stretches around the roughened region, and the secondary separation bubble occurs at the windward side of the roughness elements. The case of k = 8 is a critical case. In fact, beyond k = 8, a massive separation occurs on the upper surface, drastically affecting the aerodynamic performance.

The evolution of separation and reattachment on the airfoil surface is analyzed to further explore the unsteady load force with the dynamics of separation bubbles. Around the roughened region, clear convection motions of separation bubbles can be observed. The convection velocity is largest for the smooth airfoil, and the roughness elements retard the convection motions, which is strongest for the critical case (k = 8). These convective bubbles also affect the wake evolution downstream of the trailing edge. The asymmetric wake is analyzed by employing the scaling law by Hah and Lakshminarayana.⁵⁹ The wake profiles on the pressure side do not follow in the very-near-wake region for $k≥8$ due to the downward convective bubbles, but history effect is vanishing far away from the trailing edge, and wake profiles follow the scaling law well. The wake defect is decreasing with k, but the variation is much weaker than the wake width, showing a sharp increase from k = 6 to 8. The wake scaling analysis, to some extent, confirms the drag mechanism induced by the roughness elements.

The present simulations are performed at a specific Reynolds number and angle of attack, and it is possible that some of the underlying mechanism changes if the Reynolds number is increased, or the angle of attack is varied. It would be desirable to verify this conjecture by performing simulations at higher Reynolds number and different angles of attack, which can be potentially done in the future.

The Cray XC40 Shaheen II at KAUST was used for all simulations reported. This research was partially supported under KAUST OCRF URF/1/1394-01 and under baseline research funds of R. Samtaney. W. Gao acknowledges the financial support by Visiting Researcher Fund Program of State Key Laboratory of Water Resources and Hydropower Engineering Science (Grant No. 2020HLG02).

The authors have no conflicts to disclose.

Wei Gao: Conceptualization (equal); Formal analysis (equal); Investigation (lead); Methodology (equal); Validation (lead); Writing – original draft (lead); Writing – review & editing (equal). Ravi Samtaney: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Resources (lead); Supervision (lead); Writing – review & editing (equal). Matteo Parsani: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.