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

*Zhi Xu*

## **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.

**Keywords:** wind turbine blade, icing condition, icing characteristics, numerical simulation, heat and mass transfer mechanism

## **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:


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.

*Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

## **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 \tag{1}$$

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^\# + \rho\_1 \overrightarrow{g} \tag{2}$$

$$
\sigma^{\ddagger} = -\delta^{\ddagger} p\_1 + \tau^{\ddagger} \tag{3}
$$

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{split} \frac{D\overline{\nu}}{Dt} &= c\_{\text{b1}} (1 - f\_{\text{t2}}) \overline{\text{Sv}} + \frac{1}{\sigma \operatorname{Re}\_{\text{os}}} \left\{ \nabla \cdot \left[ (\overline{\nu} + \nu) \nabla \overline{\nu} \right] + c\_{\text{b2}} (\nabla \overline{\nu})^2 \right\} \\ &- \frac{1}{\operatorname{Re}\_{\text{os}}} \left( c\_{\text{w}} f\_{\text{w}} - \frac{c\_{\text{b1}}}{\kappa^2} f\_{\text{t2}} \right) \left( \frac{\overline{\nu}}{d} \right)^2 + \operatorname{Re}\_{\text{o}} f\_{\text{t1}} (\Delta V)^2 \end{split} \tag{4}$$

$$\overline{\mathcal{S}} = \mathcal{S} + \frac{1}{\text{Re}\_{\text{os}}} \frac{\overline{\nu}}{\kappa^2 d^2} f\_{\nu 2} \tag{5}$$

$$\mathcal{S} = \sqrt{2\Omega\_{\vec{\eta}}\Omega\_{\vec{\eta}}} \tag{6}$$

$$\boldsymbol{\Omega}\_{\dot{\boldsymbol{w}}} = \frac{1}{2} \left( \frac{\partial V\_i}{\partial \mathbf{x}\_{\dot{\boldsymbol{\beta}}}} - \frac{\partial V\_j}{\partial \mathbf{x}\_i} \right) \tag{7}$$

where *c*b1, *c*b2, *c*w, *κ*, and *σ* are the closure coefficients, *f* <sup>w</sup> is the closure function, and *f* t1 and *f* t2 are tripping functions from laminar to turbulent flow.

Energy equation

$$\frac{\partial(\rho\_1 E\_1)}{\partial t} + \nabla \cdot \left(\rho\_1 \overrightarrow{V}\_1 H\_1\right) = \nabla \cdot \left[k\_1(\nabla T\_1) + V\_i \overrightarrow{\pi}^{\ddagger}\right] + \rho\_1 \overrightarrow{\text{g}}^{\ddagger} \overrightarrow{V}\_1 \tag{8}$$

#### *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]:

*Wind Turbine Icing – Recent Advances in Icing Characteristics and Protection Technology*

$$\frac{\partial a}{\partial t} + \nabla \cdot \left(a\vec{V}\_2\right) = 0\tag{9}$$

$$\frac{\partial \left(a\vec{V}\_2\right)}{\partial t} + \nabla \cdot \left(a\vec{V}\_2\vec{V}\_2\right) = \frac{\mathbf{C}\_\mathrm{D}\mathbf{Re}\_2}{24K}a\left(\vec{V}\_1 - \vec{V}\_2\right) + a\left(\mathbf{1} - \frac{\rho\_1}{\rho\_2}\right)\frac{\mathbf{1}}{\mathbf{Fr}^2} \tag{10}$$

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

$$\mathbf{C\_D} = \left(\frac{24}{\mathrm{Re}\_2}\right) \left(\mathbf{1} + \mathbf{0.15} \,\mathrm{Re}\_2^{0.687}\right), \,\mathrm{Re}\_2 \le \mathbf{1300} \tag{11}$$

$$\mathbf{C\_D = 0.4, \,\mathrm{Re}\,\_2\mathrm{}^\circ\mathrm{1300}\tag{12}$$

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

$$
\beta = \frac{\overrightarrow{a \, V\_2} \cdot \overrightarrow{n}}{V\_{\infty}} \tag{13}
$$

#### *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]:

$$
\vec{V}\_f = \frac{\mathcal{Y}}{\mu\_f} \vec{\tau}\_{1, \text{wall}} \tag{14}
$$

$$
\sigma\_{1, \text{wall}} = C\_f \frac{1}{2} \rho\_1 V\_\infty^{\ 2} \tag{15}
$$

$$\mathbf{C}\_{f} = \left[\mathbf{3.476} - \mathbf{0.707}\ln(k\_{\mathrm{i}}/s)\right]^{-2.46} \tag{16}$$

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\_{\approx 0} \beta L \text{WC} - \dot{m}\_{\text{evap}} - \dot{m}\_{\text{ice}} \tag{17}
$$

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{aligned} \rho\_f \left[ \frac{\partial h\_f c\_f \tilde{T}}{\partial t} + \nabla \cdot \left( \vec{V}\_f h\_f c\_f \tilde{T} \right) \right] &= \left( c\_f \tilde{T}\_{2,\text{oc}} + \frac{1}{2} \vec{V}\_2^2 \right) V\_{\text{oc}} L W \text{C} \beta \\ - \mathbf{0}.5 \dot{m}\_{\text{evap}} \left( L\_{\text{evap}} + L\_{\text{subl}} \right) + \dot{m}\_{\text{ice}} \left( L\_{\text{fus}} - c\_{\text{ice}} \bar{T} \right) &+ \sigma e \left( T\_{\text{oc}}^4 - T\_f^4 \right) + \dot{Q}\_h \end{aligned} \tag{18}$$

The mass of water evaporation is given by:

$$\dot{m}\_{\text{evap}} = h\_{\text{m}} (\rho\_{\text{s}} - \rho\_{\text{os}}) = h\_{\text{m}} \left( \frac{P\_{\text{s}}}{RT\_{\text{s}}} - \frac{P\_{\text{os}}}{RT\_{\text{os}}} \right) \tag{19}$$

*Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

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

$$\rho\_f \left[ \frac{\partial h\_f c\_f \bar{T}}{\partial t} + \nabla \cdot \left( \vec{V}\_f h\_f c\_f \bar{T} \right) \right] = \left( c\_f \ddot{T}\_{2,\text{so}} + \frac{1}{2} \vec{V}\_2^2 \right) V\_{\infty} L \text{WC} \beta \text{} \tag{20.54}$$

$$-0.5 h\_{\text{m}} \left( \frac{P\_s}{RT\_s} - \frac{P\_{\infty}}{RT\_{\infty}} \right) \left( L\_{\text{evap}} + L\_{\text{subl}} \right) + \dot{m}\_{\text{ice}} \left( L\_{\text{fus}} - c\_{\text{ice}} \bar{T} \right) + \sigma \varepsilon \left( T\_{\infty}^4 - T\_f^4 \right) + \dot{Q}\_h \tag{20.55}$$

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

$$\begin{split} \rho\_f \left[ \frac{\partial h\_f c\_f \bar{T}}{\partial t} + \nabla \cdot \left( \vec{V}\_f h\_f c\_f \bar{T} \right) \right] &= \left( c\_f \bar{T}\_{2,\text{so}} + \frac{1}{2} \vec{V}\_2^2 \right) V\_{\text{os}} L \text{WCC} \beta \\\ -\mathbf{0.5St} V\_f \text{L}\_{\text{e}} &\quad \left( \frac{P\_s}{RT\_s} - \frac{P\_{\text{os}}}{RT\_{\text{os}}} \right) \left( L\_{\text{vmp}} + L\_{\text{nhl}} \right) + \dot{m}\_{\text{ixc}} \left( L\_{\text{f}\text{m}} - c\_{\text{lxi}} T \right) + \sigma \epsilon \left( T\_{\text{o}}^4 - T\_f^4 \right) + Q\_k \end{split} \tag{21}$$

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

$$\mathbf{St} = \begin{cases} \mathbf{0.5C\_f Pr\_t}^{-\frac{2}{\beta}} \\ \mathbf{0.5C\_f} / \left( \mathbf{Pr\_t} + \sqrt{\mathbf{0.5C\_f}} \mathbf{St\_k}^{-1} \right) \end{cases} \tag{22}$$

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:

$$\begin{cases} h\_f \ge 0 \\ \dot{m}\_{\text{ice}} \ge 0 \\ h\_f \bar{T} \ge 0 \\ \dot{m}\_{\text{ice}} \bar{T} \le 0 \end{cases} \tag{23}$$

### *2.1.4 Roughness model*

Roughness model is expressed as follows [26]:

$$
\left[\frac{(k\_t/c)}{(k\_t/c)\_{\text{base}}}\right]\_{V\_{\text{os}}} = 0.4286 + 0.00444139 V\_{\text{os}} \tag{24}
$$

$$
\left[\frac{(k\_\text{s}/c)}{(k\_\text{s}/c)\_\text{base}}\right]\_{L\text{WC}} = 0.5714 + 0.2457(L\text{WC}) + 1.2571(L\text{WC})^2\tag{25}
$$

$$
\left[\frac{(k\_\text{s}/c)}{(k\_\text{s}/c)\_\text{base}}\right]\_{T\_\text{s}} = 0.047T\_\text{s} - 11.27\tag{26}
$$

$$(k\_{\circ}/c)\_{\text{base}} = 0.001177 \tag{27}$$

The value of sand-grain roughness is calculated by:

$$k\_s = \left[\frac{(k\_s/c)}{(k\_s/c)\_{\text{base}}}\right]\_{V\_{\text{os}}} \left[\frac{(k\_s/c)}{(k\_s/c)\_{\text{base}}}\right]\_{L\text{WC}} \left[\frac{(k\_s/c)}{(k\_s/c)\_{\text{base}}}\right]\_{T\_s} (k\_s/c)\_{\text{base}}c \tag{28}$$

#### **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.

In-cloud measurements indicate that the liquid water content (LWC) changes from 0 to 5 g/m3 [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<sup>3</sup> 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<sup>o</sup> , and icing time was 30 min. The ice conditions were illustrated in **Table 1**.

#### **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**.

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 timemarching 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 4o , wind speed of 12 m/s, liquid water content (LWC) of 0.55 g/m<sup>3</sup> with medium volume droplet diameter (MVD) of 20 μm, ambient temperatures of 267 K and 247 K, and icing time of 420 s.

*Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

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

**Figure 1.** *Geometry and computational mesh of NACA0012 blade airfoil.*


#### **Table 1.**

*The ice conditions.*


#### **Table 2.**

*Boundary conditions.*

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.

### **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 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 *Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

**Figure 2.**

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

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.

#### **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

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

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 *Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

**Figure 4.** *The effect of ambient temperature on ice thickness.*

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

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 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

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

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 7** illustrates the relationship between the maximum ice accretion thickness and ambient temperature, and wind speed. It is obvious that the maximum

**Figure 7.** *The maximum ice thickness for various temperatures and wind speeds.*

*Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

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.

## **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.


## **Acknowledgements**

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**


## **Greek symbols**


*Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*


## **Subscripts**


## **Author details**

Zhi Xu1,2

1 College of Engineering, Northeast Agricultural University, Harbin, China

2 Heilongjiang Provincial Key Laboratory of Technology and Equipment for the Utilization of Agricultural Renewable Resources in Cold Region, Northeast Agricultural University, Harbin, China

\*Address all correspondence to: zxu@neau.edu.cn

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **References**

[1] Kusiak A, Zhang Z, Verma A. Prediction, operations, and condition monitoring in wind energy. Energy. 2013;**60**:1-12. DOI: 10.1016/j. energy.2013.07.051

[2] Willis DJ, Niezrecki C, Kuchma D, Hines E, Arwade SR, Barthelmie RJ, et al. Wind energy research: State-of-the-art and future research directions. Renewable Energy. 2018;**125**:133-154. DOI: 10.1016/j.renene.2018.02.049

[3] Cheng X, Shi F, Liu Y, Liu X, Huang L. Wind turbine blade icing detection: A federated learning approach. Energy. 2022;**254**:124441. DOI: 10.1016/j.energy.2022.124441

[4] Zhang Y, Tang N, Niu Y, Du X. Wind energy rejection in China: Current status, reasons and perspectives. Renewable and Sustainable Energy Reviews. 2016;**66**:322-344. DOI: 10.1016/j.rser.2016.08.008

[5] Jin JY, Virk MS. Study of ice accretion along symmetric and asymmetric airfoils. Journal ofWind Engineering and Industrial Aerodynamics. 2018;**179**:240-249. DOI: 10.1016/j.jweia.2018.06.004

[6] Fu P, Farzaneh M. A CFD approach for modeling the rime-ice accretion process on a horizontal-axis wind turbine. Journal of Wind Engineering and Industrial Aerodynamics. 2010;**98**: 181-188. DOI: 10.1016/j. jweia.2009.10.014

[7] Homola MC, Virk MS, Wallenius T, Nicklasson PJ, Sundsbo PA. Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades. Journal of Wind Engineering and Industrial Aerodynamics. 2010;**98**:724-729. DOI: 10.1016/j.jweia.2010.06.007

[8] Etemaddar M, Hansen MOL, Moan T. Wind turbine aerodynamic response under atmospheric icing conditions. Wind Energy. 2014;**17**:241-265. DOI: 10.1002/we.1573

[9] Villalpando F, Reggio M, Ilinca A. Prediction of ice accretion and anti-icing heating power on wind turbine blades using standard commercial software. Energy. 2016;**114**:1041-1052. DOI: 10.1016/j.energy.2016.08.047

[10] Yirtici O, Ozgen S, Tuncer IH. Predictions of ice formations on wind turbine blades and power production losses due to icing. Wind Energy. 2019;**22**:945-958. DOI: 10.1002/we.2333

[11] Ibrahim GM, Pope K, Muzychka YS. Effects of blade design on ice accretion for horizontal axis wind turbines. Journal of Wind Engineering and Industrial Aerodynamics. 2018;**173**:39-52. DOI: 10.1016/j.jweia.2017.11.024

[12] Wang Q, Yi X, Liu Y, Ren JH, Li WH, Wang Q, et al. Simulation and analysis of wind turbine ice accretion under yaw condition via an improved multi-shot icing computational model. Renewable Energy. 2020;**162**: 1854-1873. DOI: 10.1016/j. renene.2020.09.107

[13] Manatbayev R, Baizhuma Z, Bolegenova S, Georgiev A. Numerical simulations on static vertical axis wind turbine blade icing. Renewable Energy. 2021;**170**:997-1007. DOI: 10.1016/j. renene.2021.02.023

[14] Son C, Kim T. Development of an icing simulation code for rotating wind turbines. Journal of Wind Engineering and Industrial Aerodynamics. 2020;**117**: 104239. DOI: 10.1016/j.jweia.2020. 104239

*Numerical Simulation of Icing Characteristics on a Blade Airfoil for Vertical-Axis Wind… DOI: http://dx.doi.org/10.5772/intechopen.112398*

[15] Son C, Kelly M, Kim T. Boundarylayer transition model for icing simulations of rotating wind turbine blades. Renewable Energy. 2021;**167**: 172-183. DOI: 10.1016/j. renene.2020.11.070

[16] Baizhuma Z, Kim T, Son C. Numerical method to predict ice accretion shapes and performance penalties for rotating vertical axis wind turbines under icing conditions. Journal of Wind Engineering and Industrial Aerodynamics. 2021;**216**:104708. DOI: 10.1016/j.jweia.2021.104708

[17] Reid T, Baruzzi GS, Habashi WG. FENSAP-ICE: Unsteady conjugate heat transfer simulation of electrothermal deicing. Journal of Aircraft. 2012;**49**: 1101-1109. DOI: 10.2514/1.C031607

[18] Beaugendre H, Morency F, Habashi WG. FENSAP-ICE's threedimensional in-flight ICE accretion module: ICE3D. Journal of Aircraft. 2003;**40**:239-247. DOI: 10.2514/2.3113

[19] Virk MS, Homola MC, Nicklasson PJ. Relation between angle of attack and atmospheric ice accretion on large wind turbine's blade. Wind Engineering. 2010; **34**:607-614. DOI: 10.1260/ 0309-524x.34.6.607

[20] Ahn GB, Jung KY, Myong RS, Shin HB, Habashi WG. Numerical and experimental investigation of ice accretion on rotorcraft engine air intake. Journal of Aircraft. 2015;**52**:903-909. DOI: 10.2514/1.C032839

[21] Hannat R, Morency F. Numerical validation of conjugate heat transfer method for anti/de-icing piccolo system. Journal of Aircraft. 2014;**51**:104-116. DOI: 10.2514/1.C032078

[22] Clift R, Grace JR, Weber ME. Bubbles, Drops, and Particles.

New York, New York: Academic Press; 1978. ISBN: 012176950X

[23] Kim JW, Dennis PG, Sankar LN. Ice accretion modeling using an Eulerian approach for droplet impingement. In: 51st AIAA Aerospace Sciences Meeting; 07–10 January 2013. Grapevene, Texas: AIAA; 2013. pp. 1-11

[24] Beaugendre H, Morency F, Habashi WG, Benquet P. Roughness implementation in FENSAP-ICE: Model calibration and influence on ICE shapes. Journal of Aircraft. 2003;**40**:1212-1215. DOI: 10.2514/1.2318

[25] Zhang ZX, Dong ZN. Viscous Fluid Mechanics. Beijing: Tsinghua University Press; 1999 (in Chinese)

[26] Croce G, De Candido E, Habashi WG, Munzar J, Aube MS, Baruzzi GS, et al. Analytical model for spatial and temporal evolution of in-flight icing roughness. Journal of Aircraft. 2010;**47**:1283-1289. DOI: 10.2514/1.47143

[27] Gao L, Liu Y, Hu H. An experimental investigation of dynamic ice accretion process on a wind turbine airfoil model considering various icing conditions. International Journal of Heat and Mass Transfer. 2019;**133**:930-939. DOI: 10.1016/ j.ijheatmasstransfer.2018.12.181

[28] Battisti L. Wind Turbines in Cold Climates. Green Energy and Technology; 2015. DOI: 10.1007/978-3-319-05191-8

[29] Shin J, Bond TH. Results of an icing test on a NACA0012 airfoil in the NASA Lewis icing research tunnel. In: 30th Aerospace Sciences Meeting & Exhibit; 6–9 January 1992. Reno, Washington: AIAA; 1992. pp. 1-19

## **Chapter 2**
