Abstract
Development of wind energy is gaining a great interest nowadays with over 25% annual growth in any configurations and sizes. Designing a wind turbine with improved performance remains the ultimate research and development goal. Therefore, understanding the influence of the different wind turbine key parameters, e.g., tip speed ratio (TSR), twist and pitch angles, number of blades, and airfoil chord distribution, is critical. In this work, an improved blade element model (BEM) is developed by correcting the code with tip loss, Buhl empirical correction, skewed wake, and rotational effect. These corrections are necessary to extend the application of the method to the turbulent wake regime for the horizontal axis wind turbine (HAWT) configuration. Results from the developed code are compared with the NREL measured data of the twobladed Unsteady Aerodynamics Experiment (UAE) phaseVI turbine. They indicate an improved trend with the incorporation of these corrections. The performance of the threebladed 3.5kW HAWT was tested and found to follow closely the experimental trend. As the mismatch between these low fidelity analyses (BEM) and the experimental work persists, the undertaken analysis demonstrates its limitation, and it emphasizes the role of high fidelity wind tunnel and flow simulations. BEM analysis, however, shows that designing blades with proper twist and pitch angles, and targeting suitable TSR could lead to substantial gain (over 10%) in the performance.
Keywords:
BEM; TSR; Twist angle; Pitch angle; HAWT; Power coefficientBackground
Currently, significant number of approaches, theories, and models are developed to predict wind turbine performance. Although high fidelity CFD continues to be perused [13], it remains extremely costly as far as computational resources, analysis time, and required expertise. Among lower fidelity approaches, methods, and models, blade element momentum (BEM) method is the oldest and remains to be the most widely used method for predicting wind turbine performance. It was originally developed by Glauert [4] who combined blade element theory and momentum theory to analyze the airplane propeller performance. Blade element theory assumes that blades can be divided into multiple elements, which can act independently as twodimensional airfoils. The forces and moments can be calculated separately then summed to obtain the overall blade forces and moments. The other half of BEM method, the momentum theory, assumes that wind turbine harvests energy from incoming flow; thereby, the flow is subjected to pressure and moment losses. Using momentum theory, induced velocities from the momentum loss can be calculated. These induced velocities can affect the flow over blades and the forces on them. By combining the two theories and setting up an iterative process, forces and momentum on blades can be obtained.
In practice, it is those corrections that make BEM method applicable in wind turbine design. The corrections include tip loss correction [5], hub loss correction [6], Glauert [7] and Buhl empirical corrections [8], skewed wake correction [9], and ‘3D correction’ of Snel et al. [10]. Tip loss correction accounts for the influence made by vortex shedding from the blade tips into the wake. Hub loss correction accounts for the vortex shed near the hub. Glauert and Buhl empirical correction accounts for the turbulent wake phenomenon. Skewed wake correction accounts for the influence made by wind turbine when operating at yaw angles, while 3D correction of Snel et al. accounts for the lift augmentation caused by rotation.
In this work, an improved BEM method code is redeveloped that incorporates the above corrections. The developed code is validated against the NREL measurements of the UAE phaseVI turbine. The results clearly demonstrate the plausible trend of the calculation with the incorporation of the corrections. Based on the computation code, the performance of the small scale 3.5kW Horizontal Axis Wind Turbine (HAWT) was assessed along with the study of some pronounced parameters. Subsequently, the recommendations for 3.5kW HAWT design and operation condition were derived.
Methods
Mathematical formulation
The mathematical code in this work for wind turbine design is based on BEM method. It equates the losses of the air flow momentum to the load on the blades. By applying axial and angular momentum conservation s, the incremental axial/thrust force dT and angular torque dQ on blade sectors can be obtained, as illustrated below according to Manwell [11]:
where ρ is the air density, U is the mean air flow velocity, r is the local blade radius, ω is the blade angular rotational speed, a and a′ are the axial and the tangential induction factors, respectively.
According to the blade element theory [11] and by referring to Figure 1, one also has:
where σ′ is the local solidity which can be calculated by σ′ = Bc/2πr; here, B is the number of blades and c is the local blade chord length. φ is the angle of relative wind, C_{l} and C_{d} are the lift and drag coefficients, respectively.
Figure 1. Illustration of the incoming flow and angle of attack.
Thus, the total torque at the shaft (Q) is the summation of dQ for all the blades elements, and the wind turbine power is given by P = Q * ω. Eventually, the wind turbine power prediction lies in solving the axial induction factor a and the tangential induction factor a′ [11]. By combining Equations 1, 2, 3, and 4, one can solve the induction factors a and a′ as given below:
where λ_{s} is the local speed ratio which is defined as λ_{s} = ωr/U. The axial thrust coefficient C_{T} is expressed as:
By substituting the two induction factors into Equation 2, the wind turbine performance can be assessed.
Corrections for BEM method
When the BEM method was originally developed by Glauert, it was far from giving accurate and reliable results before the implementation of some important corrections. These corrections include tip and hub losses, Glauert and Buhl empirical corrections, as well as the skewed wake correction in addition to the 3D rotational correction.
Tip and hub losses
Tip and hub losses consider the influence of vortices shed from the tip and hub, which play an important role in the induced velocity distribution at the rotor. Prandtl model is the most widely used model for calculating the tip and hub loss correction factors. An approximate formula for the Prandtl tip and hub function was firstly introduced by Glauert [4] as expressed below:
Where is for tip loss factor, and is for hub loss factor, where R is the radius of blade tip and R_{hub} is the radius of blade hub. Ultimately, the correction factor can be calculated by the multiple of the tip and hub loss correction factors, as illustrated in Equation 9:
Figure 2 shows the relationship between the thrust coefficient and axial induction factor at F ≤ 1 and with and without tip and hub correction factors. A pronounced increase in the axial thrust, observed as F, increases.
Figure 2. Classical a C_{T}curve under different F values.
Due to the influence of the tip and hub correction factors, the forces derived from the momentum theory have to be modified accordingly:
Thus, the calculations for two induction factors and thrust coefficients are updated as follows:
Glauert and Buhl empirical corrections
As the induction factor is greater than 0.4, wind turbines will be under a turbulence wake. This puts an upper limit for the validity of the basic theory. Glauert developed a correction [12] to the rotor thrust coefficient based on experimental measurements of helicopter rotors with large induced velocities. Buhl [8] later developed a new relation between rotor thrust coefficient and induction factor which solved the instability caused by applying Glauert correction as illustrated below:
Figure 3 shows the Classical relationship between thrust coefficient and axial induction factor, Glauert correction relationship, Buhl empirical correction relationship, and experimental data of Moriarty et al. [6].
Figure 3. Axial thrust coefficient plot under different models at different values of F.
The widespread experimental data indicates that the thrust coefficient is not a simple function of axial induction factor in the turbulent wake region. As clearly shown from Figure 3, the Buhl empirical relationship between C_{T} and a plausibly follows the experimental data trend and could also solve the instability/mismatch problem caused by Glauert's empirical relationship at the same time.
Skewed wake correction
Although BEM method was originally proposed to solve for the axisymmetric flow, wind turbines are often running at yaw angles. This again invalidates the basic theory unless a correction is used accounting for the skewed effect. Snel and Schepers [13] derived the following correction formulation:
where χ is the wake skew angle, and ψ is the azimuth angle. However, Eggers's work [14] concluded that the correction can be large, emphasizing its inclusion.
Rotational effect
The effect of rotation was first investigated intensively for helicopter rotor. Later, the fact that aerodynamic power tends to exceed the design value for wind turbine starts attracting more attention and translated to different mechanisms. These include centrifugal pumping effect, stall delay, rotational augmentation, etc. [15]. Despite of these developments, a census approach is still lacking, particularly for the effect on the drag. Nevertheless, the 3D correction of Snel et al. received amble attention and is incorporated in present work. It provides an increase of the aerodynamic lift coefficient for the effects of rotation, which is described below:
where the ‘potential lift coefficient’ C_{l,non−pot} is defined as 2πsin(∝ − ∝_{0}), ∝_{0} is the angle of attack when lift coefficient is zero.
Summarizing, the above mathematical formulation proposes an iterative procedure for solving the axial and tangential induction factors accounting for different types of losses. As obtaining the two induction factors is indispensable for predicting the performance of wind turbines, Figure 4 illustrates the flowchart for computing these induction factors.
Figure 4. Calculation flowchart for induction factors.
The tolerance (ε) in the iterative procedures for the induction is set to be 0.001. The procedures are integrated to classical BEM method calculation, leading to the development of robust numerical procedures, i.e., extended BEM. The code follows the iterative procedure shown in Figure 4. It provides the wind turbine power generation/coefficient as an ultimate result.
Verification of the extended BEM model
The UAE phaseVI turbine [16] is used as a baseline to validate the proposed extended BEM model. Its experimental data were collected by National Renewable Energy Laboratory (NREL) in the 24 × 36m wind tunnel of NASA Ames. UAE phaseVI turbine is a twobladed and S809 airfoil type wind turbine with fixed operation state of 72 RMP rotational speed and a threedegree pitch angle. The blades has a literal sweeping chord length of 0.737 m which decreases to 0.305 m and a twist angle of 20.04° at the largest chord [15] that fades to −2.50° at the tip, as depicted in Figure 5.
Figure 5. Chord length and twist angle over blade.
The user of the code needs to provide the following necessary data: (1) the blade pitch angle, (2) sectional twist angle, (3) chord length at different radii, and (4) the lift and drag coefficients as a function of angle of attack. The aerodynamic data for S809 airfoil were taken from [15,17]. The airfoil aerodynamic data, as shown in Figure 6, were initially developed at TUDelft wind tunnel and subjected to the angleofattack shift technique for large angles of attack and empirical tool ‘stall coefficients (StC)’ for deep stall [15].
The influence of different corrections on the power coefficients of the turbine under different wind speeds were sought and compared with experimental data, as shown in Figure 7.
Figure 7. Power coefficients under different models.
It is clear from Figure 7 that the power coefficients obtained from BEM without any corrections overestimate the performance, while with tip and hub corrections, it reduces the performance at large margin. Simultaneously, the model with tip loss, Buhl empirical correction, and 3D correction of Snel et al. managed to predict the power coefficients under different wind speeds. However, the discrepancies between the results obtained from the model with tip loss and Buhl empirical correction and the experimental data can still be observed. As shown in Figure 7, the results obtained from the model strictly follow the experimental data with slight underestimation. As this mismatch becomes large in other cases, it can be attributed to two reasons: (1) Empirical correction of the tip loss and Buhl require further improvements in order to describe wind turbine and incoming wind interaction and (2) lack of accuracy in the airfoil aerodynamic data that accounts for the rotation effect that will be inherited in the 3D correction.
Application of BEM on 3.5kW HAWT
Windspot is a small size HAWT, 3.5kW threebladed turbine with a 4.05m rotor diameter. It is characterized with low cutin speed (<3m/s) and a rated power at 12m/s. It can be fitted with different set of blades, essentially at zero twist angle with a centrifugal/active pitch control system. The chord length is 0.254 m at the root, decreases linearly to 0.156 m at the blade tip.
Aerodynamic data
The Windspot airfoil was extracted from the solid model in the form of native IGES geometry. Panel method (governed by potential flow) and high fidelity CFD numerical simulations are developed and carried out to obtain the BEM required aerodynamic lift (C_{L}) and drag (C_{D}) coefficients at varying angle of attack. Figure 8 shows the geometry of the airfoil, in which the captured nodes divide the airfoil into several panels. The trailing and leading edges are further refined with additional panels.
Figure 8. Airfoil shape of Windspot.
Figure 9 shows the C_{L} result of the panel method which clearly varies linearly with the angle of attack. The panel solution is based on potential flow theory with no reference to the viscous boundary layer. These results, and when incorporated zero viscous drag, clearly will overestimate wind turbine performance.
Figure 9. Results of the C_{L}from the panel method.
High fidelity CFD simulation is also carried out for the Windspot airfoil. The flow is governed by the two dimensional, steady incompressible Navier–Stokes equations. These equations represent statements of mass and momentum conservation and are written after applying the scalar variable expansion and ensemble (over bar) averaging as:
where t is the time advancement, x_{i} is the Cartesian coordinate (i = 1, 2), u_{i} is the velocity component in x_{i} direction, ρ is the density, g_{i} is the gravitational acceleration component in x_{i} direction, τ_{ij} is the stress tensor components, p is the pressure, μ is the molecular viscosity, u_{i}^{’} is the velocity fluctuations about ensemble average velocities. The term is the Reynolds stresses and is modeled utilizing the mean (ū) velocity via the common eddy viscosity kepsilon turbulent. It is expressed as:
where k is the kinetic energy and μ_{t} is the turbulent viscosity which links k and turbulent dissipation rate such that , where f_{μ} and C_{μ} are empirical constants. Substituting Equations 21 to 20 conveniently allows summing the Reynolds stresses terms to the diffusion term (2nd right hand term in Equation 20) with an equivalent viscosity μ_{equ} = μ + μ_{t}. Therefore, closure of the above system is achieved with the integration of two additional transport equations for k and ε. The above equations are discretized on wellposed/bounded domain, forming an algebraic system of equations. The boundary conditions including incoming velocity (Dirichlet), side boundaries symmetry (Neumann), outflow at the exit, and nopenetration, noslip at the airfoil surface.
The airfoil geometry and the flow zone model were built and discretized using quadtype cells. The airfoil surface is padded with a refined cell cluster of less than 0.01% chord length initial height and at smooth growth of 1.4 to achieve one unit normalized wall distance (y+ = 1). The geometry, the baseline mesh, and the imposed boundary conditions are shown in Figure 10. The domain was subjected to the upstream uniform velocity (5 m/s which is the annual average wind speed in Masdar City [18]) and the output atmospheric pressure to be used during the iteration in the event of recirculation at the outflow boundary.
Figure 10. The meshed airfoil geometry.
Results and discussion
The Fluent iteration, double precision, and pressurebased SIMPLE solver were used to compute the aerodynamic data as a function of angle of attack with as low as 10^{−10} residuals in the continuity and the two momentum equations. Except for the pressure, a second order discretization is implemented. Mesh sensitivity studies were carried out first to assess solution dependency. Results of the 5° angle of attack case is carried out at the baseline, fine and two levels of coarseness, and these are summarized in Table 1.
Table 1. Mesh sensitivity, lift, and drag coefficients and their relative errors, and maximum y+
The mesh dependency of the solution was observed to reduce as it became finer. The baseline mesh appears to capture the solution at a reasonable level of accuracy relative to the other successive refinement and coarseness. This is also clear from the normalized value of the wall (y+) which is desired to be one unit to avoid excessive refinement. Figure 11 shows the results of the lift and drag coefficients as functions of the angle of attack for the 3.5kW Windspot airfoil from baseline mesh. The data is fitted with an appropriate polynomial for convenient integration in the BEM code.
Figure 11. CFD results of the lift and drag coefficients as functions of angle of attack.
The improved BEM is used to estimate the Windspot power coefficient, which is limited by Betz Law (<0.59) [19]. As the operation of the wind turbine is unknown, it is reasonably assumed to operate at a fixed rotational speed of 12 rad/s and 7° pitch angle. The power coefficients under different wind speeds, yet at that fixed rotational speed, are obtained and compared to experimental data [20], as depicted in Figure 12. Results of the simulation show appreciable deviation from those obtained experimentally while continue to exhibit a plausible trend.
Figure 12. Power coefficients under different wind speeds.
The influence of the pitch angle, twist angle, and TSR is also investigated. Figure 13 shows the effect of different pitch angles at a zero twist at the previously considered rotational speed (12 rad/s). A pitch angle in the range of {3,7} appears to be more optimal for a flow speed up to 8 m/s. Figure 14 depicts the influence of the twist and clearly shows that a 10° twist can lead to a measurable increase in the power coefficient which translates to a measurable increase (up to 4°) in the overall angle of attack. Higher twist resulted in near stall condition and drastic fall in the power coefficients. The influence of the tip speed ratio is also investigated at the optimal pitch and twist angles to observe their combined effect on the power coefficient as shown in Figure 15.The tip speed is varying here by the rotational speed. The results suggest that a maximum power can be obtained near TSR of 4.
Figure 13. Power coefficients under different pitch angles.
Figure 14. Power coefficients under different blade twist angles.
Figure 15. Power coefficients under different TSRs.
The wind turbine is under different operational conditions which are listed in Table 2. It is clear from Figures 13 and 14 that choosing a proper pitch angle and designing a blade with proper twist angle are important for wind turbines. Excessive choice of pitch angle combined with twist would also lower the harvested wind turbine energy. Figure 15 shows that TSR as another important parameter for wind turbine, and higher power coefficients is achieved by operating at a suitable TSR value.
Conclusion
An improved BEM code for HAWT is developed. It is incorporating four important corrections: (1) tip loss, (2) Buhl empirical correction, and (3) skewed wake correction, as well as (4) 3D correction of Snel et al. These corrections extend the application of the BEM method to the turbulent wake regime of the HAWT. The developed code was validated by comparing the results from the code against those of NREL measured data of the UAE phaseVI turbine. The accuracy of the BEM inherently depends on obtaining an accurate aerodynamic data, i.e., lift and drag coefficient as a function of angle of attack. Therefore, potential theory (via Panel method) cannot sufficiently provide such data. Furthermore, as wind turbine airfoil shapes depart from classical NACA series, high fidelity CFD analysis for turbine blade was carried out. Results of the calculation that were based on the extended BEM show a plausibly match with the experimental trend of the Windspot turbine. Furthermore, detailed performance assessment was carried out by considering the effect of the twist angle, pitch angle, and TSR. The code illustrated these effects and results show designing blades with proper twist angle that harmonized with suitable TSR and pitch angle can definitely lead to substantial gain in performance of over 10%. In the case of 3.5kW Windspot, revising the blade with 0 to 10degree twist angle and operating the wind turbine under 7degree pitch angle and TSR around 4 can lead to significant increase in power generation.
Competing interest
The authors declare that they have no competing interests.
Authors' contributions
Under IJ's supervision and as his masters degree student, SL studied and coded the BEM. Together they carried out the validation of the code and conducted the sensitivity study on 3.5kW Windspot. All authors read and approved the final manuscript.
Authors' information
IJ is currently an associate professor at Masdar Institute of Science and Technology. Before, he worked as a visiting ME professor at MIT. He also possesses extensive industrial experience of over 12 years in aerospace and automotive from Michelin and GE. SL is currently a doctorate degree student.
Acknowledgments
The financial support was received from the Masdar Institute of Science & Technology. The generous sponsorship and data support of Sonkyo Energy team (Iñigo Portillo, Christophe Lopez, Javier Vidal, and José Luis) are also highly acknowledged.
References

Janajreh, I, Qudaih, R, Talab, I, Ghenai, C: Aerodynamic flow simulation of wind turbine: downwind versus upwind configuration. Energy Convers and Mana. 51, 1656–1663 (2010). Publisher Full Text

Janajreh, I, Talab, I, Macpherson, J: Numerical simulation of tower rotor interaction for downwind wind turbine. Model Simul Eng (2010). Publisher Full Text

Hu, D, Hua, O, Du, ZH: A study on stalldelay for horizontal axis wind turbine. Renew. Energy. 31, 821–836 (2006). Publisher Full Text

Glauert, H: Airplane propellers. In: Durand WF (ed.) Aerodynamic Theory, pp. 169–360. Springer, Berlin (1935)

de Vries, O: Fluid Dynamic Aspects of Wind Energy Conversion, DTIC, VA (1979)

Moriarty, PJ, Hansen, AC: Aerodyn Theory Manual: National Renewable Energy Laboratory/TP50036881, NREL, CO (2005)

Glauert, H: A General Theory of the Autogyro. Report and Memorandum No. 1111, British ARC, UK (1926)

Buhl, ML: N R E Laboratory: A New Empirical Relationship Between Thrust Coefficient And Induction Factor for the Turbulent Windmill State, National Renewable Energy Laboratory, CO (2005)

Pitt, DM, Peters, DA: Theoretical prediction of dynamicinflow derivatives. Vertica. 5, 21–34 (1981)

Snel, H, Houwink, R, Bosschers, J: Sectional Prediction of Lift Coefficients on Rotating Wind Turbine Blades in Stall, Netherlands Energy Research Foundation, Petten (1994)

Manwell, JF, McGowan, JG, Rogers, AL: Wind Energy Explained: Theory, Design and Application, Wiley, NJ (2010)

Glauert, H: The Analysis of Experimental Results in the Windmill Brake and Vortex Ring States of an Airscrew, HMSO, London (1926)

Snel, H, Schepers, J: Joint investigation of dynamic inflow effects and implementation of an engineering method, ECN, Petten (1995)

Eggers, A Jr., Chaney, K, Holley, WE, Ashley, H, Green, HJ, Sencenbaugh, J: Modeling of yawing and furling behavior of small wind turbines. 2000 ASME Wind Energy Symposium, 19th, AIAA, Aerospace Sciences Meeting and Exhibit, 38th, pp. 1–11. Reno, NV, USA (2000)

Lindenburg, C: Investigation into Rotor Blade Aerodynamics, ECN, Peten (2003)

Simms, DA, Hand, MM, Fingersh, LJ, Jager, DW, Cotrell, JR, Schreck, S, Larwood, SM: Unsteady Aerodynamics Experiment Phase VI: Wind Tunnel Test Configurations and Available Data Campaigns, National Renewable Energy Laboratory, CO (2001)

Musial, WD, Butterfield, CP, Jenks, MD: A comparison of twoand threedimensional S809 airfoil properties for rough and smooth HAWT rotor operation, Solar Energy Research Inst, USA (1990)

Janajreh, I, Talab, I: Wind data collection and analyses at Masdar City for wind turbine assessment. Int J Ther Envi Eng. 1, 43–50 (2010). Publisher Full Text

Betz, A: Das maximum der theoretisch möglichen ausnützung des windes durch windmotoren. Zeitschrift für das gesamte Turbinenwesen. 26, 307–309 (1920). PubMed Abstract

Site Expérimental pour le Petit Eolien de Narbonne: Rapport de Test n° 19 version 2 du 27 Avril 2010 Eolienne Sonkyo Energy Windspot 3.5 kW (2010). (2010) http://www.sepenmontplaisir.fr/index.php?option=com_docman&task=doc_download&gid=29&Itemid=27 webcite. Accessed 2010