Numerical Investigation on the influence of meteorological conditions on wind turbine icing

Keeping the blade structure and shape in the design condition is the premise and guarantee for efficient operation of wing turbine, but wet and cold climate will lead to icing blade, then the blade shapes were changed and the output power of the wind turbine is reduced. A numerical frame for wind turbine icing was developed in present work, which consists of a numerical method for impingement of water droplets based on the multiple frames of reference and Euler method and a 3D heat and mass transferring method for wind turbine icing, then the effects of several meteorological parameters, including the temperature, Liquid Water Content (LWC), Mean Volume Diameter (MVD) and the inflow velocity, on ice shape of the blade were investigated in the numerical frame. The present method can provide supports for prediction of icing process and aerodynamic characteristic after icing of wind turbine blades.


Introduction
The blade is the key component of the wind turbine, keeping the blade structure and shape in the design condition is the premise and guarantee for efficient operation of wing turbine [1].However, wet and cold climate will lead to icing blade, then the blade shapes were changed, and the output power of the wind turbine is reduced.A four-year study conducted by University of Quebec showed that the average power loss of wind turbine resulting from blade icing was 26% [2].Therefore, blade icing is a key issue for safe and efficient operation of wind turbines in high and cold regions [3].
At present, the icing wind test and numerical simulations have been used to study the influence of icing on wind turbines.Due to the limitation of the size of the test section of wind tunnel, the numerical simulation presents great potential for blade icing due to its low cost, fast calculation speed and wild application range.FLUENT, LEWICE and XFOIL have been employed in the investigation of blade icing [4].And the ice distributions on the blade [6], the icing morphology [6], and its harm to wind turbines [7][8] have been proposed in previous paper.
In previous studies, the blade rotation motion was simplified to non-rotation, and replacing 3D blade was replaced by several 2D sections in calculations, resulting in inaccurate predictions.Moreover, the influence of centrifugal force on icing was not considered, which is also one of the reasons for the reduction of calculation accuracy.In present work, a numerical frame for wind turbine icing was developed, which consists of a numerical method for impingement of water droplets based on the multiple frames of reference and Euler method and a 3D heat and mass transferring method for wind turbine icing, then the effects of several meteorological parameters, including the temperature, Liquid Water Content (LWC), Mean Volume Diameter (MVD) and the air speed, on ice shape of the blade were investigated in the numerical frame.The present method can provide supports for prediction of icing process and aerodynamic characteristic after icing of wind turbine blades

Numerical methodology
The water droplets in the air were transported to the blade surface where they were impacted and collected.Then the inertia and distribution of water droplets also feedback the flow field.In order to reproduce this phenomenon in numerical simulations, the air flow field and the water droplet field need to be coupled in the calculation.However, it is necessary to decouple the two processes in order to simplify the calculation, that is, the airflow and water droplets were solved separately based on the assumption that there is no reaction to the convection field of water droplets.

Governing equation for air
In the rotation region, the reference coordinate system was located on the blade, so that the unsteady problem of blade rotation was transformed into a steady problem when the governing equation was established, and the variables remain exactly the same in the interface between the rotation and stationary regions.
The relative velocity is employed in the mass conservation equation for the rotational region, that is, the interface mass flux is calculated by the relative velocity.Thus, the governing equations can be written in general form as r r ( ) Where ρ is density,   is diffusion coefficient, r q  is source term,  ⃗ is air velocity in rotating coordinate system.
The  in Equation (1) represents variables.When  =1 and r q  =0, Equation ( 1) is the continuity equation; when   ⃗, Equation ( 1) is the momentum equation and the source term can be written as Where q  is the source term of non-rotating momentum equation,  ⃗ is the angular velocity,  ⃗ is radius vector of fluid particle in rotating coordinate system.
At the interface between the rotation and stationary regions, the vector velocity in the rotating and stationary coordinate systems can be linked as follows Where u  is the air velocity in non-rotating region.

Governing equation for water droplets
Water droplets exist as particles in the air, which is suitable for modelling by Lagrange method.However, the Euler method was used by most researchers due to the excessive computational cost.The water droplets are treated as continuous fluids in frame of Euler, After the volume fraction was introduced, the continuous and momentum equations can be written as follows Where  ⃗ is the droplet velocity, d  is the droplet density,  ⃗ is acceleration of gravity, and the inertial factor K can be defined as 18 Re 24 Where a  is aerodynamic viscosity coefficient, d is water droplet diameter, D C is drag coefficient of water droplet, and Re is relative Reynolds number.The mass flow of water droplets can be obtained by the equation as follows Where  ⃗ represents normal direction of the wall.

Icing model
Messinger model was employed in present work for calculating icing shapes.For clear ice, the mass of frozen water can be obtained by solving the mass conservation and energy conservation equations, then the icing shape can be obtained.For structured grid, the mass conservation and energy conservation equations can be written as Where the mass of impinged water im m  can be obtained from Equation (7), va m  is the mass of evaporated water, so m  is icing mass in the controlled body per unit time, in m  and out m  represents respectively the water mass of inflow and outflow controlled volume.so E  is the energy leaving the control body result from freezing, including the heat from freezing water, latent heat and sensible heat during freezing process.im H  is mass of leaving the controlled body result from evaporating, including the heat from evaporating water and latent heat during evaporating process.f Q  is the heat generated by water droplets impinged with surface of the controlled volume, c Q  is the heat from air friction, is convective heat transfer between airflow and the blade wall.
in H  and out H  Represent the heat flowing into and out the controlled volume of overflow water.The subscript nb represents the four directions of controlled volume: e, w, n, s.

Temperature
There are four cases with different temperature, the parameters had been listed in Table 1, have been calculated for investigating the effect of temperature on ice shape.The predicted ice shape was presented in Figure 1.When the temperature was -20℃, the temperature is too low to freeze the water droplets after the collision, and there is no overflow.Consequently, the numerical result presents a small icing range in chordwise, a large icing thickness and frost ice on the surface.When the temperature rise to -15℃, the numerical result presents a blunt leading edge, ice horn on up and lower surface, and frost ice.When the temperature rise to -10℃, the blunt leading edge still can be presented, but the chordwise freezing range increases, the calculation show that the overflow and the smaller icing thinkness.Furthermore, the ice morphology changes from frost ice to mix and clear ice.When the temperature rise to -5 ℃ , the numerical result presents increasing overflow, larger freezing range and smaller ice thinkness, then the ice morphology is clear ice.

Liquid Water Content
In order to investigating the influence of the liquid water content (LWC), four cases were calculated and compared with each other, the parameter was provided in Table 2.The computational ice shape was presented in Figure 2. The numerical results show that the amount of ice increases with the increase of LWC, which can be found in increasing ice thinkness and developing significantly downstream of the upper and lower limit ice position.However, the ice thickness around the stagnation point remains unchanged when the LWC is greater than a certain value.Furthermore, under the condition of mass and flow convservation, the mass of frozen water changes little and the mass of overflow water increases , then the water flows to the trailing edge by rotating linear velocity and centrifugal force, and a small part escapes from the tip.

Mean Volume Diameter
In order to investigating the influence of mean volume diameter (MVD), five cases were calculated and compared with each other, the parameter was provided in Table 3.The computational ice shape was presented in Figure 3.The numerical results present little change with MVD around the stagnation point but the thickness of the upper and lower limit ice position.The collection rate of water droplets leads to differences in ice thickness and range.The inertia of the particles increases with the increasing MVD, then the aerodynamic following becomes worse, and the impact range and value on the blade will be expanded.

Inflow velocity
In order to investigating the influence of the inflow velocity, three cases were calculated and compared with each other, the parameter was provided in Table 4.The computational ice shape was presented in Figure 4.The numerical results show that the icing range and amount increase with the increasing rotational speed.With increasing rotational speed and flow velocity, the impacting range and momentum become larger.However, the little difference can be found in icing thickness around the leading edge when the rotational speed and incoming flow increase to a certain value.Because the increasing impact velocity of water droplets leads to increasing kinetic energy, the proportion of freezing proportion is relatively reduced although the amount of water collected inceases.

Conclusion
A numerical frame was proposed in present work to predict the 3D icing shape for wind turbine.Then the effects of several meteorological parameters, including the temperature, LWC, MVD and the air speed, on ice shape of the blade were investigated in the numerical frame.The numerical result shows that: (1) The effect of temperature on icing is mainly reflected in the icing morphology, the higher temperature results in the larger icing range and clear ice.
(2) The Higher LWC results in the larger icing range and amount, but the ice thickness around the stagnation point remains unchanged when the LWC is greater than a certain value.
(3) The Larger MVD results in the larger icing range.
(4) The higher RPM results in the larger icing range and amount, but the little difference can be found in icing thickness around the leading edge when the rotational speed and incoming flow increase to a certain value

Figure 1 .
Figure 1.Influence of static temperature on icing shape.

Figure 2 .
Figure 2. Influence of LWC on icing shape.

Figure 3 .
Figure 3. Influence of MVD on icing shape.

Figure 4 .
Figure 4. Influence of inflow velocity on icing shape.

Table 1 .
The Parameter of Calculations for Temperature.

Table 2 .
The Parameter of Calculations.for LWC

Table 3 .
The Parameter of Calculations for MVD.

Table 4 .
The Parameter of Calculations for Inflow Velocity.