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 ( [subscript]CL[/subscript] ) and Drag Coefficient ( [subscript]CD[/subscript] ) 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 k-[straight epsilon] 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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Ziaul Huque 1,2 and Ghizlane Zemmouri 1,2 and Donald Harby 1,2 and Raghava Kommalapati 2,3

Recommended by Mahesh T. Dhotre

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

Received 9 January 2012; Revised 13 February 2012; Accepted 13 February 2012

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 ( _CL and _CD ) 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 omitted; refer to PDF]

Figure 2: Mesh generated around the airfoil.

[figure omitted; refer to PDF]

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.

[figure omitted; refer to PDF]

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.

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 _σa is minimized and quality parameter ^Radj2 is made as close as possible to one [12]. The _σa and ^Radj2 are defined as follows.

Let N be the number of data points and let _Np be the number of coefficients, and error _ei at any point i is defined: [figure omitted; refer to PDF] where ^fip is the actual value of the function at the design point and ^fip is the predicted value. Hence, [figure omitted; refer to PDF] [figure omitted; refer to PDF] where [figure omitted; refer to PDF] The number of data N has to be greater than the number of coefficients _Np so that the denominator of (2) is always positive and well posed. Since ^Radj2 needs to be as close as possible to 1 to represent a good fit, the terms in the numerator of (3) ^{(_σa )2} (_NP -1) should be less than or equal to the denominator ^{∑1N(_yi -y...)2} so that ^Radj2 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 _Pt of size N . From the parent set it then creates an offspring population _Qt of size N . The NSGA-II algorithm, instead of finding the non-dominated front of _Qt only, uses the combined population of _Pt and _Qt to form _Rt of size 2 N . Then, a nondominated sorting is used to classify the entire population _Rt . Although this requires more effort compared to performing a non-dominated sorting on _Qt 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, _Fi , 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 _Rt is 2 N , 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.

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 _CL Versus α of CFD simulation results with NREL experimental data.

[figure omitted; refer to PDF]

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 omitted; refer to PDF]

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

[figure omitted; refer to PDF]

3.2. Performance Trend

3.2.1. _Cp Distribution

Figures 7(a)-7(c) show the _Cp 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 _Cp distribution at the top surface of the airfoil, indicating lower pressure, and the top line represents the _Cp 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 _Cp 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) _Cp distribution of NACA 63-421 airfoil at Re=68,459 . (b) _Cp distribution of NACA 63-421 airfoil at Re=479,210 . (c) _Cp distribution of NACA 63-421 airfoil at Re=958,422 .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

3.2.2. Integrated Pressure Coefficient, _Cp

Figure 8 represents the overall integrated pressure coefficient ( _Cp ) 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, _Cp has little change as we move from a lower Reynolds number ( Re=68,459 ) to a higher Reynolds number ( Re=958,422 ). The _Cp 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 _Cp 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 _Cp ( _Cp around 1.3 or 1.4) at stall condition which are much smaller than FX 63137 and E 387 airfoils ( _Cp 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.

[figure omitted; refer to PDF]

3.2.3. Coefficient of Lift, _CL

Figure 9 shows the plot of lift coefficient for NACA 63-421 airfoil. _CL is plotted as functions of angles of attack and Reynolds number. The general trends of all the plots are similar as expected; that is, _CL increases with increasing α and Re. Some of the observations from the plots are as follows.

(i) The variations of _CL between different Re are not significant.

(ii) The differences are more significant at higher α for NACA 65-421 and E 387 airfoil.

(iii): _CL 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 _CL = around 1.8 and 1.6 at α=^12° .

Figure 9: CFD simulation results of _CL versus α of NACA 63-421 airfoil.

[figure omitted; refer to PDF]

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, _CD

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 _CD 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 _CD 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 _CD ( _CD =0.5259 ) at Re=68,459 and α=12 . Hence there is an optimum combination of α and Re for the maximum ratio of _CL by _CD for each airfoil. These optimum conditions are presented in the next sections.

Figure 10: CFD simulation results of _CD versus α of NACA 63-421 airfoil.

[figure omitted; refer to PDF]

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.

CFD normalized results
Re	A	CL	CD
0.499999	0	0.208173299	0.152746
0.499999	0.25	0.467045629	0.167201
0.499999	0.5	0.682900599	0.204313
0.499999	0.75	0.868028883	0.25002
0.499999	1	0.967890613	0.335964
1	0	0.227377477	0.098445
1	0.25	0.450145952	0.120713
1	0.5	0.657935167	0.16095
1	0.75	0.828468275	0.217986
1	1	0.937932094	0.279319
0.071428	0	0.145183592	0.597703
0.071428	0.25	0.372561069	0.630909
0.071428	0.5	0.718397603	0.688374
0.071428	0.75	0.877093255	0.822994
0.071428	1	1	1

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 [figure omitted; refer to PDF]

Drag Coefficient [figure omitted; refer to PDF]

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 _CL and _CD which indicate good capabilities for this airfoil.

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

Observations	CL	CD
R2	0.994164	0.992416
R2 Adj	0.990921	0.988203
Root Mean Square Error	0.028065	0.031312
Mean of Response	0.627276	0.381842

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 _CL versus _CD and _CL /_CD versus _CL are depicted using only 100 nondominated solutions for NACA 63-421 airfoil in Figures 11 and 12. The _CL versus _CD plot represents 100 best combinations of _CL and _CD 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 _CL /_CD versus _CL . For example, for NACA 63-421 airfoil, the optimum point is at _CL =0.717 for _CL /_CD around 40. This optimum _CL at the maximum value of _CL /_CD corresponds to an angle of attack of around 4°.

Figure 11: 2D presentation of non-dominated solutions of _CL versus _CD .

[figure omitted; refer to PDF]

Figure 12: 2D presentation of non-dominated solutions of _CL /_CD versus _CL .

[figure omitted; refer to PDF]

Similar optimizations were performed on all the airfoil shapes. The results obtained are as follows: for NACA 63-218 airfoil, the optimum _CL =0.54 and _CD =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 _CL =0.7 and _CD =0.129 ; for FX 63137 airfoil, the optimum _CL =0.976 and _CD =0.024 at α=2.10 and Re=730,996 ; for NACA 63-421 airfoil, _CL =0.717 and _CD =0.018 for Re=770,937 and α=4.06 ; for NACA 64-421 airfoil, at α=3.12 and Re=725,083 for _CL =0.68 and _CD =0.02 ; finally for NACA 65-421 airfoil, the optimum _CL =0.775 and _CD =0.023 corresponding to α=3.07 and Re=759,718 . The optimum _CL corresponding to maximum _CL /_CD 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 _CL .

Figure 13: Optimum _CL versus α for all airfoils.

[figure omitted; refer to PDF]

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

References

[1] US Government, Department of Energy, US Department of Energy on Green House Gases 2009, http://www.eia.doe.gov/bookshelf/brochures/greenhouse/Chapter1.htm

[2] Department of Energy, US Fossil Energy 2009, http://www.fossil.energy.gov/index.html

[3] C. Kong, J. Bang, Y. Sugiyama, "Structural investigation of composite wind turbine blade considering various load cases and fatigue life," Energy , vol. 30, no. 11-12, pp. 2101-2114, 2005.

[4] H. Glauert The Elements of Aerofoil and Airscrew Theory , of Cambridge Science Classic Series, Cambridge University Press, 1993., 2nd.

[5] R. E. Wilson, P. B. S. Lissaman, "Applied aerodynamics of wind power machines, PB238595," Report , no. NSF-RA-N-74-113, NTIS, Springfield, Va, USA, 1974.

[6] H. Glauert, "Airplane propellers," Aerodynamic Theory , vol. 4, Dover, 1963.

[7] M. Hand, D. Simms, L. Fingersh, D. Jager, J. Cotrell, S. Schreck, S. Larwood, "Unsteady aerodynamics experiment phase vi: wind tunnel test configurations and available data campaigns,", no. NREL/TP-500-29955, NREL, 2001.

[8] M. Alberdi, M. Amundarain, A. J. Garrido, I. Garrido, O. Casquero, M. De La Sen, "Complementary control of oscillating water column-based wave energy conversion plants to improve the instantaneous power output," IEEE Transactions on Energy Conversion , vol. 26, no. 4, pp. 1021-1032, 2011., [email protected]

[9] A. E. Marzani, F. C Ruiz, M. A. Rodriguez, M. T. Parra Santos, "Numerical modelling in wave energy conversion systems," Energy , vol. 33, no. 8, pp. 1246-1253, 2008.

[10] "UIUC Airfoil Coordinates Database," http://www.ae.illinois.edu/m-selig/ads/coord_database.html

[11] T. Goel, R. Vaidyanathan, R. T. Haftka, W. Shyy, N. V. Queipo, K. Tucker, "Response surface approximation of pareto optimal front in multi-objective optimization," in Proceedings of the 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, pp. 2230-2245, Albany, NY, USA, September 2004.

[12] R. H. Myers, D. C. Montgomery Response Surface Methodology-Process and Product Optimization Using Designed Experiment , Wiley-Interscience, 1995.

[13] K. Deb Multi-Objective Optimization using Evolutionary Algorithms , John Wiley and Sons, New York, NY, USA, 2004.

[14] M. S. Selig, B. D. McGranahan Wind Tunnel Aerodynamic Tests of Six Airfoils for Use on Small Wind Turbines , University of Illinois at Urbana Champaign, Urbana, Ill, USA, 2004.

[15] JMP, "The statistical discovery software," Version 5, SAS Institute Inc., Cary, NC, USA, 1989-2002.