REVIEW OF STUDIES ON THE CFD-BEM APPROACH FOR ESTIMATING POWER LOSSES OF ICED-UP WIND TURBINES

However, the methodology was not adequately presented and discussed for wind turbine icing. This paper reviews the results of almost all the up-to-date published papers that approached this method, summarizing the findings and federates the research in that field to conclude with concrete facts and details that advance research in this domain.

The application of computational fluid dynamics (CFD) in wind turbine design and analysis is becoming increasingly common in research on wind energy, resulting in a better knowledge of the aerodynamic behaviour of rotors. Due to the deformation of blade airfoils on account of icing, a significant drop in aerodynamic performance brings wind turbines to lose considerable portions of their productivity. Estimating power degradation due to icing via 3D simulation, although it is essential to capture the three-dimensional turbulence effects, is very costly in computational resources despite technological development; it then becomes unfeasible when it comes to different operation scenarios to estimate icing originated power losses. The Quasi-3D simulation based on the CFD-BEM method is a practical alternative for generating wind turbines' power curves. It showed effectiveness in predicting performance up to a certain level. More than few studies in the literature have adopted this approach to generate the power curve for both clean (un-iced) and iced-up wind turbines. However, the methodology was not adequately presented and discussed for wind turbine icing. This paper reviews the results of almost all the up-to-date published papers that approached this method, summarizing the findings and federates the research in that field to conclude with concrete facts and details that advance research in this domain.

…………………………………………………………………………………………………….... Introduction:-
Given the reality of global warming, the energy transition to renewable energy is gradually increasing to meet increased energy needs and reduce greenhouse gas emissions. Wind energy is an essential source of renewable energy that is increasingly used to generate electricity instead of or in addition to other sources of energy. Over the past decade, the global production capacity of wind power has grown at an estimated 28% annually [14]. However, this technology is highly dependent on the geographic location of its uses.
Cold regions are well known for their high wind energy potential and low population density, but they are also subject to frequent icing events. Icing determines higher loads, fatigue, vibrations and energy losses. To mitigate these effects, additional ancillary equipment (anti-icing and/or de-icing systems) are required [12]. Commercial and prototype level sensors, blade ice protection (anti-icing and de-icing) systems, and other solutions for icing conditions are available on the market. However, only limited information about their effectiveness has been published, and different systems may be better suited for different conditions [13].

634
Icing is a serious problem in icy climates. However, icing problems are common and have severe consequences on wind turbine production, maintenance and service life. The accumulation of ice changes the shape and the roughness of wind turbine blades, alters their aerodynamic properties, and significantly reduces the resulting energy production. In order to evaluate icing consequences and optimize the operation of Ice Protection Systems, research on wind turbine icing depends more on numerical simulations than experimental methods. Computer simulation can help predict and analyze the iced profiles of wind turbine blades, thus calculating the aerodynamic degradation due to icing and the energy necessary for de-icing under specific weather conditions. Simulation is an essential means of compliance that is not so expensive as the experimental tests. With highly accurate, robust codes, it could be possible to wean off expensive wind tunnel tests.
The elementary study for estimating the impact of ice accretion on wind turbines is the two-dimensional simulation of the blade's airfoil sections in limited arrays of meteorological and operational parameters. The several published studies cited on wind turbine icing simulation were primarily up to quite recently dependent on 2D simulations for limited parameter analysis and validation, not so much for design and optimization, nor for the optimization of ice prevention systems. As a principal objective to estimate the drop-in power production during any icing conditions, several kinds of simulation have to be conducted, and the 2D simulation will no longer fulfill the purpose. The aerodynamic performance has to be estimated in icing conditions at each blade section for every angle of attack, angle of pitch, rotational speed and other relevant parameters. A typical horizontal-axis wind turbine (HAWT) has twisted blade geometry of usually multiple forms of sections. This special form of blade geometry requires exclusively 3D simulations. However, icing affects not only the blade; it affects the entire wind turbine with various types of 3D effects, such as the effect of flowing air along the radial direction [1].Moreover, the rotation of the blades will also affect ice accretion. Therefore, there is a need for a full 3D simulation of the rotating turbine. Even a simulation for the entire wind farm is also important to consider the effect of the aerodynamic wakes behind the turbines and to account, with adequate accuracy, for production losses during specific icing conditions. Threedimensional numerical simulation of ice accretion on the entire envelope (the complete turbine with the rotating blades) is indispensable to account for the influence of the affecting parameters on the ice formation on the blades [2]. Differences in the predicted ice shapes via 2D and 3D simulations are demonstrated in a study on complex wind turbine icing events using FENSAP-ICE Simulation [3]. It has been demonstrated that higher velocity magnitude and higher ice growth were observed in 2D simulations than in 3D simulations; this behaviour was mainly associated with the absence of the flow interaction in the z-direction [4] Very few studies in literature performed fully 3D icing simulation [3][4][5][6]. These attempts are almost conducted for validation purposes. Testingdifferent possible scenarios of wind turbine operation under different icing conditions is time and computing cost-intensive. However, it is essential to mention that in all these studies, the simulation was carried out only for one blade, not for the entire turbine. The analysis was based on periodic conditions and assuming that all the turbine blades behave the same way regardless of their position. Due to the complexity of the physics of a highly timevarying phenomenon, it would be challenging to perform these kinds of simulations within the uncertainty over the mechanism of ice growth [7]. In view of that, it is highly recommended to own a reliable and powerful multi-physics simulation tool adapted to evaluate effectively icing impact on wind turbine performance and operation in cold climates. With the presence of that powerful tool, it could be possible to acquire more understanding of the conditions that lead to ice formation on the wind turbine and examine the effectiveness of the IPS used to reduce the impact of icing on wind turbines.
Recent advances in the 3D icing simulation codes developed for aeronautics can capture various 3D effects in aircraft that can also be advantageous for the simulation of the 3D effects of icing on a rotating wind turbine [2]. This paper focuses on the recent developmentsinmodelling and simulation of iced-up wind turbines'performance using advanced commercial icing tools. It provides guidelines to the CFD-BEM approach for generating the power curve for clean and iced-up wind turbines. The objective is to present the methodology of the Quasi-3D simulation to quantifying production losses under icing conditions. The present paper is structured as follows: Section 2 summarizes the CFD approach employed in the simulation of ice accretion. Section 3 presents the Blade Element Momentum (BEM) theory for the aerodynamic analysis of wind turbines. Section 4 presents the details to consider during the validation of the results. Section 5 summarizes, in a table, the academic research based on CFD-BEM mixed approach. Finally, The last section presents the main conclusion about using this approach.

Computational Fluid Dynamics for Ice Accretion Modeling
Computational Fluid Dynamics (CFD) is a commonly applied tool of fluid mechanics used to analyze fluid flows in order to generate solutions for complex physical phenomena using numerical methods. It has become today an indispensable tool for the research on wind turbine icing and de-icing. Researchers' actual interest is using Commercial CFD programs (such as ANSYS, OpenFOAM, and others) to simulate the behaviour of ice around wind turbine blades. The general objective is to evaluate the short-and long-term impact of icing on wind turbine performance and optimize the operation of the adopted icing protection techniques. Application of the CFD on aerodynamic analysis for wind turbines, including meshing and pre-processing, has been reviewed and discussed in several published works [8][9][10][11]. In this section, we will discuss the application of the CFD on ice accretion prediction for wind turbines in addition to few additional points on meshing and CFD model configuration, given that published papers on wind turbine icing CFD simulation do not provide much information about the configuration of the CFD model for the pre-processing phase of the simulation.

Calculation Domain
For CFD studies, a large computational domain should be generated along with an appropriate spatial discretization to consider all phenomena and avoid the exterior walls' negative effect (or mesh). It is also convenient to use the model in purely turbulent mode, leaving the solver to highlight the other laminar and turbulent models as well as the transition models proposed by the software [12]. Finally, regarding the accuracy of the model, it is crucial to define its boundary conditions correctly.
The computational domain for wind turbine icing simulation, being an external flow problem, describes the region around the geometry of interest (airfoil, blade or whole wind turbine) where the flow solution is required. Its shape and size depend mainly on the aerodynamic characteristics of the geometry and its effects on the surrounding flow field. For the whole wind turbine icing simulation, we should consider, for instance, the large magnitude of wakes of low energy flow left behind the turbine. In the case of an airfoil, we should similarly consider the separation of streamlines due to geometrical deformation because of icing at relatively high angles of attack. Therefore, the sides of the computational domain should be set far enough away from the geometry to ensure that the boundary conditions assigned to the outer domain do not influence the flow next to it, reducing the quality of the results. In fact, without forgetting the importance of low-cost computation, the further the borders are from the airfoil, the more space is left for the turbulence of the wake to fade out before reaching the boundary conditions imposed on the borders. We note in particular that the greater the incidence, the more turbulence there is and the more it takes a significant homothety ratio to achieve stabilization [13].
C-grid topology is typically used for simulating airfoils for both clean and iced conditions since it is more appropriate to account for turbulence due to high incidence that allows the angle of attack to be varied [14]. For icing simulation around wind turbine blades airfoil, a C-grid domain that extends 20 chord lengths in front of the airfoil and 30 chord lengths behind the airfoil would be appropriate to ensure that the entire effects of flow disturbance are captured in the fluid domain [15]. In Figure 1, the recommended dimensions, airfoil location and boundary conditions are shown. The proposed domain dimensions are based on previous explored icing simulation works that showed good agreement with experimental and numerical data [15][16][17]  636

Spatial Discretization
The spatial discretization or meshing process is an important step in numerical solutions in general. There are different mesh types (structured, unstructured, cartesian), and the correct choice will depend on the problem we want to solve. When simulating an airfoil under icing conditions, a structured mesh has been widely used due to its versatility to be optimized and regenerated when the ice starts to build up and the geometry changes. Figure  2illustrates an example of a structured mesh developed for icing simulation around NACA 0012 airfoil to be used with FENSAP-ICE. However, regardless of the implemented mesh type, an essential factor to consider for the accuracy of the simulations is sufficient refinement at the walls. Some examples are presented in Figure 3 and Figure . To take into an account the boundary layer phenomena (viscosity and turbulence effect), in general, the condition + less than 1 must be satisfied for a correct performance of the turbulence models (such as k -ω and k -ω SST) and to capture the velocity profile along with the airfoil boundary layer [8]. The evaluation of + depends on air density, freestream velocity, dynamic viscosity and the boundary layer length. It is noteworthy that air properties are defined for the atmospheric temperature considered in the simulation.

Mesh Adaption for Ice Contaminated Geometries
When Ice build-up, the airfoil geometry deforms, yielding some elements of the computational mesh around the airfoil boundary distorted and inappropriate. Remeshing, or morphing, therefore, becomes necessary using the iced 637 airfoil geometry obtained by ice accretion simulation to account for the new airfoil boundary and avoid convergence problems in the next run of the simulation or when computing the aerodynamic solution for the iced airfoil. In modelling ice accretion, the engineer must develop his mesh to be adaptive to geometry deformation as the shape of the ice is evolving. After each considered time-scale cycle of ice accretion, the grid has to be updated either by morphing or by re-meshing to conform to the new ice shape. The adaptive mesh will be able to hug the new geometry without distortion may producing degenerated elements. (See Figure ). For moving boundaries, such as mesh displacement due to ice accreting, FENSAP-ICE uses the ALE (Arbitrary Lagrangian-Eulerian) formulation [18]. Mesh adaptation techniques are provided with ANSYS FENSAP-ICE documentation and other related studies [19]. A mesh morphing-based technique for adapting and optimizing the mesh to efficiently handle ice accretion simulations was proposed by Costa, et al. [20]  Based on the ice growth shape calculated, the airfoil geometry deforms and is updated. Prior to conduct the final aerodynamic simulation (if the study is intended to estimate the aerodynamic performance of the iced airfoil), the mesh should also be updated to hug the newly generated iced boundary. Once the overall ice-accreted shape is evaluated for the desired period of time, we can look at performance analysis, in which we do final aerodynamic calculations to estimate the performance degradation due to the deformation of the blade's airfoil. Depending on the objective of the study regarding the desired ice shape surface accuracy, the grid for the overall geometry may be generated after a single run for the entire duration of ice accretion, a simplification leads to inaccuracy on the boundary generated, or alternatively, multiple runs can be performed by dividing the total icing time into shorter periods in which, the grid is updated every time to account for the periodic ice growth before conducting the next cycle of simulation [21]. The thickness of the ice is determined, and a new surface grid has to be generated to be adapted with new geometry. This has to be performed in a repetitive manner to have enough accuracy on the iced boundary reconstructed. The overall calculations estimate the total mass and form of ice over the entire duration of ice accretion. Performing a single simulation cycle for the entire period of ice accretion provides simple and approximated geometries, while a process of multiple simulations is more computationally intensive with a more accurate representation of the shape of ice resulted.
To demonstrate the influence of the choice of each of the last-mentioned schemes of simulation, we conducted two icing simulations on the airfoil NACA 64-618 of a chord length of 1.419 m situated at 95% of the blade span of the NREL 5MW, pitch-controlled virtual wind turbine. The two schemes of simulation have been used for the same icing conditions: MVD=20µm, LWC = 0.22g/m 3 , relative velocity V r =75,88 m/s, Total accretion time t=60 minutes, AOA=5.824° and atmospheric temperature T=−10℃. The k -ω SST turbulence model is used, and the Shin et al. roughness model is considered for the single-shot approach and the Beading model for the multi-shot approach. The results showed that the multi-shot scheme of simulation provided a more detailed boundary of the formed ice shape, which reflects mainly on the surface roughness effects on the subsequent simulations as well as on the aerodynamic parameter estimation (see Figure 4). The process of mesh and updating the governing parameters through multiple simulations can be performed individually by transferring the solution files for every periodic simulation or, if enabled by the solver, it can be done in automated sequences of multiple runs for the aerodynamic, trajectory and thermodynamic modules. This is to be enabled to repeat the other calculations (flow, impingement and ice growth) without having to re-mesh for every considered modular iteration. Recently published papers made use of this functionality in the advanced icing CFD software [22,23]. STAR-CCM+ software provides a 3D morphing motion technique to predict the time-scale progression of ice formation to account for the deformation of airfoil geometry by adjusting the portions of the mesh near the affected wall boundaries [23]. FENSAP-ICE provides this functionality through the -Multi-shot Ice Accretion with Automatic Mesh Displacement‖ that enables automated feedback loops with the updated parameters for the configuration and the execution into a single setup. For example, if we have 21 minutes of icing event, a single shot will use 21 minutes to estimate the ice thickness based on the initial conditions of geometry and governing parameters. Whereas, if multi-shot is adopted with seven discrete time-steps, the governing parameters and the geometry will be updated every three minutes of simulation [18]. In Figure 5, a scheme for both approaches is presented. For every shot of ice calculation, it is noteworthy that it is recommended to integrate roughness height when running flow solution (for example, in FENSAP-ICE, the -Beading model‖ is used to satisfy this condition) [24]. This option improves the computed ice shapes considering roughness effects on the shear stress and heat flux as well as the aerodynamic parameters estimated for the final resulting ice shape [18,24]. 639

Mesh convergence study
Accuracy is crucial for icing simulation to minimize dependence on costly experimental tests. The mesh generation is a critical phase in a CFD-based icing analysis, given its influence on the calculated solution. An excellent quality mesh is essential for obtaining an accurate, robust and meaningful calculation result. Meshing defects have a severe impact on the convergence, the precision of the solution and especially on the computing time. A good quality mesh is associated with the minimization of the elements presenting distortions (skewness). It is essential to achieve zero or marginal skewness near the airfoil surface [25]. Otherwise known, to ensure the continuity of any numerical solution, the accuracy is dependent mainly on the mesh used, not only in terms of its refinement and element size but also in terms of elements' orientation and aspect ratio and on a good refinement in the regions presenting a strong gradient, such as the boundary layer. A good mesh must also be sufficiently smooth and undergone a mesh refinement study. In order to minimize potential error and optimize calculation time, a mesh convergence study (also known as mesh independency/refinement study) has to be conducted to determine the optimum number of mesh elements for the aerodynamic performance simulation [23]. A Mesh refinement study involves conducting calculations with a successive mesh refinement until obtaining valid converged results for the most relevant variables. For each mesh refinement, the size of the elements will be reduced, which means that each element will be replaced by smaller elements. The results for each variable of interest are observed at each step. The convergence of an iterative process is achieved when a precision criterion is satisfied. Since the RANS equations are nonlinear, the numerical solution will not converge towards an exact solution whatever the mesh is refined. In fact, the exact solution is unknown; the error for each iteration is then based on the change of the calculated quantity between two iterations. A solution error should be verified and evaluated depending on a specified level of accuracy.
For icing simulation, our most interest reflects on the aerodynamic parameters' estimation. Therefore, a convergence study is essential to be applied on lifts (Cl) and drag (Cd) coefficients for a representative angle of attack, usually that of the operational design of the chosen wind turbine (see Figure 6). Similarly, for icing simulation convenience, it is advised to conduct this convergence study for both Cl and Cd of the iced airfoils. As an example from the wind turbine icing literature, a test of convergence on the lift coefficient has been conducted by Zanon, et al. [17]to model ice accretion on the airfoil NACA64-A17 on the NREL 5 MW reference wind turbine. Another example from wind turbine icing literature, a discretization error between grids has been tracked by Hildebrandt [15] in a mesh refinement study carried out on Cl, Cd for both clean and iced airfoils while introducing the time of simulation as an additional decisive factor for the choice of mesh. In commercial software, different approaches are used to achieve grid refinement convergence. Basically, we may opt for the adaptive meshing tool provided by the software to automate the mesh refinement until a user-specified level of accuracy is reached.
One can seek convergence study globally all over the domain or locally wherever the refinement is essential. For icing simulation on airfoils, sufficient mesh resolution should be applied near the stagnation region to capture the local characteristics of water droplets impingement accurately (see Figure 3 and Figure .). Therefore, a convergence study may be applied locally to account for airfoil deformation due to ice accretion in specific regions (around the airfoil or near the stagnation region within the impingement limits). To judge the required density implemented in the regions of interest, only a mesh convergence study can confirm that the mesh has sufficient refinement in this region. Generally speaking, smooth transitions between regions of different densities have to be respected. Rapid changes in mesh density or elements of deficient shape cause artificial disturbances and poor results. We should also be careful that the shape of the elements in these regions of concern should not stray too far from the reference shape when the reconstructed mesh has to hug the ice accreted. 640

Boundary conditions
For the icing problem, each module (aerodynamic, droplet trajectory and thermodynamic) deals with a system of partial differential equations, for which the solution requires the determination of appropriate boundary conditions to be specified to the discretized domain's boundaries. Some of the most commonly used boundary conditions are velocity-inlet, pressure-outlet, anti-slip wall, among others. However, in the icing problem, some of the boundaries will require more than one condition. For example, considering the icing problem as a multiphase flow (air as continuous phase and supercooled water droplets dispersed in the air), if we chose the air as a continuous phase, this will be automatically transferred to the droplets. This means that the droplets will be forced to go around the airfoil and avoid impingement, which leads to errors in predicting the water droplets' impingement and, consequently, the ice shape. Therefore, the airfoil should behave as a permeable wall for droplets but simultaneously as a wall for the air. In order to solve this problematic issue, some CFD software allows to couple a User Defined Function (UDF) with the calculations to account for the continuity, momentum, and energy equations for the droplets [16].

Blade Element Momentum (BEM) Theory
The Blade Element Momentum (BEM) theory is commonly used in wind turbine design to calculate the aerodynamic performance and the power curve of turbines [8,26]. These calculations can be performed for both clean and iced wind turbines in order to estimate the power losses due to icing [26].  From the momentum theory, we can obtain equations to define the thrust and torque on an annular section of the rotor as a function of the axial ( ) and tangential ( ′ ) induction factors of the flow conditions [27].

Equation 2
Where, is the air velocity and the rotational speed of the rotor.
From the blade element theory (see Figure 8), the lift and drag forces are perpendicular and parallel to a relative wind speed ( ), respectively, which is the vector sum of the wind velocity at the rotor ( 1 − ), and the wind 641 velocity due to rotation of the blade ( 1 + ′ ) [27,28]. Figure 8 shows the relationships between the angles, forces and velocities in one section of a blade. From that, it can be seen that:

Equation 5
Where, is the axial induction factor, = is the local tip speed ratio, is the angle of attack, is the angle relative to the wind, and is the pitch angle.
Regarding the above, the blade element theory also allows obtaining two more equations, which define the normal force or thrust ( ) and the tangential force or torque ( ) on the annular rotor section as a function of the flow angles at the blades and airfoil characteristics.  According to BEM theory, the blade element force causes the change of momentum of the air passes through the disk swept. As a result, it is assumed that there is no radial interaction between the flows passing contiguous disks; however, this condition is only valid if the axial flow induction factor does not vary radially, and it is uniform, i.e. the blades have uniform circulation. Ref. Section 3.5.2 of Burton, et al. [28].
Using two-dimensional aerofoil characteristics, the axial and tangential flow induction factors and ′ can be estimated iteratively [3]. The following equations are useful for carrying out the process [27,28]:

Equation 9
Where, = cos + sin , = sin − cos , is the lift coefficient, is the drag coefficient and is the chord solidity which is defined as the total blade chord length ( ) at a given radius ( ) divided by the circumferential length at that radius: To determine the complete performance characteristic of a rotor, i.e., how the power coefficient varies over a wide range of tip speed ratio, requires an iterative solution that consists of initializing and ′ to be zero, assuming = 0. Drag is excluded from the calculation of the flow induction factors because drag in attached flows is caused by skin friction and does not affect the pressure drop across the rotor [28]. Once the induction factors have been obtained, , and are calculated; then the process is repeated until convergence is obtained [27,28]. With the final results for , ′, , and , the torque can be calculated by equating equations 4 and 9 [28]:

Equation 11
Finally, the power contribution from each annulus is given by [27]:

= Equation 12
And the power coefficient is given by [27]:

Equation 13
Power coefficient Cp of the wind turbine depends, therefore, in addition to the number of blades , on the aerodynamic characteristics of the airfoils of the considered sections in the blade, their radius , their pitch angle , the angle of incidence , the upstream wind speed , and the rotational speed . The previous missioned parameters of the velocities triangle of airfoils can be reduced to a single variable, the relative velocity , that is used as an input parameter for the CFD calculations.

Improved BEMT:
Due to the assumptions mentioned at the beginning, this method has deficiencies in the study of three-dimensional flows. Therefore, some corrections such as Prandtl's tip loss factor, Spera's correction, and the Du-Selig stall delay model have been defined in order to improve BEM accuracy of prediction [8,29,30].

Prandtl's Tip Loss Factor
Air tends to flow around the tip from the lower to the upper surface, reducing lift and thus power near the tip. This is due to the suction side of a blade has lower pressure than the pressure side. This is more evident in turbines with fewer and wider blades. The correction factor (see Equation 16) was introduced into the equations discussed previously, mainly affecting the forces derived from momentum theory.

Equation 14
As can be seen in Eq. 14, this correction factor is a function of the number of blades ( ), the angle of relative wind ( ), and the position on the blade ( / ). Special attention should be paid to the angle in the inverse cosine, which should be in radians [29].

Stall Delay Model
The stall region in rotating blades begins in areas with a higher angle of attack than a static airfoil in a wind tunnel. The effects of centrifugal force help to reduce the adverse pressure gradient that causes an airfoil to stall. Aerodynamic stall can be mathematically described using laminar boundary layer theory as follows [8]:

Equation 16
Where, ∆ and ∆ are the differences in lift and drag coefficients obtained from the flow with no separation, / is the local chord length normalized according to the radial direction. Considering the effects of tip loss and stall delay, the induction factors can be written as [8]:

Viterna -Corrigan (VC) Stall Model
This model is based on the fact that HAWTs have a high aspect ratio with a characteristic flow field similar to a 2D airfoil. Thus, the lifting line theory has been used to calculate the lift coefficient before the stall angle of attack. When the flow moves into the fully developed stall region, the VC model is used to predict the lift and drag 643 coefficients for an angle of attack ( ) higher than the angle of attack in which the maximum coefficient is observed ( ), and is written as follows [8]:

Equation 20
Where, and the are:

Equation 22
Where, , and C d,s are respectively the lift and drag coefficients at an angle of attack α s of 20⁰ in which , depends on the aspect ratio, as follows [8]: Spera's Correction Spera identified an empirical relationship between the angle of attack and the thrust coefficient ( ) in the turbulent wake state, which is written as [8]:

Equation 25
Whereas if < 0.96, then the axial induction factor is solved using Equation 19.

CFD-BEM Mixed Approach (Quasi-3D simulation method)
In a 2D simulation, ice accretion is reproduced over only one airfoil section along the blade. It provides important information about the aerodynamic performance drop while the airfoil shape is deformed by ice accretion. However, this kind of study is not representative of quantifying the power production losses due to icing. Even though some references consider that the airfoil section at 85% spanwise position represents the average ice accretion on a blade power generation [21,31,32], and at that radius, arises the most significant contribution to turbines' power production [33]. Nevertheless, 2D solutions do not consider the 3D effects of flow. Fully 3D solutions simulate the blade as a whole in a 3D domain and take the 3D flow effects into account [5], except that this type of simulation is tricky to implement, especially when it comes to several scenarios of icing conditions. The Quasi-3D simulation is a rational alternative method to incorporate the simplicity of 2D analysis and the capability of the improved BEM theory to capture the 3D effects of flow [34]. By coupling the BEM theory with the solution of RANS equations via CFD simulation can be used to generate the power curve for a wind turbine. The application of the CFD-BEM approach to generate the power curves for both clean and iced WT helps estimate wind turbine production losses due to icing resulting in a less computationally expensive solution. Both methods (the Quasi-3D and the Fully-3D simulation) have a similar theoretical background, but a different order of modelling but a different representation of icing effects on wind turbine blades.Sagol [6] carried out a critical analysis using both a Quasi-3D and a Fully-3D simulation approach for NREL Phase VI wind turbine blades. The study shows that the two approaches resulted in different ice shapes and eventually different iced blade performance. The fully-3D method has more accurate predictions for a clean wind turbine blade with difficulties validating iced-up wind turbines in the absence of experimental data. Even though Fully-3D simulations have higher computational cost, 2D icing analyses, by not considering the 3D features of rotational flow, may lack accuracy [6].
In the Quasi-3D simulation method, the computed aerodynamic coefficients for the clean and iced blades using the CFD are employed together with the BEM theory to generate the wind turbine's torque, power, and Cp curves for both normal and icing conditions [26]. Firstly, estimate ( ), ( ) for clean and iced airfoils for several sections of the blade using a CFD icing simulation software. Then estimate the power with any software that uses the Blade Element Momentum (BEM) theory in its code. The iterative method is explained by Switchenko, et al. [3]. The polar needs to be extrapolated to 360• for each section in order to determine load and power calculations under various operational conditions [23]. This can be done in the polar extrapolation module using the Viterna equation

644
[23, 35,36]. Only extrapolated polar data can be used to simulate a turbine or rotor. Estimating Cl, Cd using the CFD solution of the Reynolds-averaged Navier-Stokes equations RANS is more accurate.
The aerodynamic losses due to icing are represented by the reduction in Cl and the increase in Cd. They are usually demonstrated through drawing the Cl and Cd for both clean and iced wind turbines on the same graph. However, what determines the drop in aerodynamic performance is the reduction in the lift to drag ratio (represented by Cl/Cd) due to ice-induced roughness.The power losses are represented by the degradation of the wind turbine's power curve due to ice. They are usually demonstrated through drawing the power curve P(V) or the power coefficient C P (V) as a function of the wind speed or the power coefficient C P (TSR) or the torque coefficient T(TSR) as a function of the Tip Speed Ratio, for both clean and iced wind turbines on the same graph. Figure 9 illustrates the overall process of the Quasi-3D Simulation using CFD solver to compute the aerodynamic characteristics ( ) & ( ) for each section's airfoil for both clean and iced airfoils, and the application of the BEM Theory to generate the power curves for both ordinary and icing conditions using the CFD computed aerodynamic characteristics & ( ).
The annual power production losses resulting from ice build-up on wind turbine blades can be estimated at a particular location using wind speed frequency data for that location as well as a power curve generated for both clean and iced-up wind turbines [37,38]. As an integral example for the preceding steps of the proposed methodology, Turkia, et al. [37] describe a method to develop the reduced power curves for a typical variable speed pitch controlled 3 MW iced-up wind turbine. The icing conditions were defined using a numerical weather prediction model. A relation between the standard stationary cylinder and the rotating wind turbine blade was determined using an ice accretion model described in ISO 12494 (2001) standard [39]. Ice accretion is modelled on two blade sections in three specific icing events using the VTT in-house code TURBICE, while the airflow was modelled using ANSYS FLUENT. The power curves were generated for both clean and iced conditions using FAST software. The results were validated against real-life observations. Limiting the study to the dry type of ice has reflected in underestimating the icing impact on power production.
When thinking in terms of the reliability of numerical approaches, the use of the CFD-BEM method is questionable from the point of view of Bai and Wang [8]. The BEM theory is not reliable for simulating the distribution of the aerodynamic loads on wind turbine blades, particularly for cases related to the separation of the boundary layer on the suction surface of the airfoil [8], a phenomenon considerably expected for the ice contaminated airfoils. Dynamic stall on the airfoils dramatically increases yaw moments and the cyclic blade loads [40]. The CFD-BEM approach cannot predict the detailed flow field produced by a finite number of blades. It is effective in predicting thrust but tends to overestimate the power. Therefore, for generating the power curve for an iced-up wind turbine, the CFD-BEM approach needs to improve its accuracy by considering the 3D effects around the blade [8] 645 Figure 9:-Power Loss Calculations using CFD-BEM mixed approach @Fahed Martini.

Validation options
Since numerical tools represent an increasingly used option in wind turbine icing research [35], there is necessary to compare the experimental and numerical icing simulation results for the typical wind turbinesunder icing conditions [41]. Whatever the objective of conducting ice modelling studies for wind turbines, CFD numerical solution must be validated before using its results. While almost all the available commercial engineering CFD codes that integrate icing modules were originally developed for aeronautics, it is important to carefully validate the code for the range of Reynolds numbers used for wind turbines. So far, numerical codes developed to predict ice accretion for aircraft have been recently validated more or less for wind turbine icing applications. However, there is still a need for validation against experimental results for typical wind turbines in different icing conditions and various operation scenarios. Generally speaking, one should have a conception of the type of results before conducting the simulation. Therefore, it could be convenient to start with a simple case for which we already know its results to have confidence in the owned simulation tool. In this case, several studies validated the CFD results with a 2D experimental tests [3,14,16,17,21,23,37,[42][43][44][45][46][47][48][49][50][51]. The validation process is important to demonstrate the accuracy of the model used in the simulations.
For a complex phenomenon such as icing, several upstream studies must be conducted before proposing a precise model capable of achieving a three-dimensional icing simulation on wind turbine rotating blades. Several CFDbased published studies used the NACA 0012 experimental results for validation. Some papers use the NACA64-618 since it is chosen for the outboard third of the blade for the virtual NREL offshore 5-MW baseline wind turbine. Besides, this type of airfoil (NACA series 6-Digit) has been used for large wind turbine blades [52].The first option consists of validating the resulted aerodynamic characteristics of the clean airfoil. The reference information on the aerodynamic performance of clean airfoils is available in literature either from theoretical studies or from experimental tests. Consequently, the validation of the clean airfoil is contained in the context of satisfying the searcher in the simulation tool that he uses. The iced airfoil validation consists of verifying the ice shape obtained under the same icing conditions.Two common wind turbines are used in literature for validation, mainly of the clean scenario: The two-bladed NREL-phase VI real horizontal-axis wind turbines: This turbine, composed of S809 airfoil, was experimentally characterized using wind tunnel tests at NASA's Ames Research Center in normal conditions. The characteristics information required to quantify the full-scale, threedimensional (3-D) aerodynamic behaviour of this turbine is available in the technical report provided by the National Renewable Energy Laboratory (NREL) [53]. This turbine often works with a constant rotor speed at various wind speeds, which differs from the variable-pitch and variable-speed wind turbine in operation mode nowadays [1]. Han, et al. [54] present ice accretion tests on this model.
The conventional three-bladed upwind variable-speed variable blade-pitch-to-feather-controlled virtual wind turbine known as the -NREL offshore 5-MW baseline wind turbine‖ The specifications of this representative wind turbine are available in the technical report by the National Renewable Energy Laboratory (NREL) [55]. The NREL 5-MW wind turbine blade comprises six different airfoils (Span=63m), see Table 1.

Summary of Research Survey on the Quasi-3D Simulation Studies
Results of the investigation on ice accretion modelling and simulation for HAWTs are presented in this section. The objective is to summarise the up-to-date attempts of quantifying the icing impact on wind turbines energy production via the 3D and Quasi-3D simulation. The objective of each of these studies is summarised along with other valuable information, such as the type of ice treated, the software used for both CFD and BEM theory simulation, the geometry and the wind turbine type considered, the parameter used to represent the results of the simulations, in addition to the availability of other types information such as the adopted roughness and turbulence models and the presence of experimental validation. Mainly, the two above-mentioned types of wind turbines (the NREL 5 MW and the NREL Phase VI) were used for validation, while other types were treated too.
From reviewing the simulation studies in the literature on wind turbine icing, the studies that conducted 3D and Quasi-3D simulations based on the CFD-BEM approach are listed and summarised in Table 2. The review is limited to HAWT investigations.

Conclusion:-
Several studies have been conducted on the icing simulation for wind turbines. Most of them were 2Dsimulations based on computational fluid dynamics (CFD). Despite few attempts to simulate one blade in a 3D symmetric configuration, no full-scale simulation of a rotating wind turbine under icing conditions was considered. The CFD-BEM approach is a practical alternative to generate the power curves for both clean and iced-up wind turbines in order to estimate the power drop due to icing. In this paper, a comprehensive review of the use of the CFD-BEM method in wind turbine icing research is presented in the context of estimating the power curve degradation due to ice accretion. Quasi-3D simulation studies based on the CFD-BEM mixed approach are reported relying on a methodical and concise literature survey. The literature review revealed that more than few studies used this approach to estimate the energy loss during specific icing conditions, as summarized in Table 2. However, these studies adopted different conditions for multiple types of wind turbines using different software programs. Furthermore, programs that use the BEM theory are subject to limitations in the calculations of power curves. Therefore, additional losses related to wind turbine geometry should be considered.Using the Panel method to estimate the aerodynamic characteristics (as in the case of XFOIL in QBlade) is not recommended for the iced airfoils, especially for high angles of attack. The deformation of the airfoil, mainly on the leading edge, will increase the possibility of flow separation on the airfoil, which could not be captured by the panel method.
Finally, in order to better advance the research in this field, it is recommended to unify the research on icing simulation by adopting one of the above-mentioned types of wind turbines. The NERL-Phase VI is composed of one 649 airfoil type and has the characteristics available to quantify the full-scale, three-dimensional (3D) aerodynamic behaviour. In addition, experimental results are available for validation for this wind turbine.

Data Availability Statement:
This study does not report any data.

Funding:
This research received no external funding

Conflicts of Interest:
The authors declare no conflict of interest