1 Introduction

It has been an important issue for aircraft design to reduce aerodynamic drag for a long time. This goal is particularly critical with further tension of global energy and gradual strengthening of awareness of environmental protection. The reduction in drag can significantly improve aircraft efficiency (less fuel consumption),thereby increasing the range and payload of aircraft. Laminar flow technology is the most promising avenue to reduce aircraft drag. Natural laminar flow airfoil is a direct and easy drag reduction technique. Through effective shape designing, an extensive laminar flow area can be maintained on the surface of the airfoil,thereby reducing the frictional drag of the aircraft. With the increasing use of composite material and progress of the structure and manufacture technology,the surface of aircraft becomes smoother,which is conducive to the realization of laminar flow. On the other hand,when the wing-load and the aspect-ratio of modern aircraft are larger,the flight attitude is higher,which decreases the Reynolds number (Re) based on chord and makes it easer to achieve the laminar flow.

In previous works,low-speed and relatively low Reynolds number airfoils have been widely studied. These researches have been widely used for sailplane and general aviation light plane natural laminar flow (NLF) wings. While for transonic,high Reynolds number cases,emphasis has been placed on supercritical airfoils. Almost all transonic high Reynolds number aircrafts use supercritical airfoils. The laminar-flow technology is a very efficient approach to improve the performance of supercritical airfoils. However,very little work has been done for laminar flow supercritical airfoil design. Eggleston and Poole^[1] have done some primary work on thick supercritical airfoils with NLF. However,the method is old and the precision is low to predict the transition and the shock wave. Viken andWagner^[2] have studied the limits of compressible NLF airfoils without considering the off-design points. Qiao^[3] have done some research on natural flow supercritical airfoils with low Reynolds numbers. All works concerning with high speed NLF airfoils do not consider the robust characteristic,i.e.,the sensitivity of natural laminar flow at high Mach number (Ma) condition,which is the main obstacle for the application of transonic NLF airfoils in engineering.

In this paper,the robust design method is introduced to extent the application of NLF airfoils to a transonic,high Reynolds number flow case. A more flexible parameter method is introduced to include a more wide design space. The γ-Re_θt transition model is used to analyze the transition flow. Considering the uncertainty of the flight Ma,a robust design method is introduced to improve the aerodynamic characteristic at a large Ma range. At last,an airfoil is investigated by comparing the configuration and flow condition between the turbulent flow supercritical airfoil and the optimized airfoil. 2 Numerical simulation of boundary layer transition

The development of a design (optimization) system to potentially reduce the viscous drag depends on the reliability of transition prediction methods. Langtry and Menter^{[4, 5]} have developed a local correlation-based transition model,called the γ-Re_θt model coupled with the shear-stress transport (SST) k-ω turbulent model in 2004^[6],which satisfies most of the requirements for a transition model fully compatible with a conventional computational fluid dynamics (CFD) solver.

The γ-Re_θt transition model solves the intermittent function to control the production of turbulent kinetic energy to achieve local transition. The transport equation of the intermittent function is

where P_γ is the source term of γ,E_γ is the destruction term of γ,ρ is the density,F_l is an empirical correlation that controls the length of the transition region,Fo is the transition onset function,Ft is used to disable the destruction/relaminarization source outside a laminar boundary layer or in the viscous sublayer,c_e1 = 1,c_e2 = 50,c_a1 = 2,and c_a2 = 0.06 are model constants,S is the magnitude of the strain-rate,μt is the kinematic Eddy viscosity,and is the magnitude of the vortices.

The standard transport equation to calculate Re_θt is

where ts is the time scale,Re_θt is the momentum thickness Reynold number,P_θt is the source term of Re_θt,σ_θt = 2.0 is the model constant,and F_θt is the mixing function.

The intermittent function is the probability of turbulence,which is introduced to the SST k- w model to control the production and dissipation term of the turbulent kinetic energy equation.

The transonic flow around NLR7301^[7] airfoil is analyzed to validate the precision of the model. The test of the airfoil is accomplished in the NLR wind tunnel. The result is broadly quoted to validate the numerical method of the transition. The computation condition is that Ma = 0.3 ~ 0.85 and Re = 2.0 × 10⁶. The freestream turbulence intensity T u = 0.6%,and the viscosity ratio μ_t/μ = 2.5.

Figure 1 shows that the drag coefficient (C_d) varies with the Mach number predicted by the γ-Re_θt transition model. It shows that the data agree well with experimental data. When Ma is less than 0.7,the max numerical error is less than 9%,which satisfies engineering rules. When Ma is great than 0.7,the shock wave appears and interacts with the boundary layer,and the laminar boundary layer transition occurs due to the shock waves.

Figure 2 shows that the transition position of the upper surface of the airfoil changes with Ma. The transition position is accurately predicted,and all errors fully fill the demand of engineering.

3 Free-form deformation (FFD) method based on non-uniform rational B-spline (NURBS) curve

The FFD method^{[8, 9, 10]} is a very versatile approach proposed by Sederberg and Parry which deforms the object by moving the node of the control box surrounding the object. Mathe-matically,the deformation is a mapping function X = F(x) from R³ to R³. Sederberg and Parry^[8] used the trivariate Bernstein polynomial and control frame to construct the function F(x). The FFD method can be easily applied to any kind of complex objects. In this paper,the NURBS-based FFD method is established to achieve the airfoil deformation. The displacement ΔX of any node point X(s,t,u) in the control box is defined as follows:

where P_i,j,k and ΔP_i,j,k are the original coordinate matrix and the displacement of the control point (i,j,k),respectively. The number of the control points is (s,t,u) are the local curvilinear coordinates constructed by the control box. Bⁱ⁻¹ _l−1 (s) is the (i−1)th Bernstein polynomial of degree (l − 1) defined by Equation (6) can be rewritten in matrix form as follows:

The FFD method can be implemented by the following three steps. First,construct a control box around the object and calculate the local curvilinear coordinates (s,t,u) of all nodes of the object lying in the box and the blending function B(s,t,u). Secondly,move manipulating control points of the lattice based on the design variables and define ΔP. Finally,calculate the deformation object by Eq. (6).

Figures 3 and 4 show a simple application of the FFD method to a two-dimensional (2D) airfoil. Figure 3 shows the original airfoil and the position of each point in the control volume. Figure 4 shows the comparison between the deformation of the airfoil after moving the control point and the original airfoil.

4 Construction of aerodynamic optimization design system 4.1 Optimization algorithm

In this paper,the non-dominated sorting genetic algorithm II (NSGA-II)^[11] is used to con- duct the multi-object optimization. The NSGA-II is developed based on the original non- dominated sorting genetic algorithm (NSGA). There are three main shortages for the NSGA method^[11]: (i) computational complexity; (ii) non-elitism approach; (iii) a sharing parameter needed to be specified. Deb et al.^[11] modified the NSGA approach in order to alleviate all the above difficulties. They named the algorithm as NSGA-II. The NSGA-II is a multi-objective evolutionary algorithm developed for handling complex aerodynamic optimization problems. It behaves perfectly with two- or three-objective problems. The flow chart of the design procedure is shown in Fig. 5.

4.2 Experimental design and surrogate model

The wildly used experimental design methods are uniform design,randomized design,or- thogonal design,Latin hypercube design,etc. In this paper,the Latin hypercube method is used to select samples. The surrogate model is a key technology for engineering design,which can improve the design efficiency and explore a wider design space. The Kriging model^{[12, 13]} is very efficient and has sufficient flexibility to represent nonlinear problems. Therefore,it is introduced to replace the CFD to evaluate the characteristic of the airfoil here. 4.3 Robust design

Traditional deterministic optimization techniques tend to “over-optimize” with uncertain parameters. These techniques produce good solutions at the design point but poor characteris- tics at the off-design point. Therefore,the robust design method^{[14, 15]} is introduced to provide the solutions which are insensitive to the uncertainty of input parameters. The robust design optimizes the mean and the variance of the airfoil performance simultaneity in an uncertainty range. The Mach number submits the uniform distribution within the range from 0.72 to 0.78. The optimization object is to minimize the mean and the variance of the drag coefficient over the range of the operating condition. The optimization mathematical model can be expressed as follows:

where C_l is the lift coefficient,C_ld is the design lift coefficient,x^min _i and x^max _i are the upper boundary and the lower boundary of design variables,and μ and σ are the mathematics ex- pectation and the variance of the drag coefficient which can be obtained by the Monte Carlo simulation as follows: where Ns is the Monte Carlo simulation degree,C_d(X,Ma_k) is the drag coefficient predicted by the Kriging model,and P(Ma) is the probability density distributed with Ma. In the optimization system established in this paper,it is clear that the drag characteristic in the uncertain Ma range is controlled by μ,which can only reflect the trend of the total drag characteristic. The function of the variance σ is to change the slope of the drag-Mach curve with different Ma. 5 Robust design of NLF supercritical airfoil

The NASA SC(2)-0412 airfoil is selected as the original airfoil. Assume that the Mach number obeys normal distribution from 0.72 to 0.78,the flight turbulence intensity is T u = 0.2%,and the turbulent viscosity ratio is 10. The cruise design condition is as follows:

where Clc is the cruise lift coefficient.

In this optimization case,1 000 samples are produced by the Latin hypercube method,and the object is evaluated for each sample by the numerical simulation of boundary layer transition with the γ-Re_θt transition model. Then,the surrogate model is set up based on the Kriging model. The following aerodynamic optimization mathematical model is established based on the design requirement:

where V is the design variable space,and Ns = 3 000. The parameters of NSGA-II are that the population size is 100,the crossover probability is 0.8,the mutation rate is 0.05,and the maximum iteration number is 1 000. The non-dominated solution obtained by NSGA-II is shown in Fig. 6.

After analyzing the results of the Pareto front points,four representative airfoils are selected to analyze the characteristic of the robust natural flow supercritical airfoil,as shown in Fig. 6. Airfoil Opt1 is just a single point designed airfoil with a low mean drag coefficient and a high variance. Airfoils Opt2 and Opt3 are compromises between the mean drag coefficient and the variance. Airfoil Opt4 is a robust airfoil with a high mean drag coefficient and a low variance.

Figure 7 shows the original and optimized airfoils. It is shown that when the head radius of the optimized airfoil decreases,the maximum thickness location shifts backward,and the camber of the airfoils increases. The airfoils are forward loaded after the optimization,which agrees well with the NASA result^[2]. The leading part of the upper surface of airfoil Opt2 is much more agreeable with that of Opt1,and the lower surface of airfoil Opt2 is like Opt3.

5.2 Aerodynamic characteristic analysis

Figure 8 shows the drag divergence curves of the five airfoils. From the figure,it can be seen that the drag decreases after the optimization,the drag of airfoil Opt1 suddenly increases after the design point,while the drag of robust airfoil Opt4 changes smoothly. The drag of Opt2 is lower than the drag of Opt3 with the design Mach number,while it increases suddenly at higher Ma. The drag of airfoil Opt3 is lower than the original airfoil at all the design Mach number range,but not as robust as airfoil Opt4.

To analyze the airfoil further more,the drag is divided into two parts: pressure drag and friction drag,as shown in Fig. 9. It can be seen that the divergence of the drag is controlled by the pressure drag,which reflects the development of the shock wave on the airfoil. From Fig. 9, it can be seen that the friction drag of the airfoil decreases by considering the laminar flow. When Ma is low,the friction drag is larger than the pressure drag; when Ma goes larger,the wave drag is lager than the friction drag,especially for airfoil Opt1. Airfoil Opt3 has the lowest drag among the airfoils. Airfoil Opt4 is much like NASA0412,which is essentially a robust airfoil at the design condition.

The design principle is different from the subsonic cases,where the pressure drag is low and steady,the object is just to postpone the transition and decreases the fraction drag. However, when it comes to transonic cases,the control of the shock wave and its development are the primary task. Therefore,the pressure drag is the primary drag for higherMa. It is a compromise for a laminar flow supercritical airfoil. The friction drag will concede to improve the drag divergence characteristic. 5.2 Pressure and friction coefficient distribution

The pressure distribution of original and optimized airfoils is shown in Fig. 10(a) with

It is seen that for airfoil Opt1,a longer favorable pressure gradient region exists,and it ends with a shock wave,which means a long laminar region. For airfoil Opt2,the shock wave is weaken, and the favorable pressure gradient region is shorter. The weak shock wave benefits the drag divergence characteristic of the airfoil. For airfoil Opt3,the shock is forward and weaker than the original airfoil. Therefore,the laminar flow region is also shorten. The minimum pressure point of the lower surface shifts backward to decrease the friction drag. A pressure coefficient with a weak shock wave will shorten the laminar flow area,and the friction drag will increase. However,the development of shock is controlled. Therefore,as Ma increases,the shock will not be as strong as the single point optimized one. This is the principle of the natural laminar flow supercritical airfoil,which is much different from the traditional laminar flow airfoil.

Figure 10(b) shows the friction coefficient distribution of the four airfoils at the design cruise condition. The laminar flow region of both surfaces increases after optimization. Airfoil Opt1 gives the latest transition position among the airfoils. The laminar flow area is primary on the lowest surface of all the optimized airfoils. The laminar flow area of the upper surface of the robust airfoils is shorter than airfoil Opt1. This is a concession to the robust characteristic.

The pressure distribution of the original and optimized airfoils with Ma = 0.78 is shown in Fig. 11(a). The shock wave of airfoil Opt1 is more intense than any other airfoil,and the position shifts backward. This verifies the design principle of the laminar flow supercritical airfoil. A weak shock wave exists at the lower surface,which is due to a favorite pressure gradient to the laminar flow region. The weak shock wave can control the development of the shock and the drag divergence characteristic of the airfoil.

Figure 11(b) shows the friction coefficient distribution of the five airfoils with Ma = 0.78. The laminar flow regions of both surfaces increase after optimization. Airfoil Opt1 gives a much later transition position than any other airfoil,the shock on the upper surface is so strong that the boundary layer separates. The relaminarization after the shock wave also occurs around the airfoil. The transition point moves backward as Ma increases. Too later transitions mean a small friction drag and a strong shock. Therefore,the transition point should not move too backward,and the shock will keep weak. 5.3 Pressure gradient and transition position

To analyze the design rule of robust natural laminar flow supercritical airfoils,the maximal pressure gradient is calculated. Figure 12 shows the position and value of the maximal pressure gradient. It is clear that for a robust airfoil,the position of the maximal pressure gradient should move forward at lower Ma like Opt4,i.e.,the pressure should resume earlier,because as Ma increases,the shock will move backward and become stronger. However,this will make the flow transit to turbulence earlier. Therefore,Opt2 and Opt3 are good choices.

The transition locations of the upper and lower surfaces of the original and optimized airfoils with varying Ma are shown in Fig. 13. The laminar flow region of the upper surface becomes larger as Ma increases,while the laminar flow region of the lower surface is steady. To ensure the drag divergence characteristic,the transition point on the upper surface should not be too backward.

6 Conclusions

The robust optimization design of high-speed laminar flow airfoils is investigated. By com- paring the obtained results,the following conclusions can be summarized:

(i) Natural laminar flow can be achieved on a transonic supercritical airfoil at the chord Reynolds number of 10 million. The drag of the airfoil can be greatly decreased by considering the laminar flow technology.

(ii) Robust design should be introduced to trade-off the drag and the drag divergence char- acteristic of the airfoil. The compromise between the mean value and the variance is achieved by multi object evolution strategy.

(iii) Wind tunnel experiments should be conducted to test the results and idea of the work in the future. Higher Mach number and Reynolds number condition cases should be investigated, too. The thickness and design lift coefficient will be changed to find the design philosophy of a high speed natural flow airfoil at high Reynolds numbers.