Aerothermodynamic Optimization of Aerospace Plane Airfoil Leading Edge

AbstrAct: Aiming to mitigate the aerodynamic heating during hypersonic re-entry, the aerothermodynamic optimization of aerospace plane airfoil leading edge is conducted. Lift-todrag ratio at landing condition is taken as a constraint to ensure the landing aerodynamic performance. First, airfoil profile is parametrically described to be more advantageous during the optimization process, and the Hicks-Henne type function is improved considering its application on the airfoil leading edge. Computational Fluid Dynamics models at hypersonic as well as landing conditions are then established and discussed. Design of Experiment technique is utilized to establish the surrogate model. Afterwards, the previously mentioned surrogate model is employed in combination with the Multi-Island Genetic Algorithm to perform the optimization procedure. NACA 0012 is taken as the baseline airfoil for case study. The results show that the peak heat flux of the optimal airfoil during hypersonic flight is reduced by 7.61% at the stagnation point, while the lift-to-drag remains almost unchanged under landing condition.


INTRODUCTION
Aerospace planes (ASP) encounter severe aerodynamic heating during the hypersonic phase of atmospheric re-entry.Its reusable nature dictates that it should be able to shield the underlying structure from excessive temperatures during hypersonic flight and still have good aerodynamic performance at the landing speed.Regarding the Space Shuttle, a double-delta wing configuration was adopted to optimize the hypersonic flight as well as to obtain a good lift-to-drag ratio for landing (Launius and Jenkins 2012).
Since Computational Fluid Dynamics (CFD) plays a critical role in the aerospace industry, airfoil optimization has been widely studied during the design process of a winged vehicle.Buckley and Zingg (2013) developed a weighted-integral objective function to perform multipoint aerodynamic shape optimization in which a range of operating conditions were involved.A 2-step approach was introduced to conduct the aerodynamic and structural optimization of the adaptive wing leading edge (Sun et al. 2013).Various algorithms were employed for the aerodynamic optimization.A novel global optimization algorithm based on the particle swarm one was developed and applied to a low-velocity airfoil optimization (Yang et al. 2015).Koziel and Leifsson (2014) proposed an approach utilizing the multi-objective evolutionary algorithm together with surrogate model to obtain the Pareto front of a transonic airfoil.Li et al. (2012) developed an efficient method using the response surface model and genetic algorithm to optimize the transonic airfoil.Xia and Chen (2015) performed the aerothermodynamic optimization of a hypersonic wing profile to decrease the maximum heat flux.However, research on the hypersonic Zhou C, Wang Z, Zhi J, Kretov A aerothermodynamic optimization for ASP considering the landing aerodynamic performance, is still limited.
In this study, based on the Space Shuttle re-entry case, the aerothermodynamic optimization of an airfoil leading edge was carried out to alleviate the severe aerodynamic heating during hypersonic re-entry, while aerodynamic characteristics under landing condition were simultaneously considered.First, linear superposition method of analytic function with a modifi ed Hicks-Henne type function was used for parametric modelling.CFD models for both hypersonic and landing conditions are described here.Th en, the adopted optimization approach is presented, followed by a discussion on the optimized results.

SPACE SHUTTLE RE-ENTRY DESCRIPTION
Th e atmospheric re-entry is particularly challenging, and the Space Shuttle is designed for a predefi ned schedule to survive the extreme environment (Launius and Jenkins 2012;Powers 1986).Figure 1 illustrates the nominal re-entry fl ight corridor for Space Shuttle (Sellers 2004).It was controlled to fl y at a 40° angle of attack, producing high drag, not only to slow it down to landing speed, but also to reduce re-entry heating.Th e large amount of potential and kinetic energy is dissipated as heat as the Space Shuttle enters the atmosphere.A large detached bow shock wave carries away most of the heat, with the rest transferred to the vehicle through convection and radiation (Allen and Eggers Jr 1958;Tetzman 2010).
Th e aerodynamic heating is the severest at the stagnation point.It is assumed that the stagnation zone was in chemical equilibrium.Although the gas behind the shock is most likely in a non-equilibrium state, the approximation of chemical equilibrium boundary layer is reasonable in the stagnation zone.Th e heat fl ux at the stagnation point, which is the maximum heat fl ux of the leading-edge, can be estimated by (Bian and Zhong 1986): where: c and m are constants; ρ refers to the local air density; ρ 0 is the air density at sea level; R s represents the radius of curvature at the stagnation point; V is the vehicle's velocity; V 0 = 7.9 km/s is the fi rst cosmic velocity.
Th e peak temperature of the outer surface is always close to the radiation equilibrium temperature.Th e heat fl ux arising from the aerodynamic braking should not cause the temperature on the ASP surface (T w ) to exceed the maximum permissible values for materials placed on the outer surface.As for the stagnation point, where: σ = 5.67 × 10 −8 W/(m 2 K 4 ) stands for Stefan-Boltzmann constant; ε is the emissivity of the surface, depending on the material processing and surface temperature.
Aft er a series of steep S-shaped banking turns, the vehicle lowered its nose into a shallow dive and began its approach to the landing site.Th en the nose was pulled up to fi nally slow down the vehicle to approximately 100 m/s at touch-down.Unlike commercial airliners, the Space Shuttle glides to runway with no power and has relative low lift -to-drag ratio, so it needs a big angle of attack to maintain the longitudinal fl ying quality (Powers 1986).

AIRFOIL PARAMETERIZATION
Th e airfoil profi le is regenerated through altering the value of the control points during the optimal design of 2-D airfoil, while parametric description tends to be more advantageous.Th is paper uses the linear superposition method of analytic function to fi t the airfoil profi le, which is defi ned from the baseline one, type function and corresponding coeffi cients, as: Aerothermodynamic Optimization of Aerospace Plane Airfoil Leading Edge where: y b ( x) represents the baseline airfoil profile; n and a k are the number of control points and the coefficients, respectively; is the perturbation of the baseline airfoil.
In this paper, NACA 0012 is chosen as the baseline airfoil, which is similar to the airfoil for Space Shuttle Orbiter Columbia NACA 0012-64 (Rochelle et al. 1973).These airfoils have the same leading-edge radius, and the only difference is the location of the maximum thickness.NACA 0012 airfoil can be described as:

NUMERICAL MODEL computAtionAl Fluid dynAmics model description
The flow around an airfoil is numerically simulated by solving compressible 2-D Navier-Stokes equations.The airfoil chord length is set as 5 m.C-type structured grids around the airfoil are generated by using the commercial software CFD-GEOM.Grid convergence studies are conducted, and approximately 4 × 10 5 cells are distributed in the domain.A close view of the mesh distribution is illustrated in Fig. 3.Then, the commercially available CFD-FASTRAN is employed to calculate both the heat flux at the hypersonic condition and the aerodynamic coefficients at the landing condition.
As for the hypersonic case in this paper, laminar flow model is adopted since the air is thin and the Reynolds number is small at that altitude.Radiative wall boundary condition is used while the emissivity of the airfoil wall is assumed to be 0.8.It allows for radiation heat flux at the wall according to the Stefan-Boltzmann Law (see Eq. 2).Thus, a balance is formed for heat flux at the wall between conduction to the wall and radiation from it.Regarding the landing condition, k-ε turbulent model with wall function is used to solve the problem.
Since only the leading-edge is considered in this study, a modified Hicks-Henne type function (Hicks and Henne 1978;Zhou et al. 2014) is employed to control the first quarter of the airfoil profile, which is expressed as: where: e(k) = 1n 0.5 / 1n x k , 0 ≤ x k ≤ 1, x k = 4 x, x k (k = 2, 3, 4, 5).In this paper, coefficients a 1 ~ a 5 and a 6 ~ a 10 are used to control the upper and lower surfaces of the airfoil, respectively; x k (k = 2, 3, 4, 5) are set to be [0.04,0.1, 0.3, 0.6].The constraint a 1 = a 6 is applied to maintain the continuity of the airfoil leading edge.The previously mentioned modified type function with parameters setting is illustrated in Fig. 2.

numericAl model VAlidAtion Hypersonic Condition
According to a typical re-entry trajectory of Space Shuttle (Sellers 2004), its altitude and velocity profile is shown in Fig. 4. A series of hypersonic flow simulations were conducted from 7,300 to 3,050 m/s through the re-entry stage, where the angle of attack was maintained at 40°.The heat flux variation with velocity at the stagnation point is shown in Fig. 5. Normally, the velocity for the maximum heat flux is about 80 -85% of the re-entry velocity (Bian and Zhong 1986;Sellers 2004).In this study, the first cosmic velocity is taken as the re-entry velocity.As shown in Fig. 5, the maximum heat flux occurs when the velocity decreases to approximately 6,700 m/s, i.e., (4) Zhou C, Wang Z, Zhi J, Kretov A roughly 84.8% of the re-entry velocity.Thus, the CFD model used to calculate the heating results at the stagnation point under hypersonic condition is considered to be viable.
data (Ladson 1988) for NACA 0012 at high Reynolds numbers show that the airfoil will normally stall around α = 16°, which is consistent with the numerical results obtained.The comparisons indicate that the numerical model is suitable for the landing condition problem.

Landing Condition
Simulations of aerodynamic characteristics of the NACA 0012 airfoil were carried out at Re = 3.4 × 10 7 using the previously described computational model.Angles of attack α ranging from 0° to 16° were considered.Figure 6 shows the lift coefficient variation with angle of attack.The slope of the lift curve for an airfoil at high Re can be estimated by the empirical formula (Lu 2009): where:c stands for relative thickness of the airfoil.
According to Fig. 6 and Eq. 6, the relative error of the lift curve slope between the CFD results and the empirical formula result is calculated to be less than 2%.The numerical result is very close to the theoretical one.In addition, experimental

OPTIMIZATION DESIGN optimizAtion ApproAch
The whole optimization process is shown in Fig. 7.All parts were integrated using the Isight framework (Dassault Systèmes Simulia Corp. 2012).Latin Hypercube Sampling (LHS; McKay   et al. 1979) was adopted for the design of experiment (DOE).
Variables are normally referred to as factors in a DOE study, while the values are known as levels.With the LHS technique, the design space for each factor is uniformly divided, and then these levels are randomly combined to specify sample points defining the design matrix.It provides an efficient method for generating random sample points, which are uniformly distributed over the entire design space.For each sample, MATLAB® was used to automatically generate the corresponding airfoil profile database.Then, CFD-GEOM and CFD-FASTRAN were employed for meshing and flow field calculation, respectively; hypersonic heating environment, as well as landing performance, were also obtained.Afterwards, surrogate models (Liem et al. 2015) are established based on DOE results.Specifically in this study, Response Surface Method (RSM; Park et al. 2009) surrogate models were used.The RSM is a statistical technique to explore the relation between design variables and responses.Low-order polynomials are usually applied to approximate the response of an actual analysis.A number of simulations, accomplished at the previous DOE stage, are required initially to construct a model.Then it can be used in optimization with a small computational cost, as only polynomial calculation is involved.For the current optimization problem, quadratic polynomial functions were adopted.Another set of random points in the design space were chosen to check these models.Surrogate models were continuously updated with additional sample points until the accuracy requirement was satisfied.Then, these RSM models were used to replace the numerical ones in the following optimization process.
During optimization, the Multi-Island Genetic Algorithm (MIGA; Wang et al. 2015) was employed.Genetic algorithms (GA) are widely used due to their advantage to treat complex non-linear optimizations.MIGA, a further development of GA, divides each population of individuals into several subpopulations called islands, and traditional genetic operations are performed on each island separately.Several individuals are then selected from each island and migrated to different ones periodically.The migration operation maintains the diversity of probable solutions and prevents the premature phenomena.

optimizAtion results
During this optimization study, peak heat flux at the stagnation point under the hypersonic condition was regarded as the objective function, while the lift-to-drag ratio at landing condition was treated as the constraint.Coefficients of control points in Eq. 3 were taken as design variables.The optimization problem is described as: where: K and K 0 refer to the lift-to-drag ratio of the optimal and baseline airfoil, respectively.
According to the description in the sections "Space Shuttle Re-Entry Description" and "Numerical Model", 2 typical flight conditions were considered (Table 1).
Isight was used to integrate MATLAB® code, CFD-GEOM as well as CFD-FASTRAN to conduct DOE and the optimization process.First, 150 sample points were selected using LHS to conduct CFD analyses for both hypersonic and landing conditions, and the design space is: a 1 , a 5 , a 6 , a 10 are among [−0.01,0.01], while a 2 ~ a 4 and a 7 ~ a 9 are among [−0.02,0.02].Then, RSM surrogate models were constructed for both objective and constraint functions.Another 20 random points were used to evaluate the accuracy of the surrogate model.Details are shown in Table 2, where RMSE stands for root mean square error.It is shown that the approximations for heat flux and lift-to-drag ratio are of high quality.
Afterwards, the surrogate model was used to replace the previous CFD models to carry out the optimization process.The critical parameters of MIGA are: the sub-population size is 10, the number of islands is 10, the number of generations is 30, the rate of crossover and mutation are 0.9 and 0.01, respectively, the rate of migration is 0.1, and the migration interval is 5.
The optimal results are shown in Table 3, where C d represents the drag coefficient.The characteristics of both baseline and (7) * Corresponding author: zhijin@nuaa.edu.cn(Zhijin Wang). ( optimal airfoils are presented.Compared with the baseline airfoil results, the optimal one has a less severe aerodynamic heating environment at the stagnation point under hypersonic condition and maintains the lift-to-drag ratio at the same level when landing.Specifically, the peak heat flux is reduced by about 7.61%.
The normalized leading-edge profiles of the baseline and optimal airfoils are illustrated in Fig. 8, where c represents the airfoil chord length.The optimal airfoil is flatter around the stagnation point.Specifically, the radius of curvature at the stagnation point is 0.758 m for the optimal case, while it is 0.690 m for the baseline one.The result is consistent with that of Eq. 1, i.e., the heat flux at the stagnation point is inversely proportional to the square root of the nose radius of the leading edge.
Figure 9 shows the heat flux distribution of the upper and lower surfaces of both airfoils.The maximum heat flux is reduced for the optimal case.
The optimized variables were also input to perform the CFD analysis.As shown in Table 3, the relative error between the RSM and CFD results is very small, which indicates that the surrogate model has a fairly good accuracy.The heat flux contours of the baseline and optimal airfoils are shown in Fig. 10, where the stagnation point positions are also marked.

CONCLUSION
An aerothermodynamic optimization procedure considering the landing aerodynamic performance has been developed for NACA 0012 airfoil.In the optimization study, a modified Hicks-Henne type function is first adopted to parametrically describe the airfoil leading edge.CFD models are then established and further validated to simulate the hypersonic and landing problem.An optimization approach composed of DOE, RSM and MIGA is used to obtain the optimal airfoil.It is found that the surrogate model results agree well with the CFD ones.The optimal airfoil has a lower peak heat flux at the stagnation point compared with the baseline one.Meanwhile, the lift-to-drag ratio at landing condition is nearly the same as that of the baseline airfoil.

Figure 3 .
Figure 3. Close view of mesh around the airfoil.

Figure 6 .
Figure 6.Lift coefficient variation with angle of attack.

Figure 4 .
Figure 4. Space Shuttle's altitude versus velocity for a typical re-entry.

Figure 9 .
Figure 9. Heat flux distribution of baseline and optimal airfoils.

Figure 8 .
Figure 8. Shape of baseline and optimal airfoils.

Table 1 .
Aerothermodynamic Optimization of Aerospace Plane Airfoil Leading Edge Flight conditions used for optimization.

Table 2 .
Evaluation of the surrogate model.