Airfoil properties [25].

## Abstract

The airfoil shape of horizontal axis wind turbine (HAWT) blade is optimized using genetic algorithm (GA). The algorithm is set to find the final airfoil shape with the highest gliding ratio (GR) and larger laminar boundary layer regime along the airfoil surface. The main aim is to find the best airfoil shape of higher lift coefficient with reduced drag in boundary layer from the reference airfoil shape. A 3D correction law is applied to model the effect of optimized airfoil in 3D rotational augmented situation. The thrust and power curves are generated by the blade element (BEM) and free vortex (FV) codes with 3D and loss correction. The higher power production is given when the wind turbine blades are designed using the optimized airfoil. This increment is thought to be made from the efficiency caused by the reduced separation bubbles from reduced turbulent boundary layer and 3D rotational augmentation. To validate its effectiveness in case of soiled condition, the aerodynamic parameters of airfoils are recalculated by enforcing the airfoil to undergo earlier transition, which models the leading edge roughness. The results indicate the soiled condition that does not affect the aerodynamic efficiency of the airfoil due to the positive effect of 3D rotation augmentation.

### Keywords

- optimization
- genetic algorithm
- airfoil
- wind turbine blade
- aerodynamics
- rotational augmentation

## 1. Introduction

The purpose of optimization is to find optimal solutions in scientific or engineering problems. The optimization can also be applied to many wind turbine problems. According to different objectives, constraints, algorithms, tools, and models, various types of optimization solutions are possible. In literature, Jabaraj and Iniyan [1] mentioned optimization models in the wind conversion system with different modelings such as planning model, energy supply-demand model, and forecasting model in 2000. The computational optimization algorithm is mentioned by Ba

The biological concept that survives the fittest individual in the environment among the others is applied in EC [22]. GA, which is one of ECs, runs until it finds the fitness individual with the highest fittest level from the given objective function. It considers the individual solution candidate as the gene, which is the concept of reproduction unit of Mendel [23], and the individuals are exposed to different strategies of genetics to make another generation of the solution candidate pool. The reproduction strategy includes reproduction, crossover, and mutation [24]. The airfoil shapes within given upper and lower bounds make possible solution candidates. The algorithm runs for the fittest individual in objective function f(x) = {GR + Xtr}, which means the algorithm finds the airfoil with the highest Gliding Ratio (GR) and the latest transition point (Xtr), in other words, larger laminar boundary layer regime on the airfoil surface. It finds the best airfoil of the highest lift coefficient (Cl) and lowest drag coefficient (Cd) in the generated airfoil candidate pool [25].

Two airfoils are compared to show their GR, Cl, Cd, and Power production in the wind turbine unit. As the algorithm is run with calculations of the airfoil in 2D, the correction law to consider the 3D effects and Rfoil software is applied [26, 27]. The 3D rotational effect of rotating machines has been found by many. The lift coefficient of the fan blade was found to be three times higher by Himmelskamp in [28]. The lift coefficient of a wind turbine blade was also found to be higher at the inboard sections of the blade by the experiments of Ronsten and Bruining [29, 30]. Later, correction laws for the 3D effect were tried by numerical investigations. The quasi-3D approaches by Hansen [31] and Snel [32, 33] led to the quasi-3D Navier-Stokes mode [34], which confirmed its validity by Shen and Soerensen [35].

As the rotation of the rotor was found to reduce separation and transition by the Coriolis force [35], the 3D correction terms are considered together with the optimized airfoil shape. By comparing the results corrected by the 3D correction law, the effect of optimized shape for higher GR and larger laminar boundary layer in lift coefficient and power production under 3D rotational effects can be deduced. The rotor Power and thrust curves show the combined effect of optimized airfoil on lifting efficiency in the blade unit by BEM theory [36] and FV method [37]. The thrust and power curve comparison leads us to see the effect of the blade lift efficiency increment caused by the optimized airfoil. Moreover, as the total power from the rotor is considered based on each section of blade annulus in BEM, the effects each 3D corrected aerodynamic parameter values of blade section with optimized airfoils are combined to contribute to the increase of power production. Moreover, lifting line of FV method, which calculates total external force and the lift of vortices strength, is also based on the lifting lines on the divided blade segments [37]. Its power calculation also reveals the gathered influence of increased efficiency of the optimized airfoil in sections of the blade.

The compensation for some missing correction laws in power calculations from BEM and FV codes with 3D correction law is possible with the code B-Go. The code B-Go is validated with experimental and computational results, which confirm its reliability based on several correction terms, including tip loss correction, and flow conditions such as massive flow separation takes place [38].

Another realistic rotational situation of airfoil is soiled condition [39]. The leading edge of the airfoil is exposed to debris, dirt, soil, and pollution, and so on in the real situation. The Cl parameter of airfoils is calculated with forced transition. As this contamination is known to decrease the rated and maximum power [40], the roughness sensitivity of airfoil is of importance in the generation of new airfoil. As the contamination usually forces the transition of boundary layer to the leading edge of airfoil, this study made optimized and reference airfoils to have Xtr = 0.05 on the suction side and Xtr = 0.1 on the pressure side as it is recommended in the work of [41].

This chapter illustrates the results of the airfoil and design shape of turbine blade in Section 2, followed by the aerodynamic characteristics in Section 3. The power calculation of the turbine blade with optimized airfoil with corrected BEM and B-Go is shown in Section 4. The airfoil validity in soil condition is elaborated in Section 5. The summarization of the results and their interpretation are shown in Section 6.

## 2. Airfoil and blade design

The optimized airfoil called S809gx is generated with the settled GA algorithm at the Reynolds number (Re)

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

Thickness (%) | 20.99 | 20.3 |

Max. thickness possible (%) | 38.3 | 38.7 |

Max. camber (%) | 0.99 | 0.87 |

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

Because the optimized airfoil is found from the algorithm run to finish at the higher GR and larger Xtr point value at specific angle of attack (AOA) [25], airfoil S809gx shows to have 121% higher GR value, 168% larger laminar boundary layer region on the suction side of airfoil, and 125% larger laminar boundary layer on the pressure side at AOA 7°, as depicted in Table 2. It also indicates 140% higher GR values, 400% larger laminar boundary layer region on suction side, and 162% larger laminar boundary layer on the pressure side at AOA 21.5°. Those AOA values are chosen as the representative angle for fully attached and stall separation flow around airfoil.

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

AOA (°) | 7 | 7 | 21.5 | 21.5 |

Cl | 0.899 | 0.8793 | 1.0264 | 0.9149 |

Cd | 0.011 | 0.0127 | 0.1566 | 0.1958 |

GR | 85.29 | 69.50 | 6.553 | 4.672 |

Xtr (suction side) | 0.272 | 0.162 | 0.016 | 0.004 |

Xtr (pressure side) | 0.677 | 0.540 | 1.000 | 0.616 |

The shape factor H of boundary layer [44] is plotted for both airfoils at targeted angle of attack. The optimized airfoil shows H values to be lower than 2.0 at

Moreover, the trailing edge of S809 suction side also has separation with high H value, while S809gx has the smaller H value at the trailing edge. The pressure side transition is thought to be more violent at the S809 by the extremely different H values of transition point of pressure side of each airfoil, see Figures 2 and 3, In detail, it can be seen from Figures 2 and 3, that the shape factor of the airfoil drops significantly at x/c

To check the validity of the optimized airfoil in soiled condition, boundary layer transition is forced to be 0.05 on the suction side and 0.1 on the pressure side, based on the roughness sensitivity experiment in [41], simulated in Rfoil for its 3D consideration [27].

Although GR values of airfoils are similar in different flow regimes in Table 3, the optimized one shows to have larger laminar boundary layer region over different AOA ranges in the forced transition situation, see Figure 4.

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

AOA (°) | 7 | 7 | 21.5 | 21.5 |

GR | 61 | 61 | 7.53 | 6.67 |

Xtr (suction side) | 0.05 | 0.05 | 0.0135 | 0.0025 |

Xtr (pressure side) | 0.1 | 0.1 | 0.1 | 0.1 |

Figure 2 shows the leading edge of airfoil as background of the graph. The airfoil S809gx has the larger laminar boundary layer region over all angle of attack values. It indicates that the optimized airfoil is shaped to have larger laminar boundary layer region even after transition is forced to be earlier than the normal state. Based on the assumption that the soiled condition triggers earlier boundary layer transition that occurs earlier than clean air condition [41], the optimized S809gx airfoil can be also useful under real air contamination situation [39], which will be discussed further in Section 5.

The wind turbine blade design with the optimized airfoil and the reference one is compared in Figures 5 and 6. The blade is designed with the same twist angle and chord length distribution based on Ref. [42], and the only difference is the airfoil type.

The blades designed with each airfoil are visualized in Figures 5 and 6. The airfoil distribution along the radial position with chord length distribution is based on the NREL Phase VI design guidelines [42], see Tables 4 and 5.

Radial position (m) | Chord length (m) | Twist (°) | Airfoil |
---|---|---|---|

0.508 | 0.218 | 0 | Circular |

0.66 | 0.218 | 0 | Circular |

0.883 | 0.183 | 0 | Circular |

1.008 | 0.349 | 0 | Circular |

1.067 | 0.441 | 0 | Circular |

1.133 | 0.544 | 0 | Circular |

1.257 | 0.737 | 20.04 | S809 |

1.343 | 0.728 | 18.07 | S809 |

1.51 | 0.711 | 14.29 | S809 |

1.648 | 0.697 | 11.91 | S809 |

1.952 | 0.666 | 7.98 | S809 |

2.257 | 0.636 | 5.31 | S809 |

2.343 | 0.627 | 4.71 | S809 |

2.562 | 0.605 | 3.42 | S809 |

2.867 | 0.574 | 2.08 | S809 |

3.172 | 0.543 | 1.15 | S809 |

3.185 | 0.542 | 1.115 | S809 |

3.476 | 0.512 | 0.494 | S809 |

3.781 | 0.482 | −0.015 | S809 |

4.023 | 0.457 | −0.381 | S809 |

4.086 | 0.451 | −0.475 | S809 |

4.391 | 0.42 | −0.92 | S809 |

4.696 | 0.389 | −1.352 | S809 |

4.78 | 0.381 | −1.469 | S809 |

5 | 0.358 | −1.775 | S809 |

Radial position (m) | Chord length (m) | Twist (°) | Airfoil name |
---|---|---|---|

0.508 | 0.218 | −3.00 | Circular |

0.660 | 0.218 | −3.00 | Circular |

0.883 | 0.183 | −3.00 | Circular |

1.008 | 0.349 | −3.00 | Circular |

1.067 | 0.441 | −3.00 | Circular |

1.133 | 0.544 | −3.00 | Circular |

1.257 | 0.737 | 17.04 | s809gx |

1.343 | 0.728 | 15.07 | s809gx |

1.510 | 0.711 | 11.29 | s809gx |

1.648 | 0.697 | 8.91 | s809gx |

1.952 | 0.666 | 4.98 | s809gx |

2.257 | 0.636 | 2.31 | s809gx |

2.343 | 0.627 | 1.71 | s809gx |

2.562 | 0.605 | 0.42 | s809gx |

2.867 | 0.574 | −0.92 | s809gx |

3.172 | 0.543 | −1.85 | s809gx |

3.185 | 0.542 | −1.89 | s809gx |

3.476 | 0.512 | −2.51 | s809gx |

3.781 | 0.482 | −3.02 | s809gx |

4.023 | 0.457 | −3.38 | s809gx |

4.086 | 0.451 | −3.475 | s809gx |

4.391 | 0.420 | −3.92 | s809gx |

4.696 | 0.389 | −4.35 | s809gx |

4.780 | 0.381 | −4.47 | s809gx |

5.000 | 0.358 | −4.78 | s809gx |

## 3. Aerodynamic parameters

Regarding the aerodynamic parameters like Cl, Cd, and GR, two airfoils show similar distribution over the angles of attack. However, the optimized airfoil indicates slightly increased Cl and decreased Cd. Those small advantageous differences are summed up to show increased GR.

The Cl values are corrected with the 3D correction law, mainly considering twist angle and chord per radius ratio of the blade in turbine unit. The corrected value conveys the effect of the Coriolis, centrifugal force, delay of separation, and so on in rotational augmentation [32, 33]. As the rotational effect is significant in lift force, the correction law is only applied in Cl, not Cd. The reference experiments are found in the works of [45, 46]. The calculations are done in Re=

The optimized airfoil Cl values show slightly advantageous over stall angle of attack region compared to the reference one, see Figure 7. The drag coefficient is also smaller than the reference, as shown in Figure 8. Although the airfoil was designed to have better GR value by 2D calculation at the target of angle of attack of 7°, the 3D corrected value also shows the advantage of airfoil S809gx over the different angle of attack ranges.

Cl values corrected by 3D correction law of Hansen [31] in Figure 7 are also supported by the GR calculation of software Rfoil, which improves the 2D prediction with the treated laminar and turbulent boundary layer closure problems [47, 48] in Figure 9. Based on the Rfoil validity [27], increased GR values of the airfoil S809gx especially at fully attached angle of attack range (5–13°) show positive implication for improvement of the following power production in the wind turbine unit.

## 4. Thrust and power curves

To run the code simulation for reference turbine and turbine with optimized airfoil, operation properties are set as Table 6. A fixed pitch value of 3° (equal to turbine 1) is controlled in blade distribution property in Table 5 for turbine 2. Blade 1 is designed with the reference airfoil S809 [42], and blade 2 is designed with the optimized airfoil S809gx.

Turbine 1 | Turbine 2 | |
---|---|---|

Power regulation | Stall | Stall |

Transmission | Single | Single |

V cut in (m/s) | 6.00 | 6.00 |

V cut out (m/s) | 25.00 | 25.00 |

Rotational speed (rpm) | 71.63 | 71.63 |

Outer radius (mm) | 5000.00 | 5000.00 |

Fixed pitch (°) | 3.00 | 0.00 |

Variable losses | 0.22 | 0.22 |

Blade type | Blade 1 | Blade 2 |

Although the advantage of the optimized airfoil in Cl, Cd, and GR values seems to be negligible in Figures 7–9, power production curve shows how blade lift efficiency is improved by those airfoil construction, as shown in Figure 10. The thrust forces are calculated to show the values that are similar at blades with both airfoils, see Figure 11. The similarity in value distribution of thrust implies the power increment that is caused by increased lift velocity from the blade designed with the optimized airfoil.

As the Power Production is calculated by all lifting efficiencies of each blade section [48], combined advantage of each section of blade airfoil produces largely increased power production simulation, especially in the inflow velocity range of 7–25 m/s. Considering the discrepancies between different tools of calculation, the optimized airfoil turbine produced 150% larger power than the reference one. As the optimized airfoil turbine power prediction is based on the simulation of reference one, validated with its experimental data, the discrepancy between experimental data and calculation in high velocity (15–25 m/s) should be considered more.

The thrust curves show similar value pattern with B-Go code, except the fact that thrust force is estimated to be higher at the velocity of flow stall regime, where the prediction can be misled in BEM and FV codes [38]. The B-Go code thrust calculation also supports the increment in blade velocity with the optimized airfoil as the thrust is not increased drastically in the turbine blade designed with the S809gx, see Figure 12.

The power value validation in high velocity, which would cause stall delay in blade [33], is tried with the state-of-the-art-code B-Go. The B-Go has been coded with the off design flow region where separation frequently occurs to make BEM code to be challenged in prediction [38], the power values at high velocity are predicted in spite of the discrepancy with experimental data. Although BEM and FV calculations show the value gap in the velocity of 15–25 m/s region, the B-GO codes show the similar pattern with the other codes. The increasing power values in stall region by BEM prediction show the weakness of BEM at the stall region [37, 38]. It also implies the possibility, which experimental data might have had the error in the stall region, as the experimental set up also has their limit in stall region measurements. The optimized airfoil turbine shows c.a. 150% higher power production in stall region, see Figure 13.

## 5. Performance of the optimized airfoil under soiled condition

As it is shown in Section 2, the optimized airfoil shape shows larger laminar boundary layer regimes even under a forced transition situation, which imitates the soiled condition, as shown in Figure 4. Not only smaller drag within the boundary layer but also the lift coefficient at the forced transition is benefitted because of the enlarged laminar boundary region at optimized airfoil. As it is shown in Figure 14, the Cl value difference of optimized airfoil at forced transition and normal transition situation is negligible. The influence of difference of Cl on GR, calculated by Rfoil, is also demonstrated in Figure 15. The optimized airfoil is less sensitive to changes in inflow conditions. This is not only caused by its 2D characteristics but also being supported by the 3D rotational effects, which delays flow separation and reduces the turbulent boundary layer drag [32, 33]. The optimized airfoil can be tolerable in efficiency decrement in soiled condition or other causes of earlier transition occurrences.

## 6. Conclusions

Stochastic optimization, GA, has been applied to optimize airfoil shape toward larger GR and advantageous boundary layer transition in HAWT. The resulted airfoil shows a 121% higher GR, c.a. 120–170% larger laminar boundary layer on the airfoil surface in targeted AOA (°). The Cl, Cd, and GR values of two airfoils seem to be slight in the results; however, the power production predicted by different codes shows the combined effect of optimized airfoil rotor sections that lead to 150% higher power production. The thrust curves show a similar distribution pattern, indicating the velocity of the blade designed with optimized airfoil influences the power improvement, not the thrust force. The corrected BEM code with 3D rotational augmentation and B-Go codes for the stall region are used for compensating the prediction weakness of BEM in flow separation. The airfoil validity in soil condition is simulated with a forced transition, which shows a negligible lifting coefficient decrement in the optimized airfoil. The laminar boundary layer is still broader at optimized airfoil in forced transition, which indicates that the optimized airfoil shape is useful for realistic airflow with dirt and 3D rotation.

## Acknowledgments

Busan Brain 21 project of BMC is appreciated for its funding in this research.