Numerical study of icing impact on the performance of pitch-regulated large wind turbine

This paper presents a study of the impact of icing on the performance of a pitch-regulated large wind turbine. Numerical simulations of six blade sections of the NREL 5 MW wind turbine at various free stream velocities are performed. Blade Element Momentum (BEM) method along Computational Fluid Dynamics (CFD) bases multiphase numerical simulations are used for this study. Analysis shows that the simulated parameters are in good agreement with the real conditions for each blade element during operation, except for the three-dimensional effects. The analysis of accreted ice shapes and air/droplet flow fields around the blade profile sections was carried out, and the calculation of aerodynamic performance, and energy production degradation was also performed. The tip of the blade is most affected by icing, it is characterized by the greatest changes in the aerodynamic performance. Maximum reduction in the wind turbine performance is estimated to be around 24%.


Introduction
The Covid-19 pandemic had an impact on global energy demand in 2020. Energy consumption as well as CO 2 emissions were reduced during this period. But now economic growth and global energy demand are recovering rapidly after the impact of the pandemic. According to the Global Energy Review (2021) global economic output grew by 5.9% in 2021, impacting a 6% increase in CO 2 emissions. In absolute values, CO 2 emissions in 2021 exceeded the pre-pandemic 2019 emissions level by 180 megatons. At the same time, the post-Covid era is characterized by an increase in renewable energy production, which reached an all-time high in 2021 (The International Energy Agency, 2021). Wind power generation increased by 270 TWh, which shows the largest growth among all renewable energy sources. There is now 743 GW of wind power capacity worldwide, avoiding more than 1.1 billion tons of CO 2 globally, equivalent to South America's annual carbon emissions (Global Wind Energy Council, 2021).
Many places that are promising for the construction of wind power plants are located in the high north regions (Abbey et al., 2006;Alm and Nygaard, 1995). The wind turbines in cold climate regions can produce more power than those located in the other regions as air density is higher at low temperatures (Carreno-Madinabeitia et al., 2021). But high wind speeds combined with low air temperatures, and supercooled water droplets may cause ice formation on the structures of wind turbines particularly the rotor blades. Icing can lead to degradation of blade aerodynamic characteristics, losses in power production, partial, or complete mechanical failure of the wind turbine, measurement errors, risk to human health, and life (Parent and Ilinca, 2011). Therefore, the operation of wind turbines in cold climate conditions requires additional technical solutions and a good understanding of the physics of ice accretion.
Ice mainly accumulates along the wind turbine blades and can change their aerodynamic profile shape, increases its surface roughness, which affects its aerodynamic characteristics. Icing results in reduced lift force and increased drag force of the blade (Barber et al., 2011;Bragg et al., 2005;Homola et al., 2010;Jin and Virk, 2020). Depending on the accreted ice shape, the loss in lift coefficient can be between 10 and 65% (Ibrahim et al., 2018). The total torque decreases when the lift force of the wind turbine blades decreases, which leads to a reduction in power output (Botta et al., 1998;Kraj and Bibeau, 2010;Nixon, 2020;Zanon et al., 2018). The results of field studies show that degradation of wind turbine blade aerodynamics caused by icing leads to a significant decrease in rotor speed and in some cases even shutdown of wind turbines due to the lack of sufficient torque to rotate wind turbines. Power generation losses caused by icing are strongly dependent on the duration of icing. For a 30-hour icing period, losses can reach up to 80% (Gao et al., 2021).
Ice accretion on wings and aerodynamic profiles of aircraft is widely studied (, 1952;Gulick, 1938;Jacobs, 1934;Jones and Williams, 1936;Pettit, 1959). At the same time, much fewer publications are devoted to icing of wind turbines. And this problem can still be considered insufficiently studied.
The type of wind turbine control system at high wind speeds affects the control strategy of the operation. Some studies are devoted to stall regulated wind turbines , whose blades are designed in a way that the rotor speed remains almost constant and energy production does not increase above the required limit due to the stall effect when wind speed increases. On the other hand, there are pitch regulated wind turbines (Yuan and Tang, 2017) with an active control system using where the blades are rotated around their axis to reduce the torque. Homola et al. (2012) investigated the icing of the pitch-regulated wind turbine NREL 5 MW. Shapes of the accreted ice were obtained. Aerodynamic and power performance characteristics were calculated. But the research was carried out only for an inflow velocity of 10 m/s. Han et al. (2018) studied the ice accretion along the leading edge of the blade tip. The ice accretion model was validated by comparing the simulation results with previously published wind tunnel test results. As a result of the study, the performance of the NREL 5 MW wind turbine was calculated and compared before and after icing. However, the aerodynamic performance calculation was carried out only for NACA 64-618 airfoil and for an inflow velocity of 11.1 m/s. Etemaddar et al. (2014) studied the impact of various parameters (angle of attack, relative velocity, chord length, aerodynamic profile, temperature, liquid water content, median volumetric diameter, and relative humidity) on the shape of ice formation. But CFD calculation of the aerodynamic characteristics was performed only for two airfoils.
This paper is aimed to investigate the impact of icing on the performance of pitch-regulated large wind turbine using a combined approach -BEM method and multiphase CFD modeling at different velocity conditions. The paper will investigate icing along six sections of the NREL 5 MW wind turbine blade from the root to the tip ( Figure 2). Simulations are performed for four different inflow velocities to assess the impact of icing on aerodynamic performance and power production of the wind turbine blade.

Numerical approach
In this study, several numerical simulation tools were applied to calculate the characteristics of a horizontal axis pitch-regulated wind turbine under icing conditions. The applied workflow includes the following elements: blade element momentum (BEM) code; CFD simulations with flow, water droplet and ice solvers.
BEM method. The BEM methodology is widely used for wind turbine design and calculation of the main wind power characteristics (Yirtici et al., 2019). BEM is described in detail by Glauert (1935) and Hansen (2008). This approach allows to calculate the forces acting on the wind turbine blade, and accordingly the energy output of the wind turbine, using the geometric parameters of the planar sections of the blade, and the characteristics of the incoming flow. In this paper, the BEM method is used to calculate the power production of a wind turbine based on the aerodynamic forces calculated by the CFD method.
The BEM method divides the blade into a finite number of sections, which are approximated by planar models. The method estimates the forces acting in the cross section on each element through lift and drag coefficients as a function of geometric parameters and the inflow angle. The total force values for the blade are determined by numerical integration along the blade span (Hansen, 2008).
The BEM theory assumes some simplifications: there is no aerodynamic interaction between the elements, the blade forces are determined only by the lift and drag characteristics, the free flow is perpendicular to the rotation plane. Thus, the BEM method does not determine the three-dimensional influence of blade elements on each other. But it allows to use two-dimensional modeling, which saves computational resources considerably.
The BEM theory relates the axial and tangential induction coefficients to the aerodynamic forces acting on the turbine blade. The axial induction coefficient represents the amount of air flow passing inside the rotor disk. The tangential induction coefficient determines how many rotations are created in the flow after the rotor (Ibrahim et al., 2018). The axial a and tangential a' induction factors are defined as: where f is the inflow angle, C n is the normal force coefficient, C t is the tangential force coefficient, s is the rotor solidity.
The thrust T and torque M on the element of thickness dr are defined as: where B is the number of blades, r is the radial position of the element, p n is the normal component of the vector sum of the lift and the drag, p t is the tangential component of the vector sum of the lift and the drag.
There are many implementations of BEM codes in software. In this paper, the open access software QBlade was used to perform BEM calculations (Marten and Wendler, 2013). This software is developed for design and aerodynamic calculations of wind turbine blades. New tip-and root loss models after Shen et al. (2005) were used. The new Shen tip loss model overcomes the inconsistencies of the Prandtl correction model, which overestimates tip loads compared to experimental data. Also, 3D correction after Snel et al. (1992) and foil interpolation were used in the simulation to get more accurate results.
CFD modeling. Multiphase CFD simulation codes are widely used to predict ice formation on surfaces, such as aircraft wings, or wind turbine blades, using a multiphysics analysis involving heat transfer and the CFD approach. In this paper, ANSYS-FENCEP-ICE software is used for simulation of icing of a planar blade element of a horizontal-axis pitch-regulated wind turbine. CFD modeling allows to predict ice accretion shapes, ice mass, calculate the flow field and forces acting on the blade.
The modeling process of ice formation and the study of the effect of icing on wind turbine performance consists of the following steps: flow solution, water droplet/ice crystal collision modeling, ice accretion modeling, mesh displacement for the iced airfoil, re-solving the flow for the iced airfoil (see Figure 1).
FENSAP air flow solution is a 3D finite element Navier-Stokes analysis package. Analysis of the icing problem begins with an airflow solution over a clean geometry and ends with a series of airflow solutions over an iced geometry to evaluate the performance degradation caused by ice accretion. FENSAP solves steady and unsteady compressible three-dimensional Navier-Stokes equations: where r is the density, V is the velocity vector, s ij is the stress tensor, g is the gravity vector. DROP3D is 3D finite element Eulerian water droplet impingement solver. DROP3D provides water concentration, droplet velocity vectors, water catch efficiency distributions, impingement patterns, shadow zone characteristics and impingement limits over the entire domain without the laborious iterative procedure of seeding droplets at injection points.
ICE3D is 3D finite volume ice accretion and water runback solver. This module based on fine-grain partial differential equations for the complex thermodynamics of ice formation. It yields 3D ice shape, water film thickness and surface temperature on any number of complex 3D surfaces. This paper uses a single-shot solution methodology. The flow field and droplet solutions are calculated once, and the icing solution is computed for the total icing duration as a single interval. The results of the flow field calculation are used as input to the droplet solver. The droplet solution is then passed to the ice accretion simulation module to predict the ice mass and shape for a single icing time interval. The single-shot solution is chosen for modeling the icing process because it requires less computational resources compared to the multi-shot solution.

Model features
The geometry of the NREL 5 MW wind turbine was used in the paper as a reference model for the study.  Table 1.
Six sections were chosen for this study (see Figure 2, Table 2) to reduce the use of computing resources. Most of the sections are on the outer part of the blade, since according to (Homola et al., 2012;Shu et al., 2018) ice mostly accrete closer to the blade tip section compared to the root section.   The CFD calculation was performed using the k-w SST turbulence model (Menter, 1993;Wilcox, 1988). This model is widely used, in particular for flow calculations near wind turbine blades, and provides reliable flow calculations (Menter, 1994;Moshfeghi et al., 2012;Rocha et al., 2014).

Mesh
For CFD modeling, a two-dimensional C-shaped orthogonal grid model was created. The boundaries in front, above, and below the airfoil are spaced at 15 chord lengths, and the boundary behind is spaced at 25 chord lengths (see Figure 3). The height of the first cell is calculated based on the requirement that the value of the near-wall function y+is less than one for the simulated cases (Gantasala et al., 2019;Zanon et al., 2018).
An insufficient number of cells can increase the error of the obtained results, and an excessive number of cells can unnecessarily increase the computation time. Therefore, a study of grid independence on a clean NACA 64-618 airfoil was carried out. Six meshes with different number of cells were constructed. To change the number of cells, the number of nodes along the main directions -above, below, before, and after the profile, and the number of nodes on the airfoil were changed.
The drag and lift coefficients were calculated for a Reynolds number of 6 million and an angle of attack of 0°. The difference in the results of calculating the coefficients when using meshes with different numbers of cells and the time taken to perform simulations for a given mesh were compared. The results are shown in the Figure 4. A mesh with 153,000 cells has been selected for further simulations considering a reasonable balance of accuracy of the results and optimal use of computational resources. The delta Cl and delta Cd values for this mesh are 0.03% and 0.28%, respectively. A mesh with 248,000 cells reduces this error, but at the same time requires almost twice as much time per simulation, therefore it was not chosen. Zanon et al. (2018) chose a mesh with 140 k cells for

Results
The incoming flow characteristics were calculated for six different blade sections located at radial distances r/ R = 0.19, 0.45, 0.58, 0.71, 0.84, 0.98 (sections A to F, respectively) along the blade of the NREL 5 MW wind turbine. Each section corresponds to the airfoil: A -DU40; B -DU25; C -DU21; D, E, F -NACA64-618. Four inflow velocities (cases 1-4) were chosen to calculate the blade aerodynamic characteristics and wind turbine performance -10, 15, 20, and 25 m/s. The inflow free stream velocities were chosen based on the NREL 5 MW wind turbine power curve. The maximum speed of 25 m/s was determined by the cut-out speed of the wind turbine according to the power curve.
The angle of attack, axial and tangential induction factors and relative flow velocity were calculated. Since the NREL 5 MW is pitch regulated wind turbine for each inflow velocity the pitch angle and tip speed ratio were calculated. The pitch angles are different from zero after reaching nominal power on the power curve, i.e., after a nominal wind speed of 11.4 m/s. The pitch angle was calculated for clean blade conditions. Table 3 contains the results of calculating the swept flow characteristics for Cases 1-4. For each of the six selected sections of the NREL 5 MW wind turbine blade, the simulation and calculation of the aerodynamic characteristics before and after the icing event were performed. Table 4 contains the icing simulation conditions. The option ''Sand-grain roughness -file'' was selected for the flow simulation after icing. In this case, the roughness values on the wall surface are imported from the input file, which is generated during the simulation of ice formation on the surface of the airfoil. Figure 5 shows the simulated shapes of the ice accretion on the blade profiles at the inflow velocity of 10 m/s. The results show that ice accretion on the leading edge of the airfoil occurs near the stagnation point. The stagnation point is closer to the pressure side because the angle of attack was positive. The ice formation increases as it moves from the root section of the blade (Section A) to the tip of the blade (Section F). For Section A, which is closest to the blade root (r/R = 0.19), ice formation is extremely insignificant. The maximum ice thickness according to the simulation results was 9.7Á10 24 m. The maximum ice formation is observed for the Section F closest to the blade tip (r/R = 0.98). The stagnation point is on the pressure side, and the largest area of ice formation is on this side of the blade. But partially the ice formation passes to the suction side, where it has the maximum thickness of 2.5Á10 22 m. These results prove that the tip section of the wind turbine blade is exposed to more ice formation compared to the root section. Figure 6 presents the icing simulation results for Section F for the Cases 1-4. The Section F was the most exposed to icing is the reason for selecting this section for the comparison of the results of the four cases.
The results of the cases provided a comparison of the ice formation shapes during wind turbine blade operation with the pitch control strategy for different inflow velocities. For Case 1, the ice accretion is mostly along the pressure side (inflow velocity 10 m/s) because the angle of attack is positive. The angle of attack changes to negative after reaching a nominal velocity of 11.4 m/s and simulating pitch control system conditions for a clean blade. Therefore, for cases 2-4 (inflow velocities 15, 20, and 25 m/s), the ice accretion gradually moves to the opposite side of the blade. At the same time, the ice accretion thickness does not change significantly for different cases.
Two-dimensional simulation of the flow around the airfoils was repeated after the analysis of the ice shapes to assess the influence of the ice accretion on the aerodynamic performance of the wind turbine blade sections. Figure 7 shows the relative velocity distribution fields around the airfoil for Section F (r/R = 0.98), where the greatest ice formation along the leading edge was observed. Cases 1 and 4 were analyzed before and after blade icing. For Case 1 (inflow velocity 210 m/s, angle of attack 26.5°) no significant flow changes are observed: the flow circulation around the blade is generally symmetrical, there are vortices at the leading edge caused by ice formation, as well as a slight flow separation at the trailing edge of the airfoil on the suction side. For Case 4 (inflow velocity 225 m/s, angle of attack 25.9°) a non-deep stall caused by strong flow circulation near the leading edge due to ice formations and flow separation was observed.
A comparison of the aerodynamic characteristics before and after the icing event was performed to assess the degradation of the blade aerodynamic performance. Figure 8 illustrates the change in the lift coefficient due to leading edge icing for Sections A-F under Case 1 conditions. It is obvious that a higher level of icing on the leading edge of the sections closer to the blade tip leads to a more significant change in the aerodynamic performance. The lift coefficient for Section A does not change, and for Section F it decreases by 13.39%. It is important, because Then, the changes of lift and drag coefficients (see Table 5) were analyzed in more detail for Cases 1-4 for Section F, as exposed to the most icing. Negative values indicate a decrease in the coefficient compared to the  clean profile, while positive values indicate an increase. Obviously, there is a significant degradation of the aerodynamic characteristics for the inflow velocities .20 m/s. Figure 9 demonstrates the pressure coefficient distribution for Section F before and after icing. The iced airfoil is associated with a change in the pressure coefficient distribution. The icing causes disturbances at the leading edge, which causes a reduction of the pressure difference. This results in a degradation of the blade aerodynamic characteristics.  An analysis of the impact of icing on the power production of the wind turbine was the final stage of the study. For a given wind speed, the CFD simulation results for each section are the lift and drag coefficients. These coefficients are used to calculate the lift and drag forces acting on selected blade sections. The contribution to the torque of each section is then calculated. The total blade torque value is found by integrating the values for the sections along the blade span using BEM method. Finally, the power extracted by the wind turbine is determined by multiplying the torque of one blade by the number of blades and by the rotation speed. Figure 10 shows the results of the power calculation for wind turbine under icing conditions. For an inflow free stream velocity of 10 m/s, the reduction in power production is not significant. But with increasing inflow velocity the degradation of performance gradually increases. The maximum performance reduction is 24%.

Conclusion
In this paper, the authors investigated the impact of icing on the performance of a pitch-regulated wind turbine using a combined approach -Blade Element Momentum code and CFD modeling. The NREL 5 MW reference wind turbine was used as a test case for icing simulations. The study included simulations for six airfoils (Sections A-F) from root to tip section along the blade span. Four cases with different inflow velocities were considered to evaluate the impact of icing on wind turbine performance and create a power curve after icing. An air temperature of 210°C and a duration of 1 hour were used as the icing simulation conditions. The two-dimensional simulation allowed to simulate the flow conditions that are specific to wind turbine rotation, except for the three-dimensional influence of the sections on each other. Relative velocity, inflow angle, axial and tangential induction factors were calculated for each section.  Figure 9. Pressure coefficient distribution for Section F before and after icing.
The results showed that the sections of the blade close to the blade tip were the most affected by icing. For section A (r/R = 0.19), it was noticed that ice formation was very less. The maximum ice thickness according to the simulation results was 9.7Á10 24 m. The largest ice formation was observed for section F (r/R = 0.98). The stagnation point was on the pressure side, the largest area of ice formation was on this side of the blade, but the ice had the greatest thickness on the suction side 22.5Á10 -2 m. For cases 2-4, the ice shapes move to the upper side of the airfoil due to the change in angle of attack, but the ice thickness is almost the same for different cases. The icing of the leading edge of the blade causes changes in the aerodynamics of the flow -swirls appear near the leading edge, with increasing inflow velocity the flow separation increases, and stall appears.
Calculation of aerodynamic characteristics demonstrated that the greatest losses are observed in the sections close to the tip of the blade due to large ice formation. Also, with the increase of inflow velocity the aerodynamic characteristics degradation raised. Accordingly, this influenced energy performance -energy losses increased with increasing inflow velocity and reached a maximum of 24% at an inflow velocity of 25 m/s. This study was conducted using two-dimensional simulation and the Blade Element Momentum theory, which does not take into account the aerodynamic interaction between the blade elements. Therefore, the next stage of this study may be to conduct simulation using three-dimensional rotor model and compare obtained results with the results of this work.

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.