Abstract
Most icing studies use 2D or quasi-3D models, ignoring the effects of blade span and tip icing. In this paper, a new model for predicting 3D vertical axis wind turbine blades icing by supplementing the blade span is established, the icing characteristics of blades under multiple working conditions are simulated step by step from the time scale, and the unsteady icing laws are compared and analyzed. Ice properties and surface temperature distribution determine the corresponding deicing power. The analysis results show that the unfrozen water film flows in the spanwise direction of the blade, resulting in a smaller icing limit on the pressure surface than the two-dimensional model, which indicates that the chordwise anti-icing structure can be set in a small range. At the same time, the annular ridge ice formed at the end of the blade slows down the reduction rate of the lift-drag ratio of the blade, effectively reducing the influence of the leading edge ice shape on the aerodynamic performance of the blade. Finally, the critical ice melting power density required to reach equilibrium temperature at the leading edge of the blade is reduced by about 11% compared to the trailing edge due to the increased thermal resistance of the gas-to-blade matrix interface by the ice layer.
1. Introduction
The low temperature in cold regions results in high air density, which is suitable for the development of wind energy [1]. More than 30% of the wind turbines installed in cold regions in the world face the problem of icing. The icing changes the aerodynamic shape of the blades and reduces the lift-to-drag ratio of the blades. The annual loss of power generation due to icing is as high as 50% [2]. Vertical Axis Wind Turbine (VAWT) can capture wind energy in all directions without complex yaw mechanism. It is widely used in cities with high population density, especially the problems of ice throwing, noise, additional load, health, and safety caused by icing are particularly prominent [3].
Icing is a natural phenomenon, and the main research methods include theoretical analysis, numerical simulation, and ice-wind tunnel experiments [4]. Due to the limited size of the experimental section, Li et al. used an isometric model to conduct an ice wind tunnel experiment. Temperature changes in the natural environment play a major role in icing. When the ambient temperature is between 0°C and -10°C, the water droplets hit the leading edge of the blade, delaying freezing to form a high-roughness ice shape [5]. When the ambient temperature is lower than -10°C, the water droplets freeze instantaneously to form streamlined ice similar to the blade profile. At this time, the icing only increases the thickness of the airfoil and has little effect on the aerodynamic performance [6].
The numerical method can efficiently evaluate the performance of the VAWT under different icing conditions, reduce the experimental cost, and improve the analysis efficiency [7]. Among them, the 2D model is based on a single airfoil to explore the icing situation of the blade. That is, take the shape of a certain section of the blade and build a model on a two-dimensional plane. Based on the 2D model, Shu et al. found that static icing occurs mainly on the leading edge of the blade [8]. Baizhuma et al. added centrifugal force terms and Coriolis force terms to the flow control equations to account for the effect of rotation on icing, which occurs on the entire blade [9]. The quasi-3D model superimposes the icing conditions of a single airfoil to form a complete blade icing. That is, the blades are intercepted at certain intervals to obtain multiple blade sections, each section is calculated and analyzed separately, and the results of each section are connected in series to characterize the icing condition of the entire blade, which is mostly used for large horizontal axis wind turbine blades. The research results of Yirtici et al. show that the power loss of a 30 kW horizontal axis wind turbine can reach 20% for one hour when it freezes [10].
It is worth noting that all previous studies were based on 2D or quasi-3D models. At the same time, the icing simulation process adopts the one-way steady-state icing method. However, 2D or quasi-3D models ignore the spanwise flow characteristics of the unfrozen water film on the blade surface and cannot fully reflect the icing state of the blade in the 3D flow field [11, 12]. Secondly, during the icing process, the flow field and the trajectory of water droplets determine the icing shape, and the change of the aerodynamic shape of the blade after icing will affect the former. The flow field, water droplet field, and temperature field are continuously interacting with each other, and the icing of wind turbine blades is essentially the result of a bidirectional fluid-solid-thermal multifield coupling synergy.
In order to improve the ability of VAWT to cope with cold climates, the icing mechanism under the synergistic effect of multiple fields was explored, and the thermal power distribution on the surface of the anti-icing blade was optimized. In this paper, a 3D VAWT blade icing model is established to supplement the blade spanwise effect on icing. FENSAP-ICE is used to predict the icing characteristics of the unfrozen water film flowing in the blade spanwise direction in the 3D flow field. In the temperature range of -10~0°C with the highest icing frequency, the icing characteristics of blades under multiple working conditions are simulated. Comparative analysis of the aerodynamic performance changes of the icing blades in the three-dimensional flow field. Finally, the anti-icing power of the corresponding area is determined according to the three-dimensional blade icing characteristics.
The detailed of the current method is given in Section 2. Section 3 verifies the accuracy of the flow field, water droplet field, and temperature field in the model. Section 4 applies this method to discuss the blade icing characteristics, aerodynamic performance changes, and the deicing power density in corresponding regions under different working conditions. The conclusions will be followed in Section 5.
2. Numerical Method
In this section, the numerical calculation method of static VAWT blade icing is described, and the icing model establishment process is shown in Figure 1.
The flow state of the external flow field in FENSAP-ICE is calculated by solving the flow governing equation, calculating the velocity, pressure, and temperature at each grid node, and obtaining the surface heat flow and shear force. On the basis of the external flow field, the water droplet impact characteristics and icing shape are calculated. Among them, the multishot method is used for the calculation of unsteady icing of blades. The icing process is divided into three intervals according to the time scale. The calculation result is used as the initial condition of the next interval to realize the transfer of information in the continuous process [13].
The angle of attack of the fixed tip speed ratio VAWT varies with the azimuth angle, as shown in Equation (1). In this study, the flow fields at 0° and 8° attack angles were calculated, and the transient icing characteristics of the VAWT when rotated to different azimuth angles were compared and analyzed.
The aerodynamic performance of the blade is more sensitive to geometrical changes in the leading edge of the airfoil; so, other components such as struts are not considered. The fluid domain boundary does not interfere with the flow state around the blade, and the minimum dimension is 20 times the chord length [14]. The computational domain includes the blade model established by the NACA0018 airfoil, and the span length is 0.5 m. The fluid domain is divided into structural meshes by C-type segmentation, and the blade end meshes are refined by O-type segmentation to ensure the orthogonality of the end meshes, as shown in Figures 2 and 3. When the incoming wind speed is 10 m/s and the angle of attack is 4°, the grid independence is verified by changing the grid height and node number of the first layer.
As shown in Figure 4, when the number of grids exceeds , the changes of lift coefficient and drag coefficient tend to be stable with the increase of grid number [15]. The viscosity calculation is carried out in the boundary layer near the wall of the blade. The wall function is set to 5, the mesh growth rate is 1.2, the thickness of the first layer mesh is 0.15 mm, and the total number of meshes after division is 718,775.
2.1. Airflow Solution
Wind turbine blades work in a low-speed flow field, and it is approximately considered that air is an incompressible fluid [16]. The viscosity of air will cause the separation of the gas boundary layer of the blade, and the flow state will change from laminar flow to turbulent flow. The one-equation Spalart-Allmaras turbulence model with low free turbulence is used to calculate the solid-wall turbulence problem, a transport equation is used to solve the dynamic eddy viscosity, and the effect of turbulence is expressed in the momentum equation [17]. The momentum equation of three-dimensional flow can be represented by three nonlinear equations, and its vector form is shown in formula (2):
In the calculation of viscous laminar boundary layer, the Prandtl number and surface recovery factor can be calculated based on the dynamic viscosity, thermal conductivity, and specific heat of air at constant pressure, which are used for the calculation of blade surface heat transfer [18]. The dynamic viscosity and thermal conductivity are calculated according to formula (3) and formula (4):
where , is the air static temperature, , and .
The icing surface roughness increases with the heat transfer rate. As shown in equations (5)-(7), the NASA icing empirical correlation is used to calculate the ice surface roughness [19]: where is the chord length, is the free flow velocity, is the liquid water content, and =0.001177.
Finally, the NASA icing surface roughness is calculated as shown in equations (8) and (9):
2.2. Droplet Field Solution
The Euler method is used to establish the motion equation of the water droplet in the three-dimensional flow field. The water droplet is regarded as a continuous phase. The continuity equation and momentum equation of the water droplet phase are solved by the finite volume method, as shown in formulas (10) and (11). Among them, it is assumed that the initial velocity of the water droplet is equal to the incoming flow velocity, and the volume of the water droplet does not affect the flow field quality of the surrounding flow. The water droplet volume fraction and water droplet velocity distribution of each node in the spatial grid are introduced, and finally, the water droplet impact limit and water droplet collection efficiency on the blade surface are obtained.
As shown in Figure 5, during the movement of supercooled large water droplets, under the combined action of pressure and surface tension, the droplets change from spherical to oblate disk shape. The drag coefficient of the droplets continues to increase during the movement, eventually leading to the breakup of the droplets. The small water droplets formed by crushing are easily affected by the flow field and deviate from the initial motion trajectory, increasing the impact limit of the water droplets. In addition, after the water droplet hits the blade surface, it bounces into the flow field due to the reaction force of the blade surface and hits the blade surface again under the action of the airflow. This process will reduce the water droplet collection efficiency on the blade surface.
The Mundo model determines the splashing and bouncing of droplets through the parameters defined by the Weber number and the Oh number. It is suitable for large droplets with a median volumetric diameter of more than 40 μm, and splitting occurs when the critical Weber number defined by Pilch and Erdman is exceeded. The critical Weber number is as follows Equations (12) and (13) show: where is the viscosity of the water drop, is the size of the water drop, and is the surface tension of the water drop.
The water droplet collection efficiency is calculated as Equation (14): where is droplet volume fraction, and is droplet velocity.
2.3. Icing Thermodynamic Model
The icing of wind turbine blades is an unsteady process. The icing shape on the blade surface is determined by the collection efficiency of water droplets. There is a continuous interaction between the air flow field, the motion trajectory of water droplets, and the icing shape. The flow field and the trajectory of water droplets determine the icing state. During the icing process, blade geometry changes, which in turn lead to changes in surface heat flux, water droplet collection efficiency, and shear stress. The control volume elements on the blade surface were taken for analysis, and the mass transfer was shown in Figure 6.
In the energy conservation equation, the kinetic energy of impinging water droplets, latent heat of fusion, evaporation, and sublimation, as well as heat of radiation and convective cooling, is added as a source term. The mass and energy conservation equations are shown in formulas (15) and (16):
The flow characteristics of the backflow water are represented by the velocity of the surface water film. The water film is subjected to shear force, and the velocity at in the microelement is shown in formula (17). Since the water film has a certain thickness, the overall motion speed of the water film is represented by the average speed, as shown in formula (18): where and are the mass of water droplets flowing into the control volume element. , , and are the mass of water droplets flowing out of the control volume element. , , and are the energy flowing into the control volume element. , , , , and are the energy lost by the control volume element.
The thickness of the icing on the blade surface is calculated by equation (19):
3. Validations
In this section, aerodynamic solution, droplet impingement, and ice accretion modules are validated with existing results.
The VAWT blade uses NACA0012 as the reference airfoil. The span length is 0.5 m, and other airfoil data are quoted from reference [19]. The surface pressure distribution and icing shape of the blade are calculated by the method in this paper and compared with the experimental data of Manatbayev et al. [19]. The distribution of the pressure coefficient on the blade surface is shown in Figures 7 and 8.
The figure compares two cases where the angle of attack is 0° and 10°. When the angle of attack is 0°, the upper and lower surfaces of the blade have the same pressure distribution, which is in good agreement with the reference experimental data. When the angle of attack is 10°, the result obtained on the pressure surface of the blade is slightly lower than the reference data, but the overall flow field around the blade can be accurately calculated.
The accuracy of the drop field and temperature field is verified by comparing the final icing shape. The specific data are shown in Table 1. The verification results are shown in Figure 9.
The icing limit of the blade pressure surface calculated by this method is similar to the reference data, and the difference is only at the leading edge of the blade, which is caused by the flow of the water film on the blade surface. Therefore, the leading edge ice thickness is smaller, but the suction surface ice thickness is larger. The verification results show that this method can calculate the icing shape well.
4. Results and Discussion
4.1. Three-Dimensional Static Blade Icing Characteristics
In this section, the results of three-dimensional wind turbine blades at 0° and 8° attack angles are calculated to compare and analyze the transient icing characteristics of the VAWT rotating to different azimuth angles.
The blade adopts NACA0018 airfoil with a chord length of 0.5 m. The ambient temperature will cause the delayed freezing or instantaneous freezing of the supercooled water droplets hitting the blade surface. In the temperature range of 0 ~ 10°C, the icing characteristics of the blade at -4°C, -7°C, and -10°C are calculated parameters. The settings are shown in Table 2.
In the calculation process of case 1 and case 2, the mesh distortion becomes negative volume, and the unstructured mesh with better geometrical adaptability is used instead. Under the condition of case 1, the ice shapes of the two types of grid models were compared for 10 minutes, and the accuracy verification results are shown in Figure 10. The icing limit calculated by the unstructured grid is slightly higher than that of the structured grid, but the overall ice shape is the same. Therefore, the mesh type does not have an impact on the calculation results of the cases in this study.
The shape of the icing on the blade surface is shown in Figures 11 and 12. The water film flowing in the spanwise direction of the blade causes the icing limit on the pressure surface to be smaller than that calculated by the two-dimensional model.
(a) ,
(b) ,
(c) ,
(a) Case1
(b) Case2
(c) Case3
At -4°C, the ice accumulated on the suction surface to form an obvious upper ice angle, and the sudden roughness would cause the separation of the gas boundary layer, resulting in flow loss. At an angle of attack of 8°, the water film flows on the pressure surface of the blade to form streamlined ice similar to the blade profile. At -7°C, the temperature difference between the blade surface and the ambient temperature increases, resulting in faster convective heat transfer and icing rate, lowering the height of the upper ice angle, and increasing the icing thickness. At the same time, a lower ice angle appears on the pressure surface. When the temperature continued to decrease to -10°C, the increase in temperature difference further increased the freezing rate of water droplets hitting the blade surface, and no obvious ice angle appeared.
At -4°C, the freezing ratio of supercooled water droplets decreases, decreases, increases in the energy conservation equation, and higher ice ridges cannot be formed at the tip of the blade, but the icing limit is increased. At -10°C, the freezing ratio is high, increases, decreases, the ice ridge height is higher, and the ice thickness increases significantly. The reduction of the water film flow leads to a decrease in of the next control volume microelement, which in turn limits the further increase of the ice ridge height, and thus reduces the icing limit. At -7°C, the is replenished by the backflow of unfrozen supercooled water droplets, resulting in the highest ice ridge height.
As shown in Figure 13, the tip icing limit decreases with decreasing temperature, and the ridge ice height first increases and then decreases with decreasing temperature. At -7°C, the ice ridge height is the highest.
4.2. Three-Dimensional Icing Blade Aerodynamic Analysis
In this section, the changes in lift coefficient, drag coefficient, and lift-drag ratio of 2D and 3D icing blades are compared to evaluate the effect of icing on the aerodynamic performance of the blade, and the differences between the 2D and 3D models are illustrated.
Figure 14 shows the variation of blade lift coefficient with time at different temperatures. The lift coefficient of the unfrozen 3D model is 52.5% lower than that of the 2D model, which is due to the periodic vortex shedding at the blade tip [20]. The lift coefficient decreases with the increase of icing time, and icing at -4°C causes the most lift coefficient loss. After freezing for 30 minutes, the lift coefficient of the two-dimensional model decreased by 46.7%, and the lift coefficient of the three-dimensional model decreased by 29%. This may be because the ice ridges formed after the icing of the blade ends perturb the flow field and act as a spoiler to slow down the reduction rate of the three-dimensional blade lift coefficient.
As shown in Figures 15 and 16, with the increase of icing time, the overall blade drag coefficient shows an increasing trend, and the overall lift-drag ratio shows a decreasing trend. Freezing at -4°C produces the highest drag coefficient. After freezing for 30 minutes, the drag coefficient of the two-dimensional model increased by 10.1%, and the drag coefficient of the three-dimensional model increased by 13.4%.
When icing for about 10 minutes, a small amount of icing forms Kruger flaps on the leading edge of the blade, and the windward attack angle of the blade decreases, which reduces the flow loss of the gas boundary layer of the blade, resulting in a decrease in the blade drag coefficient and an increase in the lift-to-drag ratio. With the increase of icing time, icing has little effect on the drag coefficient, but has a greater effect on the lift coefficient, the lift-drag ratio decreases and gradually becomes stable, and the aerodynamic performance of the blade decreases.
The calculation results of the three-dimensional model and the two-dimensional model have the same trend, and both reach the minimum lift-to-drag ratio at the freezing temperature of -4°C. The maximum lift-to-drag ratio is reached at the freezing temperature of -10°C. The difference is that in the three-dimensional flow, the separation vortex generated at the end causes the blade lift to be 49.9% lower than the two-dimensional model. The flow state at the end is shown in Figure 17.
(a) Case 1
(b) Case 2
(c) Case 3
The difference between the maximum and minimum lift-drag ratios of the two-dimensional model is 65.4%, while the lift-drag ratio of the three-dimensional model is only 26.7% different. That is to say, the three-dimensional blade lift-to-drag ratio is not sensitive to the change of icing temperature. In the temperature range of -10~0°C, the icing at the end of the blade reduces the influence of the leading edge ice shape on the aerodynamic performance of the blade. To sum up, the 3D model takes into account the effect of end icing on the overall performance. Therefore, the 3D model is more accurate than the 2D model.
Table 3 shows the blade icing parameters under different working conditions. Case 8 uses uniform water droplet size for calculation, and case 9 uses unidirectional icing process for calculation.
As shown in Figure 18(a), under the condition of , the leading edge of the blade forms a streamlined ice similar to the airfoil profile, and the increase of LWC and surface roughness will intensify the convective heat transfer on the blade surface. Compared with Figure 18(b), it can be seen that the increase of LWC will significantly change the shape of blade icing and increase the icing thickness.
(a)
(b)
In Figure 18(b), case 8 has a very high icing limit, which is because the liquid water in the environment adopts a single droplet size distribution, lacks large-sized water droplets, and does not consider the impact characteristics of large water droplets. Therefore, the droplet size and size distribution are crucial for calculating the chordwise icing characteristics of the blade. Case 9 calculated with a single icing process did not observe the formation of ice angle on the pressure surface and suction surface of the leading edge of the blade, which could not reflect the result of the backflow water flow and the coupling between the icing ice shape and the flow field.
Figure 19 shows the ice formation at the blade end. Under the influence of unsteady flow, the icing extends along the spanwise direction of the blade, and the icing ice is ridge-like. Under the influence of the airflow, the small water droplets formed by the breaking of large water droplets in the incoming air deviate from the initial motion trajectory, which reduces the collection efficiency of water droplets on the blade surface, resulting in a smaller icing range of case 7 than case 8, but the height of the ice ridges is similar. Therefore, the droplet size and size distribution did not affect the icing height at the tip of the blade.
The height behind the ice ridge in case 7 using multishot calculation decreases slowly, and this is because the unfrozen water droplets are formed by flowing at the end of the blade. Case 9 using a single icing calculation cannot accurately simulate the flow characteristics of backflow water and the coupling of ice shape and flow field.
4.3. Anti-Icing Power Analysis
In this section, the anti-icing heat load is applied to the blade section according to the abovementioned icing characteristics, so as to achieve an efficient energy ratio and a reasonable temperature distribution on the blade surface.
As shown in Figure 20(a), the blade is divided into front, middle, and rear parts, and the leading edge of the blade is subdivided into walls 1~4. The ice mass ratio of each region satisfies the formula (20):
(a)
(b)
where is the ice mass of the corresponding area, and is the summary ice mass of the blade.
Figure 20(b) is the percentage of ice mass in the blade zone. The icing mainly occurs in the 4% area of the leading edge of the blade, and the ice mass accounts for up to 90%. There is a small amount of icing on the trailing edge of the blade. Therefore, deicing should be carried out mainly on the leading edge of the blade.
The wind turbine blades are deiced by heating, and the critical ice-melting power density is used to describe the thermal balance relationship on the blade surface [21]. The flow field in the working environment of the wind turbine blade is a low-speed flow field, and the frictional heat generated by the air viscosity cannot be ignored. Therefore, this paper improves the original model and adds the frictional heat source term for analysis. The external heating structure is adopted, which does not involve the heat conduction of the blade base material.
is the minimum power density at which the temperature of the contact surface between the blade and the ice layer is , and the temperature of the interface between the ice layer and the airflow is less than 0°C, as shown in formulas (21) and (22):
where is the thermal power density released by cooling to after the phase transition of supercooled water droplets, and and are constants, 0.95 and , respectively. is the convective heat transfer coefficient, is the specific heat of air, is the ice growth rate per unit area, is the specific heat capacity of water, is the specific heat capacity of ice, and is the heat of fusion.
The critical ice melting power density of each region of the blade is shown in Figure 21. Compared with unfrozen blades, icing is equivalent to increasing the thermal resistance between the blades and the low-temperature airflow, reducing radiation and convective heat losses. The ice thickness of the leading edge of the blade is large, the thermal resistance of the ice layer is higher, and a small critical ice melting power density is required to melt the ice.
(a) No supercooled water droplets
(b) There are supercooled water droplets
As shown in Figure 21(a), when there are no supercooled water droplets in the air, at this time , the critical ice melting power density required for the leading edge of the blade to reach the equilibrium temperature is reduced by at most 11% compared to the trailing edge. As shown in Figure 21(b), when there are supercooled water droplets in the air, the energy brought by the freezing and impact of supercooled water droplets needs to be considered. The critical ice melting power density required to reach equilibrium temperature at the leading edge of the blade is reduced by up to 75% compared to the trailing edge.
5. Conclusion
This study provides a new method for numerical calculation of icing on static VAWT blades, explores the effects of icing, and analyzes the corresponding deicing power in different regions of the blade. The following conclusions can be drawn: (1)Static VAWT icing occurs at the leading edge of the blade, and the water film flows in the spanwise direction, which reduces the icing limit of the three-dimensional blade pressure surface. At -4°C and -7°C, the leading edge forms an angular ice shape, and the blade ends form a semicircular ice ridge. The height of the ice ridge first increases and then decreases with the decrease of temperature, and the end icing limit decreases with the decrease of LWC content and temperature.(2)VAWT icing significantly reduces the aerodynamic performance of the blade, and icing at -4°C causes the most power loss. After 30 minutes of icing, the lift coefficient decreases by 29%, the drag coefficient increases by 13.4%, and the lift-to-drag ratio decreases from 6.37 when it is not icing to 3.99. Compared with the two-dimensional model, the icy ridge at the end of the blade forms a spoiler, which slows down the reduction rate of the lift-drag ratio of the blade. When freezing for about 10 minutes, a small amount of icing will form a flap on the leading edge of the blade, which will slightly improve the aerodynamic performance of the blade. The actual three-dimensional model blade icing aerodynamic performance degradation is not serious.(3)The icing quality of the leading edge of the three-dimensional blade accounts for the highest proportion of 90%, which is the main deicing energy consumption area. The ice layer increases the thermal resistance between the blades and the airflow and slows down the convective heat transfer rate. The deicing power density required for the blade surface to reach the equilibrium temperature is smaller than that of the middle and trailing edges.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This work is financially supported by the Guangxi Natural Science Foundation (2018GXNSFAA050077), Project Director of Guangxi Manufacturing System and Advanced Manufacturing Technology Key Laboratory (20-065-40-006Z), and Guilin University of Electronic Science and Technology Graduate Education Innovation Project Funding Project (2021YCXS016), which is gratefully acknowledged by the authors.