## 1. Introduction

The airfoil is the geometrically shaped structure for mechanical force generation from the relative movement between the airfoil and surrounding airflow of the airfoil structures [1]. For wind turbines, the airfoil shape of the blades influences the turbine power production. The lifting efficiency of the blades determines the effectiveness of rotor rotation to cause productive energy conversion from wind kinetics to rotor rotation, which leads to higher electricity generation from the drive unit.

From the late eighteenth century, the curving surface geometry was discovered to be advantageous for lifting efficiency in windmill by Smeaton [2]. In the 1880s, Lilienthal discovered the specific shape from the bird’s wings, which inspired the airplane invention by Wright brothers [3]. The research of Prandtl and Tietjens revealed the benefit of the thick airfoil through their mathematical skills and wind tunnel tests in 1917 [4]. The US National Advisory Committee for Aeronautics (NACA)-generated airfoil groups, “NACA airfoil families” in the 1930s, and the series have been widely used even in these days [5]. In contradiction to the mathematical methods to calculate the pressure distribution of airfoil, Jacobs proposed the airfoil design which causes the desired pressure distribution. The laminar flow of airfoil was expanded to cause higher L/D ratio and smaller drag [6]. Later, different types of airfoils for various airplane design and off-design requirements were continuously designed.

The first wind turbine blades were also designed by the airfoils from aeronautic applications. However, in the 1980s, the airfoils specially dedicated for wind turbines were begun to be made due to the defects of aeronautic airfoils applied in a wind turbine. The sensitivity roughness effect on the leading edge arose to be the required element for wind turbine airfoil. The airfoil series for stall-regulated, variable-pitch control wind turbine was developed by NREL in 1984, incorporated with SERI and Airfoils [7]. The wind turbine-dedicated airfoils with the thickness from 15 to 40% of the chord were also made by the team of the Delft University of Technology with the design objective of low sensitivity to roughness, Gurney flaps, and trailing-edge wedge consideration [8]. The airfoils from Risø were designed to have high aerodynamic efficiency and slender blade shape [9]. The airfoil design using numerical optimization for tip region of the blades was researched by Grasso [10]. As mentioned in these studies, the higher aerodynamic efficiency, insensitivity to roughness effect, structural stability and smooth post-stall exhibition, etc., are required for wind turbine airfoil design. To accomplish these objectives, boundary layer consideration of the wind turbine airfoil can be advantageous as it was proven from the laminar airfoil by Jacobs [6].

The boundary layer of the airfoil is exerted by additional pressure generated by the curvature shape of airfoil compared to the constant pressure on boundary layer made of the plate with zero incidences. The pressure distribution on the edge of the boundary layer is same with the pressure distribution on the wall in the plate. However, due to streamline curvature of airfoil surface, the pressure gradients and compensation for the centrifugal force of the streamline flow are generated inside the boundary layer. Furthermore, the transition point of the boundary layer on the airfoil is determined by the outer flow and its pressure difference generated by the curvature shape of the surface [11].

To generate the airfoil shape which has the advantage for pressure distribution in the boundary layer and transition points, genetic algorithm (GA) optimization was used in this study. As all airfoils are designed for higher aerodynamic performance, GA objective functions therefore had two objectives—higher transition points of the larger laminar boundary layer and higher gliding ratio (GR). The airfoil S809 of NREL airfoil series for the wind turbine was chosen as a reference. The shape of insensitiveness to the roughness effect of the airfoil S809 could be maintained in the optimized airfoil. The final evaluations of turbine performance were done with the sample of stall-regulated wind turbine of NREL phase VI, which consists of the same airfoil-type composition [12].

The B-spline parameterization was used for the airfoil description, and the *y* points of the spline were considered to be the variables. The values of boundary layer parameters and GR of the airfoil were calculated by the flow solver XFOIL. The power performances of turbine unit with blades of the optimized airfoil were calculated by using blade element method (BEM) of the software QBlade. The CFD simulations from OpenFOAM^{®} were performed to visualize the improved aerodynamic aspects.

Section 2 explains the GA airfoil optimization method, Section 3 presents the aerodynamic and boundary layer results of the optimized airfoils with improved power production of turbine unit, Section 4 visualizes the airflow of the optimized airfoil with the reference, and Section 5 concludes this chapter.

## 2. Genetic Algorithm optimization for airfoil

The GA algorithm is based on the principle of the survival of the fittest and natural selection, observed by Darwin [13]. As the various bird breaks were developed for different foods that they can survive with, the airfoil was set to be shaped to survive at the condition of the highest GR and transition point. To put the airfoil in mathematical form, the B-spline was used, and its variables were set as the *y* point of control points in MATLAB^{®} (
**Figure 1**
).

The smoothness and number of variables were set according to the previous research for effective GA operation within a given computation time [14]. The airfoil is described with the B-spline. The corresponding equations and the example of description figure are the following. The *x, y* points are for the B-spline control points to compose *P* matrix in *coefs* for determining the smoothness of B-spline with *k.* When the constant *a* is designated, the *knots* are defined. The *x* points of B-spline are defined with MATLAB function *linspace* which divides the *x*-axis space according to the defined airfoil chord:

To achieve higher aerodynamic efficiency and larger laminar boundary layer region at the same time, the multi-objective function was formulated to evaluate both parameters of the airfoil. The objective function code was written as

Then, *TXtr* is the transition point of the top of airfoil, and *BXtr* is the transition point of the bottom part of airfoil surface. The GR values, *TXtr* and *BXtr* values, were given from the results of XFOIL calculation [15]. The given Re number was 10^{6}, and the angle of attack was 7°. The algorithm flow diagram is shown in
**Figure 2**
.

## 3. The optimized airfoil

The airfoil shape with higher GR and transition points were given after the convergence of default set of population and running operation in GA. The thickness was almost the same as the reference, which was desirable thickness for the stall-regulated wind turbine, but the curvature shape was slightly changed (
**Figure 3**
and
**Table 1**
).

S809 | S809gx | |
---|---|---|

Thickness (%) | 20.99 | 20.3 |

Max. thick. pos. (%) | 38.3 | 38.7 |

Max. camber (%) | 0.99 | 0.87 |

Max. camber pos. (%) | 83.3 | 43.6 |

However, this slight curvature shape change generated different GR and boundary layer transitions (
**Figures 4**
–
**8**
). As the gliding ratio was set to be higher at the angle of attack 7°, there was some region (off-design) where GR values of the optimized airfoil are lower than the reference. But the angle of attack range of increasing GR ratio was larger at the optimized airfoil (design point). This means that the optimized one has the smaller range of separation or stall occurrence, which is desirable for the stall-regulated wind turbine [7] (
**Figure 4**
).

Although the maximum GR values of the reference and optimized airfoils were similar, the advantage of the optimized is represented by the drag distribution in
**Figure 5**
. Significantly, reduced drag coefficients are found in
**Figure 5**
, which advocate the lower drag of the optimized curvature shape, especially at the stall regime. The reduced drags can be explained with the enlarged laminar boundary layer region of the optimized (
**Figures 6**
–
**8**
).

The transition point distributions of the reference and optimized airfoils are compared in all angle of attack regime. The top and bottom mean the upper and lower sides of the airfoil. Including the angle of attack 7°, where the optimization was calculated, the transition points of the optimized S809gx showed the higher values than the reference S809, which means a larger laminar boundary layer. This tendency was found in all flow regime—fully attached, separation-transition, and dynamic stall regime. The boundary layer region difference can be also found by the skin-friction coefficient (Cf) (
**Figure 9**
). In the separation-transition flow regime where the tip-speed ratio (TSR) = 5 and the incoming velocity is 8.3 m/s, transition point comparison of two airfoils at the top and bottom parts is presented in
**Table 2**
. The Cf value at the leading edge of the top of the airfoil showed much higher value than the reference one. The Cf value hill range was also higher at the reference than the optimized one which indicates that the optimized one has the smaller shear wall stress on the top.

Although the values were relatively smaller than the top, the bottom part Cf distribution shows the increasing and decreasing tendency change at the transition point. The optimized airfoil shows the changing point to be located to the right side than the reference. If the separation occurs at the point where shear stress gets zero, the optimized airfoil could have the smaller adverse pressure gradient region than the reference one at the bottom side.

The improved aerodynamic efficiency of the optimized airfoil through laminar boundary layer enlargement affected the power increment of rotor turbine with the optimized airfoil. The wind turbine simulation of NREL phase VI blades with the airfoil S809 and S809gx is done with the blade incorporation. The settings for power performance comparison of two turbines are represented in
**Table 3**
(
**Figure 10**
). The wind turbine with optimized airfoil shows higher annual yield and power production at a given velocity condition (**Figure 11**). The increased power production means the effective rotor rotation of the wind turbine with the optimized airfoil, which was affected by lifting efficiency of the blades from the aerodynamically upgraded airfoil composition. The 3D rotation effect to reduce the root vortex could be another reason of increased rotor efficiency of the rotor as explained in Ref. [16].

## 4. CFD simulation

To visualize the improved aerodynamic behavior in the optimized airfoil compared to the reference, the CFD simulations were performed.

The meshes for two airfoils were generated with Gmsh [17] (
**Figure 12**
). The Gmsh tool offers a certain tool for mesh refinement so that the validity of the mesh results can be assured. As the six decimals of accuracy are used, the changes are only appreciated by the fourth digit with the maximum of 0.1% [18].

The software OpenFOAM^{®} uses SimpleFOAM solver with Spalart-Allmaras turbulence model for Reynolds-averaged Navier–Stokes (RANS) equations as the governing equation. It was developed for the aerospace flow problems including wall-bounded flow for boundary layers under the adverse pressure gradients.

The transport equation with the working variable
*S*). The right-hand side of the equation also includes the third term, which is the destruction term. The faster-decaying motion in the outer part of the boundary layer is expressed with the function (*f _{w}
*). The detailed derivation and explanation about equations are found in Ref. [19].

The solver visualized the following airflow with ParaView^{®}. The velocity distribution at the outer flow and boundary layer flow distributions are visualized through the aforementioned solution method. The flow time was set to be 1000 with the time interval of 50. The dynamic stall regime angle of attack was at around 22°; the outer flow TSR was 1.5 when incoming velocity was 27.7 m/s.

In general, the reference airfoil shows larger stall area with more laminar separation bubble occurrences than the optimized one. At the time point *T* = 650, the difference between two airfoils was clearly visualized. At *T* = 650, the slope of the stall area is more smooth at the optimized airfoil with the smaller total area (
**Figure 13**
).

The laminar separation bubbles are only found in the reference when the pressure gradient is less drastic at the optimized one. This is supported by the pressure visualization in
**Figure 14**
. The adverse pressure gradient difference was higher at the reference, and the area with the minus pressure range is largely found at the reference. This milder stall and separation effects explain the smaller drag at the optimized airfoil in dynamic stall regime.

## 5. Conclusion

The airfoil optimization using GA for higher aerodynamic efficiency and larger laminar boundary layer is achieved. Furthermore, the turbine power performance increment and air flow visualization of reduced stall from the optimized airfoil are proven via simulation. The higher GR and transition point values due to the reduced drag from expanded laminar boundary layer region are presented for the optimized airfoil. The boundary layer characteristics such as diminished skin-friction coefficients and higher transition point values over all flow regimes are found to be the reason behind the improved aerodynamics of the optimized airfoil. The results present the contribution of specifically airfoil shape for aerodynamic efficiency and modified boundary layer distribution on HAWT performance improvement including an effective lifting of the blades and rotor rotation.