Abstract
The goal of this paper is to employ a multiobjective genetic algorithm (MOGA) to optimize the shape of a well-known wind turbine airfoil S809 to improve its lift and drag characteristics, in particular to achieve two objectives, that is, to increase its lift and its lift to drag ratio. The commercially available software FLUENT is employed to calculate the flow field on an adaptive structured mesh using the Reynolds-Averaged Navier-Stokes (RANS) equations in conjunction with a two-equation SST turbulence model. The results show significant improvement in both lift coefficient and lift to drag ratio of the optimized airfoil compared to the original S809 airfoil. In addition, MOGA results are in close agreement with those obtained by the adjoint-based optimization technique.
1. Introduction
With recent emphasis on emission-free renewable energy, wind energy has taken a center stage in recent years with exponential growth in deployment of wind-turbines worldwide. Among wind-turbines, horizontal-axis-wind-turbines (HAWTs) are mostly deployed for power generation in Megawatt range. It is well established that the power generated by a HAWT is a function of the number of blades, the of the blade airfoil section, and the tip speed ratio λ (= rotational speed of the blade at tip/wind speed in free stream). Thus, one of the goals of the design of a wind turbine blade is to maximize its . As a result, there has been significant effort devoted in recent years to shape optimization of wind turbine blade to achieve high . In last two decades, aerodynamic shape optimization has become an important tool in aircraft design [1–4]. The focus of this paper is on the aerodynamic shape optimization of airfoil sections used in wind turbine blades since they affect their aerodynamic performance which in turn influences the amount of power a wind turbine can generate [5]. In modern wind-turbines, thick airfoils such as NACA-63XXX and NACA-64XXX are frequently used; however, new airfoil families are increasingly being developed because of multiple requirements of aerodynamics performance at rated power conditions and off-rated power conditions as well as strong structural properties [6]. National Renewable Energy Laboratory (NREL) has developed a family of airfoils for HAWT applications [7] since 1984.
The present paper focuses on the optimization of most well-known NREL airfoil, known as the S809 airfoil. This airfoil is 21% thick laminar flow airfoil whose design and experimental results are given in [8]. NREL Phase II, Phase III, and Phase VI HAWT blades are composed of S809 airfoil from root to tip [9]. Under class 3 to 4 wind conditions, S809 is subjected to low Mach number (almost incompressible) flow with Reynolds numbers in the range of one to two million. Laminar separation can occur on the suction surface for angles of attack ranging from zero to 5.13 degrees. Turbulent trailing edge separation occurs when angle of attack increases [10].
This paper presents shape optimization of S809 airfoil using a multiobjective genetic algorithm (MOGA). The commercially available software FLUENT is used for calculation of the flow field using the Reynolds-Averaged Navier-Stokes (RANS) equations in conjunction with a two-equation SST turbulence model. Using MOGA, globally optimal S809 airfoil shape is obtained which maximizes both and for a given wind speed, rotational speed, and pitch setting. The results show that the aerodynamics characteristics of optimized S809 are significantly improved. For the purpose of validation, the results are compared with those obtained by Ritlop and Nadarajah [10] using the adjoint equation-based optimization method.
2. Brief Description of Genetic Algorithm and Airfoil Parameterization
2.1. Single Objective Genetic Algorithm (SOGA)
Genetic algorithms are a class of stochastic optimization algorithms inspired by the biological evolution. In GA, a set or generation of input vectors, called individuals, is iterated over, successively combining traits (aspects) of the best individuals until a convergence in the objective values (namely or ) is achieved. In the context of the present paper, the individuals are airfoils which are generated in each generation of GA. These airfoils in a given generation are randomly generated with some constraints so that their thickness, camber, and other geometric properties do not vary a great deal from the original S809 airfoil. The population of airfoils in a given generation is selected by the user; a population of 20 airfoils is selected in this work. Execution of GA employs the following steps [11, 12].(1)Initialization: randomly create individuals (airfoils).(2)Evaluation: evaluate the fitness of each individual (airfoil) by calculating the fitness or objective function (e.g., lift coefficient or lift to drag ratio).(3)Natural selection: remove a subset of the individuals. Often the individuals that have the lowest fitness value are removed; although culling, the removing of those individuals with similar fitness is sometimes performed.(4)Reproduction: pick pairs of individuals to produce an offspring. This is often done by roulette wheel sampling; that is, the probability of selecting some individual for reproduction is given by
A crossover function is then performed to produce the offspring. Generally, crossover is implemented by choosing a crossover point on each individual and swapping alleles or vector elements at this point as illustrated in Figure 1.(5)Mutation: randomly alter some small percentage of the population.(6)Check for convergence: if the solution has converged, return the best individual observed. If the solution has not yet converged, label the new generation as the current generation and go to step 2. Convergence is achieved after a certain number of generations when the fitness value does not change for a number of consecutive generations (generally 3 to 5).
2.2. Multiobjective Genetic Algorithm (MOGA)
For many design problems, it is desirable to achieve, if possible, simultaneous optimization of multiple objectives [13]. These objectives, however, are usually conflicting, preventing simultaneous optimization of each objective [14]. Therefore, instead of searching for a single optimal solution, a multiobjective genetic algorithm (MOGA) is necessary to find a set of optimal solutions (generally known as the set of Pareto-optimal solutions). A solution belongs to the Pareto set if there is no other solution that can improve at least one of the objectives without degrading any other objective. For Pareto-optimal solutions, any individual inside the Pareto set dominates any individual outside the set while any individual in the set is not dominated by another individual in the solution set. An individual dominating another individual in the Pareto sense means that this individual is either better than or the same as the other individual for all objectives and at least there is one objective or fitness function for which this individual is strictly better than the other individual. A solution is said to be Pareto optimal if and only if there does not exist another solution that dominates it. The set of all Pareto- optimal solutions is called the Pareto-optimal set. The MOGA used to find the Pareto-optimal solutions to the airfoil optimization problem in this study is widely known as NSGA-II. It has the following three features: (1) it uses an elitist principle, (2) it uses an explicit diversity preserving mechanism, and (3) it emphasizes nondominated solutions in a population [15]. The basic philosophy behind the NSGA-II algorithm is as follows.
First, the population is initialized based on the constraints. The population is then sorted based on nondomination into various fronts, the first front being completely the nondominant set in the current population and the second front being dominated by the individuals in the first front only and so on for the other fronts. Individuals (airfoils) in each front are assigned rank (fitness) values based on the front to which they belong. Individuals in first front are assigned a given fitness value and individuals in the second front are assigned another fitness value and so on. In addition to the fitness value, a new parameter called crowding distance is calculated for each individual. The crowding distance is a measure of how close an individual is to its neighbors. Large average crowding distance usually results in better diversity in the population. Parents are selected from the population by using binary tournament selection based on the rank and the crowding distance. An individual selected in the rank is lesser than the other if the crowding distance is greater than the other. The selected population generates offsprings from crossover and mutation operators. The population with the current offsprings is sorted again based on nondomination and only the best individuals are selected, where is the population size. The selection is based on rank and on the crowding distance on the last front. The implementation procedure of NSGA-II is briefly as follows [16].(1)At 0th generation: a random parent population of airfoils of size is created based on constraints; it is sorted based on the nondomination into various fronts, the first front being the completely nondominant set in the population and the second front being dominated by the individuals in the first front only and so on for the other fronts. Once the nondominated sorting is complete, the crowding distance is assigned. Since the individuals are selected based on the rank and crowding distance, all the individuals in the population are assigned a crowding distance value. Crowding distance is assigned frontwise. Once the individuals are sorted based on nondomination and with crowding distance assigned, the selection is carried out using a crowded-comparison-operator. is then sent to the selection, recombination, and mutation operators to create offspring population of size .(2)At th generation: a combined population of size is formed and is sorted according to nondomination. Then again the individuals in are divided into the best nondominated set , the next-best nondominated set , and so on. If the size of is smaller than , all members of go to , with the remaining members chosen from , until the size of is . It is accomplished as described for the 0th generation. Then the new population is sent to selection, crossover, and mutation operators to create a new population of size .(3)Termination: the procedure is terminated when convergence criterion is met. The convergence is considered achieved when the two objective values do not change from one generation to the other; generally, the algorithm termination condition is applied when no improvements are observed after a number of consecutive generations.
The java code package utilized in this study is called jMetal. It is a Java-based framework for multiobjective optimization using metaheuristics. It is easy to use and is flexible and extensible [17].
2.3. Airfoil Parameterization
The airfoil shapes are parameterized using Bezier curves. Bezier construction is a curve fitting method for constructing free-form smooth parametric curves which are widely used in CAE design data structure modelling and computer graphics application [18]. A Bezier curve is defined by a set of Bezier control points. Each curve can be expressed as polynomial equations containing the information of Bezier control points. The number of control points required to parameterize a curve depends on the shape of the curve.
Each airfoil is divided into top and bottom boundary curves by the airfoil chord joining its leading edge and trailing edge. Considering the complexity of S809 airfoil shape, 12 control points are used for parameterization. For an airfoil curve, two points are fixed since they represent the leading and trailing edge of the airfoil. The intermediate points are allowed to move within the specified boundaries. A maximum thickness constraint of 19%–22% of chord is imposed on the airfoil. The constraints applied to the Bezier control points are shown in Table 1.
3. The Shape Optimization of NREL S809 Airfoil
This section presents the optimization process for S809 airfoil using the multiobjective genetic algorithm (MOGA). An optimization procedure is established by coupling the MOGA code coupled with the mesh generation code “ICEM” and the CFD solver “FLUENT” as shown in Figure 2.
The individuals in each generation of MOGA are represented by a set of control points, which generate the airfoil shape through the Bezier curve. The mesh around the airfoil shape is generated using the grid generation software ICEM, which is used to create a two-dimensional structured or unstructured mesh as an input to CFD solver FLUENT. FLUENT is used to calculate the flow field for given flow conditions. Using the flow field data, FLUENT calculates the lift coefficient and the drag coefficient which are used to calculate the objective values for a given airfoil shape. Using the information about the objective values for all the airfoils in a given generation, MOGA is applied to create a next generation of airfoils and the process is repeated to obtain the Pareto front following the MOGA procedure outlined above. From the Pareto front approximation, Pareto front-based solution for the best objective values is obtained. The airfoil shape that corresponds to the best objective values is the final shape of the optimized airfoil.
3.1. Implementation of MOGA
NSGA-II [16] and the jMetal [17] multiobjective GA software packages are employed. We choose 20 individuals (airfoils) for each generation. The crossover rate of 0.9 is considered. The mutation rate is determined to be 1/24. jMetal MOGA framework offers multiple operators; here, we employ the simulated binary crossover (SBX) operator and the polynomial mutation operator for crossovers and mutations, respectively. The selection employs the binary tournament operator.
The multiobjective optimization algorithm is performed with two objective functions. The first objective is to minimize 10/ and the second objective is to minimize . The goal is to find the Pareto front for these two objective functions. When the value of both objective functions does not change for a number of consecutive generations (normally 3 to 5), the solution is considered converged to a Pareto-optimal approximation of and . The airfoil shape that corresponds to the optimal objective values is the final shape of the optimized airfoil.
3.2. Mesh Generation
The commercially available software “ICEM” is used to generate a structured C-mesh around the NREL S809 airfoil. A reply file is scripted to automatically generate mesh around different airfoils in a given generation. The reply file is edited to be able to generate mesh based on different airfoil shapes. Figure 3 shows a typical C-mesh around an airfoil. In this mesh, there are 58460 quadrilateral cells. Far field boundary is set at 20 chord lengths. Adaptive meshing is employed. The mesh is transported to the CFD solver FLUENT for calculation of the flow field.
(a)
(b)
3.3. Flow Field Computations
The commercial software FLUENT is employed to calculate the lift coefficient and the drag coefficient of an individual airfoil in a generation. Because of low Mach number , the convergence is very slow. It is due to the well-known fact that the time-marching flow solvers designed for compressible flows do not perform well at very low subsonic Mach numbers due to large condition number of the associated eigenvalue matrix [19]. To avoid the slow convergence, we multiply the free stream velocity by four which increases the Mach number but still keeps it within the low subsonic flow regime. To keep the Reynolds number unchanged, chord length of the airfoil is divided by four. First we validate our calculations by comparing them with the numerical results of Ritlop and Nadarajah [10] and the experimental data of Ramsay et al. [20]. The simulations are performed at Reynolds number of 1 million and Mach number of 0.044. Nine different angles of incidence are selected for validation, which are 0.0°, 2.1°, 4.1°, 6.1°, 8.2°, 10.1°, and 11.2°. The turbulence model employed is a two-equation k-omega SST turbulence model. Tables 2 and 3 and Figures 4 and 5 show that our computed results for lift and drag coefficient agree well with the experimental values [20] as well as the computations of Ritlop and Nadarajah [10]. These comparisons validate the numerical methodology used in our calculations.
It can be seen from Figure 4 that while the two sets of computational results for lift coefficient agree quite well with each other for angle of attack up to 11 deg., the experimental results agree with the computations only up to an angle of attack of 8 deg. There is separation observed in the experiment for AOA > 8 deg. which is not found in the calculations by two different investigators. The separation near the aft region of the airfoil in the experiment for AOA > 8 deg. is the reason for discrepancy between the computed and experimental values of the lift coefficient. In Figure 5, the computed and experimental values of the drag coefficient are compared. The drag is a very difficult quantity to compute as well as to measure experimentally. Nevertheless, the two sets of computations are in acceptable agreement. The discrepancy between the computed and experimental values of the drag coefficient is again due to the reasons explained above.
A journal file is written for autorunning of FLUENT in the optimization process. Temperature and static pressure are defined at standard sea level condition and are taken as 288.16 K and 101325 Pa, respectively. Both values are quite reasonable for a wind turbine whose maximum altitude does not exceed a few hundred meters. Density is taken as and laminar viscosity is taken as kg/m·s.
Wind-turbines generate power due to rotation of the blades. Blade Element Momentum (BEM) theory is used to determine the generated power [21]. The BEM theory is based on Glauert’s propeller theory [22] which has been modified for application to wind-turbines [23]. In the present paper, we are simply interested in determining the relative velocity faced by an airfoil of the wind turbine as shown in Figure 6 from [23]. The mathematical expressions for axial velocity and blade velocity at radius are given by the axial induction factor and the radial or rotational induction factor , respectively. The relative velocity can be expressed as [23] where is the angular velocity of wind turbine blade, r is the radial position of airfoil section, and is the free stream velocity.
Using the published data in [9], the effective angles of attack and the relative velocity for free stream velocity , rotational speed rpm, and pitch setting of 5 deg. are obtained as shown in Table 4.
4. Results and Discussion
As mentioned before, in the application of multiobjective genetic algorithm (MOGA), we optimize the S809 airfoil at three locations of the blade, 23%, 58%, and 98% locations from the center of the rotor which correspond to the root, mid, and tip section of the blade, respectively. We consider the free stream wind velocity of 6.7 m/s, rotational speed of 72 rpm, and pitch setting of 5 deg. We set two objectives: minimize 10/ and . The airfoil shape that results in the best possible lowest values of both objectives from the Pareto front approximation is considered the shape of the optimized airfoil. Figures 7, 8, and 9 show the evolution process of the airfoil shape and variation in lift to drag ratio at various generations of MOGA at 23%, 58%, and 98% span location, respectively. Figure 10 shows the dominant solutions in the Pareto front for the airfoil sections at three span locations. Figures 11, 12, and 13 show the comparison of airfoil shapes and pressure coefficients between the MOGA optimized airfoils and the original airfoils at three different span locations.
(a) 23% span airfoil
(b) 58% span airfoil
(c) 98% span airfoil
(a)
(b)
(a)
(b)
(a)
(b)
A comparison between the present optimized airfoil shapes using MOGA and those obtained by Ritlop and Nadarajah [10] using the adjoint method is shown in Figure 14. Tables 5 and 6 show the comparison of and for the present MOGA optimized airfoils with those optimized by Ritlop and Nadarajah [10] using the adjoint method. For the lift coefficient, the present results compare quite well with those in [10]. For the ratio , the present results give a slightly higher value than those given in [10]. These results show that MOGA can give optimized airfoils that meet the two objectives of minimizing 10/ and ; the results are close to those obtained by the adjoint method.
(a) 23% span airfoil
(b) 58% span airfoil
(c) 98% span airfoil
5. Conclusions
In this paper, a multiobjective genetic algorithm (MOGA) has been employed to optimize the shape of a well-known wind turbine airfoil S809 to improve its lift and drag characteristics, in particular to achieve two objectives, that is, to increase its lift and its lift to drag ratio. The commercially available software FLUENT is employed to calculate the flow field on an adaptive structured mesh using the Reynolds-Averaged Navier-Stokes (RANS) equations in conjunction with a two-equation SST turbulence model. The results show significant improvement in both the lift coefficient and lift to drag ratio of the optimized airfoil compared to the original S809 airfoil. In addition, the present results are in close agreement with those obtained by the adjoint-based optimization technique [10].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.