Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind Turbine under Various Icing Conditions

Abstract

The phenomenon of icing on wind turbines gives rise to significant liability concerns in regions characterized by cold and humid climates, especially those with extreme climatic conditions. Accordingly, investigating the icing characteristics is essential for the safety operation of wind turbines. In this chapter, an icing model coupling water film flow with water film evaporation considering airfoil surface roughness is developed to investigate the effect of icing conditions on the icing characteristics of a blade airfoil for vertical-axis wind turbines by numerical simulation. The mechanism of heat and mass transfer under various icing conditions is explored. The results show that the simulated and experimental ice shapes on the airfoil agree well. The ice shape contour fluctuates along the airfoil surface at higher ambient temperature due to water film flow and heat flux variation. A large area of airfoil surface is covered by ice accretion at high wind speed due to an increase in driving force acting on water film and convective cooling between water film and air. The maximum ice thickness changes more significantly at wind speeds of 2–7 m/s than that at wind speeds of 7–12 m/s. This contributes to theoretical basis for exploring anti/de-icing method in wind turbines.

1. Introduction

Wind energy is widely used in electricity generation due to huge reserve and cleanness [1, 2]. Owing to the fact that wind turbines are installed in humid and cold regions, ice may accumulate on the wind turbine blades. It causes liability issues such as aerodynamic performance degradation, increased noise, decreased lifetime, ice shedding, and safety risks [3, 4]. Therefore, investigating the icing characteristics on wind turbines is crucial to explore anti/de-icing methods.

In the last few decades, numerical simulation methods can provide the insight of flow and heat transfer physics in detail, which is not easily possible using icing wind tunnel [5]. Numerical simulation is therefore reliable-efficient ways to investigate the icing characteristics of wind turbines under various icing conditions. For instance, Fu et al. [6] built the model to investigate the effect of angular speed on the icing characteristics at different locations of blades for a horizontal-axis wind turbine by numerical simulation. It was found that ice accumulation increased at the tip of blades with an increase in angular speed. Homola et al. [7] conducted numerical simulation to investigate the influences of ambient temperature and medium volume droplet diameter (MVD) on the icing characteristics. It was found that ice horns were observed at higher ambient temperature, and ice accretion increased with an increase in MVD. Etemaddar et al. [8] employed numerical approach to analyze the effect of ice accretion on the aerodynamic performance at various icing conditions. It was found that wind power loss increased by up to 35% when ice accumulated on the blades in below rated wind speed. Villalpando et al. [9] described the ice accretion shapes on different positions of a blade under icing conditions by numerical simulation. Yirtici et al. [10] employed numerical simulation to predict the influence of ice accretion on the power losses under various icing conditions. Jin and Virk [5] compared the ice accretion shapes on symmetric and asymmetric airfoils under icing conditions and reported that the ice accretion shapes were more streamlined on symmetrical airfoils as compared to asymmetric airfoils. Ibrahim et al. [11] evaluated the influence of blade thickness on the icing characteristics, and reported that the ice accretion thickness at the location of flow stagnation increased significantly with an increase in blade thickness due to higher droplet collection efficiency. Wang et al. [12] developed an improved multi-shot icing model to investigate the effect of ice accretion shapes on the aerodynamic performance on various positions of airfoils under yaw condition, and found that ice horn at 94% span position resulted in aerodynamic performance degradation due to flow separation. Manatbayev et al. [13] developed an icing model considering rotating effect to analyze the effect of ice accretion shapes on the aerodynamic performance for vertical axis wind turbines at various angles of attack and found that the glaze ice had a significant effect on aerodynamic performance degradation. Son et al. [14, 15] presented a more accurate method with moving reference frame and surface roughness effect to predict maximum ice accretion thickness and growth direction on different positions of airfoils. Baizhuma et al. [16] proposed an icing model coupling sliding mesh technique with multiple reference frame to analyze the effect of icing characteristics on the aerodynamic performance for rotating wind turbines under various icing conditions.

From the literature review, it is found that icing conditions including ambient temperature and wind speed have a great effect on the icing characteristics of wind turbine blades. Currently, there are some limitations in the above research as follows:

1. Most of researchers developed models to investigate the water film flow and icing process on clean airfoils under various icing conditions by numerical simulation. However, less attention has been paid to the effect of the water film flow and evaporation on the icing characteristics of rough airfoil surfaces.

2. The mechanism of the icing distribution on a rough airfoil of wind turbines is not clear.

Therefore, the main purposes of this study are as follows: (a) to propose the relationship between water evaporation and water film flow considering surface roughness effect; (b) to compare simulated and experimental ice shapes in order to validate the feasibility of numerical method; and (c) toinvestigate the heat and mass transfer mechanism of the ice accretion on a blade airfoil of wind turbines under different icing conditions.

2. Modeling

2.1 Mathematical model

2.1.1 Airflow model

The external flow and heat transfer characteristics of wind turbine blades is investigated according to conservation of mass, momentum, and energy, which is expressed as follows [17]:

Continuity equation

$$\frac{\partial {\rho }_1}{\partial t}+\nabla \cdot \left({\rho }_1\overrightarrow{V_1}\right)=0$$E1

Momentum equation

$$\frac{\partial \left({\rho }_1\overrightarrow{V_1}\right)}{\partial t}+\nabla \cdot \left({\rho }_1{\overrightarrow{V}}_1\overrightarrow{V_1}\right)=\nabla \cdot {\sigma }^{\mathit{ij}}+{\rho }_1\overrightarrow{g}$$E2

$${\sigma }^{\mathit{ij}}=-{\delta }^{\mathit{ij}}{p}_1+{\tau }^{\mathit{ij}}$$E3

The Spalart-Allmaras model is employed at low Reynolds number to determine reliable wall shear force and heat transfer values for near-wall regions by numerical approach, which is expressed as follows [11, 18, 19]:

$$\begin{array}{l}\frac{D\overline{\nu }}{\mathit{Dt}}=c_{\mathrm{b}1}\left(1-f_{\mathrm{t}2}\right)\overline{S}\overline{\nu }+\frac{1}{\sigma {Re}_{\infty }}\left\{\nabla \cdot \left[\left(\overline{\nu }+\nu \right)\nabla \overline{\nu }\right]+c_{\mathrm{b}2}{\left(\nabla \overline{\nu }\right)}^2\right\}\\ -\frac{1}{{Re}_{\infty }}\left(c_{\mathrm{w}}{f}_{\mathrm{w}}-\frac{c_{\mathrm{b}1}}{{\kappa }^2}{f}_{\mathrm{t}2}\right){\left(\frac{\overline{\nu }}{d}\right)}^2+{Re}_{\infty }{f}_{\mathrm{t}1}{\left(\mathit{\Delta V}\right)}^2\end{array}$$E4

$$\overline{S}=S+\frac{1}{{Re}_{\infty }}\frac{\overline{\nu }}{{\kappa }^2{d}^2}{f}_{\nu 2}$$E5

$$S=\sqrt{2{\mathrm{\Omega }}_{\mathit{ij}}{\mathrm{\Omega }}_{\mathit{ij}}}$$E6

$${\mathrm{\Omega }}_{\mathit{ij}}=\frac{1}{2}\left(\frac{\partial V_i}{\partial x_j}-\frac{\partial V_j}{\partial x_i}\right)$$E7

where $c_{\mathrm{b}1}$, $c_{\mathrm{b}2}$, $c_{\mathrm{w}}$, $\kappa$, and $\sigma$ are the closure coefficients, $f_{\mathrm{w}}$ is the closure function, and $f_{\mathrm{t}1}$ and $f_{\mathrm{t}2}$ are tripping functions from laminar to turbulent flow.

Energy equation

$$\frac{\partial \left({\rho }_1{E}_1\right)}{\partial t}+\nabla \cdot \left({\rho }_1{\overrightarrow{V}}_1{H}_1\right)=\nabla \cdot \left[k_1\left(\nabla T_1\right)+V_i{\tau }^{\mathit{ij}}\right]+{\rho }_1\overrightarrow{g}{\overrightarrow{V}}_1$$E8

2.1.2 Droplet flow model

The Eulerian two-fluid approach is employed to get droplet flow distribution on the blade airfoil for wind turbine according to continuity and momentum equations with droplet volume fraction, which is expressed as follows [20, 21]:

$$\frac{\partial \alpha }{\partial t}+\nabla \cdot \left(\alpha {\overrightarrow{V}}_2\right)=0$$E9

$$\frac{\partial \left(\alpha {\overrightarrow{V}}_2\right)}{\partial t}+\nabla \cdot \left(\alpha {\overrightarrow{V}}_2{\overrightarrow{V}}_2\right)=\frac{C_{\mathrm{D}}{Re}_2}{24K}\alpha \left({\overrightarrow{V}}_1-{\overrightarrow{V}}_2\right)+\alpha \left(1-\frac{{\rho }_1}{{\rho }_2}\right)\frac{1}{{\mathrm{Fr}}^2}$$E10

The drag coefficient is determined by Reynolds number, which is given by [22]:

$$C_{\mathrm{D}}=\left(\frac{24}{{Re}_2}\right)\left(1+0.15{Re}_2^{0.687}\right),{Re}_2\le 1300$$E11

$$C_{\mathrm{D}}=0.4,{Re}_2\ge 1300$$E12

Meanwhile, droplet collection efficiency is given by [23]:

$$\beta =\frac{\alpha \overrightarrow{V_2}\cdot \overrightarrow{n}}{V_{\infty }}$$E13

2.1.3 Water film and ice accretion model

The water film velocity is considered as the linear film velocity distribution normal to the wall where zero velocity is imposed, which is expressed as follows [18, 24]:

$${\overrightarrow{V}}_f=\frac{y}{{\mu }_f}{\overrightarrow{\tau }}_{1,\text{wall}}$$E14

$${\tau }_{1,\text{wall}}=C_f\frac{1}{2}{\rho }_1{V_{\infty }}^2$$E15

$$C_f={\left[3.476-0.707ln\left(k_s/s\right)\right]}^{-2.46}$$E16

The mass conservation equation consists of the mass transfer caused by the incoming water droplets, water evaporation, and ice accretion, which is written by:

$${\rho }_f\left[\frac{\partial h_f}{\partial t}+\nabla \cdot \left({\overrightarrow{V}}_f{h}_f\right)\right]=V_{\infty }\mathit{\text{βLWC}}-{\dot{m}}_{\text{evap}}-{\dot{m}}_{\mathrm{ice}}$$E17

The energy conservation equation consists of the heat transfer caused by the incoming water droplets, phase change, ice accretion, convective, and radiative heat fluxes, which is calculated by:

$$\begin{array}{l}{\rho }_f\left[\frac{\partial h_f{c}_f\tilde{T}}{\partial t}+\nabla \cdot \left({\overrightarrow{V}}_f{h}_f{c}_f\tilde{T}\right)\right]=\left(c_f{\tilde{T}}_{2,\infty }+\frac{1}{2}{\overrightarrow{V}}_2^2\right)V_{\infty }\mathit{\text{LWCβ}}\\ -0.5{\dot{m}}_{\text{evap}}\left(L_{\text{evap}}+L_{\text{subl}}\right)+{\dot{m}}_{\mathrm{ice}}\left(L_{\mathrm{fus}}-c_{\mathrm{ice}}\tilde{T}\right)+\mathit{\sigma \epsilon }\left(T_{\infty }^4-T_f^4\right)+{\dot{Q}}_h\end{array}$$E18

The mass of water evaporation is given by:

$${\dot{m}}_{\text{evap}}=h_{\mathrm{m}}\left({\rho }_{\mathrm{s}}-{\rho }_{\infty }\right)=h_{\mathrm{m}}\left(\frac{P_{\mathrm{s}}}{{\mathit{RT}}_{\mathrm{s}}}-\frac{P_{\infty }}{{\mathit{RT}}_{\infty }}\right)$$E19

Eq. (19) is substituted into Eq. (18), which is obtained as follows:

$$\begin{array}{l}{\rho }_f\left[\frac{\partial h_f{c}_f\tilde{T}}{\partial t}+\nabla \cdot \left({\overrightarrow{V}}_f{h}_f{c}_f\tilde{T}\right)\right]=\left(c_f{\tilde{T}}_{2,\infty }+\frac{1}{2}{\overrightarrow{V}}_2^2\right)V_{\infty }\mathit{\text{LWCβ}}\\ -0.5h_{\mathrm{m}}\left(\frac{P_{\mathrm{s}}}{{\mathit{RT}}_{\mathrm{s}}}-\frac{P_{\infty }}{{\mathit{RT}}_{\infty }}\right)\left(L_{\text{evap}}+L_{\text{subl}}\right)+{\dot{m}}_{\mathrm{ice}}\left(L_{\mathrm{fus}}-c_{\mathrm{ice}}\tilde{T}\right)+\mathit{\sigma \epsilon }\left(T_{\infty }^4-T_f^4\right)+{\dot{Q}}_h\end{array}$$E20

The mass of water evaporation caused by water film flow considering roughness effect is written by:

$$\begin{array}{l}{\rho }_f\left[\frac{\partial h_f{c}_f\tilde{T}}{\partial t}+\nabla \cdot \left({\overrightarrow{V}}_f{h}_f{c}_f\tilde{T}\right)\right]=\left(c_f{\tilde{T}}_{2,\infty }+\frac{1}{2}{\overrightarrow{V}}_2^2\right)V_{\infty }\mathit{\text{LWCβ}}\\ -0.5\mathrm{St}{V}_f{{\mathrm{L}}_{\mathrm{e}}}^{-\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(\frac{P_{\mathrm{s}}}{{\mathit{RT}}_{\mathrm{s}}}-\frac{P_{\infty }}{{\mathit{RT}}_{\infty }}\right)\left(L_{\text{evap}}+L_{\text{subl}}\right)+{\dot{m}}_{\mathrm{ice}}\left(L_{\mathrm{fus}}-c_{\mathrm{ice}}\tilde{T}\right)+\mathit{\sigma \epsilon }\left(T_{\infty }^4-T_f^4\right)+{\dot{Q}}_h\end{array}$$E21

Stanton number is different for laminar and turbulent regimes according to empirical relationship [25], which is expressed as follows:

$$\mathrm{St}=\left\{\begin{array}{l}0.5C_f{{Pr}_{\mathrm{t}}}^{-\frac{2}{3}}\\ 0.5C_f/\left({Pr}_{\mathrm{t}}+\sqrt{0.5C_f}{{\mathrm{St}}_{\mathrm{k}}}^{-1}\right)\end{array}\right.$$E22

A series of functional equations are added to close the above equations due to unknown variables such as the water film thickness, the interface equilibrium temperature, and the ice accretion mass, which are expressed as follows:

$$\left\{\begin{array}{l}{h}_f\ge 0\\ {\dot{m}}_{\mathrm{ice}}\ge 0\\ h_f\tilde{T}\ge 0\\ {\dot{m}}_{\mathrm{ice}}\tilde{T}\le 0\end{array}\right.$$E23

2.1.4 Roughness model

Roughness model is expressed as follows [26]:

$${\left[\frac{\left(k_s/c\right)}{{\left(k_s/c\right)}_{\text{base}}}\right]}_{V_{\infty }}=0.4286+0.0044139V_{\infty }$$E24

$${\left[\frac{\left(k_s/c\right)}{{\left(k_s/c\right)}_{\text{base}}}\right]}_{\mathit{LWC}}=0.5714+0.2457\left(\mathit{LWC}\right)+1.2571{\left(\mathit{LWC}\right)}^2$$E25

$${\left[\frac{\left(k_s/c\right)}{{\left(k_s/c\right)}_{\text{base}}}\right]}_{T_{\mathrm{s}}}=0.047T_s-11.27$$E26

$${\left(k_s/c\right)}_{\text{base}}=0.001177$$E27

The value of sand-grain roughness is calculated by:

$$k_s={\left[\frac{\left(k_s/c\right)}{{\left(k_s/c\right)}_{\text{base}}}\right]}_{V_{\infty }}{\left[\frac{\left(k_s/c\right)}{{\left(k_s/c\right)}_{\text{base}}}\right]}_{\mathit{LWC}}{\left[\frac{\left(k_s/c\right)}{{\left(k_s/c\right)}_{\text{base}}}\right]}_{T_s}{\left(k_s/c\right)}_{\text{base}}c$$E28

2.2 Geometry and computational mesh

The chord length and span length of the blade airfoil are 0.5334 and 0.1 m, respectively. The computational domain was meshed by structured elements, as shown in Figure 1. The near-wall refinement method with inflation factor (1.1) was employed to generate the nonuniform grids on the blade airfoil. The mesh cells of 41,130 were used to guarantee high computational precision and improve computational efficiency due to less than 0.5% relative error between 41,130 and 89,356 mesh cells in numerical approach.

Figure 1.

Geometry and computational mesh of NACA0012 blade airfoil.

In-cloud measurements indicate that the liquid water content (LWC) changes from 0 to 5 g/m³ [27, 28]. The actual droplet size distribution in-cloud is represented by a single variable called the median volume diameter (MVD), which varies from 10 μm (in freezing fog) to 5000 μm (in freezing rain) [28]. In addition, the ambient temperature and wind speed have a great effect on the icing characteristics of wind turbine blades. According to the book reported by Battisti [28], under operating conditions of wind turbines, the ambient temperature changes in the range of 233–323 K, and wind speed varies in the range of 2–25 m/s. In this study, the icing conditions were selected as follows: liquid water content (LWC) of 0.6 g/m³ with medium volume diameter (MVD) of 50 μm, ambient temperatures of 263–268 K, and wind speeds of 2–12 m/s. An angle of attack was 12^o, and icing time was 30 min. The ice conditions were illustrated in Table 1.

Parameters	Value
Wind speed (m/s)	2–12
Ambient temperature (K)	263–268
Liquid water content (g/m^3)	0.6
Medium volume droplet diameter (μm)	50

Table 1.

The ice conditions.

2.3 Boundary condition and solution method

The far-field and wall boundary conditions were employed in external flow field and solid surface of a blade airfoil for wind turbines in numerical simulation, respectively. The coupled boundary condition was used to deal with fluid-solid interface. The detailed parameters were displayed in Table 2.

Type	Values
Far-field	Pressure: 101325 Pa; Temperature: 263–268 K; Speed: 2–12 m/s; LWC: 0.6 g/m^3
Wall	No slip, coupled condition

Table 2.

Boundary conditions.

Numerical approach was employed to investigate the icing characteristics on a blade airfoil of wind turbines under various icing conditions by FNESAP-ICE. Spatial discretization was performed by finite element method with a streamline upwind term. The matrix equations linearized by the Newton approach was iteratively solved using the implicit Gear scheme and a generalized minimum residual procedure. The time-marching scheme is explicit, and the automatic time step was employed to obtain the optimal time step for each grid/film speed combination, resulting in considerable computational cost savings. The ALE (Arbitrary Lagrangian Eulerian) formulation was used for the mesh displacement due to ice accretion in time or body motion.

2.4 Model validation

Numerical method is validated in terms of ice shapes on a NACA 0012 airfoil in the NASA Icing Wind Tunnel [29]. The ice conditions were as follows: angle of attack of 4^o, wind speed of 12 m/s, liquid water content (LWC) of 0.55 g/m³ with medium volume droplet diameter (MVD) of 20 μm, ambient temperatures of 267 K and 247 K, and icing time of 420 s.

It is found that the simulated and experimental ice shapes on the airfoil agree well at ambient temperatures of 267 and 247 K, as shown in Figure 2. Ice horns are observed at the leading edge at ambient temperature of 267 K, as displayed in Figure 2a. The impinging water will immediately freeze into ice at the leading edge due to high droplet collection efficiency. On the other hand, the impinging water will partly form unfrozen water film due to the fact that the latent heat is not removed at relatively high ambient temperature. The runback water flows along the surface of blade airfoil where convective cooling rate is high, freezing to form ice horns. It is found that the rime ice appears at the leading edge when the ambient temperature is 247 K, as illustrated in Figure 2b. Compared with the glaze ice, ice horns are not observed at the leading edge under rime ice condition. This is attributed to the fact that all water droplets impinging at the leading edge will immediately freeze into ice instead of water film formation at relatively low ambient temperature, and ice accretion grows gradually in the direction of the impinging water droplets, forming smooth and streamlined ice shape.

Figure 2.

The simulated and experimental ice shapes on the airfoil (a) Ambient temperature of 267 K and (b) ambient temperature of 247 K.

3. Results and discussion

3.1 The effect of ambient temperature on icing characteristics

To investigate the influence of ambient temperature variation on ice accretion shape and rate, numerical analyses were performed at different ambient temperatures of 263–268 K. Figure 3 illustrates the ice shapes on the blade airfoil at different ambient temperatures for wind speed of 7 m/s. The collision efficiency is not affected by the ambient temperature variation, but the rate and shape of ice accretion change with the change in ambient temperature due to surface heat flux variation. At lower ambient temperature, the rate of ice accretion grows rapidly and the ice shape is smooth due to higher droplet freezing fraction. In contrast, at higher ambient temperature, most of droplets do not freeze into ice due to the fact that the latent heat released by solidification is not removed by air, forming unfrozen water film to flow along the airfoil surface caused by airflow effect. Therefore, the growth rate of ice accretion becomes slow, and the ice shape is not smooth at higher ambient temperature as compared to lower ambient temperature.

Figure 3.

The effect of ambient temperature on ice shapes (a) 268 K, (b) 265 K, and (c) 263 K.

Figure 4 shows the ice accretion thickness distribution on the blade airfoil at different ambient temperatures. With the ambient temperature decreasing from 268 to 263 K, it is observed that the maximum ice accretion thickness increases from 4.55 to 4.73 mm. This is mainly due to the fact that the droplet freezing fraction increases with the ambient temperature decreasing. The growth rate of ice accretion becomes faster at ambient temperature of 263 K than that at ambient temperature of 268 K. This is explained by the fact that most latent heat released by solidification is removed by air at lower ambient temperature. In addition, the ice accretion thickness fluctuates along the airfoil surface at higher ambient temperature due to water film flow and heat flux variation. In a word, the discrepancy in icing characteristics caused by various ambient temperatures can be explained by the insight of coupling water film flow with heat transfer in detail.

Figure 4.

The effect of ambient temperature on ice thickness.

3.2 The effect of wind speed on icing characteristics

Figure 5 shows the effect of wind speed on the ice shapes at ambient temperature of 263 K. It is found that the maximum ice accretion thickness is observed at the leading edge. This is due to the fact that the largest amount of impinging water can be collected at the leading edge, freezing to form the thickest ice layer at ambient temperature of 263 K. With the wind speed gradually increasing, the ice accretion at the leading edge grows along the direction, which is normal to the blade surface or tangential to trajectory of the blade airfoil. This can be attributed to the fact that the amount of the impinging water increases as the wind speed increases, forming thicker ice layer due to convective heat transfer enhancement caused by higher wind speed. In addition, a wider ice shape contour covered on the blade airfoil is observed at higher wind speed. The impinging water partly forms unfrozen water film due to the fact that the latent heat released by solidification cannot be immediately removed by air. Thus, the unfrozen water film driven by airflow moves toward the further downstream where the convective cooling rate enhances due to higher wind speed, resulting in freezing into ice layer.

Figure 5.

The effect of wind speed on ice shapes at 263 K.

Figure 6 shows the ice accretion thickness on the blade airfoil for different wind speeds at ambient temperature of 263 K. It is found that ice accretion distribution is wider at wind speed of 12 m/s than that at wind speed of 2 m/s. When the wind speed increases in the range of 2–12 m/s with steps of 5 m/s, the maximum ice accretion thickness varies from 0.783 to 8.16 mm. This is explained by the fact that more impinging water droplets accumulate at the leading edge with an increase in wind speed, resulting in forming thicker ice accretion caused by solidification. Meanwhile, unfrozen water film flows along the surface of blade airfoil due to shear force of air. Owing to the fact that the driving force acting on water film is proportional to wind speed, a wider water film distribution on the surface of blade airfoil is obvious at higher wind speed. Further, the convective heat transfer rate between air and water film enhances due to high wind speed, leading to an increase in ice accretion area. This explains the effect of various wind speeds on the icing characteristics.

Figure 6.

The effect of wind speed on ice thickness at 263 K.

Figure 7 illustrates the relationship between the maximum ice accretion thickness and ambient temperature, and wind speed. It is obvious that the maximum ice thickness varies more significantly at wind speeds of 2–7 m/s than that at wind speeds of 7–12 m/s. This discrepancy can be attributed to the heat and mass transfer mechanism of water droplets on the blade airfoil. At lower wind speed, an increase in wind speed increases the mass rate of impinging water droplets at the leading edge, which is more dominant than water film evaporation caused by wind speed. Meanwhile, owing to the fact that convective heat transfer rate on the surface of water film enhances in terms of Stanton number, the water film will freeze to form thicker ice layer. At higher wind speed, the evaporation rate of water film and driving force acting on water film will become more significant. Therefore, the growth rate of the maximum water film become slower at high wind speed as compared to low wind speed. This explains the difference in ice growth rate at various wind speeds. The maximum ice thickness reduces slightly with an increase in ambient temperature. This is explained by the fact that the latent heat caused by solidification is not removed immediately at higher ambient temperature, leading to a decrease in ice accretion thickness. Additionally, investigating the maximum ice accretion thickness will lay a foundation for exploring anti-icing energy demand on the blade airfoil of wind turbines.

Figure 7.

The maximum ice thickness for various temperatures and wind speeds.

4. Conclusions

In this study, an icing model coupling water film flow with water film evaporation considering airfoil surface roughness has been developed to investigate the effect of icing conditions (ambient temperature and wind speed) on the icing characteristics (ice shape and ice thickness) of a blade airfoil for wind turbines by numerical simulation. The mechanism of heat and mass transfer for wind turbine icing under various icing conditions has been explored. The following conclusions have been obtained.

1. The simulated and experimental ice shapes on the airfoil agree well under different icing conditions.

2. With the ambient temperature decreasing from 268 to 263 K, the maximum ice accretion thickness increases from 4.55 to 4.73 mm due to an increase in the freezing fraction.

3. The ice shape contour fluctuates along the airfoil surface at higher ambient temperature due to water film flow and heat flux variation.

4. With the wind speed increasing in the range of 2–12 m/s, the maximum ice accretion thickness varies from 0.783 to 8.16 mm due to an increase in the droplet collection efficiency.

5. A large area of airfoil surface is covered by ice accretion at high wind speed due to an increase in driving force acting on water film and convective cooling between water film and air.

6. The maximum ice thickness changes more significantly at wind speeds of 2–7 m/s than that at wind speeds of 7–12 m/s.

Acknowledgments

This study was supported by the National Youth Foundation of China (Grant No. 52106228).

Conflict of interest

The author declares no conflict of interest.

Appendices and nomenclature

c	chord length (m)
cice	specific heat capacity of ice (J/kg K)
cf	specific heat capacity of water film (J/kg K)
CD	drag coefficient
d	characteristic length (m)
E	internal energy (J/kg)
Fr	Froude number
g→	gravity vector (m/s2)
hf	water film thickness (m)
k	thermal conduction coefficient (W/m K)
ks	surface roughness (m)
K	inertia parameter
L	latent heat (J/kg)
Le	Lewis number
LWC	liquid water content (g/m3)
ṁ	mass rate per unit area (kg/m2 s)
n→	surface normal vector
Ω	pressure (Pa)
Q̇h	heat flux (W/m2)
Re	Reynolds number
S¯	vorticity (1/s)
St	Stanton number
t	time (s)
u	temperature (K)
T˜	equilibrium temperature (°C)
V→	velocity (m/s)
y	distance normal to the wall (m)

Greek symbols

α	droplet volume fraction
β	droplet collection efficiency
ε	emissivity
μ	dynamic viscosity (Pa s)
ν¯	kinematic viscosity (m2/s)
ϕ	density (kg/m3)
τ	shear stress (kg/m s2)

Subscripts

1	air
2	droplet
evap	evaporation
fus	fusion
f	water film
ice	ice
subl	sublimation
∞	free stream