*2.3. 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 θ.

**Figure 2.** ε-θ correlation in literature [5,21] and proposed correlation in black point-dashed line.

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 *Dend* = θ/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].

**Figure 3.** Damage modeling. (**a**) Leading-edge (LE) damage; (**b**) TE damage.

## *2.4. 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.

**Figure 4.** CFD validation: (**a**) Illustration of the adopted computational domain; (**b**)validation of the numerical setup in respect to data from [5] at Re <sup>=</sup> 1.5 <sup>×</sup> <sup>10</sup>6.

An unstructured triangularmeshis used. The airfoil's boundarylayerismodelled 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 <sup>×</sup> 105 elements and 500 elements along the airfoil surface, a medium mesh with 2.8 <sup>×</sup> 105 elements and 650 elements along the airfoil's surface, and a fine mesh with 3.6 <sup>×</sup> 105 elements and 750 elements along the airfoil's surface. The lift (Cl) and drag (Cd) 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 <sup>×</sup> 106. 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.

## *2.5. 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].
