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

1 Department of Mechanical Engineering, Prairie View A&M University, P.O. Box 519, Mail Stop 2525, Prairie View, TX 77446, USA 2 Center for Energy and Environmental Sustainability, Prairie View A&M University, P.O. Box 519, Mail Stop 2500, Prairie View, TX 77446, USA 3 Department of Civil and Environmental Engineering, Prairie View A&M University, P.O. Box 519, Mail Stop 2510, Prairie View, TX 77446, USA


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

Approaches
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 (C L and C D ) 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.

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

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.
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 v = 1 m/s (for Re = 68, 459), v = 7 m/s (for Re = 479, 210), and v = 14 m/s (for Re = 958, 422).The Realizable k-epsilon turbulence model and a second-order upwind solution method were used to get more accurate results.

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 σ a is minimized and quality parameter R 2 adj is made as close as possible to one [12].The σ a and R 2 adj are defined as follows.Let N be the number of data points and let N p be the number of coefficients, and error e i at any point i is defined: where f p i is the actual value of the function at the design point and f p i is the predicted value.Hence, where The number of data N has to be greater than the number of coefficients N p so that the denominator of (2) is always positive and well posed.Since R 2 adj needs to be as close as possible to 1 to represent a good fit, the terms in the numerator of (3) (σ a ) 2 (N P − 1) should be less than or equal to the denominator N 1 (y i − y) 2 so that R 2 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.

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 P t of size N. From the parent set it then creates an offspring population Q t of size N.The NSGA-II algorithm, instead of finding the non-dominated front of Q t only, uses the combined population of P t and Q t to form R t of size 2N.Then, a nondominated sorting is used to classify the entire population R t .Although this requires more effort compared to performing a non-dominated sorting on Q 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, F i , 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 R t is 2N, not all fronts may be accommodated in N 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.

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.
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.the airfoil indicating higher pressure.As the angle of attack increases from 0 to 12 for any Re, the area under the C p 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.(iii) C L 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 α It is obvious from the previous observations that different airfoils behave differently with angle of attack and Reynolds numbers.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.

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: International Journal of Chemical Engineering   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 C L and C D which indicate good capabilities for this airfoil.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 C L .

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 C L , 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.

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

Figure 3 :
Figure 3: The flow domain with the boundary conditions.

Figure 4 :
Figure 4: Comparison of C L Versus α of CFD simulation results with NREL experimental data.
Drag, C D .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 v = 7 m/s (Re = 479, 210) to v = 1 m/s (Re = 68, 459), the C D curves increase drastically for all six airfoils, while if we increase the velocity from v = 7 m/s to v = 14 m/s (Re = 958, 422), the C D 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 C D (C D = 0.5259) at Re = 68, 459 and α = 12.Hence there is an optimum combination of α and Re for the maximum ratio of C L by C D for each airfoil.These optimum conditions are presented in the next sections.

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 C L versus C D and C L /C D versus C L are depicted using only 100 nondominated solutions for NACA 63-421 airfoil in Figures11 and 12

Figure 12 :Figure 13 :
Figure 12: 2D presentation of non-dominated solutions of C L /C D versus C L .
Pressure Coefficient, C p .Figure8represents the overall integrated pressure coefficient (C p ) 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, C p has little change as we move from a lower Reynolds number (Re = 68, 459) to a higher Reynolds number (Re = 958, 422).The C p 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 C p 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.Figure 5: Colour contour of static pressure around NACA 64-421 airfoil at α = 6 and Re = 479, 210.
L increases with increasing α and Re.Some of the observations from the plots are as follows.(i)The variations of C L between different Re are not significant.(ii) The differences are more significant at higher α for NACA 65-421 and E 387 airfoil.

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

Table 2 :
Quality parameters of response surface of NACA 63-421 airfoil.