Direct numerical simulation of transitional ﬂow past an airfoil with partially covered wavy roughness elements

freestream airfoil length ( C ), of attack ( AoA are ﬁxed, i.e., Re c = 5 × 10 4 , AoA = 10 ◦ The leading edge roughness elements by streamwise sinusoidal-wave geometry that is uniform in the spanwise direction with ﬁxed varying wave numbers ( k ) from 0 to 12. validation of the smooth baseline case ( k = 0), roughness effects drag coefﬁcients. coefﬁcient monotonically coefﬁcient


Accepted to Phys. Fluids 10.1063/5.0123437
We perform direct numerical simulation (DNS) of flow past the NACA (National Advisory Committee for Aeronautics) 0012 airfoil with partially covered wavy roughness elements near the leading edge.The Reynolds number Re c based on the freestream velocity (U 0 ) and airfoil chord length (C), and the angle of attack (AoA) are fixed, i.e., Re c = 5 × 10 4 , AoA = 10 • .The leading edge roughness elements are characterized by streamwise sinusoidal-wave geometry that is uniform in the spanwise direction with fixed height whereas varying wave numbers (k) from 0 to 12. Based on validation of the smooth baseline case (k = 0), the roughness effects on the aerodynamic performance are evaluated in terms of the lift and drag coefficients.The drag coefficient decreases monotonically with k, while the variation of the lift coefficient with k is similar to the "drag crisis" phenomenon observed in cylinder flows.The sharp variations of lift and drag coefficients from k = 6 to 8 indicate that k = 8 is a critical case, beyond which massive separation occurs and almost covers the airfoil's suction side and dominates the airfoil aerodynamic performance.
To further reveal the underlying mechanism, the flow structures, pressure, skin friction coefficients, shear layer transition onset, and boundary layer separation are quantified to investigate the roughness effects.The roughness elements strongly modify the separation behavior, whereas they have little effect on the transition onset.The unsteady interactions and convections of separation bubbles downward the trailing edge also change the wake evolution.Based on the scaling of the asymmetric wake profiles, we find that the wake defect is gently decreasing with k, but the increase of the wake width is much stronger.This scaling analysis gains a better insight into the roughness effects on the momentum deficit flow rate, which confirms the drag mechanism with different roughness wave numbers.a) wei.gao@kaust.edu.sa

I. INTRODUCTION
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. 3 .The pioneering work of quantifying the roughness effect in pipe flows was led by Nikuradse 4 , and Moody 5 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, and 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 10 , and yielding significant departures from Nikuradse's initial proposal.Cheng, Pullin, and Samtaney 11 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. 12 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 13 , irregular ice 14 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 Accepted to Phys.Fluids 10.1063/5.0123437portion 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 15 .Based on representative ice geometries, Bragg, Broeren, and Blumenthal 16 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 17 .However, there are minimal numerical investigations of the impact of accretions on airfoils (refer to the review paper by Stebbins et al. 18 ).
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 20 .
Serson, Meneghini, and Sherwin 21 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, Piomelli, and Lehmkuhl 20 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, Haeri, and Joseph 22 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, Mughal, and Ash-Accepted to Phys.Fluids 10.1063/5.0123437worth 23 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, Miller, and Lung 24 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. 18 , 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 Re c = 5 × 10 4 , 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 § II.The DNS code is validated in § III.
The results are presented and analyzed in § 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 § V.

II. PHYSICAL AND NUMERICAL SETUP A. Physical setup
The artificial wavy roughness elements are placed near the leading edge on the suction side, as shown in Fig. 1.The tuple (s 0 , y n ) 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.05L u , 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 ten percent of the upper surface length (i.e., 0.1L u ).The roughness height is set to h = 2 × 10 −3 C, and the roughness wave number (RWN) is k = 0.1L u /λ .
The location of the roughened region follows the experimental work of Fujino, Yoshizaki, and Kawamura 25 , which is also adopted by Bagade, Bhumkar, and Sengupta 26 .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, Yoshizaki, and Kawamura 25 , but slightly smaller than the cylindrical element (h = 3 × 10 −3 C) as adopted by Kumar, Piomelli, and Lehmkuhl 20 .The angle of attack (AoA) and Reynolds number are fixed, i.e., AoA = 10 • and Re c = 5 × 10 4 .This Reynolds number is smaller than the realistic (high speed) flight (Re c ∼ 10 6 ), however, it is applicable to the small-scale autonomous aircraft, wind turbines, etc 27 .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, 12, to investigate the effects of different roughness geometries, in which k = 0 corresponds to the smooth case.

B. Numerical details
The incompressible Navier-Stokes equations in the generalized curvilinear coordinates read where (ξ 1 , ξ 2 , ξ 3 ) = (ξ , η, ζ ) denote the generalized curvilinear coordinates, and U m and F m i are given by The triplet (x 1 , x 2 , x 3 ) = (x, y, z) denotes the Cartesian coordinates with corresponding velocity components denoted as (u 1 , u 2 , u 3 ) = (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 energyconservative fourth-order finite difference scheme proposed by Morinishi et al. 28 .Specifically, the convective term is discretized in the skew-symmetric form to minimize the aliasing error 29 .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. 32 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 33 2. Sketch of the numerical setup and computational domain.
All the simulations are performed with a mesh N ξ ×N η ×N z = 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 units ∆η + 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 39 .The total simulation time is approximately 60C/U 0 , and only the last 45C/U 0 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.

III. CODE VALIDATION
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.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, Ohtake, and Muramatsu 40 .However, there is some discrepancy near the pressure plateau, which is usually used to identify the transition onset and separation point in experiments 41 .
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 0 ) are in good agreement with the experimental results from Ohtake 42 , and the reverse flow near the leading edge, i.e., x/C = 0.1 is well captured.For the RMS velocity profiles (u /U 0 ), the near-wall gradient of u /U 0 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.

IV. NUMERICAL RESULTS AND DISCUSSIONS
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.

A. Lift and drag coefficients
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.Reynolds number increases beyond the sub-critical regime [43][44][45] .The "drag crisis" phenomenon is attributed to either transition or dynamics of separation bubbles 45 .We will further discuss the behavior of C L in the following sections.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 behaviour remains quite important since it is the principal indication of near-wall dynamics including near-wall separation and reattachment.

B. Mean flow field
Based on the above 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 onesided velocity-derivative method 46 .
From the C p profile shown in Fig. 7(a), we observe that the magnitude of the pressure coeffi- the pressure side, except for x/C > 0.9, which would affect the wake development downstream of the trailing edge (see § 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. 47 (compressible DNS with Mach number of 0.4).With the increase of 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. 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 /U 2 0 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 49 and experimental investigations 50 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 (p w,rms ) and wall shear stress (τ w,rms ) along the suction side of the airfoil for different roughness wave numbers.For k ≤ 6, both p w,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 p w,rms and τ w,rms for x/C ≤ 0.85.Meanwhile, τ w,rms is an order of magnitude smaller than p w,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 the same as the contribution of wall pressure to lift and drag coefficients, as shown in Fig. 6.The mean level of p w,rms on the suction side for k ≥ 8 is larger than smaller k cases, and produces the sharp increase of C L,rms and C D,rms from k = 6 to 8 as shown in Fig. 5.

C. Unsteady separation and reattachment
In the previous section, 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.tiple vortical structures along the suction side. 33The evolution of separation and reattachment, although shown only for a short duration of 5C/U 0 , 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 is 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 51 .It is found that the convective velocities for different roughness wave numbers are V c = (0.6 ± 0.06)U 0 , which is similar to the convective velocity of the reattaching shear layer in backward-facing step flows (V c = 0.6U 0 ) 52 .For the smooth airfoil case (k = 0), V c = 0.66U 0 is larger than the NACA0012 case with the same Re c but AoA = 5 • , i.e., V c = 0.56U 0 as reported by Zhang and Samtaney 33 , which confirms the argument that larger AoA can produce larger V c Accepted to Phys.Fluids 10.1063/5.0123437 to the interactions of secondary bubbles (merging and break-up of bubbles with decreased |C f |) 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 51 .The interactions of instantaneous separation bubbles are also reported by Gao et al. 31 for the NACA0018 case (AoA = 5 • , Re c = 10 5 ), 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 § IV D for further discussions).
Figure 11 shows the ratio of reversed flow for all k cases (time of near-wall reversed flow (C f (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.

D. Wake development
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 § IV B) 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 Accepted to Phys.Fluids 10.1063/5.0123437self-similarity wake profile fails for the asymmetric wake 55,58,59 .Following the scaling approach proposed by Hah and Lakshminarayana 59 , 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; U d = U 0 −U c 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 c , is equal to U d /2 on the pressure and suction sides, respectively (see Fig. 12).Hah and Lakshminarayana 59 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 59 .
12. Sketch of the mean velocity profile in the wake region.
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 scaling approach, the maximum velocity defect, (U 0 − U c )/U 0 and wake width, l w = l s + l p 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, Accepted to Phys.Fluids 10.1063/5.0123437 13.Scaling of the wake defect velocity at (a) x/C = 1 and (b) x/C = 2 for different roughness wave numbers, k.From bottom to top, k is increasing from 0 to 12, shifted by 1, 2, ..., 5 along the ordinate; ---, scaled wake profile from Eq. (7).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 60 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).

V. CONCLUSIONS
We perform a direct numerical simulation of flow past the NACA0012 airfoil at 10 degrees angle of attack and Reynolds number 5 × 10 4 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 is 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, Accepted to Phys.Fluids 10.1063/5.0123437tions, one cannot become an eternal learner.Without becoming an eternal learner, we cannot possibly learn even a fraction of the ocean of knowledge.Ravi, you were an eternal learner and an inspiration to all of us.You left a huge gap, but you also gifted us with things we must carry on.
And we will do so.

Accepted to
for DNS of flow past the NACA0018 airfoil at Re c = 5 × 10 4 .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 L z = 0.8C, as recommended by Zhang and Samtaney 33 , 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) = (U 0 , 0, 0), U 0 = 1 at Accepted to Phys.Fluids 10.1063/5.0123437 the inlet, the convective boundary condition ∂ u/∂t + U B ∂ 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.

Figure 3 - 1 FIG. 3 .
Figure3shows the distribution of time-and spanwise-averaged pressure coefficient C p for the

Figure 5 FIG. 5 .
Figure 5 shows the mean and RMS values of the lift and drag coefficients for different roughness wave numbers.In what follows, ∆C L 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 (∆C L /C 0 L = 1%) but a sharp decrease from k = 6 to k = 8 (∆C L /C 6 L = 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 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 ∆C D /C 0 D = 21.1%.Finally, a sharp increase from k = 6 to k = 8 is observed (∆C D /C 6 D = 85.1%), similar to the lift coefficient.The RMS of the lift and drag coefficients, C L,rms and C D,rms , show similar trend with k.For k ≤ 8, they grow with k with much stronger variations, i.e., C 6 L,rms /C 0 L,rms = 4.3 and C 6 D,rms /C 0 D,rms = 11.7.A sharp increase from k = 6 to k = 8 is also observed for both C L,rms and C D,rms , i.e., C 8 L,rms /C 6 L,rms = 4.1 and C 8 D,rms /C 6 D,rms = 2.1.For k ≥ 8, the variations in the mean and RMS of lift and drag coefficients are not large.

3 FIG. 7 .
FIG. 7. (a) Distribution of the pressure coefficient, C p around the airfoil, and (b) the skin friction coefficient, C f , on the suction surface for different roughness wave numbers, k.The shaded region denotes the roughened region.The horizontal dashed line C f = 0 in (b) is shown for convenience: zero crossing of this line indicates separation and reattachment.

Figure 10 FIG. 9 .
Figure10shows the space-time diagram of the instantaneous spanwise-averaged skin friction coefficient for different roughness wave numbers, which also represents the evolution of the mul-

FIG. 11 .
FIG. 11.Ratio of the time of near-wall reversed flow (C f (x, z,t) < 0) to the total simulation time for different roughness wave numbers, k.(a) k = 0, (b) k = 4, (c) k = 6, (d) k = 8, (e) k = 10, ( f ) k = 12.The vertical dashed lines denote the mean separation and reattachment lines detected from the mean C f -profiles, as shown in Fig. 7.