1 Introduction

Large scale wind turbine techniques have become a favorable form of clean energy transformation during the last two decades. The aerodynamic characteristics of wind turbines play crucial roles in the capture and conversion of wind energy. Engineering aerodynamic methods, mainly based upon the blade element momentum (BEM) theory or lifting line model^[1], have been used for a long time. These methods are used in the applications such as rapid calculation of aerodynamic loadings and wind turbine blade design. In this century, with the worldwide growth of wind turbine industry, both the sizes and power of wind turbines increased quickly every year. This makes the accurate prediction of wind turbine aerodynamic performance much more challenging than ever before. Those traditional measures are no longer suitable for both aerodynamic analysis but also design of the mega-Watt wind turbines today. In many cases, engineers and researchers from the wind turbine industry are often facing the following aerodynamic problems.

Firstly, large scale wind turbines all work in the extremely complex atmospheric environment. In some areas, they may have to survive extreme meteorological conditions, such as strong gusts or typhoons. Besides rotating with the blades, the flow around the wind turbines shows significant unsteadiness and high turbulence. With the unsteady inflow, the dynamic stall phenomenon is inevitable and makes it difficult for an accurate temporal prediction of aerodynamic blade loading. Furthermore, the rotation of blades also produces both the centrifugal and Coriolis forces. They bring about the local radial flow along the blades. Without deliberate correction of the rotational effects, classical BEM-based engineering analyses often under-predict the dynamic power output of wind turbines.

Secondly, the interaction between the wind turbine flow and the cascades of atmospheric turbulence structures is dominant at all heights of rotating wind turbine blades. Since the mechanism is not thoroughly investigated yet, proper modelling techniques are still required for depicting the influence from turbulence on phenomena, such as three-dimensional dynamic flows and blade boundary layer separations.

Thirdly, most traditional analysis and design methods for both airfoil profiles and wind turbine blades are borrowed from the techniques for high-speed aeronautical airfoils and wings. But the common working conditions of wind turbines are limited with the much lower wind speed. Correspondingly, most state-of-arts wind turbine blades and their airfoil profiles show the higher thickness, different flow separation patterns, leading-edge roughness sensitivity, and other distinctive characteristics. Here, the aerodynamics of low-speed incompressible flow is of great importance for the performance of wind turbines. As a consequence, many widely used wind turbine aerodynamic analysis and design techniques have to be modified for better predictions.

In contrast, the computational fluid dynamics (CFD) methods, which are established both on the active researches of high-efficiency numerical computation algorithms and on the rapid development of computer engineering, show significant advantages over the traditional theoretical and empirical analysis methods^[2]. Numerous academic efforts have been exerted in these fields. Progress and achievements are so remarkable and are out of the coverage of this article. It has been proved that in comparison with theoretical methods, the modern CFD methods can be used for the numerical simulation of flow field with extremely complex geometries. Today, numerical implementations based on the multi-block structured grids or unstructured meshes are both available. For the flow field around rotating wind turbine blades, the three-dimensional unsteady Reynolds-averaged Navier-Stokes (RANS) equation system can be solved by high-order numerical algorithms for the accurate prediction of flow characteristics of wind turbines. Effects from turbulent structures and the three-dimensional flow separation can be captured with properly chosen turbulence models.

In this paper, the CFD methods are used for solving the fundamental aerodynamic problems of large scale wind turbines. The recent progresses on this research are also briefly reviewed.

2 Aerodynamic aspects of large scale wind turbine airfoils

For the three-dimensional unsteady turbulent flow problems, the numerical methods available today are mainly categorized into three branches as follows: the direct numerical simulation (DNS), the large eddy simulation (LES), and the RANS solution. They roughly cover different scales of turbulent flow structures. For a convincible numerical solution of high Reynolds flow, both the DNS and LES methods have critical basic demands on the spatial resolutions of computational meshes^[3]. Today, the results from the LES, the detached eddy simulation (DES) methods, and the hybrid LES/RANS methods have already been reported for the numerical prediction of wind turbine flow fields^[4-5]. However, this article shows that, taking into account both the required accuracy and the computation costs, the RANS solver is still the most feasible way for understanding the key aerodynamic aspects of large scale wind turbines.

Take the example of the flow around an S809 airfoil^[6]. This airfoil has a considerable thickness. The RANS prediction for its lift and drag features is observed with the serious dispersion among several typical turbulence models. This is true especially at the stall or post-stall stage. At the large angle of attack (AOA), there occurs a large area of flow separation near its leeward surface. The numerical comparison of the lift coefficients is shown in Fig. 1. Here, the $k$-$\omega$ shear stress transport (SST) model also shows a considerable discrepancy from the experimental results, although it may be considered as the best turbulence model for prediction of separation flows today. During our investigation, it is found that a good prediction of flow transition point as well as the correct simulation of the laminar separation bubble is quite important in the overall separation flow simulation of this airfoil. From the sensitivity analysis of local parameters in the $k$-$\omega$ SST model, an SST-based parameter correction can provide much better stall/post-stall lift results for S809-like mid-/high-thickness airfoils. In this way, considerable improvements are obtained both for prediction of the stall AOA and for the maximum lift coefficient (see Fig. 2) .

Two-dimensional aerodynamic aspects of airfoils are the bases of the overall performance of wind turbine blades. They are closely related with each other in many cases. As shown in Fig. 3, not a coincidence, the highest prediction error of the sectional tangential force $\Delta C_{\rm T}$ of both the cut of r/R=47% span of the National Renewable Energy Laboratory (NREL) phase VI blade and an isolated S809 airfoil are located at the same range of the AOA, say, near 20°. Here, $r$ is the spanwise coordinate, and $R$ is the blade span. Therefore, the transition and flow separation of wind turbine airfoils at the large AOA are worthy of investigation, ahead of the full three-dimensional blade analysis^[6].

Meanwhile, in order to obtain a better prediction of the transition location, the Michel's relation is used, and the computation error of skin friction is considerably reduced. This also brings a better simulation result for the total drag (summed by the friction plus pressure drag), as shown in Fig. 4. With this correction, a more accurate simulation of dynamic stall for pitching S809 airfoil is obtained. Compared with the full turbulence simulation based on the default $k$-$\omega$ SST model, the current method shows a better ability of capturing the limiting force cycle of such a pitching airfoil (see Fig. 5) .

3 Efficient numerical flow simulation for large-scalewind turbine blades

The complexity of wind turbine working conditions results in the typical three-dimensional flow rotation, unsteadiness, dynamic stall and considerable turbulence effects. In those beginning days, the CFD-based methods are not more advantageous than the classical analysis and model-based prediction methods. In 2000, the blind code tests of phase VI blades coordinated by the NREL of United States showed a great discrepancy among all computational results from the experiment data. The output power of wind turbines at low wind speeds spread across 25% to 175% from the experimental reference. At large wind speeds, the dispersion of the results reached 35% to 275% ^[7].

In order to attack against the above problems, it is critical to choose proper numerical algorithms and to take into account the three-dimensional dynamic response and rotation effects of wind turbine flows. At all testing conditions, the wind speeds fall in the range of incompressible flow. If a fully coupled compressible Navier-Stokes solver is used for such a low-speed flow case, the reduction in the computational efficiency or accuracy would be expected. Moreover, the prediction of dynamic stall is also the care of wind turbine engineers. Therefore, a feasible strategy for the numerical computation of wind turbine aerodynamics is to develop high-accuracy numerical algorithms for low-speed incompressible flows with proper preconditioning and flow separation-capable turbulence models. Then, the solution of the RANS equation system can be used to capture the flow structures around wind turbine blades, and the key characteristics of wind turbines can be obtained.

In the current researches, the aerodynamic performance of an NREL phase VI blade is thoroughly investigated by the CFD method within a finite volume framework. Figure 6 is the mesh topology as well as the comparison of low-speed shaft torque (LSSTQ). It indicates that, with the $k$-$\omega$ SST model, a considerable improvement can be obtained over those typical CFD results in the original blind tests^[8]. This is true especially in the wind speed range beyond 14m/s. The grid independence is verified on two meshes, one with 22 million cells, the other one with 6 million, respectively.

At the mid-speeds (about 10$-$13m/s), the flow near the blade surfaces will show the partial flow separation along the span-wise direction. This will cause the occurrence of the three-dimensional separation pattern and appearance of spiral points on the surface, as shown by the limiting streamlines near the blade surface in Fig. 7. Such a flow structure generally means an abrupt span-wise change of aerodynamic loading. Since no efforts in the above computation have been exerted yet for the prediction of transition, this may cause the discrepancy of torque values from experiments in this speed range (see Fig. 6) .

Besides the computation for the loadings of an isolated blade pair, its wake effects are also investigated by the CFD methods. For a full comparison with experiments, a one-to-eight scaled model is used. Figure 8 gives the comparison of vorticity in the blade wake between the CFD simulation and particle image velocimetry (PIV) measurements. An agreement has been found here for the vorticity magnitude or the spatial distribution of wake vortices.

Based upon the above work, the CFD method is further used for the flow simulation of a 1.5 MW NH-1500 blade. Its performance analysis has been compared with the results of the free vortex wake (FVW) method. The partial mesh topology is shown in Fig. 9. Based on the computation, the radial distribution of the axial induction velocity is captured. The comparison from both the CFD and the FVW results in Fig. 10 shows that the axial induction factor reaches 1/3 near the blade disk, and the axial velocity component near the wake boundary is far greater than the inside one. The conclusions agree very well with the classical theories. The vorticity in Fig. 11 gives the axial dissipation of wake vortices and the radial vortical expansion in the near wake, which agrees considerably with the prediction from the actuator disk theory.

4 Multi-object optimal design of large scale wind turbine blades

In addition to the satisfactory aerodynamic performance, there are also requirements on the structure and weights of wind turbine blades. Based on the genetic algorithm, an efficient multi-object optimization is now available for the integrated blade design, with both the power output and the material weight being considered^[9-10]. This method takes the basis of the rapid non-dominated sorting genetic algorithm. A more efficient constraint implementation replaces the penalty method^[11-12]. Elite control and dynamic distribution evaluation techniques are introduced to improve the distribution of population for higher efficiency in cases with arbitrary numbers of objectives.

This method is used for the optimized design of a 5 MW blade^[13] and a 1.5 MW NH-1500 blade. The objectives are the maximum of the power coefficient $C_{\rm p}$ at the design wind speed and the minimum of the blade weight. Take the example of the 5 MW blade. The Pareto optimal solutions of the 5 MW blade are shown in Fig. 12. Each point in this set stands for one possible optimal solution^[13]. Their monotonic increasing distribution divides the whole domain into two sections, i.e., the ideal solution I and the feasible solution II. This monotonicity also explains that with the increase in the blade weight, the two objectives contradict with each other. Therefore, in problems like the wind turbine blade design, there is no single solution which satisfies all of the objectives.

The above solutions bring us three typical solution states shown in Fig. 13 (also labelled in Fig. 12) . The comparison among them shows that the chord lengths increase with the power coefficient $C_{\rm p}$ while the twists and thickness change reversely with $C_{\rm p}$. Therefore, the output of power will be higher with the increase in the blade solidity or decrease in the thickness. For example, the case II (see Fig. 12) in the solution set of 5 MW blade with lower weight and manufacturing costs can produce the larger power coefficient compared with the NREL.

The 1.5 MW NH-1500 blade also comes from the optimal design of the above algorithm. Here, the objectives include the maximization of annual power output, minimization of blade weight and minimization of axial thrust. The distribution of both thickness and chord length is shown by the curves labelled by $C_{\rm p}=0.486$ in Fig. 14. For comparison, another optimal solution from the Perato set is also shown in Fig. 14. Although its power coefficient $C_{\rm p}$ reaches 0.5, due to a lower thickness, it is not chosen from the structural viewpoint. Finally, both the $C_{\rm p}$ values come from the evaluation of GH Bladed.

5 Integrated aerodynamic performance of NH-1500 blade

The above NH-1500 blade was manufactured in 2011. Its aerodynamic performance has been evaluated by the FVM and the CFD methods. A one-to-sixteen scaled model has been implemented for wind tunnel validation tests. The comparison of the power coefficient between computational and experimental results can be seen in Fig. 15. For most values of the tip-speed-ratio (TSR) $\lambda$, the results from the FVW methods are slightly higher than the experimental ones. In the design TSR range, the CFD data agree quite well with the wind tunnel results. Its prediction can be verified in this case.

The $C_{\rm p}$ results at the design condition from different methods are listed in Table 1. Here, the scaled model experiments give a power coefficient of 0.492, which is a considerable improvement over its contemporary commercial products. Considering the influence of scaling, the power output of a full scale blade should be even higher than this one. Furthermore, the error between the CFD and experiment results at this TSR value is only $2.6%$, which proves the ability of current methods for the accurate prediction of blade aerodynamic performance.

6 Conclusions

This article introduces the recent progresses on large scale wind turbine aerodynamics and the application of the CFD methods in this field, which is supported by the National Basic Research Program of China (973 Program). The results of the numerical flow simulation of mega-Watt wind turbine blades, multi-objective blade optimal design techniques, and the aerodynamic design and analysis of a 1.5 MW blade are discussed in detail.

Acknowledgements The authors would like to thank the financial support from the 973 Program and the National Natural Science Foundation of China. We would also like to thank Prof. Jingping XIAO, Dr. Li CHEN, and other colleagues for providing the experimental data from wind tunnel tests.