Multi‐site aerodynamic optimization of wind turbine blades for maximum annual energy production in East Iran

Aerodynamic optimization of a horizontal axis wind turbine blade is performed for maximizing total annual energy production in three different sites located in the east part of Iran. While most researches have focused on the optimization of wind turbine blades for a specific location or air velocity, the main objective in practical problems is to find the best single design for multiple locations giving the highest energy output. The current paper addresses this issue and shows that multi‐site optimization of wind turbine blades is the key to maximizing energy extraction at multiple sites while reducing the manufacturing costs by proposing one optimum blade for all sites. In order to achieve this goal, the known Riso wind turbine blade is used as the base geometry and is optimized using a modified blade element momentum (BEM) theory and genetic algorithm. The objective function is maximizing sum of energy production at the selected sites. Blade twist and chord length distributions are chosen as design parameters. By running various optimization cases, it is shown that the multi‐site optimum blade has the highest sum of annual energy production by 4.45% increase compared with the Riso wind turbine blade.


| INTRODUCTION
Increased demand for energy is one of the main concerns for countries in recent years. Finite sources of fossil fuels, along with environmental problems, are becoming more expensive and that is why renewable energies attracted more attention than ever before. However, there is still a long way to go toward harvesting more from them, as in 2018, only about 4% of the world's energy consumption was supplied from modern renewables used to generate electricity, including wind, geothermal, solar, biomass, and waste. Among those renewables, the wind contributed more than 51%. Although the trends show an increase in renewables and mainly wind energy capacity for the future, being costly is one of the main barriers toward the global extension of its usage, which reduces competitiveness with fossil fuels. In an example, although high wind potential exists in Iran, still more than 98.5 percent of electricity was supplied via fossil fuels in 2018. 1 One practice in keeping wind energy costs down is to design site-specific wind turbines, whose rotor blades are optimized for specific local wind conditions that extract the highest possible amount of kinetic energy from wind. However, this is in contrast with the current trend in the wind turbine industry as manufacturers prefer to design and construct a set of standard wind turbines, and customers have to select the most suitable one among those. 2 Obviously, the selected WT may not be the most efficient one for that specific site; in other words, manufacturers sacrifice performance for lower design and construction costs.
This paper studies the possibility and efficiency of multi-site-specific optimization of HAWT blades in the east part of Iran. Design and manufacturing costs, site preparation and installation cost, maintenance cost, and financing cost form the overall cost of a WT. 3 On the other hand, energy production is the main output of this system. Therefore, multi-site optimization contributes to the efficiency increment firstly through reducing costs by making similar blades for multiple sites and therefore reducing the manufacturing costs, and secondly, by maximizing the energy extraction.
Iran is the 17th largest country in the world, with an area of 1,648,000 [km 2 ], and significantly varying climatic conditions across the country. On the other hand, the wind energy industry is just at the beginning of its way, and more importantly, when it comes to the manufacturing of technologically sophisticated parts like the blades, Iran is facing technical challenges. 4 The highly varying climate across the country and challenges in the wind turbine manufacturing process will justify the importance of introducing a single wind turbine that works satisfactorily across multiple sites.
The main attempts on WT blade optimization can be categorized into single-objective and multi-objective optimizations. The single-objective optimization has been the topic of various researches, for either maximum AEP, 5,6 maximum power, 7,8 minimum cost of energy, 9 and maximum output torque. 10 In other works, the multi-objective optimization is used, some of which take maximum AEP and minimum blade loads (thrust and flap-wise root moment), 11 maximum AEP, minimum blade mass and flap-wise root moment, 12 maximum power and minimum blade tip displacement, 13 maximum AEP, minimum blade cost and mass, 14 maximum power and minimum noise, 15 maximum AEP, and minimum blade mass. 16 More wind turbine optimization methodologies can also be found in. 17 While most researches have focused on the optimization of wind turbine blades for a specific location or air velocity, the question that arises is: what would happen if the optimum blade is located on other sites with different wind conditions? In the current work, the possibility of designing an optimum blade for maximizing total annual energy production in multiple sites is investigated. Initially, three different sites in the eastern part of Iran are selected with different wind speeds during the year, that is, Binalood, Khaf, and Zabol. Then, separate optimum blades are designed for each site (named as single-site optimums). Next, double-site optimums are introduced, and, finally, triple-site optimization is carried out trying to maximize the sum of AEP of all sites, and the results from all optimization cases are compared against each other. All mentioned optimization cases are initially done for blade twist angle distribution, and then both twist angle and chord length distributions are optimized; by comparing those, readers will gain insight into the effect of each optimization parameter in overall blade performance. Some improvements are also made to the BEM code, to make its validation against field experiments of Riso WT more accurate, especially at higher wind speeds.

| MATHEMATICAL MODELING
The required mathematical description for aerodynamic calculations, along with the Riso WT characteristics, are described, including BEM theory, its corrections and validation, Riso WT specifications, Weibull distribution for wind modeling, and AEP calculation.

| Blade element momentum theory
Blade element momentum theory is a well-known method for performance analysis of wind turbines under specified geometrical and meteorological conditions and is documented in numerous references, like 18 or 19 . BEM method combines blade element and momentum theories. The momentum theory balances the forces at the blade under the conservation of linear and angular momentums in a control volume approach, while the blade element theory analyses the forces on a blade section as a function of blade specific geometry. Thus, BEM theory relates blade shape to the forces of wind and of significant consideration among them here, the rotor power extraction. This work uses the classical BEM theory described in 19 together with some corrections that will be explained in more detail.
The rotor introduces two velocities in the field. An axial velocity component in the opposite direction of wind, whose magnitude is equal to aU 0 , and a tangential velocity component in the opposite direction of rotor rotation, equal to a � Ωr in the rotor plane at some radial distance r along F I G U R E 1 Wind turbine blade section nomenclature at some radial position r | 2171 AZIZI And JAHAnGIRIAn the blade; where a and a ′ are known as axial and tangential induction factors, respectively (see Figure 1). Assuming uniform and steady wind conditions upstream of the wind turbine, each blade section experiences normal and tangential velocities as: Among different references, there is a discussion on how the induction factors in momentum and blade element theories are linked when combining them. In the momentum theory, the wind turbine is suggested to be an actuator disk with an infinite number of blades and azimuthally averaged loading, while in the blade element theory, the exact number of blades is inserted directly into the equations and induction factors are local to each blade. The discrete number of blades in momentum equations are accounted for via a tip loss fac-torF tip , and in blade element theory via a solidity factor of , both defined as: where is the angle of the relative wind to the plane of blade rotation or tan = U n U t . It should be recalled that there are a variety of models for tip loss factor in the literature, and here the classical model of Prandtl is utilized.
There is controversy on how to apply the tip loss factor; whether it should be applied to the momentum change, the flow rate, the induction factors, or a combination of them. Here, the suggestion of Glauert on the application of the tip loss factor is implemented, it is applied such that a ∞ = a b F tip and a � ∞ = a � b F tip , where a ∞ ,a � ∞ are azimuthally averaged values, and a b ,a ′ b are local to each blade element. Starting from momentum and blade element theories, and after some manipulations, one obtains two formulae below for axial and tangential induction factors when momentum and blade element theories are combined (dropping sub- where c n is the nondimensional projection of forces in the normal direction, that is, c n = C l cos + C d sin and similarly c t = C l sin − C d cos .

| BEM theory corrections
In this part, some common corrections to the BEM theory are introduced and applied to the current code to make WT power curve prediction more accurate. Firstly, the high thrust correction model of Spera is applied to overcome the notorious momentum theory breakdown problem. When the axial induction factor from Equation (5) exceeds a value of a c = 0.34, the following formula is replaced for axial induction factor calculation 19 : Again, as suggested in, 19 instead of Equation (5), the following equation by Harman is used for tangential induction factor calculation, where r = rΩ U 0 is local speed ratio: Moreover, a relaxation factor of 0.3 ensures convergence of BEM loop iterations on each blade segment.
Similar to the tip loss factor, a hub loss factor can be defined as: And the total loss factor is defined as: BEM method requires aerodynamic coefficients of used airfoils in blade construction, usually obtained from 2d wind tunnel measurements, here denoted by C l and C d for lift and drag coefficients of the profiles, respectively. Since these coefficients are nonlinear in origin, BEM equations will become nonlinear, therefore.
Going back to the 1950s, Himmelskamp recognized that rotating airfoils experience lift coefficients higher than stationery, or 2d sections, as happens in case of rotating WT blades. 20 The Himmelskamp effect, in short, is stall delay at high angles of attack. From then, corrections for including this effect have been the discussion topic of various researches, such as. 21 For the Riso WT, aerodynamic coefficients are given in annex xviii 22 ; although corrected for the actual Reynolds number of Riso WT, these coefficients are still 2d. Hence, there must be a correction applied to take into account the Himmelskamp effect, usually known as a 3d correction.
Wind turbine blades typically experience large angles of attack, beyond 2d wind tunnel measurements, so there must be a method of extrapolating aerodynamic coefficients to the high angles of attack. One such popular method is the Viterna-Corrigan method 23 that relies on flat plate assumption to obtain poststall lift and drag coefficients. Although aerodynamic data for Riso WT blade airfoils are available for incidence angles from 0 to 90 [deg.] in annex xviii, this work uses the Viterna-Corrigan (V-C) model for implementing a 3d correction model as described in paragraphs below.
For the Riso WT, there are two 3d correction models proposed and validated against experimental data in the two articles by Martinez et al 21 and Lanzafame et al. 24 Unfortunately, in the work of Martinez et al, they refer to a formula in, 25 which is proved to be incorrect in the second edition of the same book, 18 so their BEM code results will be referenced here as "incorrect formula." Although the model of Lanzafame et al fits well with experimental data at lower wind speeds, for higher speeds, a new modified 3d model proves to match better to the data. Inspired by the work of Martinez et al, a modified 3d correction model is applied to 2d airfoil data of Riso WT. Following, 23 the V-C model defines: where C l s and C d s are lift and drag coefficients at an angle of attack s , respectively, and C d max is the maximum drag coefficient at the angle of attack = 90 • . As stated in the V-C model, C d max depends on the blade aspect ratio and equals to: Martinez et al suggested s be set at somewhere in the fully developed stall regime, at the angle of minimum lift coefficient, but the original V-C model states that s is taken at the onset of stall, usually s ≃ 15 • . This paper follows the latter, that is, s = 12 • just below the stall angle for all sections along the Riso WT blade. On the rest of computations, a similar approach to Martinez et al method is taken: That is, for angles of attack below s = 12 • , the 2d data from annex xviii 22 database are used, while for angles higher than s , the lift and drag coefficients are averaged of 2d and the V-C model coefficients, ie A similar equation for the drag coefficient is also used. During validation tests, a stall delay model, such as the model of Snel, 26 is included in the above equations, but the results did not show much improvement, hence neglected. During BEM iterations, aerodynamic coefficients for some slight negative angles of attack are also needed, which can be obtained from XFOIL solver. 27

| Baseline wind turbine
In this study, a namely 100 [kW] HAWT with a diameter of 19 [m] is chosen for optimization and related comparisons; the known Riso wind turbine is well described in the publicly available and valuable ECN annex xviii 22 database. It is a three-blade stall-regulated HAWT with upwind rotor, tapered and twisted blades, a hub height of 30 [m], and the compromising airfoils start from NACA 63-225 at the root (blade radius of 2.7 [m]) to NACA 63-212 at the tip (blade radius of 9.5 [m]). Cut-in wind speed is 4 [m/s] and for cutout, although not clearly stated in annex xviii, a feasible assumption of 25 [m/s] is utilized. The field measured power curve of Riso WT is presented in annex xviii. The validation of some existing BEM codes against measurements has been a discussion topic in various publications, such as Méndez and Greiner, 5 Ceyhan, 7 Martinez et al, 21 and Lanzafame et al. 24

| BEM code validation
Riso WT blade geometrical parameters are given for 9 sections along the blade in annex xviii. 22 However, for the BEM code to be precise, the number of blade sections must be increased. So, 49 blade elements are selected along the blade and geometrical properties, and aerodynamic coefficients are interpolated using a spline interpolation method from the original 9 sections. The plot of current 3d corrected BEM code results vs experimental data and other BEM codes available in the literature, as described in Section 2.3, is shown in Figure 2. Also, the power curve without the 3d correction model is plotted and is tagged as "2d Data." The convergence criteria for induction coefficients calculations are set such that the sum of changes in both axial and tangential induction (10) C l = 1∕2C d max sin 2 + K l cos 2 sin

Number of blade elements
Notably, all calculated powers in this work are mechanical; but it can be easily converted via the efficiency factor of gearbox and generator to electrical power for the interested readers.

| Annual energy production calculation
The annual energy yield of a WT at specific wind conditions is calculated using 2 : where V 1 , V 2 denote cut-in and cut-out wind speeds, respectively, and, 8760 is the number of hours in a year.P(u) is power at speed u and Pr (u) is the probability of occurrence of u. In many books, this equation is numerically integrated and expanded using the trapezoidal rule, but since there is already a Simpson integration code written in the current BEM code, the AEP is also calculated by Simpson integration rule for more accuracy. Generally, for representing the wind regime in a specific region, the Weibull probability density function is used, given by: where A [m∕s] is the Weibull scale parameter and k is the nondimensional shape parameter, both specific to each site and must be obtained from actual wind measurement data by statistical methods.
In Figure 3, the annual energy production variation with the computational time for different numbers of blade elements is plotted for the site of Khaf. According to that, 49 blade elements are found to be reasonable; since AEP error is less than 0.4% compared with the results of 269 blade elements calculations, with the computational time, only about 17.3% of corresponding 269 blade elements. By repeating the calculations for other sites and obtaining similar results, the rationality of selecting 49

| SITE SELECTION
The eastern part of Iran is selected as the study area since the wind potential there is relatively good, and there are many disperse small towns that are far from power plants.
Weibull parameters for three different towns in this area are gathered and used for AEP calculation: The first town is Binalood located in the center of Khorasan Razavi province, in the northeast of Iran, famous for strong winds during all year. The next town is Khaf, also located in the south most of Khorasan Razavi province, and is known for its ancient windmills or Asbads (a colloquial Persian translation for windmill), which still are in use today. The third town is Zabol, located in Sistan and Baluchistan province in the southeast of Iran, just below South Khorasan province. Zabol is globally known for its "120-day winds," a continuous dust storm that blows all the summer long. The earliest world recorded windmill design goes back to 1300 AD, 28 a vertical axis windmill in Sistan within the same province, and interestingly some were still in use in recent years. Today all three towns mentioned above own a wind farm and produce electricity. Iran wind map at 40 [m] height is also shown in Figure 4 with those three towns marked on it. The Weibull parameters at the hub height of Riso WT; 30 [m], and also site altitude and corresponding air density along with some statistical indices are given in Table 1. Air density is calculated from site altitude using the formulae presented in. 29 The Weibull distributions are plotted and compared in Figure 5; although all three sites are located in the east part of Iran, they have different wind speed distributions during a year, which make the design of a single optimum WT blade a daunting task. For Zabol, Weibull parameters are initially given at 10 [m] height in, 30 so Equations (18)  where U z is the wind speed at height z, U r is the measured wind speed at height z r , and a is power-law exponent, usually accepted for most onshore applications a = 0.15. For Weibull shape factor extrapolation, the following empirical formula by Justus et al can be used 31 : where k z is the shape parameter at height z, k r is the shape parameter at height z r , and the two constants z ref = 10[m], c = 0.088.

| GENETIC ALGORITHM OPTIMIZATION
Initially pioneered by John Holland in the 1960s-the genetic algorithm is a numerical optimization and parameter search method inspired by Darwin's evolution theory or "survival of the fittest." GA can find the global optimum while dealing with a large number of design variables. This is mainly due to the fact that it improves a set of solutions instead of one as in the case of traditional optimizers. It also has the diversity property that tends to find the global optimum among many possible local optima. Furthermore, GA does not require any knowledge of derivatives and is suitable for problems with a large number of design variables and nonlinearities that are usually involved in aerodynamics applications; while the traditional optimizers will become inefficient in such problems. Hence, GA has proved to be a robust technique for wind turbine applications. 5,7,35,36 MATLAB release R2016b optimization toolbox is utilized for GA implementation in this work, and the BEM code in FORTRAN is used as WT performance computation engine, coupled to MATLAB GA via MEX files. The computation procedure follows the flowchart shown in Figure 6. Design variables are 9 section twist angles for twist optimization cases, and 9 section chord lengths along with 9

| Design parameters bounds
For twist angles, the interval −20 • ≤ (r) ≤ +40 • is chosen and for chord length, the interval is selected, where c is a constant and c 0 (r) is the initial Riso WT chord distribution. This is done to prevent abrupt changes in chord distribution. When AEP is set as the objective function, GA tends to decrease the blade area near the root and increases it at the tip, which is in contrast with blade structural aspects. Thus, following, 5 a value of c = 0.1 is found to yield logical results.

| Design parameters constraints
Section twist angles are selected such that they decrease from root to tip, or (r i+1 ) ≤ (r i ); to ensure that all groundless solutions are omitted during GA optimization. But there is no constraint on the root section twist angle. For section chord lengths, similarly, this constraint is set as c(r i+1 ) ≤ c(r i ) . Another constraint is applied on section chords to make the optimized blade area is equal to (or smaller than) initial Riso WT blade area:� c(r)dr = 1 where A 0 is the initial Riso blade area. This prevents too much change in section chord particularly to keep the blade weight unchanged and prevent unwanted structural changes which may lead to blade structural failure.

| GA objective function
Because the MATLAB GA tool can only function as a minimizer, the objective function is set to the minus of AEP of the WT in [MWh]. For multi-site optimizations, sum of AEP at the sites is selected as the objective function.

| GA validation
The GA tool is validated against a similar attempt on WT twist and chord distribution optimization by Méndez and Greiner, 5 with the objective function of maximizing AEP for    Figure 7 against the Riso characteristics. As shown, both twist and chord distributions are consistent enough with those in. 5

| RESULTS
The GA optimization process is run initially on the three towns of Binalood, Khaf, and Zabol, separately, and then, on all possible combinations of them: Binalood-khaf, Binalood-Zabol, Khaf-Zabol, and Binalood-Khaf-Zabol; making a total of 7 optimization cases. When taking into account the twist angle, and both twist angle and chord length distributions as design parameters, the number of cases is increased to 14. A number of 5 random GA runs are performed for each case to confirm the GA solution.

| Twist angle distribution optimization
The first 7 cases are those regarding the optimization of the twist angle distribution. Final AEP values in [MWh] for all twist optimization cases are given in Table 2. While columns represent different optimum blades, each row corresponds to a specific single or multiple-site; thus, the performance of initial Riso WT and all other optimized blades are evaluated and compared on 7 different site complexes. Detailed optimization results are also compared with the initial Riso WT, and corresponding results are shown in rows tagged as "%Riso." Italic values show AEP of Riso WT at each site, and diagonal bolded values represent optimization results of each case, and other values, although not related to any optimization, are given for performance comparison of each blade at other sites for interested readers. For example, by reading values on column "2.Khf Opt.," row "Bnl," the reader finds that Khaf optimum blade produces 175.93 [MWh] energy if located in Binalood and that is 21.59% less than that of initial Riso WT at Binalood; next row shows that Khaf optimum blade produces 448. 16 [MWh] in Khaf, which is 12.08% higher than Riso WT at Khaf, and so on. As tabulated, by optimizing the twist angle distribution, a total amount of 967.51 [MWh] energy are produced during the year at the three selected sites, which is 1.73% higher than the Riso WT. In Figure 8, two sample GA convergence diagrams are shown. In Figures 9 and 10, the optimum twist angle distribution for single-site and multi-site optimums is plotted, respectively. While twist angles for Binalood and Zabol are in the same range, Khaf has a higher twist angle range, which makes sense as mean wind speed in Khaf is higher. Also, in multi-site optimums, except the Binalood-Zabol case, twist angles are more toward Khaf optimum twist angle. Since AEP increases with the cube of wind speed, and the mean speed in Khaf is higher, GA multi-site optimum is moved toward it, whose contribution in total AEP is higher.

4.45
Bold values represent main optimization results for each case.

| Twist angle and chord length distributions optimization
Similar to the twist angle distribution optimization case above, Table 3 compares the total AEP values for all multisite optimization cases of twist angle and chord length distributions. Rows represent different sites, and each column corresponds to an optimized blade; for a detailed description of data presented, please refer to the description of Table 2 under Section 5.1. As can be seen, the multi-site optimization of the twist angle and chord distributions has increased the total AEP by 4.45% in comparison with the Riso WT, a value that has not been achieved by any of the other optimum blades in this study.  Figure 11 shows two sample GA convergence plots for twist and chord optimization runs. Figures 12 and 13 compare twist angle distribution for single-site and multi-site optimums, respectively. Similarly, chord length distribution is plotted and compared in Figures 14 and 15, with each blade area printed on it. Again, similar to the twist angle optimums, the multi-site optimum twist and chord distributions-except Binalood-Zabol case-are more toward Khaf results; because of the higher share of Khaf in total AEP as its mean speed is higher.

| Wind power density
An interesting study here would be to find out the amount of extractable energy from wind in each site at different wind speeds.  One can multiply the amount of power per area at each wind speed with the probability of that wind speed from the Weibull distribution, and the result gives the distribution of power per square meter of WT swept area at different wind speeds, usually called the power density. The plot of wind power density vs wind speed is drawn in Figure 16 for Khaf; the vertical axis represents the amount of wind power per area multiplied by probability. Total available wind power at different wind speeds is also shown and is named as total WPD. Then, from that value, the highest possible amount extractable by the most efficient WT belongs to the Betz optimal WT, whose efficiency equals 16/27. Power extraction of Riso WT and blades with optimum twist and chord distributions is also shown. While the Betz limit is theoretical, in practice, stall-regulated wind turbines have efficiencies of almost half of the Betz limit. This figure clearly shows that while the Khaf single-site optimum blade tends to extract more energy from the higher wind speeds, the triple-site optimum has more focused on lower speeds in which the other two locations (Binalood and Zabol) have higher probabilities.

| CONCLUSIONS
The performance of multi-site optimization of a HAWT blade was evaluated for three specific sites of Binalood, Khaf, and Zabol, with different wind speed distributions in the east part of Iran. The blade element momentum theory was used as the turbine output power computation engine, and the genetic algorithm was utilized as the optimizer. Total annual energy production at the sites was defined as the objective function. Twist angle and chord length distributions along the Riso WT blade were optimized for all single, double, and triple combinations of those sites. It was shown that performing optimization for one specific site, or even two sites did not guarantee that the optimized blade will perform satisfactorily in other sites; while the total annual energy production at the three sites with the Riso blade was 951.01 [MWh], the multi-site blade with optimum twist angle distribution gave a sum of AEP equal to 967.51 [MWh] with 1.73% increase. Additionally, the blade with optimum twist angle and chord length distributions increased the sum of AEP at the sites to 993.33 [MWh], by 4.45% improvement with respect to Riso WT. Thus, it was proved that among all other optimum blades, the multi-site optimum blade has the highest sum of AEP at the selected sites, and also, the chord length distribution is more effective than the twist angle for achieving higher energy output.