Numerical Investigation of Aerodynamic Performance of Wind Turbine Airfoils with Ice Accretion

One of the emerging problems in modern computational fluid dynamics is the simulation of flow over rough surfaces, wind turbine blades with ice on its surface for instance. An alternative method to numerically simulate rough surfaces instead of using a grid with y+ < 1 criterion is to use rough wall functions (RWFs) that models the flow behavior in case of the presence of roughness. This work aims to investigate different rough wall function models to find out the model that can provide the most accurate results with the lowest computational resources possible. This aim was achieved by comparing coefficients of lift 5 and pressure resulting from CFD simulations with wind tunnel results of an airfoil with actual ice profiles collected from the site. After implementing new RWFs in OpenFOAM framework and validating the results with published experimental results, the comparison showed that momentum RWF provided the best agreement between simulation and experimental results while using only 25% of the number of cells used with smooth RWF. The conclusions of this work will be applied in a simulation code within OpenFOAM framework to simulate airflow fields of wind turbines with ice accretion. 10


Introduction
For the last few decades, energy generation from wind has increased rapidly as a renewable source of energy. Because of this rapidly increasing demand for wind energy, wind turbine manufacturers have been trying to increase the size of the turbine rotors. By increasing rotor diameter and tower height, the single turbine can now harvest more power from wind in almost the 15 same horizontal area occupied by old, small size wind turbines. However, some of the best wind energy sites like in Europe and North America suffer from icing atmospheric conditions. These conditions can decrease the annual energy production by up to 40% as shown by Sailor et al. (2008). This drop in blade aerodynamic performance occurs due to several detachments and reattachment areas on the surface due to the presence of a rough ice surface.
Since the 1930s, different experiments were conducted in order to investigate the effect of roughness on fluid flow in general. 20 One of the pioneers in this filed was Nikuradse (1933) who investigated turbulent flow in pipes with different relative roughness values and Reynolds numbers between 10 4 to 10 6 . He noticed that the boundary layer follows the log-law as in the case of smooth surface. However, in case of rough surface, there is a clear velocity shift (∆u). A few years later, Schlichting (1937) studied the internal flow of a square-section channel with one rough wall. This rough wall had spherical segments, cones and angular shapes. He found that the velocity shift depends on four different geometrical parameters of roughness elements, 25 namely: element's height, the area projected on the surface, area projected in flow direction and the average distance between elements. Accordingly, he derived an equation to calculate the friction coefficient on the plate surface. Later, Moody (1944) provided an engineering method to calculate friction losses in pipes due to roughness.
For external flow applications, Blanchard (1977) and Hosni et al. (1991) provided measurements of heat transfer performance of a heated flat plate with hemispherical roughness elements. Also, Hosni et al. (1993) studied the effect of roughness shape 30 on the boundary layer. Achenbach (1977) studied heat transfer of a roughened circular cylinder at different roughness heights and Reynolds numbers.
From computational point of view, it is known that to correctly model the viscous sub-layer in computational fluid dynamics (CFD) methods, the height of the grid's first cell should fulfill the y + ≤ 1 criteria. Additionally, to avoid high aspect ratio cells, large number of cells over the surface should be generated. Accordingly, one method to simulate rough surfaces is to generate 35 a fine grid that fits the rough surface. Wang and Zhu (2018) simulated ice accretion on the NREL Phase VI wind turbine blade.
In their numerical setup, the used a grid y + ≈ 1, which resulted in a grid with total of 7.3 million cells for a half-cylindrical domain for one blade. This gives a good idea about how much computational cost is needed to perform such simulations in the standard way. However, this cost can be much more increased for larger turbine blades manufactured nowadays.
Another approach is to modify the mathematical models in a way that the behaviour of flow over rough surfaces is grasped 40 without the need for a fine grid. The growing use of CFD in fluid flow simulations in both research and industrial applications has caused an increased interest in this approach and more new models have been developed. The challenge is to find the proper mathematical models which are often based on the results of the aforementioned experiments. Chen and Patel (1988) introduced a two-layer model and a wall function in the k-model to simulate rough surfaces. Hellsten and Laine (1997) provided an extension for the k-ω-SST turbulence model to simulate the behaviour of flow, so did Wilcox et al. (2006) with 45 the k-ω model. Instead of providing changes to turbulence models, rough wall functions (RWFs) were developed to simulate the behaviour of turbulent flow near walls. In general, wall functions must provide adequate handling for flow change from its stream velocity to stagnation on the wall according to the no-slip condition. Thus RWFs should provide additional models that account for roughness. Recently, Da Silva et al. (2011) applied a new ν t wall function base on the work of Suga et al. (2006). Also, Knopp et al. (2009) and Chedevergne and Aupoix (2017) provided k and ω wall functions that are capable to adapt for 50 the presence of roughness on the surface.
In our work here, we aim to find a computational model that is able to simulate the airflow around an iced wind turbine airfoil.
This model should provide the most accurate results using the least computational resources possible. This aim is approached by comparing the CFD simulation results of different RWFs with experimental results of actual wind turbine ice accretion profiles. The studied ice profiles were collected from the site and tested in a wind tunnel after being molded and attached to the 55 airfoil. In this work, three different RWFs were implemented in OpenFOAM ® v6 (Greenshields (2017)) and compared : Knopp et al. (2009) RWF , which will be referred in this work as DLR-RWF, Da Silva et al. (2011) RWF, which will be referred as Momentum-RWF, and Chedevergne and Aupoix (2017) RWF fitted to Colebrook's experimental results, which will be referred as Colebrook-RWF.
This article is structured as follows, Sec. 2 explains the mathematical model of each of the used RWFs. After that, Sec. 3 60 explains the preparation of the rough ice for the simulation and the corresponding computational mesh generation. Sec. 4 provides the results of each case studied in this work. The last section provides discussion and conclusions of this work and how these results can be applied to general rough surfaces simulations.

Mathematical Models
In this section, a brief description of the RANS models used in flow field simulations are introduced as a background for the studied RWFs. Subsequently, the details of the implemented RWFs are also explained given.
In RANS models, Navier-Stokes equations are solved as time-averaged while extra one equation is solved (like in Splart-Allmaras) or two equations (like in k-ω SST) or more equations are solved to close the system of equations required to solve the flow field (Versteeg and Malalasekera (2007)).

70
In this work, only steady-state RANS turbulence models, namely Spalart-Allmaras (Spalart and Allmaras (1992)) and Menter kω SST (Menter (1994)), were used to simulate the turbulent flow field over the iced airfoils. In general aerodynamic simulation problems, both models show good results, especially in the pre-stall AoA's. However, each model has its own advantages. Since Spalart-Allmaras (SA) turbulence model is a one-equation model, it is usually more stable and requires less computational power. On the hand, k-ω SST can reach better results due to solving more turbulence parameters, but of course with more 75 computational cost.
In this work, these two different turbulence models were specifically chosen to fit for different wall functions. As will be explained in the next section, one of the used RWFs defines a relationship for turbulent kinematic viscosity ν t . This type will be used only with SA model since it mainly solves ν t equation to simulate the flow field. The other two RWFs adapts k and ω values near the wall for the presence of roughness. Accordingly, these two RWFs will be used only with k-ω SST turbulence 80 model. Nikuradse (1933) found that the boundary layer log-law behavior in flow over both rough and smooth surfaces are the same except for a velocity shift (∆u). The difference between most of the proposed RWFs in various literature is how to calculate this shift. Schlichting (1937) worked out the parameters on which the velocity shift depends on, the four parameters: average 85 elements height (K avg ), the area projected on the surface (A p ), area projected in flow direction (A s ), and the average distance between elements (L avg ) (as shown in Fig. 1) are contracted together to one dimensionless parameter called equivalent sand roughness height (K s ) that can be calculated by the equations presented in Dirling (1973)

Rough wall functions
(2) Figure 1. Roughness element geometrical parameters So, once the roughness profile is known, K avg , L avg , and D avg can be calculated. Assuming roughness elements take certain geometrical shape for approximation, A p and A s can be calculated.In the current study, roughness elemnts are assumed to take conical shapes. Accordingly, A p /A s = πD avg /2K avg and Equations 1 and 2 can be implemented to calculate K s . The final value of the velocity shift in all RWFs used in this work depends on the calculated value of K s . that represents this velocity shift as a function of K s according to Cebeci (2004) ideas in the equations Finally, turbulent viscosity can be calculated by: ). 105 Since this wall RWF deals mainly with ν t , it is used with Spalart-Allmaras turbulence model to simulate rough surface.

DLR RWF
In their work, Knopp et al. (2009) followed the ideas of Aupoix and Spalart (2003) in adding the velocity shift effect to the wall function. However, instead of applying the shift on turbulent viscosity as in Aupoix and Spalart (2003), they used their procedure to adapt k and ω turbulence parameters. The new k and ω RWF states that Also the same like DLR RWF, this one is used with k-ω SST turbulence model.

Profiles and grids of CFD simulations
To enable the usage of the aforementioned RWF's, rough ice profiles have to be smoothened to give an equivalent smooth 125 surface. By using this smooth surface and the K s value calculated from Equ. 1, the RWF's can be used in the simulation cases.
This section explains how these profiles should be prepared and how the computational grid should be generated.

Wind turbine airfoil with ice accretion
To investigate the different wall functions described in Sec. 2, they were applied on a wind turbine airfoil with two different ice accretion profiles. The two ice profiles were collected from site, scanned and then machined as a leading edge attachment 130 for wind tunnel measurements.
Both ice profiles take the horn-ice form Fig. 2 which usually occur under severe icing conditions and also both of them extend to 7.5-10 % of chord length. However, it can be noticed that profile 1 is smoother than profile 2. Also it can be noticed that profile 1 takes a more aerodynamic shape than profile 2 since profile 2 has two horns which forms a stagnation area between them.

Roughness parameters calculations
Since all RWFs treat the rough wall as smooth wall with a velocity shift as explained in Sec. 2.2, the actual rough surface of the ice profile is replaced with another equivalent smooth surface. The new smooth surface will be used to generate the computational grid around the profile and will be numerically treated as rough surface, i.e. a velocity shift will be added to the smoothened surface. To calculate this new surface, the rough surface was smoothened with cubic splines and . Similar to the 140 calculation procedures of roughness parameters indicated in the DIN-EN-ISO 4287 standard (DIN (2010)), the rough surface should be smoothed as shown in Fig. 3a, using a cubic spline in this work, to find an average, smooth surface. After that, this smooth surface could be mapped to the x-axis as shown in Fig. 3b to ba able to analyze the roughness parameters.
Knowing the distance between roughness elements and height of elements, average sand roughness height explained in Fig.1 can be calculated using Eq. (1). By analyzing the laser scanned ice surfaces in different literature like Hawkins et al. (2017) 145 and McClain et al. (2018), it can be assumed that roughness elements take conical shapes. Accordingly, A p = πD 2 avg /4 and A s = 0.5K avg D avg . The above analysis gives results in K s = 1 × 10 −3 m and 1 × 10 −2 m for profiles 1 and 2 respectively.

Grid generation
In order to use rough wall functions, the height of the first cell center should be large enough to cover the roughness element.
Along with converting the rough surface into a smoother one, the resulting grid is much coarser than the grid required by 150 smooth wall functions which require y + (1) < 1 to be able to correctly simulate the boundary layer. Accordingly, the studied approach in this work requires less computational resources. In case of the two ice profiles studied in this work, to properly generate a grid that fulfills the condition of y + (1) < 1, is found to require number of cells ≈ 4 × 10 5 cells in 2D simulations.
On the other hand, when using rough wall functions with the proper first cell height and roughness smoothing, it is found to  To provide a comparison between the results of fully resolved surface roughness and using RWFs, the next section will also provide a comparison between C l values resulted from using fine and coarse grids, shown in figures 5a and 6a respectively, and coarse grids, shown in figures 5b and 6b, used with RWFs.

160
It should be noticed that the first cell height in case of RWF simulations is relatively larger that in case of fully resolved mesh to fit the criteria explained by Liu (2014). This criteria mandates that the height of the cell center of the first cell near the wall should be larger than the roughness height. Accordingly, a mesh independence test was done to select the first cell height values indicated in Tab 1.  For simulation results, validation case will be introduced to give evidence of correct implementation of the new codes. Then, the results of the simulation cases shown in Tab.1 will be shown.

Implementation and validation
To apply the aforementioned RWF, equations 7-13 were implemented within OpenFOAM framework. The Momentum RWF is 170 already available in OpenFOAM (named nutURoughWallFunction). However, the other two wall functions, namely DLR and Colebrook RWF, are not available. These two wall functions were implemented and validated against experimental results published by Achenbach (1977). In this paper, the author analyzed the flow field around a circular copper cylinder with pyramidal roughness elements to study both fluid flow and heat transfer. This experiment was done on a 0.5 m length, 0.15 m diameter copper cylinder at different roughness heights and Reynolds numbers.

175
It is known that RANS simulations of flow filed around a cylinder is already unstable and hard to predict, especially when it comes to the prediction of separation location. However, this case was selected to identify any false implementation of RWFs in OpenFOAM that may lead to unstable numerical solution. While analyzing the results of these validation cases, one must keep in mind that such a complicated simulation case of flow over a cylinder at high Reynolds number using steady-state RANS simulation is difficult. Accordingly, one should not expect good agreement between numerical and experimental simulation on all locations over the cylinder. These low expectations comes from the fact that all RANS models fail to accurately predict the flow separation location even for airfoils with relatively 195 high AoA. For cylinder case, the situation is even harder. Since the aim of this validation process is to make sure that the RWFs have been correctly implemented and are working properly and not to judge the mathematical behavior of these RWFs, it can be concluded that the RWFs follows the trend of C p distribution over the cylinder. 4.2 Case 1: profile 1 at Re = 2.6 × 10 6 : In this case, profile 1 was tested at Re = 2.6 × 10 6 , which is a relatively low Reynolds number compared to the other two cases.  This case is exactly like the previous case except for being tested at Re = 5.1 × 10 6 . This high Reynolds number is challenging for RWFs since it leads to more violent separation and hence are harder to be predicted.

215
This case shows good agreement (Fig.10)  In this case, the flow over profile 2 was simulated. It can be noticed from Fig.2 that profile 2 has not only two ice horns, but also it has a rougher surface , i.e. has a higher K s value. The effect of this complex shape can be obviously seen   Fig.13b, Colebrook and DLR RWFs shows large deviations over the whole airfoil while Momentum RWF shows better agreement. In figures 13c and 13d good agreement between all models and experimental results can be noticed. The results shown in the last sections show that the flow field has different behaviors with using each RWFs. To find out which of the three RWFs resulted in the most accurate results, error analysis should be done to find out which RWF had less deviation from the experimental results. C p distribution was chosen as the main criteria to compare since as explained earlier, the C p distribution gives a better understanding of the shape of the flow field over the body which is very important to simulate ice accumulation. In this work, the average absolute error between C p calculated from pressure measurements in wind tunnel and 240 the corresponding C p calculated from simulations using the The standard error of estimates: where σ avg is the standard error, C p,exp and C p,sim are coefficients of pressure of experiments and simulations respectively, and N is the number of experimental points.
As shown in the Figures 14a and 14b, it can be noticed that the higher the AoA value, the more error between simulation and 245 wind tunnel results except for some results at AoA = 4 o in case 1. While in Fig. 14a, we can see that the error values are not increasing with the increase of AoA in the same rates as shown in the previous two figures. These error results shows that all RWFs has a limited capabilities in simulating the detachment of the flow and accordingly give results deviated from actual results. Also, in case 3, where profile 2 was simulated, the flow should be highly unsteady which is not suitable for RANS simulations. However, the Momentum RWF managed to predict the C p distribution better than the other models except for In this work, three different rough wall functions were tested on iced wind turbine airfoils. The aim of this comparison was to find out which RWF will be the most suitable one to simulate ice accretion on wind turbine blades exposed to icing atmospheric conditions. To be able to apply these RWFs, DLR and Colebrook RWFs were implemented to OpenFOAM v6 CFD framework 255 along with the existing Momentum RWF. After that, two ice profiles that collected from wind site, molded to airfoil, and tested in the wind tunnel were smoothed by a cubic spline to find the equivalent smooth surface. Also, roughness parameters described in Sec. 3.2 were calculated and used to calculate the velocity shift value ∆u.
In this section, general remarks and discussions of the results are introduced. Then the conclusions from the outcomes of this work will be highlighted.

General remarks and discussion
Regarding the rough ice surface shown in Fig. 2, one should keep in mind that the ice formation process (or roughness formation in general) is a stochastic phenomenon. Which means, if the same airfoil exposed to the same atmospheric conditions, the ice profile resulting will not be exactly the same. Hence, to find the real smooth surface required for the simulations, a large number of ice profiles of the same atmospheric conditions should be studied and averaged to find the real average surface.

265
Also, each rough surface will give different attachment and reattachment bubbles locations, and might have different overall pressure distribution. However, the scope of this work is only to prove that the RWF approach results in good results compared to experiments. on the different turbulence parameters. From the geometry of the two ice profiles shown in Fig.2, it can be noticed that both profiles are slightly inclined downwards. This inclination forms a relatively big separation bubble behind the profile on the lower surface as shown in Fig. 15. That is why we can see deviations between simulations and experimental results of C p distribution in this region (i.e. the lower surface between x/c = -0.1 to 0.1) in most of the studied cases.
It can be noticed from the results that agreement of C l curves between simulations and experimental results doesn't necessarily 275 mean that flow field was accurately simulated. For example, in case 3, DLR RWF shows better C l agreement than other models at AoA = 4 o while from C p distribution, DLR RWF showed larger deviations from experimental results. The same happened with case 2 at AoA = 0 o and 4 o . The overall C l at these AoAs in the end had good agreement because the deviations on the upper and lower surfaces compensate each other which can be misleading in this case. To accurately access the C l , it should be studied together with the C p distribution curves.

280
It is noteworthy that all the RWFs used are only algebraic equations that express the behavior of the flow near that wall.
Accordingly, the three different RWFs have the same computational cost. However, the overall computational cost depends only on number of cells in the used computational grid and the turbulence model itself as discussed in Sec. 2.

Conclusions
The simulation cases shows that DLR RWF provides the best agreement of cases 1 and 2 simulations with experiments com-285 pared to the other models. While in case 3, Momentum RWF provides the least errors between experiments and simulations.In some cases, the simulations show fair agreement at locations that witness violent separation like concave areas between ice and airfoil. This fair agreement is expected due to using steady state RANS simulations.
From this work, we can conclude that in simulations of rough, iced airfoil profiles, the Momentum RWF should be used to for the simulation of the performance drop of wind turbine blades exposed to icing atmospheric conditions with minimum computational effort.  Competing interests. The authors declare no conflict of interest.
Disclaimer. TEXT