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 k-... and SST k-ω 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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Fernando Villalpando 1 and Marcelo Reggio 1 and Adrian Ilinca 2

Recommended by Guan Yeoh

1, École Polytechnique de Montréal, Campus de l'Université de Montréal, 2500 Chemin du Polytechnique, Montréal, QC, H3T 1J4, Canada
2, Laboratoire de Recherche en énergie éolienne, Université du Québec à Rimouski, 300 Allée des Ursulines, Rimouski, QC, G5L 3A1, Canada

Received 23 September 2010; Revised 10 February 2011; Accepted 15 February 2011

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×1⁰⁶ . The flow parameters at the inlet are a velocity V=40 m/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.

[figure omitted; refer to PDF]

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 ^y+ 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.

[figure omitted; refer to PDF]

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 k-... and SA models were applied showed no further variation, while slight differences remained following the application of SST k-ω . 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.

Grid convergence.

(a) Lift coefficient

[figure omitted; refer to PDF]

(b) Drag coefficient

[figure omitted; refer to PDF]

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 k-... model, the SST k-ω (Shear-Stress-Transport k-ω ) model, the SA (Spalart-Allmaras) model, and the RSM (Reynolds stress model).

3.1. RNG k-...

The RNG model incorporates modifications over the standard k-... [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 k and the dissipation rate of the turbulent kinetic energy ... are defined as follows: [figure omitted; refer to PDF] [figure omitted; refer to PDF]

The turbulent kinetic energy production term is defined as follows: [figure omitted; refer to PDF]

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 (S ), reducing the destruction of epsilon, and, as a consequence, k 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 (_αk , _α... ). The constant values for this model are _η0 =4.38 , β=0.012 , _Cμ =0.0845 , _C1... =1.42 , and _C2... =1.68 .

The RNG k-... 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 k-... 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 [figure omitted; refer to PDF]

The production and destruction terms, _Gν and _Yν , respectively, are defined as follows: [figure omitted; refer to PDF] where S is the mean strain rate, y is the wall distance, and _f[upsilon]2 and _fw are damping functions. The constant values are _Cb1 =0.1355 , _Cb2 =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 k-ω

The original k-ω [15] model is reputed to be more accurate than k-... 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 k-... in the far field and the standard k-ω 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 _μt and the turbulent Prandtl numbers _σk and _σω . The transport equation for k and the specific dissipation rate ω are the following: [figure omitted; refer to PDF]

The production terms of k and ω are defined, respectively, as follows: [figure omitted; refer to PDF] _Gk being calculated with (3), [figure omitted; refer to PDF] and the dissipation terms are given by [figure omitted; refer to PDF]

The cross diffusion term that blends the two models is defined as follows: [figure omitted; refer to PDF] ^β* , β , _αω , and _F1 are damping functions, and the constant value is _σw,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 [figure omitted; refer to PDF]

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 k-ω 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 k revealed that the main effect of this function is to reduce k (k 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 k-ω , RNG k-... , 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 k-ω model.

Lift and drag coefficients.

(a) Lift coefficient

[figure omitted; refer to PDF]

(b) Drag coefficient

[figure omitted; refer to PDF]

At a low angle of attack (up to 8°), the four models show good agreement with the experimental data. In the 8° to 1^6° range, lift and drag coefficients begin to show weak oscillations and numerical results that differ from the experimental results. At high angles (2^4° and over), the four models give similar results, but overpredicted. In this range, the SST k-ω 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 k-... model predicts the highest values of both coefficients and a maximal lift angle that is far from the experimental value. Between 2^2° and 2^4° , both coefficients change abruptly.

The SST k-ω 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 1^8° , 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 1^6° for each turbulence model. The RNG k-... and RSM models predict the smallest zones that coincide with the highest _Cl and the lowest _Cd predicted at this angle. In the case of the RNG k-... model, the separation point is almost at the same place as the RSM, but the streamlines are barely perturbed, which may be because the k-... 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 k-ω models predict the largest zones that are of similar size. However, the SST k-ω 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 k-ω and RSM-RNG k-... models at the trailing edge is also present in the _Cp curves shown in Figure 7.

Streamlines after the application of various turbulence models at 1^6° .

(a) SA model (separation at x/c=0.52 )

[figure omitted; refer to PDF]

(b) RNG k-... model (separation at x/c=0.80 )

[figure omitted; refer to PDF]

(c) SST k-ω model (separation at x/c=0.50 )

[figure omitted; refer to PDF]

(d) RSM model (separation at x/c=0.88 )

[figure omitted; refer to PDF]

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 _Cp 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 omitted; refer to PDF]

Figure 7: Pressure coefficient at 1^6° following application of various turbulence models.

[figure omitted; refer to PDF]

Figure 7 shows the results for an angle of 1^6° . In this case, the pressure side agrees very well with the experimental values. On the suction side, however, the minimum _Cp 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 _Cl . From Figures 4 and 7, we see that the RNG k-... model leads to the highest peak, which corresponds to the highest value of _Cl , while the SA leads to the lowest peak and the lowest value of _Cl . At the trailing edge, it can be seen that _Cp depends on the size of the circulation zone. The smallest zone (RNG k-... ) has _Cp values close to zero and the largest zone (SST k-ω ) has the highest _Cp 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 2^4° and using the SST k-ω 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.

Relative mass imbalance contours.

(a) Minimal mean

[figure omitted; refer to PDF]

(b) Maximal mean

[figure omitted; refer to PDF]

The reasonable approximation of the _Cl and _Cd 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 k-... and SST k-ω 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 2^8° .

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 k-... model gives the highest prediction of lift and maximal lift angle. At 1^6° , 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 k-ω models are very similar over the 0° to 1^6° 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 k-ω 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 k-ω 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 k-ω 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.

References

[1] World Wind Energy Report 2009 , World Wind Energy Association, Germany, 2009.

[2] M. B. Bragg, A. P. Broeren, L. A. Blumenthal, "Iced-airfoil aerodynamics," Progress in Aerospace Sciences , vol. 41, no. 5, pp. 323-362, 2005.

[3] H. Seifert, F. Richert, "Aerodynamics of iced airfoils and their influence on loads and power production," in Proceedings of the European Wind Energy Conference, Dublin Ireland, October 1997.

[4] C. Mayer, A. Ilinca, G. Fortin, J. Perron, "Wind tunnel study of electrothermal de-icing of wind turbine blades," The International Journal Society of Offshore and Polar Engineers , vol. 17, pp. 182-188, 2007.

[5] S. L. Yang, Y. L. Chang, O. Arici, "Navier-Stokes computations of the NREL airfoil using a k-ω turbulent model at high angles of attack," Journal of Solar Energy Engineering, Transactions of the ASME , vol. 117, no. 4, pp. 304-310, 1995.

[6] Y. Chen, Y. Zhiquan, L. Deyuan, H. Fupeng, "Numerical simulation of large angle-of-attack separated flows over airfoils of HAWT rotors," Wind Engineering , vol. 30, no. 1, pp. 35-42, 2006.

[7] Fluent Inc. Fluent 6.3.26 User's Guide , Fluent Inc., 2005.

[8] Fluent Inc. Gambit 2.2.30 User's Guide , Fluent Inc., 2004.

[9] C. Bak, P. Fuglsang, J. Johansen, I. Antoniou, "Wind tunnel tests of the NACA 63-415 and a modified NACA 63-415 airfoil,", no. Risø-R-1193, Risø National Laboratory, Roskilde, Denmark, 2000.

[10] W. P. Jones, B. E. Launder, "The prediction of laminarization with a two-equation model of turbulence," International Journal of Heat and Mass Transfer , vol. 15, no. 2, pp. 301-314, 1972.

[11] V. Yakhot, S. A. Orszag, "Renormalization group analysis of turbulence. I. Basic theory," Journal of Scientific Computing , vol. 1, no. 1, pp. 3-51, 1986.

[12] M. Wolfstein, "The velocity and temperature distribution of one-dimensional flow with turbulence augmentation and pressure gradient," International Journal of Heat and Mass Transfer , vol. 12, pp. 301-318, 1969.

[13] P. Spalart, S. Allmaras, "A one-equation turbulence model for aerodynamic flows," AIAA Paper , no. 92-0439, 1992.

[14] P. R. Spalart, S. R. Allmaras, "One-equation turbulence model for aerodynamic flows," Recherche Aerospatiale , no. 1, pp. 5-21, 1994.

[15] D. C. Wilcox, "Reassessment of the scale-determining equation for advanced turbulence models," AIAA Journal , vol. 26, no. 11, pp. 1299-1310, 1988.

[16] F. R. Menter, "Two-equation eddy-viscosity turbulence models for engineering applications," AIAA Journal , vol. 32, no. 8, pp. 1598-1605, 1994.

[17] M. M. Gibson, B. E. Launder, "Ground effects on pressure fluctuations in the atmospheric boundary layer," Journal of Fluid Mechanics , vol. 86, no. 3, pp. 491-511, 1978.

[18] B. E. Launder, "Second-moment closure: present ... and future?," International Journal of Heat and Fluid Flow , vol. 10, no. 4, pp. 282-300, 1989.

[19] B. E. Launder, G. J. Reece, W. Rodi, "Progress in the development of a reynolds-stress turbulence closure," Journal of Fluid Mechanics , vol. 68, no. 3, pp. 537-566, 1975.

[20] D. Wilcox, "Turbulence modeling: an overview," AIAA Paper , no. 2001-0724, 2001.