Abstract
The very first step in the simulation of ice accretion on a wind turbine blade is the accurate prediction of the flow field around it and the performance of the turbine rotor. The paper addresses this prediction using RANS equations with a proper turbulence model. The numerical computation is performed using a commercial CFD code, and the results are validated using experimental data for the 3D flow field around the NREL Phase VI HAWT rotor. For the flow simulation, a rotating reference frame method, which calculates the flow properties as time-averaged quantities, has been used to reduce the time spent on the analysis. A basic grid convergence study is carried out to select the adequate mesh size. The two-equation turbulence models available in ANSYS FLUENT are compared for a 7 m/s wind speed, and the one that best represents the flow features is then used to determine moments on the turbine rotor at five wind speeds (7 m/s, 10 m/s, 15 m/s, 20 m/s, and 25 m/s). The results are validated against experimental data, in terms of shaft torque, bending moment, and pressure coefficients at certain spanwise locations. Streamlines over the cross-sectional airfoils have also been provided for the stall speed to illustrate the separation locations. In general, results have shown good agreement with the experimental data for prestall speeds.
1. Introduction
Prediction of the aerodynamic loads on a horizontal axis wind turbine (HAWT) is both an important and a complex process that takes place during the design stage. It is important because it is directly related to crucial characteristics of the wind turbine, such as its power curve, structural loads, and noise generation. The power curve determines the energy output of the wind turbine, and therefore its cost effectiveness, structural loads determine the size of the various wind turbine components and materials to be used, and noise generation is a consideration for the location of a wind turbine and should be kept to a reasonable level. The complexity of a flow simulation around a wind turbine blade to determine aerodynamic loads comes from large angles of attack, a variety of Reynolds numbers resulting from very long blades, and rotational effects. Further complicating the issue is the atmospheric boundary layer, variable wind speeds along the blade, and interaction of the rotor with the nacelle and the tower.
Although 2D and quasi-3D methods like blade element momentum analysis may predict aerodynamic loads up to a certain level, 3D analysis is required to capture all the turbulence features that are inherently three dimensional. As computational cost drops thanks to technological development, the use of computational fluid dynamics (CFD) for wind turbine design and analysis is becoming increasingly widespread and results in a better understanding of the aerodynamic phenomena on the rotor flow field.
Studies on wind turbine aerodynamics in the literature focus on different aspects of the science, like performance analysis, fluid-structure interaction, acoustics, and icing. Investigation into performance analysis is primarily aimed at estimating aerodynamic loads, and the effect of various parameters on these loads. Duque et al. [1] explored the ability of various methods to predict wind turbine power and aerodynamic loads. Results showed that all the methods, namely blade element momentum (BEM), vortex lattice, and reynolds averaged navier stokes (RANS), perform well for prestall regimes. They found that the RANS code OVERFLOW, although not perfect, gives better predictions of power production for stall and poststall regime modeling than other methods. Another study, based on CFD analysis (Sorensen et al. [2]), showed that performance and wind turbine load predictions are very accurate, except in the stall region. A commercial wind turbine company, Siemens, analyzed their own large-scale wind turbine using a commercial CFD code, ANSYS-CFX, with transition and fully turbulent models [3]. Transition models improve drag prediction but overestimate the lift compared to fully turbulent models. Benjanirat et al. [4] used CFD with various turbulence models (Boldwin-Lomax, Spalart-Allmaras, and k-ε), with and without wall corrections, on an experimental wind turbine. The k-ε model with wall correction yielded the best results when compared with experimental data. A more recent study, conducted by Sezer-Uzol and Long [5] using a generic CFD code, PUMA2, have shown that time accurate inviscid results are also compatible with the experimental data.
Predicting the effects of tower, nacelle, and anemometer on the rotor flow field have also been investigated in several works. Smaili and Masson [6] and Zahle and Sorensen [7] both concluded that CFD is an effective tool for evaluating the influence of the nacelle, even for different positions and alignments of the wind turbine blade. A successful study of the interaction between rotor and tower [8] emphasizes the importance of CFD in solving complex flow phenomena. All these studies show that, even though CFD requires more resources for unsteady problems, it is cheaper than conducting full-scale or scaled wind turbine experimental analysis but yet provides sufficiently accurate results. By improving its capability to simulate widely separated flows through better turbulence modeling near the wall and in the flow field, CFD tools are becoming more and more accurate for the aerodynamic and aeroelasticity analysis of wind turbine blades.
Even though the final objective of the current study is to perform icing simulation on wind turbine blades, the characteristics of the flow field and the integrated loads of the clean blade are nevertheless required. To obtain these data, a grid study is performed to represent the flow field accurately with a minimum grid size. Three different grids, with 1.6 M, 1.9 M, and 2.2 M cells, are tested using the commercial CFD code (ANSYS-FLUENT [9]) at a prestall speed of 7 m/s and the results compared with experimental data. Once the grid convergence analysis has been conducted, various two-equation turbulence models are tested to find the one that best fits the experimental NREL Phase VI HAWT results. Thus, for each turbulence model, standard k-ε, RNG k-ε, realizable k-ε, and SST k-ω simulations are performed at a wind speed of 7 m/s, characterized by a reduced stall region on the rotor blade, to obtain pressure distributions and integrated aerodynamic loads on the blade. Results are compared with experimental data to find the most suitable turbulence model. Once that model has been selected, new simulations are performed throughout the operational wind speed range of the wind turbine.
2. Methodology
2.1. Experimental Data
For this study, we selected one of the unsteady aerodynamics experiments (UAEs), the NREL Phase VI wind turbine, conducted by the National Renewable Energy Laboratory (NREL) in the NASA-Ames wind tunnel at Moffett Field, California, in 2000 [10, 11]. The NREL organized the data used in this paper. The advantage of using a wind tunnel of such gigantic proportions (24.4 m by 36.6 m, or 80 ft by 120 ft) to perform a full scale 10 m diameter wind turbine test is that the blockage effect is not significant.
The general characteristics of the NREL Phase VI wind turbine are given in Table 1. The S809 airfoil is used for the entire blade. Linear chord and nonlinear twist distributions are shown in Figure 1. The blade root is completely circular up to 0.883 m, and, from this point up to 1.2573 m, the shape smoothly transforms from a circle into an airfoil configuration. The 3D CAD model of the blade was generated using the commercial software Rhinoceros 4.0 [12]. Figure 2 shows the geometry.
(a)
(b)
(a)
(b)
From the experimental data, the upwind, 3° pitch, and nonyaw configurations were selected as validation data. The experimental parameters to be validated are the pressure coefficients at pressure tap locations (30%, 46.6%, 63.3%, 80%, and 95% span), the low-speed shaft torque, and the bending moment on the blade throughout the operational wind speed range. The experimental data are provided for a 30 second time period, and time-averaged data are used for the comparison.
For the selected data, the wind turbine rotation direction is counter-clockwise when viewed from upwind. The cone angle is 0°. The rotor speed is 72 rpm, and the pitch angle is 3°. The pitch angle is defined at 75% of the span, 0.75, and the pitch axis is given as 30% of the chord, 0.3.
2.2. Numerical Study
Today, CFD tools, whether commercial or developed inhouse, are routinely used for investigating the aerodynamics and fluid-structure interactions around numerous configurations. Of the commercial tools available, ANSYS FLUENT 12.0.16 is used in the current study. The package is a finite volume-based solver, which allows both structured and unstructured grids to discretize the computational domain [9]. The software allows transient calculations as well as steady-state computations to be performed. Parallel computation capabilities are also offered to handle large meshes and to reach solutions in a reasonable time.
For this study, an incompressible, steady, Reynolds-averaged Navier Stokes (RANS) model is applied to solve the problem on a “rotating reference frame.” The basic idea behind the rotating reference frame is the assumption that it is the flow field that rotates, and not the rotor, which means that an unsteady flow field turns into a steady flow with respect to the rotating reference frame. This approach simplifies the problem in terms of boundary conditions (no sliding mesh is required), computational cost, and postprocessing results.
To reduce the complexity of the problem and the time spent on computational analysis, we made a number of hypotheses prior to numerical analysis. Initially, effects of the tower and the nacelle are ignored to reduce mesh size and the complexity of the problem. Although they may affect the performance and flow field of the rotor, they are neglected in many studies when they are not the focus of the research [1–4]. Moreover, on the assumption that the flow field is 180° axisymmetric, one blade is modeled instead of two, and a “rotational periodic” boundary condition is applied, in order to reduce the computational cost. As explained above, since the blockage effect of the wind tunnel is minimal, the wind tunnel walls are not taken into account in the numerical domain. Finally, although the experimental data are time dependent, numerical analysis is conducted at steady state with the rotating reference frame model, which predicts time-averaged quantities. As a result, the time-averaged quantities of the experimental data and the numerical results can be compared.
Two-equation turbulence models have been widely used to simulate the flow field in engineering applications. As the name implies, these models have two independent transport equations, one for turbulent kinetic energy, and the other for turbulent dissipation or specific dissipation rate. Two equation models are complete, which means that no additional equations are needed to model the turbulence, and that they both depend on the Boussinesq assumption [13]. Details of turbulence models that are applied in this study are briefly explained below.
2.2.1. Standard k-ε Turbulence Model
The transport equations of the k-ε models are based on the turbulent kinetic energy, k, and the dissipation rate, ε. The simplest, the standard k-ε model, proposed by Launder and Spalding [14], is based on the assumption that flow is fully turbulent. This model gives better results for fully turbulent flows. The transport equations for the standard k-ε model are as follows: The turbulent viscosity is calculated using k and as follows: The model coefficients, which are empirically determined, are given as in [14]:
2.2.2. RNG k-ε Turbulence Model
A more developed k-ε turbulence model, RNG k-ε, which is based on renormalization group theory [15], has correction terms for swirling flow, low Reynolds number flow, and flow with high velocity gradients. The transport equations of the RNG k-ε model are very similar to those of the standard k-ε model, except as shown below: where The model coefficients, which are derived analytically by RNG, are given as in [15]: As can be seen from (5), the coefficient of dissipation term is modified to provide better adaptability of the model to rapid strained flows [9]. Moreover, the effective viscosity term in both transport equations improves the model for low Reynolds numbers and near-wall regions. For swirled flows, the RNG k-ε model of ANSYS FLUENT [9] has an optional correction model that calculates turbulent viscosity as a function of swirl strength. For this study, the swirl correction is enabled as the flow field rotates.
2.2.3. Realizable k-ε Turbulence Model
The realizable k-ε model, proposed by Shih et al. [16], is a new turbulent viscosity model that accounts for rotation and strain in the flow. The transport equations of the model are the following: where The realizable k-ε model is significantly different from other k-ε models, in that the turbulent viscosity coefficient depends on mean strain and rotation rates, turbulent kinetic energy, and energy dissipation, the use of which ANSYS FLUENT recommends for better prediction of the turbulent viscosity [9]. The model coefficients of the realizable k- εmodel are given as in [9]:
2.2.4. SST k-ω Turbulence Model
This model, which was presented by Menter [17], is a combination of the standard k-ω model and the transformed k-ε model. It benefits from the advantages of both these turbulence models in different flow regions, activating the standard k-ω model near the wall and the k-ε model away from the surface. The SST k-ω model has been found to be more accurate than other turbulence models [9] where the term is
The model coefficients are given as in [9]:
2.2.5. Computational Domain
The geometries of the wind turbine and the domain for the numerical analysis are generated using Rhinoceros [13]. Figure 3 illustrates the computational domain. The inlet and outlet boundaries are placed at 3 times and 5 times the diameter from the wind turbine, respectively. Since an axisymmetric flow field is required for using the rotating reference frame and the rotational periodic boundary condition, a semicylindrical domain has been generated. To improve accuracy, a second small semicylindrical domain around the blade has also been generated, in order to refine the mesh in this region, as illustrated in Figure 3.
An unstructured mesh was chosen for the domain discretization. The surface mesh was generated using GAMBIT, formerly the companion software of FLUENT. For the surface mesh, the size functions of GAMBIT that enable the user to control the grid distribution over the surface were used. The number of nodes was maximized near the leading edge, where the large pressure gradients exist. “Meshed” size functions provided smooth transactions of mesh size from edges to face, whereas “curvature” size functions provided better representation of curved surfaces by keeping the mesh size to a minimum, as seen in Figure 4(a). To resolve the boundary layer, 20 prism layers normal-to-surface elements on the blade were generated. The initial height selected was 10-5 m, which guarantees y+ < 10 on the blade surface. Three different grids with mesh sizes of 1.6 M, 1.9 M, and 2.2 M were used. The surface mesh and prism layer are shown in Figures 4(a) and 4(b), respectively.
(a)
(b)
2.2.6. Boundary Conditions
For the inlet boundary condition, the velocity normal to the inlet boundary and the turbulence intensity, which are provided by experimental data, were imposed. The outer cylindrical domain is also treated as an inlet, and the same boundary conditions are applied. For the outflow boundary, the pressure outlet boundary condition is enforced by applying atmospheric pressure, as the flow is far from the wind turbine. For the inner surfaces, as explained above, rotational periodic boundary conditions are applied. Finally, the blade surface was treated as a no-slip wall boundary condition; that is, a zero velocity is imposed.
3. Results
In this section, results of the comparison of 3D numerical studies and the UAE data [11] are presented, starting with the results of the grid optimization study. To find the most convenient turbulence model, various two-equation turbulence models available in ANSYS FLUENT were applied, and the pressure coefficients were compared with the UAE data. The models that are available in the software are standard k-ε, RNG k-ε, realizable k-ε, standard k-ω, and SST k-ω. Initially, the results of these comparisons are presented for a 7 m/s wind speed and preselected spanwise locations, as explained in the Section 2. As will be shown below, the k-ω SST turbulence model fits the UAE data better, and so this model was selected for the rest of simulations at higher wind speeds. To further verify the numerical study, moments on the wind turbine blade, that is, low speed shaft torque and root flap moment, were compared. Finally, the numerical tool’s ability to predict separated flow was shown by comparing prestall and stall velocity simulations.
3.1. Grid Study
To ensure the accuracy of the results and to keep the computational cost to a minimum, a basic grid convergence study was performed by generating three different grids: a coarse grid with 1.6 M elements, a medium grid with 1.9 M elements, and a fine grid with 2.2 M elements. One noticeable difference between the grids is in the number of elements on the leading edge of the blade, as can be seen in Figure 5.
(a)
(b)
(c)
The computed results are compared against experimental data at three different spanwise locations: root, midspan, and blade tip, as illustrated in Figures 6(a), 6(b), and 6(c), respectively. As shown, the coarse grid solution, represented by the cross-symbol, overpredicts the pressure coefficient for all spanwise locations, whereas medium and fine grids provide solutions that agree well with the experimental data. After midchord, where the pressure gradients are relatively low, the coarse grid leads to a reasonable prediction of the pressure coefficient on the blade. In spite of the improved behavior in this region, it is clear that, globally, the coarse grid does not yield satisfactory results. In terms of the solutions obtained using medium and fine meshes, it can be seen that the pressure coefficient calculated with the fine mesh is only slightly better than that computed with the medium-sized grid. Because of this minor difference, we conclude that the medium-sized grid can be selected for the rest of the calculations while still keeping the computational cost relatively low.
(a)
(b)
(c)
3.2. Comparison of Turbulence Models
The results obtained with the above-referenced turbulence models and medium-sized grid are compared at 7 m/s, which is the pre-stall, and therefore stable, region over the operational range. The pressure coefficient distribution over the chord is presented for three spanwise locations as root 0.3, mid-span 0.63, and tip 0.95. The pressure coefficient is obtained by (14), where represents the wind speed, is the local rotational speed at that spanwise location, and indicates the local chord. Comparison of the calculations obtained with the various turbulence models against experimental data are shown in Figures 7(a)–7(c) At first glance, with the exception of the standard k-εmodel, all the turbulence models seem to predict the pressure coefficient well at the pressure side and at the suction side. The standard k-ε model overestimates the pressure at the leading edge.
(a)
(b)
(c)
At the root of the blade (Figure 7(a)), where the rotational speed is relatively very low, the compatibility of all the turbulence models seems good. However, if the results are evaluated in detail at the leading edge and tip of the section, differences among the models become clearer. As seen in detail in Figure 7(a), the pressure coefficient predicted with the SST k-ω model is the best when compared with the experimental data, and the computed pressure coefficient is the worst when using the standard k-ε model.
At the mid-span location (Figure 7(b)), the calculations with the standard k-ε model further deviate from the experimental data at the leading edge of the blade. The detailed comparisons of the pressure coefficients are shown at both the suction and pressure surfaces of the leading edge. The results of the other models are close at this location.
The ability of a turbulence model to accurately simulate flow becomes evident at the tip of the blade, where the rotational speed is the highest. As seen from the details in Figure 7(c), the predictions provided by the SST k-ω model are closer to the experimental data, and those obtained with the standard k-ε model deviate the most. Moreover, it is worth noting that the differences between the computations performed with the realizable k-ε and RNG k-ε models are very small, and these solutions appear to be the same in Figure 7(a)–7(c). Although there is no separation at all at this location at 7 m/s, the prediction capability of all the turbulence models decreases visibly compared to that at the root and mid-span locations.
As a result, among the turbulence models, the best overall compatibility with experimental data is shown by the SST k-ω model for all spanwise locations. Similar conclusions have been drawn by Villalpando et al. [18, 19], who analyzed the flow over an NACA 63415 airfoil and showed that although all the turbulence models predict the lift, drag, and pressure coefficients well, the k-ω SST model provides a better estimate of the vortex shedding patterns of the flow.
3.3. Comparison of Moments on the Blade
Using a medium-sized grid and the SST k-ω model, a number of simulations were performed to predict the integrated loads on the rotor through the operational wind speed interval of the wind turbine at 7 m/s, 10 m/s, 15 m/s, 20 m/s, and 25 m/s. Moments on the blade, namely, low speed shaft torque (LSST) and root flap bending moment (RFBM), are compared against UAE data over the operational wind speeds, as shown in Figures 8(a) and 8(b). Figure 8(a) reveals that the numerical study underpredicts the LSST by up to 20% compared with the UAE data, except at the 10 m/s stall speed, where the LSST is underpredicted by about 40%. In spite of the large deviation, numerical torque exhibits the same trend as the UAE data. The rather poor compatibility of the LSST data is attributed to the drag, that is, the dominant contribution to the LSST, especially at low wind speeds and for low CFD prediction capability for widely separated flows, as mentioned in Section 1. In contrast to LSST, RFBM comparison shows fairly good compatibility with the UAE data, as shown in Figure 8(b).
(a)
(b)
Streamlines over the blade and the velocity field at spanwise locations 0.3, 0.466, 0.633, 0.8, and 0.95 are given for 10 m/s in Figure 9, which shows that the flow is widely separated at mid-board locations. Separation over the blade affects the accuracy of the prediction of forces adversely, as mentioned above. Comparison of the pressure coefficients against experimental data at 10 m/s is shown in Figures 10(a)–10(e). It can be seen that the pressure coefficients are not as well predicted as in the 7 m/s test case. The effect of separation can be clearly seen at the suction surfaces 0.466, 0.633, and 0.8, where there are separation bubbles. In fact, at each spanwise location, prediction of the pressure coefficient at the suction side is less accurate than that obtained on the pressure side.
(a)
(b)
(c)
(d)
(e)
4. Concluding Remarks
Performance of the NREL Phase VI experimental wind turbine rotor has been simulated using the commercial CFD tool ANSYS FLUENT. The simulation conditions selected were an upwind turbine configuration with 0° yaw and 3°pitch alignment. For the simulations, a rotating reference frame model is used. A grid with 1.9 million elements was selected through a grid optimization study. For handling turbulence, the two-equation models available in ANSYS FLUENT were used for 7 m/s, and the k-ω SST model was found to be the most appropriate.
The resultant moments obtained through the operational wind speed range (7, 10, 15, 20, and 25 m/s) were compared against experimental data. The comparison of moments shows good agreement, especially for the RFBM. Although the LSST results deviate from the experimental data in the 20% range, the trend similar to that of the experimental data. Comparisons of pressure coefficients for stall speed (10 m/s) were presented for five different spanwise locations. The adverse effects of separation on prediction capability have been shown. However, as mentioned above, the capability of the CFD method to simulate highly separated flows remains poor. Moreover, the major deviation from the experimental data for the leading edge, caused by the large pressure gradients, is a problem that can be addressed by increasing the grid density in this region.
Although the current results can only be seen as preliminary, the flow field data and the velocity and pressure distributions obtained will serve as initial data for the multiphase analysis of air and water that is required for icing simulation. In a subsequent work, calculations at various wind speeds will be conducted on an iced wind turbine blade, and the performance of the blade will be compared with that of the clean blade presented in this study.
Symbols and Abbreviations
| : | Turbulent kinetic energy per mass (J/kg) |
| : | Atmospheric pressure (Pa) |
| : | Static pressure (Pa) |
| : | Radial location (m) |
| : | Radius of the wind turbine (m) |
| : | User-defined source terms for turbulent kinetic energy, dissipation rate, and specific dissipation rate, respectively |
| : | Mean rate-of-strain tensor () |
| : | Time (s) |
| : | Flow velocity component (m/s) |
| : | Reynolds stress tensor |
| : | Wind speed (m/s) |
| : | Turbulence dissipation rate () |
| : | Dynamic viscosity (Pa-s) |
| : | Turbulent viscosity (Pa-s) |
| : | Kinematic viscosity () |
| : | Density of the air (kg/m3) |
| : | Specific dissipation rate () |
| : | Angular speed (rad/s) |
| HAWT: | Horizontal axis wind turbine |
| LSST: | Low speed shaft torque |
| NREL: | National renewable energy laboratory |
| UAE: | Unsteady aerodynamics experiment. |
Acknowledgment
E. Sagol expresses her gratitude to WesNET for their financial support and to Scott Schreck of the NREL for the experimental data.