2. Genetic Algorithm optimization for airfoil

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

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 advan-

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

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

airplane design and off-design requirements were continuously designed.

tageous as it was proven from the laminar airfoil by Jacobs [6].

the curvature shape of the surface [11].

86 Flight Physics - Models, Techniques and Technologies

airfoil-type composition [12].

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:

$$P = [\mathbf{x}; y] \tag{1}$$

$$\text{coefs} = [PP] \tag{2}$$

$$k \text{nots} = -2 \ast a : 2 \ast a + n \tag{3}$$

Figure 1. Airfoil B-spline parameterization with upper and lower bounds.

Figure 2. Genetic algorithm flow diagram.

Airfoil Boundary Layer Optimization Toward Aerodynamic Efficiency of Wind Turbines http://dx.doi.org/10.5772/intechopen.70895 89

$$k = \text{length}(\text{knots}) - \text{length}(\text{coefs}) \tag{4}$$

$$Bsp = spmk(k \text{nots}, \text{coef}) \tag{5}$$

$$\mathbf{x} = \left[ \text{linspace}\left(1, 0, \frac{n}{2}\right), \text{linspace}\left(0, 1, \frac{n}{2}\right) \right] \tag{6}$$

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

$$f = \left(-\left(T\mathbf{X}tr + B\mathbf{X}tr + GR\right)\right) \tag{7}$$

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<sup>6</sup> , and the angle of attack was 7�. The algorithm flow diagram is shown in Figure 2.

#### 3. The optimized airfoil

Figure 2. Genetic algorithm flow diagram.

88 Flight Physics - Models, Techniques and Technologies

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).

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


Table 1. The optimized airfoil S809gx and the reference S809.

Figure 4. Gliding ratio (Cl/Cd) distribution of the airfoil at the root, middle, and tip parts of the blades.

laminar boundary layer. This tendency was found in all flow regime—fully attached, separationtransition, 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 tipspeed 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.

Airfoil Boundary Layer Optimization Toward Aerodynamic Efficiency of Wind Turbines http://dx.doi.org/10.5772/intechopen.70895 91

Figure 5. Drag coefficient distribution of airfoils.

laminar boundary layer. This tendency was found in all flow regime—fully attached, separationtransition, 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 tipspeed 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

Figure 4. Gliding ratio (Cl/Cd) distribution of the airfoil at the root, middle, and tip parts of the blades.

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

S809 S809gx

has the smaller shear wall stress on the top.

Figure 3. The optimized airfoil S809gx and the reference S809.

90 Flight Physics - Models, Techniques and Technologies

Table 1. The optimized airfoil S809gx and the reference S809.

Figure 6. Transition point (Xtr) of the airfoils in fully attached flow regime.

Figure 7. Transition point (Xtr) of the airfoils in transition-separation flow regime.

Figure 8. Transition point (Xtr) of the airfoils in dynamic stall flow regime.

Figure 9. Skin-friction coefficient (Cf) of the top and bottom surfaces of the airfoils.


Table 2. The optimized airfoil S809gx and the reference S809.

Figure 8. Transition point (Xtr) of the airfoils in dynamic stall flow regime.

92 Flight Physics - Models, Techniques and Technologies

Figure 9. Skin-friction coefficient (Cf) of the top and bottom surfaces of the airfoils.


Table 3. Settings for the simulation of turbines consist of the airfoil S809 and S809gx.

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.

Figure 10. The blade designs with airfoil S809 (left) and S809gx (right).

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.

Airfoil Boundary Layer Optimization Toward Aerodynamic Efficiency of Wind Turbines http://dx.doi.org/10.5772/intechopen.70895 95

Figure 11. Annual energy production [kWh] of NREL phase VI turbine with the airfoil S809 and S809gx.

Figure 12. Generated mesh of the reference airfoil S809.

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

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

maximum of 0.1% [18].

under the adverse pressure gradients.

94 Flight Physics - Models, Techniques and Technologies

Figure 10. The blade designs with airfoil S809 (left) and S809gx (right).

The transport equation with the working variable υ~ is mentioned in Ref. [19]. In model assumption, the eddy viscosity is considered as the capability of turbulent flow to transport momentum. The production term is assumed to be increased linearly with the magnitude of the vorticity (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 (fw). The detailed derivation and explanation about equations are found in Ref. [19].

Figure 13. Velocity visualization with LIC, the airfoil S809 (top), and S809gx (bottom) at T = 650.

Airfoil Boundary Layer Optimization Toward Aerodynamic Efficiency of Wind Turbines http://dx.doi.org/10.5772/intechopen.70895 97

Figure 14. Pressure visualization, the airfoil S809 (top), and S809gx (bottom) at T = 650.

Figure 13. Velocity visualization with LIC, the airfoil S809 (top), and S809gx (bottom) at T = 650.

96 Flight Physics - Models, Techniques and Technologies

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.
