Assessment of Turbulence Models for Flow Simulation around a Wind Turbine Airfoil

Abstract

This investigation focuses on the application of the computational fluid dynamics tool FLUENT to the study of flows over the NACA 63–415 airfoil at various angles of attack. With the aim of selecting the most suitable turbulence model to simulate flow around ice-accreted airfoils, this work concentrates on assessing the prediction capabilities of various turbulence models on clean airfoils at the large angles of attack that cause highly separated flows to occur. The study was undertaken by conducting simulations with the one-equation Spalart-Allmaras (SA) model, the two-equation RNG 𝑘-𝜖 and SST 𝑘-𝜔 models, and the Reynolds stress model (RSM). Domain discretization was carried out using general quadrilateral grids generated with GAMBIT, the FLUENT preprocessing tool. Comparisons were made with available experimental data.

1. Introduction

Today, the worldwide concern for the environment can be expressed as a huge need for low entropy production during energy transformation processes. This essentially means low-polluting, low-wastage, and low energy-degradation methods and technologies. At the forefront of these is wind technology, and wind power plants constitute the world's fastest growing energy source [1]. Testament to this is the proliferation of wind farms that has occurred in recent years, and the planned installation of many others all around the world. However, in certain regions, wind energy extraction faces specific challenges imposed by climatic conditions.

In cold climates, the formation of ice on the wind turbine blades is one of the main concerns, as ice-induced roughness on the blades reduces lift and increases drag. This results in production losses that should be quantified, in order to obtain a picture of how much there is to be gained from the development of technologies designed to prevent ice accretion in cold regions.

The phenomenon of ice accretion and performance losses on wind turbine blades is a subject which is of particular interest to the wind turbine research community (see, e.g., [2–4]). The problem is that ice accretion on the blades will develop in a severe performance penalty, because large separated flow regions will develop at low angles of incidence (as they do on clean airfoils at high angles of incidence), followed by a significant lift reduction and an increase in drag. Few studies have analysed airfoils with large separated flows at Reynolds numbers characteristic of wind turbines [5, 6]. In fact, simulating the flow phenomenon around wind turbine blades is a challenging task, because that flow is three dimensional, unsteady, and turbulent. In wider spectra, it is useful to assess the ability of turbulence models to lead to an accurate prediction of lift and drag when separated flows occur.

In this investigation, only geometries contained in a two-dimensional space are addressed, and the flow is considered to be incompressible. This study concentrates on the aerodynamics of airfoils at low and high angles of attack.

To conduct the current computations, FLUENT [7], a finite-volume flow solver was applied, while the meshing process was achieved using the accompanying preprocessor, GAMBIT [8]. An initial grid with approximately 85 000 cells was generated. This mesh was then refined near the airfoil by applying FLUENT, which resulted in a typical grid of about 140 000 cells.

The 2D steady-state results on the NACA 63-415 airfoil are compared to data provided by the RISO National Laboratory [9]. This recognized research center carries out wind tunnel testing and also performs CFD calculations by applying their own code, EllipSys2D, to study and improve the efficiency of wind turbine blades.

2. Numerical Model

The numerical counterpart of the RISO experimental test is carried out on a 0.60 m chord, NACA 63–415 airfoil at Re=1.6×10⁶. The flow parameters at the inlet are a velocity 𝑉=40m/s and a turbulence intensity of 1%.

The meshed domain used for the calculations is shown in Figure 1. The top, bottom, and left boundaries were placed at a distance of 12 chords from the airfoil, while the right boundary was placed at 20 chords.

Figure 1 

A global view of the mesh.

Focusing on the flow simulation over iced airfoils, for which the characteristics and the calculation of the flow near the airfoil are even more crucial than for clean airfoils, as in the case of ice accretion prediction, we chose to address the viscous sublayer. Thus, a grid with a dimensionless distance 𝑦⁺ close to 1 was generated.

The airfoil pressure and suction sides were meshed with 312 uniformly distributed nodes, and 100 nodes with a growth ratio of 1.06 were placed between the airfoil and the far-field boundaries. Then, the mesh was refined near the airfoil to obtain the established wall distance, leading to a final mesh with approximately 140 000 cells. An idea is given in Figure 2 of a typical resulting mesh (detail showing the grid after being adapted in FLUENT).

Figure 2 

Near-wall view of the mesh.

A grid independence study was performed to verify that the solution will not change subsequent additional refinements. The study was based on simulations performed at an angle of incidence of 8°, with the lift and drag coefficients taken as the control parameters. Figures 3(a) and 3(b) show the computed values of the lift and drag parameters according to the various turbulence models. At the second grid level (80 000 cells), the aerodynamic parameters calculated after the RNG𝑘-𝜖 and SA models were applied showed no further variation, while slight differences remained following the application of SST𝑘-𝜔. The results obtained with the RSM continued to change up to the third grid level (140 000 cells). These calculations led us to conclude that it is reasonable to consider that the grid with 140 000 cells guarantees results on all four turbulence models, independently of grid size.

(a) Lift coefficient

(b) Drag coefficient

Figure 3 

Grid convergence.

3. Turbulence Models

Before going into the details of the flow simulation around the NACA 63–415 airfoil, let us recall the basic notions underlying the turbulence models used in this investigation, namely, the RNG (renormalization group) variation of the 𝑘-𝜖 model, the SST𝑘-𝜔 (Shear-Stress-Transport 𝑘-𝜔) model, the SA (Spalart-Allmaras) model, and the RSM (Reynolds stress model).

3.1. RNG 𝑘-𝜖

The RNG model incorporates modifications over the standard 𝑘-𝜖 [10] model, which improve its performance. These modifications were derived using the mathematical technique derived from what is called “renormalization group” theory [11]. The transport equations for the turbulent kinetic energy 𝑘 and the dissipation rate of the turbulent kinetic energy 𝜖 are defined as follows:𝜕+𝜕𝜕𝑡𝜌𝑘𝜕𝑥_𝑖𝜌𝑘𝑈_𝑖=𝜕𝜕𝑥_𝑗𝛼_𝑘𝜇_eff𝜕𝑘𝜕𝑥_𝑗+𝐺_𝑘−𝜌𝜖+𝑆_𝑘, (1)𝜕+𝜕𝜕𝑡𝜌𝜖𝜕𝑥_𝑖𝜌𝜖𝑈_𝑖=𝜕𝜕𝑥_𝑗𝛼_𝜖𝜇_eff𝜕𝜖𝜕𝑥_𝑗+𝐶_1𝜖𝜖𝑘𝐺_𝑘−𝐶_2𝜖+𝐶_𝜇(𝑆(𝑘/𝜖))³1−𝑆(𝑘/𝜖)/𝜂₀1+𝛽(𝑆(𝑘/𝜖))³𝜌𝜖²𝑘+𝑆_𝜖. (2)

The turbulent kinetic energy production term is defined as follows:𝐺_𝑘=−𝜌𝑢^′_𝑖𝑢^′_𝑗𝜕𝑈_𝑗𝜕𝑥_𝑖. (3)

The main modification is the addition of a second coefficient in the third term on the right side of (2). This coefficient becomes negative in regions of large strain rate (𝑆), reducing the destruction of epsilon, and, as a consequence, 𝑘 is also reduced, yielding lower turbulent viscosity values.

Other modifications are a differential equation to calculate the turbulent viscosity and the equation to calculate the inverse effective Prandtl numbers (𝛼_𝑘,𝛼_𝜖). The constant values for this model are 𝜂₀=4.38, 𝛽=0.012, 𝐶_𝜇=0.0845, 𝐶_1𝜖=1.42, and 𝐶_2𝜖=1.68.

The RNG𝑘-𝜖 model is only suited for zones far from the wall and cannot address the viscous sublayer directly. In order to be able to do this, the Enhanced Wall Treatment (EWT) option was devised. Based on the turbulent Reynolds number, the EWT model specifies a near-wall zone, where the eddy viscosity is calculated using the one-equation model of (Wolfstein [12]). Outside this zone, the eddy viscosity is calculated, as usual, by the RNG𝑘-𝜖 model. Then, a blending equation for the eddy viscosity is used to couple the two zones.

3.2. SA

The SA [13, 14] is a one-equation model that solves an empirical transport equation for the eddy viscosity 𝜈. This equation models the production, transport, diffusion, and destruction of the modified turbulent kinematic viscosity. One of its advantages is the simplicity with which the free stream and wall boundary conditions can be imposed. The transport equation for 𝜈 is𝜕𝜌𝜈+𝜕𝜕𝑡𝜕𝑥_𝑖𝜌𝜈𝑈_𝑖=𝐺_𝜈+1𝜎_𝜈𝜕𝜕𝑥_𝑗𝜈𝜕𝜈𝜇+𝜌𝜕𝑥_𝑗+𝐶_𝑏2𝜌𝜕𝜈𝜕𝑥_𝑗²−𝑌_𝜈+𝑆_𝜈. (4)

The production and destruction terms, 𝐺_𝜈 and 𝑌_𝜈, respectively, are defined as follows:𝐺_𝜈=𝐶_𝑏1𝜌𝜈𝑆+𝜅²𝑦²𝑓_𝜐2𝑌𝜈,_𝜈=𝐶_𝑤1𝜌𝑓_𝑤𝜈𝑦², (5) where 𝑆 is the mean strain rate, 𝑦 is the wall distance, and 𝑓_𝜐2 and 𝑓_𝑤 are damping functions. The constant values are 𝐶_𝑏1=0.1355, 𝐶_𝑏2=0.622, and 𝜈=2/3.

The SA model is able to resolve the viscous sublayer when the mesh is fine enough near the walls.

3.3. SST 𝑘-𝜔

The original 𝑘-𝜔 [15] model is reputed to be more accurate than 𝑘-𝜖 in the near-wall layers. It has been successful for flows with moderate adverse pressure gradients, but it has a 𝜔 equation which is very sensitive to the values of 𝜔 in the free stream. The SST [16] corrects this problem by solving the standard 𝑘-𝜖 in the far field and the standard 𝑘-𝜔 near the walls. To improve its performance for adverse pressure flows, the SST considers the effects of the transport of the turbulent shear stress in the calculations of the turbulent viscosity 𝜇_𝑡 and the turbulent Prandtl numbers 𝜎_𝑘 and 𝜎_𝜔. The transport equation for 𝑘 and the specific dissipation rate 𝜔 are the following:𝜕+𝜕𝜕𝑡𝜌𝑘𝜕𝑥_𝑖𝜌𝑘𝑈_𝑖=𝜕𝜕𝑥_𝑗𝜇𝜇+^𝑡𝜎_𝑘𝜕𝑘𝜕𝑥_𝑗+𝐺_𝑘−𝑌_𝑘+𝑆_𝑘,𝜕+𝜕𝜕𝑡𝜌𝜔𝜕𝑥_𝑖𝜌𝜔𝑈_𝑖=𝜕𝜕𝑥_𝑗𝜇𝜇+^𝑡𝜎_𝜔𝜕𝜔𝜕𝑥_𝑗+𝐺_𝜔−𝑌_𝜔+𝐷_𝜔+S_𝜔. (6)

The production terms of 𝑘 and 𝜔 are defined, respectively, as follows:𝐺_𝑘=min𝐺_𝑘,10𝜌𝛽^∗𝑘𝜔, (7)𝐺_𝑘 being calculated with (3),𝐺_𝜔=𝛼_𝜔𝜈_𝑡𝐺_𝑘, (8) and the dissipation terms are given by𝑌_𝑘=𝜌𝛽^∗𝑌𝑘𝜔,_𝜔=𝜌𝛽𝜔². (9)

The cross diffusion term that blends the two models is defined as follows:𝐷_𝜔=21−𝐹₁𝜌𝜎_𝑤,21𝜔𝜕𝑘𝜕𝑥_𝑗𝜕𝜔𝜕𝑥_𝑗, (10)𝛽^∗, 𝛽, 𝛼_𝜔, and 𝐹₁ are damping functions, and the constant value is 𝜎_𝑤,2=1.168.

3.4. RSM

The RSM [17–19] solves a transport equation for each term of the Reynolds stress tensor and an equation for the dissipation rate. In 2D, five transport equations are solved, in addition to the continuity and momentum equations. The model takes into account effects that one- and two-equation models are not able to consider, for example, the anisotropy of the turbulence. However, the Reynolds stress transport equations require the modeling of complex terms, which may compromise the model’s accuracy𝜕𝜌𝜕𝑡𝑢^′_𝑖𝑢^′_𝑗+𝜕𝜕𝑥_𝑘𝜌𝑈_𝑘𝑢^′_𝑖𝑢^′_𝑗𝜕=−𝜕𝑥_𝑘𝜌𝑢^′_𝑖𝑢^′_𝑗𝑢^′_𝑘+𝑝𝛿_𝑘𝑗𝑢^′_𝑖+𝛿_𝑖𝑘𝑢^′_𝑗+𝜕𝜕𝑥_𝑘𝜇𝜕𝜕𝑥_𝑘𝑢^′_𝑖𝑢^′_𝑗−𝜌𝑢^′_𝑖𝑢^′_𝑘𝜕𝑈_𝑗𝜕𝑥_𝑘+𝑢^′_𝑗𝑢^′_𝑘𝜕𝑈_𝑖𝜕𝑥_𝑘+𝑝𝜕𝑢^′_𝑖𝜕𝑥_𝑗+𝜕𝑢^′_𝑗𝜕𝑥_𝑖−2𝜇𝜕𝑢^′_𝑖𝜕𝑥_𝑘𝜕𝑢^′_𝑗𝜕𝑥_𝑘+𝑆_RMS. (11)

To resolve the viscous sublayer with the RSM, the Low-Re Stress-Omega was selected. This option models the pressure strain term (fourth term on the right) of the Reynolds stress equations using the LRR approach (Launder-Reece-Rodi) [19] and solves the 𝜔 equation [20]. This option gives the RSM model characteristics similar to those of the 𝑘-𝜔 model.

The setup of the RSM model requires more attention than the others. With the default options and second-order discretization for the momentum and for the turbulence equations, the model showed early convergence problems at low angles of attack, and, in certain cases, the viscosity ratio reached unrealistic values (beyond 10⁵). After several tests, we empirically determined the setup to be used for the simulations; at low (8° or less) and high (26° or more) angles of attack, the simulations were performed with the default setup. For the rest of the calculations (10°–24°), we disabled the Shear Flow Corrections option. All the simulations were performed with a second-order discretization for the momentum and a first-order discretization for the Reynolds stresses.

The Shear Flow Corrections option enables a function that modifies the values of the dissipation term in the 𝜔 equation, the first pressure strain term, and the dissipation rate of the Reynolds stresses. Studying the influence of the turbulent kinetic energy 𝑘 revealed that the main effect of this function is to reduce 𝑘 (𝑘 is calculated with the diagonal of the Reynolds stresses).

4. Numerical Simulations

An initial set of computations was carried out using the steady-state formulation with the SA, SST𝑘-𝜔, RNG𝑘-𝜖, and RSM models for the 0°–28° angle of attack range.

Figures 4(a) and 4(b) show the computed results using averaging to estimate the aerodynamic coefficients when oscillations occurred. In these figures, RISO Laboratory data are also shown, consisting of wind tunnel records, as well as numerical data calculated with the EllipSys2D program using the SST𝑘-𝜔 model.

(a) Lift coefficient

(b) Drag coefficient

Figure 4 

Lift and drag coefficients.

At a low angle of attack (up to 8°), the four models show good agreement with the experimental data. In the 8° to 16^∘ range, lift and drag coefficients begin to show weak oscillations and numerical results that differ from the experimental results. At high angles (24^∘ and over), the four models give similar results, but overpredicted. In this range, the SST𝑘-𝜔 model gives the smoothest curve.

The SA model gives the best prediction of the maximal lift angle; however, the lift value is overpredicted. Also, this model gives the smoothest curves for the coefficients, except at high angles of attack.

The RNG𝑘-𝜖 model predicts the highest values of both coefficients and a maximal lift angle that is far from the experimental value. Between 22^∘ and 24^∘, both coefficients change abruptly.

The SST𝑘-𝜔 model shows similar values to the SA model, except in the model's prediction of maximal lift. At high angles of attack, this model gives the smoothest curve.

The RSM model does not provide better results than the alternatives. It shows a drastic change in the coefficient behavior at 18^∘, which means that the model faces problems with handling large recirculation zones. Moreover, the simulations take much more time, since more equations are solved.

Figure 5 shows the computed streamlines at the trailing edge at an angle of attack of 16^∘ for each turbulence model. The RNG𝑘-𝜖 and RSM models predict the smallest zones that coincide with the highest 𝐶_𝑙 and the lowest 𝐶_𝑑 predicted at this angle. In the case of the RNG𝑘-𝜖 model, the separation point is almost at the same place as the RSM, but the streamlines are barely perturbed, which may be because the 𝑘-𝜖 models overpredict the turbulence kinetic energy production, and momentum exchange allows the flow to remain close to the walls (Figures 4(a) and 4(b)). The SA and SST𝑘-𝜔 models predict the largest zones that are of similar size. However, the SST𝑘-𝜔 and RSM models predict unstable regions with vortex separation, which may explain the difficulty the models have in reaching convergence at this angle. The similarity between the SA-SST𝑘-𝜔 and RSM-RNG𝑘-𝜖 models at the trailing edge is also present in the 𝐶_𝑝 curves shown in Figure 7.

(a) SA model (separation at 𝑥/𝑐=0.52)

(b) RNG𝑘-𝜖 model (separation at 𝑥/𝑐=0.80)

(c) SST𝑘-𝜔 model (separation at 𝑥/𝑐=0.50)

(d) RSM model (separation at 𝑥/𝑐=0.88)

Figure 5 

Streamlines after the application of various turbulence models at 16∘.

Regarding the pressure coefficient, all the models give curves with similar shapes and values close to the experimental data. Figure 6 shows the pressure coefficient distribution at 0°. At this angle, the four models fully agree with the values given by RISO, except near the leading edge and at the stagnation point, where 𝐶_𝑝 is higher than the experimental value. Concerning the lift and drag, the four models give a good approximation of lift, and only the RSM gives a good estimate of drag, while the rest of models overpredict it.

Figure 6 

Pressure coefficient at 0° following application of various turbulence models.

Figure 7 

Pressure coefficient at 16∘ following application of various turbulence models.

Figure 7 shows the results for an angle of 16^∘. In this case, the pressure side agrees very well with the experimental values. On the suction side, however, the minimum 𝐶_𝑝 is underpredicted, which means that a higher suction acts on the upper side of the airfoil. There is an obvious relation between the peak found at this point and the overprediction of 𝐶_𝑙. From Figures 4 and 7, we see that the RNG𝑘-𝜖 model leads to the highest peak, which corresponds to the highest value of 𝐶_𝑙, while the SA leads to the lowest peak and the lowest value of 𝐶_𝑙. At the trailing edge, it can be seen that 𝐶_𝑝 depends on the size of the circulation zone. The smallest zone (RNG𝑘-𝜖) has 𝐶_𝑝 values close to zero and the largest zone (SST𝑘-𝜔) has the highest 𝐶_𝑝 values.

In order to understand the situation from a different perspective, a relative mass imbalance, defined as the mass imbalance on each cell divided by the mass of the cell, was introduced. Since this index can be positive or negative (mass source or mass sink), its absolute value was used to calculate the mean value over the whole domain.

For a flow having an angle of incidence of 24^∘ and using the SST𝑘-𝜔 model, the mean value was monitored during 10 000 iteration steps. Figures 8(a) and 8(b) show the cells having a relative mass imbalance of more than ±0.5% at the iteration where the mean is at its highest and at its lowest value. This representation shows regions surrounding the airfoil in which a steady-state solution was not reached, suggesting that transient calculations need to be considered.

(a) Minimal mean

(b) Maximal mean

Figure 8 

Relative mass imbalance contours.

The reasonable approximation of the 𝐶_𝑙 and 𝐶_𝑑 values obtained by all the models may be explained by the fact that the relative mass imbalance remains low. It also indicates that steady-state solutions are a reasonable starting point for the transient calculations.

5. Conclusions

This paper has presented a study of the prediction of two-dimensional flows over a wind turbine airfoils. The methodology was based on the application of the commercial program FLUENT.

The intention was to provide wind turbine researchers with an assessment of turbulence modeling for the calculation of the fundamental lift and drag parameters for flows over wind turbine blade sections operating under icing conditions. Four turbulence models, the one-equation Spalart-Allmaras (SA) model, the two-equation RNG𝑘-𝜖 and SST𝑘-𝜔 models, and the Reynolds stress model (RSM), were applied to predict flows over a clean NACA 63-415 airfoil at angles of attack ranging from 0° to 28^∘.

From these results, we estimate that the RSM is the least appropriate for solving this type of flow problem. It has shown to be much more sensitive in the presence of moderate to large recirculation zones, and it predicts high viscosity ratio values, requiring the readjustment of the setup options each time to prevent this problem from occurring. The calculations conducted with this model require more computational time with no noticeable improvement, when compared with the three alternatives for handling turbulence.

The RNG𝑘-𝜖 model gives the highest prediction of lift and maximal lift angle. At 16^∘, the model also predicts the smallest recirculation zone and the separation point that is closest to the trailing edge.

The results obtained with the SA and the SST𝑘-𝜔 models are very similar over the 0° to 16^∘ range. In this work, these models always provided lower and less oscillatory aerodynamic coefficients and did so for the entire range of angles of attack tested. This behavior is not related to a faster or slower convergence rate.

The SA and SST𝑘-𝜔 models were also found to be the most adequate for simulating the current flows. As seen in Figures 4(a) and 4(b), they are shown to have a better agreement with the solution shape given by the experimental data. Nevertheless, in this study, we noted that the SA model has a tendency to predict stable recirculation zones. Also, in this model, turbulence decay is only based on wall distance and velocity gradients and cannot handle turbulence dissipation, as the SST𝑘-𝜔 model does. This implies that turbulence intensity remains constant in the wake. Although the other models are prone to showing similar nonphysical recirculation zones, this is the most unfavorable case.

Recirculation zones are present in both clean and iced wind turbine airfoils. In the case of ice accretion simulation, this zone impacts the water impingement calculation over the suction side. Thus, a good approximation of flow in this zones is important in numerical studies of wind turbine flow.

Therefore, despite the fact that the SA gives better predictions near maximal lift, we conclude that the SST𝑘-𝜔 model is the best suited to simulating flow around both clean and iced wind turbine airfoils.

Acknowledgments

The authors gratefully acknowledge the financial support of the NSERC, through a strategic grant, and of FQRNT, through a team project grant.