Abstract
Wind energy in icing and low-temperature climate has a huge growth potential, but rotor icing effects on turbine dynamics and lifetime are not well known and simulations with iced rotor are not required in current IEC 61400-1 turbine design standard. In this article, simulations with iced rotor are compared to measured mechanical loads. The dynamic behaviour of the wind turbine was simulated with FLEX5 aeroelastic code for Senvion MM92 2 MW wind turbine. Simulations with typical iced airfoil lift and drag coefficients, aerodynamic and mass imbalances for iced rotor were performed and compared to measured iced turbine loads. Resulting iced turbine simulation parameters can be used in defining new design load cases for cold climate turbines. The most representative simulation parameter combination was achieved with a symmetric aerodynamic penalty applied on all blades and an asymmetric rotor mass imbalance of 166 kg ice load on two blades and 83 kg ice load on one blade.
Introduction
Wind power capacity in cold climate, defined by IEA Wind Task 19 Group (2011) as weather conditions of atmospheric icing and low-ambient temperatures which expose wind turbines to conditions outside their normal design limits, is growing rapidly. By 2013, 69 GW (24% of global wind power capacity) were installed and another 50 GW of wind energy will be built in cold climates by 2017 (Navigant Research, 2013). Building wind power to cold climate sites is attractive because of favourable wind conditions and low population density (Lehtomäki and Wallenius, 2015; Navigant Research, 2013). However, icing of wind turbine blades remains one of the main challenges for cold climate sites. Icing of blades has become more significant while turbine size has increased. The most common way of ice accretion is via in-cloud icing. The risk of in-cloud icing increases drastically with elevation and therefore icing risk of blades is much larger for a modern 3 MW turbine with 120 m hub height compared to old turbines with less than 100 m hub heights (Rissanen and Lehtomäki, 2015).
Ice accretion on blades causes production losses, risk of ice shedding and vibrations which can be harmful for turbine component lifetime. However, ice on the blades will decrease airfoil lift and increases drag coefficients. Decreased power production usually means decreased rotation speed and mechanical loads. It has been unclear how combination of increased vibrations and decreased power production effect together to the turbine lifetime and how ice accretion should be taken into account in turbine design. Current IEC 61400-1 standard (IEC 61400-1:2005(E), 2005) does not propose specific design load cases for iced turbines.
In the past, several dynamical load simulation studies for iced wind turbines have been made (Battisti, 2005; Frohboese, 2007; Ganander, 2003; Mayr et al., 2006; Rissanen and Uski, 2005; Salam et al., 2015). The biggest challenge in previous studies has been the absence of model validation. Load measurements and high-quality icing observations have not been available from the same site and simulation model validation has been difficult. In some of these studies, icing effects to turbine aerodynamics has been estimated using a general ISO standard for iced structures (ISO 12494:2001, 2001). This ISO standard is designed for stationary structures, not for a rotating wind turbine blade. Mass imbalance has usually been based on a worst-case scenario where one blade is clean and two blades are iced. This assumption might be appropriate for ultimate load analysis, but it is highly conservative for fatigue analysis. Large mass imbalance leads to very high side-to-side tower loads and this was concluded to be the main reason for increased loads in many of the previous studies. Ganander and Rissanen simulated different mass and aerodynamical imbalances and results were compared to 1P component of measured power for a 600 kW turbine (Ganander, 2003; Rissanen and Uski, 2005). The aim of the studies was to find a simple way to simulate ice-induced aerodynamic imbalances. However, the validation method was not accurate enough and the lifetime of components was not calculated. Frohboese used pitch error to model aerodynamic imbalance of iced turbines (Frohboese, 2007). In most of these studies, tower side-to-side loads (perpendicular to wind speed direction) were concluded to be the most sensitive for increased loads for iced turbines.
Iced airfoil aerodynamic properties can be estimated using existing wind tunnel measurements. However, imbalances (aerodynamical and mass) due to uneven ice build-up and ice shedding cannot be estimated in a wind tunnel or with numerical analysis, and therefore, the connection between the amount of ice, imbalances and mechanical loads cannot be made without full-scale field measurements. TechnoCentre Éolien (TCE) has a test site in Rivière-au-Renard, Gaspé, where Senvion performed strain gauge load measurements parallel to extensive weather and icing observations by TCE. (Lehtomäki et al., submitted)
The motivation of this work was to find validated simulation parameters to represent realistic iced wind turbine behaviour under normal power production operation. These parameters can be later used to define design load cases for other cold climate wind turbines. The target was to find aerodynamic properties (lift and drag) for iced airfoils and aerodynamic and mass imbalances for rotor that produce similar effects to turbine fatigue loads and power production compared to measurements. Measured icing cases were not analysed separately in Lehtomäki et al. (submitted). The goal was to find average airfoil properties and rotor imbalance that correspond with typical icing events.
Iced airfoil aerodynamics
There is a limited amount of detailed information available for ice effects on wind-turbine-specific airfoil aerodynamics and associated ice mass effects. Almost all studies examined are performed in icing and/or wind tunnels with typical two-dimensional (2D) aircraft airfoils. In total, eight different references (referred later as icing airfoil references) were used (Addy, 2000; Addy et al., 2003; Bragg et al., 2007; Broeren, 2002, 2010; Busch et al., 2007; Hochart, 2008; Tammelin et al., 1998). The most extensive study made at NASA focused on detailed ice accretion effects of three typical aviation airfoils for various atmospheric conditions in two separate wind tunnels. With this extensive study (and several other individual icing wind tunnel measurement studies), some first statistics for nearly 200 individual icing cases and their effects on airfoil penalties can be summarized. As can be seen from NASA ice accretion tests (Addy, 2000), in-cloud icing conditions were tested in an icing wind tunnel for the following parameter ranges:
Duration (t) = 1–40 min
Flow speed (v) = 67–146 m/s
Liquid water content (LWC) = 0.1–0.66 g/m3
Mean volume diameter (MVD) = 15–56 µm
Temperature (T) = −4°C … −20°C
It should be highlighted that the aircraft industry is typically interested in extreme event analysis for flight safety reasons and not on typical icing events. This means that the parameters used in icing wind tunnel studies are typically more harsh (large water content and droplet size) than typically observed in nature during low level in-cloud icing events (Lehtomäki and Rissanen, 2014). However, icing event duration in these studies have been much shorter compared to typical icing events near ground level. It can be assumed that longer and milder icing events results similar ice shapes. Typically smaller and thinner airfoils are more prone to icing than thicker airfoils (Homola et al., 2009). This might indicate that aerodynamic penalty effects on thick airfoils in the turbine blade root section may be overestimated if assuming similar aerodynamic behaviour as for thinner iced aviation airfoils. On the outer third of wind turbine blades, where most of the power and loads are generated, same type of airfoils are used as in aviation industry.
Iced airfoils were visually divided into four categories (start, light, moderate and extreme icing) to ease comparison of results and derive statistics. Figure 1 presents some examples of the different categories and the effects to the lift coefficient. Figure 1 also represents some typical icing durations and the approximated range of ice mass per kilogram on the airfoil. However, special caution should be observed with the icing durations in Figure 1 as the other icing parameters (often LWC) were typically larger than for typically observed during ‘in nature’ in-cloud icing events. It can be clearly seen from Figure 1 that airfoil lift decreases more and more as more dramatic ice shapes accumulate. Another clear trend is that airfoil penalties are not constant over all angles of attack (α), but that more airfoil penalty is introduced by increasing angle of attack. This means that the iced airfoil reaches stall much earlier than clean airfoils.

Top: ice accretion shapes from various icing categories. Bottom: respective lift coefficient (CL) curves of clean and iced cases. x-axis is the angle of attack and y-axis is the lift coefficient (Addy, 2000).
Based on the four categories, lift and drag coefficients for clean and iced airfoils were extracted from the eight icing airfoil references. The relative change from clean to iced performance was investigated. Iced airfoil penalty is presented as a penalty factor. The penalty factor, thus, is the relationship between clean and iced airfoil performance so that penalty factor = iced/clean and can applied to any clean airfoil lift (CL) and drag (CD) coefficients to represent iced airfoil performance. Used reference cases are listed in Table 1.
Airfoil lift and drag coefficient measurement studies used in average lift and drag calculations.
Figures 2 and 3 present the averages of iced penalty factors per ice category. In Figure 2, the purple dashed line represents the average of all ice categories combined. The continuous black line is a second-order polynomial fit, providing a high correlation coefficient, to the average of all ice categories. The best-fit function with a more modest correlation coefficient for CD was linear due to large scatter of data. Similarly, to the average of all ice categories, best-fit functions were fitted for other categories and results summarized in Table 2. The best-fit functions (continuous black lines in Figures 2 and 3) for the average of all ice categories can be considered as typical airfoil penalty effects in icing conditions for aviation airfoils.

Average iced airfoil lift penalty per ice category and the best-fit function.

Average iced airfoil drag penalty per ice category and the best-fit function.
Iced penalty factor as function of angle of attack (α) for different ice categories.
Model and simulations
The modelled turbine was a Senvion MM92 CCV with 2.05 MW rated power, 80 m hub height and 92 m rotor diameter. It is a cold climate version of MM92 turbine, but there is no blade ice protection system installed in turbine which was used in model validation. Simulations were made with an aeroelastic wind turbine simulation software Flex5 and ‘black box’ turbine model provided by the turbine manufacturer. The ‘black box’ feature enables third parties to simulate turbines with Flex5 code without having access to the source code, turbine controller specifications or detailed information about turbine structural specifications or rotor aerodynamics. Flex5 uses a reduced modal representation of blades and tower and has 28 degrees of freedom (DOFs) in total. For the aerodynamic behaviour, the combined blade element and momentum (BEM) theory is applied. The control parameters, as well as subsystems (e.g. pitch device), are implemented in the Flex5 model. In this analysis, only 21 DOFs were activated and the supervisory controller was deactivated.
In order to find realistic iced turbine behaviour of the simulation model with minimized amount of different simulation cases, realistic values for parameters were defined based on literature and old studies. Table 3 describes the simulation model parameters that were modified during analysis. ‘Parametric study’ analysis method was selected to find out parameter combinations that produces similar fatigue loads and power loss with measurements. The parametric study simulation strategy means that one or two parameters are modified at a time and others are kept constant. Each simulation consists of wind speeds 3–24 m/s with 1 m/s increment (22 wind bins) and 6 turbulent wind fields for each wind bin for IEC turbulence intensity class A. A Weibull–Rayleigh wind distribution and average wind speed of 8.5 m/s are assumed. The simulations are thus consistent with IEC 61400-1 ed3 design load case 1.1 (DLC1.1) resulting to a total of 132 10-min simulations per parametric study case. In total, three wind and air parameters (1–3) and five turbine parameters (4–8) were modified in the parametric simulation study. All possible parameter combinations were not simulated.
List of modified iced turbine simulation model parameters.
Ice mass
For blade ice mass (simulation parameter 4 in Table 3), two different ice mass distribution formulas were used: the ‘industry standard’ GL ice mass formula (Germanischer Lloyd, 2010) and a new linearly increasing ice mass developed by VTT as below
where M is the mass distribution on the rotor blade (kg/m), A is a variable to be verified with turbine load measurements,
ISO 12494 has also been used to evaluate ice mass and distribution on blades. ISO 12494 ice distribution equation is a function of icing duration, object width and wind speed. With this equation, amount of ice can be huge if icing event duration is long. Ice distribution also differs from GL and linear distributions because the equation includes object width. Based on this equation, a thick blade root collects more ice than the blade tip. However, based on numerical analysis, ice mass per length increases towards blade tip (Homola, 2011). Therefore, ISO 12494 ice distribution was not simulated. Figure 4 illustrates the three different blade ice mass models used in the analysis.

Examples of blade ice mass distribution.
For turbine fatigue loads, mass imbalance is more dominant than actual ice mass. Ice mass models were used to define maximum amount of ice on blades and because of uneven ice shedding, worst possible imbalance is when one blade is clean and two blades are fully iced. Also, lesser imbalances with two fully iced and one partially iced blade were simulated.
Analysis methods
Fatigue loads caused by ice were compared to clean, baseline operation using damage equivalent loads (DELs). Because DELs depend on material properties specified by the slope of the S-N curve, 1 Hz DEL was calculated using Miner’s rule. Change in component lifetime was calculated using a standard IEC method (IEC 61400-1:2005(E), 2005) of weighting wind speed–specific loads with wind speed probability distribution which in this case was the Weibull–Rayleigh distribution with a shape factor k = 2 and average wind speed of 8.5 m/s. Wöhler exponent of m = 3 was assumed for steel tower and m = 10 was assumed for blade root. First, ice was assumed to be present on the blades all year round (8760 h) which is highly conservative assumption and for 1 month a year (750 h) which is closer to reality but still on the conservative side compared to 160 h measured instrumental icing in TCE test site (Lehtomäki et al., submitted). Absolute lifetime values were not analysed. Iced turbine results were compared to clean, non-iced baseline results with identical ambient conditions for selected wind, turbulence and wind seed bin. By comparing iced turbine values to clean values (similarly as for iced airfoil performance evaluation), we obtain a relative value indicating the change in performance as follows
With formula (2), it is easy to estimate whether the iced turbine simulation results had a higher (over 1.0) or lower (below 1.0) impact of the turbine behaviour and performance. Relative value of 1.0 indicates normal operation. Relative fatigue loads were calculated with the same method as in Lehtomäki et al. (submitted) for validation data. Coordinate system used in analysis is presented in Figure 5.

Mechanical load measurement coordinate system of test turbine, and in red circles, the used signals for analyses, GL coordinate system (Germanischer Lloyd, 2010) (see colour version of this figure online).
Validation data
TCE owns and operates a R&D test site in Rivière-au-Renard, Gaspé. The wind turbine model was made by the turbine manufacturer Senvion, and the model has been validated with measurement data from TCE site in ice-free conditions. The load measurement camping was from December 2011 to May 2012. Icing events were detected with speed deficit of heated and unheated anemometers. If the 10-min average of a nacelle mounted unheated cup anemometer showed a 20% reduction in wind speed compared to the heated ultrasonic anemometer, the 10-min period was tagged as an icing event. During the load measurement campaign, a total of 10 icing events were recorded and the duration of these icing events varied from 5 to 77 h. Due to the fact that the icing event duration was defined with anemometers, and since ice erosion increases with wind speed (Mazin et al., 2001) causing the rotating blade to recover faster from ice than a stationary object, there is likely ice-free turbine operation included in these icing events. Ice growth is also faster for rotating blades, and it is possible that all icing events were not detected by anemometers. Changes in fatigue loads of different components based on measurements are presented in Table 4 (Lehtomäki et al., submitted).
Measured P90 of relative power and component lifetime during 10 icing cases.
Average measured fatigue loads were not changed much because uniformly decreased aerodynamic lift and increased drag of all blades also decreases fatigue loads. However, periods of increased fatigue loads were also found, especially in tower root side-to-side loads (Mxf). The winter 2011/2012 was relatively mild (Clausen et al., 2013) in Gaspé, and this study was made with only one turbine type. Icing effects to turbine lifetime is turbine manufacturer and model specific and even controller settings have large effects on dynamical behaviour of a turbine during icing events. In extreme conditions, heavily iced wind turbines cannot operate and lifetime of components will not be largely affected. Because of these uncertainties, a conservative analysis approach using the 90th percentile values (P90) of all measured loads during icing events was selected to determine simulation target values.
Maximum estimated ice mass on a blade at TCE site was 192 kg (Lehtomäki et al., submitted). Therefore, moderate ice masses were simulated in addition to huge GL ice mass.
Results
Relative change in lifetime was calculated for every component and every simulation case using method described in chapter 3. A total of seven different signals listed in Table 4 were compared between measurements and simulations and a perfect match was not found. The most representative simulation case was selected by similar behaviour in simulations and measurements and not only matching actual lifetime values. A simulation case with smallest average mean absolute error of all simulated seven signals was not representative because the operation point of turbine was not similar to measurements. The same lifetime values can be received with different rotation speed, but with different icing parameters and production losses. In this case, power loss was smaller than in measurements. This means that imbalances have also been smaller than measured. Therefore, it was important to first match produced power with target values. After the simulated power production matched measurements, similar behaviour was searched from rest of the signals.
A selection of simulation results are presented in Table 5. Simulation case with 166 kg ice mass on two blades and 83 kg in one blade was selected as the most representative case. Other parameters in this case were ‘average of all’ aerodynamic penalty in every blade and 1.3 kg/m3 air density. If selecting criteria would have been lowest mean absolute error of all seven signals, the best case would have been ‘average of all’ aerodynamic penalty in one blade and two blades are ice free (aero-penalty imbalance in one blade). However, as there was only 4% power loss, the simulated turbine operating point was different compared to measurements. In addition, simulations showed an increase in side-to-side (Mxf) fatigue loads compared to fore-aft loads (Myf) and these results were not as clear as in target values. Average aero penalty in two blades produces increased fatigue loads in both tower side-to-side and fore-aft directions, which was also different compared to validation data. Behaviour with only pitch offset without mass imbalance did not match measurements as there was no increase seen in tower root side-to-side (Mxf) loads.
Examples of simulation results.
Table 5 also presents fatigue loads for 50% of GL standard ice mass in two blades. There is a huge 200% increase in tower base side-to-side fatigue loads. Although this may be a realistic design load case for ultimate load analysis, it is not applicable in fatigue analysis. In reality, the turbine controller would stop the turbine quickly due to excess tower vibrations. Table 6 shows how different rotor ice imbalances affect the Mxf ultimate loads. No substantial increases were found in measurements and thus not reported in table below but an inter-comparison between different simulation results is possible. Using the GL blade, ice model produces very high Mxf ultimate loads. For ultimate load analysis, the new linear ice mass formula (1) is multiplied by a factor of two for representing a worst-case scenario. The new proposed linear blade ‘IEC’ ice mass formula (1) produces similar results as the reduced GL × 0.5 blade ice model.
Ultimate loads for tower Mxf assuming blades 1 and 2 are iced, blade 3 is ice free.
Conclusion and discussion
Wind turbine fatigue loads during iced rotor normal power production operating mode was analysed by comparing measured mechanical loads to simulation cases with Flex5 aeroelastic simulation code. Simulated loads were compared to measured loads during icing events at TCE test site in Quebec, Canada, for a 2-MW Senvion MM92 turbine. Icing events were defined by measured wind speed difference between heated and unheated anemometers, and duration of icing event is the duration of instrumental icing. This duration can differ significantly from blade icing. Due to the selected icing event analysis method, normal operation with clean blade periods were included in measured icing events and target values were diluted increasing the method uncertainty. Higher fatigue loads would have been seen in single icing events and in shorter periods within events. All measured icing cases were merged together and individual icing events were not analysed. Aim of this work was to find average simulation parameters for iced turbine load analysis and the presented approach is suitable for that purpose.
In order to assess icing impacts on a more broad manner, conservatism was added intentionally to counteract uncertainties regarding short measurement campaign length, and milder than average winter season. In addition, measurements were performed only on one turbine type and results can be different with different turbine. Lifetime fatigue loads for iced turbine are turbine-type dependant and especially the turbine controller has huge effects to the loads during icing events. Also, tower height and flexibility can effect to side-to-side mass imbalance loads. Due to these uncertainties, conservative P90 values of measured fatigue loads were selected to target values for simulation.
Perfect match between simulations and measurements was not found. As the turbine performance measurement targets for simulations were defined statistically, there really never exists exactly this type of turbine behaviour in icing conditions. In reality, iced turbine behaviour changes in time (due to uneven ice accretion and ice shedding between blades) and in simulations, all parameters remain constant. Instead of minimizing mean absolute error, similar behaviour and operating point of turbine were searched from simulations. This means that most important signals were power and tower base side-to-side Mxf as power indicates the amount of aerodynamic penalty resulting to a shift in turbine operating point and Mxf as the most affected signal due to ice. The validated final simulation model parameters are conservative producing higher loads that what was measured, so that the proposed parameters could be globally utilized for other turbine types and models. Only tower top, tower base and blade root loads for single turbine were analysed in this study. Icing effect to other components or to other turbine types needs to be analysed separately.
Average effect of icing to tower or blade root lifetime is relatively small, except for tower base side-to-side loads. On one hand, uneven ice accretion and shedding on the turbine rotor will increase vibrations and decrease component lifetime. On the other hand, ice on the blades decreases airfoil lift and increases drag resulting static forces on the rotor, decreased power production and rotation speed which increase turbine lifetime. In extreme case, turbine cannot rotate at all with ice on the blades and fatigue loads are very low compared to even rotating ice-free turbine. Also, the wind turbine is more prone to side-to-side vibrations because of the rotor’s low aerodynamic damping compared to fore-aft direction (Hansen and Martin, 2008). Side-to-side vibration is mainly caused by mass imbalance of the rotor. In measurements, side-to-side loads were predominant and similar behaviour was selected as a target for simulations.
The final simulation model (applying an iced airfoil penalty factor to all blades and implementing an ice rotor mass imbalance to the rotor) is a simplification of reality, but it is intended to be an engineering approach to a very complex phenomenon of turbine rotor icing. The intention on this engineering iced turbine model is to be a global model applicable to any turbine model or type for making it possible to turbine designer to take into account cold and icing climate affects more realistically than before. The final simulation model parameters were used in proposing new icing design load cases for the new IEC 61400-1 ed4.
Footnotes
Acknowledgements
The authors would like to thank TechnoCentre Éolien for measurement data and Senvion for load measurement data and simulation model.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.
Funding
This work was funded by Norden IceWind and TEKES IcedBlades (part of IPSPoC) projects.
