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 (, 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 (, , and ) 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 m/s (for ), m/s (for ), and m/s (for ). 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 is made as close as possible to one [12]. The and 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:
where is the actual value of the function at the design point and is the predicted value. Hence,
where
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 needs to be as close as possible to 1 to represent a good fit, the terms in the numerator of (3) should be less than or equal to the denominator so that 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 (Pcross): 1.00,(iv)Distribution parameter (for crossover): 20,(v)Mutation probability (Pmut): 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 of size . From the parent set it then creates an offspring population of size . The NSGA-II algorithm, instead of finding the non-dominated front of only, uses the combined population of and to form of size 2. Then, a nondominated sorting is used to classify the entire population . Although this requires more effort compared to performing a non-dominated sorting on 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 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 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 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.

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.

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 () to a higher Reynolds number (). 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.

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 .(iv)Both FX 63-137 and E 387 indicate reaching stall condition at around .(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 .

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 m/s () to m/s (), the curves increase drastically for all six airfoils, while if we increase the velocity from m/s to m/s (), 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 () at and . 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.

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

Drag Coefficient

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.

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 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 and at and ; for E387 airfoil, the best combination is at and which give and ; for FX 63137 airfoil, the optimum and at and ; for NACA 63-421 airfoil, and for and ; for NACA 64-421 airfoil, at and for and ; finally for NACA 65-421 airfoil, the optimum and corresponding to and . 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.