Assessing wind turbine energy losses due to blade leading edge erosion cavities with parametric CAD and 3D CFD

Wind turbine leading edge erosion is a complex installation site-dependent process that spoils the aerodynamic performance of wind turbine rotors. This gradual damage process often starts with the formation of pits and gouges leading ultimately to skin delamination. This study demonstrates the application of open source parametric CAD functionalities for the generation of blade geometries with leading edge erosion damage consisting of pits and gouges. This capability is key to the development of high-fidelity computational aerodynamics frameworks for both advancing knowledge on eroded blade aerodynamics, and quantifying energy losses due to erosion. The considered test case is an offshore 5 MW turbine featuring leading edge pit and gouge damage in the outboard part of its blades. The power and loads of the nominal and damaged turbines are determined by means of a blade element momentum theory code using airfoil force data obtained with 3D Navier-Stokes computational fluid dynamics. An annual energy loss between about 1 and 2.5 percent of the nominal annual energy yield is encountered for the considered leading edge damages. The benefits of adaptive power control strategies for mitigating erosion-induced energy losses are also highlighted.


Introduction
Wind turbines (WTs) operate in harsh environmental conditions. Over time, rain, hail and airborne abrasive particles such as sand, erode rotor blades, particularly at the leading edge (LE). Alterations of the LE geometry due to erosion reduce significantly the blade aerodynamic performance and, thus, rotor power, with some field data and numerical studies pointing to reductions of the annual energy production (AEP) from between 2 and 3 percent, up to 10 percent [1][2][3]. LE erosion is a complex cross-disciplinary problem, involving meteorology, material science, chemistry and aerodynamics. As for aerodynamics, the performance loss mechanisms vary significantly as erosion progresses. At the early stages, when the surface geometry alterations amount to fractions of a millimeter, erosion results in increased surface roughness, which shifts the laminar-to-turbulent transition towards the LE, yielding increased drag forces [4] with respect to the design intent. As erosion progresses and the geometry alterations amount to several millimeters or more, however, the physical phenomena accounting for aerodynamic performance loss become even more complex, involving both increments of the local wall shear stress, and increased growth rate of the airfoil boundary layers due to the energy dissipated in erosion pits and gouges [5]. At advanced LE delamination stages, this faster boundary layer growth may also result in lower lift values at moderate to high angles of attack (AoA) [6]. In the presence of significant LE erosion, the amount of airfoil performance reduction due to surface roughness and larger alterations of the LE geometry may result from as yet unknown nonlinear interactions occurring on/between these two geometric scales. This LE erosion study focuses on the impact of sizeable geometry alterations of blade LE on WT power and energy yield, namely erosion-induced pits and gouges observed at most erosion stages.
The assessment of AEP losses due to lift reduction and drag increment caused by large LE geometry alterations due to erosion is receiving increasing attention, since the availability of reliable estimates of the resulting loss of revenue would enable wind farm operators to optimize maintenance and repair operations, which are particularly costly in the offshore case. Sareen et al. [2] carried out extensive experimental campaigns to assess the aerodynamic performance loss of the DU 96-W-180 airfoil subject to varying LE erosion levels, from early erosion stage pits and gouges to advanced delamination. These analyses relied on photographic records of operational multi-MW rotor blades and eroded blades under repair collected over 10+ years. Two-dimensional (2D) Navier-Stokes (NS) computational fluid dynamics (CFD) and the NREL WT hydro-aero-servo-elastic code FAST were used in [7] to assess the impact of a particular LE delamination damage on the AEP of the NREL 5 MW reference turbine. The study reported an AEP reduction of 8 percent. Wang et al. [5] developed a geometric model of pitting erosion for the 2D NS CFD analysis of the aerodynamic performance of WT airfoils featuring this damage type. Each erosion pit was modeled as a semi-circular cavity. They carried out parametric analyses of the dependence of the performance loss of the NREL S809 airfoil on the depth, surface density, surface extent and location of the erosion pits, finding critical damage levels above which airfoil performance did not undergo significant further reduction.
Recently, a novel AEP loss estimation system (ALPS) was proposed for delivering in a few seconds the AEP loss of WTs and wind farms using photographic data of operational turbine blades as input [6]. ALPS relies on aerodynamic databases of damaged airfoils generated with NS CFD, artificial neural networks exploiting the databases to rapidly determine force data of eroded blade sections not contained in the database, and blade element momentum theory (BEMT) for rotor performance analysis. LE delamination patterns of the NREL 5MW turbine blades were considered in [6]; the database had 6042 damaged airfoils, whose performance was determined for 16 AoAs in two weeks wall-clock time with 2D NS CFD using 20 16-core computer cluster nodes. For a typical offshore wind farm site in North West England and a fairly general LE delamination damage distributed non-uniformly along the blades, ALPS predicted an AEP loss of about 8 percent.
This study presents new development work of the ALPS technology, aiming at a) improving the knowledge of the energy loss mechanisms associated with pit and gouge erosion patterns, b) developing new parametric CAD approaches to damaged geometry generation, and c) further improving the ALPS portability and user-friendliness to enable assessing the performance of rotors whose blades feature complex and realistic LE erosion patterns. Key objectives are to: i) demonstrate the flexibility of a novel commercial CAD-free framework for the parametric generation of realistic LE erosion patterns including arbitrarily distributed and sized pits and gouges, and ii) analyze with three-dimensional (3D) NS CFD and BEMT the sensitivity of the performance of the NREL 5 MW turbine to the radial extent of the considered erosion pattern. The key novel contribution of this study is the direct 3D simulation of erosion cavity aerodynamics, enabled by an automated parametric CAD functionality for generating eroded blade portions, and the use of 3D NS CFD. Sections 2 and 3 discuss, respectively, the classification and CAD reconstruction of LEs featuring erosion cavities. Sections 3 and 4 present the CFD set-up and a validation study, respectively. The comparison of the power and regulation curves, and AEP estimates of the NREL 5 MW turbine with clean and realistically eroded LEs, are discussed in Section 5, followed by a study summary.

Survey and parametrization of LE erosion cavities
The erosion of WT blade LEs is a gradual process [2,8] that, once the damage exceeds the levels of high surface roughness, results in the formation of distributed pits (stage 1); as erosion progresses, pit density and size increase, leading to the formation of larger pits or gouges (stage 2); from this time, the gouges grow and/or coalesce, causing LE delamination (stage 3). Analysis of photographic images of WT blade  [2,3,8] and on commercial web sites show that this 3D damage type often has radial patterns which are roughly periodic. In this study, therefore, a parametrized model of a blade slice or strip with pitted LE is developed. The lift and drag data of the strip are computed using 3D NS CFD, and the lift and drag curves are used by a BEMT code for determining the rotor performance at the radial positions of the damaged airfoil with the pitting pattern of the strip under consideration. Each strip features five chordwise cavity rows. The pits of each row are semispherical, geometrically identical and equally spaced along the airfoil curvilinear coordinate. The parameters defining each cavity row are: 1) pit diameter d, 2) curvilinear length of the row on the suction side Lss, measured from the LE, 3) curvilinear length of the row on the pressure side Lps measured from the LE, 4) pit spacing D, namely the curvilinear distance between two adjacent pits, and 5) spanwise position zl of the pit row, taken to be the position of the pit centers.
To make a suitable choice of the pitting damage parameters for the study below, seven photographic images of real pitted LE portions were analyzed, associating to each quasi-periodic damage pattern a set of approximate values of the five geometric parameters above, normalized by d or the mean chord c. The surveyed pitting damages were categorized into three classes: small, for 0.001< d/c<0.002, medium, for 0.002<d/c<0.006, and large, for d/c>0.006. The results of the survey are reported in Tab. 1, in which s denotes the mean spanwise distance of adjacent pit rows.
The data of Tab. 1 highlight that: a) small pits present a fairly uniform distribution over the eroded LE, with isotropic spacing of about 2d to 5d in both chordwise and spanwise directions; b) medium size pits also have comparable surface density in the two directions; c) large pits have notably larger density in the spanwise than in the chordwise direction, with the former varying between 2d to 20d+. In some of the considered cases, the extent of the erosion damage appears to be comparable or slightly larger on the suction rather than the pressure side. This may appear counterintuitive, but may occur because the rain-induced erosion damage also depends on the magnitude and direction of the relative speed of the droplets when they impact on the blade, and these parameters are weather-and airfoil design-dependent. Indeed, simulations of rain-and sand-induced LE erosion also point to possible larger erosion extent on the suction side, although the erosion depth is higher on the pressure side [9,10].

Parametrized generation of blade surface with eroded leading edge
The 3D surface of each periodic strip of the eroded blade, required for the CFD grid generation, is built with a CAD functionality obtained by linking an in-house MATLAB script to the open source Computational Geometry Algorithms Library (CGAL) [11]. This functionality enables the automated parametric generation of a large number of LE geometries with erosion cavities, making use of the geometric parameters d, Lss, Lps, D and zl introduced in Section 2. The periodic strip geometry is defined in a 3D Cartesian system with origin at the LE at the spanwise location of one of the two periodic boundaries. The x-and y-axis are in the chordwise and spanwise directions, respectively.
The 3D geometry of each periodic blade strip with eroded LE requires the xy coordinates of the nominal airfoil and the values of the five damage parameters, and is accomplished following the steps of the pseudo-code below: 1. Load xy coordinates of nominal airfoil (example in Fig. 1

CFD set-up
The aerodynamic forces on the nominal and damaged strips are obtained with incompressible 3D NS ANSYS FLUENT simulations. In all analyses, the turbulence closure is accomplished with the k-ω shear stress transport (SST) model. For the simulations of nominal strips, Menter's γ-Reϑ transition model is used to account for the effects of the laminar-to-turbulent transition of the boundary layers. In all cases, a freestream turbulence intensity of 5% and a turbulent-to-laminar viscosity ratio of 1 are used. At the airfoil surface, a no-slip condition is enforced, and no wall functions are used.
Hybrid unstructured grids including hexahedra, prisms and tetrahedra are used. The C-shaped far field boundary is placed at about 30 chords from the airfoil in all analyses. A velocity inlet far field boundary condition (BC) is applied on the front, lower and upper parts of the far field boundary, whereas a pressure outlet BC is applied on the rear part. Periodicity BCs are applied on the two lateral sides of the physical domain.

Mesh refinement and validation
The analysis of the NACA 64-618 airfoil at Reynolds number (Re) of 6M is considered herein. The reference (fine) grid features NA=471 nodes on the airfoil surface at any spanwise location, and a total of NT=725,698 elements. The 3D airfoil has spanwise extent Δz=0.025c, discretized with 9 grid nodes. The minimum wall distance is 4.0x10 -6 c, which gives nondimensionalized minimum wall distance y + of about 1 for all considered AoAs. For assessing the solution sensitivity to mesh refinement, flow simulations are also performed using a coarser grid with NA=236 and NT=125,250. The minimum wall distance in the coarse grid is twice the fine grid value. Views of the LE region and the whole airfoil region of the fine grid cross section are provided, respectively, in the left and right images of Fig. 3.

Results
The developed parametric CAD functionalities for the CFD analysis of turbines whose blades have LE erosion cavities are demonstrated by determining the power, load and regulation curves of the NREL 5 MW offshore turbine with nominal and eroded blades. Three levels of LE erosion damage are assessed: a moderate one affecting the 15% outermost blade portion (damage A), a medium one affecting the 30% outermost portion (damage B), and an extended one affecting the 43% outermost portion (damage C). The blades feature the NACA 64-618 airfoil from the tip to 70% rotor radius (RT), and the DU 93-W-210 airfoil from 0.70RT to 0.57RT.

CFD analyses of nominal and eroded blade strips
The analyses of the nominal and eroded NACA 64-618 and DU 93-W-210 airfoils are performed at Re=7M and Re=6M, the same values to which the aerodynamic force data used by the NREL 5 MW turbine refer. At these Reynolds numbers, the adopted minimum wall distance of 4.0x10 -6 c of the grid nodes closest to the blade surface gives a y + value of about 1 in all analyses below.
The fine grid data of the nominal NACA 64-618 airfoil presented in Section 5 are used for the nominal turbine analysis below. For the nominal DU 93-W-210 airfoil, a similar grid with NA=471 nodes on the airfoil at all spanwise locations and a total of NT=736,795 elements is used. The 3D model of this airfoil has Δz=0.025c with 8 grid nodes covering this length. The eroded geometry of the two blade portions featuring the considered airfoils results from stacking a radially periodic pitted strip having 5 pit rows. The damage parameters of the strips used for the two airfoils are shown in Tab. 2. The grid of the eroded strip of the NACA 64-618 airfoil has Δz=0.025c and a total of 2,275,870 elements, while that of the DU 93-W-210 airfoil has Δz=0.025c and 2,425,963 elements. A cross section of the grid in and around an erosion pit of the NACA 64-618 airfoil is shown in Fig. 5.    Figure 6b compares the drag polars of the same profiles, highlighting notable cd increments of the eroded geometry. One finds that for the aforementioned AoA range, the drag coefficient increases by up to 80% (from 0.0106 to 0.0192). Figures 6c and 6d

Comparative performance analysis of turbines with nominal and eroded blade geometries
The power P, torque Q and thrust T of the nominal 5 MW WT and the WT with LE erosion geometry defined by damage C (LE erosion cavities from tip to 0.57RT) are plotted against the wind speed U in the left plot Fig. 7, where the solid line-curves refer to the nominal turbine, and the dotted line-curves refer to the damaged turbine. The rated value of P is about 5.3 MW, as the analyses refer to the aerodynamic power with no account of electrical and mechanical losses. The rotor speed Ω, the blade pitch β, and the tip-speed ratio (TSR) λ of the same two turbines are plotted against U in the right plot of Fig. 7. All curves referring to the nominal turbine are determined with the FAST code, and the power control of the nominal turbine is hardwired in the code and described in [6]. The curves of the damaged turbine are instead determined using the NREL BEMT code AERODYN and an external power control strategy implemented in MATLAB [6]. The adaptive control of the damaged turbine accounts for the alterations of the optimal TSR and power coefficient CP to mitigate power losses due to erosion, and, like the nominal controller, uses the blade pitch to limit the power at wind speed above rated. The cl and cd data for the blade sections at radii lower than 0.57RT are those provided by NREL, whereas those for the remainder of the blade are provided by the CFD analyses presented above for both nominal and damaged WTs. The key effects of blade damage C appear to be that: i) the rated power is achieved at U=12 m/s, which is higher than the rated speed of 11.4 m/s of the nominal WT and indicates significant power losses as U increases from cut-in to rated speed, and ii) the blade pitch of the damaged turbine decreases more slowly from rated to cut-out speed so as to increase AoA levels and compensate for overall lower lift and higher drag of the eroded blade sections. The λ-CP curves of the nominal WT and those featuring the erosion damage levels A, B and C are depicted in Fig. 8. As expected, the reduction of the maximum CP due to LE pitting increases with the radial extent of the erosion. However, it is noted that the CP reductions incurred moving from the nominal WT (CP=0.468) to the damage A WT (CP=0.458), and from the damage A WT to the damage B WT (CP=0.448) are about twice that incurred passing from the damage B WT to the damage C WT (CP=0.442). This is because the higher-radius blade sections contribute more to the overall rotor power, and, therefore, erosion at lower radii impacts less on CP. For the examined LE damages, fairly small TSR variations occur. However, this result is likely to be damage-dependent. Therefore, the erosioninduced variability of both CP and λ should be considered, in general, to implement loss-mitigating adaptive control strategies.  To provide a more tangible measure of the impact of the considered LE erosion damages on the AEP revenue, this metric is computed for the nominal and the three damaged turbines for a typical offshore site on the northwest coast of England, with wind frequency distribution defined by a Rayleigh distribution with mean wind speed of 9.3 m/s. The results of this analysis are reported in Tab. 3. One sees that the AEP loss of 0.92% due to damage A doubles when doubling the extent of the LE erosion cavities radially inwards (damage B), whereas a further similar increment of the damaged area in the same direction (damage C) yields a much more moderate increment of the AEP loss.

Conclusions
A new parametric CAD functionality, utilizing the open source CGAL software for generating surfaces starting from point clouds, has been developed and applied to the parametric geometry generation of wind turbine blade strips with LE affected by erosion cavities. For the first time, this functionality has enabled the assessment of the airfoil performance reduction due to erosion cavities making use of highfidelity 3D Navier-Stokes CFD. Making use of 3D Navier-Stokes simulations to obtain eroded blade force data, and BEMT codes to determine rotor performance and loads, it has been found that the AEP loss of the NREL 5 MW turbine featuring 3 levels of LE erosion cavities experiences AEP losses between 1% and 2% for the examined damage types and extents. This work brings the ALPS technology [6] for assessing wind farm revenue losses due to LE erosioninduced AEP reductions and optimizing blade inspection and repair costs closer to industrial deployment. Ongoing work includes generating force databases of 3D airfoil geometries with LE affected by both erosion cavities and delamination, and using artificial neural networks for rapidly predicting the AEP of wind farms whose turbines have blades featuring general size and pattern of surface defects.