Large-scale wind turbine blade design and aerodynamic analysis

Incorporating controlled elitism and dynamic distance crowding strategies, a modified NSGA-II algorithm based on a fast and genetic non-dominated sorting algorithm is developed with the aim of obtaining a novel multi-objective optimization design algorithm for wind turbine blades. As an example, a high-performance 1.5 MW wind turbine blade, taking maximum annual energy production and minimum blade mass as the optimization objectives, was designed. A 1/16-scale model of this blade was tested in a 12 m × 16 m wind tunnel and the experimental results validated the high performance. Moreover, both the computational fluid dynamics (CFD) method and a free-vortex method (FVM) were applied to calculating the aerodynamic performance, which was consistent with the experimental data. For completeness, the CFD and FVM were used to analyze the wake structure, and good and consistent results were obtained between them.

Because of the issues raised by environmental concerns and fossil fuel exhaustion, wind power has been widely accepted as one of the most important renewable energy resources. However, there are still many key problems to be solved in the fundamental research on wind energy application, such as the integrated design of wind turbines, experimental validation and aerodynamic performance prediction.
Wind turbine design cannot be seen as an independent process. There are close links between the components of a wind turbine. This is related to many optimization objectives such as the maximum power output of the rotor, minimum cost of electricity (COE), and minimum tower thrust, as well as meeting the constraints of control, noise and manufacturing [1][2][3][4][5]. Because of its complexity, wind turbine design is mainly based on the cyclical design processes consisting of optimization, verification and modification, which not only depend highly on the designer's experience, but also modifies the optimum aerodynamic configuration of wind turbine blades. Consequently, the performance of the wind turbine is not always ideal. The method of objective weighting [6,7], based on a single-objective idea, is usually employed to solve the multi-objective optimization problem of wind turbine design. This method deals with optimization where the multi-objective problem is converted into a single-objective problem with emphasis on one particular optimum solution at a time. It is not a real multiobjective optimization method and cannot be generally used to solve complex multi-objective problems. It greatly restricts the integrated system design. Therefore, a highly effective and reliable multi-objective optimization algorithm is needed to cope with the multi-variable, multi-objective and multi-constraint optimization problem of the wind turbine.
The aerodynamic performance of the wind turbine is usually studied by two means, i.e. experiment and calculation. Wind tunnel testing can be executed under controlled test conditions and therefore is generally used to determine the wind turbine aerodynamic performance [8]. Recently, a new flow visualization technique, nanoparticle-based planar laser scattering, has revealed the details of the turbulent boundary layer, complex vortex and wake with high precision [9,10], demonstrating great potential for capturing the wind turbine flow structure. The classic blade element momentum (BEM) theory is widely applied to wind turbine aerodynamic calculation, to which the effects of centrifugal force and gravity have recently been introduced for improvements in determining the load of wind turbines [11].
Among wind turbine aerodynamic performance prediction methods, computational fluid dynamics (CFD) is a more general method based on rigorous mathematical derivation and better physical features with fewer assumptions. Another advantage of CFD is that the airfoil aerodynamic data, which are necessary input data for blade element momentum (BEM) theory, is not required. The CFD methods are used in calculations both for wind turbine aerodynamic loading [12][13][14] and for aerodynamic interaction between wind turbines [15][16][17]. Various CFD methods, such as RANS [18], DES [19], LES [20], transition prediction [21], overset grids [22], hybrid schemes of CFD and BEM [23], have been developed and used in wind turbine aerodynamics, providing an unprecedented detailed view and understanding of the flow around and behind the wind turbine.
The vortex method is an alternative in wind turbine aerodynamics because of its higher efficiency and lower cost than the CFD method. A prescribed wake approach [24], which is coupled with a dynamic stall model and a threedimensional rotational effect model, has been used in the Blind Comparison organized by the National Renewable Energy Laboratories [25]. Although computationally efficient, prescribed wake methods are strictly limited in many complex conditions due to lack of sufficient and reliable data for the wake description. In contrast, the free-vortex method (FVM), at least in principle, has fewer fundamental limitations and can be used to model more challenging operation conditions. Recently, various FVM models have been proposed for wind turbine aerodynamic analysis [26][27][28], showing that FVM models can offer more accurate results than BEM theory, which is still widely used in the wind energy industry, and are less expensive than CFD methods. However, they still rely strongly on the accuracy and availability of airfoil aerodynamic data.
In this study, a modified NSGA-II algorithm is used for multi-objective optimization of wind turbine blade design, in particular using a 1.5 MW rated wind turbine as an example. For validation of the aerodynamic design, a wind tunnel experiment has been conducted using a 1/16-scale model of the 1.5 MW blade. A full FVM and a CFD solver have been developed for calculating the blade aerodynamic performance and wake structure. This CFD solver, based on an unstructured finite volume method, employs both a rotating reference frame to simulate the blade rotation and fully implicit coupled schemes to reach quick global convergence.

Blade design
A multi-objective optimization algorithm was developed for wind turbine blade design to achieve high performance. This algorithm was then applied to designing the NH1500 blade especially for a three-bladed, pitch-controlled horizontal-axis wind turbine with a rated power output of 1.5 MW. The wind turbine parameters are shown in Table 1.

The modified NSGA-II
The integrated design of a large-scale wind turbine blade with multi-objective optimization is highly complicated. It is related to many influencing factors such as constraints, optimization variables and objectives, among which conflicts may exist. Moreover, the solutions of objective functions are also highly complicated multimodal optimization problems. To obtain appropriate multi-objective optimization algorithms for wind turbine design, the following three issues must be considered.
(i) The algorithm can handle the optimization problems for any number of objectives theoretically. It is not sensitive to the correlations and conflicts among the objectives.
(ii) The algorithm can effectively solve the constrained problems with any number of constraints, in which the penalty function method should be avoided.
(iii) The algorithm should have good convergence and robustness in solving the specified multimodal problems of wind turbine design.
In response to these specific requirements, a modified NSGA-II algorithm was developed based on a fast and non-dominated sorting genetic algorithm (NSGA-II) [29]. NSGA-II is an efficient sorting algorithm and judges which individual is superior or inferior only according to the dominance relationship without any conversion of the objectives [30]. The algorithm also employs the crowding distance approach to replace the sharing function method, which can control the distribution of individuals without specifying the sharing parameter. For constraint handling, an efficient method named the proposed constraint-handling approach is added to the NSGA-II algorithm [31]. This approach avoids the penalty function method by introducing a virtual value of constraint violation, which can be obtained through transforming the constraints into dimensionless equations, and provides a general method for all of the constraints. NSGA-II has displayed good convergence and robustness and become one of the benchmarks of multi-objective optimization algorithms.
Both the controlled elitism strategy [32] and the dynamic crowding distance strategies [33] have been incorporated into a new algorithm to improve the lateral diversity and uniform diversity of non-dominated solutions. This modification can avoid local convergence of NSGA-II while it retains all its advantages. In addition, to decrease the computational burden and the complexity of algorithmic operation, the simulated binary crossover (SBX) operator and polynomial mutation [34] are also employed in the modified NSGA-II in this study.

Optimization objectives and constraints
An ideal wind turbine blade design is to reach minimum cost of energy under the condition of multiple objectives and constraints. However, the cost of the wind turbine involves many factors. Here, to simplify the optimization, only two conflicting design objectives are chosen for the NH1500 wind turbine blade: (i) Maximum annual energy production in the given wind farm.
(ii) Minimum blade mass under the given structural style. The NH1500 blade pursues the best aerodynamic performance in the given conditions. Consequently the maximum annual energy production is the primary goal in this paper. The modified BEM theory is used in the calculation of the maximum annual energy production. The modifications include the tip loss, hub loss, and correction of the axial induced velocity factor at the large thrust state [35].
The NH1500 blade is based on the preferred structure. The blade structural style is shown in Figure 3, which is primarily made up of double webs [36]. UD GFRP is used in the spar and trailing edge parts, while the 45° GFRP through the leading edge and skin as well as sandwich parts, is filled with PVC foam. For the NH1500 blade, the design loads are quickly evaluated by a set of simplified extreme conditions which are chosen based on IEC61400-1 [37]. According to the design loads, a structural design method is Figure 1 Blade section structure schematic. adopted to ascertain the optimal station and thickness for both spar and trailing edge through the theory of a thinwalled beam element. Just as with the station and thickness, the blade can meet the requirements of strength and stiffness while it consumes a minimum of material.
The determination of the constraints in wind turbine blade design is a semi-empirical process. From the aspect of aerodynamic layout, the NH1500 design only needs to supply some essential restrictions, for which the chord, twist and thickness of the blade decrease from the maximum thickness station to the tip. To cover as many blade shapes as possible and retain blade smoothing, 13 design variables and the natural boundary spines of fourth order are applied in this study. Meanwhile, an approach of smoothing prior to optimization is used in the design to meet the manufacturing requirement and avoid the loss of aerodynamic performance.

Performance of NH1500 blade
The aerodynamic layout of the NH1500 blade is shown in Figure 2, where the position of the blade spar can clearly be seen. The length of the blade is 41.5 m, and the mass is approximately 5.6 tons, which is similar to current commercial 1.5 MW blades. It has been shown that every single objective cannot simultaneously reach the optimum in multi-objective optimization, but a compromise among the objectives is needed [38]. The main concern in this design is the aerodynamic performance of the blade at an acceptable blade mass. The NH1500 rotor power and shaft torque, calculated using GH-BLADED software, are shown in Figure  3. Between cut-in wind speed and 8 m/s, the rotor speed  varies linearly tracking the optimal rotor power coefficient, which can be up to 0.49 according to the GH-BLADED calculation. The rated power is achieved at a wind speed of 10.4 m/s, and pitch control is then activated to maintain constant power output above the rated wind speed.

Experiment
The experimental investigation of the NH1500 blade aerodynamic performance was conducted in the 12 m×16 m low speed wind tunnel at the China Aerodynamics Research and Development Center (CARDC). In the test section, the flow velocity vector diverges less than 2° from the central line. The tested rotor was equipped with three 1/16-scale models of the NH1500 prototype blade (Figure 4) to keep minimum blockage in the tunnel. The blockage factor for this experiment was 11.1%.
In wind tunnel testing, the similarity parameters of the tip speed ratio λ, the Reynolds number Re, and the tip Mach number M could not be the same as the full-scale turbine and for the wind tunnel model. Yet the tip speed ratio is a determining factor of the performance, and hence this similarity must be satisfied. The tunnel flow velocity was settled and then the rotor speed was adjusted to keep λ as consistent as possible with the full-scale turbine. The similarity parameter M was simultaneously satisfied. The torque balance and stress balance were used to measure the torque and force of the turbine model. After the raw data of the wind turbine were obtained, the measured aerodynamic performance was corrected to take into account the solid and wake blockage through the so-called wall pressure signature method.

CFD method
The tip speed of the NH1500 full-scale blade is about 75 m/s, equivalent to a Mach number of about 0.22. Therefore, it is reasonable to assume the density of air as a constant in this simulation, ignoring compressibility.
The governing equations solved in this study are the where Ω is the control volume, S the boundary of Ω, and with the contra-variant velocity V defined as For a rotor in un-yawed oncoming flow, it is convenient to define a reference frame rotating with the rotor to conduct the simulation in a steady manner. In such a rotating reference frame, the Coriolis acceleration and the centrifugal acceleration can be considered in source terms by where   is the angular velocity of the rotating frame, and r  the position vector.
The viscous fluxes v  F are related to laminar viscosity and turbulent viscosity. The turbulent viscosity in this study is calculated from the shear stress transport (SST) k-ω turbulence model. This is recognized as one of the best linear eddy viscosity turbulence models providing a relatively good prediction of turbulence in both attached and separated flows.
The spatial and temporal terms in the equations are discretized using the second-order upwind scheme and implicit time integration, respectively. The variable gradient is computed using the Green-Gauss cell-based method. The pressure gradient and the face mass flux are calculated using the coupled algorithm with fully implicit discretization.

FVM method
A free vortex model (FVM) was developed based on vortex theory. In this method, the blade was modeled as a series of elements, which were represented as the line of bound vor-ticity lying along the blade quarter chord line. The control point of every element was also laid at the blade quarter chord line. The trailing vorticity was modeled as a discrete series of sequential, finite and straight line vortex filaments extending downstream from the trailing edge of the blade element boundary. The shed vorticity could also be modeled to account for time-varying bound vorticity. The wake was allowed to distort freely under the influence of the local velocity field. If the incoming flow was steady and not yawed, the contribution of shed vorticity could be ignored. The wake was cut off after a length of two times the rotor diameter to save computational time.
The induced velocities were determined by means of the Biot-Savart law. To avoid the unreasonable induced velocity when the control point was infinitely close to a vortex element, a modified vortex core model was applied based on the Lamb-Ossen [39] model. The vortex core size was a function of the vortex element movement time taking account of the dissipation effect due to air viscosity [40,41]: where ζ is the wake age angle, ω the rotor angular velocity, r 0 the core radius of the wake vortex at ζ=0, α L Lamb's constant (α L =1.25643),  the kinematic viscosity of air,  v the wake vortex intensity and a 1 a constant (a 1 =0.1) [42]. r 0 in this study equaled a tenth of the local chord for bound vortex and 0.05R for trailing and shed vortices. In the present work, a pseudo implicit technique [43] has been introduced to improve the stability of free-vortex iteration. The wake node position was determined by a relaxation factor between the previous position and the next distorted under the influence of the local velocity field: where w is the relaxation factor and here w=0.5. The subscript denotes the wake node position and the superscript denotes the iterative step in eq. (6). In this work, the initial wake geometry was assumed to be a rigid wake. The FVM was then applied to the whole flow field to generate a new wake structure. Based on this, the procedure was repeated until convergence of the wake geometry was achieved.

Results and discussion
The calculation of both CFD and FVM was carried out for the full-scale 3-bladed NH1500 wind turbine. Figure 5 shows the comparison of the variation of power coefficient C P with tip-speed ratio (TSR) from the experimental data, the CFD calculation and the FVM calculation in un-yawed flow. It can be seen from the figure that the maximum C P was achieved at the optimum TSR of around 9.5 for the experiment and both the calculations. The maximum C P was 0.492 for the experiment, 0.505 for the CFD calculation and 0.528 for the FVM calculation, respectively. This demonstrates that the NH1500 blade is aerodynamically highly efficient. Not only was the experimental maximum C P slightly lower than the calculated maximum C P , but also all the power coefficients from the experiment were lower than those from the both calculations generally. This was reasonable because the Reynolds number for the 1/16-scale model used in the wind tunnel test was certainly lower than that for the full-scale blade at the same tip speed, resulting in lower lift and higher drag. As for the both calculations, the result from the FVM was higher than that from the CFD, since the former only took into account the viscosity in the vortex core but not in the whole flow fields, leading to a higher maximum C P of 0.528, which was too close to the well-known Betz limit 0.593. Nevertheless, it was assured that the maximum C P for the full-scale rotor was higher than the experimental value 0.492 considering the effect of the Reynolds number, and therefore the maximum result of 0.505 from the CFD may be more convincing. Figure 6 shows the finally convergent wake geometry from the FVM. The wake vortex filaments were allowed to develop freely downstream during the simulation. It can be  clearly seen from the figure that the tip vortex and the root vortex rolled up from a certain wake age angle. Figure 7 shows the contours of the axial induced factor from both the CFD and FVM at the wind speed of 8 m/s (TSR=9.36). It is interesting that the agreement between the contours from the two different methods was good and the axial induced factor at and around the blade was close to 1/3, which corresponded to an aerodynamically optimum wind turbine according to classic momentum theory. The axial induced factors around the wake outer boundary and center were much smaller than the mid-part of the wake. It was obvious that the axial induced factor increased from upstream to downstream, decreasing the axial flow velocity, which was consistent again with momentum theory. The comparison shows that both the CFD and FVM methods could well describe the wind turbine wake structure. Figure 8 shows the locations of tip vortex cores from the CFD and FVM calculations, in which the result from CFD is illustrated by vorticity contours and the result from FVM  by solid black circles. The wake vortical strength is visualized through color bands. Three tip vortex cores are clearly allocated in the figure. The vorticity becomes progressively weaker downstream because of viscid dissipation. From the locations calculated by both CFD and FVM, it was quite clear that the wake expanded in the radial direction because of the lower axial velocity downstream (Figure 7). The agreement between the vortex core locations from CFD and FVM was noted to be very good.

Conclusions
A high performance wind turbine blade for a 1.5 MW turbine has been designed using a modified NSGA-II algorithm and the aerodynamic performance has been validated by both experiment and calculation. This study also focused on the wind turbine wake structure analysis through CFD and FVM calculations.
(i) The modified NSGA-II had a good performance in convergence and robustness when handling multi-objective, multi-variable and multi-constraint optimization problems, providing a general approach for wind turbine blade design.
(ii) Both the experiment and calculation showed that the power coefficient of the designed NH1500 blade can achieve a maximum value of more than 0.492, demonstrating that the blade is aerodynamically highly efficient.
(iii) Both CFD and FVM show significant potential for wind turbine flow field calculations. Their results were well consistent with theoretical analysis. They all caught the tip and root vortices and agreed well in the locations of the tip vortex cores.