Aerodynamic Shape Optimization of a Missile Using a Multiobjective Genetic Algorithm

The aim of this paper is to demonstrate the effects of the shape optimization on the missile performance at supersonic speeds. The N1G missile model shape variation, which decreased its aerodynamic drag and increased its aerodynamic lift at supersonic flow under determined constraints, was numerically investigated. Missile geometry was selected from a literature study for optimization in terms of aerodynamics. Missile aerodynamic coefficient prediction was performed to verify and compare with existing experimental results at supersonic Mach numbers using SST k-omega, realizable k-epsilon, and Spalart-Allmaras turbulence models. In the optimization process, the missile body and fin design parameters need to be estimated to design optimum missile geometry. Lift and drag coefficients were considered objective function. Input and output parameters were collected to obtain design points. Multiobjective Genetic Algorithm (MOGA) was used to optimize missile geometry. The front part of the body, the main body, and tailfins were improved to find an optimum missile model at supersonic speeds. The optimization results showed that a lift-to-drag coefficient ratio, which determines the performance of a missile, was improved about 11-17 percent at supersonic Mach numbers.


Introduction
Although there are many variants of missiles specialized for different purposes from short-range cruise types to intercontinental ballistic ones, the flight performance criteria used to measure their effectiveness are common: range, speed, and maneuverability. The primary factors affecting the range and speed of a missile can be specified as environmental conditions, launch/propulsion system, and aerodynamic properties. Doubtlessly, the major goal in aerodynamic design of a missile is catching a large lift coefficient (C L ) and a small drag coefficient (C D ) or as an equivalent description a high lift-todrag ratio which is determined by the shape and size of the nose and body, as well as the shape, size, and location of the canard, wing, and tailfin.
Obviously, the first two approaches have high operational costs due to the consecutive reproduction requirements in actual flight tests and supersonic flow rate requirements in wind tunnel tests. Conversely, the third approach has high computational complexity due to the requirements to solve the coupled nonlinear partial differential equations resulting from the interactions of fluid with surfaces. Fortunately, advanced computer programs such as Computational Fluid Dynamics (CFD) software packages make it possible to solve fluid flow problems in complex geometries within a reasonable time by using numerical solution techniques. However, in the aviation engineering field where shaped surfaces are investigated at subsonic or supersonic flow, some assumption loses its validity due to the separation, reattachment, eddies, vortex formation, and vortex shedding phenomenon arising from adverse pressure gradient, boundary layer growth, and circulation. At that point, the choice of the turbulence model is crucial to accurately predict the boundary flow separation and shock boundary layer interaction [1]. Solving the decomposed RANS equations instead of Navier-Stokes equations reduces the computational requirements and makes it possible to simulate practical engineering flows for complex models [2]. The comparative studies from the related literature that utilize the experimental data and simulations results show that the solutions obtained by RANS models are appropriate to use in a shape optimization process of a missile as explained as follows.
Nguyen et al. [3] accomplished a two-phase study composed of optimization and validation steps to improve the aerodynamic characteristics of a missile. The experimental data and Missile DATCOM results were utilized to attain the optimum geometry, and the ANSYS Fluent model was implemented to verify the optimized configuration. Another two-step study was carried out by Vidanović et al. [4]. In the first step, CFD simulation results were obtained for the AGARD-B model with a generic-shaped nose, and the results were validated by using experimental data obtained from the model with the same structure. In the second step, the simulation model was introduced to acquire numerical predictions for a nongeneric nose configuration. Sahu [5] studied the flow control of the finned projectile to create asymmetric pressure distribution and provide aerodynamic control changing the flow field in the aft finned region of the projectile. Ocokoljić et al. [6] performed a comparative study on a guided missile that intends to modify the front part and is aimed at improving the aerodynamic attributes. The outcomes revealed that experimentally obtained aerodynamic loads are in compliance with CFD-simulated aerodynamic coefficients. Ageev and Pavlenko [7] proposed a study to decrease the aerodynamic drag at supersonic speeds for the body of revolution. In order to perform the optimization process, ANSYS DesignXplorer and IOSO optimizer were used. Drag coefficient and volume were determined as objective function and constraint, respectively. The front part of the body was improved to make it slightly blunted, and part of volume of the front part of the body was transferred to the back part to form the back face. The results showed that the aerodynamic drag was decreased about 20% using the best variation when the framework of RANS was compared with the Sears-Haack body. Vidanović et al. [8] focused on the external design of N1G and AGARD-B models at various angles of attack (AoA) and for different Mach numbers within the supersonic flow regime. To predict the lift and drag coefficients, the CFD model and experimental study were implemented.
Körpe and Kanat [9] proposed a study related to aerodynamic optimization of a UAV wing. Xfoil was used to predict drag and lift coefficients, and the results were compared with XLFR5 and ANSYS Fluent. The sequential quadratic programming was used for the optimization problem. Riddle et al. [10] solved the shape optimization problem of a missile by using a genetic algorithm method. The predictions of aerodynamic coefficients were obtained by using both AERODSN routine and Missile DATCOM software packages. Runduo and Xiaobing [11] studied the aerodynamic shape optimization subject. To solve this multiobjective optimization problem, they used nondominated sorting genetic algorithm (NSGA-II) and real-coded genetic algorithm (RGA) methods. Similarly, Yang et al. [12] made an effort on an aerodynamic shape optimization approach to modify the canards and tailfins of a guided missile by aiming to maximize the range. For this purpose, he interlinked the real-coded adaptive range genetic algorithm method with trajectory analysis. He and Agarwal [13] performed shape optimization to improve lift and drag characteristic for the NREL S809 airfoil using the genetic algorithm. The optimization results were compared with an adjoint-based optimization technique.
In this study, the N1G missile model was selected to validate aerodynamic coefficients and perform aerodynamic shape optimization. CFD solution was performed by using ANSYS Fluent to predict and validate aerodynamic coefficients at high AoA and supersonic Mach numbers. SST k -omega, realizable k-epsilon, and Spalart-Allmaras turbulence models were used to perform CFD solution at 4°and 6°AoA. Experimental results that are available, e.g., the study of Vidanović et al. [8], were verified and compared with CFD solution. It was inferred that SST k-omega gave a closer result to the experimental study. Optimization processes were performed using MOGA at 4°AoA and supersonic Mach numbers 1.4, 2, and 2.5. Lift and drag coefficients were selected as objectives that are an important factor in terms of the aerodynamic performance. The results of the optimization problem solution showed that lift and drag coefficients were satisfactorily improved when compared with the baseline geometry model.

Material and Method
2.1. Missile Geometry and Mesh Generation. The investigated missile geometry was selected as the N1G model which is a body-tail-fin configuration. The solid model of the missile geometry was formed in Design Modeler in ANSYS. The dimension of the missile is given as seen in Figure 1. A three-dimensional computational domain and mesh generation were formed using Design Modeler and Mesh in ANSYS, respectively. The three-dimensional mesh generation was required because the use of periodicity was not suitable due to the four-fin arrangement and angle of attack. The computational domain was modeled to be a cylinder with 10 body lengths upstream from the tip of the missile nose, approximately 20 model body lengths downstream from the model. The radius of the cylinder is 10 model body lengths.
The unstructured hybrid mesh was generated using the environment of Mesh in ANSYS. The mesh generation of a missile is based on capturing the shock wave since the CFD analyses have been performed at supersonic flow. Therefore, in order to resolve the boundary layer of the missile and capture a viscous sublayer to the wall function that relaxed numerical efforts, twenty-five layers of prismatic cells were generated around the missile body and fins. Tetrahedral mesh elements were generated for the remaining part of the computational domain. The cross-sectional views of the mesh generation are presented for the missile body and tail fins in Figure 2.
In order to show mesh independency, the number of the mesh elements was generated between approximately 100000 to 4.2 million. Growth rates were changed with respect to mesh density. For finer mesh, the prismatic cell generated on the boundary layer was formed with twenty-five layers and a 1.1 growth rate. The size of the mesh and its growth rate were also changed to obtain course and finer mesh. The mesh setting selected proximity and curvature to capture 2 International Journal of Aerospace Engineering the flow field around the missile body. A patch conforming method was used to generate tetrahedral mesh for the solution domain. After the CFD solution was performed for each case and finer mesh elements, it was observed that 2008827 elements were enough to obtain a good agreement with experimental data. It also ensured the convergence of the calculated aerodynamic coefficients and residuals. The missile body surface was selected to be a wall type due to no-slip condition. The computation was initially performed on a course mesh with a nondimensional distance from the body about y + ≈ 20-40. The nondimensional distance was y + ≈ 0:8-7 for finer mesh. Figure 3 represents the y + value for the missile body. In this figure, y + values are changing between 0.8 and 7 for the whole missile. The mesh density was chosen  3 International Journal of Aerospace Engineering according to convergence criteria of the computed aerodynamic coefficients and residuals, and sufficient amount of solution time and computer memory were also taken into account. Moreover, mesh independency was ensured for finer grids that are between about 2 million and 4.2 million elements [14]. Figure 4 shows mesh independency for C D values at 2.5 Mach number.

Flow Solver.
In this study, the ANSYS Fluent flow solver was used to obtain drag and lift coefficients and the flow field around the missile geometry. The Fluent code that uses a finite volume method solves conservation equations. The Navier-Stokes equations are applied to the control volume, and the equations can be simplified to provide convenience. Hence, RANS equations, which take into account the viscous effects in a simpler way, are used for the solution of the flow [15].
The density-based, implicit, steady, compressible solver was used to compute the flow field. SST k-omega, realizable k-epsilon, and Spalart-Allmaras turbulence models were used to solve missile aerodynamics since they are suitable to solve complex aerodynamic shape at supersonic Mach number and high AoA. The implicit formulation with the Roe-FDS flux type was chosen for solution, and the least square cell-based one was selected for gradient. The secondorder upwind was selected for flow. The change of aerodynamic coefficients and flow residuals was tracked to determine convergence during the computations. The convergence takes place between 1100 and 1200 iterations. When the change of the aerodynamic coefficient value became less than 1% during 100 iterations and the flow residuals reached to 10 -5 , the computation run was finished. Solution convergence is given in Figure 5 that shows residuals of solution and C D value versus iterations, respectively.
The governing equations of continuity, momentum, and energy are solved by the fluid solver. The entire system of governing equations is represented by [16]  . τ is the viscous stress tensor, and ρ, E, and p are the density, total energy, and pressure, respectively. q is the heat flux vector.    [8]. Table 1 represents the experimental results, which were performed, and Fluent results and discrepancies of the results.
The CFD solution results observed that these three turbulence models which are SST k-omega, realizable k-epsilon, and Spalart-Allmaras were in good agreement with the experimental results. However, the SST k-omega turbulence model gave reasonable results at supersonic flow and high AoA. In addition, the SST k-omega model is convenient to solve near the wall region. The discrepancy between the CFD solution and experimental results may stem from experimental error or uncertainty. A force measurement system and nonlinear behaviour of the force balance can cause this error. Asymmetries of the model or dynamic model motion can be the possible reason in force measurements. Other possible reasons may be mechanical balance design, deterioration of a strain gage, support system of the missile body, or the not fully adjusted angle of attack of the missile [4].
Drag and lift coefficients versus Mach numbers at 4°AoA are presented in Figures 6 and 7, respectively. Drag and lift coefficients versus Mach numbers are also presented at 6°A oA in Figures 8 and 9, respectively. The lift-to-drag ratio is presented to show how to change the performance of the missile at high AoA. Figure 10 represents the lift-to-drag ratio versus high AoA at supersonic Mach numbers. From Figure 10, it can be observed that the performance of the missile increases up to 12°AoA and it has approximately the same value between 22°and 30°AoA for each Mach number. However, the ratio decreases after 12°AoA because flow separation may occur on the missile body.

Multiobjective Optimization Design
Aerodynamic shape optimization is a crucial task for air vehicles. All external parts of the air vehicle are important in terms of aerodynamic performance. In optimization design, a single design objective is not frequently sufficient to face the optimization problem. Real-world optimization problems are generally suited to be modeled using multiple objectives. So, the designer prefers to solve a multiobjective optimization problem.
3.1. Objective Functions. The aerodynamic performance of the missile can be determined by computing lift and drag coefficients. In this study, the objective functions have been determined as lift and drag coefficients. The aim of the optimization study is to increase the lift coefficient and decrease the drag coefficient.

Design
Variables. The selected missile model consists of three main parts which are the nose, body, and tailfins. These parts have an effect on the aerodynamic characteristics of the missile model. The important task in achieving optimum solution is to select design variables. In our study, 12 parameters that are body length (L b ), body radius (R b ), nose curve radius (R N ), nose front radius (R F ), tailfin position (X F ), fin span (L FS ), leading-edge (γ le ) and trailing-edge (γ te ) sweepbacks, root chords (L FR ) and tip chords (L FT ) of the tailfin, and root (T r ) and tip (T t ) thickness of the tailfin were chosen to find the optimum size and shape of the missile geometry. Figure 11 represents the design variables for the missile body, nose, and tailfin.   International Journal of Aerospace Engineering

Constraints.
The optimization results should meet the requirement of the performance and provide some shape limits. The first constraint is the position of the tailfin. Summation of the tailfin length and body length should be equal or smaller than the total length.
The second constraint is the lift and drag coefficients. The lift coefficient of the optimum missile geometry should be higher than that of the baseline missile geometry while the drag coefficient should be smaller than that of the baseline missile geometry.
Subject to g i ðxÞ ≤ 0, i = 1, 2, ⋯, m where x = ½x 1 , x 2 , ⋯, x j T is the vector of design variables, f ðxÞ is the multiobjective vector, f n ðxÞ is the objective function, g i ðxÞ and h i ðxÞ are the constraints, and x i l and x i u are the lower and upper bounds of the design variables, respectively.  Figure 10: Lift-to-drag ratio versus AoA at supersonic Mach numbers.  A multiobjective optimization problem can be converted to a single-objective one. A nondimensional form must be generated by converting each cost term in the total cost expression. A multiobjective cost function can be expressed as follows: where f n ðn = 1, 2, ⋯, nÞ indicates the cost function. The value of each cost term decreases into an interval of ½0 1. In order to obtain a nondimensional form, each cost function must be divided by the referenced values or baseline missile values. The defined equation is presented as follows: The designers should cope with some tradeoffs during design. The parameters that are based on relative weighting of requirements should be optimized to decide the tradeoffs between conflicting requirements. The weightings are multiple with cost expression, and their values are determined according to the importance of the objectives [18]. The weighted and nondimensional cost function can be expressed as follows: Baseline, lower, and upper values of the design variables are given in Table 2.  Mach number

Optimization Algorithm.
In this study, the Multiobjective Genetic Algorithm (MOGA) that is a variant of the popular NSGA-II (nondominated sorting genetic algorithm-II) was used to solve the optimization problem by means of ANSYS DesignXplorer. Genetic algorithms are widely used to generate high-quality solutions for optimization, and they are a stochastic search approach (derivative methods are not applied) that is based on randomized operators. They have been applied to numerous problems in engineering and science since they are convenient to solve complex problems such as external geometry size optimization problems. In addition, GA solves optimization problems in a reasonable amount of time. This is crucial in point of waste of time and cost for industrial application. GA gives quite accurate solution results for optimization problems when previous studies are examined. Therefore, the use of GA is considered to be a more effective method and practical for this study and similar applications [18]. In an engineering problem, GA is suited to solve a multiobjective optimization problem. It is possible to find a diverse set of solution searching different regions of a solution space thanks to GA. In a crossover process, GA may utilize structures of good solutions to create nondominated solutions in undiscovered parts of the Pareto front that stands for a set of optimal solutions in a space of objective functions in multiobjective optimization problems. Hence, GA is the most popular for multiobjective optimization design [17]. MOGA uses a rank-based fitness assignment method. A fitness function is a specific type of objective function that determines the optimum of a solution in the genetic algorithm so that a particular chromosome can be sorted against other chromosomes. A design point that is determined by fitness function is used to decide whether to incorporate the design for the next iteration in the selection process. We generated a multiobjective optimization model using the genetic algorithm. The optimization processes are represented in Figure 12.

10
International Journal of Aerospace Engineering 3.6. Optimization Problem Solution. Optimization problem solution was carried out providing objective function under specified constraints and selected design variables. The calculation of the objective functions was performed at 1.4, 2, and 2.5 Mach numbers and 4°AoA using RANS with the SST k-ω turbulence model in Fluent. Optimization was carried out by means of ANSYS DesignXplorer. The geometry was formed in Design Modeler, and the specified design variables were selected to be transferred to the parameter set block of the ANSYS Workbench. The formed geometry was transferred to ANSYS Mesh, where a computational grid that is associated with geometry was generated at the different boundary regions. The generated computational grid was then transferred to ANSYS Fluent. The calculation was performed to obtain lift and drag coefficients with respect to preset initial and boundary conditions. These coefficients were transferred to the parameter set block of the ANSYS Workbench. In this way, input and output parameters were gathered together in the parameter set block. Input and output parameters are automatically transferred from the parameter set block to ANSYS DesignXplorer. There are three steps in DesignXplorer to perform the optimization process. In the first step, parametric investigation is carried out. Input data can be generated using various methods. In this study, a central composite design was applied. In the design of experiments, 281 design points were generated for 12 input parameters. This means that 281 different missile models were constructed and CFD solution was performed for each design point. In the second step, a response surface (approximating function) which depends on the obtained data was generated. Response surface construction was formed using the Kriging surface that enables the best solution of the optimization problem. After these calculations, optimization was performed using MOGA.

11
International Journal of Aerospace Engineering The configuration of the optimization process is generating 100 samples and 100 samples per iteration. Three candidate solutions were obtained. Each candidate solution was verified, and one of them was chosen and presented as an optimum result. The convergence of the optimization problem was achieved after 876 evaluations. The optimum body was then generated, and CFD solution was performed to verify optimization results.

Optimization
Results. It was observed that the lift-to-drag ratio shows satisfactory improvement after the optimization process was performed. The results of the optimization problem are given in Table 3. Design variables are also presented for the optimum missile model in Table 4. The results showed that the body radius is more effective in reducing the drag coefficients with respect to the nose curve radius and nose front radius. In addition, fin thickness, fin root and tip chord lengths, and leading-edge angle are more effective on the lift coefficient when compared with other design variables. Figures 13-15 show the pressure contours for baseline and optimum missile model comparison at 1.4, 2, and 2.5 Mach numbers and 4°AoA, respectively.
It can be observed that a low-pressure area decreases the optimum missile model when compared to the baseline geometry at 1.4, 2, and 2.5 Mach numbers. This means that the drag coefficient decreases for the optimum model. In addition, the pressure that occurs on tailfins of the leading edge of the baseline model at 1.4 Mach number is higher than that of the optimum model. This shows that tailfins of the optimum model are improved. At the same time, lower drag force and larger lift force on the tailfins represent larger control force. The higher pressure occurs on the lower area of the nose of the optimum model when compared to the baseline  13 International Journal of Aerospace Engineering model. This result shows that after optimization, the lift force occurring in lower part of the nose increases. In addition, pressure coefficient distribution is presented to understand the increase in the lift-to-drag ratio of the missile in Figures 16 and 17. The fin span and fin length increase while fin thicknesses decrease, so the fin area increase for the optimum missile model. This means that the lift force increases and the drag force decreases. It is observed that the length of the missile is affected by the drag reduction and pressure occurring on the tip of the nose for optimum shape is bigger than that for baseline shape for Mach 2 and 2.5. This is the reason for increasing the lift force on the front of the missile due to AoA. However, pressure difference on the fin is higher at Mach 1.4 when compared with Mach 2 and 2.5. It is shown that the increase in the lift-to-drag ratio is mostly due to fin design at Mach 1.4. Pressure distribution differences between the lower and upper surfaces are small since the solution process is performed at 4°AoA which is a small value for the missile; however, experimental data is available at 4°AoA for the selected missile. Consequently, the contribution of the optimum model to the lift force is more than that of the baseline model after optimization processes.

Conclusion
In this paper, missile aerodynamic analysis was performed using SST k-omega, realizable k-epsilon, and Spalart-Allmaras turbulence models at 1.4, 2, 2.5, and 4 Mach num-bers. The missile geometry model was referenced from literature which was implemented in the experimental study by Vidanović et al. [8]. It was concluded that CFD results showed good agreement with the experimental study. However, the SST k-omega turbulence model gave more accurate results since it provides superior solution performance for the thin boundary layer, recirculation, and separation. To show mesh independency, mesh was generated from course to finer mesh elements (100000 to 4.2 million). However, the results of the CFD solutions showed negligible change between 2 and 4.2 million mesh elements. In addition, the lift-to-drag ratio was also calculated to show the missile performance at high AoA. The ratio increased up to 12°AoA, and after that, it decreased since drag forces increase due to flow separation. In the optimization process, design variables were determined to perform optimization. The missile body and fin parameters need to be estimated to design the optimum missile geometry. Lift and drag coefficients were considered objective function for the optimization process. Input and output parameters were collected, and 281 design points were generated. CFD solution was then performed for each design point. The Multiobjective Genetic Algorithm (MOGA) was used to optimize missile geometry. The front part of the body, the main body, and tailfins were improved for each selected Mach number. The results of the optimization solution problem showed that the lift-to-drag ratio increases by 17.17%, 16.05%, and 11.89% for Mach numbers 1.4, 2, and 2.5, respectively. It was concluded that the body 15 International Journal of Aerospace Engineering radius is more effective in reducing the drag coefficients with respect to the nose curve radius and nose front radius. In addition, fin thickness, fin root and tip chord lengths, and span are more effective on the lift coefficient when compared with other design variables.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

16
International Journal of Aerospace Engineering