Fluctuations of angle of attack and lift coefficient and the resultant fatigue loads for a large Horizontal Axis Wind turbine

Unsteady loads are a major limiting factor for further upscaling of HAWTs considering the high costs associated to strict structural requirements. Alleviation of these unsteady loads on HAWT blades, e.g. using active ﬂ ow control (AFC), is of high importance. In order to devise effective AFC methods, the unsteady loading sources need to be identi ﬁ ed and their relative contribution to the load ﬂ uctuations experienced by blades needs to be quanti ﬁ ed. The current study investigates the effects of various at- mospheric and operational parameters on the ﬂ uctuations of a and C L for a large HAWT. The investigated parameters include turbulence, wind shear, yawed in ﬂ ow, tower shadow, gravity and rotational imbal- ances. The study uses the DTU's aeroelastic software HAWC2 . The study identi ﬁ es the individual and the aggregate effect of each source on the aforementioned ﬂ uctuations in order to distinguish the major contributing factors to unsteady loading. The quanti ﬁ cation of contribution of each source on the total fatigue loads reveals > 65% of ﬂ apwise fatigue loads is a result of turbulence while gravity results in > 80% of edgewise fatigue loads. The extensive parametric study shows that the standard deviation of C L is 0.25. The results support to design active load control systems by highlighting the magnitude of C L and a variations experienced by HAWTs, and thus the dC L that needs to be delivered by an AFC system. of of unsteadiness on uctuations


Introduction
Previous studies have shown that increasing the size of wind turbines is a driving factor towards more environmentally friendly and cost-efficient wind energy [1,2]. However, as wind turbines become larger, the rotating blades subsequently become longer. Longer blades necessitate the use of more slender and flexible designs, which lead to larger tip deflection in steady state, near rated wind speed. Larger deflection for the slender blades may be also a result of structural mode coupling due to unsteady, spatially inhomogeneous, wind fields and loads. This imposes strict requirements, leading to thicker blade laminates with higher structural strength and stiffness in order to withstand the inherently higher loads, due to the effects of gravity and non-uniform inflow conditions, and to limit the resultant larger deflections, ultimately leading to more expensive blades.
In order to refrain these undesirable high costs, the sources of unsteady loads on the blades need to be identified and their resulting fatigue damage needs to be determined. Furthermore, the resultant load fluctuations of each source need to be quantified. Identification of the load fluctuations helps to effectively design the active load control systems in order to alleviate the identified fluctuations [3]: the dC L that needs to be delivered by an AFC system will be known.
The atmospheric boundary layer imposes a range of operating conditions on wind turbine blades during their lifetime. These include unsteady fluctuations in wind speed and direction as well as gradients of mean velocity in both vertical and lateral directions. Previous studies have shown atmospheric characteristics such as turbulence [4e6] and stability [7e9] are important for wind turbine unsteady loading. Additionally other effects including wind shear [10e12], yawed inflow [13], gravity, tower shadow, mass and aerodynamic imbalances and wake effects [14e16] also have a significant contribution to load fluctuation experienced by HAWT rotor blades.
Unsteady loading can lead to structural resonance and fatigue damage and finally structural failure and reduction of lifetime for wind turbine blades [17]. Part of the unsteady loading are caused by fluctuations in the angle of attack and consequently fluctuations of forces and moments on blades. Therefore, quantification of the fluctuations in angle of attack and lift coefficient under various loading conditions is important for wind turbine blades.
Though previous efforts [18,19] mainly focused on the investigation of mean value of a and C L , few studies addressed the fluctuations of a and C L . Moreover, such studies either investigated the effect of fluctuation of a on the turbine output power [20] or only studied the aggregate effect of all sources of unsteadiness on the fluctuations of C L and fatigue damage [21]. Therefore, the quantification of the individual effect of unsteady loading sources on the fluctuations of a and C L and the associated fatigue damage has received minimal attention to this point. The current study intends to investigate the influence of various unsteady loading cases on the fluctuations of a and C L for a large HAWT in a systematic and holistic approach. Additionally, a comparison of lifetime equivalent fatigue loads of the blade root bending moment corresponding to each source of unsteadiness is performed.
As such, the present study can support the identification of the major contributors to unsteady loads on wind turbine blades.
Moreover, quantification of the fluctuations in a and C L can provide guidelines for the design of active load control systems for wind turbine blades [3,22,23].
The current study uses the BEM-based aeroelastic wind turbine design code 'HAWC2' from DTU Wind Energy. Although highfidelity CFD calculation are also commonly employed to study the wind turbines [24e26], aeroelastic codes offer the advantage of being computationally much cheaper while offering a satisfactory level of agreement with higher-fidelity modeling and experiments (see subsection 2.1). This grants us with the opportunity to study a wide range of parameters at a reasonable computational cost.
The paper starts with a methodology section where the investigated wind turbine and the employed aeroelastic code are described in subsection 2.1. The simulations settings and test matrix is presented in subsection 2.2 and the employed method for calculation of fatigue loads are discussed in subsection 2.3. The results and discussion are presented in section 3 where the fluctuations of angle of attack and lift coefficient are discussed in subsection 3.1 and the resultant fatigue loads are elaborated in subsection 3.2. The conclusions are provided in section 4.

Wind turbine and aeroelastic code
The 10 MW reference wind turbine (DTU-10MW-RWT) developed by DTU is employed in the current study. DTU-10 MW-RWT is a 3-bladed upwind horizontal-axis variable-speed pitch-regulated yaw-controlled wind turbine with schematic shown in Fig. 1. The turbine is selected as a representative of future trend of large wind turbines. Geometrical and operational characteristics of this turbine are described in Table 1. Further details can be found in Ref. [27]. The geometrical and operational characteristics of this turbine (aerodynamics, structure and controller) are used as inputs in the simulations. The simulations are performed using HAWC2 (Horizontal Axis Wind turbine simulation Code 2nd generation) software. HAWC2 is an aero-servo-hydro-elastic wind turbine design core developed by DTU Wind Energy at Risø Campus in Roskilde, Denmark 'www. hawc2.dk' [28e30]. The software includes 3 modules to account for the wind, the structure and the aerodynamics in order to simulate the turbine while a controller module is included in the model of the specific turbine which will be used in any simulation.

Nomenclature
The wind module includes two components: 1) a deterministic component in order to set the average wind speed with wind steps and the wind shear; 2) a stochastic component uses the Mann turbulence model [31] in order to generate a full 3D anisotropic turbulent inflow wind characteristics. The selected model is based on an algorithm to generate 2D/3D turbulent atmospheric boundary layer data for different wind speeds and turbulence intensities. The model is commonly used in wind energy and wind engineering applications. It is important to note that in BEM-based models, the effect of turbulence on instantaneous blade section aerodynamic loads is indirectly modeled. This means that when the velocity fluctuates, a new induced velocity is calculated which results in a modified angle of attack. Subsequently, the airfoil lift and drag forces for the new angle of attack are 'updated' from the airfoil polars provided to the code. Potential flow solution is used in order to account for the effect of the tower shadow.
Structural module is based on a multi-body formulation of the turbine [32]. More specifically, every body (e.g. one blade) is simulated as a set of Timoshenko beam elements (6 degrees of freedom) and every HAWT component (e.g. rotor, tower, nacelle) has its own reference frame and consists of one or more bodies. The coupling between the bodies and components is based on algebraic constraints (e.g. fixed relative position, joints, etc.) which accounted for large rotation and translation in the coupling point [30].
Aerodynamic module is based on a blade element momentum (BEM) theory [33]. BEM is a theory that combines momentum theory with blade-element-theory in order to calculate the loads on the blades of a turbine/propeller. BEM in HAWC2, however, is extended, compared to the classical BEM, with the following features (also see Table 2): Highly-loaded rotor: erroneous predictions of BEM at high induction values are corrected based on Glauert's empirical correction [33]. Tip loss: in order to account for the finite number of blades, the corrections by Prandtl's [34] are included.
Non-uniform induction: in order to account for non-uniform induction over each annular blade element, the rotor disc is discretized in a stationary polar grid so that variations of induced velocities over the span and azimuth can be captured. This improves the results for sheared, skewed and turbulent flows [35]. Dynamic inflow: this model allow to take into account the significant dynamics of transient rotor wake after a sudden change in the blade loading, e.g. change in pitch or yaw angle [36,37]. Skewed inflow: two considerations are included in order to account for the skewed inflow (for a yawed or tilted rotor): 1) the reduction in the mean induction as a function of skewness angle which is done using the Glauert method [33] 2) the azimuthal variation of induction which is caused by the skewed wake which is done based on the method proposed by Coleman [38]. Large blade deflections: in order to account for this, two considerations are required: 1) the blade forces will not be normal to the rotor plane; 2) the effective diameter of the turbine will decrease. Corrections are implemented to account for this. Dynamic stall: the oscillations in angle of attack will result in unsteady lift and drag forces on the blades to be different than the static values. The corrections are implemented using a modified Beddoes-Leishmann model for wind turbines [39e41]), and thus capture the oscillations due to the unsteady inflow and/or fluctuations due to turbulence affecting the unsteady blade section loads. It is noted that recent studies [42,43] have suggested that a more accurate estimation of unsteady loads due to dynamic stall should consider synchronous oscillations in angle of attack and inflow velocity, even though it could be also be argued that quite complex dynamic stall models already exist [44] and that future efforts should address the coupling of the different aerodynamic models instead [45]. Nacelle and tower drag: aerodynamic drag for the tower and nacelle are calculated based on solution of potential flow. The solution provides information regarding the stagnation region in front of the tower and nacelle and consequently the respective velocity field is updated. The employment of the modified Beddoes-Leishmann model ensures the solution includes the effect of this sudden change in the velocity field in front of the two bluff bodies on the unsteady loads on blades. 3D rotational effects on airfoil polars: the 2D airfoils polars for DTU-10MW-RWT input into HAWC2 are corrected for 3D rotational effects [46]. Validation studies have been performed for HAWC2 against higher-order models and/or experiments in order to ensure the accuracy of the findings. The studied cases include free inflow and even more challenging inflow situation of partial wake of upstream turbine where good agreement has been observed in all studied cases. Some of the studies are briefly listed below in chronological order: Blade root bending moments and flow angle on blades for a 2 MW wind turbine operating at the Tjaereborg wind farm in Table 1 Geometrical and operational characteristics of DTU-10MW-RWT [27]. Denmark were compared against experimental data by Larsen et al. [47] in 3 different inflow conditions with 3% TI: free inflow, 1=3 and 2=3 wake of upstream turbine. Blade root bending moments and load distribution along span on blades for NREL-5MW-RWT in half-wake of an upstream turbine were compared against 3D CFD results by Barlas et al. [48]. Load and axial induction distribution along span, power and thrust force for NREL-5MW-RWT in high-shear atmospheric boundary layer were compared against a free-wake vortex code and 3D CFD by Madsen et al. [35]. Blade root bending moments, power and thrust force for NREL-5MW-RWT at different wind speeds are compared against 3D CFD by Heinz [49]. Tangential and normal forces along the blade span for NREL-5MW-RWT in the half-wake of an upstream turbine at 10 m=s with 10% TI were compared against 3D CFD results by Larsen et al. [50] (see Fig. 2).

Parameter
The blade root bending moments, deflection distribution along span, thrust force and power versus wind speed for DTU-10MW-RWT were compared against a free-wake vortex code in Ref. [51] during the 'INNWIND.EU' project.
Several other validations studies are performed in different inflow conditions by Larsen et al. [52], Barlas et al. [53], Ostachowicz et al. [54] and Pirrung et al. [55] where HAWC2 results were compared against higher-order models and/or experiments and satisfactory agreement was observed.

Simulation settings and test matrix
In order to identify the effects of sources of unsteadiness on fluctuations of a and C L , fatigue damage corresponding to the blade root bending moments, a series of 34 aeroelastic simulation cases were carried out. Each simulation scenario consisted of 11 individual cases with wind speed ranging from 4 À 24m=s (bin width of 2m=s). Upper and lower limits of this range correspond to cut-in and cut-out velocities, respectively. The wind shear is defined using a 'power law' wind profile with an exponent (a ¼ 0:2) according to IEC standard [56]. Equation (1) shows the expression for the wind profile where U, U t , z and z t are mean wind speed, mean wind speed at the reference height, height and the reference height, respectively. Table 3 presents the structural, controller and wind inputs to the simulations. Simulation and sampling periods were as prescribed by the design requirements for wind turbines IEC À 61400 À 1 edition-3 standard [56] (see Table 4).
The simulation matrix describing the 34 tested cases is shown in Table 5. The considered parameters are tower shadow (TS), wind shear (SH), turbulence intensity (TI) and yawed inflow (Y). The values for the yaw, wind shear and turbulence show the yaw angle ( ), wind shear power law exponential and turbulence intensity (%) respectively. The test cases are selected to enable systematic analysis of each of the aforementioned parameters. Case-01 considers only the effect of loads due to gravity, mass and aerodynamic imbalances. A star sign ('Ã') is used to show the presence of tower shadow in columns. The 'y' sign shows case-30 as the reference case Fig. 2. Normal force along the blade span calculated using HAWC2 and 3D CFD (both performed by Larsen et al. [50]) for NREL-5MW-RWT in the half-wake of an upstream turbine (on the right side of the rotor disc) at 10 m=s with 10% TI with blade at (a) 270 + : out of the wake; (b) 90 + : in the wake. for all the comparisons. Case-30 corresponds to inputs according to a Class-A wind turbine [56]. The distribution of turbulence intensity versus wind speed was calculated based on reference turbulence intensity (I ref ), mean wind speed (U) and turbine hub height (Equation (2)) as prescribed by Ref. [56]. I ref is the expected value of hub-height turbulence intensity at a mean wind speed of 15 m=s averaged over 10 min.

Calculation of fatigue loads
Fatigue equivalent loads for different wind speeds were calculated using rainflow counting algorithm and Palmgren-Miner's summation rule using Equation (3) [57e60]. F EQ is equivalent fatigue load, n EQ is the number of equivalent cycles, which is 600 for 10-min fatigue cycles with 1 Hz and 10 7 for lifetime fatigue cycles. F i and n i are the fatigue loads ranges and cycles, respectively. m is the W€ ohler exponent characterizing the material properties and is set to 10 in the present study [60]. The lifetime equivalent fatigue loads (LEFL) can be calculated using Equation (4) [60] where F EQ ÀTOT is LEFL, F EQ ðUÞ is EFL as a function of mean wind speed calculated using n EQ ¼ 10 7 . The value of n EQ is the typical value commonly used as the number of load cycles that a wind turbine experiences during a lifetime of 20 years. PðUÞ is the wind speed probability density function (PDF) described later in this section.
The fatigue damage (D) at each wind speed was calculated by dividing the fatigue equivalent loads for 1 Hz fatigue cycles at the specified wind speed by the sum of fatigue equivalent loads at the whole range of wind speeds from cut-in (U in ) to cut-out (U out ) for 1 Hz fatigue cycles (Equation (5)). This does not include the wind speed probability density function. The effect of wind speed probability density function was introduced using Equation (6) and normalized using the sum from cut-in (U in ) to cut-out (U out ) using Equation (7).
In order to calculate the lifetime equivalent fatigue loads (LEFL) of the blade root bending moments, the probability density function of wind speed is required. This is determined employing the Weibull distribution to describe the wind climate. The turbulence seeds at each wind speed was calculated based on the IEC standard [56]. The employed Weibull distribution has a shape factor (k) of 2.03 and a scale factor (A) of 11.9 (see Fig. 3). The Weibull distribution was fitted to 23 years of wind data with one measurement per hour and a return rate of 91:5% taken from the off-shore 10 mheight meteorological mast at station 321 Europlatform provided by Koninklijk Nederlands Meteorologisch Instituut (KNMI) in the Netherlands. The wind data was scaled up to the DTU-10MW-RWT hub height (119 m) using logarithmic law. The station is located 40 km away from the Dutch coast at latitude 51 + 59 0 56 00 north and longitude 3 + 16 0 34 00 east. The availability of wind data for such a long period of 23 years for this offshore location ensures that the employed wind data is a good representation of the wind climate of a potential offshore installation location for such a large HAWT. The load time series consisted of 11 time series covering the cutin to cut-out velocities (4 À 24 m=s with steps of 2 m=s). I ref was the respective value for a class A wind turbine, i.e. 0.16.

Fluctuations of angle of attack and lift coefficient
The current parametric study was conducted based on the The '*' sign shows that tower shadow is considered for that case.
results of the HAWC2 simulations as listed in Table 5 to elucidate the effect of sources of unsteadiness on fluctuations in a and C L in the rotor plane for various wind speeds. The study has excluded the root sections with r R < 0:3 as these inboard stations have a negligible contribution to aerodynamic loads felt at the blade root because of their somewhat cylindrical cross section and small effective velocity, and as such are not shown. Additionally, as shown by CFD results for such large HAWTs such as NREL-5MW-RWT and DTU-10MW-RWT, e.g. Lin et al. [61], Bak et al. [62] and Ostachowicz et al. [54], the root cross-flow effects are very much limited to r=R < ¼ 0:3, therefore their impact on radial sections with r=R > 0:3 which are investigated in the present study are expected to be negligible.
As previously mentioned, the wind turbine used in the simulations is the DTU-10MW-RWT [27]. At this point, it is necessary to introduce the various airfoils used in this turbine and their design characteristics in order to better understand the operating conditions based on a, for different radial positions (r). Therefore, the relevant data are provided in Table 6. The distribution of the design a and C L for the airfoils of the DTU-10MW-RWT over the span is illustrated in Fig. 4.
Prior to the outcomes of the parametric study, the reference case (see Table 5) is examined in order to identify the wind speed where the highest a is experienced for the largest portion of the blade.
Variations of a versus span shown in Fig. 5 illustrates that below rated wind speed, 11.4 m=s, the average a slightly increases from inboard to outboard stations while this reverses after rated wind speed when the turbine starts to pitch. This is in agreement with the results from various modeling methods; including panel method, lifting line and lifting surface; from Ref. [19]. Variations of a versus wind speed shown in Fig. 6 shows the maximum a is experienced at velocities slightly less than the rated wind speed (before the turbine starts pitching). The highest a was experienced at 8 m=s, then 10 m=s. For all the inboard and outboard sections the value rapidly decreases at lower wind speeds. The rate of reduction of mean value of a (from maximum value at 8 m=s) increases from inboard to outboard stations. This is in agreement with results of    [21] where a similar trend was found. This preliminary study has been conducted for the other unsteady loading cases (from Table 5) and similar results were observed. Therefore, the focus of the study has been set at 8 m=s where the highest a was experienced on the blade. This is because the blade sections are approaching stall, and hence dynamic stall and hysteresis effects become more pronounced and may lead to larger amplitude loading [37,63].
The mean value of a and C L at 8 m=s with fluctuations (standard deviation) shown as shaded area are shown in Fig. 7 versus span for the reference case. It is seen that the mean values have a small drop in the midspan with slightly higher values at the tip. This is in agreement with the results of [19], and consistent with the design a spanwise distribution shown in Table 6.
The result shows that fluctuations of a and C L are higher near the root and lower near the tip. This is simply due to lower rotational velocity near the root. Accordingly the same fluctuation in wind speed will lead to larger fluctuations for a and C L at the root than the outboard region. Similar trend was observed by Refs. [21,64]. The fluctuations of a start from ±4 near the root and decrease to ±2 near the tip. This corresponds to DC L ¼ ±0:4 near the root and DC L ¼ ±0:25 near the tip. In order to compare the operating a and C L versus the corresponding airfoils polars for the local blade section and to identify the sections that experience stall during operation at 8m=s (highest a), a values are plotted against the polars in Fig. 8. This shows that sections from midspan to tip hardly experience stall during operation while sections with r R < 0:35 happen to operate in stall mode for about 1% of their operational time. The probability of occurrence of stall along the blade span for different load cases is illustrated in Fig. 9a. The results show that, for the studied turbine, increasing I ref from 8% to 20% (case-07 to 10) gradually increases the stall probability at the inboard stations from approximately 0 up to 2:2% while the effect diminishes towards the midspan. Insignificant differences are observed for outboard stations. In addition, comparing case-09 and 30, tower shadow and wind shear exhibit almost no influence on the stall probability. Moreover, comparing case-30 and 32 where the only difference is 10 + yaw angle for the latter case, yawed inflow is also found to have a negligible effect of the stall probability over the span. These findings are in agreement with previous research from Ref. [21].
Instantaneous values of C L versus azimuth for r ¼ 80 m is plotted in Fig. 9b to give an overview of the C L in each blade rotation. Similar behaviour is also observed by Ref. [20]. The peaks at about 0 + or 360 + in the figure correspond to where blade is at its vertically highest position and the trough at 180 + corresponds to where the blade is at its lowest vertical position parallel to the tower height. The near sinusoidal variation is a result of the wind shear that results in a gradual increase of wind speed with the altitude of the blade. Consequently the blade experiences higher a when it reaches the highest vertical position and it gradually decreases when it comes to the lowest vertical position.
The mean values a and C L for various unsteady loading cases were compared (see Fig. 10). Case 01 from Table 5 was defined as having only gravity and imbalances as the sources of unsteadiness. Since this was computed with a (space and time) homogeneous wind field, the remaining operational cases, including tower shadow, wind shear, yaw and turbulence, were compared to identify the individual effects. It is seen that almost all the sources of unsteadiness reduce the mean value of a and C L along the span.
Moreover, the effect is larger as multiple sources are considered, compared to single source cases.
In contrast, the scenario is totally different for the fluctuations of a and C L (see Fig. 11) where all the sources of unsteadiness increase the fluctuations of a to various amounts. The values are normalized with the relevant values for case-01 from Table 5 to make the increments more self-explaining. The normalized values are shown in Fig. 12. The most important conclusions stemming from this comparison are: Tower shadow has negligible effect on the fluctuations. Yaw angle (of 10 + ) increases the fluctuations by a factor of 2 almost uniformly along the span. Wind shear increases the fluctuations by a factor of 2 near the root. This factor increases almost linearly to 4 for r R ¼ 0:3 À 0:7 and stays constant outwards to tip. Turbulence (I ref ¼ 0:16) will result in a significant increase in fluctuations. The level of increase starts from a factor of 10 near the root, increases to a factor of 12 around the midspan (where the mean values are the lowest) and decreases to a factor of 8 at the blade tip.
A graphical comparison of the effect of various sources of  representative value for a Class-A HAWT. This finding can support an estimation of the effect that a potential active load control system will be required to counteract, although the magnitude of required DC L also depends on some other parameters, e.g. the spanwise extent of actuation employment.
Instantaneous a for a single rotation of the blades is plotted over the rotor plane (see Fig. 14) for the reference case, representative of a Class-A wind turbine. Though wind shear results in higher a in the upper half plane, the influence of turbulent eddies of various sizes is observed as regions of low and high a in the rotor plane. Therefore, turbulence was found to be the main reason for large DC L and unsteadiness in loads over the entire blade, consistent with Fig. 12. Additionally, turbulence results in local unsteadinesses over the span. Thus, load control systems acting locally and distributed radially on the blade might be very successful. This outcome can suggest distributed flow control [65] methods as a promising candidate for load control on wind turbine blades.
A similar study was done for various unsteady loading cases and the following key points are inferred from the results: Tower shadow results in a region of abrupt low a as a result of the stagnation region.
Yaw results in a combined vertical and horizontal gradient in a. The horizontal gradient in a can be attributed to the skewed wake effect [66] which results in variations of streamwise velocity (and alpha) where larger magnitude of variation is on horizontal (left and right) azimuthal positions. While the vertical gradient in a can be caused by the advance-retreating blade effect [13,37,67,68] which results in variations of tangential or in-plane velocity (and alpha) where larger magnitude of variation is on vertical (top and bottom) azimuthal positions.
As expected, wind shear results in a vertical gradient in a where approximately 2 + higher a is experienced in the upper half. Local fluctuations in a were observed due to turbulent eddies of different sizes passing the rotor plane while the blade is rotating as well (rotational sampling).
A combined effect is observed (see Fig. 14) where all unsteady sources are considered, as representatives of Class-A wind turbine.

Turbulence and wind shear
An investigation of the effect of turbulence intensity and wind shear on LEFL for the blade root bending moments was also carried out, for flapwise (Mx), edgewise (My) and torsional (Mz) moments. Cases 1,6e10 and 4,15e19 and 24,27e31 from Table 5 showed (see     15) an almost linear increase in normalized LEFL with increasing turbulence intensity for flapwise, edgewise and torsional moments, compared to the reference case. Moreover, wind shear was found to have a significant effect on flapwise, edgewise and torsional LEFL. Wind shear resulted in an almost constant increase in edgewise and torsional LEFL at different turbulence intensities, while the effect of wind shear on the flapwise LEFL was smaller at higher turbulence intensities.
Additionally, when I ref increased from 0.00 to 0.16 (case-01 and case-09 from Table 5), the normalized flapwise LEFL increased from 0.34 to 1.00. The comparison implies that within the studied range, approximately 66% of flapwise LEFL is caused by turbulence and the other 34% is a result of combination of gravity, mass and aerodynamic imbalance, tower shadow and wind shear. The increase in the normalized edgewise LEFL with increasing turbulence (to I r ef ¼ 16%) is about 14%. This highlights the lower impact of wind conditions for the edgewise LEFL and the greater impact of gravity and imbalances which contributed to over 86% of the edgewise LEFL.
The effect of turbulence on fatigue damage is found to be very significant. For case-1, a constant distribution of fatigue damage versus wind speed was found. However, additions of turbulence (cases 6e9) change this behavior and result in more fatigue damage at higher wind speeds, which is consistent with previous studies [69]. This effect is found to be larger for flapwise than edgewise blade root bending moments. Additionally, Fig. 15 shows the presence of turbulence results in a significant jump in magnitude of flapwise LEFL.
Fatigue damage distribution for the reference case (a representative of a Class-A wind turbine) versus wind speed and a comparison of LEFL values for flapwise (Mx), edgewise (My) and torsional (Mz) moments are shown in Fig. 16 and Fig. 17. The magnitude of flapwise LEFL is found to be primarily a result of turbulence. Secondary effect is due to gravity and imbalances and tertiary due to wind shear (see Fig. 18a). The case was different for the edgewise LEFL where the magnitude of LEFL is mainly caused by gravity and imbalances with turbulence as the second dominant source (Fig. 18b).
A similar behavior is observed for the fatigue damage resulting from all the individual blade root bending moments versus wind speed, which is a gradual increase in fatigue damage from cut-in velocity to cut-out velocity. However, when the effect of the wind speed PDF was introduced into fatigue damage, the highest fatigue damage for all the blade root bending moments was found near the rated wind speed. This is in agreement with results of [7,8,21,70] where they found similar behavior for blade root bending moments and the corresponding fatigue damage.

Overall contribution
A comparison of the contribution of sources of load unsteadiness for flapwise (Mx) and edgewise (My) blade root bending   moments on LEFL was carried out (see Fig. 18). Fig. 18a elucidated the main contribution of turbulence for flapwise LEFL while gravity, imbalances and wind shear were the other prominent effects. The flapwise LEFL caused by aerodynamic imbalance (pitch and twist off-set of the blades) increases by 8% when I ref increases to 0.16. Fig. 18b shows gravitational loads are the main contribution for edgewise LEFL while the weaker contributions are from turbulence and wind shear. The main conclusions are stated below: Gravity and imbalances result in near 18% of flapwise LEFL but more than 80% of edgewise LEFL Tower shadow has a negligible effect for both LEFL with less than 0:5% contribution Wind shear results in approximately 8% of flapwise LEFL and 4% of edgewise LEFL Yaw (±10 + ) results in near ±2:5 À 5% of flapwise LEFL and less than ±1 À 2% of edgewise LEFL.   The conclusions stated above imply that alleviation of turbulence-induced local unsteadiness with active load control systems can directly decrease the flapwise LEFL. Lower flapwise LEFL may allow for designing lighter HAWT blades and could lead to a decrease in the edgewsie LEFL as a result of decreasing the blade mass.

Conclusion
This study investigated the effect of different sources of unsteadiness (i.e. gravity, imbalances, tower shadow, wind shear, yaw and turbulence) on loads and moments of HAWT blades, particularly on the fluctuations of angle of attack and lift coefficient. Quantification of these values can support the design of relevant active load control mechanisms to alleviate unsteady loading. Moreover, the corresponding lifetime equivalent fatigue loads were compared. The systematic study is based on 34 aeroelastic simulations of the DTU-10MW-RWT using DTU's HAWC2 software. The main conclusions are: The standard deviation of a and C L are almost uniform along the blade with slightly higher values near the root in the presence of turbulence.
Turbulence dramatically increases the standard deviation of a and C L . For Class-A wind turbine with reference turbulence intensity of 0.16, standard deviations are 10 À 14 times compared to the case with no turbulence.
Wind shear and yaw increase the standard deviation of a and C L for 2 À 4 times, respectively. Tower shadow has a negligible effect. For a Class-A wind turbine, the standard deviation of C L are in the range of DC L ¼ À0:25 to þ0:25. This finding supports the design of active load control devices for wind turbines by clarifying the dC L that needs to be delivered by an AFC system. Fluctuations were found to be mainly local due to turbulent eddies, which indicates high effectiveness for active load control systems distributed along the span. More than 65% of flapwise fatigue loads is caused by turbulence. Gravity and wind shear are the second and third main contributors to flapwise fatigue loads. More than 80% of edgewise fatigue loads is caused by gravitational loads. Turbulence and wind shear are the second and third main contributors to edgewise fatigue loads.