Study of Ice Accretion on Horizontal Axis Wind Turbine Blade Using 2D and 3D Numerical Approach

In order to optimize the wind turbine operation in ice prone cold regions, it is important to better understand the ice accretion process and how it affects the wind turbine performance. In this paper, Computational Fluid Dynamics (CFD) based 2D and 3D numerical techniques are used to simulate the airflow/droplet behaviour and resultant ice accretion on a 300 kW wind turbine blade. The aim is to better understand the differences in the flow behaviour and resultant ice accretion between both approaches, as typically the study of ice accretion on the wind turbine blade is performed using simple 2D simulations, where the 3D effects of flow (air & droplet) are ignored, which may lead to errors in simulated ice accretion. For 2D simulations, nine sections along a 300 kW wind turbine blade are selected, whereas for 3D study, a full-scale blade is used. The obtained results show a significant difference in the ice accretion for both approaches. Higher ice growth is observed in 2D approach when compared with the full-scale 3D simulations. CFD simulations are carried out for three different icing conditions (typical, moderate and extreme), in order to estimate the extent of differences the different operating conditions can introduce on the ice accretion process in the 2D and the 3D simulations. Complex ice shapes are observed in case of extreme ice conditions, which affect the aerodynamic performance of the blade differently from typical and moderate ice conditions.


I. INTRODUCTION
The primary appeal for the wind energy in cold climate regions are the higher air densities and wind speeds, which provide significant benefits to power production. However, the low temperatures and the risk for atmospheric icing associated with them poses additional challenges for wind turbine systems, when compared with the warmer climate [1]. According to the WindPower Monthly report, 127GW of installed wind power capacity was located in cold climate regions at the end of 2015, with an expected total installed capacity reaching 186GW in 2020, which would be a remarkable 30% share of the global forecasted wind capacity [2]. However, atmospheric ice accretion can be one of the primary The associate editor coordinating the review of this manuscript and approving it for publication was Huiqing Wen . hazards for the design and safety of wind turbine [3], [4]. This highlights the need to better understand the ice accretion process on wind turbine blades with the aim to improve safety and to reduce the Capital Expenditure (CAPEX) and Operational Expenditure (OPEX) related to wind turbine operations in the cold regions.
Atmospheric ice accretion on wind turbine blade mainly occurs due to the impingement of super-cooled water droplets, which may freeze on blade surface [5]. The accreted ice changes the blade geometry and increases its surface roughness, thus altering the aerodynamic performance of wind turbine blade [6], [7]. This primarily results in the decrease of lift force and increase of drag force, which leads to reduction in the aerodynamic performance of blade and overall power production. The accreted ice shape along the blade's leading edge depends upon many variables, such as the geometry of wind turbine blade and operating conditions, most important being the relative wind velocity, atmospheric temperature, droplet size, droplet distribution spectrum and liquid water content [8].
The ice accretion process and its effects on wind turbine's aerodynamic performance can be studied using both the experimental and the numerical approaches. Along with the field measurements, the advanced multiphase CFD simulations have become an important tool for simulating and determining the performance of wind turbine blade under icing conditions [9]- [14]. In 2011, Barber and Jafari [15] analyzed the ice accretion process on the wind turbine blades and its effects on the turbine performance using both CFD and the field measurement data. They found that the ice accretion at higher altitude sites can cause an AEP lose up to 17% and the outer section (95 − 100% of the blade length) are the most severely iced. Moreover, Lamaraoui and Benoit [16] found that the majority of power production loss due to atmospheric icing is primarily happening at the 80% of the blade section close to tip. Thus, the tip section of the blade is typically used in studying the atmospheric ice accretion and its effects on the wind turbine performance.
The existing codes for ice accretion such as LEWICE and CANICE [17] were developed strictly for aviation with redundant features analyses, while FENSAP-ICE is a secondgeneration, state-of-the-art CFD icing code capable of 2D and 3D icing simulations for a large variety of applications. CFD simulations of atmospheric icing on wind turbine blades can be performed using 2D blade ''planes'' and 3D full blade approach. FENSAP has been tested for icing on rotating components such as helicopter rotors and engines. This highlights the possibility to use this code for wind turbine rotor blade also. Modern CFD tools such as LEWICE [18]- [20], CANICE-WT [21], TURBICE [22], and FENSAP-ICE [23], [24] are capable for modelling the icing using both 2D and 3D techniques of wind turbine. The differences of these software including operational environment, prolonged exposure of the wind machine to icing, and temporal variations in the external meteorological parameters. Etemaddar [25] have studied the effects of atmospheric and system parameters on the 2D icing simulated with the LEWICE code for different sections along the blade and estimated aerodynamic coefficients using FLUENT after ice accretion, following with the validation of obtained results against experimental results of NREL 5MW wind turbine. Switchenko et al. have performed a series of 2D and 3D simulations of complex wind turbine icing events in Canada, using the FENSAP-ICE and compared with the experimental data [26]. Study of ice accretion on wind turbine blade shows that smaller wind turbines can be more adversely affected by the atmospheric icing than the larger ones.
Currently there are large number of wind farms are being built at high-altitude mountainous areas. When compared to the European sites, the icing characteristics of the Chinese sites such as the accreted ice shapes, ice density and its adhesion to the surface is quite different, mostly due to different climatic characteristics [12]- [14]. Thus, studying the ice accretion in such regions presents a unique opportunity. Chongqing University (CQU) is one of the major research institute, which has carried out various studies on the atmospheric ice accretion and its effect on the wind turbine operations in China, especially due to the access to the Xufeng Mountain natural icing test station [15], [16], [18], [19]. The researchers at CQU have carried out the CFD simulation of airflow behaviour along full scale wind turbine blade using 3D approach using ANSYS-FLUENT and have found that higher wind speeds lead to more abrupt flow separations, especially for glaze ice [27]. They have not performed any icing simulations, rather used ice template shapes along wind turbine blade and simulated the airflow behaviour to estimate the aerodynamic performance. In order to extend the results from the previous CQU study, the focus of this paper is to better understand the differences in the ice accretion process on the CQU 300kW wind turbine using both 2D and 3D numerical techniques. For the 2D analysis, nine sections have been selected along the blade, whereas for 3D analysis the entire blade is used in full scale.

II. WIND TURBINE SPECIFICATIONS
The work in this paper is performed using 300kW wind turbine blade. This wind turbine is installed at Xufeng Mountain natural icing test station, which is situated in the Southwest Hunan Province of China (27.1258 • N, 110.5606 • E). The wind turbine is located on a complex terrain at the average elevation of 1400 m.a.s.l. at the top of the Xufeng Mountain, which has an annual average precipitation of 1810 mm that belongs to the typical micro-topography icing prone region [28]. Figure 1 shows the terrain and the installation location of CQU 300kW wind turbine, where typically more than 90 days annually of ice cover and high humidity are observed.

III. DESIGN OF EXPERIMENT
Typical procedure for the numerical modelling of atmospheric ice accretion on wind turbines is to perform the CFD simulation on the 2D blade profile sections (airfoils), primarily to estimate its aerodynamic performance. The main advantage of 2D approach is the simplicity and computational efficiency of the setup. However, the main drawback of the 2D simulations is its inability to provide information about the flow in the z-direction, which can be crucial if more ''physics-rich'' processes are expected to occur that may influence the results. Therefore, this study investigates the potential differences in the ice accretion process on the CQU 300 kW wind turbine using both 2D and the 3D CFD simulations. The design of experiment section consists of following three parts: geometry, mesh and the CFD setup.

A. GEOMETRY MODEL
For simplicity, it is assumed that the structural model is limited to a surface representation of the wind turbine blade VOLUME 8, 2020 FIGURE 1. Overview of Xufeng mountain nature icing test station and 300kW wind turbine. [26]. and the surface in question is assumed to be composed from a collection of airfoil shapes that are lofted in the blade axis directions. In order to extend the previous CQU research, [28] the focus is on the full scale blade geometry and the chosen nine cutting planes, situated from root to tip in the ascending order 14.6 m long blade is attached to the hub with a radius of 3 m, which gives the total rotor radius of 17.6 m. The on-site observations at the Xufeng Mountain icing test station show that the outer section of the wind turbine blade near tip often has the largest ice thicknesses, thus, the majority of the chosen planes in this 2D study are located close to the tip section of the blade. The blade is composed of 15 airfoils which are shown in Figure 2. The main geometric parameters for these nine chosen sections are summarized in Table 1.

B. MESH
Numerical mesh is generated using Pointwise for 2D airfoils and the ANSYS Workbench for 3D blade. Fine meshes are made in order to accurately determine the boundary layer characteristics (shear stresses and heat fluxes) of the turbine blade. Mesh sensitivity study was carried out using coarse, medium and fine meshes to accurately determine the boundary layer characteristics (shear stress and heat flux). During mesh sensitivity analysis, number of mesh elements and y+ value less than 1 for first cell layer was selected based upon the heat flux calculations, where a numerical check was imposed that the heat flux computed with the classical formulae dT/dn should be comparable with the heat flux computed with the Gresho's method. Mesh sensitivity study showed that the effect of mesh size on droplet solution was negligible, however some flow quantities including convective heat flux on the blade surface was sensitive to the mesh size. 2D airfoil meshes use the C-type structured numerical grid while 3D mesh of full-scale blade is generated using unstructured tetrahedral mesh. For 2D mesh the total number of hexahedral elements, ranges from 75761 to 84056 and approximate 4.88 million grid cells were used for fullscale 3D analysis. The examples of these meshes are shown in Figures 3 and 4.

C. CFD MODEL
CFD based multiphase numerical simulations are carried out using ANSYS FENSAP-ICE, which uses the finite element method Navier-Stokes analysis with Eulerian water droplet impingement solver to analyses the steady flow field around the turbine blade. Atmospheric ice accretion on a wind turbine blade can be numerically simulated by means of integrated thermo-fluid dynamic models, which involve the fluid flow simulation, droplet behavior, surface thermodynamics and phase changes. Airflow behavior is simulated by solving the nonlinear partial differential equations for the conservation of mass, momentum and energy.
where ρ is the density of air, v is the velocity vector, subscript α refers to the air solution, T refers to the air static temperature in Kelvin, σ ij is the stress tensor and E and H are the total initial energy and enthalpy, respectively. The sand grain roughness for the iced surface is calculated using following Shin and al. roughness model.  Two-phase flow (air and water droplets) is numerically simulated using the Eulerian approach, where the super cooled water droplets are assumed to be spherical. The Eulerian two-phase fluid model consists of the Navier-Stokes equation with the water droplets continuity and momentum equation. The water droplet drag coefficient is based on the empirical correlation for the flow around the spherical droplets described by Clift and Weber [29] ∂α ∂t where α is the water volume fraction, Vd is the droplet velocity, CD is the droplet drag coefficient and Fr is the Froude number. Surface thermodynamics is calculated using VOLUME 8, 2020 the mass and energy conservation equations, considering the heat flux due to convective and evaporative cooling, heat of fusion, viscous and kinetic heating. (8) Equation 9 expresses the conservation of energy: The coefficients ρ f , c f ,c s , σ, ε,L evap , L fusion are physical properties of the fluid. The reference conditions T ∞ , V ∞ , LWC are the airflow and droplets parameters. Equation 10 express the Jones (glaze) formula: where d is the droplet diameter in microns, V d is the droplet impact velocity (m/s) andT wall is the wall temperature (Celsius). ALE (Arbitrary Langrangian Eulerian) formulation is used for the mesh displacement during ice accretion. This approach adds the grid speed terms to the Navier-Stokes equations to account for the mesh velocity [30].

IV. RESULTS AND DISCUSSION
Ice accretion along wind turbine blade changes its geometric shape, which effects the flow behaviour along pressure and suction sides of the blade and results a change in its aerodynamic coefficient and overall performance. In this study, CFD based numerical analysis are carried out to simulate the airflow behavior, droplet behaviours and ice shapes using CQU clean and iced single airfoils (2D) and full scale blade (3D). Figure 5 shows the velocity vectors and the pressure coefficient distribution for the full-scale blade in the 3D numerical simulations, while Figures 6 shows the comparison of velocity magnitudes and the pressure coefficients on the nine chosen planes in the 2D and the 3D CFD simulations, respectively. The 3D planes in Figure 6 are extracted from full-scale 3D simulations of the blade using plane-cutting option in CFD post processing. A clear difference in the airflow behaviour is observed in both cases.

A. AIRFLOW BEHAVIOR COMPARISON
The first immediate observation when comparing the results between 2D and 3D simulations is that the 2D  simulations have significantly higher values of the velocity magnitudes. This implies that the flow separation in the 2D case occurs more abruptly than in the 3D case, which results into sharp pressure gradients. Moreover, there is a significant change in the vorticity patterns, occurring behind the trailing edge when comparing the results for the velocity magnitudes between 2D and 3D case. In the 3D case, the vorticity is present in all nine planes of interest, with the sections, closer to the root of the blade having a non-zero value of the velocity magnitudes in said vortices, which gradually reaches values close to zero in the section closer to the tip. In comparison, the results from the 2D simulation show more extreme vorticity in sections closer to the root (Plane 2, 4 and 6) with velocity magnitudes being close to zero, while the sections closer to the tip do not have vorticity in the wake at all. This suggests that the flow velocity magnitude in the 3D case has a significant non-zero component in the z-direction along the blade, which can cause complex fluid-structure interactions and as a result a more complicated flow pattern. The same sort of interaction is absent from the 2D cases simply by virtue of them having no extrusion in the z-direction. The primary explanation for this difference is the asymmetry of the blade in the 3D, mainly due to twist angle between different sections of the blade, such as Plane 8 have no vortex in 2D case but in 3D case it shows vorticity pattern.
Thus, based on 2D and 3D separation study of aerodynamic characteristics, the comparison of pressure coefficient (Cp) for both 2D-airfoil and 3D-plane cases have carried out in Figure 7 for clean and iced conditions. The good agreement between single airfoil and single blade cutting planes are found in most cases, especially in leading edge. Reversing point are carried out more forward in single blade's most cutting planes, due to the lift and drag forces have effects the flow behaviours which affect Cp, however, the uncommon trends is found in the 3D mid section ( Plane 6, 8,10 and 12) because of flow effects and wake effects. Compared with 3D simulation results, 2D results have smaller Cp range because the windward and leeward have different lift and drag force due to the shadow regions and aerodynamic changes. Icing types leads a big difference in 2D simulation, however for the 3D simulation the difference not obvious. It may be the reason that the flow behaviours are difference from 2D and 3D due to the higher and lower artificial viscosity respectively. Figure 8 shows the droplet collision efficiency comparison for the moderate, typical and extreme icing conditions in the 2D and the 3D simulations. Since the droplet impingement simulations have been performed with the Langmuir D distribution, which is a seven-bin distribution, only the spectrum-averaged results are presented, in order to streamline the comparison.

B. DROPLET BEHAVIOR COMPARISON
Since the operating conditions for the different icing types in Table 2 are constant with the exception of the MVD and LWC values, the differences in the droplet behaviour will be primarily influenced by the droplet inertia, which depends on the airflow velocity magnitude distributions and the MVD values as a function of parameter Vd 2 where 'd' is the droplet diameter. Thus, the preceding discussion on the differences in the velocity magnitudes between the 2D and the 3D numerical simulations, is also relevant for the droplet collision efficiencies comparison, as the higher velocity magnitudes in the 2D case will directly result in the increase in the droplet inertia, which will in turn will increase the overall and local collision efficiencies, maximum impingement angles and the droplet impact velocities. This effect in terms of the local collision efficiencies and the maximum impingement angles is clearly visible in Figure 8, where the results from 2D simulations displays higher values of the local collision efficiencies and the maximum impingement angles, when compared with the 3D results. This will directly result in more accreted ice masses, thicker ice shape, covering larger area and higher accreted ice densities, as they are a function of the droplet inertia in the Jones (glaze) icing density formulation, in the 2D simulations. These conclusions are applicable for all three tested icing typestypical, moderate and severe, as the 2D simulations show consistently higher values of the local collision efficiencies and the maximum impingement angles for all planes and at all icing types. The only thing which will change with the change in icing type is the ''ratio'' of the droplet collision efficiencies between the 2D and the 3D  simulations, however, due to non-linear dependence of the droplet collision efficiencies on the droplet inertia the exact ''ratio'' will also be highly non-linear. Tables 3 and 4 show the comparison of accreted ice densities and the maximum ice thicknesses for the nine planes of interest in the 2D and the 3D simulations, while Figure 9 shows the ice accretion on full scale blade for three different icing conditions, whereas Figure 10 provides the comparison of the accreted ice shapes after 60 minutes of ice accretion. The primary observation from the Table 3 is that for the typical icing conditions the ice density values are higher in the 3D simulations than in 2D. For the rest of the cases the values of accreted ice densities are equal at 840 kg/m 3 . The difference in values for the typical icing conditions is peculiar, considering that from the comparison of the droplet collision efficiencies it follows that in 2D cases the droplet collision efficiencies and the maximum impingement angles are higher than in the 3D simulations, thus implying that the droplet inertia is higher in 2D simulations and by extension, the droplet impact velocities should also be higher, which would result in higher accreted ice densities. However, the higher values of the accreted ice densities for the typical icing conditions in the 3D cases suggest that the surface temperature must be lower in the 2D cases possibly as a result of a difference in the heat fluxes between the 2D and the 3D simulations. When it comes to the maximum ice thicknesses from Table 4 and the accreted ice shapes from Figure 10, it is obvious that the CFD simulations in the 2D cases accrete significantly more ice than the 3D simulations for all tested icing conditions. This coupled with the fact that the velocity magnitude and pressure coefficient distributions imply that the airfoils in 2D may experience higher magnitude of lift and drag forces than 3D simulations, and may results in the airfoils in the 2D simulations being more significantly affected by the accreted ice.

V. CONCLUSION
Within this study, a series of multiphase CFD simulations on a 300 kW wind turbine blade have been performed using 2D and 3D numerical simulations. The 2D simulations were performed on the nine blade sections of interest, while the 3D simulations were performed on the entire wind turbine blade in full scale. The primary objective of this paper was comparison of the multiphase flow behaviour between the 2D and the 3D CFD simulations, primarily in connection with the atmospheric ice accretion. The obtained results show that in 3D simulations there is a significant velocity component and flow interactions in the z-direction, which are obviously absent from the 2D simulations. As a result of this flow interaction in the third dimension, the velocity magnitudes are reduced in the 3D simulations when compared to the 2D simulations. This, in turn, affects the ice accretion process, as the higher velocity magnitudes in the 2D cases result in the higher droplet inertia, collision efficiencies and the maximum impingement angles, which results in more ice mass accreted along with the thicker and larger ice shapes present in the 2D simulations. Coupled with the airflow velocity magnitudes, the result suggest that the airfoils in the 2D simulations are more aerodynamically loaded than the full blade in 3D, and the higher ice accretion values will affect the aerodynamic performance of the airfoils more significantly in the 2D simulations when compared to the 3D results. The primary reason behind such differences is believed to be the asymmetry of the wind turbine blade in 3D, primarily due to the twist angle, which causes the complex flow patterns to arise in the first place. MUHAMMAD SHAKEEL VIRK received the bachelor's degree in aerospace engineering from NUST Pakistan, and the master's degree in aerospace engineering and the Ph.D. degree in computational fluid dynamics, U.K. He is currently a Professor of cold climate technology with the University of Tromsø (UiT), Norway. He established the research activities and infrastructure related to atmospheric icing, cold climate technology, computational fluid dynamics (CFD), and wind energy. With his experience of renewable energy and CFD based numerical modeling of complex turbulent flow behaviors and heat transfer, he is also having a good understanding of the aircraft design and aspects related to high-speed aerodynamics from his aerospace engineering background. He is also leading the Arctic Technology and Icing Research Group, UiT. He is having vast academic and research network and actively working with many universities/research centers of Europe, North America, and Asia. He is also engaged in various applied research projects co-funded by industry and various national and European funding agencies. He has successfully managed and participated as principal and co-investigator in various international research projects funded by European Union, Norwegian Research Council, and industry. He is the author of 140 scientific publications peer-reviewed journal publications and international conference proceedings. He has received numerous scientific prizes, including best paper awards and the EU Arctic Award 2018 from European Union, for one of his research project focused on wind energy in cold climate.