Uncertainty Quantification of the Effects of Blade Damage on the Actual Energy Production of Modern Wind Turbines

Wind turbine blade deterioration issues have come to the attention of researchers and manufacturers due to the relevant impact they can have on the actual annual energy production (AEP). Research has shown how after prolonged exposure to hail, rain, insects or other abrasive particles, the outer surface of wind turbine blades deteriorates. This leads to increased surface roughness and material loss. The trailing edge (TE) of the blade is also often damaged during assembly and transportation according to industry veterans. This study aims at investigating the loss of AEP and efficiency of modern multi-MW wind turbines due to such issues using uncertainty quantification. Such an approach is justified by the stochastic and widely different environmental conditions in which wind turbines are installed. These cause uncertainties regarding the blade’s conditions. To this end, the test case selected for the study is the DTU 10 MW reference wind turbine (RWT), a modern reference turbine with a rated power of 10 MW. Blade damage is modelled through shape modification of the turbine’s airfoils. This is done with a purposely developed numerical tool. Lift and drag coefficients for the damaged airfoils are calculated using computational fluid dynamics. The resulting lift and drag coefficients are used in an aero-servo-elastic model of the wind turbine using NREL’s code OpenFAST. An arbitrary polynomial chaos expansion method is used to estimate the probability distributions of AEP and power output of the model when blade damage is present. Average AEP losses of around 1% are predicted mainly due to leading-edge blade damage. Results show that the proposed method is able to account for the uncertainties and to give more meaningful information with respect to the simulation of a single test case.


Introduction
Wind turbine damage has in recent years gained interest from industry and academia in an effort to keep aging wind parks around the globe productive. According to Rempel [1], in the early days of the wind energy industry there was the general misconception that once the blade is in operation, no further maintenance is required. This has changed, partly due to a considerable number of field reports that have started to surface in recent years highlighting extreme and worrying examples of early blade deterioration. For instance, Rempel states that blades as young as three years of age can show signs of wear and that blades of 87 out of 111 wind turbines in a wind farm off the shores of Denmark had to be dismantled and brought to shore after less than five years in operation due to severe leading-edge (LE) damage, as shown in Røndgaard [2]. It should be noted that blade damage is undoubtedly not the only source of uncertainty that affects the power production of a wind farm. Other common sources of uncertainty are related to environmental conditions, with uncertainties in wind speed and turbulence intensity being the main ones. As these are not the topic of the present study, which focuses specifically on the effects of blade damage, they are not included in the uncertainty quantification; however, in order to ensure that the It should be noted that blade damage is undoubtedly not the only source of uncertainty that affects the power production of a wind farm. Other common sources of uncertainty are related to environmental conditions, with uncertainties in wind speed and turbulence intensity being the main ones. As these are not the topic of the present study, which focuses specifically on the effects of blade damage, they are not included in the uncertainty quantification; however, in order to ensure that the study is up to the present simulation standards, they are accounted for using the standard procedures of the International Electrotechnical Commission (IEC).

Materials and Methods
In this section, the main details regarding the numerical modelling tools and choices that were made are given. Firstly, the numerical methods used to estimate the uncertainty associated to blade damage are introduced. The hypotheses regarding the input uncertainties are then explained. Finally, in the following subsections the numerical tools used in the required deterministic model evaluations are detailed.

Stochastic Approach
The stochastic approach exploited in the present work falls into the class of the arbitrary polynomial chaos (aPC) as implemented by Oladyshkin and Novak [13]. This technique falls into the field of the study of aleatory uncertainty, which only accounts for deviations of boundary condition and geometrical parameters. The present approach does not include the contribution of the limits of the numerical approach adopted. The deviation or the effect of such limitation have been considered as negligible. CFD has been validated and run according to best practices, including grid independence study. This approach has the advantage of providing stochastic results (or PDFs) without the need to change the algorithm of the numerical tools employed in the simulations. These kinds of approaches are generally known as non-invasive methods as reviewed by Iaccarino [14] and more detailed in Carnevale [15] and Ahfield [16]. The PDF of a specific quantity of interest is extracted by reproducing a surface response obtained by a certain number of simulations (or deterministic realization). The boundary conditions for these simulations are set to reproduce the PDF representing the aleatory parameter. The process of selecting appropriate boundary conditions is known as sampling. The sampling process is usually obtained by means of selecting the boundary condition using the Monte Carlo method filtered by the proper PDF. The approach as described implies a large number of simulations and it is not reliable for application where CFD solvers are used for each single deterministic prediction. This would require a high computational cost to complete the simulation campaign.
A strategy to overcome this limitation consists of a clever choice of the boundary conditions resulting in a limited number of simulations. The convolution of this boundary conditions is representative of a specific PDF. This approach is known in literature as the probabilistic collocation point (PCM). The PCM are obtained as quadrature points of a linear system built on the basis consisting in a set of polynomials (polynomial chaos, PC). The choice of these polynomials corresponds to make a strong assumption on how the response surface is determined. The surface response will be as the weighted functions corresponding to a specific PDF. Mathematic foundations can be found in Tatang et al. [17]. This particular approach has been successfully applied to CFD simulations in Carnevale et al. [15,18] and Salvadori et al. [19]. The particular approach proposed allows weaker hypothesis to be considered on the PDF of the aleatory parameter. The aPC only demands the existence of a finite number of moments and does not require the complete knowledge or even the existence of a probability density function. This approach has also been employed in Ahlfield et al. [16], where the stochastic behavior physical parameters are characterized by discontinuity and Gibbs phenomena. The aPC extends chaos expansion techniques by employing a global polynomial basis.
Let's consider a generic aleatory variable ξ propagating on a specific output of interest Y = f (ξ), where f is a general unknown stochastic model (or PDF); it can be expressed as a d-order expansion: According to the general theory of PCM the characteristic statistical quantities of Y(ξ) can be evaluated by the coefficient c i , and the momentum and variance are expressed as follows: The peculiarity of the aPC approach is related to the strategy adopted to determine the orthonormal basis of polynomial P (i) . These polynomials have been determined by the moment-based approach detailed in Oladyshkin et al. [13]. Once the aPC, which represents an orthonormal basis, has been identified, the collocation points are obtained by means of a quadrature procedure. Given an aleatory variable ξ ± σ associate with a PDF f (ξ), the more general expression of its quadrature is In the previous equation, the left-hand side is the stochastic representation of the aleatory variable ξ associated with the PDF f (ξ). The right-hand side is its expansion on the basis P(ξ k ), where the ω(ξ k ) is the weighting term (in this context we can consider ω(ξ k ) = 1), R M (Y) is the remainder approaching zero as d-order of the expansion increases and the collocation points ξ k are such that the formula k=0 ω(ξ k )P(ξ k ) = 0 is satisfied for the moment µ(ξ) and the µ(ξ) ± σ.

Probability Density Functions
The random input variables are introduced in the model using PDFs. Although this is not specifically required by the adopted aPC method, which is on the other hand able to operate on any kind of available data, in the present study PDFs were assumed based on an expert's opinion due to the lack of publicly available information regarding the studied parameters. In fact, the PDFs are based on the assumptions of Bortolotti et al. [20], who also attempt to deal with input uncertainties in aero-servo-elastic wind turbine models. Two beta functions are used for both LE Erosion Factor ε and TE Damage Factor τ. They are appropriately scaled to match the support these variables are defined upon. The values of the PDFs are reassumed in Table 1. The adopted PDFs are assumed as representative of cases where medium-low blade damage is present or of sites with challenging environmental conditions where regular maintenance is performed.

Blade-Damage Model
The first stage of the modelling process consists of modelling the blade damage itself. Blade damage is modelled through shape-modification of selected airfoils along the wind turbine's blades. Trailing-edge damage is reproduced by a simple truncation of the airfoil's trailing-edge. The amount of TE truncation with respect to the airfoil's cord is expressed as the above-introduced TE Damage Factor τ. Leading-edge damage is instead modelled through a more complex shape modification. This, which is caused by leading-edge delamination, is based on two main parameters, the maximum erosion depth θ and the chord-wise coverage of the damaged area ε. Both the influence of ε and of θ are studied in this research. However, in order to quickly estimate the global influence of leading-edge damage with respect to trailing-edge damage and to keep the analysis synthetic with only two random input variables, θ and ε were related through an empiric correlation. This assumption is supported by existing studies, where these two variables seem to be related to each other. However, Gaudern et al. [21] and Sareen et al. [5] found two very different ε-θ curves, as shown in Figure 2. As both curves are found by field examination of the blades and given that there is no clear way of assessing which of the two curves is more accurate for this study, a mean curve is proposed here (also shown in Figure 2). We can now refer only to ε as the aforementioned LE Erosion Factor, as this value now also uniquely determines θ.
damage is modelled through shape-modification of selected airfoils along the wind turbine's blades. Trailing-edge damage is reproduced by a simple truncation of the airfoil's trailing-edge. The amount of TE truncation with respect to the airfoil's cord is expressed as the above-introduced TE Damage Factor τ. Leading-edge damage is instead modelled through a more complex shape modification. This, which is caused by leading-edge delamination, is based on two main parameters, the maximum erosion depth θ and the chord-wise coverage of the damaged area ε. Both the influence of ε and of θ are studied in this research. However, in order to quickly estimate the global influence of leadingedge damage with respect to trailing-edge damage and to keep the analysis synthetic with only two random input variables, θ and ε were related through an empiric correlation. This assumption is supported by existing studies, where these two variables seem to be related to each other. However, Gaudern et al. [21] and Sareen et al. [5] found two very different ε-θ curves, as shown in Figure 2. As both curves are found by field examination of the blades and given that there is no clear way of assessing which of the two curves is more accurate for this study, a mean curve is proposed here (also shown in Figure 2). We can now refer only to ε as the aforementioned LE Erosion Factor, as this value now also uniquely determines θ. Once a value of ε is selected, the damaged airfoil shape is generated with a purposely developed Matlab ® tool. The LE of the airfoil is moved inward by a maximum depth of θ. Similarly to what was done by Schramm et al. [22], the leading-edge was flattened. The height of the flattened area is imposed to be ℎ = 2. Damage extends up to ε on the suction side of the airfoil and up to 1.3ε on the pressure side, as done in [5]. This is also motivated by the fact that wind turbine airfoils are designed to operate with a positive angle of attack (AoA), and therefore, the pressure side of the airfoil is more exposed to wear. The depth of delamination at the end of the damaged area is equal to = /3.
The TE and LE damage models are shown in Figure 3. The models are also described in further detail in [10]. The present model is a simplified version of the real LE damage pattern adequate for a parametric study like the present one, which cannot therefore reproduce all the features of a real, three-dimensional damaged blade. The model, however, is in line with the proposals of other authors [22,23] and also qualitatively reproduces the damaged shapes obtained from computational models [8,24], as seen in experiments [5]. Once a value of ε is selected, the damaged airfoil shape is generated with a purposely developed Matlab ® tool. The LE of the airfoil is moved inward by a maximum depth of θ. Similarly to what was done by Schramm et al. [22], the leading-edge was flattened. The height of the flattened area is imposed to be h = 2θ. Damage extends up to ε on the suction side of the airfoil and up to 1.3ε on the pressure side, as done in [5]. This is also motivated by the fact that wind turbine airfoils are designed to operate with a positive angle of attack (AoA), and therefore, the pressure side of the airfoil is more exposed to wear. The depth of delamination at the end of the damaged area is equal to D end = θ/3. The TE and LE damage models are shown in Figure 3. The models are also described in further detail in [10]. The present model is a simplified version of the real LE damage pattern adequate for a parametric study like the present one, which cannot therefore reproduce all the features of a real, three-dimensional damaged blade. The model, however, is in line with the proposals of other authors [22,23] and also qualitatively reproduces the damaged shapes obtained from computational models [8,24], as seen in experiments [5].  The lift and drag coefficients of the airfoils are calculated using CFD. The numerical set-up was 196 used by the authors and has been presented in detail in [10]; however, the main parameters will be 197

CFD Setup
The lift and drag coefficients of the airfoils are calculated using CFD. The numerical set-up was used by the authors and has been presented in detail in [10]; however, the main parameters will be reassumed herein. The ANSYS ® FLUENT ® (Version 18.2) solver is used to calculate the 2D polars. A Reynolds-averaged Navier-Stokes (RANS) approach is used. The Navier-Stokes equations are solved in a coupled manner with second order upwind spatial discretization. Turbulence closure is achieved with the k-ω Shear Stress Transport (SST) model. A bullet-shaped computational domain is used, as with this shape open-field conditions can be modelled with only one inlet and one outlet boundary condition. In order to ensure that the boundary conditions do not influence the results, the computational domain is 74 chord lengths long and 40 chord lengths wide, as shown in Figure 4a.

CFD setup 195
The lift and drag coefficients of the airfoils are calculated using CFD. The numerical set-up was 196 used by the authors and has been presented in detail in [10]; however, the main parameters will be 197 reassumed herein. The ANSYS ® FLUENT ® (Version 18.2) solver is used to calculate the 2D polars. A 198 Reynolds-averaged Navier-Stokes (RANS) approach is used. The Navier-Stokes equations are 199 solved in a coupled manner with second order upwind spatial discretization. Turbulence closure is

206
An unstructured triangular mesh is used. The airfoil's boundary layer is modelled with a 207 quadrilateral inflation layer from the blade surface. A total amount of 46 prismatic layers are used. 208 To ensure grid independence, three meshes were tested with varying number of elements. A coarse mesh 209 with 1.3 × 10 5 elements and 500 elements along the airfoil surface, a medium mesh with 2.8 × 10 5 elements 210 and 650 elements along the airfoil's surface, and a fine mesh with 3.6 × 10 5 elements and 750 elements 211 along the airfoil's surface. The lift (Cl) and drag (Cd) coefficients are calculated with CFD between 20° 212 and 30° of AoA; values for AoA higher and lower than this are extrapolated using Viterna's method 213 [25]. A total roughness height of 0.4 mm is imposed on the airfoil's nose trough an equivalent sand-214 grain roughness height, estimated through the simple correlations provided in [26]. This roughness 215 height is selected based on the observations of several authors [5, 6,21] and models medium to 216 advanced pitting and gauging of the LE. As the focus of the LE damage model is on advanced stages 217 of damage, a constant value of roughness was considered suitable across all the LE-damaged cases. 218 An unstructured triangular mesh is used. The airfoil's boundary layer is modelled with a quadrilateral inflation layer from the blade surface. A total amount of 46 prismatic layers are used. To ensure grid independence, three meshes were tested with varying number of elements. A coarse mesh with 1.3 × 10 5 elements and 500 elements along the airfoil surface, a medium mesh with 2.8 × 10 5 elements and 650 elements along the airfoil's surface, and a fine mesh with 3.6 × 10 5 elements and 750 elements along the airfoil's surface. The lift (C l ) and drag (C d ) coefficients are calculated with CFD between 20 • and 30 • of AoA; values for AoA higher and lower than this are extrapolated using Viterna's method [25]. A total roughness height of 0.4 mm is imposed on the airfoil's nose trough an equivalent sand-grain roughness height, estimated through the simple correlations provided in [26]. This roughness height is selected based on the observations of several authors [5, 6,21] and models medium to advanced pitting and gauging of the LE. As the focus of the LE damage model is on advanced stages of damage, a constant value of roughness was considered suitable across all the LE-damaged cases.
It is important to point out that CFD is by its nature deterministic, i.e., the same simulation is expected to give the same results if the same settings are used. In this sense, it does not add any source of uncertainty in the analysis. On the other hand, it is true that using different numerical settings to solve the same test case could lead to different results. On this basis, it is very important that the CFD approach is robust and validated with experiments whenever possible. In the present study, in particular, the numerical set-up was validated with respect to available experimental data from [5]. The clean and eroded data is obtained from the DU96W-180 airfoil that was tested in clean and damaged configurations ("stage 5" erosion in [5]) for a Reynolds number of 1.5 × 10 6 . Figure 4b demonstrates good agreement between the experimental values and CFD predictions, with limited differences that can be attributed to the unspecified wind tunnel turbulence level and to the surface finish of the reference model.

Aeroelastic Setup
The lift and drag coefficients of the damaged airfoils are used in an aero-servo-elastic model of the DTU 10 MW RWT [11]. The model is developed within NREL's open-source simulation code OpenFAST (NREL, CO, USA) [12]. The aerodynamic module, AeroDyn is based on blade element momentum (BEM) theory. As in all BEM codes, the wake is modelled with a series of concentric annuli, upon which a momentum balance is imposed. The blades are modelled trough lift and drag coefficients. Corrections for high induction (Glauert correction), blade tip and root losses, tower shadow, skewed flow and dynamic stall are included [27]. The coefficients of the dynamic stall model are tuned based on the lift, drag and moment coefficients of the damaged airfoils. Blade damage was considered from 70% of the rotor span outwards. The reasoning behind this choice has to do with the fact that the LE damage phenomena considered are mainly related to erosion, which is most influent where the local blade inflow velocities are highest. Other authors also applied damage from 70% of the rotor span outwards [6]. The lift and drag coefficients of the damaged airfoils are applied uniformly to the entire damaged area.
Fully flexible blades and tower are modelled with the structural dynamics module ElastoDyn. The modal formulation allows for a fairly accurate computation of the structural dynamics with very low computational cost. The Delft Research Controller (DRC) [28] is used in this study. This open-source baseline controller is able to regulate torque and pitch. Constant-torque operation is selected above rated wind speed. The control parameters are tuned based on the report of [29].

General DLC Setup
The DTU 10 MW RWT is a state-of-the-art reference rotor, developed in recent years as a benchmark for researchers and industry in the field of wind energy. It features a 178-meter diameter rotor with aerodynamic features like gurney flaps that help this conceptual turbine reach a rated power of 10 MW at a wind speed of 11.4 m/s. The tower height is 119 m and the nominal revolution speed is 9.6 rpm, which equates to a tip speed just shy of 90 m/s. The complete definition of the turbine and all of its parameters can be found in Bak et al. [11]. To estimate the AEP of the turbine, a power-production design load case (DLC) is simulated. This is done through sixty-six 10-minute simulations with wind speeds at a hub height between 4 and 24 m/s. Six turbulent seeds per wind speed are simulated, in compliance with the minimum requirements of the IEC 61400-1 [30]. The wind fields also feature wind shear and misaligned flow with respect to the rotor plane. By simulating several cases, uncertainties regarding atmospheric conditions are dealt with, and their influence is accounted for in this study.
It is important to note that turbulence affects power production and other key turbine figures in a complicated manner, as this depends both on the interaction between the controller and the incoming wind speed and on the complex blade boundary-layer phenomena amongst other things. The interaction between large turbines and the turbulent atmospheric boundary layer is out of the interest of the present study and has been evaluated in detail by Churchfield et al. and Nandi et al. [31,32]. Moreover, as other authors have pointed out when studying a similar multi-MW wind turbine in an aero-servo-elastic modelling framework [20], six turbulent realizations are enough to guarantee good convergence on the AEP statistics.
AEP is calculated using a Rayleigh wind-speed probability density function with a mean of 10 m/s as specified by IEC class IA, which is the design class of the DTU 10MW. The AEP obtained using a Rayleigh distribution with a mean wind speed of 8.5 m/s (corresponding to IEC class IIA) will also be briefly analyzed as this could be more representative of the impact of blade damage on sites with lower mean wind speeds.

Results
The aPC resulting collocation points are qualitatively shown in Figure 5 and detailed in Table 2. For each point, the corresponding damaged airfoil geometry is generated and CFD calculations were Energies 2020, 13, 3785 8 of 17 performed as described in Section 2.3. With the resulting airfoil data, aero-servo-elastic BEM simulations were performed as described in Section 2.4. be briefly analyzed as this could be more representative of the impact of blade damage on sites with 273 lower mean wind speeds. 274

Results 275
The aPC resulting collocation points are qualitatively shown in Figure 5 and detailed in Table 2.  276 For each point, the corresponding damaged airfoil geometry is generated and CFD calculations were 277 performed as described in Section 2.3. With the resulting airfoil data, aero-servo-elastic BEM 278 simulations were performed as described in Section 2.4. 279 280 Figure 5. Arbitrary polynomial chaos (aPC) resulting collocation points' plot in ε-τ space.

Aerodynamic performance 283
In this section the aerodynamic performance under uncertainties is discussed. In Figure 6 the 284 mean variation in power coefficient (CP) with respect to the clean reference turbine is shown. The 285 standard deviation and associated probability contours are also shown. The CP mean value is lower 286 than the nominal one for all the wind speed bins except for the 4 m/s one. In this wind speed bin, the 287 average gain is about 1%. The reasons that cause such gains are related mainly to the TE damage; 288 however, this gain in performance, while conceptually interesting, is weakened by two factors. First, 289 at 4 m/s the power is about 60 times lower than the nominal one, and thus the effect on the AEP will 290 be minimal. This can be seen clearly in Figure 7. Secondly, there is a high dispersion in the CP values 291 and therefore the expected value is hard to predict. The high dispersion is due to the extremely 292 different response from the damaged airfoils. Both gain and power losses at this wind speed occur. 293 The time averaged AoA from 30% of the blade span to tip goes from 0° to 5°. This allows some of the 294 damaged airfoils to operate with favourable lift and drag forces with respect to others. More details 295 about this behaviour are given below. 296 The highest value for the mean decrease in CP is of −2.6% at 10 m/s. At this wind speed the 297 reduction in CP can exceed −12%. Moreover, from 8 m/s to 12 m/s, mostly only power losses occur. In 298 this wind speed range, a significant part of the total turbine's energy is produced; therefore, power 299

Aerodynamic Performance
In this section the aerodynamic performance under uncertainties is discussed. In Figure 6 the mean variation in power coefficient (C P ) with respect to the clean reference turbine is shown. The standard deviation and associated probability contours are also shown. The C P mean value is lower than the nominal one for all the wind speed bins except for the 4 m/s one. In this wind speed bin, the average gain is about 1%. The reasons that cause such gains are related mainly to the TE damage; however, this gain in performance, while conceptually interesting, is weakened by two factors. First, at 4 m/s the power is about 60 times lower than the nominal one, and thus the effect on the AEP will be minimal. This can be seen clearly in Figure 7. Secondly, there is a high dispersion in the C P values and therefore the expected value is hard to predict. The high dispersion is due to the extremely different response from the damaged airfoils. Both gain and power losses at this wind speed occur. The time averaged AoA from 30% of the blade span to tip goes from 0 • to 5 • . This allows some of the damaged airfoils to operate with favourable lift and drag forces with respect to others. More details about this behaviour are given below.
The highest value for the mean decrease in C P is of −2.6% at 10 m/s. At this wind speed the reduction in C P can exceed −12%. Moreover, from 8 m/s to 12 m/s, mostly only power losses occur. In this wind speed range, a significant part of the total turbine's energy is produced; therefore, power reductions in this region will eventually lead to a significant reduction in AEP. Finally, for wind speeds higher than 14 m/s, shown in the grey-shadowed region in Figure 6, the damage effects are no longer visible, as from this wind speed onwards a lower pitch-to-feather regulation is able to compensate for the aerodynamic losses.
The power output per wind speed bin is shown in Figure 7. Upon examination of this figure, it is apparent that the blade damage has a greater impact on power output between 8 m/s and 12 m/s, confirming what was seen in the relative trends of Figure 6. At 4 m/s, however, as previously pointed out, the mean power output is only 174 kW, higher than the 172 kW of the nominal case. Due to the little power produced, this difference as well as the high standard deviation of ±7 kW (±4%) are not visible in the plot, further highlighting how such variation has little impact on the overall performance. In order to better understand the global results, each wind speed bin can be examined more in detail. Energies 2020, 13, x FOR PEER REVIEW 9 of 18 reductions in this region will eventually lead to a significant reduction in AEP. Finally, for wind 300 speeds higher than 14 m/s, shown in the grey-shadowed region in Figure 6, the damage effects are 301 no longer visible, as from this wind speed onwards a lower pitch-to-feather regulation is able to 302 compensate for the aerodynamic losses. 303 304 Figure 6. Variation in power coefficient, mean value (μ), standard deviation (σ) and probability.

305
The power output per wind speed bin is shown in Figure 7.

316
The response surfaces reporting the differences in CP for the wind speed bins that show the most 317 relevant differences are shown in Figure 8. For the wind speed bin of 4 m/s the response surface 318 slightly overestimates the CP of the nominal geometry. Such behavior is shown in Figure 8a around 319 the ε = 0, τ = 0 point. On the other hand, the response surface prediction gives good results at 8 m/s 320 and 10 m/s where the CP variation predicted for the nominal geometry is zero as expected. 321

305
The power output per wind speed bin is shown in Figure 7. Upon examination of this figure, it 306 is apparent that the blade damage has a greater impact on power output between 8 m/s and 12 m/s, 307 confirming what was seen in the relative trends of Figure 6. At 4 m/s, however, as previously pointed 308 out, the mean power output is only 174 kW, higher than the 172 kW of the nominal case. Due to the 309 little power produced, this difference as well as the high standard deviation of ±7 kW (±4%) are not 310 visible in the plot, further highlighting how such variation has little impact on the overall 311 performance. In order to better understand the global results, each wind speed bin can be examined 312 more in detail. 313

316
The response surfaces reporting the differences in CP for the wind speed bins that show the most 317 relevant differences are shown in Figure 8. For the wind speed bin of 4 m/s the response surface 318 slightly overestimates the CP of the nominal geometry. Such behavior is shown in Figure 8a around 319 the ε = 0, τ = 0 point. On the other hand, the response surface prediction gives good results at 8 m/s 320 and 10 m/s where the CP variation predicted for the nominal geometry is zero as expected. 321 The response surfaces reporting the differences in C P for the wind speed bins that show the most relevant differences are shown in Figure 8. For the wind speed bin of 4 m/s the response surface slightly overestimates the C P of the nominal geometry. Such behavior is shown in Figure 8a around the ε = 0, τ = 0 point. On the other hand, the response surface prediction gives good results at 8 m/s and 10 m/s where the C P variation predicted for the nominal geometry is zero as expected.

323
In the 4 m/s wind speed bin, an increase in CP for several combinations of ε and τ can be noted. 324 To explain this unexpected trend, one can consider the collocation point pairs γ 7 & γ 2 (same ε and 325 the highest and the lowest τ, respectively) and γ10 & γ3 (same τ and the highest and the lowest ε, 326 respectively). Therefore, looking at the pair γ2 & γ7 the influence of τ is highlighted, while looking 327 at the pair γ10 & γ3 the influence of ε is highlighted. Point γ7 shows the highest increase in CP (about 328 10%), while γ2 shows a mild decrease in CP, about −1.5%; thus, as shown in Figure 8, power increases 329 as tau increases. The other γ-pair shows the opposite behavior, for γ10, the power coefficient 330 decreases by 12%, while γ3 shows an increase in the power coefficient of about 3%, and thus, power 331 In the 4 m/s wind speed bin, an increase in C P for several combinations of ε and τ can be noted. To explain this unexpected trend, one can consider the collocation point pairs γ 7 & γ 2 (same ε and the highest and the lowest τ, respectively) and γ10 & γ3 (same τ and the highest and the lowest ε, Energies 2020, 13, 3785 10 of 17 respectively). Therefore, looking at the pair γ2 & γ7 the influence of τ is highlighted, while looking at the pair γ10 & γ3 the influence of ε is highlighted. Point γ7 shows the highest increase in C P (about 10%), while γ2 shows a mild decrease in C P , about −1.5%; thus, as shown in Figure 8, power increases as tau increases. The other γ-pair shows the opposite behavior, for γ10, the power coefficient decreases by 12%, while γ3 shows an increase in the power coefficient of about 3%, and thus, power decreases as tau decreases. To better understand the trends, the lift and drag coefficients for the FFAW3-241 airfoil (i.e., the airfoil present in the damaged part of the blade) for the four damage levels are shown in Figure 9 with respect to the reference configuration. In general, lift decreases and drag increases for all of the damaged configurations as expected. Focusing on the mean AoA recorded for the various damaged configurations at 4 m/s in Figure 9a, it is clear how the mean AoA increases for all of the damaged cases. This is due to the lower lift of the damaged cases. A new operational equilibrium point in the BEM code is then reached, with a lower induction and thus a higher AoA.

323
In the 4 m/s wind speed bin, an increase in CP for several combinations of ε and τ can be noted. 324 To explain this unexpected trend, one can consider the collocation point pairs γ 7 & γ 2 (same ε and 325 the highest and the lowest τ, respectively) and γ10 & γ3 (same τ and the highest and the lowest ε, 326 respectively). Therefore, looking at the pair γ2 & γ7 the influence of τ is highlighted, while looking 327 at the pair γ10 & γ3 the influence of ε is highlighted. Point γ7 shows the highest increase in CP (about 328 10%), while γ2 shows a mild decrease in CP, about −1.5%; thus, as shown in Figure 8, power increases 329 as tau increases. The other γ-pair shows the opposite behavior, for γ10, the power coefficient 330 decreases by 12%, while γ3 shows an increase in the power coefficient of about 3%, and thus, power 331 decreases as tau decreases. To better understand the trends, the lift and drag coefficients for the 332 FFAW3-241 airfoil (i.e., the airfoil present in the damaged part of the blade) for the four damage levels 333 are shown in Figure 9 with

340
As a consequence of the increased AoA, lift and drag forces slightly increase and, more 341 importantly, are more tangentially and axially oriented. The new force composition generates more 342 torque and more power for some of the combinations of ε and τ. As shown in Figure 10, the same 343 phenomena are occurring for all the damaged configurations: a change in the lift and drag coefficients 344 As a consequence of the increased AoA, lift and drag forces slightly increase and, more importantly, are more tangentially and axially oriented. The new force composition generates more torque and more power for some of the combinations of ε and τ. As shown in Figure 10, the same phenomena are occurring for all the damaged configurations: a change in the lift and drag coefficients leads to a different BEM equilibrium point with different induction and AoA along the entire area of the blade affected by damage. However, increasing ε also significantly increases drag, leading to lower performance and offsetting the benefit of a higher AoA, despite the change of orientation of the forces. For instance, in γ10, the highest increases in drag are observed, exceeding 30% at an AoA of around 2 • .
In Figure 11 the average AoA for the nominal and four damage combinations for all the wind speed bins is shown. For all the damaged combinations, the highest average AoAs are predicted in the 8 m/s and 10 m/s wind speed bins. At 8 m/s mean wind speed the average AoA for the nominal case at 78% blade span is about 6.9 • , while the damaged cases work at an even higher AoA due to decreased induction, as previously discussed. In these wind speed bins, there is no power increase in any combination of ε and τ. From the analysis of Figure 9, the higher the AoA, the wider the difference is in lift and drag coefficients. This ultimately leads to the power losses observed in Figures 6-8, with peaks that exceed −10% at 8 m/s and −12% at 10 m/s, respectively. It is also interesting to note that ε is the main cause of performance decrease and has a more pronounced effect than τ. This is due to the fact that LE damage causes a reduction in the stall AoA of the airfoil, which strongly influences high-AoA operation and a more pronounced increase in drag than TE damage.
The probability distributions found from the evaluation of the computed response surfaces at 4 m/s, 8 m/s and 10 m/s are shown in Figure 12. At 4 m/s, the variation in C P is most affected by uncertainties.
The peak is located at 1% of variation in C P , but the resulting distribution is fat-tailed. Indeed, the standard deviation is ±4.1% and the probability to lose or gain C P are about 40% and 60%, respectively. At 8 m/s and 10 m/s, the distributions are strongly asymmetric and have lower standard deviations with respect to the 4m/s case and are equal to ±1.7% and ±2% at 8 m/s and 10 m/s, respectively. In both cases, the probability for a C P gain is zero and losses always occur. They both have a marked left tail, but a higher dispersion at 10 m/s is found. The probability peak is clearly located on the right of the mean value at −1.5% and −1.7% for 8 m/s and 10 m/s, respectively. case at 78% blade span is about 6.9°, while the damaged cases work at an even higher AoA due to 355 decreased induction, as previously discussed. In these wind speed bins, there is no power increase in 356 any combination of ε and τ. From the analysis of Figure 9, the higher the AoA, the wider the difference 357 is in lift and drag coefficients. This ultimately leads to the power losses observed in Figure 6-8, with 358 peaks that exceed −10% at 8 m/s and −12% at 10 m/s, respectively. It is also interesting to note that  359 is the main cause of performance decrease and has a more pronounced effect than . This is due to 360 the fact that LE damage causes a reduction in the stall AoA of the airfoil, which strongly influences 361 high-AoA operation and a more pronounced increase in drag than TE damage. 362 case at 78% blade span is about 6.9°, while the damaged cases work at an even higher AoA due to 355 decreased induction, as previously discussed. In these wind speed bins, there is no power increase in 356 any combination of ε and τ. From the analysis of Figure 9, the higher the AoA, the wider the difference 357 is in lift and drag coefficients. This ultimately leads to the power losses observed in Figure 6-8, with 358 peaks that exceed −10% at 8 m/s and −12% at 10 m/s, respectively. It is also interesting to note that  359 is the main cause of performance decrease and has a more pronounced effect than . This is due to 360 the fact that LE damage causes a reduction in the stall AoA of the airfoil, which strongly influences 361 high-AoA operation and a more pronounced increase in drag than TE damage. 362  Figure 12. At 4 m/s, the variation in CP is most affected by uncertainties. 367 The peak is located at 1% of variation in CP, but the resulting distribution is fat-tailed. Indeed, the 368 standard deviation is ±4.1% and the probability to lose or gain CP are about 40% and 60%, 369 respectively. At 8 m/s and 10 m/s, the distributions are strongly asymmetric and have lower standard 370 deviations with respect to the 4m/s case and are equal to ±1.7% and ±2% at 8 m/s and 10 m/s, 371 respectively. In both cases, the probability for a CP gain is zero and losses always occur. They both 372 have a marked left tail, but a higher dispersion at 10 m/s is found. The probability peak is clearly 373 located on the right of the mean value at −1.5% and −1.7% for 8 m/s and 10 m/s, respectively. 374

Annual Energy Production (AEP)
The uncertainties in AEP estimation are discussed in this section. The AEP was calculated according to IEC 61400-1 standard turbine classes. A Weibull wind speed distribution with shape factor of 2 and average values of 8.5 m/s and 10 m/s were used to model sites of IEC wind class II and IA. In particular, class IA is chosen as this is the design class of the DTU 10 MW RWT and class IIA is chosen as representative of medium wind speed sites, where such a turbine might also be installed. The availability factor was assumed to be 1. This assumption is justified by the fact that relative variations are mainly analyzed in the present study. The variation in AEP for the two wind distributions is shown in Figure 13. Both response surfaces well predict the trends around ε = 0, τ = 0, showing no variation in AEP in that point. The LE erosion, ε, is the main driver for AEP reduction, as decreases are mostly along the ε axis. The trailing edge damage, τ, has a minor influence in AEP, as clearly visible in Figure 13. Moreover, the trailing edge damage contribution seems to be dependent on the erosion level. For instance, if one considers the six combinations of ε and τ (where ε = 0, 4, 8 and τ = 0, 3) for wind class IIA shown in Table 3, the point ε = 4, τ = 0 gives a variation in AEP of −1.87%, while the point ε = 4, τ = 3 gives a variation of −2.24%. Therefore, for ε = 4, the trailing edge damage increases losses by 0.37%. By performing the same consideration for ε = 8, trailing edge damage increases losses by 0.82%. This means that the TE contribution to losses increases as ε increases. Table 3. Annual Energy Production (AEP) reduction for some of the computed ε-τ combinations. The highest variation in AEP predicted by the response surface is −10.35% at ε = 8, τ = 3 for class IIA, as seen in Figure 13a. The highest simulated AEP reduction is −6.21% for γ10 (class IIA, Figure 13a). For wind class IA, the highest variations in AEP are lower than the ones predicted for class IIA and amount to −8.56% in ε = 8, τ = 3 and −5.02 in γ10, as shown in Figure 13b. Such differences are due to the Weibull wind speed PDFs. The probability for the machine to work in the bins range from 8 to 12 m/s, where the highest losses in power occur, are 36% and 31% for IIA and IA, respectively. This difference is the main cause of different variations in AEP for the two classes.   The highest variation in AEP predicted by the response surface is −10.35% at ε = 8, τ = 3 for class 396 IIA, as seen in Figure 13a. The highest simulated AEP reduction is −6.21% for γ10 (class IIA, Figure  397 13a). For wind class IA, the highest variations in AEP are lower than the ones predicted for class IIA 398 and amount to −8.56% in ε = 8, τ = 3 and −5.02 in γ10, as shown in Figure 13b. Such differences are 399 due to the Weibull wind speed PDFs. The probability for the machine to work in the bins range from 400 8 to 12 m/s, where the highest losses in power occur, are 36% and 31% for IIA and IA, respectively. 401 This difference is the main cause of different variations in AEP for the two classes.

403
Finally, we consider the probability distributions AEP variation shown in Figure 14. As is also 404 the case for the previously shown distributions, the PDFs are obtained by sampling the response 405 surfaces 250,000 times. The mean and standard deviations of the PDFs are −1.21% and ±1.04% for 406 class IIA and −0.98% and ±0.84% for class IA. Such mean reductions are indeed significant on a multi-407  Finally, we consider the probability distributions AEP variation shown in Figure 14. As is also the case for the previously shown distributions, the PDFs are obtained by sampling the response surfaces 250,000 times. The mean and standard deviations of the PDFs are −1.21% and ±1.04% for class IIA and −0.98% and ±0.84% for class IA. Such mean reductions are indeed significant on a multi-MW scale turbine and are in line with the finding of Eisenberg et al. [33] but seem to be lower than the values indicated by most of the present research [5, 6,8]. Both distributions show a clear peak, with the mode of the PDFs below 1% AEP loss in both cases. For both the IEC 61400-1 IA and IIA scenarios, the left tails of the distributions are long, reaching values of 6-8% AEP reductions, coherently with the response surface shown in Figure 13. The probability associated to values of AEP reduction in the order of 3-8%, which most authors indicate, is almost insignificant in the present test case. It is important to stress that these results depend on the assumed PDFs, which are, as discussed, based on published literature and appear reasonable based on the authors' experience. Moreover, as Fiore and Selig have suggested [34], larger turbines seem to be impacted less by LE damage phenomena such as erosion. However, results suggest that the commonly forecasted reductions might be based on heavy-damage scenarios, which, whilst not unrealistic, have low probability of occurrence.
Energies 2020, 13, x FOR PEER REVIEW 14 of 18 Eisenberg [33], the turbine's power output does not experience any significant variation for wind 422 speeds above rated when the blades are damaged and therefore the higher the mean wind speed, the 423 lower the variation in AEP. These results clearly depend on the IEC class that was chosen. Lowering 424 the average wind speeds even further (IEC Class III), the turbine is expected to spend less time at 425 rated power, and therefore, AEP losses are expected to further decrease for the present test case. 426 Although low wind speed sites have recently been exploited for wind turbine installation, specially 427 designed machines with low specific power are being installed in such sites, resulting in machines 428 that are able to spend significant time at rated power even in these sites. As noted in [33], a utility-429 scale machine will spend 40% to 60% of its time at rated power, where blade damage has no effects. 430 In addition, although the main cause of LE erosion is related to the rotational tip velocity, it can be 431 argued that in lower wind speed sites, less transport of abrasive particles will arise, therefore leading 432 to less erosion. 433    The wind class IIA shows higher standard deviation and higher left tail length. As previously mentioned, which is due to the fact that in the class IIA scenario the turbine operates at rated power for a shorter period of time with respect to the class IA scenario. In fact, as also pointed out by Eisenberg [33], the turbine's power output does not experience any significant variation for wind speeds above rated when the blades are damaged and therefore the higher the mean wind speed, the lower the variation in AEP. These results clearly depend on the IEC class that was chosen. Lowering the average wind speeds even further (IEC Class III), the turbine is expected to spend less time at rated power, and therefore, AEP losses are expected to further decrease for the present test case. Although low wind speed sites have recently been exploited for wind turbine installation, specially designed machines with low specific power are being installed in such sites, resulting in machines that are able to spend significant time at rated power even in these sites. As noted in [33], a utility-scale machine will spend 40% to 60% of its time at rated power, where blade damage has no effects. In addition, although the main cause of LE erosion is related to the rotational tip velocity, it can be argued that in lower wind speed sites, less transport of abrasive particles will arise, therefore leading to less erosion.

Conclusions
This study proposes the use of an uncertainty quantification approach to the modelling of the effects of blade damage on the performance of multi-MW wind turbines. The proposed approach aims at overcoming some of the issues associated with the evaluation of a single test case. In fact, treating blade damage as a random phenomenon, bias due to a specific test case of a combination of bladedamaging factors can be avoided and more general conclusions can be drawn. The entire process is simulated numerically. First, geometric shape modifications are applied to the airfoils that constitute the turbine's blade. Lift and drag coefficients are calculated using CFD. The newly found coefficients are then applied to an aero-servo-elastic model of the wind turbine. Uncertainties are propagated through the model using an arbitrary polynomial chaos method.
Results show that LE damage has the larger influence on power and AEP losses. For the selected test case, TE damage has little impact, except for when the turbine is operating at very low wind speeds, where a slight performance increase is observed due to TE damage. Focusing on absolute values, maximum average power reductions are observed at 8 m/s and 10 m/s mean wind speeds and are of 2.2% and 2.6%, respectively. The most unfavorable damage combinations simulated showed a decrease in AEP of up to over 6%. By looking at the probabilistic framework, however, the configurations with the highest probability of occurring based on the input PDFs show AEP reductions of below 1% in both IEC classes I and IIA. Indeed, mean AEP reductions of just below 1% for class IA and just above 1% for class IIA are estimated. These values, whilst significant, seem to be notably lower than what is commonly forecasted in published literature that, however, is strictly site-or turbine-dependent. It is important to point out that the results of the present study do not indicate that published literature values are unrealistic (even though sometimes a too large span coverage of erosion is considered), however, for the present test case, representative of modern turbine size and design trends, such values seem to have very low probabilities of occurrence. Indeed, AEP decreases exceeding 10% are noted in the present study. Blade damage is an issue that should still be taken very seriously by the industry, due to its structural implications that were not investigated in the present work; however, the impact on AEP does not seem to be as pronounced as early research indicated. A great deal of factors could cause these discrepancies, which could be due to the radial damage extension considered and size and hence the Reynolds number of the turbine, which are not investigated herein and therefore remain an open issue, where additional research would definitely be beneficial. As pointed out by other authors, LE damage seems to have a lower effect on larger wind turbines. Although this is not the focus of the present work, the results of this study, if put into perspective with other published literature that reports higher AEP decreases on smaller turbines, seem to confirm this.
Moreover, as previously discussed, the results strongly depend on the input PDFs. The presented method can be however adapted to different input PDFs, which are hopefully more extensively supported by field data. Nevertheless, the present assumptions can be considered realistic for medium-low damaging environments or for blades where regular maintenance schedules are planned. It is also important to point out that these results are valid strictly only for the present test case. A selection of different study cases might influence the results significantly, as, in the authors' experience, LE damage affects different airfoils to different degrees. Finally, the LE damage model also influences the results. Although the model is calibrated and tested with respect to experimental data and is adequate for the present parametric framework, it is hard, if not impossible, to accurately reproduce small, stochastic features that might influence the sectional efficiency significantly.
In conclusion, even considering these factors, it is apparent that the present statistical approach is able to give designers a better picture of the impact of blade damage. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflicts of interest.