Optimization of Wind Turbine Airfoil Using Nondominated Sorting Genetic Algorithm and Pareto Optimal Front

Abstract

A Computational Fluid Dynamics (CFD) and response surface-based multiobjective design optimization were performed for six different 2D airfoil profiles, and the Pareto optimal front of each airfoil is presented. FLUENT, which is a commercial CFD simulation code, was used to determine the relevant aerodynamic loads. The Lift Coefficient (𝐶_𝐿) and Drag Coefficient (𝐶_𝐷) data at a range of 0° to 12° angles of attack (𝛼) and at three different Reynolds numbers (Re=68,459, 479, 210, and 958, 422) for all the six airfoils were obtained. Realizable 𝑘-𝜀 turbulence model with a second-order upwind solution method was used in the simulations. The standard least square method was used to generate response surface by the statistical code JMP. Elitist Non-dominated Sorting Genetic Algorithm (NSGA-II) was used to determine the Pareto optimal set based on the response surfaces. Each Pareto optimal solution represents a different compromise between design objectives. This gives the designer a choice to select a design compromise that best suits the requirements from a set of optimal solutions. The Pareto solution set is presented in the form of a Pareto optimal front.

1. Introduction

According to the US Department of energy, the combustion of fossil fuels results a net increase of 10.65 billion tones of atmospheric carbon dioxide every year [1] which deteriorates the environmental balance. Fossil fuels also give out sulfur dioxide into the air, which, after reacting with the moisture in air, produces sulfuric acid and leads to acid rain. Furthermore, depletion of these nonrenewable sources of energy is taking place at a rapid pace because of the increasing demands of energy, with the modernization of our society. Estimates from the US Department of Energy predict that the years of production left in the ground for oil are 43 years, gas 167 years, and coal 417 years [2]. Therefore, it is critical that we start looking for some renewable sources of energy that can be used as alternatives to fossil fuels. Renewable energy resources can play a key role in producing local, clean, and inexhaustible energy to supply the growing demand for electricity, heat, and transportation fuel. Wind energy as a source of energy to produce electricity is favoured widely as an alternative to fossil fuels. It is plentiful, renewable, widely distributed, and clean and produces no greenhouse gas emissions. Wind turbines convert kinetic energy from the wind into mechanical energy which can be used to generate electricity. When the wind blows and flows around the blades of the wind turbines, which have essentially airfoil cross sections, it generates lift forces which makes the blades spin. These blades are connected to a drive shaft that turns an electric generator to produce electricity which therefore can be sent through a cable down the turbine tower to a transmission line. The blades of a wind turbine rotor are generally regarded as the most critical component of the wind turbine system [3]. The aerodynamic profiles of wind turbine blades have crucial influence on aerodynamic efficiency of wind turbines. Even minor alterations in the shape of the profile can greatly alter the power curve and noise level. Some of the important design parameters include the number of blades, blade solidity, blade taper, and twist as well as tip-speed ratio. The aerodynamic theory of the wind turbines gradually developed, starting with the simple one-dimensional momentum analysis [4] of the actuator disc to the more commonly used BEM theory.

The BEM theory is based on the assumption that the flow at a given annulus does not affect the flow at adjacent annuli [5]. This allows the rotor blade to be analysed in sections, where the resulting forces are summed over all sections to get the overall forces of the rotor. The theory uses both axial and angular momentum balances to determine the flow and the resulting forces at the blade. BEM methods are very fast and reliable in the design process, nevertheless these are limited due to their two dimensional nature. These codes require tabulated data for the lift, drag and moment distributions versus the angle of attack to calculate the blade aerodynamic loads. Furthermore, empirical corrections are necessary to account for rotational effects near the root and three dimensional flows around the tip region [6].

Over the last several years, the wind turbine community has started to look at CFD methods to complement wind tunnel [7] and in field tests on the understanding of the complex flow physics around rotating wind turbine blades. CFD codes can be very useful on the calculation of aerodynamic coefficients required by engineering methods and on the explicit determination of loads since no corrections are necessary as in BEM method.

The overall goal of this study is to perform a response surface-based multiobjective optimization of selected 2D airfoil profiles using Elitist Nondominated Sorting Genetic Algorithm (NSGA). In order to achieve this overall goal, several specific objectives were determined. The first specific objective was to identify several airfoil profiles with their geometric coordinates. The second objective was to perform CFD simulations around the airfoils. Simulations for each airfoil were performed for several values of Re and 𝛼. The third objective was to determine response surfaces for lift and drag coefficients as a function of Re and 𝛼. The fourth and final objective was to perform the optimization using genetic algorithm to determine a set of nondominated solution for each airfoil. Based on the optimization results, designers can opt choose multiple airfoils for a single blade depending on Re and 𝛼 variation along the length of the blade after appropriate twist and taper is applied. It is also worth to mention other interesting works [8, 9] that are reported where numerical models are used to design turbines specifically in the wave energy area.

2. Approaches

2.1. CFD Modelling

In this work, we consider the flow around six different airfoil shapes (NACA 63-218, E387, FX63-137, NACA 63-421, NACA 64-421, and NACA 65-421) at 3 different Reynolds number (Re=68,459, Re=479,210, and Re=958,422) for a range of 0° to 12° angles of attack. These airfoils are created from a set of vertices generated from the University of Illinois at Urbana Champagne (UIUC) airfoil database [10]. These vertices are connected using a smooth curve, creating the surface of the airfoil. A flow domain is created surrounding the airfoil and this domain is split for meshing purposes.

The CFD data of the 15 simulated cases for each airfoil were used to generate a response surface. The response surfaces were fit using standard least-square regression with quadratic polynomial using JMP. These response surfaces are obtained between design variable (Re and 𝛼) and objective functions (𝐶_𝐿 and 𝐶_𝐷) for each airfoil profile. All the design variables and the objective functions are normalized between 0 and 1 based on their maximum and minimum values in order to determine the response surface.

2.2. Grid Description

Grid generation is the most important step in the CFD simulations. The quality of the grid plays a direct role on the quality of the analysis, regardless of the flow solver used. Additionally, the solver will be more robust and efficient when using a well-constructed mesh. In this work, structured grids were generated using the commercial code GAMBIT. Figure 1 shows a 2D mesh of the entire domain using a map scheme with around 50,000 quadrilateral elements (the number varies slightly for different airfoils) while Figure 2 shows the blowup of the mesh generated around the airfoil.

Figure 1 

2D mesh of the entire domain using a map scheme with around 50,000 quadrilateral elements.

Figure 2 

Mesh generated around the airfoil.

In order to have a stable solution, the generated grids had the least number of elements with high aspect ratios. To be able to resolve adequately the boundary layer along the airfoil wall, grid points were clustered near the wall. The grids were also clustered near the trailing edge in order to catch the flow separation.

2.3. Boundary Conditions

Boundary conditions specify the flow and thermal variables on the boundaries of the physical model. They are, therefore, a critical component of the CFD simulations, and it is important that they are specified appropriately. In this work, 3 different types of boundary conditions were used: no-slip boundary condition over the airfoil surface, inlet boundary condition for free stream flow, and pressure outlet. The outlet boundary of the domain was set to a constant pressure value. It was set to be atmospheric pressure. The object in the computational domain (i.e., the airfoil surface), around which the flow was simulated, was set to be no-slip boundary (wall). The no-slip boundary condition sets the stream wise velocity to zero. The velocity Inlet boundary condition was used to define the flow velocity at the flow inlet. Figure 3 shows the different assigned boundary conditions for all the CFD simulations as follows.

Figure 3 

The flow domain with the boundary conditions.

In Gambit, the boundary conditions were declared (i.e., wall, velocity inlet, and pressure outlet), but actual values for these boundaries were defined in fluent. For velocity inlet, we used 3 different velocities for each airfoil at every angle of attack. We set 𝑣=1 m/s (for Re=68,459), 𝑣=7 m/s (for Re=479,210), and 𝑣=14 m/s (for Re=958,422). The Realizable 𝑘-epsilon turbulence model and a second-order upwind solution method were used to get more accurate results.

2.4. Response Surface Methodology

The response surface method fits an approximate function to a set of experimentally or numerically evaluated design data points [11]. There are various response surface approximation methods available in the literature [11, 12], with the polynomial-based approximations being the most popular. In this technique, an appropriate ordered polynomial is fitted to a set of data points, such that the adjusted RMS error 𝜎_𝑎 is minimized and quality parameter 𝑅²_adj is made as close as possible to one [12]. The 𝜎_𝑎 and 𝑅²_adj are defined as follows.

Let 𝑁 be the number of data points and let 𝑁_𝑝 be the number of coefficients, and error 𝑒_𝑖 at any point 𝑖 is defined: 𝑒_𝑖=𝑓^𝑎_𝑖−𝑓^𝑝_𝑖, (1) where𝑓^𝑝_𝑖 is the actual value of the function at the design point and 𝑓^𝑝_𝑖 is the predicted value. Hence,𝜎_𝑎=∑^𝑁_𝑖=1𝑒²_𝑖𝑁−𝑁_𝑃,𝑅²_adj𝜎=1−^2𝑎𝑁_𝑝−1∑^𝑁_𝑖=1𝑦_𝑖−𝑦², (2)(3) where∑𝑦=^𝑁𝑖=1𝑦_𝑖𝑁_𝑝. (4) The number of data 𝑁 has to be greater than the number of coefficients 𝑁_𝑝 so that the denominator of (2) is always positive and well posed. Since 𝑅²_adj needs to be as close as possible to 1 to represent a good fit, the terms in the numerator of (3) (𝜎_𝑎)²(𝑁_𝑃−1) should be less than or equal to the denominator ∑^𝑁₁(𝑦_𝑖−𝑦)² so that 𝑅²_adj will always be positive. In this study, the response surface method is applied with two objectives, namely, to generate response surface from the CFD simulation results and to approximate the global Pareto optimal front by representing one objective in terms of others. Both aspects will be discussed in upcoming sections.

2.5. Optimization Approach

The methodology used for generating Pareto optimal front is a multiobjective evolutionary algorithm (MOEA). The specific algorithm used is the Elitist Nondominated Sorting Genetic Algorithm (NSGA-II) [13]. All genetic algorithm codes use some form of sorting scheme to get nondominated solutions. Nondominated solutions are the best solutions. Among the nondominated solutions one cannot be said to be better than the other. NSGA-II uses an explicit diversity-preserving mechanism. The starting point is the identification of constraints, performance criteria, design variables, and allowable range of design variables. The response surfaces generated from the results of the CFD simulations are incorporated into the NSGA-II code. In running the code input values for several parameters must be provided. These parameters and their best values, as suggested by Deb [11, 13] after their extensive parametric study, are as follows:(i)Population size: 100,(ii)Generations: 250,(iii)Crossover probability (P_cross): 1.00,(iv)Distribution parameter (for crossover): 20,(v)Mutation probability (P_mut): 0.250,(vi)Distribution parameter (for mutation): 200.

Where the population size is the size of the non-dominated solutions and the generations are equivalent to the number of iterations. Crossover probability, mutation probability, distribution parameter for crossover, and distribution parameter for mutation are used to create the offspring population from the parent population. The crossover probability is mainly responsible for the search aspect of the genetic algorithm while mutation probability keeps the diversity in the population. The distribution parameter for crossover controls the diversity of the children solutions obtained after crossover while distribution parameter for mutation controls the spread of the solutions after mutation.

In NSGA-II algorithm, the code first creates a parent population of 𝑃_t of size 𝑁. From the parent set it then creates an offspring population 𝑄_t of size 𝑁. The NSGA-II algorithm, instead of finding the non-dominated front of 𝑄_t only, uses the combined population of 𝑃_t and 𝑄_t to form 𝑅_t of size 2𝑁. Then, a nondominated sorting is used to classify the entire population 𝑅_t. Although this requires more effort compared to performing a non-dominated sorting on 𝑄_t alone, it allows a global non-domination check among the offspring and parent solutions. Once the non-dominated sorting is over, the new population is filled by solutions of different non-dominated fronts, 𝐹_𝑖, one at a time. The filling starts with the best non-dominated front and continues with solutions of the second non-dominated front and so on. Since the overall population size of 𝑅_t is 2𝑁, not all fronts may be accommodated in 𝑁 slots available in the new population. All fronts which could not be accommodated are simply deleted. When the last allowed front is being considered, there may exist more solutions in the last front than the remaining slots in the new population. Instead of arbitrarily discarding some members from the last front, a niching strategy [13] is used to choose the members of the last front based on crowding distance. The solutions kept in the population are those which have the largest crowding distance thus keeping the diversity of the solution. This new set of solutions is now the parent set for the next generation. The procedure is then repeated till the best non-dominated set is obtained.

3. Results and Discussions

3.1. Flow Field Overview

The CFD validation is a very important part of computational fluid dynamics. It is used to evaluate the accuracy of CFD results. We compared the CFD results of Lift Coefficient of E387 airfoil for 5 different angles of attack with National Renewable Energy Laboratory (NREL) experimental data [14] for Re=479,210 as shown in Figure 4. The CFD results of Lift Coefficient show good agreement with NREL experimental results.

Figure 4 

Comparison of 𝐶𝐿 Versus 𝛼 of CFD simulation results with NREL experimental data.

Colour contours of static pressure and velocity vectors by velocity magnitude of NACA 64-421 airfoil at an angle of attack of 6° and at Re=479,210 are shown in Figures 5 and 6. The static pressure plot clearly shows the higher pressure (indicated by red, yellow, and green colours) on the bottom surface and lower pressure (indicated by blue colour) on the top surface. In the velocity vector plot, higher velocity corresponding to lower pressure and lower velocity corresponding to higher pressure can be clearly observed.

Figure 5 

Colour contour of static pressure around NACA 64-421 airfoil at 𝛼=6 and Re=479,210.

Figure 6 

Colour velocity vectors by velocity magnitude around NACA 64-421 airfoil at 𝛼=6 and Re=479,210.

3.2. Performance Trend

3.2.1. 𝐶_𝑝 Distribution

Figures 7(a)–7(c) show the 𝐶_𝑝 distribution for NACA 63-421 for the three Reynolds numbers. In each figure and for each angle of attack (uniform colour), the bottom line represents the 𝐶_𝑝 distribution at the top surface of the airfoil, indicating lower pressure, and the top line represents the 𝐶_𝑝 distribution on the bottom surface of the airfoil indicating higher pressure. As the angle of attack increases from 0 to 12 for any Re, the area under the 𝐶_𝑝 curve increases indicating larger pressure difference between the bottom and the top surfaces. Similar trend is observed for different Re with the same angle of attack. These are expected trends for any airfoil.

(a)

(b)

(c)

Figure 7 

(a) 𝐶𝑝 distribution of NACA 63-421 airfoil at Re=68,459. (b) 𝐶𝑝 distribution of NACA 63-421 airfoil at Re=479,210. (c) 𝐶𝑝 distribution of NACA 63-421 airfoil at Re=958,422.

3.2.2. Integrated Pressure Coefficient, 𝐶_𝑝

Figure 8 represents the overall integrated pressure coefficient (𝐶_𝑝) as a function of angle of attack (𝛼) of NACA 63-421 airfoil at the three different Reynolds numbers. As expected, as we increase the angle of attack, the overall pressure coefficient increases for all six airfoils. However, within the same airfoil, 𝐶_𝑝 has little change as we move from a lower Reynolds number (Re=68,459) to a higher Reynolds number (Re=958,422). The 𝐶_𝑝 of NACA 63-218 airfoil increases continuously as we increase the angle of attack which indicates that it has not reached the stall condition yet, while the 𝐶_𝑝 plot of the other airfoils starts to flatten at around 11° to 12° of angle of attack which indicates that it is close to its stall condition. In addition, NACA 63-218, NACA 63-421, NACA 64-421, and NACA 65-421 airfoils have small integrated 𝐶_𝑝 (𝐶_𝑝 around 1.3 or 1.4) at stall condition which are much smaller than FX 63137 and E 387 airfoils (𝐶_𝑝 around 1.6 or 1.8). Thus, we can conclude that the stall conditions could vary significantly between various airfoil profiles.

Figure 8 

Integrated pressure coefficient of NACA 63-421 airfoil at different Re.

3.2.3. Coefficient of Lift, 𝐶_𝐿

Figure 9 shows the plot of lift coefficient for NACA 63-421 airfoil. 𝐶_𝐿 is plotted as functions of angles of attack and Reynolds number. The general trends of all the plots are similar as expected; that is, 𝐶_𝐿 increases with increasing 𝛼 and Re. Some of the observations from the plots are as follows.(i)The variations of 𝐶_𝐿 between different Re are not significant.(ii)The differences are more significant at higher 𝛼 for NACA 65-421 and E 387 airfoil.(iii)𝐶_𝐿 for NACA 63-218 at all Re did not reach the stall conditions; that is, the stall condition will be reached at much higher than 𝛼=12^∘.(iv)Both FX 63-137 and E 387 indicate reaching stall condition at around 𝛼=12^∘.(v)Both FX 63-137 and E 387 show smaller variation with Re and reach higher values of 𝐶_𝐿 = around 1.8 and 1.6 at 𝛼=12^∘.

Figure 9 

CFD simulation results of 𝐶𝐿 versus 𝛼 of NACA 63-421 airfoil.

It is obvious from the previous observations that different airfoils behave differently with angle of attack and Reynolds numbers.

3.2.4. Coefficient of Drag, 𝐶_𝐷

Drag Coefficient as a function of angle of attack of NACA 63-421 at the three different Reynolds numbers is shown in Figure 10. As the velocity goes down from 𝑣=7 m/s (Re=479,210) to 𝑣=1 m/s (Re=68,459), the 𝐶_𝐷 curves increase drastically for all six airfoils, while if we increase the velocity from 𝑣=7 m/s to 𝑣=14 m/s (Re=958,422), the 𝐶_𝐷 curves do not change a lot. As expected, for lower Re and larger 𝛼, the higher is the Drag Coefficient. For instance, NACA 64-421 airfoil has the highest 𝐶_𝐷 (𝐶_𝐷=0.5259) at Re=68,459 and 𝛼=12. Hence there is an optimum combination of 𝛼 and Re for the maximum ratio of 𝐶_𝐿 by 𝐶_𝐷 for each airfoil. These optimum conditions are presented in the next sections.

Figure 10 

CFD simulation results of 𝐶𝐷 versus 𝛼 of NACA 63-421 airfoil.

The CFD simulation results for the 15 cases are shown in Table 1 of NACA 63-421 airfoil where the design variables and objective functions are given in normalized form.

Table 1 

Normalized design variables and objective function values from CFD simulations of NACA 63-421.

3.3. Response Surface Approximation

The CFD data of 15 cases were used to generate a response surface for each of the two objective functions for each airfoil shape. The response surfaces were fit using standard least-square regression with quadratic polynomial using JMP [15]. The following response surfaces for each of the objective functions were obtained as functions of the two design variables of NACA 63-421 airfoil:

Lift Coefficient𝐶_𝐿()+()−()−()−().=0.1306+0.1646∗Re1.2649∗𝛼0.0766∗Re∗Re0.1732∗Re∗𝛼0.3843∗𝛼∗𝛼 (5)

Drag Coefficient 𝐶_𝐷()+()+()−()+().=0.6795−1.6805∗Re0.1358∗𝛼1.1322∗Re∗Re0.2220∗Re∗𝛼0.2329∗𝛼∗𝛼 (6)

The quality of the response surface of this airfoil is shown in Table 2. The response surface for the entire objective had very high adjusted coefficient of both 𝐶_𝐿 and 𝐶_𝐷 which indicate good capabilities for this airfoil.

Table 2 

Quality parameters of response surface of NACA 63-421 airfoil.

3.4. Optimization Results

In order to obtain the Pareto optimal solutions, the two response surface equations were incorporated in the NSGA-II code with the input parameters as mentioned in the previous section. After the simulation, the code generates a file containing 100 non-dominated solutions created during the final iteration. Non-dominated values are the best values according to the desired maximization of the objective functions among the entire population. For better understanding 2D plots of 𝐶_𝐿 versus 𝐶_𝐷 and 𝐶_𝐿/𝐶_𝐷 versus 𝐶_𝐿 are depicted using only 100 nondominated solutions for NACA 63-421 airfoil in Figures 11 and 12. The 𝐶_𝐿 versus 𝐶_𝐷 plot represents 100 best combinations of 𝐶_𝐿 and 𝐶_𝐷 corresponding to the best combinations of Re and 𝛼. Therefore, the designer can choose any of these combinations and get good results. But, to be able to get one optimum result of these 100 combinations, we can plot 𝐶_𝐿/𝐶_𝐷 versus 𝐶_𝐿. For example, for NACA 63-421 airfoil, the optimum point is at 𝐶_𝐿=0.717 for 𝐶_𝐿/𝐶_𝐷 around 40. This optimum 𝐶_𝐿 at the maximum value of 𝐶_𝐿/𝐶_𝐷 corresponds to an angle of attack of around 4°.

Figure 11 

2D presentation of non-dominated solutions of 𝐶𝐿 versus 𝐶𝐷.

Figure 12 

2D presentation of non-dominated solutions of 𝐶𝐿/𝐶𝐷 versus 𝐶𝐿.

Similar optimizations were performed on all the airfoil shapes. The results obtained are as follows: for NACA 63-218 airfoil, the optimum 𝐶_𝐿=0.54 and 𝐶_𝐷=0.026 at 𝛼=4.58 and Re=777,284; for E387 airfoil, the best combination is at 𝛼=2.90 and Re=721,398 which give 𝐶_𝐿=0.7 and 𝐶_𝐷=0.129; for FX 63137 airfoil, the optimum 𝐶_𝐿=0.976 and 𝐶_𝐷=0.024 at 𝛼=2.10 and Re=730,996; for NACA 63-421 airfoil, 𝐶_𝐿=0.717 and 𝐶_𝐷=0.018 for Re=770,937 and 𝛼=4.06; for NACA 64-421 airfoil, at 𝛼=3.12 and Re=725,083 for 𝐶_𝐿=0.68 and 𝐶_𝐷=0.02; finally for NACA 65-421 airfoil, the optimum 𝐶_𝐿=0.775 and 𝐶_𝐷=0.023 corresponding to 𝛼=3.07 and Re=759,718. The optimum 𝐶_𝐿 corresponding to maximum 𝐶_𝐿/𝐶_𝐷 as a function of angle of attack for all the airfoils was plotted in Figure 13. From this figure, the designer can choose the airfoil shape corresponding to the angle of attack dictated by the twist angle he/she is using. For example, if the twist angle is 4° at any blade section, the designer can use NACA 63-421 airfoil shape for that section of the blade in order to get the optimum 𝐶_𝐿.

Figure 13 

Optimum 𝐶𝐿 versus 𝛼 for all airfoils.

4. Conclusions

The Pareto optimal front in multiobjective optimization problem is useful to visualize the tradeoffs among different objectives. In addition to identify compromise solutions, this also helps the designer set realistic design goals. The goal of the current research is focused on the determination and optimization of wind turbine airfoil performance. For this purpose, six different two-dimensional airfoil profiles were studied over two important design variables. The NSGA-II approach of optimization and response surface methodology has been used to generate Pareto optimal front. The optimum 𝐶_𝐿, corresponding to the best combination of 𝛼 and Re, is different from one airfoil shape to another. In summary, the proposed systematic approach is very useful for optimization of designs with many variables and objectives and can be practically applied in many disciplines.

Acknowledgment

This work is supported by the National Science Foundation (NSF) through the Center for Energy and Environmental Sustainability (CEES), a CREST Center, award no. 1036593.