AIRFOIL SHAPE OPTIMIZATION USING BÉZIER CURVE AND GENETIC ALGORITHM

. There are different types of airfoil used in many applications such as energy production, aerospace, mixing of fluid products. Design optimization studies are still being carried out on the airfoil type structures. The airfoil section is the most important factor affecting the quality and efficiency of the performed work. The aim of this study is the optimization of the airfoil shape to generate more lift than the original airfoil shape creates. For this purpose, Bézier curves are used to generate the airfoil polar points, XFOIL is used as a flow solver and MATLAB is used to create optimization codes using the genetic algorithm. The results show that the created optimal airfoil shape produces more lift than the original airfoil shape. In this study, design optimization studies are supported by flow analysis using ANSYS Fluent.


Introduction
The pioneering research carried out in the airfoil parameterization and optimization can be summarized as follows. A hybrid evolutionary-adaptive directional local search method was used by Lim and Kim (2019) for convergence enhancement. They used XFOIL and a two-dimensional structured grid RANS solver for flow simulation. The CST method and B-spline method were used for airfoil shape parameterization. Hansen (2018) used class-shape-transformation technique for parameterization and XFOIL as flow solver. Derivative-free Covariance Matrix Adaptation Evolution Strategy (CMA-ES) was used for optimization. Jeong and Kim (2018) optimized thick airfoil using genetic algorithm and they used Akima curve fitting method for airfoil parameterization. They improved the lift-to-drag ratio by 30%~40% compared to the baseline airfoil. Lu et al. (2018) proposed a new airfoil parameterization method, called IGP method. They used thin airfoil theory, it means, the camber and the thickness are defined separately and there are two independent optimization problems. Due to independence, the computational cost decreased. Ziemkiewicz (2017) studied parametric modeling for airfoil shape description. Tandis and Assareh (2017) used genetic-based bees algorithm (GBBA). Crossover and neighborhood searching operators were used in this method and they were derived from genetic algorithm and bees algorithm, respectively. It helps to increase the convergence speed. Koreanschi et al. (2017) applied a genetic algorithm to the airfoil which has morphed upper surface. The genetic algorithm results were compared with the artificial bee colony and a gradient method. Adaptive upper surface was manufactured using carbon fiber composite materials and positioned between 20% and 65% of the airfoil chord. There are two ailerons in the study, one of them is rigid and the other one is flexible. Yang et al. (2018) presented a study based on Bézier curve parameterization and radial basis function interpolation. They used genetic algorithm for aerodynamic optimization. Sun et al. (2015) studied airfoil inverse design. For this purpose, they used PARSEC parameterization and ANN. They expect that under the given aerodynamic conditions, ANN gives the airfoil geometry. Reddy et al. (2016) studied multi-element winglets (3 elements) using radial basis function response surface approximation coupled with a genetic algorithm. They analyzed each configuration via OpenFOAM. They found that the multielement winglet concept increases the aerodynamic performance, but, if more winglet elements added in a single wing-tip, the added elements become thinner and more delicate. For this reason, both aerodynamic and structural analysis must be performed in the optimization process. Fincham and Friswell (2015) studied on a morphing system called Fishbone Active Camber (FishBAC) morphing system. This system allows the shape of the airfoil adaptation in order to produce optimum performance in a wide range of flight conditions. They used genetic algorithm for optimization and radial basis function interpolation for smooth shape changes. Della Vecchia et al. (2014) studied on the coupling PARSEC parameterization and evolutionary algorithm to optimize the airfoil shape. The procedure is to find Nash equilibrium and decide the direction of the airfoil shape modification. Mukesh et al. (2014) used PAR-SEC geometry representation method, Panel technique and genetic algorithm to optimize the NACA 2411 airfoil geometry in order to increase the lift coefficient and tested the optimized airfoil in the wind tunnel. Salunke et al. (2014) discussed Bezier curve, PARSEC techniques and combination of Bezier-PARSEC techniques which cover wide range of airfoil. Ebrahimi and Jahangirian (2014) used adaptive parameterization and genetic algorithm. Melin (2013) described the airfoils by a set of parametric Bezier curves. Timnak and Jahangirian (2018) studied on the optimization of the airfoil using PARSEC parameterization and genetic algorithm. Ribeiro et al. (2012) studied on the optimization of the airfoil using genetic algorithm and artificial neural networks together. They found that using ANN reduces computational time. Kharal and Saleem (2012) determined the airfoil geometry using Bézier-PARSEC parameterization and used three different neural networks, Feed-forward backpropagation, generalized regression, and radial basis neural network, for learning. The feed-forward backpropagation neural network yielded better results. Sripawadkul et al. (2010) compared the five airfoil shape parameterization techniques: Ferguson's curves, Hicks-Henne bump functions, B-Splines, PARSEC and Class/Shape function transformation. They considered Parsimony, Intuitiveness, Orthogonality, Completeness and Flawlessness characteristics to compare the parameterization techniques. Derksen and Rogalsky (2010) studied on Bezier-PARSEC (BP) parameterization which accelerates the convergence. Four separate Bézier curves are defined in BP parameterization. The class function/shape function transformation (CST) method was represented by Kulfan (2008). Khurana et al. (2008) studied on PARSEC parameterization coupling with the particle swarm optimization and artificial neural network. PARSEC parameters were used as inputs and coefficient of lift is used as the output of the ANN. Mengistu and Ghaly (2008) studied the genetic algorithm coupled with an ANN which uses a backpropagation algorithm. They used NURBS (Non-Uniform Rational B-splines) for airfoil shape parameterization.
The objective of this paper is to present a procedure for the optimization of the airfoils and validate the obtained results using another CFD software. For this purpose, a practical approach is presented to airfoil shape optimization.
In this study, the Bézier curve is used to make the design parametric. The utilized parameterization method has six control points and these control points are used for the optimization process. In order not to increase the computational time during the optimization study, the control points number haven't been increased above six. Two of these control points are fixed at the origin and the tip of the airfoil. The other control points are changed to define the optimized airfoil profile during the optimization study. To conduct these parameterization and optimization processes, MATLAB and XFOIL programs are utilized simultaneously. ANSYS Fluent is used for validation purposes in the study.
This study is organized as follows. In section 1, the methodology for the airfoil shape optimization using Bézier curves and genetic algorithm, and the implementation of the procedure is explained. The mesh generation process is described. In section 2, the results of CFD simulations are given and discussed. Finally, the findings of the study are given in the conclusions part.

Airfoil shape optimization
In the aerodynamic design approach, there are two kinds of design methods; conventional aerodynamic design method and inverse design method (Sun et al., 2015). In the first design method, aerodynamic features are obtained using CFD analyses or experiments for the given airfoil first, and then the airfoil geometry is optimized. These steps are repeated until the result is satisfactory. In the second design method, the airfoil geometry is obtained according to given aerodynamic features. In this study, the conventional aerodynamic design method is used. The optimization process requires many evaluations and makes changes on the airfoil profile during the run time. Evolutionary algorithms are population-based optimization algorithms and they are based on Darwin's natural evolution theory. Natural selection and survival of the fittest are the main principles of this theory. Evolutionary algorithms select the optimum solutions by imitating the natural selection processes (Messac, 2015).

Genetic algorithm
One of the evolutionary algorithms is Genetic Algorithms (GAs). A population of solutions is used by GA in the optimization process to reach the global optimum. The GA makes modifications on the population during the optimization. Considering the termination criteria, the parents are chosen by the genetic algorithm from the current population of the solution, and finally, GA reaches the optimal solution or Pareto frontier (Messac, 2015). The flowchart of the used optimization process in the present work can be seen in Figure 1. CP represents the control points determined using Bézier curves and can be seen in Figure 2 for the original and optimized airfoil. p0 is the vector representation of the original airfoil's control points. p1 contains the control points of the new population matrix. Cl indicates the lift coefficient of the airfoil. The algorithm is performed for 20 generations of 40 individuals. The crossover and mutation probabilities are set to 75% and 20%, respectively. The termination criteria of this study are the number of generations. With the help of the constraints, only appropriate solutions are evaluated by the flow solver. The objective function and constraints are defined under the operating conditions of Re = 10 6 , Ma = 0.4: : here; t, c, upper y and lower y are thickness, chord length, y coordinate vector of the airfoil's upper surface, and y coordinate vector of the airfoil's lower surface, respectively. p0 and p1 contain the optimization variables which are the control points defined by the Bézier curves.
In the present work, the idea is to reach the optimized airfoil using modification of the control points created by the Bézier curve. The abscissa and ordinate of the control points are the only design variables in the study. To achieve the optimization goal, the genetic algorithm is employed. A genetic algorithm is an iterative optimization method and a new airfoil profile is created at each iteration step if the constraints are assured. The created profile is solved using XFOIL and Cl values are saved as well. At each step, elite results are selected and kept for the next generation. This process continues until the termination criteria are satisfied.

Bézier curve
Airfoils can be described by point clouds or by mathematical functions. In the present study, the Bézier curve is used for airfoil geometry parameterization. Control points is defined using point clouds of the airfoil (Fazil & Jayakumar, 2011). The detailed definition and the formula of the Bézier curve can be found (Rogers & Adams, 1990). A parametric Bézier curve is defined by: (4) where B i is the vertices of a Bézier polygon,

( )
, n i J t is the Bernstein basis function and n is the degree of the defining Bernstein basis function. The summation of the basis functions is equal to one.
Derivatives need to be considered when combining two Bézier curves to ensure the continuity of the slope and curvature. The first and second derivatives of the Bézier curve are given by: The derivatives of the basis function are obtained by: , , First and second derivatives at the beginning and the ends of a Bézier curve (t = 0 and t = 1): Original and optimized airfoils and control points can be seen in Figure 2. There are 97 points on the airfoil surface and the number of control points is 6. Control points of the original airfoil and the optimized airfoil can be seen in Table 1. Control points are used by the genetic algorithm as design variables to create airfoil coordinates in the optimization process. XFOIL is an open-source flow solver developed by (Drela, 1989) and rapidly analyze airfoil profiles and determine the airfoil characteristics. XFOIL code combines the panel method and the integral boundary layer formulation. The computational cost of the XFOIL is small. For this reason, XFOIL is preferred for the flow solver during the optimization process (Mauclère, 2009;Morgado et al., 2016).
XFOIL provides the lift and drag coefficient for the given angle of attack (AoA), Reynolds number, and Mach number (Anitha et al., 2018). No value is written to the output file if there is no convergence when the program completes the solution procedure (Mauclère, 2009). Lift and drag coefficients are obtained for the sequence of alpha in XFOIL program using the generated airfoil profiles by the genetic algorithm. The number of panels used in the XFOIL flow solver is 160. In the process of obtaining the optimized NACA 4415 airfoil shape, MATLAB and XFOIL flow solver are used simultaneously.

Numerical procedure 1.2.1. Basic formulations
In this study, a dimensionless analysis is employed. Reynolds number is defined as: where, ρ − fluid density, ϑ − fluid velocity, c -chord length, µ − dynamic viscosity and ν − kinematic viscosity. In most studies, lift (Cl) and drag (Cd) coefficients are the two used parameters to evaluate the aerodynamic performance and defined as: where, L -lift force, D -drag force and S -airfoil area. The coefficients of lift and drag are solved using Reynolds Averaged Navier-Stokes (RANS) equations of mass and momentum which are denoted as: Conservation of mass: where, u is velocity vector, ( 2 p u −∇ + µ∇ ) represents the internal forces and F ρ represents the external forces. In the numerical analyses, y + is used as the dimensionless distance from the wall to the first node. y + values are dependent on the mesh resolution. The turbulent boundary layer's different regions are: laminar sub-layer 5 y + < ; transition or buffer layer (5 30) y + < < ; and turbulent or log-layer 30 y + > . It is important that the first cell adjacent to the wall should not be located in the buffer zone. In the viscous sublayer region 5 y + < and the loglayer region 30 y + > , u + is defined as: (20) where y is the first cell height, u + is the dimensionless velocity, u τ is the friction velocity, u is the velocity, ω τ is the wall shear stress, 0.41 κ ≈ (Von Kármán constant) and 5 B ≈ (White, 2011).

Model specifications, mesh generation and mesh independence study
The numerical results of this study are obtained based on the original and optimized NACA 4415 airfoil profiles. The chord length of the airfoil (c) is 0.1 m. In this study, the Reynolds number is selected as 10 6 and the Mach number (Ma) is taken as 0.4. As a turbulence model, Spalart-Allmaras with a standard set of coefficients is selected. The pressure-velocity coupling is handled using the SIMPLE scheme. Second-order upwind discretization scheme is utilized for all RANS equations. The residual convergence level of the simulations is between 10 -5 and 10 -6 depending on the airfoil's angle of attack. Boundary conditions have been set as pressure far-field along the computational domain and wall along the airfoil.
Mesh details in the computational domain are given in Figure 3. The computational accuracy and computational time depend on mesh quality. To generate successful mesh grids in the computational domain, it is divided into two parts. There is high-density mesh around the airfoil at the near-field area and the radius of this zone is 5c. Far-field is defined at the distance of 12.5c from the origin of the computational domain. The sizing option of the overall mesh is proximity and curvature. The edge sizing and face sizing method are used to refine the mesh quality in the near-field area and around the airfoil. Inflation layer is used along the airfoil wall with a number of layers as 25 and a growth rate as 1.2 to provide the y + < 5 requirement. y + values are about 1 for all ANSYS Fluent analyses. According to the mesh statistics on ANSYS Fluent, the numbers of elements are nearly 130,000. The boundary conditions and mesh generation methods presented in this section were used for both original and optimized airfoils in the analyses.
The results of the mesh independence study are given in Table 2. The analyses are run on the computer with Intel Core i7-6700HQ, 2.6 GHz CPU and 16 GB RAM. 2 nd grid is selected for the study in order not to increase the mesh element number and the computational time since CPU times are dependent on the number of elements. The selected number of elements is sufficiently enough to compute the airfoil flow. The y + value is below 1 for the selected grid at 0° attack angle.

Results and discussions
In this study, NACA 4415 airfoil is employed and the optimization study is conducted to increase the airfoil performance. As a parameterization technique and optimization method, Bézier curve and GA are selected, respectively. As a result of the optimization study, the camber of the optimized airfoil is increased compared to the original airfoil profile. Thereby, the lift coefficients have been improved. XFOIL software is used while processing the optimization study to find the optimal solution. XFOIL and ANSYS Fluent give information about aerodynamic coefficients of the lift and drag of the airfoil profile. The validation study of the XFOIL simulations is conducted at the various angle of attacks using the Spalart-Allmaras turbulence model in ANSYS Fluent. The significant part of the study is the investigation of the airfoil profile's aerodynamic behavior. The comparison of the pressure and velocity contours in the computational domain for both original and optimized airfoils can be seen in Figure 4. These contours are outcomes of the ANSYS Fluent. The flow analyses are performed between -5° and 20° by increasing the angle of attack by 5° in each step. The pressure and velocity con- There is a pressure difference between the upper and lower surfaces of the airfoil, so lift force occurs on the airfoil. As can be seen, the pressure difference is greater for the optimized airfoil profile, so the lift coefficients are increased compared to the original airfoil profile. Lift coefficients increase for both airfoils up to 15° of attack angle. After this degree of attack angle, the lift force attenuates and this situation is called a stall.
The stagnation points at which the flow velocity is nearly zero can be seen in the vicinity of the leading edge in the velocity contours in Figure 4. The flow velocity is greater on the upper airfoil surface than the lower airfoil surface. This situation is expected when the pressure distribution is taken into account.
Coefficient of lift with respect to the angle of attack and coefficient of drag for both original and optimized NACA 4415 airfoils can be seen in Figure 5. Considering Figure 5, it can be seen that there is a good agreement between the results of XFOIL and ANSYS Fluent software.
Both analyses show a similar lift curve. The optimized airfoil shape creates more lift force than the original airfoil shape as expected, and the enhancement can be seen from both XFOIL and ANSYS Fluent coefficient of lift results and the corresponding numerical values are presented in Table 3 and Table 4.
The ratio of lift to drag coefficients can be seen in Table 3 and Table 4 and, the change of these ratios is expressed as the percentage difference between the original and optimized airfoils for both analyses results. XFOIL results show the improvements in terms of increasing the coefficient of lift to coefficient of drag ratio, but the ratio decreases for 5° and 10°, however, this is a possible result since the study is single objective optimization study.
ANSYS Fluent results show the improvements in terms of increasing the coefficient of lift to coefficient of drag ratio, but this ratio decreases by 10° and above. However, this situation is possible because the objective of the study is to increase the lift coefficient, and the objective has been reached in the study.

Conclusions
An optimization framework is generated combining parameterization, computational fluid dynamics, and aerodynamic optimization in order to increase the aerodynamic performance of the NACA 4415 airfoil. This procedure can be applied to any airfoil profile. Geometric parameterization has a key role in the optimization process in terms of the evaluation of different airfoil profiles. The results show that the proposed optimization algorithm has enhanced the lift coefficient of the airfoil. And, this procedure can be integrated into the industrial airfoil design. Bézier curve is one of the basic parametric design methods, and it gives the opportunity to optimize the airfoil using control points. In this study, Bézier curve is used for parameterization of the airfoil geometry. And a MATLAB code is written to optimize the airfoil profile in order to increase the airfoil performance. Besides, XFOIL code is used to carry out the optimization and CFD study together. ANSYS Fluent is used to validate the obtained results created by the optimization code. CFD analyses are conducted for the original and optimized NACA 4415 airfoils at different angles of attack on ANSYS Fluent software. Mesh independence study is carried out to consider the computational time and simulation accuracy, and 5 y + < requirement is assured to obtain the accurate simulation results. The coefficients produced by the flow simulation and optimization study are in very good agreement, and even though the study is single objective optimization study, the Cl/Cd ratio has been improved for some AoA in the analyses. In conclusion, significant outcomes have been obtained in terms of increasing the lift coefficient considering the study is the single objective optimization study.
Finally, a multi-objective optimization study can be applied to the same airfoil and observed the outcomes as future work and the optimized airfoil profile can be examined in terms of determining the lift and drag coefficients experimentally.