Experimentally validated three-dimensional computational aerodynamics of wind turbine blade sections featuring leading edge erosion cavities

Wind turbine blade leading edge erosion reduces the lift and increases the drag of the blade airfoils. This occurrence, in turn, reduces turbine power and energy yield. This study focuses on the aerodynamic analysis of large and sparse erosion cavities, observed in intermediate to advanced erosion stages, whose size and surface pattern do not lend themselves to experimental and numerical analysis by means of distributed roughness models alone. Making use of three-dimensional Navier-Stokes computational fluid dynamics enhanced by laminar-to-turbulent transition modeling, and geometrically resolving individual erosion cavities, the study validates this simulation-based approach for predicting the aerodynamics and performance loss of blade sections featuring the aforementioned erosion cavities against available experimental data. It is found that the considered cavities can trigger transition, indicating the necessity of both resolving their geometry in the simulations and also modeling distributed surface roughness, of typically lower level, as this latter affects the properties of boundary layers and, if sufficiently high, may trigger transition over the entire spanwise length affected. The energy yield loss of a utility-scale turbine due to the considered erosion pattern is found to be between 2.1% and 2.6% using measured and computed force data of the nominal and eroded outboard blade airfoil. A parametric analysis of the cavity geometry suggests that the geometry of the cavity edge has a much larger impact on aerodynamic performance than the cavity depth. show 10A. c the section normal to the cavity axis, and the chordwise section containing the cavity axis. All results highlight the complex 3D character of the cavity flow, and the left and central subplots also point to the lack of flow symmetry discussed above. The j v j contours of the spanwise section (left subplot) also highlight the thicker boundary layer on the right-hand side of the cavity due to the blockage associated with part of the cavity fluid leaving the cavity in this region. It has been found that the cavity flow asymmetry is observed also using symmetry BCs on both lateral boundaries of the domain. In light of the symmetry of the airfoil model and the boundary conditions, the lack of flow symmetry may look surpris-ing. However, small numerical perturbations of the numerical solution, such as those resulting from the CFD grid not being symmetric about the midspan plane, may result in a symmetric flow state not being achievable. In the context of the continuous governing conservation laws, this cir-cumstance may correspond to the existence of a supercritical pitchfork bifurcation, 34 whereby, above a critical Reynolds number and/or AoA, the symmetric solution becomes unstable and the nonsymmetric solution with significant levels of vortical cavity flow in the plane normal to the cavity axis becomes instead stable. In the authors' view, the nonsymmetric vortical solution is more representative of real erosion cavities, featuring irregular shapes which would prevent symmetric solutions even ignoring the desymmetrizing perturbations due to rotational effects and blade twist and taper.

consequently, significant turbine and wind farm energy losses. 2 Large energy losses due to roughness and geometry alterations of the LE surface region can also occur due to insect, 3,4 dust, 5 and ice accretion. 6,7 Alterations of the LE surface geometry due to erosion result in drag increase and lift reduction of the airfoils making up the rotor blades, and both occurrences reduce rotor torque. The loss is exacerbated by the fact that erosion is particularly severe in the outboard blade region, where most energy capture takes place. The consequent wind power loss reduces the turbine annual energy production (AEP) by an amount depending on both the erosion severity and wind characteristics at the installation sites, with some field data and numerous numerical studies pointing to AEP reductions from between 2% and 5% 8-10 up to 10%. 2,11 The assessment of AEP losses due to LE erosion is paramount to wind farm operators, as these energy losses cause a loss of revenue. Thus, the availability of reliable erosion-induced AEP losses would enable wind farm operators to optimize maintenance and repair operations, which are particularly costly in the offshore case. Partly prompted by this industrial requirement, a novel AEP loss estimation system (ALPS) 11 was recently developed. Making use of artificial neural networks and a force database of eroded airfoil with LE delamination generated using highperformance computing-enhanced two dimensional (2D) computational fluid dynamics (CFD), ALPS can return in a few seconds WT AEP losses due to LE erosion using surface state data of operational turbine blades as input.
AEP losses due to LE erosion may be alleviated in several ways. One approach consists of using blade LE protection technologies such as coatings, tape, and shields. 1 This technology has the potential of significantly alleviating the LE erosion problem, but it still faces significant challenges, such as the adhesion of the protective system and the blade material, defects of the protective system due to its application and/or operation, and possible deviations of the blade airfoil shapes from their design intent due to the protective system. A recently proposed alternative approach to mitigating AEP losses due to erosion is to monitor rain falls at the installation site and reduce the rotor speed during the rainfall events that have the highest impact on LE erosion. 12 The approach relies on the frequency of these events being very low and the possibility of preserving the integrity of the LE surface by accepting a relatively modest AEP loss due to using suboptimal rotor speeds. The deterioration of blade section aerodynamic performance due to increased LE surface roughness, which represents the initial stage of LE erosion, has been addressed for a long time by designing WT airfoils whose performance sensitivity to LE roughness is minimal. [13][14][15] One of the key effects of surface roughness is to lead to premature laminar-to-turbulent transition of the airfoil boundary layer (BL). This occurrence depends also on whether the average roughness height k is notably smaller than the BL thickness or is comparable to it. 16,17 One approximate method to assess the impact of surface roughness on transition is to use a Reynolds number Re k based on k and the wall friction velocity u τ and to compare this value to a critical threshold Re k ð Þ cr . 18 The model predicts that the considered roughness level triggers transition if Re k > Re k ð Þ cr .
Many previous research efforts to assess the impact of LE erosion on WT airfoil performance have focused mainly on the detrimental effect of relatively small, often distributed, LE roughness. Furthermore, field observed larger erosion patterns have been rarely used in experimental and numerical investigations. However, there is some evidence that blade maintenance is often performed when the erosion level is well beyond the level of higher distributed surface roughness. This study aims at investigating the aerodynamic losses of WT airfoils at erosion stages between a moderate increase of surface roughness and severe LE delamination.
The performance loss mechanisms vary as erosion progresses. At the early stages, when the surface geometry alterations amount to fractions of a millimeter and surface roughness increases, the laminar-to-turbulent transition moves towards the LE, 17 increasing drag and decreasing lift due to thicker BLs and higher pressure drag. As erosion progresses, however, resulting in the formation of larger and sparser cavities, the phenomena accounting for performance loss become more complex, involving also fluid dynamic interactions on scales larger than that of the BL thickness. On the experimental side, Sareen et al. 2 carried out extensive experimental campaigns to assess the performance loss of the DU 96-W-180 airfoil subject to varying LE erosion levels, from early erosion stage pits, through intermediate stage pits and gouges to advanced LE delamination.
These analyses were conducted by using an airfoil model with LE deliberately damaged with patterns inspired by photographic records of operational multi-MW rotor blades and eroded blades under repair collected over 10+ years. Gaudern 19 performed similar field data-based experimental investigations that differed from that of Sareen et al. 2 primarily for using profiled shape reproducing simplified versions of observed damages rather than damaging directly the airfoil model.
On the computational side, Wang et al. 20 used 2D Reynolds-averaged Navier-Stokes (RANS) CFD to investigate the aerodynamic performance of WT airfoils featuring LE erosion pits. Each cavity was modeled as a semicircular cavity, and parametric analyses of the dependence of the aerodynamics performance loss on the depth, surface density, surface extent, and location of the erosion pits were carried out, determining critical damage levels above which the airfoil performance did not decrease further in a significant way. Campobasso et al. 11 used 2D RANS CFD to investigate the impact of severe LE delamination on both airfoil and WT rotor performance. The analyses were validated against wind tunnel measurements of the DU 96-W-180 airfoil with LE delamination, 2 and a good agreement of numerical results and measured data was obtained.
Han et al. 9 also used 2D RANS CFD to assess the impact of LE erosion and contamination on WT power and energy yield. The severe LE delamination pattern they examined was based on blade surface damage data of a commercial 660-KW turbine. Castorrini et al. 10 developed a scripted and modular CAD functionality and used three-dimensional (3D) RANS CFD to assess the impact on large and sparse LE erosion cavities on the power curve and AEP of the National Renewable Energy (NREL) 5-MW reference turbine, 21 but the study focused on computational geometry and WT overall performance analysis.
This study focuses on the aerodynamic analysis of relatively large and sparse LE erosion cavities, whose size and surface pattern may make difficult or impossible their experimental and numerical analysis by means of distributed roughness models alone. The key objectives are (a) validating the 3D RANS approach for predicting the aerodynamics of WT blade sections featuring relatively large and sparse erosion cavities; (b) estimating the AEP loss due to this erosion damage, typical of intermediate erosion stages; (c) assessing numerically the role of laminar-toturbulent transition and additional distributed surface roughness on eroded airfoil aerodynamics; and (d) analyzing the cavity flow and its dependence on key geometric parameters of the cavity geometry. The article is structured as follows: Section 2 reports the definition of the physical airfoil model used in the wind tunnel tests of Sareen et al. 2 and that of the digital model used herein for the 3D RANS analysis of one of those tests.
Section 3 briefly presents the CFD code, summarizes the settings used for the reported analyses, and provides in greater detail the formulation of the turbulence and transition models along with approach to incorporate the impact of distributed surface roughness in these models. Section 4 presents the adopted CFD grid topology, provides the physical domain definition, and presents the key results of the mesh sensitivity analysis.
Comparison of computed and measured force data, result sensitivity to distributed roughness and transition modeling, cavity flow patterns, and the impact of cavity flows on the overall performance of the eroded airfoil are discussed in Section 5, and the key findings of this project are summarized in Section 6. gouges are modeled as semispherical cavities with center on the airfoil surface. The edge cavity is kept sharp without chamfer. The position of the cavity centers of the digital model is also obtained by assuming a normally distributed curvilinear coordinate of the cavities in the chordwise direction and a uniform distribution in the spanwise direction. It is assumed that on the upper side the normal distribution of the curvilinear coordinate of the cavity centers achieves a value of 3 standard deviations σ at 10% c, whereas the 3σ value is achieved at 13% c on the lower side.

| ERODED AIRFOIL GEOMETRY AND ANALYZED WIND TUNNEL TEST
Since the number of cavities on both sides has been found to be very small for a value of the curvilinear coordinate s larger than 5.5% c, it has been decided to discard all cavities for s=c > s cav =c ¼ 0:0724 on both sides. The final digital model has a total of 41 pits and 27 gouges. The number of cavities on the two sides of the digital model is reported in Table 1 Figure 2A, whereas the geometry obtained by repeatedly copying and translating in the spanwise direction the baseline geometry with span equal to 10% that of the physical model is depicted in Figure 2B.
The 3D digital airfoil model and the physical domain of the CFD study are generated by using a MATLAB script, a Python script, and the  radius are provided in the MATLAB script output file. This operation results in the creation of quasi-semispherical erosion cavities in the airfoil LE surface region.
As mentioned above, the span length of the digital model used in the CFD study herein has been reduced to 10% that of the physical model used in the wind tunnel tests. To verify the validity of this modeling choice, digital models with span lengths of 5%, 10%, 15%, and 20% that of the physical model were generated, and their aerodynamic performance was analyzed by means of CFD simulations for angles of attack of À 5 , 0 , and 5 . It was found that for all three flow regimes, the lift and drag forces of the models with 10%, 15%, and 20% differed negligibly from each other, which highlighted the choice of the 10% span length to be the best trade-off of solution reliability and computational cost.

| COMPUTATIONAL AERODYNAMICS APPROACH
In this study, the computational aerodynamics of the considered eroded airfoil model is investigated solving the 3D RANS equations by means of the ANSYS FLUENT CFD code, with the discretization scheme and integration method summarized in Section 3.1. The adopted turbulence closure and transition model are discussed in Section 3.2, whereas the method used to include in the analysis the effects of distributed surface roughness is presented in Section 3.2.1. The steady pressure-based RANS equations are solved in all cases, and the turbulence closure is achieved by using the two-equation k À ω SST model, as discussed in Section 3.2. Some of the simulations below also use a two-equation model to account for the effects of laminar-toturbulent transition of the airfoil BLs. The space discretization of all conservation laws uses the second-order upwind scheme. The COUPLED solver, whereby the continuity and momentum equations are solved in a strongly coupled fashion and all other transport equations are solved in a loosely coupled fashion, is adopted for the integration of the governing equations.

| Turbulence and transition modeling
In the CFD analyses presented below, turbulence and transition modeling rely on the Langtry-Menter four-equation transitional SST model, [26][27][28] referred to as the γ À Re θ SST model hereafter. The method couples the SST turbulence model 29 with so-called local correlation-based transition modeling (LCTM). LCTM does not attempt to directly simulate the physics of laminar-to-turbulent transition, but rather relies on empirical correlations that relate local boundary layer quantities, such as momentum thickness θ and other local flow field properties such as pressure gradient and turbulence intensity, to locate the transition onset.
As explained in Menter et al., 26 one of the four equations of the γ À Re θ SST model is a transport equation for the turbulent intermittency γ.
This variable is coupled to the SST turbulence model and is used to turn on and off the production term of the turbulent kinetic energy k in the associated transport equation downstream of the transition point. A second transport equation is solved for the so-called "transition onset momentum thickness Reynolds number" f Re θt . This is done to capture the nonlocal influence of the turbulence intensity, which changes due to the decay of k in the freestream, as well as changes in the freestream velocity outside the boundary layer. This additional transport equation ties the empirical correlations used by the model to the onset criteria in the intermittency equation and allows using the model for general geometries.
The transport equations for γ and f Re θt are, respectively, where the symbols ρ, v, μ, and μ t denote, respectively, density, velocity vector, molecular, and eddy viscosity and σ f and σ θt are model constants provided in Langtry and Menter. 28 The production term P γ in Equation (1) is responsible for increasing intermittency at the transition onset location and maintaining it in turbulent regions of the BL, whereas the destruction term E γ enables modeling of BL relaminarization. The expression of both source terms is reported in Langtry and Menter. 28 The expression of the production term P θt in Equation (2) is in which t and c θt are a suitable time-scale and a model constant, respectively, and the variable Re θt is determined with an empirical correlation depending on the local turbulence intensity and Thwaites' pressure gradient coefficient λ θ , defined as with djvj=ds being the flow acceleration in the streamwise direction s. 28 The blending function F θt in Equation (3) tends to 0 outside the BL, activating the source term P θt and forcing the transported variableRe θt to approach the value Re θt obtained with the empirical correlation. In the BL, F θt tends to 1, thus switching off the production of γ and enabling diffusion ofRe θt between the BL edge and the wall boundary.
At each iteration of the solution process, the local value ofRe θt is used to determine the critical Reynolds number Re θc , the threshold above which the intermittency starts to increase in the boundary layer. This is accomplished by using another empirical correlation of the form Re θc ¼ fðRe θt Þ. 27 When the local value ofRe θt exceeds the threshold Re θc , the model triggers production of intermittency via the P γ source term in Equation (1).
The production of intermittency attempts to model the transition process by forcing the source terms of the k-equation of the SST k À ω model to remain low before transition and progressively allowing it to increase to its fully turbulent form as transition starts or progresses. The transport equations of k and the specific dissipation rate ω of this turbulence model 29 are, respectively, where σ k , σ ω , and σ ω2 are model constants and F 1 is a blending function of the k À ω and k À ϵ turbulence models. The changes brought about by LCTM in the SST transport equations are limited to (a) the production termP k and the destruction termD k in Equation (4), which are made to depend also on γ, and (b) the blending function F 1 in Equation (5), due to required additional robustness of this function when used in the framework of the γ À Re θ SST model. 26 The boundary conditions for the k and ω equations at wall and inlet boundaries are described in Menter et al. 29 At wall boundaries, a zero normal flux condition is used for both γ and f Re θt . At inlet boundaries, γ is set to 1, and f Re θt is set to the value of Re θt determined by using the empirical correlation linking this variable to turbulence intensity and using the value of this quantity at the considered boundary.
The turbulent viscosity is computed using the equation where a 1 ¼ 0: 31. 29 In all analyses of this study, however, a 1 has been set to 0.28, as this value improved significantly the agreement between computed and measured airfoil force data in all validation tests. Other studies have also reported alterations of this coefficient prompted by similar observations. 7

| Distributed surface roughness
The presence of distributed surface roughness affects the field variableRe θt , due to the increase of the momentum thickness, and also the wall boundary, due partly to larger wall viscous stress. The computed value ofRe θt is modified according tõ where f is a function of the geometric roughness height K s ; the value ofRe θ,rough is then used in place ofRe θ as input for the empirical correlations of the γ À Re θ SST model.
At rough wall boundaries, an automatic wall treatment 30,31 that blends the analytical solutions of the laminar sublayer and logarithmic layer is used to obtain a profile of the nondimensionalized velocity component parallel to the wall and the specific dissipation rate ω applicable from the outer edge of the logarithmic region to the wall. At smooth walls, when the height of the near-wall cell center is not sufficiently small to give a nondimensionalized minimum wall distance y þ of order 1, the friction velocity u τ and ω are determined from the analytical expressions of these two variables in the laminar sublayer or the logarithmic layer. Wall function estimates are not used if y þ ≈ 1. In the case of a rough wall, as surface roughness increases, the viscous sublayer is progressively destroyed. This is accounted for by artificially increasing y þ by half the dimensionless roughness height K þ s , that is, by setting This causes u τ and ω at the wall being determined from their logarithmic profiles when the roughness height is larger than the height of the nearwall cell. The law of the wall used to compute the friction velocity at rough walls differs from that of smooth wall for a constant offset ΔB depending on K þ and has form

| PHYSICAL DOMAIN, MESH TOPOLOGY, AND MESH SENSITIVITY ANALYSIS
All CFD simulations of this study are three-dimensional, and their physical domain is the region between the DU 96-W-180 airfoil surface and a closed C-shaped far field boundary at a distance of 30c from the airfoil in all directions. The spanwise extent of the physical domain is 0:186c, which corresponds to 10% of the span of the physical model of Sareen et al. 2 At both spanwise ends, the domain is closed by two plane boundaries extending from the airfoil to the far field boundary.
All grids are hybrid, consisting of prismatic, tetrahedral, pyramidal, and hexahedral elements, and are built using the ANSYS Meshing grid generation software. In the LE region, the near-wall grid consists of a thin inflation layer made up of prismatic grid elements. This layer is surmounted by an unstructured grid buffer layer made up of tetrahedra and pyramids next to its upper edge. From here, the grid has a structured pattern until the far field boundary. The grid from the remaining surface of the airfoil model to the far field boundary is fully structured. A view of the grid around the airfoil and one of the physical domains in the airfoil plane are reported in the left and right images of Figure 3A, respectively. The unstructured grid buffer in the LE region is used for both the nominal and eroded airfoil geometries, as visible in Figure 3B and 3C, respectively.
An enlarged view of the eroded airfoil LE, featuring erosion pits and gouges, is reported in Figure 3D. The surface mesh triangular elements are the bases of the prismatic grid elements used in the inflation layer in the LE region, a choice also adopted for the grid of the nominal airfoil. All erosion cavities of the damaged airfoil model have sharp circular edges, and the inflation layer covering the airfoil also covers their semispherical surface.
To assess the dependence of the CFD solutions on the grid refinement, a coarse, medium, and fine refinement have been used to determine the force coefficients of both the nominal and eroded airfoil models. In all grids, the distance of the first grid nodes off the airfoil surface from the surface itself is 4 Â 10 À6 c. For the considered case of Re  Table 3. The parameter N H is the number of grid points used to discretize the circular edge of each erosion cavity, and each entry of the column provides the minimum and maximum values of N H . The maximum values of y þ in the cavities have been found to be smaller than 1 for all three grids. It is also noted that the values of N Z used for the eroded airfoil grids, reported in Table 3, are notably larger than those used for the nominal airfoil grids, reported in inflation layer in the eroded LE region and the structured inflation layer over the remainder of the airfoil. This disparity arises from the necessity of increasing notably the surface mesh refinement in the aforementioned unstructured inflation layer and would result in the average cell size in the spanwise direction differing by about 1 order of magnitude at the interface between unstructured and structured inflation layers if the N Z values of the nominal airfoil grids were also used for the eroded airfoil grids.
As indicated above, the Reynolds number used for all simulations of this study is 1.5 M. For the mesh sensitivity analyses presented in this section, the turbulence intensity I and the turbulent/laminar viscosity ratio μ T =μ at the far field boundaries are set to 5% and 1, respectively. A velocity inlet boundary condition (BC) is applied on the C-shaped part of the far field boundary, whereas a pressure outlet BC is applied on the rear part of the far field boundary. Periodicity BCs are used on the lateral boundaries, normal to the airfoil span. The no-slip condition is enforced at the airfoil surface, and the wall is considered smooth in all mesh sensitivity analyses herein. The impact of turbulence on the mean flow is accounted for by using the k À ω SST model for the turbulence closure. 32 The effects of laminar-to-turbulent transition of the airfoil boundary layers are modeled using the γ À Re θ transition model. 26,28 All simulations reported in this article were run using 2000 iterations, using an underrelaxation factor of 0.5 for the solution of all equations. This choice resulted in the residuals of all conservation laws decreasing by 3 to 6 orders of magnitude in all reported simulations. Figure 4A reports the curves of the lift coefficient c l of the nominal airfoil computed using the coarse, medium, and fine refinement grids against the AoA α in the range À8 ≤ α ≤ 14 . It is observed that no significant differences exist between the c l curves of the medium and fine grid predictions over the considered range of AoA. Conversely, notable differences between the coarse and medium refinement grid predictions exist for α > 6 , where the coarse grid solution underpredicts c l with respect to both the medium and fine grid results. Smaller differences between the coarse and medium refinement grid solutions also exist for α < À 4 , where the magnitude of the negative lift is underpredicted by the coarse grid analysis.
The curves of c l against drag coefficient c d (drag polars) of the same airfoil geometry computed using the three grids are provided in Figure 4B. Also in this case, very good agreement of the medium and fine refinement grid predictions is observed over the considered range of AoA. The differences between the coarse and medium grid predictions are instead significant for α > 6 and α < À 4 , where the coarse grid analysis not only underestimates the magnitude of c l with respect to the medium and fine grid solutions but also overestimates c d . Figure 4C presents the curves of the ratio c l =c d of the nominal airfoil against α obtained with the coarse, medium, and fine grid analyses. Only very small differences exist between the c l =c d curves of the medium and fine grid predictions over the considered range of AoA. Conversely, significantly larger differences between the coarse and medium refinement grid predictions exist for α > 6 and, to a minor extent, α < À 4 , where the coarse grid solution underpredicts the magnitude of c l =c d with respect to both the medium and fine grid results. These differences are due to the coarse grid underestimating the magnitude of c l and overestimating c d at high and low α with respect to the other two grid analyses. The results above indicate that the medium refinement grid provides a grid-independent solution, and for this reason, all other nominal airfoil analyses presented and discussed in this study have been performed using this grid. Figure 5A reports the curves of the lift coefficient c l of the eroded airfoil computed using the coarse, medium, and fine refinement grids for À8 ≤ α ≤ 14 . Qualitative trends are similar to those observed for the nominal airfoil mesh refinement analysis. Also for the eroded airfoil geometry, differences between the coarse and medium refinement grid predictions exist for α > 6 , but these are smaller than those encountered in the nominal airfoil analysis. At the lowest AoA values, the coarse and medium grid c l values are very close to each other unlike the two corresponding grid predictions were in the case of the nominal airfoil. No significant differences exist between the c l curves of the medium and fine grid predictions over the entire range of AoA, similarly to the case of the nominal airfoil mesh refinement analysis.
The drag polars of the eroded airfoil geometry computed using the three grids are provided in Figure 5B. Very good agreement of the medium and fine refinement grid predictions is observed over the considered range of AoA. Some differences between the coarse and medium grid predictions are instead visible at the largest and, to a minor extent, smallest values of α. These differences are due to the coarse grid analysis slightly underestimating the c l magnitude and overestimating c d in these AoA ranges. Figure 5C compares the coarse, medium, and fine grid curves of c l =c d versus α for the eroded airfoil geometry. Negligible differences exist between the c l =c d curves of the medium and fine grid predictions over the considered range of AoA. Relatively small differences between the coarse and medium refinement grid predictions exist for α > 6 and, to a minor extent, α < À 4 , where the coarse grid analysis underpredicts the magnitude of c l =c d with respect to both the medium and fine grid results. These differences are due primarily to the coarse grid underestimating the magnitude of c l at high and low α with respect to the other two grid analyses. These results indicate that the medium refinement grid provides a grid-independent solution, and for this reason, all other eroded airfoil analyses presented and discussed in this study have been performed using this grid.
As discussed in Section 2, the assumptions and the analyses leading to setting the aspect ratio of the digital model of the eroded airfoil to 10% the aspect ratio of the model tested in the wind tunnel were partly motivated by limiting computational costs by avoiding the use of excessive grid cell counts. In line with that choice, it is noted that, as the medium refinement grid yielding a mesh-independent solution for the analyzed eroded airfoil model has nearly 8.4 M elements, using the aspect ratio of the tested model would require grids with 84 M or more elements, which was impractical for the study reported herein, due to the large number of parametric analyses performed over a fairly wide AoA range.

| RESULTS
In the domain is substantially larger than the wind tunnel, and the measured level of I cannot be enforced at the same distance where it was measured in the tunnel. All presented simulations refer to Re = 1.50 M, which is also the value for which the measured force data below are provided.

| Comparison with measured data, and impact of distributed roughness and transition
Measured and computed curves of c l against α and c d against α are compared in Figure 6A and 6B, respectively, whereas measured and computed drag polars and curves of c l =c d against α are compared in Figure 6C and 6D, respectively. The AoA range for which measured c l data are available is smaller than that for which measured c d data are available, 22 as observed by comparing Figure 6A and Figure  The c l curves of Figure 6A highlight that both experiments and simulations predict a negligible impact of the considered erosion pattern at low AoAs, for À3 < α < 3 . Both methods predict a significant c l drop at high AoAs, with reductions in excess of 0.2 at α ≈ 8 . As discussed below, part of the reason for this reduction is premature transition of the damaged airfoil SS BL, resulting, in turn, in thicker SS BL on the rear part of the airfoil SS and anticipated flow separation over the nominal airfoil working at the same nominal AoA. At the lowest AoAs considered, the experimental measurements do not show a significant detrimental impact of erosion on the lift force. Conversely, the simulations predict reductions of the eroded airfoil jc l j of up to about 0.1. This trend difference may be due to some discrepancies between the digital and physical erosion patterns on the PS, as the airfoil performance at low AoAs is significantly affected by the PS surface state. It is also noted that for the airfoil featuring LE delamination, which is a more severe erosion stage than that associated with the erosion cavities considered in this study, the experimental data of Sareen et al. 2 show a jc l j reduction at low AoAs. Moreover, that reduction is predicted fairly well by CFD 11 simulations of the same experiment.
The fact that the uncertainty affecting the definition of the simpler LE delamination geometry is likely to be smaller than that affecting the geometry definition of the LE erosion cavities supports the assumption that the present overprediction of the lift loss at low AoAs is due to some discrepancies between the eroded PS geometry of the physical model and that used for the CFD simulations.
The curves of Figure 6B also highlight a fair agreement of the measured and computed differences of the c d values of the eroded and nominal airfoils. One also notes that the CFD simulations reproduce well the experimentally observed trend of c d increase of the damaged airfoil as AoAs increase, although the computed drag coefficients are slightly higher than the measured ones at the highest AoAs considered. The overall agreement of measured and computed drag polars reported in Figure 6C, as well as that of the measured and computed profiles of c l =c d against α reported in Figure 6D, is fairly good. It is particularly encouraging that the experimentally measured reductions of the maximum c l =c d due to the considered erosion cavities and the AoA at which the maximum c l =c d occurs are predicted fairly well by the presented CFD analyses.
To examine in quantitative terms the measured and computed differences of force data, Figure 7A reports the c l difference between the eroded and nominal airfoils against the AoA, whereas Figure 7B reports the c d difference between the same airfoils. The curve labeled "exp." refers to the experimental data, whereas the curve labeled "dam. rough CFD" refers to the reference CFD data reported in Figure 6. These curves confirm overall fairly good agreement of measured data and reference CFD simulations. Furthermore, they both indicate a c l reduction of nearly 0.14 and a c d increase of about 0.0061 at α ¼ 6 ∘ , which is a representative operating condition for WT blade sections featuring this airfoil.
At α > 6 , however, the reference CFD simulation overpredicts the experimental level of c l reduction and c d increase.
To try and separate the impact of distributed surface roughness and the larger and sparse geometry perturbations due to the considered erosion cavities on the patterns of c l and c d variations examined above, three additional CFD simulation sets have been performed, namely, (a) one using the eroded geometry without distributed roughness at the LE but using transition modeling as the reference simulations, labeled "dam.
smooth CFD" ; (b) one using the eroded geometry without distributed roughness and without transition modeling, labeled "dam. FT CFD"; and (c) one using the nominal geometry without distributed roughness and without transition modeling, labeled "nom. FT CFD." The difference between these three force data sets and the smooth nominal airfoil CFD data set is also plotted in Figure 7. Some important observations can be made. Firstly, the drag increase and the lift reduction predicted by the transitional analysis are overall lower when no distributed roughness is included, as expected: at α ¼ 6 ∘ , a typical AoA for blade sections featuring this airfoil, the c l reduction and the c d increase predicted without distributed roughness are about 0.07 and 0.003, respectively. As shown below, the higher losses of the transitional analysis with LE roughness are due to the fact that at low AoA, roughness increases the BL thickness, decreasing lift and increasing profile drag; at higher AoA, the considered distributed roughness level also yields premature transition, indicating that at these AoA levels, the Reynolds number Re k based on the roughness height k exceeds the critical threshold Re k ð Þ cr . 18 The minimum c l and c d differences between the two transitional simulations occur at the maximum AoA of 10 considered in this analysis, indicating that at this value the flow separation occurring after stall is not significantly affected by LE roughness.
A second key note is that the drag increase and the lift reduction predicted by the fully turbulent simulations of both the nominal and eroded airfoils without roughness are significant and differ relatively little from each other over the considered AoA range. This indicates that if the distributed LE roughness is sufficiently large to trip BL transition at the LE for all AoAs (Re k < Re k ð Þ cr ), one may avoid both modeling laminar-toturbulent transition and resolving geometrically sparse and larger erosion cavities, significantly reducing computational costs of CFD simulations.
The results of Figure 7 also show that the aerodynamic losses predicted by the transitional simulation of the eroded airfoil without roughness become significant only for α ≥ 6 but remain notably smaller than those of the transitional simulations of the eroded airfoil with roughness for α > 2 and are comparable to these for lower AoAs. This suggests that when the LE of WT blades featuring the considered erosion cavity pattern has a level of distributed LE roughness that starts shifting BL transition towards the LE (Re k < Re k ð Þ cr ) only from certain AoAs up, the aerodynamic analysis needs to model BL transition and surface roughness and also resolves geometrically the erosion cavities. This is important to avoid underestimating the aerodynamic performance loss due to the considered type of LE erosion.
The relative patterns of the curves in Figure 7 for α > 6 are affected not only by the BL state but also by the development of flow separation on the airfoil SS. As the AoA increases above 6 , this flow separation grows more rapidly for the transitional analysis with LE roughness and is instead more slow for the fully turbulent cases, with the smooth LE transitional result lying in between. These observations, supported by the profiles of static pressure coefficient discussed below, also highlight the importance of modeling transition and roughness and resolving cavities when the LE distributed roughness is not sufficient to trip transition at the LE for all AoAs.  Figure 8C. Figure 8B and 8D is structured as Figure 8A and 8C, respectively, but they refer to the nominal airfoil with transitional and fully turbulent BL and the eroded airfoil with LE erosion cavities assuming smooth surface and fully turbulent BL. The transitional solution profiles of the nominal airfoil are reported in all four subplots of Figure 8 to enable a more direct cross-comparison of all five solutions. In all plots, solid and dashed lines refer to the airfoil SS and PS, respectively. At α ¼ 0 ∘ , the SS and the PS of the nominal airfoil experience transition between 50% and 60%, and between 60% and 70% chord, respectively. As visible in Figure 8C, the SS adverse pressure gradient at the position where transition occurs is relatively modest, and, therefore, detrimental impact of LE cavities and roughness do not alter the mean position where transition occurs for the two eroded airfoil set-ups ( Figure 8A). On the other hand, the erosion cavities increase the mean surface viscous stress in the front part of the airfoil, and the addition of distributed LE roughness results in the SS boundary layer becoming turbulent from onset. One also sees that a significant level of spanwise variability of the surface viscous stress exists from the LE to about 60% of the chord when the airfoil features erosion cavities. The μ cp profiles of the three transitional analyses in Figure 8C differ very little from each other, and this explains why the c l reduction at α ¼ 0 ∘ predicted by the two transitional analyses of the eroded airfoil in Figure 7A is very small. The c d increase of the same two analyses in Figure 7B is instead due to the higher surface viscous stress on the PS.
At α ¼ 0 ∘ , the surface viscous stress estimates on both the SS and the PS of the nominal and eroded airfoils predicted by the fully turbulent analyses are very close to each other and both higher than those of the transitional estimate for the nominal airfoil, as seen in Figure 8B. The closeness of the two fully turbulent c f profile predictions indicates that at this AoA, the geometric resolution of the erosion cavities does not affect significantly the aerodynamic performance. This is also supported by the fact that the spanwise variability of c f of the eroded airfoil is very small. The thicker BLs associated with fully turbulent BLs from the LE also result in slightly lower aerodynamic loading, as deduced by crosscomparing the c p profiles in Figure 8D.
The contour plots of the intermittency γ on a surface at a constant distance of 0.5 mm from the airfoil over 20% chord from the LE for α ¼ 0 ∘ are reported in the subplots of Figure  right to the airfoil PS. Nonzero values of γ indicate the onset of transition and close to 1 the BL is taken to be turbulent. The contour plots of Figure 9A and 9B highlight that at the considered low AoA, no significant turbulent activity exists in the first 20% chord of the nominal airfoil, which is consistent with the conclusions drawn inspecting the c f profiles of Figure 8A that transition occurs from 50% chord on both sides. The γ contours of the LE region with erosion cavities in Figure 9C and 9D reveal that the cavities promote turbulence, with γ reaching the upper limit of 1 behind many of the larger cavities on the airfoil PS. As seen in the corresponding c f profiles of Figure 8A, however, this enhanced turbulent activity is not sufficient to move transition to the LE on either side. The inclusion of LE surface roughness promotes turbulence even further on both sides, as seen in Figure 9E and 9F. Due to this, the PS BL of the airfoil with LE cavities and distributed roughness becomes fully turbulent from the onset, as already seen in the corresponding c f profile of Figure 8A. Figure 10, which has the same structure of Figure 8, compares the skin friction coefficient c f and the static pressure coefficient c p of the nominal and eroded DU 96-W-180 airfoils at α ¼ 6 ∘ obtained using the same physical set-ups considered in the analysis reported in Figure 8 for the weight of distributed LE roughness in triggering transition at the LE is comparable or possibly even higher than that of larger and sparser erosion cavities. The μ cp profiles of Figure 10C indicate significantly higher differences among the predicted loading than observed at α ¼ 0 ∘ . The highest and lowest loading are obtained by the nominal and eroded airfoil with added distributed roughness, respectively, and the loading of the airfoil with cavities and without distributed roughness lies in between. The lower loading of the eroded airfoil with roughness is partly due to a flow separation on the SS close to the trailing edge (TE) starting at the considered AoA. This flow separation is less pronounced for the eroded airfoil without roughness and is absent for the nominal airfoil. The aforementioned TE flow separation in the eroded airfoil case is caused by the thicker and less energetic SS BL due to premature transition, and these characteristics are exacerbated by the addition of distributed roughness.
Similarly to the α ¼ 0 ∘ case, also at α ¼ 6 ∘ the surface viscous stress estimates on both the SS and the PS of the nominal and eroded airfoils predicted by the fully turbulent analyses are very close to each other and both higher than those of the transitional estimate for the nominal airfoil, as seen in Figure 8B. The closeness of the two fully turbulent c f profiles confirms that also at this AoA, the geometric resolution of the erosion cavities does not affect significantly the aerodynamic performance. This is also supported by the fact that the spanwise variability of c f of the eroded airfoil is very small. The thicker and less energetic BLs due to LE transition also result in lower aerodynamic loading, as deduced by crosscomparing the c p profiles in Figure 10D. Cross-comparison of the μ cp profiles of Figure 10C and 10D shows that the profiles of the transitional simulation of the eroded airfoil with distributed roughness, on one hand, and the fully turbulent analyses of the nominal and eroded airfoil, on the other, are very close, confirming that, once the SS BL is turbulent from the LE, the geometric resolution of the considered cavity pattern and extent does not significantly alter the prediction of the performance reduction.
The γ contour plots of the transitional flow simulations on a surface at a constant distance of 0.5 mm from the airfoil over 20% chord from the LE for α ¼ 6 ∘ are reported in the subplots of Figure 11. Similarly to the case of Figure 9, the top, middle, and bottom rows of images refer to the nominal airfoil, the eroded airfoil with LE erosion cavities and without distributed surface roughness, and the eroded airfoil with LE erosion cavities and distributed surface roughness of 279 μ m, respectively. The images on the left refer to the airfoil SS, and those on the right to the airfoil PS. Nonzero values of γ indicate the onset of transition and close to 1 the BL is taken to be turbulent. The contour plots of Figure 11A and 11B highlight that, like at α ¼ 0 ∘ , also at α ¼ 6 ∘ no significant turbulent activity exists in the first 20% chord of the nominal airfoil, which is consistent with the corresponding c f profiles of Figure 10A, showing that transition occurs from about 40% chord on the SS and about 70% chord on the PS. The γ contours of the LE region with erosion cavities in Figure 11C show that at the considered AoA, which is typical of rated wind speed operation, the role of erosion cavities in promoting turbulent activity and transition is substantial, with γ reaching and maintaining the upper limit of 1 behind the majority of the larger cavities on the airfoil SS. As seen in the c f profiles of Figure 10A, this enhanced turbulent activity is sufficient to move the SS BL transition to the LE. The γ contours of Figure 11C also highlight a significant variation of γ in the spanwise direction, which correlates with the larger cavity positions on the surface. This variability also yields large variations of c f in the spanwise direction, which explains the relatively large standard deviation of the SS c f observed in Figure 10A. The inclusion of LE surface roughness promotes turbulent activity even further on the SS, as seen in Figure 11E. As a consequence, the SS BL of the airfoil with LE cavities and distributed roughness has a nearly constant value of γ ¼ 1 over the entire span, yielding, in turn, lower standard deviation of c f , as seen in Figure 10E.

| Estimates of annual energy production losses
To estimate the AEP loss level associated with the computed and measured aerodynamic performance loss of the DU 96-W-180 airfoil highlighted in Figure 6, an AEP analysis of the NREL 5-MW reference turbine 21 is considered. The turbine aerodynamic performance is determined using the NREL AeroDyn blade element momentum theory code, 33 assuming steady wind with no vertical shear and modeling all wind turbine components as rigid. In these analyses, the NACA 64-618 airfoil, which is used from 70% rotor radius to blade tip and has thickness of 18%, like the DU 96-W-180 airfoil, has been replaced with the latter airfoil. In the AEP loss analysis reported below, the AEP of two variants of the  Figure 12, with the profiles of normal force per unit length F n shown in Figure 12A, and those of tangential force per unit length F t shown in Figure 12B. All results refer to a close-to-rated wind speed of 11 m/s and rotor angular speed of 12.1 RPM. Also in this case, the computed force coefficients of the nominal and eroded DU 96-W-180 airfoils of Figure 6 have been used. Both plots indicate notable reductions of both force components on the outboard blade portion affected by LE erosion. These reductions, also occurring at lower wind speeds, account for the aforementioned C P reduction due to erosion. For wind speeds above the rated wind speed, the tangential and normal force components of the nominal and eroded blades also differ. However, the turbine power can be kept constant and equal to the rated value relying on the blade pitch control that curtails the power to the rated value fairly independently of LE erosion. 11 The site considered for the AEP analysis has a Rayleigh wind frequency distribution characterized by mean wind speed of 9.36 m/s. AEP estimates have been considered, all using the same wind frequency distribution and differing for the power curves computed with AeroDyn. Two power curves have been obtained using the computed estimates of the force coefficients of the nominal and eroded DU 96-W-180 airfoils reported in Figure 6, and two power curves have been obtained using the measured force coefficients of the same airfoils. The AEPs thus  Table 4. One sees that the order of magnitude of the AEP losses based on computed and measured force coefficients is comparable, as the former amounts to about 2.1% and the latter to about 2.6%.
It is noted that the DU 96-W-180 force data used for these analyses refer to Re ¼ 1:5 M, whereas the outboard blade of this turbine operates at Re ≥ 7 M. At the relatively low local AoAs of these blade sections, however, the c l and c d sensitivity to Re is modest, and the presented AEP loss levels are deemed realistic.

| Sensitivity to cavity shape and depth, and cavity flow patterns
The presented numerical analysis of the experiments reported in Sareen et al. 2 and Sapre 22 is based on the assumption of a semispherical shape with sharp edges of the erosion cavities applied to the tested airfoil. However, no specific information on this aspect has been reported. In order to investigate the uncertainty affecting the detailed geometry definition of the pits and gouges applied to the physical model of the eroded DU 96-W-180 airfoil used in the wind tunnel tests, a parametric computational analysis of the impact of different cavity shapes has been undertaken.
Five different cavity shapes have been analyzed, namely, the semispherical cavity used in all analysis of Section 5. The side view of all considered cavity types is shown in the top part of Figure 13, whereas the bottom part provides the curves of lift coefficient c l , the drag coefficient c d , and the ratio c l =c d versus the AoA α of the airfoils featuring cavities with Shapes 1-5 at the LE. Inspections of these results indicate that no lift loss is incurred for any of the considered cavity types for À2 < α < 2 and that the drag increase in the same range is very small for all cylindrical cavities and marginally higher for all semispherical cavities. As α increases, however, both the lift reduction and the drag increase of all five damaged airfoils become more severe. The smallest impairment of the aerodynamic performance is obtained with the cylindrical holes, and the lift and drag degradation with this cavity type appears to be unaffected by both the cavity depth and the chamfer. respectively. These observations, made on the basis of this particular parametric analysis, point to the fact that the impact of erosion cavities on airfoil performance depends significantly on both the overall internal geometry of the cavity and the geometry close to the airfoil surface. This is indicated by the fact that the same chamfer notably increases the aerodynamic performance impairment due to the considered semispherical cavities but does not alter the force coefficients of the airfoil featuring cylindrical cavities. Cross-comparison of the force coefficients of airfoil featuring the considered cylindrical cavities also indicates negligible dependence of the aerodynamic performance loss on the cavity depth. ing the cavity laterally and at the rear edge constitutes an obstruction for the more energetic fluid flowing above the cavity. When the main and the cavity streams meet, the main BL becomes turbulent and thickens. A similar occurrence was also noted behind the forward facing step forming in the case of LE delamination. 11 Cross-comparing the flow field patterns around the cavity edge of the five considered cavities also shows that the largest blockage resulting from the flow leaving the cavity is that of the chamfered semispherical cavity, indicating that sharper cavity edges may be less detrimental for aerodynamic performance. However, the fact that the same chamfer has negligible impact on the performance of the airfoil with cylindrical cavities indicates that the slope of the entire lateral wall of the cavity (all the way down to the cavity floor) has a significant impact on the airfoil performance. This may be due to the fact that cavities with lateral walls normal to the airfoil surface confine more effectively the recirculating cavity fluid, reducing the perturbation to the airfoil BL around the cavity. In operation, the shape of the cavity edges is likely to depend on the structural properties of the LE material at the erosion time. the section normal to the cavity axis, and the chordwise section containing the cavity axis. All results highlight the complex 3D character of the cavity flow, and the left and central subplots also point to the lack of flow symmetry discussed above. The jvj contours of the spanwise section (left subplot) also highlight the thicker boundary layer on the right-hand side of the cavity due to the blockage associated with part of the cavity fluid leaving the cavity in this region. It has been found that the cavity flow asymmetry is observed also using symmetry BCs on both lateral boundaries of the domain. In light of the symmetry of the airfoil model and the boundary conditions, the lack of flow symmetry may look surprising. However, small numerical perturbations of the numerical solution, such as those resulting from the CFD grid not being symmetric about the midspan plane, may result in a symmetric flow state not being achievable. In the context of the continuous governing conservation laws, this circumstance may correspond to the existence of a supercritical pitchfork bifurcation, 34 whereby, above a critical Reynolds number and/or AoA, the symmetric solution becomes unstable and the nonsymmetric solution with significant levels of vortical cavity flow in the plane normal to the cavity axis becomes instead stable. In the authors' view, the nonsymmetric vortical solution is more representative of real erosion cavities, featuring irregular shapes which would prevent symmetric solutions even ignoring the desymmetrizing perturbations due to rotational effects and blade twist and taper.

| CONCLUSIONS
A thorough numerical investigation into the impact of LE erosion cavities on the aerodynamic performance of wind turbine airfoils has been presented. Three-dimensional RANS CFD has been used to study the aerodynamics and the performance of the nominal DU 96-W-180 airfoil and one of its eroded variants considered in the wind tunnel experimental measurement campaign of Sareen et al. 2 The considered variant is representative of intermediate to advanced LE erosion of utility-scale wind turbines. The eroded LE simulations have both resolved the geometry of large and sparse erosion cavities and modeled lower level distributed surface roughness to assess the impact of these erosion features on laminar-to-turbulent transition and overall airfoil performance combined and in isolation. Overall good agreement between measured and computed force data of the nominal and the eroded airfoil geometries has been found. Making use of a variant of the NREL 5-MW turbine obtained by replacing the NACA 64-618 with the DU 96-W-180 airfoil, an AEP loss between 2.1% and 2.6%, based on both experiments and simulations, has been observed.
It has been found that the resolved isolated cavities can trigger LE transition of the SS BL from an AoA of about 6 , starting to decrease lift and increase drag with respect to the smooth nominal airfoil values. The addition of distributed surface roughness with roughness height smaller than the cavity characteristic size worsens the performance loss. When modeling the BLs as fully turbulent, the geometric resolution of the erosion cavities has a small impact on the aerodynamic performance at all considered AoAs, confirming indirectly that the key role of the resolved larger cavities is to trigger transition. Consequently, geometric resolution of large and sparse LE erosion cavities and consideration of laminar-toturbulent transition in the AoA range at which the outboard part of WT blades operates (2 < α < 6 ) are key to more reliable estimates of turbine performance losses. Modeling distributed surface roughness, whose roughness height is likely to be smaller than the typical dimension of sparse erosion cavities, is also important, as this feature affects the properties of BLs and, if sufficiently large, may trigger transition over the entire spanwise length affected. A simulation-based parametric study into the impact of the geometry of LE erosion cavities on the aerodynamic performance impairment, aiming to estimate the weight of the uncertainty on the erosion cavity geometry in the analysis of the wind tunnel data at hand, shows that the depth of cylindrical holes has a very small impact on the performance of the eroded airfoil. It also highlights that the aerodynamic loss depends in a strongly coupled fashion on both the overall cavity shape and the geometry of the cavity around its intersection with airfoil surface. Semispherical cavities are found to induce greater aerodynamic losses than cylindrical cavities, seemingly due to the more effective confinement of the cavity fluid associated with cylindrical cavities, which thus reduce flow perturbations of the airfoil BL. Cavity chamfers increase the performance reduction caused by semispherical cavities, whereas they have very limited impact when applied to cylindrical holes.
This parametric analysis provides a good starting base for wider and more general investigations into wind turbine blade performance losses due to blade erosion.
All CFD grids and boundary conditions used in this study have been made available at https://www.sites.google.com/site/mscampobasso/ research-data for use by other researchers.