Numerical simulation of rime ice accretion on a three-dimensional wind turbine blade using a Lagrangian approach

Abstract

The accreted ice on wind turbine blades significantly deteriorates the blade aerodynamic performance and consequently the power production. The existing numerical simulations of blade icing have mostly been performed with the Eulerian approach for two-dimensional (2D) blade profiles, neglecting the possible three-dimensional (3D) rotating effect. This paper conducts a numerical simulation of rime ice accretion on a 3D wind turbine blade using the Lagrangian approach. The simulation results are validated through previously published experimental data. The icing characteristics along the blade radial direction are then investigated in detail. Significant radial airflow along the blade is observed, which demonstrates the necessity of 3D simulation. In addition, more droplets are found to impinge on the blade surface near the tip region, thereby producing severer ice accretion there. The accreted ice increases almost linearly along the blade radial direction in terms of both ice mass and maximum ice thickness.

Graphical abstract

Keywords

ice accretion / wind turbine / numerical simulation / Lagrangian approach

1 Introduction

Ice accretion on cold surfaces has become a concern for both academic research and engineering applications [1–7]. Wind turbine icing is one of the most representative issues, and it is typically prominent in cold climate regions. A series of adverse effects induced by blade icing issues has been demonstrated, such as power loss, aerodynamic performance degradation, and instrument or controller errors [8]. Ice accretion has a significant impact on the surface roughness and the blade shape, and this effect reduces power production [9]. For example, the Annual Energy Production loss caused by blade icing under extreme conditions is approximately 17% in Bernese Jura, Switzerland [10]. Clearly, it is important to study the characteristics of the blade icing to minimize wind turbine power loss.

When the super-cooled water droplets collide with a wind turbine blade surface and freeze upon impingement, ice accretion develops. The type of ice formed is affected by temperature, inflow velocity, liquid water content of air as well as droplet size [11–13]. In general, ice shapes are mainly categorized into two different groups: rime ice and glaze ice [14]. Rime ice appears when the temperature is below −6 °C and water droplets freeze after impingement [15]. In contrast, glaze ice is formed when the temperature ranges from −5 to 0 °C. After impingement, a fraction of water droplets freezes very rapidly and the remaining droplets run along the blade surface under the action of airflow and centrifugal force, and then freeze gradually outside the impingement area [8]. Rime ice is opaque and granular, whereas glaze ice appears transparent with horn and feather shapes [16]. It is accepted that due to the complexity of runback water motion and surface thermal energy exchange, glaze ice is comparatively harder to be predicted [17]. Characteristics of blade icing have been extensively studied using numerical simulations in previous research [9,18–21]. Virk et al. [18] studied the effect of attack angle on ice accretion near the blade tip using numerical simulation. Less severe icing is observed at lower angles of attack. In a study by Etemaddar et al. [20], the impacts of all atmospheric and geometric parameters on ice shapes as well as ice loads were estimated, and that work provided a systematic reference for related research. However, most numerical simulations have only focused on specific blade sections, instead of the distribution of ice along the blade radius. Hu et al. [9] studied the distributions of mass and thickness of ice affected by different icing parameters. The approximately linear growth trend from the blade root to the tip was found to be generally consistent with the results of previous experiments and field observations.

The 3D airflow has a certain impact on the impingement distribution and correspondingly the icing characteristics [22]. However, most of the current numerical simulations on blade icing only considered the 2D blade profile, neglecting the possible 3D rotating effect. Fu and Farzaneh [4] simulated the process of 3D rime ice accretion on a rotating blade using the Eulerian multiphase flow model. 3D ice shapes and the ice load distributions along the blade radius were obtained. Hu et al. [23] proposed a 3D icing model to simulate rime ice accretion on a rotating NREL Phase VI blade. The overall icing characteristics and power degenerations were investigated. These 3D numerical predictions all use the Eulerian approach, instead of the Lagrangian particle tracking approach, to simulate droplet motion, although the latter approach has been widely adopted in 2D simulations of blade icing.

In contrast to the continuum treatment in the Eulerian method, the Lagrangian technique tracks the trajectories of discrete particles by force analysis, with a clear mechanism and fewer assumptions. The study of particle behavior from a microscopic perspective is also applied in other fields such as sand transport [24,25], snow drifting [26,27] and cracking particles [28–31]. The 3D Lagrangian simulation has difficulties in the establishment of the model and complex calculations of particle trajectory and icing shape, as well as in carrying out the calculations due to the large computation cost. Therefore, further research using the 3D Lagrangian model is still needed.

This paper uses a Lagrangian method to simulate rime ice accretion on a 3D rotating wind turbine blade. A description of the icing model is given in detail. Then the reliability of numerical simulation is validated through the previous experiments. A comprehensive analysis of the airflow and the icing characteristics along blade radial direction is then conducted.

2 Numerical method

In this study, the numerical simulation of rime ice accretion mainly includes three steps: 1) calculation of airflow around the blade; 2) droplet trajectories and impingement calculation; 3) simulation of ice accumulation on the blade surface. A thermodynamic analysis is not considered due to the very rapid freezing of water droplets. The flow field is simulated with the commercial software ANSYS Fluent. Then the solution is used for determinations of droplet trajectories and impingement distributions based on the Lagrangian approach, which is implemented in MATLAB. Finally, the ice mass and ice shape are calculated by the ice accretion code.

2.1 Flow field calculation

The airflow around the rotating blade is intrinsically unsteady when viewed from the stationary frame. A moving reference frame (MRF) is employed to consider the blade rotation. This approach transforms the unsteady flow behavior into steady-state conditions.

The airflow is simulated as a steady and incompressible flow. The Reynolds-averaged Navier−Stokes equations for the MRF of the incompressible flow are written as follows [32]:

(1)\mathrm{\nabla }\cdot {\mathit{u}}_{\mathrm{r}}=0,

(2)\mathrm{\nabla }\cdot \left({\rho }_{\mathrm{a}}{\mathit{u}}_{\mathrm{r}}{\mathit{u}}_{\mathrm{a}}\right)+{\rho }_{\mathrm{a}}\left(\text{ω}\times {\mathit{u}}_{\mathrm{a}}\right)=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \mathit{\sigma }+\mathrm{\nabla }\cdot \mathit{\tau },

where u_r is the relative velocity of air, ρ_a is the density of air, u_a is the absolute velocity of air, ω is the angular velocity, p is the static pressure, and σ and τ represent viscous stress and Reynolds stress, respectively.

The air relative velocity u_r can be transformed using the following relation:

(3){\mathit{u}}_{\mathrm{r}}={\mathit{u}}_{\mathrm{a}}-\text{ω}\times \mathit{r},

where r is the position vector.

Airflow is predicted by solving the above governing equations for mass and momentum. Then the obtained solution is used for the calculation of water droplet motion.

2.2 Simulation of droplet trajectory

The Lagrangian and Eulerian methods are the main methods used to model particle motion [4]. Both approaches have been embedded into icing simulation solvers, such as LEWICE [1], TURBICE [33], and FENSAP-ICE [3]. In the Eulerian method, droplets are regarded as a continuum and the volume fraction concept is used. The droplet solutions are obtained by solving mass and momentum equations, which is similar to the simulation method of air phase. Nevertheless, the continuity assumption inevitably weakens the authenticity of droplet motion, and does not deal with, for example, droplet trajectory intersections [34]. In contrast, the Lagrangian method treats the droplets as discrete particles and tracks the droplet trajectories by force analysis, which has a clearer mechanism with fewer assumptions, and better accuracy as a consequence [35,37]. In this study, the Lagrangian method is used to simulate droplet trajectories and impingements.

2.2.1 Governing equation of droplet motion

The droplets are assumed herein to be spherical with equal size. In general, the principal forces acting on a water droplet are aerodynamic drag, gravity, and buoyancy. In the case of icing simulation, the effects of gravitational and buoyancy forces are negligible [38,39]. The gravity and buoyancy are much smaller than the aerodynamic drag force. In addition, the effect of droplet motion on the surrounding flow field is also neglected due to the low water concentration [34]. Based on the above assumptions, the governing equation of droplet motion can be expressed as

(4)m_{\mathrm{d}}\frac{{\mathrm{d}}^2\mathit{r}}{\mathrm{d}{t}^2}={\mathit{F}}_{\mathrm{D}}=-0.5C_{\mathrm{D}}{A}_{\mathrm{d}}{\rho }_{\mathrm{a}}\left|{\mathit{u}}_{\mathrm{d}\mathrm{r}}-{\mathit{u}}_{\mathrm{r}}\right|\left({\mathit{u}}_{\mathrm{d}\mathrm{r}}-{\mathit{u}}_{\mathrm{r}}\right),

where C_D is the drag coefficient, F_D is the drag force, m_d, A_d, and u_dr are the mass, the upwind area, and the relative velocity of water droplet, respectively.

The above governing equation can be rearranged as

(5)\frac{{\mathrm{d}}^2\mathit{r}}{\mathrm{d}{t}^2}=\frac{C_{\mathrm{D}}{R}_{\mathrm{e}}}{24}\frac{18{\mu }_{\mathrm{a}}}{{\rho }_{\mathrm{d}}{D}_{\mathrm{d}}^2}\left({\mathit{u}}_{\mathrm{d}\mathrm{r}}-{\mathit{u}}_{\mathrm{r}}\right),

where ρ_d is the water droplet density, μ _a is the air viscosity, D_d is the droplet diameter, and R_e is the relative Reynolds number defined by

(6)R_{\mathrm{e}}=\frac{\left|{\mathit{u}}_{\mathrm{d}\mathrm{r}}-{\mathit{u}}_{\mathrm{r}}\right|D_{\mathrm{d}}{\rho }_{\mathrm{a}}}{{\mu }_{\mathrm{a}}}.

The drag coefficient, C_D, is usually determined by the empirical formula. For a spherical particle in this study, the empirical drag coefficient is expressed as [40]:

(7)C_{\mathrm{D}}=\{\begin{array}{c}\frac{24}{R_{\mathrm{e}}}(1+0.1\text{5}{R}_{\mathrm{e}}^{0.6\text{87}}),R_{\mathrm{e}}⩽1000,\\ 0.44,R_{\mathrm{e}}>1000.\end{array}

2.2.2 Solving procedure

Based on the above airflow solution, the droplet trajectories and impingements can be obtained by integrating the motion equations of each droplet. The air velocity is calculated using linear interpolation from surrounding cell nodes. For each integration step, the motion equation is solved to get the new droplet position and velocity until impingement occurs, which is implemented by the fourth-order Runge−Kutta method. The specific solving procedures are as follows.

1) Through flow field calculation, the air velocities around the rotating blade are obtained.

2) The water droplets with equidistant spacing are released from the seeding location. The initial velocity of the droplet is set as the local air velocity, which can be determined by linear interpolation based on the solution of nearby airflow.

3) Equation (5) is integrated by the fourth-order Runge−Kutta scheme to gain the new droplet position and velocity of the first time step.

4) The air velocity at the new droplet position is calculated using linear interpolation.

5) Equation (5) is solved as step 3) to obtain the new droplet position and velocity of this integration step.

6) Steps 4) and 5) are repeated until the droplet impinges on the blade surface.

2.3 Determination of ice shape

It should be noted that this study focuses on rime ice accretion simulation. Rime ice prediction can be regarded as a validation of the numerical method to some extent. Thus, rime ice is selected for simulation in this study.

As for rime ice accretion, due to rapid heat dissipation, water droplets freeze immediately upon impingement. Therefore, thermodynamic analysis is omitted herein. The total collected water mass during the icing process is identical to the accumulated ice mass [4,23].

The whole blade surface is divided into some segments, and the corresponding area of each segment is marked as Δs_i, as shown in Fig.1(a). The segment is refined near the leading edge to catch the high gradient of collected water mass there. The number of droplets impinging in each segment is noted as n_i. The ice mass, m_i, namely the total collected water mass in each segment, can be calculated from

(8)m_i=\frac{n_i}{N}\cdot A\cdot V_0\cdot LWC\cdot t,

where N is the total number of water droplets at the seeding location, A is the area occupied by the seeding location, V₀ is the wind speed, LWC is the liquid water content, and t is the icing time.

For each separate segment, the ice thickness, d_i, is easily obtained by:

(9)d_i=\frac{m_i}{{\rho }_{\mathrm{i}\mathrm{c}\mathrm{e}}\mathrm{\Delta }{s}_i},

where ρ_ice denotes the ice density.

The above ice thickness is approximately considered as the ice thickness in the middle section of the blade to create the final shape of ice. For the convenience of observation, a 2D view of the middle cross-section is illustrated in Fig.1(b). Rime ice grows along the normal direction to the blade surface. The ice thickness develops perpendicularly to the corresponding segment, resulting in a discrete ice surface. To obtain a continuous surface of the final ice shape, the endpoints of the adjacent segments are then averaged, as is expressed by

(10){\left(x,y,\mathit{\text{z}}\right)}_i^{\prime }=(x,y,\mathit{\text{z}}{)}_i+\frac{1}{2}\left({\mathit{n}}_i\cdot d_i+{\mathit{n}}_{i+1}\cdot d_{i+1}\right),

where n_i and n_{i + 1} denote the unit normal vector of the ith and (i + 1)th surface segment, respectively.

The final ice surface is generated by connecting these surface points. The formed ice shape may present undesired protuberances due to the division of surface segments. Therefore, a simple curve smoothing treatment without changing the overall ice shape is adopted. The principle of the above algorithm for the determination of the ice shape refers to a subroutine of LEWICE, namely NWFOIL [1].

It is noted that ice accretion is a transient process during which the ice shape varies over time [41]. This causes variations of the airflow in the proximity of the blade, and the changed flow in turn affects the droplet impingement characteristics. For simplification, the whole icing process is simulated as a one-way steady computation without considering the growing ice shape on the surrounding airflow. Therefore, the airflow and droplet impingement solutions of the clean blade are adopted for the total icing period. This assumption can be considered reasonable for slight ice accretion with relatively small ice thickness as studied herein.

The flow chart of the overall calculation is shown in Fig.2.

3 Numerical model setup

3.1 Wind turbine blade model

Han et al. [17] tested the scaled ice accretion on a model wind turbine blade in an Adverse Environment Rotor Test Stand (AERTS). An AERTS S809 blade was designed and tested [17]. It is used as the research object of the current study. The radius of the blade is 1.372 m with an aerodynamic section of 1.016 m, and the length of the chord is 0.267 m, as illustrated in Fig.3. Here, r represents the local radius of the blade section.

3.2 Computation domain, grid and solver settings

A semi-cylindrical domain is used to save computation time, with a radius of 6 m (greater than 4R). Only a single blade is located in the middle of the domain, extending 6 m to both the upstream and downstream directions, as illustrated in Fig.4. The inlet is specified as a velocity inlet with uniform wind speed. The turbulent viscosity ratio and turbulent intensity are set as 10 and 5%, respectively. The outlet is specified as a pressure outlet and atmospheric static pressure is adopted herein. A periodic boundary is used at the bottom of the computational domain. The blade surface is defined as a no-slip wall, with the same rotation speed as the MRF. A symmetric boundary condition is applied to the lateral surface, which enables the boundary to act like a free-slip wall and to avoid the real wall effects.

An unstructured grid with tetrahedral cells is created to properly approximate the blade surface geometry complexity. The grid is refined near the blade surface. In addition, a grid-sensitivity analysis is conducted to show the grid-independence of the results. Fig.5 shows the variation in the blade aerodynamic torque with the total number of cells. For the grid with a cell number greater than 2 million, the simulation results are nearly grid-converged. Therefore, a grid with approximately 2.05 million tetrahedral cells is finally adopted, as shown in Fig.6.

The standard k−ε turbulence model is adopted in this simulation. Compared with other turbulence models, the standard k−ε turbulence model has a lower requirement of grid resolution near the boundary due to the adoption of a wall function, which significantly reduces the cell number in the near-wall region [42]. The scalable wall function is selected based on the complexity of the unstructured grid near the blade surface. The y⁺ value is found to be about 50–120. When residuals decrease to 10⁻⁴, the computation is considered to reach convergence. The detailed computational settings are shown in Tab.1, which are used for the flow field calculation of all simulation cases in this paper.

3.3 Computational parameters of droplet motion

For the Lagrangian approach, it is critical to determine the seeding location at an upstream distance where particles are to be released [43]. To save computational resources and meanwhile maintain sufficient accuracy, a series of tentative simulations are conducted to select an appropriate starting position in the free stream. The criterion for the selection of the appropriate starting position is that the simulated result of ice shape would not change if the distance from the position to the blade becomes longer, while it would change if the distance becomes shorter. About 285000 droplets are uniformly arranged at 1.5 m (greater than R) upstream of the wind turbine blade. The radius of the seeding location is 1.5 m. The initial velocity of the droplet is set as the local air velocity calculated by linear interpolation from adjacent grid points. A time-step of 3 × 10⁻⁵ s is selected in the trajectory integration to locate the exact impingement position.

4 Validation

In the experimental study of Han et al. [17], the LEWICE code was employed to validate the capacity of the AERTS facility. The rime ice shapes at the 99% radial section matched extremely well with 2D LEWICE predictions both in terms of ice thickness and the range of ice [17]. However, only 2D inflow, and icing conditions for the specific radial section, were provided.

To obtain the required free stream wind speed for given inflow conditions at the specified radial section, the blade element momentum (BEM) theory [44] is employed for an inverse calculation. The relationship between these inflow parameters can be described by the triangle illustrated in Fig.7.

The wind speed V₀ is obtained from

(11)V_{\text{0}}=\frac{V_{\mathrm{r}\mathrm{e}\mathrm{l}}\sin \phi }{1-a},

where V_rel is the relative wind speed at this section, a is the axial induction factor, and φ is the inflow angle. The axial induction factor a and the tangential induction factor a^{\prime } are derived as

(12)a={\left\{1+\frac{4F{\sin }^2\phi }{\sigma \left(C_{\mathrm{l}}\cos \phi +C_{\mathrm{d}}\sin \phi \right)}\right\}}^{-1},

(13)a^{\prime }={\left\{-1+\frac{4F\sin \phi \cos \phi }{\sigma \left(C_{\mathrm{l}}\sin \phi -C_{\mathrm{d}}\cos \phi \right)}\right\}}^{-1},

where F is the Prandtl’s tip-loss factor, σ is the local solidity, C_l is the lift coefficient, and C_d is the drag coefficient. The Glauert correction is introduced for a larger axial induction factor [44]. Given the inflow parameters of the experiment, the required free stream wind speed can be easily obtained.

Three test cases of Han et al. [17] are chosen for validation of the numerical simulation. Detailed icing conditions at the 99% radial section are given in Tab.2. The blade collective pitch angle is set as 0°. The relative wind speeds at this section for these cases are all within 6% of the desired values, i.e., 50 m/s [17]. The other icing parameters are identical to the experimental conditions adopted by Han et al. [17]. The droplet diameter is set to 20 μm, corresponding to a particle mass of 4.2 × 10⁻⁹ g. A constant ice density of 700 kg/m³ [45] is used for rime ice in the study.

As illustrated in Fig.8, the predicted ice shapes are compared against those obtained experimentally by Han et al. [17]. The simulated results are in good agreement with experimental data for all cases in the ice thickness and the ice-covered area. This indicates the reliability of the numerical prediction.

5 Results and discussion

A comprehensive analysis of icing characteristics along the blade radial direction is carried out, including the airflow behavior, droplet motion behavior, and final ice accretion. The simulated icing conditions are also chosen based on the study by Han et al. [17], as presented in Tab.3. Four sections are selected along the blade radius, which are located at r/R = 90%, 75%, 60%, 45%.

5.1 Airflow

Fig.9 shows the pressure coefficient contours and limiting streamlines on the blade suction surface. To get a better insight into the features of the airflow in the vicinity of the blade, the results of a lower wind speed (V₀ = 8 m/s) and a higher wind speed (V₀ = 16 m/s) are also presented herein. The pressure coefficient C_p is defined in the following formula:

(14)C_{\mathrm{p}}=\frac{P-P_0}{0.5{\rho }_{\mathrm{a}}\left[V_0^2+{\left(r\omega \right)}^2\right]},

where P is the static pressure and P₀ is the free stream atmospheric pressure.

It can be observed that the absolute magnitude of the pressure coefficient reduces gradually along the radial direction, particularly near the leading edge. For a speed of 12 m/s, the radial flow takes place over a large proportion of the blade surface, which indicates a 3D effect at this wind speed. The three-dimensionality is predominant, especially in the root region, and gradually weakens along the blade radius. Furthermore, as the wind speed increases, the separation line moves toward the blade tip and the radial flow extends to a larger area. For better observation, the streamlines at the 45% radial sections of the above three wind speeds are shown in Fig.10. The separation point moves from trailing to the leading edge as the speed increases. A higher wind speed corresponds to stronger flow separation. And extensive separation vortices appear at the speed of 16 m/s.

The 3D effect with the airflow along the blade radial direction can be significant. The influence of the 3D characteristics on the droplet impingement is indeed not negligible, and leads to a certain decrease in the droplet collection efficiency [22]. For wind turbine blade icing, especially the numerical predictions, the 2D BEM methodology [44] is widely adopted to simulate flow behavior in the vicinity of the blade, due to the ease of implementation. This approach evaluates the inflow conditions at specific blade radial sections, neglecting the impact of the possible 3D radial flow. Compared with 2D BEM, the 3D numerical simulation takes full account of the 3D rotating effect, and can better capture the details of the airflow [46]. As a consequence, the predictions of droplet impingement and resultant ice accumulation tend to be more accurate.

5.2 Droplet motion

To better understand the droplet motion behavior at different radial sections, the droplet trajectories are analyzed in detail, as illustrated in Fig.11. Significant differences in droplet trajectories are observed as the local radius changes. In addition to the variation of the travel direction, droplets tend to impinge on the blade surface near the tip region. Fig.12 shows a comparison of the collected droplet distributions around the four sections. The total number of impinged droplets in the 90% region is more than twice that in the 45% region. This trend is mainly due to the variation of airflow behavior along the blade, causing a higher wind speed and smaller attack angle at the blade tip.

5.3 Ice distribution

Fig.13 illustrates the ice shapes of corresponding radial sections. Ice mainly accretes at the leading edge and the pressure side of the blade surface due to the inflow direction. Moreover, the ice thickness increases significantly along the radial direction, which is generally consistent with previous studies [19,23]. Fu and Farzaneh [4] attributed this trend to the blade tapered structure with smaller airfoil size near the tip and the blade rotating movement. As for the untapered and untwisted blade in this study, the increase in relative wind speed along the blade radius caused by the rotating movement is the main cause. As aforementioned, the higher relative wind speed allows the blade to collect more droplets near the tip region, giving rise to a larger ice growth rate. Therefore, severer ice accretion usually forms at this location. The conclusion is consistent with that in a previous experiment [15].

For a more intuitive comparison, the ice mass and the maximum ice thickness at different radial sections are normalized by those at the 45% radial section, as shown in Fig.14. An approximately linear increasing trend of the accreted ice along the blade radial direction can be obtained, in terms of both the ice mass and the maximum ice thickness. Furthermore, the growth rate of the ice mass is greater than that of the maximum ice thickness, which can be explained by the increase of the ice-covered area along the blade. The increments between the 45% and 90% radial sections reach up to 113.3% in the maximum ice thickness and 136.7% in the ice mass, respectively.

6 Conclusions

This paper conducts a numerical study of rime ice accretion on a 3D rotating wind turbine blade using the Lagrangian approach. The reliability and accuracy of the numerical method are successfully proved by comparisons with experimental data. The icing characteristics along the blade radial direction are analyzed in detail. The main conclusions can be drawn as below.

1) The 3D effect with the radial airflow along the blade can be prominent at specific high wind speeds, which may not be accurately predicted with 2D simulations. And it is highly recommended that 3D numerical simulation be applied in wind turbine blade icing for further study.

2) The airflow separation point moves from the trailing to the leading edge as the speed increases. A higher wind speed corresponds to stronger flow separation. And extensive separation vortices appear on the suction at the speed of 16 m/s.

3) Significant differences in droplet trajectories are observed as the local radius changes. In addition, droplets tend to impinge on the blade surface near the tip region.

4) The higher relative wind speed near the tip region allows the blade to collect more droplets, thereby producing severer ice accumulation there. The increase of the accreted ice at the 90% radial section corresponds to 113.3% in ice thickness and 136.7% in mass compared with their values at the 45% radial section.

In the future, glaze icing issues on wind turbine blade affected by the 3D rotating effect could be simulated under the Lagrangian reference frame. Additionally, the influences of dynamic change in ice shape will also be taken into account.

References

[1] RuffG ABerkowitzB M. User’s Manual for the NASA Lewis Ice Accretion Prediction Code (LEWICE), 1990

[2] BraggM BHutchionTMerretJOltmanRPokhariyalD. Effect of ice accretion on aircraft flight dynamics. In: Proceedings of the 38th Aerospace Sciences Meeting and Exhibit. Reno, NV: American Institute of Aeronautics and Astronautics, 2000

[3] Beaugendre H, Morency F, Habashi W G. FENSAP-ICE’s three-dimensional in-flight ice accretion module: ICE 3D. Journal of Aircraft, 2003, 40(2): 239–247

[4] 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(4−5): 181–188

[5] Zhang Z, Jiang X, Sun C, Hu J, Huang H, Gao D W. Influence of insulator string positioning on AC icing flashover performance. IEEE Transactions on Dielectrics and Electrical Insulation, 2012, 19(4): 1335–1343

[6] Zhang Z, Zhang Y, Jiang X, Hu J, Hu Q. Icing degree characterization of insulators based on the equivalent collision coefficient of standard rotating conductors. Energies, 2018, 11(12): 3326

[7] Jin J Y, Virk M S. Study of ice accretion and icing effects on aerodynamic characteristics of DU96 wind turbine blade profile. Cold Regions Science and Technology, 2019, 160: 119–127

[8] Yirtici O, Tuncer I H, Ozgen S. Ice accretion prediction on wind turbines and consequent power losses. Journal of Physics: Conference Series, 2016, 753: 022022

[9] Hu L, Zhu X, Hu C, Chen J, Du Z. Wind turbines ice distribution and load response under icing conditions. Renewable Energy, 2017, 113: 608–619

[10] Barber S, Wang Y, Jafari S, Chokani N, Abhari R S. The impact of ice formation on wind turbine performance and aerodynamics. Journal of Solar Energy Engineering, 2011, 133(1): 011007

[11] Liu Y, Hu H. An experimental investigation on the unsteady heat transfer process over an ice accreting airfoil surface. International Journal of Heat and Mass Transfer, 2018, 122: 707–718

[12] Jin J Y, Virk M S. Study of ice accretion along symmetric and asymmetric airfoils. Journal of Wind Engineering and Industrial Aerodynamics, 2018, 179: 240–249

[13] GaoLLiuYZhouWHuH. An experimental study on the aerodynamic performance degradation of a wind turbine blade model induced by ice accretion process. Renewable Energy, 2019, 133: 663−675

[14] Makkonen L. Models for the growth of rime, glaze, icicles and wet snow on structures. Philosophical Transactions of the Royal Society of London Series A: Mathematical, Physical, and Engineering Sciences, 1776, 2000(358): 2913–2939

[15] GaoLLiuYHuH. 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

[16] ReidTBaruzziGOzcerISwitchenkoDHabashiW. FENSAP-ICE simulation of icing on wind turbine blades, part 1: performance degradation. In: 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. Grapevine, TX: American Institute of Aeronautics and Astronautics, 2013

[17] Han Y, Palacios J, Schmitz S. Scaled ice accretion experiments on a rotating wind turbine blade. Journal of Wind Engineering and Industrial Aerodynamics, 2012, 109: 55–67

[18] Virk M S, Homola M C, Nicklasson P J. Relation between angle of attack and atmospheric ice accretion on large wind turbine’s blade. Wind Engineering, 2010, 34(6): 607–613

[19] Virk M S, Homola M C, Nicklasson P J. Atmospheric icing on large wind turbine blades. International Journal of Energy and Environment, 2012, 3(1): 1–8

[20] Etemaddar M, Hansen M O L, Moan T. Wind turbine aerodynamic response under atmospheric icing conditions. Wind Energy, 2014, 17(2): 241–265

[21] Ibrahim G M, Pope K, Muzychka Y S. Effects of blade design on ice accretion for horizontal axis wind turbines. Journal of Wind Engineering and Industrial Aerodynamics, 2018, 173: 39–52

[22] Zhu X, Hu L, Chen J, Shen X, Du Z. Calculation of collection efficiency on NREL Phase VI blade. Journal of Energy Resources Technology, 2018, 140(7): 071202

[23] Hu L, Zhu X, Chen J, Shen X, Du Z. Numerical simulation of rime ice on NREL Phase VI blade. Journal of Wind Engineering and Industrial Aerodynamics, 2018, 178: 57–68

[24] Kang L, Guo L. Eulerian−Lagrangian simulation of aeolian sand transport. Powder Technology, 2006, 162(2): 111–120

[25] Huang N, Xia X, Tong D. Numerical simulation of wind sand movement in straw checkerboard barriers. European Physical Journal E, 2013, 36(9): 99

[26] Huang N, Sang J, Han K. A numerical simulation of the effects of snow particle shapes on blowing snow development. Journal of Geophysical Research, 2011, 116(D22): D22206

[27] Okaze T, Niiya H, Nishimura K. Development of a large-eddy simulation coupled with Lagrangian snow transport model. Journal of Wind Engineering and Industrial Aerodynamics, 2018, 183: 35–43

[28] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[29] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computational Methods in Applied Mathematics, 2010, 199: 2437–2455

[30] Zhang Y, Zhuang X. Cracking elements method for dynamic brittle fracture. Theoretical and Applied Fracture Mechanics, 2019, 102: 1–9

[31] Zhang Y, Yang X, Wang X, Zhuang X. A micropolar peridynamic model with non-uniform horizon for static damage of solids considering different nonlocal enhancements. Theoretical and Applied Fracture Mechanics, 2021, 113: 102930

[32] Abdelsalam A M, Ramalingam V. Wake prediction of horizontal-axis wind turbine using full-rotor modeling. Journal of Wind Engineering and Industrial Aerodynamics, 2014, 124: 7–19

[33] Makkonen L, Laakso T, Marjaniemi M, Finstad K J. Modelling and prevention of ice accretion on wind turbines. Wind Engineering, 2001, 25(1): 3–21

[34] WidhalmMRonzheimerAMeyerJ. Lagrangian particle tracking on large unstructured three-dimensional meshes. In: The 46th AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV: American Institute of Aeronautics and Astronautics, 2008

[35] VilledieuPTrontinPGuffondDBoboD. SLD Lagrangian modeling and capability assessment in the frame of ONERA 3D icing suite. In: 4th AIAA Atmospheric and Space Environments Conference. New Orleans: American Institute of Aeronautics and Astronautics, 2012

[36] Sang W, Shi Y, Xi C. Numerical simulation of icing effect and ice accretion on three-dimensional configurations. Science China. Technological Sciences, 2013, 56(9): 2278–2288

[37] Li Y, Sun C, Jiang Y, Yi X, Xu Z, Guo W. Temperature effect on icing distribution near blade tip of large-scale horizontal-axis wind turbine by numerical simulation. Advances in Mechanical Engineering, 2018, 10(11): 1–13

[38] PootsG. Ice and Snow Accretion on Structures. Hoboken: Wiley, 1996

[39] Shu L, Liang J, Hu Q, Jiang X, Ren X, Qiu G. Study on small wind turbine icing and its performance. Cold Regions Science and Technology, 2017, 134: 11–19

[40] BourgaultYHabashiW GDompierreJBoutaniosZBartolomeoW D. An Eulerian approach to supercooled droplets impingement calculations. In: The 35th Aerospace Sciences Meeting and Exhibit. Reno: American Institute of Aeronautics and Astronautics, 1997

[41] PedersenMSørensenH. Towards a CFD model for prediction of wind turbine power losses due to icing in cold climate. In: The 16th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery. Honolulu, HJ, 2016

[42] Homola M C, Virk M S, Wallenius T, Nicklasson P J, Sundsbø P A. Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades. Journal of Wind Engineering and Industrial Aerodynamics, 2010, 98(12): 724–729

[43] WirogoSSirirambhatlaSLeaderT. An Eulerian method to calculate the collection efficiency on two and three dimensional bodies. In: The 41st Aerospace Sciences Meeting and Exhibit. Reno: American Institute of Aeronautics and Astronautics, 2003

[44] HansenM O L. Aerodynamics of Wind Turbines. Oxford: Earthscan Publication, 2008

[45] Seifert H, Westerhellweg A, Kröning J. Risk analysis of ice throw from wind turbines. Boreas, 2003, 6(9): 2006–01

[46] Bai C, Wang W. Review of computational and experimental approaches to analysis of aerodynamic performance in horizontal-axis wind turbines (HAWTs). Renewable & Sustainable Energy Reviews, 2016, 63: 506–519

RIGHTS & PERMISSIONS

Higher Education Press