A Complete Procedure for Predicting and Improving the Performance of HAWT's

A complete procedure for predicting and improving the performance of the horizontal axis wind turbine (HAWT) has been developed. The first process is predicting the power extracted by the turbine and the derived rotor torque, which should be identical to that of the drive unit. The BEM method and a developed post-stall treatment for resolving stall-regulated HAWT is incorporated in the prediction. For that, a modified stall-regulated prediction model, which can predict the HAWT performance over the operating range of oncoming wind velocity, is derived from existing models. The model involves radius and chord, which has made it more general in applications for predicting the performance of different scales and rotor shapes of HAWTs. The second process is modifying the rotor shape by an optimization process, which can be applied to any existing HAWT, to improve its performance. A gradient- based optimization is used for adjusting the chord and twist angle distribution of the rotor blade to increase the extraction of the power while keeping the drive torque constant, thus the same drive unit can be kept. The final process is testing the modified turbine to predict its enhanced performance. The procedure is applied to NREL phase-VI 10kW as a baseline turbine. The study has proven the applicability of the developed model in predicting the performance of the baseline as well as the optimized turbine. In addition, the optimization method has shown that the power coefficient can be increased while keeping same design rotational speed.


Introduction
The potential of wind as an alternative energy to provide the world with power in line with a strategic plan to achieve the sustained resources development and environmental protection has made it target for many industrials and researches sectors. Nowadays, it is widely seen the tendency of the world to move toward sustainable energy approaches. Wind turbines are one of the main renewable energy technologies. However, there is a range of challenges associated with the rapid development of their technologies. Thus, the expansion of wind power requires development in three mainly fields: number, size and efficiency of wind turbines to balance the increasing demand of energy [1]. The increasing number and size have their own common and particular challenges, whereas the efficiency of wind turbines depends on different criteria such as rotor shape, drive unit and control type of the turbine. Different design methods were embedded with several optimization techniques, which enables multi-objective functions subjective to various constraints were used for improving the performance of the wind turbine under variable operating conditions [2,3,4,5]. The NREL phases for wind turbines are one of the most important targets for a vast number of studies. They are normally chosen because of their well-documented wind tunnel database that includes pressure distributions, separation boundary locations, polar characteristics and flow-visualization data [6]. The NREL phase -VI of 2-blades and 10.058 m rotor diameter stall regulated turbine which has the 20.95 % thick S809 airfoil shape [7], has been used for the present study. This turbine is originally designed with incorporation of Prandtl tip-loss and Corrigan and Schillings post-stall models to embed the three dimensionality effects [8,9], and the blade geometry was optimized for maximum gross annual energy production (GAEP) [8].
The present study aims to develop methods for optimization and performance predicting of HAWT's. MATLAB computer programs for both have been developed. The optimization part is based on modifying inputs design operating conditions and rotor shape of the TMASO program [10]. The performance prediction involves the BEM theory, Prandtl tip-loss and Glauert axial induction factor corrections, profile loss and a developed post-stall models of EOLO and Viterna. The EOLO model helps to predict performance under the rotational augmentation of the aerodynamics characteristics [11], whereas the Viterna model is particularly evident in high winds where aerodynamic stall occurs as the blade experiences high angles of attack [12].

Performance prediction
The analysis for prediction the performance of the stall-regulated NREL UAE phase-VI 10kW is mainly based on the BEM theory and it considers the post-stall region. A Matlab program is built to analyse any defined wind turbine blade shape. The information needed for the performance analysis are rotor radius R, number of blades B, number of sections, chord length distribution c (r) , tip pitch angle β o and twist angle distribution β t(r) . In addition to the aforementioned, the operation conditions, which include wind velocity range v 1 from 5m/s to 25m/s with an interval of 1m/s and the associated angular speed ω, which depend on the control strategy. The angular speed is set here as a linear function of wind velocity until v 1 reaches 7.2m/s, then it is set to be constant at 7.5rad/s, since the control strategy of the NREL 10kW turbine is stall-regulated.
In order to calculate forces and power, BEM states that induction factors have to be evaluated as they reveal the magnitude and the direction of the air flow, which is needed for finding the relative velocity. In the BEM theory a and a are not constant and dependent on λ r . In the prediction procedure, the induction factors are calculated by iteration with below equations: where σ is the local solidity (σ = Bc/2πr), φ is the relative angle between the relative velocity and the rotor plane of rotation, which is equal to the summation of the cumulative pitch β and the attack angle α, (φ = β + α). C L is the coefficient of lift, which is found along the blade sections while iteration for each α until the convergence. F is the Prandtl's tip loss correction factor, which is needed due to the pressure difference between the suction and the pressure sides of the blade that makes the air tend to flow around the tip, reducing lift force and hence power production near the tip, which is calculated as [13]: Because BEM theory does not yield reliable results for axial induction factors greater than 0.5, the Glauert's correction of calculation is used in this range.
where this equation is valid for a > 0.4, which corresponds to C F > 0.96. C F can be calculated as: Furthermore, as the value of a becomes greater than 0.2, the simple BEM theory will break down, thus we used the below correction value of a instead [14]: where, In order to predict torque and power with BEM, C L and C D need to be calculated. They rely mainly on the angle of attack which is in turn a function of the wind velocity, angular speed and the control strategy of the turbine. In addition, Re and the centrifugal pumping (or the 3-D radial flow) are important parameters that effect C L and C D . Re distribution for a rotational speed of 72rpm of the UAE phase-VI turbine is selected to be Re = 1.00 * 10 6 , which is the average value for all wind velocities from 5m/s to 25m/s along the whole blade, [15].
The prediction is done mainly by considering three different ranges of α. Since the control strategy of the turbine is a stall-regulated one, it is expected that the stall will appear then dominate the blade sections as the wind velocity reaches high values because the ω will remain constant at this range. In this study, v 1 is used rather than α to define the ranges, since it has only one value for each run, whereas α changes along the blade section. The first part of predicting C L and C D is based on curve fitting equations of Delft 2-D data up to stall angle of 9.32 degrees [12]. These data should predict results, which are in excellent agreement with measurements over the first part of the power curve from 5m/s to 8m/s for a tip pitch angle of 2,3 or 4 degrees. In the prediction procedure at the tip pitch angle of 3 degrees, 2-D Delft data provides better results by implementing the equations up to a wind velocity of 9m/s. The second part, after reaching the stall angle, a reasonable prediction agreement is achieved by implementing post stall UAE-derived airfoil data and Viterna equations over the wind velocity range of 9m/s to 25m/s. UAE-derived airfoil data of Cl and Cd were averaged, which are used up to the flat plate angle of 20 degrees, where the Viterna equations were inserted, which are guided by flat plate theory [12]. Experiments have proven that the flat plate equations agree with UAE derived data averaged along the blade at angles of attack of 20 degrees or higher. The poor prediction of the Viterna model in the velocity range between 9 − 20m/s is because of the averaged experimental data cannot help in prediction of post-stall region due to the 3-D complexities; such as producing a standing vortex in the on-board region as a result of interaction between the mean flow and the radial pumping one. Although, attention is paid when implementing Viterna equations after leading edge separation, it is taken into account that the equations are influenced by the aspect ratio of the blade, whose selection dictates the maximum drag coefficient and governs the predicted power at high wind speeds over 15m/s, it still giving poor prediction. Hence, the Viterna model is used for the velocity range of 20−25m/s, and better results are achieved by adjusting the constant of C Dmax while keeping the aspect ratio (AR) of 14.
C Dmax = 1.11 + constant * AR Airfoil stall angle and blade averaged values of C L and C D at α = 20 degrees are used as inputs for the Viterna equations, because these values C Lstall,average = 1.24 and C Dstall = 0.44 conform a glide ratio for the initial conditions, which agree with flat plate theory [12]. To bridge the gap from 2-D airfoil data to the Viterna method in the wind velocity range of 10 − 19m/s or rather in the angle of attack range from 9.21 • to nearly 20 • , the EOLO model is implemented in this post-stall region. This model is based on the flat plate equations and only valid in the transition zone between the "high lift, stall development regime, dynamic stall" and the flat plate, fully stalled regime [16]. This correlates to an angle of attack range of roughly 10 • to 20 • . For the fully stalled regime (20 ≤ α ≤ 45) the flat plate equations are: C Lmax and C Dmax are inputs for the flat plate theory equations and equal to C L,α = 45 • and C D,α = 90 • [16]. For the S809 airfoil, the value of 2.3 is assumed as maximum drag coefficient. In order to follow the flat plate theory, value of C Lmax at 45 • would be 1.24. Instead of assuming one value for C Lmax , the EOLO model uses an increased lift coefficient to take into account the radial flow along the blades that provides an increase in energy with the same flow [16]. After some analysis, a reasonable agreement is found by a nonlinear equation, which considers experimental data, After estimating C L and C D , then evaluating a and a for all range of v 1 , all required components for calculating the force, torque and power of the turbine with BEM are available. Tip loss is already considered by Prandtl's tip loss correction factor. Profile loss along the blade, which includes an appearing drag at higher wind velocities and angles of attack, is multiplied by the power [17]. Figure 1 shows the torque prediction in comparison with six further prediction models and the experimental data. Best agreement with experimental data is perceptible by the present performance analysis method, which is named "calc" in the figure. From 5m/s to 9m/s all prediction models differ only slightly. A discrepancy in torque prediction of nearly 1kN m at wind velocity of 10m/s to 19m/s is noticeable The model EOLO provides in this region among the six prediction models the smallest differences to experimental data, therefore it is used in the analysis. From 20m/s to 25m/s only the slope of the model of Viterna & Corrigan matches with experimental data. Hence, in the analysis, the trend of the Viterna model was shifted down by modifying the constant to be in the range of the experimental data. In figures 2 and 3 axial induction factor and angle of attack over wind velocity are shown, respectively. They are calculated by the post-stall model and aero-elastic code PHATAS and compared at three different spanwise locations. Aerodynamic modelling of PHATAS is based on the BEM theory as well [18]. Differences between the models in figure 2 are minimal. Shapes of the prediction models are mostly the same. Maximum deviation appears at 10m/s, which is a about 0.04 in prediction of the axial induction factor. This is due to the post-stall model used in the present study which starts at 10 m/s. In figure 3 α slopes of the prediction models differ little at a spanwise location of r/R = 0.3. At spanwise locations of r/R = 0.6 and r/R = 0.8 PHATAS model predicts angles of attack, which are up to 5 • higher. This could be justified because validation of axial induction factor is at each spanwise location only a little higher than the prediction of the model of PHATAS. Since α is dependent on ϕ, and ϕ is dependent on the axial and tangential induction factor and on λ (r) , a higher value of the axial induction factor leads to a lower value of ϕ. This reduction is increased at greater spanwise locations by the factor r. As a result, α decreases with the increase of both wind velocity and span.
In figure 4, angle of attack is plotted over r/R at different wind velocities. In addition to the validation, the performance-prediction method of LSWT and measuring method of local flow angles (LFA) are also shown. The LFA is measured with five-hole probes on stalks, whose length    Calculated angle of attack and prediction procedure by LSWT and measuring method of LFA of the NREL UAE phase-VI rotor [18].
is four fifths of the chord length and they are attached to the leading edge of the blade [19]. LFA results from the vector addition of v 1 and the local tip speed u. At low wind velocities prediction of α differ hardly between the calculated and LSWT. One reason for the increased difference of nearly 10 • at 25m/s between both models is that LSWT considers the induced effects of the blade configuration and those from the span-wise distribution of trailing vorticity, which is distinctive at higher wind velocities [12]. The much larger values of α of the LFA method can be explained due to a direct use of the local flow angles without converting to the angle of attack because of uncertainty over how the flow angle changes between the position of the five-hole probe and the leading edge of the airfoil [19].

Shape optimization of stall-regulated WT
The torque matched aerodynamic shape optimization approach (TMASO) will be used [10]. The same approach of changing the rotor shape to reach maximum power coefficient while keeping the torque matching with the drive unit is kept. Consequently, the rotational speed is kept constant at 71.6rpm. The number of blades 2, sectional airfoil S809 and blade radius 5 are kept the same. One initial angular speed, one initial wind velocity and design angle of attack need to be known as inputs for the TMASO. As a result of a constant rotational speed at 71.6rpm, the initial angular speed is 7.5rad/s (ω = 2πn/60). Initial designed wind velocity is calculated by (v 1,o = ω o R/λ o ). The value of λ o is taken as λ o = 5.2, which gives the maximum power coefficient C P max of the NREL UAE phase-VI rotor at β o = 3 • and the design wind velocity v 1,o = 7.2m/s. Modification of the design angle of attack α d is suggested, thus, a range of angles within a range of 3 • to 9 • were tested. C P max and maximum glide ratio GR max are found at angle of attack of 6 • . The initial design of the blade, which includes chord length and twist angle distributions are calculated by Schmitz method, is used as input. Airfoil characteristics C L , C D and GR for each iteration of α and Reynolds number are found simultaneously by XFOIL. The objective is to find the maximum power coefficient, which is hence, the objective function is formulated as: As the rotor coupled to the drive, T rotor should match the torque of the drive unit T drive (n). The difference between them should be minimized so that torque matching can be achieved. For this purpose, the following constraint is used during the optimization.
T drive (n) is calculated by BEM before and during rated operating. The upper and lower bounds of ω are limited to ensure keeping the optimized rotational speed in the range of the initial one (generally wind turbine generators are set to operate at a defined rotational speed, which represent the capability limit that cannot be exceeded, and it is correlated either directly or via gear-box to the grid frequency). Bounds of the velocity are more stretched. The procedure optimizes the wind speed and rotational speed together with the independent rotor blade shape decision parameters. The optimization is done by using the gradient-based function f mincon with the active-set algorithm in MATLAB. f mincon can search for the minimum of a constrained non-linear multi-variable function. A new velocity is found, which indicates a new operation condition at new value of λ and hence new values of P w , P and T rotor , which should be equal to T drive . The optimized rotor shape parameters β (r) and c (r) are calculated by using Schmitz method. Achievement of TMASO is an increment of the power coefficient, torque and power. The optimized rotational speed is kept exactly the same as the initial one, whereas the optimized tip speed ratio increases from 5.2 to 5.45 and wind velocity decreases from 7.2 to 6.85m/s at the design angle of attack of 6 • . A special treatment is considered for the chord length, where the improvement in performance always tends towards increasing chore values. Hence, the maximum possible of the chord is always expected to be provided by the optimization procedure. Therefore, an upper limit for the chord is done to have a realistic optimal design. The chord length and twist angle along the blade are shown in figure 5.

Modified Post-stall model
As indicated above, the Viterna model with some modifications is performing better in the higher wind velocity range v 1 ≥ 20, whereas the EOLO model can predict the performance within the stall region 9 < v 1 < 20. They both serve to predict the performance of the baseline turbine of NREL phase-VI 10kW. In order to make EOLO perform for different scales and shapes of turbines it has to be modified. A semi-empirical equation,which is derived from EOLO and the existing experimental data with additional observations is developed in the present study to overcome the performance prediction of any existing wind turbine. All parameters defining the turbine rotor, such as rotor radius, chord, the number of blades are included in the developed model.
where, the suffix i refers to reference wind turbine parameters, which are the radius R i = 5 m and the mean chord c i = 0.5606. These constants used to calculate the maximum lift coefficient as, then, use it in flat plate equations 8 to get C L and C D . Figures 6,7,8 below reveal the advantage of using the model over the equation (the model refer to the modified model and the equation is the equation used to estimate EOLO, equation 8) for the performance prediction of different turbine characteristics. The main criteria used to measure the advantage are limit values of the post-stall region and trends of the curves. Limits values are simply characterized by the starting value of the stall at v 1 < 10, which is trusted since it is found by 2D Delft equations and the end value at v 1 = 20, which is calculated by the Viterna model, which is trusted for wind velocity at that range since the airfoil will behave like a flat plate as explained above. The second criterion is the trend of the curve, where the standard trend is taken from the available experimental data, which shows a continued decrement of the torque after v 1 = 10 till reaching the end value of v 1 = 20. In addition, the smoothness of the curve can be another indication for the feasibility of the model. It is required to mention here that B is not set in the model since it can be predicted rather well by both the model and the equation. Nevertheless, the model still shows better trends especially for the torque in figure 6. Figure 7 shows the torque versus the wind velocity of the model in comparison with the equation for a different chord ratios, which is the ratio of any wind turbine mean chord c by the baseline turbine mean chord (c i ), which is here c i = 0.5606. It is found that the model is applicable until c/c i = 1.75, and it can predict power and torque in much better way than the equation. Hence it is included in the model, equation 17. Figure 8 shows the effect of the rotor radius on the torque for both model and equation. The model is applicable until radius value of 40 m, and for all of the radius ranges, it performs better than the equation. Consequently, the torque increases with the increment of the rotor radius, and hence the wind power increases too (for the new turbine area), and the angular velocity is kept constant, thus the power coefficient decreases. Since, the model is set based on the information obtained from the baseline turbine; for instance, in equation 18 it is set in terms of wind velocity rather than angle of attack, the relation between α and v 1 has to be checked for better prediction at different radius values. As the radius increases, the tangential velocity increases, and hence the angle of attack decreases. Figure 9 proves that fact, where it reveals that within the tested radius range there is a small change (reduction) of the mean angle of attack as the wind velocity increases, hence, the stall will delay as the rotor radius increases. Thus, the range of wind velocity defining the post-stall range has to be shifted as the radius increases. It is depicted in figure 9 with dotted lines, for the same mean angle of attack where the stall starts about α mean = 10, which is associated with wind velocities v 1 = 9.5, 10, 10.5, 11 of the radius R = 5, 10, 20, 40, respectively. The stall delay with the radius increase is depicted in figure 10, where the locus of the different radii torques form an exponential curve.   Figure 10. Effect of radius on the torque over wind velocities.

Performance prediction of the stall-regulated shape optimized WT
In the performance prediction of the shape optimized turbine, the original chord length and the twist angle distributions of the baseline turbine are replaced by equations from figure 5 of the optimized blade shape. The BEM method with the modified stall model are used in the analysis to resolve the increment of the chord length in the shape optimized turbine. In figure 11, the power prediction of the shape optimized turbine is higher than the experimental data of the baseline turbine and the velocity is decreases from v 1,o = 7.2m/s to the shape optimized turbine velocity of v 1 = 6.85m/s. In the post-stall region 9 > v 1 > 20m/s, the rated power is increased from P N,o = 10 to P N = 14.5kW . This is due to the increase in the torque from T N,o = 1.33[N.m] to T N = 2[N.m], while the angular speed is kept constant at ω o = 7.5rad/s. Hence, the post-stall model that includes the increase of the chord is performing well as the mean chord increment is limited to a value less than 2 of the baseline turbine mean chord.  Figure 12. Power coefficient increment of the shape optimized turbine as compared to the baseline turbine experimental data [16].
In figure 12, the calculated power coefficient is plotted over tip speed ratio. It is recognizable in comparison with the experimental data of the baseline turbine that C P is raised for the optimized blade shape. At the baseline turbine design tip speed ratio of λ o = 5.2 the C P is increased from about 0.35 to 0.54 at the new optimized tip speed ratio of λ = 5.45. As a result, the objective to increase the power coefficient by an optimized blade shape is achieved.

Conclusions
The developed prediction method of the stall-regulated wind turbine agrees well with the experimental data, hence, it is possible to rely on it for investigations. However, for a general performance prediction of stall-regulated wind turbines, it needs a method for predicting the lift and drag coefficients in addition to consideration of turbine characteristics and operational conditions The developed performance prediction method, which includes the change in radius and chord length, can predict the performance for a relatively high range of wind turbines scales, number of blades and chord lengths. The post-stall model parameters have to be adjusted before applying to any wind turbine. The gradient-based optimization method TMASO performs well for a large scales wind turbines. The limitation of the chord length is necessary to keep it from a rapid increase when optimizing of the rotor shape. In general, the developed methods provide a complete set of tools for performance prediction and shape optimization of different wind turbines operated with stall-regulated control strategy. In general, these developed methods serve in changing the rotor shape while keeping its angular speed constant