# Aerodynamic Design Optimization of Transonic Natural-Laminar-Flow Airfoil at Low Reynolds Number

Authors: Jing Li; Cong Wang; Huiting Bian
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The position and size of laminar separation bubble on airfoil surfaces exert a profound impact on the efficiency of transonic natural-laminar-flow airfoil at low Reynolds number. Based on the particle swarm algorithm, an optimization methodology in the current work would be established with the aim of designing a high and robust performance transonic natural-laminar-flow airfoil at low Reynolds number. This methodology primarily includes two design processes: a traditional deterministic optimization at on-design point and a multi-objective of uncertainty-based optimization. First, a multigroup cooperative particle swarm optimization was used to obtain the optimal deterministic solution. The crowing distance multi-objective particle swarm optimization and the non-intrusive polynomial chaos expansion method were then adopted to determinate the Pareto-optimal front of uncertainty-based optimization. Additionally, theγ−Re¯θt transition model was employed to predict the laminar-turbulent transition. Regarding to the established optimization methodology, a propeller tip airfoil of solar energy unmanned aerial vehicle was finally designed. During optimization processes, the minimized pressure drag was particularly chosen as the optimization objective, while the friction drag increment served as a constraint condition. The deterministic results indicate that the optimized airfoil has a good ability to control the separation and reattachment positions, and the pressure drag can be greatly reduced when the laminar separation bubble is weakened. The multi-objective results show that the uncertainty-based optimized airfoil possesses a significant robust performance by considering the uncertainty of Mach number. The findings evidently demonstrate that the proposed optimization methodology and mathematical model are valuable tools to design a high-efficiency airfoil for the propeller tip.


---

## Body

## 1. Introduction

The solar powered unmanned aerial vehicle (UAV), generally flying at an altitude of 10–40 km, where it is with low atmospheric density, has attracted wide attention due to its significant practical engineering application in recent years [1, 2]. Nevertheless, note that the low atmospheric density existing in high altitude environment could result in the Reynolds number range of the propeller airfoils falling within 105–106 [3]. Subsequently, the aerodynamic performance of airfoil could be deteriorated at such low Reynolds number (LRN), and then, the propeller efficiency would be plummeted [4, 5]. As such, inspired by the importance and urgent need of low-speed flight at LRN, comprehensive studies have been dedicated to improve the aerodynamic efficiency of its propeller [6, 7]. Zhang et al. carried out the optimization design of E387 and PSU 94–097 airfoils at Re=2.0×105 and Re=4.0×105, respectively [8]. The designed airfoils showed better aerodynamic characteristics predicted through the Xfoil code than the original ones. Ma et al. established a hierarchical multi-objective optimization platform based on Xfoil code [9]. Their results indicated that three new airfoils optimized from LRN high-lift airfoils have excellent aerodynamic characteristics in both on-design and off-design states. Wang et al. accomplished the multi-objective optimization design for SD7073 airfoil and proposed a significant design concept that is to change boundary layer properties via controlling the transition positions [10]. Chen et al. discovered the low lift-to-drag ratio for the obtained tandem airfoil at Re = 220,000 that could be overcome under the lower lift coefficient, and it particularly demonstrates a manifestation at a high angle of attack [11].These above new optimized natural-laminar-flow (NLF) airfoils presented the outstanding performance in different aspects for low-speed flight at LRN, but they can scarcely be implemented in the propeller tip. The reason for this is that the velocity nearby the blade tip would be high subsonic or transonic, while the velocity of its root region is rather low [12, 13]. Up to now, there are very few transonic LRN airfoils that are available for the blade tip of solar powered UAV. Undoubtedly, it gives rise to more difficulties in designing transonic LRN airfoil compared with low-speed NLF one. Moreover, LSB will inevitably appear on one or both surfaces of the airfoil, which are mainly responsible for increasing sensitivity of aerodynamic performance as the Reynolds number decreases to and even below 500,000 [14–16]. For this, it is necessary to employ the active flow control methods to completely eliminate the laminar separation [17, 18]. Furthermore, it has a high probability of developing shock waves on the upper airfoil surface at transonic speed flight. These phenomena would subsequently produce the strong adverse pressure gradients and compressibility effect that greatly influence the size and position of LSB and enhance the design complexity [19, 20]. For another, the fluctuations of the flight conditions, leading to the instabilities of the sizes and locations of separation bubbles, also cause more difficulties [21].According to these previous studies [12–21], it is badly in need of putting forward an effective optimization methodology to design the transonic NLF airfoil at LRN towards achieving direct application in the blade to enhance the propeller’s efficiency. The particle swarm optimization (PSO) algorithm has been widely used in aerodynamic optimization because of its excellent optimal-searching ability [22–24]. Khurana et al. built an adaptive mutation-PSO (AM-PSO) method in use for redesigning airfoil at flight envelopes encompassing low-to-high Mach numbers [25]. Their results illustrated that the AM-PSO can seek out the optimal airfoil with minimal computational resources. Tao et al. proposed an improved PSO algorithm on the basis of centroidal Voronoi tessellations (CVTs) and fitness distance ratio (FDR) [26]. The optimization of single point and robust showed that the optimized airfoil and wing both have outstanding characteristics at on-design point. Masoumi et al. developed a hybrid algorithm combining the big bang-big crunch and PSO, which can be applied for reverse and direct optimization approaches [27]. And this proposed algorithm has proper accuracy and convergence speed that facilitates to reach the desired airfoil geometry. However, it should be mentioned that most of the above PSO algorithms would fall into a local optimum. The main reason is the evolution information for individual guidance only coming from the single group. Under such conditions, the evolutionary ability of algorithm could be declined and thereupon gives rise to an algorithm premature in local optimum [28–30]. Of particular interest is the high-order nonlinear feature of the aerodynamic design space that certainly requires the algorithm to have a strong global optimization capability. Moreover, massive computational consumption of aerodynamic characteristics needs the algorithm to have a faster convergence speed [31, 32].With respect to the significant role of single-objective PSO algorithms in deterministic optimization, the high robustness airfoil should be taken into consideration, in which the performances have less sensitivity in response to flight conditions’ uncertainty. Usually, the mean and variance of the characteristics were treated as the indicators for the uncertainty design. However, the mean and variance always were converted into a single-objective optimization by a linear-weight combination [33, 34]. So, the optimization results would be influenced by the treatment of weight coefficients. In this study, the mean and variance are treated as two objectives to get over the disadvantage of the weight method. In line with previous study, the multi-objective particle swarm optimization (MOPSO), employed in many multi-objective aerodynamic optimization, tried to get the Pareto-optimal front [35]. Azabi et al. proposed an interactive MOPSO with the benefits of reducing computational time and avoiding the complexity of postdata analysis [36]. Huang et al. constructed an improved MOPSO with principal component analysis methodology, enhancing the convergence of Pareto front [37]. The level of non-dominanted Pareto front was effectively improved, and the multi-objective requirements of a wide-body aircraft were satisfied. However, the number of the Pareto-optimal front will grow dramatically for multi-objective aerodynamic design. The reasonable dominant ranking and selecting for Pareto-optimal sets are essential factors in enhancing the capability of the MOPSO [38].In present work, an optimization methodology based on PSO algorithm was proposed. It aims to design a high and robust performance transonic NLF airfoil at LRN. This article is organized as follows: in Section2.1, the PSO methodology is built up, and the algorithm performances are validated by the aerodynamic design of a supercritical airfoil; the numerical simulation and airfoil parameterization methods are given in Sections 2.2 and 2.3, respectively; in Section 3, deterministic optimization and multi-objective robust design are conducted to get an optimal airfoil to meet the design requirements; and in Section 4, several important design conclusions are summarized by the comparison of original and newly optimized airfoil.


## 2. Methodology

### 2.1. Optimization Methodology Based on the PSO Algorithm

#### 2.1.1. Multigroups Cooperative Particle Swarm Optimization Algorithm (MCPSO)

To avoid algorithm premature, a MCPSO algorithm based on the biological systems symbiosis is put forward in this section. The populations are divided into several subgroups. Every subgroup obtains the evolution information from the group itself and the other symbiotic subgroups. The symbiotic relationships among subgroups are described by a master-slave structure. The slave subgroups and the master one separately undertake the global and local searches.The basic principles of this algorithm are conducted as follows: the populations are divided into four subgroups, marked as the 1st, 2nd, 3rd, and 4th subgroups. The 1st subgroup is chosen as the master subgroup, which is mainly in charge of local search, and the others are specified as the slave subgroups for the global search. The 1st subgroup is to adopt the standard PSO, while the others are to adopt the improved PSO algorithms accomplished in our prior studies [39, 40]. The speed and position of the 1st subgroup is updated according to the following equation:(1)v⟶it+1=ωv⟶it+c1r1p⟶i1t−x⟶it+c2r2p⟶g1t−x⟶it+fc3r3p⟶g2,⋯t−x⟶it,x⟶it+1=x⟶it+v⟶it+1,where v⟶i and x⟶i are the speed and position of the particle, respectively; p⟶i1 and p⟶g1 are the local and global best particles of the 1st subgroup, respectively; p⟶g2,⋯ is the global best particle of the other subgroups; c1 and c2 are the learning factors; and f∈0,1 is the migration factor.The ability to search the optimum is verified by designing the supercritical airfoil NASA-SC (2) 0712. The on-design point isMa∞=0.76,CL=0.70,andRe=3.0×106, and the objective is minimizing the total drag at this cruise state. The following optimized mathematical model is built up for this design problem:(2)objective:minCD,s.t.:CL=0.7,thicknessmax≥0.12,where thicknessmax is the maximum thickness of the airfoil.In this testing aerodynamic design, the characteristics are obtained by solving the Reynolds-averaged Navier–Stokes (RANS) equations [41]. The airfoil parameterization method adopted in this design will be displayed in Section 2.3. There are five design variables for the upper or lower surfaces, ten design variables in all for each airfoil. The cell number of the grids is 29,000.The standard PSO is marked as GPSO. The population size of the two algorithms (GPSO and MCPSO) counts 60 particles, and the convergence condition is the maximum iteration which is set as 50. The damping strategy is utilized for the boundary treatment method [42].Aerodynamic characteristics for the original and optimized airfoils at on-design point are shown in Table1. The drag coefficient of MCPSO is smallest and decreases 32.8 counts by comparing with the original one. The lift-to-drag ratio increases by 24.5%.Table 1
Aerodynamic characteristics for the original and optimized airfoils at on-design point.
CharacteristicsNASA SC (2)-0712GPSOMCPSODrag coefficient0.0166780.0135050.013394Lift-to-drag ratio41.9751.8352.26Maximum thickness12%12%12%The optimization results of GPSO and MCPSO are shown in Figure1. The convergent curves indicate that MCPSO has less probability of falling into the local optimum than GPSO. One thing should be mentioned that the fitness of MCPSO could attain 0.13505 in the 25th iteration, while such value only could be achieved by GPSO in its 40th iteration. Furthermore, the GPSO falls into the local optimum at the iteration end, compared with the better fitness (0.013394) of the MCPSO. Seen from the surface pressure distributions presented in Figure 1(b), the strong shocks disappear on the upper surface of both two optimized airfoils. Notice that the suck peak of optimized surface pressure distribution is higher than that of the original one. The results indicate that the convergence speed of MCPSO is faster than that of GPSO, and the better optimal solution could be found in the design space quickly and accurately.Figure 1
Optimization results: (a) convergent curves; (b) comparisons of pressure distributions and airfoils.
(a)(b)


#### 2.1.2. MOPSO Algorithm Based on Crowding-Distance Sorting

The MOPSO proposed in the current study is based on the crowding-distance sorting, and the update and control of Pareto sets are as follows:(1)
Assume that the maximum size of the Pareto sets named asAt is M, and the size of At in the t iteration is m1m1≤M. There will be m2 new Pareto solutions in the following iteration, which puts together with At into another temporary Pareto sets named as At′.(2)
If the two solutions have the same objectives, randomly eliminate one of the duplicate solutions in the temporary Pareto setsAt′.(3)
The size ofAt′ will be m3m3≤m1+m2, when deleting the dominated solutions in the temporary Pareto sets At′.(4)
Sort the populations ofAt′ in the ascending order according to the crowing distance, marked as At″. The crowing-distance computation is as follows:(3)idistance=∑m=1Nobj∑i=2Nset−1fi+1m−fi−1mfmaxm−fminm,whereidistance refers to the crowing distance of the ith population; fi+1m refers to the mth objective function value of the i + 1th individual; fmaxm and fminm are the maximum and minimum values of the mth objective functions, respectively; and Nobj and Nset are the numbers of objectives and set At′. Herein, the crowing-distance computation requires sorting the population corresponding to each objective function value in ascending order.(5)
The sorting setsAt″ will be chosen as the new At′ in next iteration if m3≤M. Otherwise, the exact top M population members from the At″ will be chosen as the new At′.In general, the optimization results of MOPSO can be verified by generational distance (GD) and spacing metric (SP). The GD can represent the distance between the non-dominated set and true Pareto-optimal front, of which the formula is as follows:(4)GD=∑i=1ndi2n,where n is the size of nondomination set and di is the smallest distance between the ith individual and true Pareto-optimal front. Notice that the smaller the value of GD is, the better the convergence is.The SP is used to calculate the distribution uniformity, and its formula is as follows:(5)SP=1n−1∑i=1n∑i=1n−1din−1−di2.A smaller SP denotes a better uniformity of the sets. Thus, three typical multi-objective functions labeled by ZDT1, ZDT2, and ZDT3 are to verify the performance:ZDT1:(6)minf1x=x1minf2x=gx1−x1gxgx=1+9∑i=2nxin−1n=30,x∈0,1.ZDT2:(7)minf1x=x1minf2x=gx1−x1gx2gx=1+9∑i=2nxin−1n=30,x∈0,1.ZDT3:(8)minf1x=x1minf2x=gx1−x1gx−x1sin10πx1gxgx=1+9∑i=2nxin−1n=30,x∈0,1.Figure2 shows the Pareto-optimal fronts of these three functions obtained from MOPSO. The newly designed fronts are extremely consistent with the real ones. Moreover, as shown in Table 2, the GDs of MOPSO are better than those of NSGA-II [43] for all three functions. As for SPs, it is hard to distinguish between the two algorithms to find which one performs better.Figure 2
The comparison of the real and Pareto-optimal fronts obtained from MOPSO: (a) ZDT1; (b) ZDT2; (c) ZDT3.
(a)(b)(c)Table 2
The GD and SP obtained from MOPSO and NAGA-II in terms of three typical functions.
Test functionsGDSPNSGA-IIMOPSONSGA-IIMOPSOZDT18.146e − 23.396e − 49.539e − 27.539e − 2ZDT22.334e − 22.127e − 36.722e − 21.269e − 1ZDT31.037e − 15.957e − 34.613e − 15.613e − 2


#### 2.1.3. Uncertainty Quantification by Polynomial Chaos Expansion (PCE)

The main difficulties for uncertainty quantification are the excessive computational cost and unsatisfactory accuracies of the mean and variance. For instance, the widely used Monte Carlo simulation (MCS) [44] method basically requires at least 104 samples to get sufficiently accurate estimate of the mean and variance, which certainly causes a much high cost of computation. In the aspect of efficient and accurate uncertainty quantification, the PCE is a promising method for accurate estimation of mean and variance of the objectives [45, 46]. A non-instrusive PCE method is adopted in this part, which is a question-answer system and unnecessary to modify the original governing equations [47].This formula assumes that the random variableθ obeys a standard normal distribution, and thus, p polynomial order of a general random process Xξθ can be expressed as follows:(9)Xξθ=a0H0+∑i1=1∞ai1H1ξi1θ+∑i1=1∞∑i2=1i1ai1,i2H2ξi1θ,ξi2θ+⋯=∑j=1NP+1bjψjξ,where Hn is the Hermite polynomial, ξi1,…,ξin are the vectors consisting of n Gaussian variables, NP+1=p+n!/p!n!, and n is the number of random variables.Then, coefficientsbj can be solved by the following formula:(10)bj=Xξ,ψjξψjξ,ψjξ=∫Xξψjξωξdξψjξ,ψjξ=∑i1=1k⋯∑iN=1kωi1⋯ωiNXxi1,…,xiNψjxi1,…,xiNψjξ,ψjξ,where the number of collocation points is indicated by M=kn and ωi and xi are the Gaussian–Hermite quadrature coefficient and point, respectively.According to the orthogonality of Gaussian–Hermite polynomials, the mean and variance ofXξθ are presented as follows:(11)μXξθ=b1,σ2Xξθ=∑j=2NP+1bj2ψj2.The accuracy and efficiency for PCE and MCS approaches are compared through the following function:(12)fX=20e−x+2.72/22π+e−x−2.72/2222π2,where X is assumed to be normally distributed with μ=2.0 and σ=0.8. The random variables are assumed as ξi∈N0,12, and xi is expressed as xi=2.0+0.8ξi.The real mean and variance of this function would be calculated by the following equations:(13)μY=EfX=3.520985,σY=EfX−μY=0.507313.Figure3 plots the comparisons of estimated mean (μY) and variance (σY) between PCE and MCS methods. The sample numbers for MCS fetch two values, i.e., 100 and 500, that are utilized to appraise its estimation accuracy of the stochastic process. The results demonstrate that the mean and variance will become convergent with the increasing order P. And the estimation accuracy of mean is slightly higher than that of variance, particularly for the P value within 2–4. Additionally, the estimated μY and σY by the PCE method both fall within 0.5 percent error bands, when the polynomial order P is up to 4. Meanwhile, the estimated accuracy of μY agrees well with that of MCS using 500 samples. But the estimated accuracy of σY becomes closer to that of MCS using 500 samples with the polynomial order P increasing from 4. Thereupon, our results demonstrate that the PCE method is reasonable for the uncertain quantification with respect to its reliable accuracy and efficiency. And the same conclusion was also derived by the former study that verified the estimated results of aerodynamic characteristics by PCE [48]. Figure 3(c) shows the time cost of PCE and MCS approaches varying with the error of variances. When the error of variance is at the same level, the time cost of the MCS approach is much more than that of the PCE approach. Therefore, the PCE method is chosen to estimate the mean and variance of the objectives in present work for its efficiency.Figure 3
Comparison of estimated mean and variance between PCE and MCS methods: (a) meanμY; (b) variance σY; (c) time cost.
(a)(b)(c)


### 2.2. Numerical Simulation Method for Aerodynamic Characteristics

Theγ−Re¯θt transition model [49, 50] is combined with the local variable of SST k−ω turbulence model to predict transition. The proposed model incorporated intermittent function with transition experience relationship. So, it can control the distribution of intermittent functions in boundary layers through the relevance functions of transition momentum thickness Reynolds number and afterward dominates turbulence through the intermittent functions. Moreover, since the momentum thickness Reynolds number is related to the strain rate, the calculation of momentum thickness can be prevented, thereby enhancing the prediction accuracy of separation transition [50].The transition model validation is performed on NLR7301 airfoil applying the same boundary conditions adopted in the experiments. The simulation condition isMa=0.747andRe=2×106. The airfoil was tested at the NLR pilot closed-circuit wind tunnel [51]. There are 760×179 structured grid for the transition model solver, as shown in Figure 4. As displayed in Figure 5, the transition position (TP) and surface pressure distribution obtained from the transition model and experiment data illustrate good consistency.Figure 4
The computation structured grid: (a) grid around the airfoil; (b) leading edge grid.
(a)(b)Figure 5
The aerodynamic characteristics of the NLR 7301 obtained from CFD and experiment: (a) AoA effect on TP of upper surface; (b) AoA effect on TP of lower surface; (c) surface pressure distribution when AoA is equal to −2; (d) surface pressure distribution when AoA is equal to 2.
(a)(b)(c)(d)Also, the validation of capturing LSB is accomplished on E387 airfoil at low speed andRe=2×105. The airfoil performances predicted by transition model are compared to University of Illinois Urbana-Campaign (UIUC) wind-tunnel data [52] in Figure 6. The results displayed in Figure 6(a) demonstrate that the lift-drag curve is in good agreement with the experimental one. As presented in Figure 6(b), both the location and size of LSB obtained from CFD are coincident with the experiment data.Figure 6
The aerodynamic characteristics of the E387 airfoil obtained from CFD and experiment: (a) lift-drag curves; (b) locations of the upper surface flow features.
(a)(b)Consequently, these phenomena make us to conclude that the numerical simulation technology of transition boundary layer is accurate and reliable for the aerodynamic design optimization of transonic NLF airfoil at LRN.


### 2.3. Airfoil Parameterization Method

The disturbing basis function, based on the class/shape function transformation (CST) [53], is introduced in airfoil parameterization. The parameterization method constructed here can not only minimize the errors caused by fitting the initial parameters but also preserve the advantages of the original CST parameterization method.The original CST method is shown in the following equation:(14)ζψ=ψ0.51−ψ∑i=0nAiKiψi1−ψn−i+ψζT,where ψ and ζψ are the x-coordinates and y-coordinates of the airfoil, respectively; ζT is the y-coordinate of the trailing edge point; Ai is the design variables; and Ki=Cni=n!/i!n−i!.The disturbanceΔAi is added to the initial design variable Ai, and then, equation (14) can be written as(15)ζψ′=ψ0.51−ψ∑i=0nAi+ΔAiKiψi1−ψn−i+ψζT=ζψ+Δζψ,where ζψ′ and ζψ are the separable y-coordinates of new and original airfoils and Δζψ are the disturbances of y-coordinates, which can be denoted as(16)Δζψ=ψ0.51−ψ∑i=0nΔAiKiψi1−ψn−i.According to the aforementioned equations, the newy-coordinates can be treated as two parts, i.e., the original y-coordinates (ζψ) and disturbances (Δζψ). Thus, equation (16) can be denoted as the disturbing basis function based on CST.


## 2.1. Optimization Methodology Based on the PSO Algorithm

### 2.1.1. Multigroups Cooperative Particle Swarm Optimization Algorithm (MCPSO)

To avoid algorithm premature, a MCPSO algorithm based on the biological systems symbiosis is put forward in this section. The populations are divided into several subgroups. Every subgroup obtains the evolution information from the group itself and the other symbiotic subgroups. The symbiotic relationships among subgroups are described by a master-slave structure. The slave subgroups and the master one separately undertake the global and local searches.The basic principles of this algorithm are conducted as follows: the populations are divided into four subgroups, marked as the 1st, 2nd, 3rd, and 4th subgroups. The 1st subgroup is chosen as the master subgroup, which is mainly in charge of local search, and the others are specified as the slave subgroups for the global search. The 1st subgroup is to adopt the standard PSO, while the others are to adopt the improved PSO algorithms accomplished in our prior studies [39, 40]. The speed and position of the 1st subgroup is updated according to the following equation:(1)v⟶it+1=ωv⟶it+c1r1p⟶i1t−x⟶it+c2r2p⟶g1t−x⟶it+fc3r3p⟶g2,⋯t−x⟶it,x⟶it+1=x⟶it+v⟶it+1,where v⟶i and x⟶i are the speed and position of the particle, respectively; p⟶i1 and p⟶g1 are the local and global best particles of the 1st subgroup, respectively; p⟶g2,⋯ is the global best particle of the other subgroups; c1 and c2 are the learning factors; and f∈0,1 is the migration factor.The ability to search the optimum is verified by designing the supercritical airfoil NASA-SC (2) 0712. The on-design point isMa∞=0.76,CL=0.70,andRe=3.0×106, and the objective is minimizing the total drag at this cruise state. The following optimized mathematical model is built up for this design problem:(2)objective:minCD,s.t.:CL=0.7,thicknessmax≥0.12,where thicknessmax is the maximum thickness of the airfoil.In this testing aerodynamic design, the characteristics are obtained by solving the Reynolds-averaged Navier–Stokes (RANS) equations [41]. The airfoil parameterization method adopted in this design will be displayed in Section 2.3. There are five design variables for the upper or lower surfaces, ten design variables in all for each airfoil. The cell number of the grids is 29,000.The standard PSO is marked as GPSO. The population size of the two algorithms (GPSO and MCPSO) counts 60 particles, and the convergence condition is the maximum iteration which is set as 50. The damping strategy is utilized for the boundary treatment method [42].Aerodynamic characteristics for the original and optimized airfoils at on-design point are shown in Table1. The drag coefficient of MCPSO is smallest and decreases 32.8 counts by comparing with the original one. The lift-to-drag ratio increases by 24.5%.Table 1
Aerodynamic characteristics for the original and optimized airfoils at on-design point.
CharacteristicsNASA SC (2)-0712GPSOMCPSODrag coefficient0.0166780.0135050.013394Lift-to-drag ratio41.9751.8352.26Maximum thickness12%12%12%The optimization results of GPSO and MCPSO are shown in Figure1. The convergent curves indicate that MCPSO has less probability of falling into the local optimum than GPSO. One thing should be mentioned that the fitness of MCPSO could attain 0.13505 in the 25th iteration, while such value only could be achieved by GPSO in its 40th iteration. Furthermore, the GPSO falls into the local optimum at the iteration end, compared with the better fitness (0.013394) of the MCPSO. Seen from the surface pressure distributions presented in Figure 1(b), the strong shocks disappear on the upper surface of both two optimized airfoils. Notice that the suck peak of optimized surface pressure distribution is higher than that of the original one. The results indicate that the convergence speed of MCPSO is faster than that of GPSO, and the better optimal solution could be found in the design space quickly and accurately.Figure 1
Optimization results: (a) convergent curves; (b) comparisons of pressure distributions and airfoils.
(a)(b)


### 2.1.2. MOPSO Algorithm Based on Crowding-Distance Sorting

The MOPSO proposed in the current study is based on the crowding-distance sorting, and the update and control of Pareto sets are as follows:(1)
Assume that the maximum size of the Pareto sets named asAt is M, and the size of At in the t iteration is m1m1≤M. There will be m2 new Pareto solutions in the following iteration, which puts together with At into another temporary Pareto sets named as At′.(2)
If the two solutions have the same objectives, randomly eliminate one of the duplicate solutions in the temporary Pareto setsAt′.(3)
The size ofAt′ will be m3m3≤m1+m2, when deleting the dominated solutions in the temporary Pareto sets At′.(4)
Sort the populations ofAt′ in the ascending order according to the crowing distance, marked as At″. The crowing-distance computation is as follows:(3)idistance=∑m=1Nobj∑i=2Nset−1fi+1m−fi−1mfmaxm−fminm,whereidistance refers to the crowing distance of the ith population; fi+1m refers to the mth objective function value of the i + 1th individual; fmaxm and fminm are the maximum and minimum values of the mth objective functions, respectively; and Nobj and Nset are the numbers of objectives and set At′. Herein, the crowing-distance computation requires sorting the population corresponding to each objective function value in ascending order.(5)
The sorting setsAt″ will be chosen as the new At′ in next iteration if m3≤M. Otherwise, the exact top M population members from the At″ will be chosen as the new At′.In general, the optimization results of MOPSO can be verified by generational distance (GD) and spacing metric (SP). The GD can represent the distance between the non-dominated set and true Pareto-optimal front, of which the formula is as follows:(4)GD=∑i=1ndi2n,where n is the size of nondomination set and di is the smallest distance between the ith individual and true Pareto-optimal front. Notice that the smaller the value of GD is, the better the convergence is.The SP is used to calculate the distribution uniformity, and its formula is as follows:(5)SP=1n−1∑i=1n∑i=1n−1din−1−di2.A smaller SP denotes a better uniformity of the sets. Thus, three typical multi-objective functions labeled by ZDT1, ZDT2, and ZDT3 are to verify the performance:ZDT1:(6)minf1x=x1minf2x=gx1−x1gxgx=1+9∑i=2nxin−1n=30,x∈0,1.ZDT2:(7)minf1x=x1minf2x=gx1−x1gx2gx=1+9∑i=2nxin−1n=30,x∈0,1.ZDT3:(8)minf1x=x1minf2x=gx1−x1gx−x1sin10πx1gxgx=1+9∑i=2nxin−1n=30,x∈0,1.Figure2 shows the Pareto-optimal fronts of these three functions obtained from MOPSO. The newly designed fronts are extremely consistent with the real ones. Moreover, as shown in Table 2, the GDs of MOPSO are better than those of NSGA-II [43] for all three functions. As for SPs, it is hard to distinguish between the two algorithms to find which one performs better.Figure 2
The comparison of the real and Pareto-optimal fronts obtained from MOPSO: (a) ZDT1; (b) ZDT2; (c) ZDT3.
(a)(b)(c)Table 2
The GD and SP obtained from MOPSO and NAGA-II in terms of three typical functions.
Test functionsGDSPNSGA-IIMOPSONSGA-IIMOPSOZDT18.146e − 23.396e − 49.539e − 27.539e − 2ZDT22.334e − 22.127e − 36.722e − 21.269e − 1ZDT31.037e − 15.957e − 34.613e − 15.613e − 2


### 2.1.3. Uncertainty Quantification by Polynomial Chaos Expansion (PCE)

The main difficulties for uncertainty quantification are the excessive computational cost and unsatisfactory accuracies of the mean and variance. For instance, the widely used Monte Carlo simulation (MCS) [44] method basically requires at least 104 samples to get sufficiently accurate estimate of the mean and variance, which certainly causes a much high cost of computation. In the aspect of efficient and accurate uncertainty quantification, the PCE is a promising method for accurate estimation of mean and variance of the objectives [45, 46]. A non-instrusive PCE method is adopted in this part, which is a question-answer system and unnecessary to modify the original governing equations [47].This formula assumes that the random variableθ obeys a standard normal distribution, and thus, p polynomial order of a general random process Xξθ can be expressed as follows:(9)Xξθ=a0H0+∑i1=1∞ai1H1ξi1θ+∑i1=1∞∑i2=1i1ai1,i2H2ξi1θ,ξi2θ+⋯=∑j=1NP+1bjψjξ,where Hn is the Hermite polynomial, ξi1,…,ξin are the vectors consisting of n Gaussian variables, NP+1=p+n!/p!n!, and n is the number of random variables.Then, coefficientsbj can be solved by the following formula:(10)bj=Xξ,ψjξψjξ,ψjξ=∫Xξψjξωξdξψjξ,ψjξ=∑i1=1k⋯∑iN=1kωi1⋯ωiNXxi1,…,xiNψjxi1,…,xiNψjξ,ψjξ,where the number of collocation points is indicated by M=kn and ωi and xi are the Gaussian–Hermite quadrature coefficient and point, respectively.According to the orthogonality of Gaussian–Hermite polynomials, the mean and variance ofXξθ are presented as follows:(11)μXξθ=b1,σ2Xξθ=∑j=2NP+1bj2ψj2.The accuracy and efficiency for PCE and MCS approaches are compared through the following function:(12)fX=20e−x+2.72/22π+e−x−2.72/2222π2,where X is assumed to be normally distributed with μ=2.0 and σ=0.8. The random variables are assumed as ξi∈N0,12, and xi is expressed as xi=2.0+0.8ξi.The real mean and variance of this function would be calculated by the following equations:(13)μY=EfX=3.520985,σY=EfX−μY=0.507313.Figure3 plots the comparisons of estimated mean (μY) and variance (σY) between PCE and MCS methods. The sample numbers for MCS fetch two values, i.e., 100 and 500, that are utilized to appraise its estimation accuracy of the stochastic process. The results demonstrate that the mean and variance will become convergent with the increasing order P. And the estimation accuracy of mean is slightly higher than that of variance, particularly for the P value within 2–4. Additionally, the estimated μY and σY by the PCE method both fall within 0.5 percent error bands, when the polynomial order P is up to 4. Meanwhile, the estimated accuracy of μY agrees well with that of MCS using 500 samples. But the estimated accuracy of σY becomes closer to that of MCS using 500 samples with the polynomial order P increasing from 4. Thereupon, our results demonstrate that the PCE method is reasonable for the uncertain quantification with respect to its reliable accuracy and efficiency. And the same conclusion was also derived by the former study that verified the estimated results of aerodynamic characteristics by PCE [48]. Figure 3(c) shows the time cost of PCE and MCS approaches varying with the error of variances. When the error of variance is at the same level, the time cost of the MCS approach is much more than that of the PCE approach. Therefore, the PCE method is chosen to estimate the mean and variance of the objectives in present work for its efficiency.Figure 3
Comparison of estimated mean and variance between PCE and MCS methods: (a) meanμY; (b) variance σY; (c) time cost.
(a)(b)(c)


## 2.1.1. Multigroups Cooperative Particle Swarm Optimization Algorithm (MCPSO)

To avoid algorithm premature, a MCPSO algorithm based on the biological systems symbiosis is put forward in this section. The populations are divided into several subgroups. Every subgroup obtains the evolution information from the group itself and the other symbiotic subgroups. The symbiotic relationships among subgroups are described by a master-slave structure. The slave subgroups and the master one separately undertake the global and local searches.The basic principles of this algorithm are conducted as follows: the populations are divided into four subgroups, marked as the 1st, 2nd, 3rd, and 4th subgroups. The 1st subgroup is chosen as the master subgroup, which is mainly in charge of local search, and the others are specified as the slave subgroups for the global search. The 1st subgroup is to adopt the standard PSO, while the others are to adopt the improved PSO algorithms accomplished in our prior studies [39, 40]. The speed and position of the 1st subgroup is updated according to the following equation:(1)v⟶it+1=ωv⟶it+c1r1p⟶i1t−x⟶it+c2r2p⟶g1t−x⟶it+fc3r3p⟶g2,⋯t−x⟶it,x⟶it+1=x⟶it+v⟶it+1,where v⟶i and x⟶i are the speed and position of the particle, respectively; p⟶i1 and p⟶g1 are the local and global best particles of the 1st subgroup, respectively; p⟶g2,⋯ is the global best particle of the other subgroups; c1 and c2 are the learning factors; and f∈0,1 is the migration factor.The ability to search the optimum is verified by designing the supercritical airfoil NASA-SC (2) 0712. The on-design point isMa∞=0.76,CL=0.70,andRe=3.0×106, and the objective is minimizing the total drag at this cruise state. The following optimized mathematical model is built up for this design problem:(2)objective:minCD,s.t.:CL=0.7,thicknessmax≥0.12,where thicknessmax is the maximum thickness of the airfoil.In this testing aerodynamic design, the characteristics are obtained by solving the Reynolds-averaged Navier–Stokes (RANS) equations [41]. The airfoil parameterization method adopted in this design will be displayed in Section 2.3. There are five design variables for the upper or lower surfaces, ten design variables in all for each airfoil. The cell number of the grids is 29,000.The standard PSO is marked as GPSO. The population size of the two algorithms (GPSO and MCPSO) counts 60 particles, and the convergence condition is the maximum iteration which is set as 50. The damping strategy is utilized for the boundary treatment method [42].Aerodynamic characteristics for the original and optimized airfoils at on-design point are shown in Table1. The drag coefficient of MCPSO is smallest and decreases 32.8 counts by comparing with the original one. The lift-to-drag ratio increases by 24.5%.Table 1
Aerodynamic characteristics for the original and optimized airfoils at on-design point.
CharacteristicsNASA SC (2)-0712GPSOMCPSODrag coefficient0.0166780.0135050.013394Lift-to-drag ratio41.9751.8352.26Maximum thickness12%12%12%The optimization results of GPSO and MCPSO are shown in Figure1. The convergent curves indicate that MCPSO has less probability of falling into the local optimum than GPSO. One thing should be mentioned that the fitness of MCPSO could attain 0.13505 in the 25th iteration, while such value only could be achieved by GPSO in its 40th iteration. Furthermore, the GPSO falls into the local optimum at the iteration end, compared with the better fitness (0.013394) of the MCPSO. Seen from the surface pressure distributions presented in Figure 1(b), the strong shocks disappear on the upper surface of both two optimized airfoils. Notice that the suck peak of optimized surface pressure distribution is higher than that of the original one. The results indicate that the convergence speed of MCPSO is faster than that of GPSO, and the better optimal solution could be found in the design space quickly and accurately.Figure 1
Optimization results: (a) convergent curves; (b) comparisons of pressure distributions and airfoils.
(a)(b)


## 2.1.2. MOPSO Algorithm Based on Crowding-Distance Sorting

The MOPSO proposed in the current study is based on the crowding-distance sorting, and the update and control of Pareto sets are as follows:(1)
Assume that the maximum size of the Pareto sets named asAt is M, and the size of At in the t iteration is m1m1≤M. There will be m2 new Pareto solutions in the following iteration, which puts together with At into another temporary Pareto sets named as At′.(2)
If the two solutions have the same objectives, randomly eliminate one of the duplicate solutions in the temporary Pareto setsAt′.(3)
The size ofAt′ will be m3m3≤m1+m2, when deleting the dominated solutions in the temporary Pareto sets At′.(4)
Sort the populations ofAt′ in the ascending order according to the crowing distance, marked as At″. The crowing-distance computation is as follows:(3)idistance=∑m=1Nobj∑i=2Nset−1fi+1m−fi−1mfmaxm−fminm,whereidistance refers to the crowing distance of the ith population; fi+1m refers to the mth objective function value of the i + 1th individual; fmaxm and fminm are the maximum and minimum values of the mth objective functions, respectively; and Nobj and Nset are the numbers of objectives and set At′. Herein, the crowing-distance computation requires sorting the population corresponding to each objective function value in ascending order.(5)
The sorting setsAt″ will be chosen as the new At′ in next iteration if m3≤M. Otherwise, the exact top M population members from the At″ will be chosen as the new At′.In general, the optimization results of MOPSO can be verified by generational distance (GD) and spacing metric (SP). The GD can represent the distance between the non-dominated set and true Pareto-optimal front, of which the formula is as follows:(4)GD=∑i=1ndi2n,where n is the size of nondomination set and di is the smallest distance between the ith individual and true Pareto-optimal front. Notice that the smaller the value of GD is, the better the convergence is.The SP is used to calculate the distribution uniformity, and its formula is as follows:(5)SP=1n−1∑i=1n∑i=1n−1din−1−di2.A smaller SP denotes a better uniformity of the sets. Thus, three typical multi-objective functions labeled by ZDT1, ZDT2, and ZDT3 are to verify the performance:ZDT1:(6)minf1x=x1minf2x=gx1−x1gxgx=1+9∑i=2nxin−1n=30,x∈0,1.ZDT2:(7)minf1x=x1minf2x=gx1−x1gx2gx=1+9∑i=2nxin−1n=30,x∈0,1.ZDT3:(8)minf1x=x1minf2x=gx1−x1gx−x1sin10πx1gxgx=1+9∑i=2nxin−1n=30,x∈0,1.Figure2 shows the Pareto-optimal fronts of these three functions obtained from MOPSO. The newly designed fronts are extremely consistent with the real ones. Moreover, as shown in Table 2, the GDs of MOPSO are better than those of NSGA-II [43] for all three functions. As for SPs, it is hard to distinguish between the two algorithms to find which one performs better.Figure 2
The comparison of the real and Pareto-optimal fronts obtained from MOPSO: (a) ZDT1; (b) ZDT2; (c) ZDT3.
(a)(b)(c)Table 2
The GD and SP obtained from MOPSO and NAGA-II in terms of three typical functions.
Test functionsGDSPNSGA-IIMOPSONSGA-IIMOPSOZDT18.146e − 23.396e − 49.539e − 27.539e − 2ZDT22.334e − 22.127e − 36.722e − 21.269e − 1ZDT31.037e − 15.957e − 34.613e − 15.613e − 2


## 2.1.3. Uncertainty Quantification by Polynomial Chaos Expansion (PCE)

The main difficulties for uncertainty quantification are the excessive computational cost and unsatisfactory accuracies of the mean and variance. For instance, the widely used Monte Carlo simulation (MCS) [44] method basically requires at least 104 samples to get sufficiently accurate estimate of the mean and variance, which certainly causes a much high cost of computation. In the aspect of efficient and accurate uncertainty quantification, the PCE is a promising method for accurate estimation of mean and variance of the objectives [45, 46]. A non-instrusive PCE method is adopted in this part, which is a question-answer system and unnecessary to modify the original governing equations [47].This formula assumes that the random variableθ obeys a standard normal distribution, and thus, p polynomial order of a general random process Xξθ can be expressed as follows:(9)Xξθ=a0H0+∑i1=1∞ai1H1ξi1θ+∑i1=1∞∑i2=1i1ai1,i2H2ξi1θ,ξi2θ+⋯=∑j=1NP+1bjψjξ,where Hn is the Hermite polynomial, ξi1,…,ξin are the vectors consisting of n Gaussian variables, NP+1=p+n!/p!n!, and n is the number of random variables.Then, coefficientsbj can be solved by the following formula:(10)bj=Xξ,ψjξψjξ,ψjξ=∫Xξψjξωξdξψjξ,ψjξ=∑i1=1k⋯∑iN=1kωi1⋯ωiNXxi1,…,xiNψjxi1,…,xiNψjξ,ψjξ,where the number of collocation points is indicated by M=kn and ωi and xi are the Gaussian–Hermite quadrature coefficient and point, respectively.According to the orthogonality of Gaussian–Hermite polynomials, the mean and variance ofXξθ are presented as follows:(11)μXξθ=b1,σ2Xξθ=∑j=2NP+1bj2ψj2.The accuracy and efficiency for PCE and MCS approaches are compared through the following function:(12)fX=20e−x+2.72/22π+e−x−2.72/2222π2,where X is assumed to be normally distributed with μ=2.0 and σ=0.8. The random variables are assumed as ξi∈N0,12, and xi is expressed as xi=2.0+0.8ξi.The real mean and variance of this function would be calculated by the following equations:(13)μY=EfX=3.520985,σY=EfX−μY=0.507313.Figure3 plots the comparisons of estimated mean (μY) and variance (σY) between PCE and MCS methods. The sample numbers for MCS fetch two values, i.e., 100 and 500, that are utilized to appraise its estimation accuracy of the stochastic process. The results demonstrate that the mean and variance will become convergent with the increasing order P. And the estimation accuracy of mean is slightly higher than that of variance, particularly for the P value within 2–4. Additionally, the estimated μY and σY by the PCE method both fall within 0.5 percent error bands, when the polynomial order P is up to 4. Meanwhile, the estimated accuracy of μY agrees well with that of MCS using 500 samples. But the estimated accuracy of σY becomes closer to that of MCS using 500 samples with the polynomial order P increasing from 4. Thereupon, our results demonstrate that the PCE method is reasonable for the uncertain quantification with respect to its reliable accuracy and efficiency. And the same conclusion was also derived by the former study that verified the estimated results of aerodynamic characteristics by PCE [48]. Figure 3(c) shows the time cost of PCE and MCS approaches varying with the error of variances. When the error of variance is at the same level, the time cost of the MCS approach is much more than that of the PCE approach. Therefore, the PCE method is chosen to estimate the mean and variance of the objectives in present work for its efficiency.Figure 3
Comparison of estimated mean and variance between PCE and MCS methods: (a) meanμY; (b) variance σY; (c) time cost.
(a)(b)(c)


## 2.2. Numerical Simulation Method for Aerodynamic Characteristics

Theγ−Re¯θt transition model [49, 50] is combined with the local variable of SST k−ω turbulence model to predict transition. The proposed model incorporated intermittent function with transition experience relationship. So, it can control the distribution of intermittent functions in boundary layers through the relevance functions of transition momentum thickness Reynolds number and afterward dominates turbulence through the intermittent functions. Moreover, since the momentum thickness Reynolds number is related to the strain rate, the calculation of momentum thickness can be prevented, thereby enhancing the prediction accuracy of separation transition [50].The transition model validation is performed on NLR7301 airfoil applying the same boundary conditions adopted in the experiments. The simulation condition isMa=0.747andRe=2×106. The airfoil was tested at the NLR pilot closed-circuit wind tunnel [51]. There are 760×179 structured grid for the transition model solver, as shown in Figure 4. As displayed in Figure 5, the transition position (TP) and surface pressure distribution obtained from the transition model and experiment data illustrate good consistency.Figure 4
The computation structured grid: (a) grid around the airfoil; (b) leading edge grid.
(a)(b)Figure 5
The aerodynamic characteristics of the NLR 7301 obtained from CFD and experiment: (a) AoA effect on TP of upper surface; (b) AoA effect on TP of lower surface; (c) surface pressure distribution when AoA is equal to −2; (d) surface pressure distribution when AoA is equal to 2.
(a)(b)(c)(d)Also, the validation of capturing LSB is accomplished on E387 airfoil at low speed andRe=2×105. The airfoil performances predicted by transition model are compared to University of Illinois Urbana-Campaign (UIUC) wind-tunnel data [52] in Figure 6. The results displayed in Figure 6(a) demonstrate that the lift-drag curve is in good agreement with the experimental one. As presented in Figure 6(b), both the location and size of LSB obtained from CFD are coincident with the experiment data.Figure 6
The aerodynamic characteristics of the E387 airfoil obtained from CFD and experiment: (a) lift-drag curves; (b) locations of the upper surface flow features.
(a)(b)Consequently, these phenomena make us to conclude that the numerical simulation technology of transition boundary layer is accurate and reliable for the aerodynamic design optimization of transonic NLF airfoil at LRN.


## 2.3. Airfoil Parameterization Method

The disturbing basis function, based on the class/shape function transformation (CST) [53], is introduced in airfoil parameterization. The parameterization method constructed here can not only minimize the errors caused by fitting the initial parameters but also preserve the advantages of the original CST parameterization method.The original CST method is shown in the following equation:(14)ζψ=ψ0.51−ψ∑i=0nAiKiψi1−ψn−i+ψζT,where ψ and ζψ are the x-coordinates and y-coordinates of the airfoil, respectively; ζT is the y-coordinate of the trailing edge point; Ai is the design variables; and Ki=Cni=n!/i!n−i!.The disturbanceΔAi is added to the initial design variable Ai, and then, equation (14) can be written as(15)ζψ′=ψ0.51−ψ∑i=0nAi+ΔAiKiψi1−ψn−i+ψζT=ζψ+Δζψ,where ζψ′ and ζψ are the separable y-coordinates of new and original airfoils and Δζψ are the disturbances of y-coordinates, which can be denoted as(16)Δζψ=ψ0.51−ψ∑i=0nΔAiKiψi1−ψn−i.According to the aforementioned equations, the newy-coordinates can be treated as two parts, i.e., the original y-coordinates (ζψ) and disturbances (Δζψ). Thus, equation (16) can be denoted as the disturbing basis function based on CST.


## 3. Aerodynamic Design Problem

### 3.1. Optimization Design Definition

As mentioned above, the size and position of LSB and transition location have a great impact on aerodynamic characteristics of the transonic NLF airfoil at LRN. So, these two factors should be considered during the design process. The classical laminar bubble separation theory [54–56] indicates that the adverse pressure gradient of the upper surfaces will promote the laminar separation, and the transition will take place in free shear layer above the airfoil surface, of which the location is downstream of the SP. Consequently, the LSB would become smaller with the forward SP and RP, and the pressure drag would decrease. Conversely, the friction drag would increase by the forward TP. It is of particular interest to point out that the friction drag decreases dramatically by the backward TP. Nevertheless, an extra length of laminar flow region would result in the excessive adverse pressure gradient at the pressure recovery section, thus yielding a larger LSB on upper surface. Accordingly, the pressure drag sharply increases because of the larger LSB. Regarding to the above analysis, it should keep a balance between these two factors in the aerodynamic design of transonic NLF airfoil at LRN.The flight condition of the propeller tip airfoil, utilized in a solar energy UAV, is transonic at LRN. As for this, the on-design point is assumed asMa=0.72andRe=1.1×105. For the original airfoil SD-8000, its CFD results show that there is a large LSB on the upper surface and the proportion of pressure drag is 75%. Therefore, the decrease in pressure drag would correspondingly lead to a reduction in total drag. And it can be achieved by weakening the LSB. However, the friction drag would increase when the TP moves forward.


### 3.2. Deterministic Optimization Results and Discussion

The following optimized mathematical model is built up for the deterministic design problem:(17)objective:minCDp,s.t.:CL=0.5,thicknessmax≥0.089,CDvoptimized−CDvoriginalCDvoriginal≤0.1,where CL is the lift coefficient, CDp is the pressure drag coefficient, thicknessmax is the maximum thickness of airfoil, and CDvoptimized and CDvoriginal are the friction drag coefficients of the optimized and original airfoils, respectively. The objective of this model is to reduce the pressure drag coefficient. There are some constraints, including the growing proportion of friction drag coefficient less than 10%, the maximum thickness of airfoil greater than 8.92%, and the lift coefficient being 0.5.There are five design variables for the upper or lower surfaces, a total of ten design variables for each airfoil. The cell number of the grids is760×179 [57], of which the distributions for the leading edge and trailing edge are demonstrated in Figures 7(a) and 7(b). With respect to the grid convergence for the original airfoil presented in Figure 7(c), when the circumferential grid number is more than 600, the variations of the total and pressure drags are within 0.0008.Figure 7
(a) The grid distributions of the leading edge; (b) the grid distributions of the trailing edge; (c) variation of drag with circumferential grid number.
(a)(b)(c)In the deterministic optimization, the algorithm population in each iteration is set as 32 particles, and the maximum iteration step is set as 60. Table3 lists the aerodynamic characteristics for both the original and optimized airfoils. It can be summarized that the total drag and pitching moment coefficients separately decrease by 48 counts and 26.6%, while the friction drag coefficient and lift-to-drag ratio increase by 7.4% and 28.7% at on-design point, respectively. The friction drag increases because of the TP moving forward, and the pressure drag reduces because of the smaller size of LSB, as displayed in Figure 8(b). From the geometric characteristics enumerated in Table 4, it can be noticed that there is a backward location of the maximum thickness on optimized airfoil. Its maximum camber is smaller than that of the original one. And the same situation happens to the leading edge radius, of which the value for optimized airfoil is also smaller.Table 3
Aerodynamic characteristics for the original and optimized airfoils at on-design point.
CLCDK = CL/CDCDpCDvCmXtriOriginal airfoil0.50.022222.50.01680.0054−0.07850.598Optimized airfoil0.50.017428.70.01160.0058−0.05760.488Figure 8
(a) The original and optimized airfoils; (b) surface pressure distributions of original and optimized airfoils at on-design point.
(a)(b)Table 4
Geometric characteristics for the original and optimized airfoils.
Maximum thicknessLocation of maximum thicknessMaximum camberLocation of maximum camberLeading edge radiusOriginal airfoil0.0890.2740.01970.4910.0078Optimized airfoil0.0890.3060.01680.6970.0035Figure8 displayed two airfoils and surface pressure distributions at on-design point. Compared with the original airfoil, the optimized one with a flat upper surface has a smaller camber before the 60% chord and larger camber after the 60% chord. There is a stronger inverse pressure gradient roughly existing in 5% to 20% chord for the optimized pressure distribution, and then, this gradient becomes weaker for 50% to 75% chord. Moreover, the suction peak of the optimized airfoil obviously becomes larger. As marked in Figure 8, the SP and RP of the optimized airfoil are moved forward at the position of 15% and 60% chord, respectively. Figure 9 shows the Mach number distribution contours and the streamlines of the original and optimized airfoils at on-design point. The separation bubble on upper surface moves forward and becomes smaller compared with that of the original one. We assume that the transition will take place in free shear layer above the airfoil surface, of which the location is the downstream of SP. The forward TP results in a slight increase of friction drag. Simultaneously, the small LSB results in a considerable decrease of pressure drag. Moreover, the height of the bubble is also smaller. It can be concluded from Figures 8 and 9 that the strong inverse pressure gradient on the leading edge could lead to the forward SP, and the early recover pressure could result in the forward RP. The forward SP and RP minimize the LSB to reduce the pressure drag. Nevertheless, the forward LSB may be conducive to the formation of a forward transition position TP, subsequently causing the increase of friction drag. The lift-to-drag ratio curves of the original and optimized airfoils are shown in Figure 10. For the optimized airfoil, there is a low drag region over the range of lift coefficient from 0.3 to 0.52. It is caused by the sharp decrease in the pressure drag, as clearly presented in Figure 10.Figure 9
The Mach number contours and streamlines at on-design point: (a) original airfoil; (b) optimized airfoil.
(a)(b)Figure 10
The lift-to-drag ratio curves of the original and optimized airfoil (Re=1.1×105,CL=0.5): (a) lift to total drag; (b) lift to pressure drag.
(a)(b)It can be seen from Figure10 that the drag of optimized airfoil is larger than that for the original one, when CL = 0.2. By comparing the streamline distributions of the off-design point (CL = 0.2) for original airfoil in Figure 11(a), it is obviously noticed that there is a large flow separation in wake region of the optimized airfoil, which will increase the pressure drag in Figure 11(b). Furthermore, its skin friction coefficient increases in the trailing edge of the optimized airfoil in Figure 11(c), which will significantly give rise to the growth of the friction drag. Therefore, both the increases of pressure and friction drags result in increasing the total drag.Figure 11
Off-design point (CL = 0.2): (a) the streamlines of original airfoil; (b) the streamlines of optimized airfoil; (c) the skin friction distributions.
(a)(b)(c)Figure12 shows the comparison of aerodynamic characteristics between the original and optimized ones. The pressure drag of optimized airfoil shows a considerable decrease over the Mach number ranging from 0.70 to 0.76. The friction drag of optimized airfoil shows a slight increase, which meets the constraints at on-design point. Note that the maximum increase of the friction drag is only up to 5 counts. It can be summarized from the drag divergence curves that the total drag becomes decreased significantly with the substantial decrease of the pressure drag. The optimized airfoil has the smaller drag coefficients over the range of Mach numbers from 0.70 to 0.76. The SP and RP varying with the Mach numbers are displayed in Figure 12(d). Over the Mach number ranging from 0.70 to 0.76, the LSB of the optimized airfoil is smaller than that of the original one and thus leads to a larger reduction of the pressure drag. However, when Ma > 0.75, the RP cannot be detected on the original airfoil, leading to the sharp rise of the pressure drag. Therefore, the varying trend of the pressure drag is mostly influenced by the size and position of LSB over the studied Mach number range.Figure 12
The comparison of aerodynamic characteristics between the original and optimized airfoils: (a) pressure drag divergence curves; (b) friction drag divergence curves; (c) drag divergence curves; (d) SP and RP.
(a)(b)(c)(d)Figure13(a) shows the pressure distributions of original and optimized airfoils varying with the Mach number. It can be seen that the obvious flat pressure distribution appears around the position of separation bubble. The results also present that the strong inverse pressure gradient on the leading edge leads to a forward SP on the optimized airfoil and an inverse pressure gradient around the RP. These characteristics indicate that the laminar separation of optimized airfoil coincides with the classical LSB theory. The SP of original airfoil moves forward when the Mach number becomes larger, but the RP is difficult to be detected on the upper surface. These characteristics are consistent with the trailing edge separation bubble theory.Figure 13
The surface pressure distribution varying with Mach number: (a) original airfoil; (b) optimized airfoil.
(a)(b)


### 3.3. Multi-objective Robust Results and Discussion

Although the optimized airfoil of deterministic has an outstanding performance at on-design point, the pressure and total drags both have the sharp increases, when Ma > 0.75. Therefore, the fluctuation of flight conditions would cause a sudden drop in the airfoil’s efficiency. A multiobjective robust design is carried out in this section aiming to reduce the performances’ sensitivity to the uncertainties of flight conditions. The PCE method is employed to obtain the mean and variance of the objectives. The four Gaussian–Hermite integration points are introduced in the nonintrusive PCE method to calculate the PCE coefficients, and then, a 3-order PCE model is set up for the robust design. The population of MOPSO is set as 40, and the maximum iteration is 20.The Pareto optimal solution is chosen to deal with the mean and variance as two objectives rather than converting them into a single objective. We assume Mach number to be normally distributedMa∞∈N0.72,0.022. The following optimized mathematical model is built up for the multiobjective design problem:(18)objectives:minμCDp,σCDp,s.t.:CL=0.5,thicknessmax≥0.089,μCDv−μCDvoriginalμCDvoriginal≤0.1,where μCDp and σCDp are the mean and the variance of the pressure drag, respectively. μCDv and μCDvoriginal are the separate mean friction drags of the optimized and original ones.We focus on the marked individuals of the Pareto-optimal front named as optimization_2, as shown in Figure14. The deterministic optimized airfoil named as optimization_1 also plots in Figure 14. The optimization_2 has a smaller camber before the 40% chord compared with optimization_1. However, the camber distributions for these two airfoils are aligned with each other after 40% chord.Figure 14
(a) Pareto-optimal front; (b) optimized airfoils.
(a)(b)Table5 lists the aerodynamic force coefficients of the optimization ones at on-design point. It can be summarized that the total and pressure drag coefficients of optimization_2 are smaller than those of optimization_1. The lift-to-drag ratio of optimization_2 is better than that of optimization_1. Table 6 displays the geometric characteristics for two optimized airfoils. Compared with these of optimization_1, the positions of both the maximum thickness and camber of optimization_2 move forward slightly. And its leading edge radius becomes smaller demonstrated in Table 6.Table 5
Aerodynamic force coefficients for the optimized airfoils at on-design point.
CLCDK = CL/CDCDpCDvCmXtriOptimization_10.50.017428.70.01160.0058−0.05760.488Optimization_20.50.017129.20.01150.0056−0.05520.476Table 6
Geometric characteristics for the optimized airfoils.
Maximum thicknessLocation of maximum thicknessMaximum camberLocation of maximum camberLeading edge radiusOptimization_10.0890.3060.01680.6970.00351Optimization_20.0890.3040.01700.6960.00345Figure15 shows the comparisons of aerodynamic characteristics between deterministic and robust optimizations. When Ma < 0.72, optimization_2 has the larger pressure drag coefficients than the optimization_1. However, the contradictory trend is found for them, when Ma > 0.72. With a tradeoff of the pressure drag, the varying trend of pressure drag becomes gentler over the range of Mach number from 0.70 to 0.76. The friction drag coefficients of optimization_2 are smaller than the those for optimization_1. From the drag divergence curves of the optimized airfoils, we can see that the rapid change of drag calms down through the robust design. For the performance of optimization_2, although it has a trade-off of the pressure drag, its friction drag acquires a slighter drop. Therefore, the airfoil has a wide range of lower drag. The SP and RP varying with the Mach numbers are displayed in Figure 15(d). Two optimizations have the similar LSB length; however, the SP and RP of optimization_2 are both forward compared with those of optimization_1.Figure 15
The comparison of aerodynamic characteristics between the deterministic and robust optimizations: (a) pressure drag divergence curves; (b) friction drag divergence curves; (c) drag divergence curves; (d) SP and RP.
(a)(b)(c)(d)Figure16 shows the aerodynamic characteristics of optimized airfoils at on-design point. There is a stronger inverse pressure gradient occurring in 5% to 20% chord on optimization_2, and then, a weaker inverse pressure gradient can be found from 40% to 60% chord. Furthermore, the suction peak of the optimization_2 becomes larger than that of optimization_1. For clarity, the SPs and RPs are also marked in Figure 16. The SP and RP of optimization_2 would move forward at the positions of 11% and 58% chord, respectively. As mentioned before, the identical circumstance also happens here; that is, the strong inverse pressure gradient on the leading edge and the forward SP both lead to a forward TP. The position and length of LSB also can be seen from surface friction distributions, as shown in Figure 16. The separation bubble of optimization_2 on upper surface moves forward and becomes smaller than that of optimization_1. Figure 17 displays the aerodynamic characteristics of optimized airfoils at Re=1.1×105andCL=0.5. There is an obvious lower drag region of optimization_2 over the range of lift coefficient from 0.1 to 0.55. Additionally, with the Mach number increasing, the suction peak becomes smaller and the pressure recovery is postponed.Figure 16
The aerodynamic characteristics of optimized airfoils at on-design point: (a) surface pressure distributions; (b) surface friction distributions.
(a)(b)Figure 17
The aerodynamic characteristics of optimized airfoils atRe=1.1×105andCL=0.5: (a) lift-to-drag coefficient curves; (b) the surface pressure distribution of optimization_2 varying with Mach number.
(a)(b)The pitching moment curves of the original and optimized airfoils are shown in Figure18. It can be seen that the pitching moments of the optimized airfoils are smaller than that of the original one. Additionally, as shown in Figures 8, 14, and 16, both the maximum camber and after load of the optimized airfoils are also smaller than those of the original one, and thereby, they all will result in the decrease in the pitching moment.Figure 18
The comparison of moment coefficient between the original airfoil and optimized airfoils: (a) varying with Mach number; (b) varying with AOA.
(a)(b)Regarding to all above discussions, the mulit-objective robust design could produce an excellent airfoil. The optimized surface pressure distribution not only keeps the low drag at on-design point but also catches the robustness of low drag to overcome the fluctuation of flight condition.


## 3.1. Optimization Design Definition

As mentioned above, the size and position of LSB and transition location have a great impact on aerodynamic characteristics of the transonic NLF airfoil at LRN. So, these two factors should be considered during the design process. The classical laminar bubble separation theory [54–56] indicates that the adverse pressure gradient of the upper surfaces will promote the laminar separation, and the transition will take place in free shear layer above the airfoil surface, of which the location is downstream of the SP. Consequently, the LSB would become smaller with the forward SP and RP, and the pressure drag would decrease. Conversely, the friction drag would increase by the forward TP. It is of particular interest to point out that the friction drag decreases dramatically by the backward TP. Nevertheless, an extra length of laminar flow region would result in the excessive adverse pressure gradient at the pressure recovery section, thus yielding a larger LSB on upper surface. Accordingly, the pressure drag sharply increases because of the larger LSB. Regarding to the above analysis, it should keep a balance between these two factors in the aerodynamic design of transonic NLF airfoil at LRN.The flight condition of the propeller tip airfoil, utilized in a solar energy UAV, is transonic at LRN. As for this, the on-design point is assumed asMa=0.72andRe=1.1×105. For the original airfoil SD-8000, its CFD results show that there is a large LSB on the upper surface and the proportion of pressure drag is 75%. Therefore, the decrease in pressure drag would correspondingly lead to a reduction in total drag. And it can be achieved by weakening the LSB. However, the friction drag would increase when the TP moves forward.


## 3.2. Deterministic Optimization Results and Discussion

The following optimized mathematical model is built up for the deterministic design problem:(17)objective:minCDp,s.t.:CL=0.5,thicknessmax≥0.089,CDvoptimized−CDvoriginalCDvoriginal≤0.1,where CL is the lift coefficient, CDp is the pressure drag coefficient, thicknessmax is the maximum thickness of airfoil, and CDvoptimized and CDvoriginal are the friction drag coefficients of the optimized and original airfoils, respectively. The objective of this model is to reduce the pressure drag coefficient. There are some constraints, including the growing proportion of friction drag coefficient less than 10%, the maximum thickness of airfoil greater than 8.92%, and the lift coefficient being 0.5.There are five design variables for the upper or lower surfaces, a total of ten design variables for each airfoil. The cell number of the grids is760×179 [57], of which the distributions for the leading edge and trailing edge are demonstrated in Figures 7(a) and 7(b). With respect to the grid convergence for the original airfoil presented in Figure 7(c), when the circumferential grid number is more than 600, the variations of the total and pressure drags are within 0.0008.Figure 7
(a) The grid distributions of the leading edge; (b) the grid distributions of the trailing edge; (c) variation of drag with circumferential grid number.
(a)(b)(c)In the deterministic optimization, the algorithm population in each iteration is set as 32 particles, and the maximum iteration step is set as 60. Table3 lists the aerodynamic characteristics for both the original and optimized airfoils. It can be summarized that the total drag and pitching moment coefficients separately decrease by 48 counts and 26.6%, while the friction drag coefficient and lift-to-drag ratio increase by 7.4% and 28.7% at on-design point, respectively. The friction drag increases because of the TP moving forward, and the pressure drag reduces because of the smaller size of LSB, as displayed in Figure 8(b). From the geometric characteristics enumerated in Table 4, it can be noticed that there is a backward location of the maximum thickness on optimized airfoil. Its maximum camber is smaller than that of the original one. And the same situation happens to the leading edge radius, of which the value for optimized airfoil is also smaller.Table 3
Aerodynamic characteristics for the original and optimized airfoils at on-design point.
CLCDK = CL/CDCDpCDvCmXtriOriginal airfoil0.50.022222.50.01680.0054−0.07850.598Optimized airfoil0.50.017428.70.01160.0058−0.05760.488Figure 8
(a) The original and optimized airfoils; (b) surface pressure distributions of original and optimized airfoils at on-design point.
(a)(b)Table 4
Geometric characteristics for the original and optimized airfoils.
Maximum thicknessLocation of maximum thicknessMaximum camberLocation of maximum camberLeading edge radiusOriginal airfoil0.0890.2740.01970.4910.0078Optimized airfoil0.0890.3060.01680.6970.0035Figure8 displayed two airfoils and surface pressure distributions at on-design point. Compared with the original airfoil, the optimized one with a flat upper surface has a smaller camber before the 60% chord and larger camber after the 60% chord. There is a stronger inverse pressure gradient roughly existing in 5% to 20% chord for the optimized pressure distribution, and then, this gradient becomes weaker for 50% to 75% chord. Moreover, the suction peak of the optimized airfoil obviously becomes larger. As marked in Figure 8, the SP and RP of the optimized airfoil are moved forward at the position of 15% and 60% chord, respectively. Figure 9 shows the Mach number distribution contours and the streamlines of the original and optimized airfoils at on-design point. The separation bubble on upper surface moves forward and becomes smaller compared with that of the original one. We assume that the transition will take place in free shear layer above the airfoil surface, of which the location is the downstream of SP. The forward TP results in a slight increase of friction drag. Simultaneously, the small LSB results in a considerable decrease of pressure drag. Moreover, the height of the bubble is also smaller. It can be concluded from Figures 8 and 9 that the strong inverse pressure gradient on the leading edge could lead to the forward SP, and the early recover pressure could result in the forward RP. The forward SP and RP minimize the LSB to reduce the pressure drag. Nevertheless, the forward LSB may be conducive to the formation of a forward transition position TP, subsequently causing the increase of friction drag. The lift-to-drag ratio curves of the original and optimized airfoils are shown in Figure 10. For the optimized airfoil, there is a low drag region over the range of lift coefficient from 0.3 to 0.52. It is caused by the sharp decrease in the pressure drag, as clearly presented in Figure 10.Figure 9
The Mach number contours and streamlines at on-design point: (a) original airfoil; (b) optimized airfoil.
(a)(b)Figure 10
The lift-to-drag ratio curves of the original and optimized airfoil (Re=1.1×105,CL=0.5): (a) lift to total drag; (b) lift to pressure drag.
(a)(b)It can be seen from Figure10 that the drag of optimized airfoil is larger than that for the original one, when CL = 0.2. By comparing the streamline distributions of the off-design point (CL = 0.2) for original airfoil in Figure 11(a), it is obviously noticed that there is a large flow separation in wake region of the optimized airfoil, which will increase the pressure drag in Figure 11(b). Furthermore, its skin friction coefficient increases in the trailing edge of the optimized airfoil in Figure 11(c), which will significantly give rise to the growth of the friction drag. Therefore, both the increases of pressure and friction drags result in increasing the total drag.Figure 11
Off-design point (CL = 0.2): (a) the streamlines of original airfoil; (b) the streamlines of optimized airfoil; (c) the skin friction distributions.
(a)(b)(c)Figure12 shows the comparison of aerodynamic characteristics between the original and optimized ones. The pressure drag of optimized airfoil shows a considerable decrease over the Mach number ranging from 0.70 to 0.76. The friction drag of optimized airfoil shows a slight increase, which meets the constraints at on-design point. Note that the maximum increase of the friction drag is only up to 5 counts. It can be summarized from the drag divergence curves that the total drag becomes decreased significantly with the substantial decrease of the pressure drag. The optimized airfoil has the smaller drag coefficients over the range of Mach numbers from 0.70 to 0.76. The SP and RP varying with the Mach numbers are displayed in Figure 12(d). Over the Mach number ranging from 0.70 to 0.76, the LSB of the optimized airfoil is smaller than that of the original one and thus leads to a larger reduction of the pressure drag. However, when Ma > 0.75, the RP cannot be detected on the original airfoil, leading to the sharp rise of the pressure drag. Therefore, the varying trend of the pressure drag is mostly influenced by the size and position of LSB over the studied Mach number range.Figure 12
The comparison of aerodynamic characteristics between the original and optimized airfoils: (a) pressure drag divergence curves; (b) friction drag divergence curves; (c) drag divergence curves; (d) SP and RP.
(a)(b)(c)(d)Figure13(a) shows the pressure distributions of original and optimized airfoils varying with the Mach number. It can be seen that the obvious flat pressure distribution appears around the position of separation bubble. The results also present that the strong inverse pressure gradient on the leading edge leads to a forward SP on the optimized airfoil and an inverse pressure gradient around the RP. These characteristics indicate that the laminar separation of optimized airfoil coincides with the classical LSB theory. The SP of original airfoil moves forward when the Mach number becomes larger, but the RP is difficult to be detected on the upper surface. These characteristics are consistent with the trailing edge separation bubble theory.Figure 13
The surface pressure distribution varying with Mach number: (a) original airfoil; (b) optimized airfoil.
(a)(b)


## 3.3. Multi-objective Robust Results and Discussion

Although the optimized airfoil of deterministic has an outstanding performance at on-design point, the pressure and total drags both have the sharp increases, when Ma > 0.75. Therefore, the fluctuation of flight conditions would cause a sudden drop in the airfoil’s efficiency. A multiobjective robust design is carried out in this section aiming to reduce the performances’ sensitivity to the uncertainties of flight conditions. The PCE method is employed to obtain the mean and variance of the objectives. The four Gaussian–Hermite integration points are introduced in the nonintrusive PCE method to calculate the PCE coefficients, and then, a 3-order PCE model is set up for the robust design. The population of MOPSO is set as 40, and the maximum iteration is 20.The Pareto optimal solution is chosen to deal with the mean and variance as two objectives rather than converting them into a single objective. We assume Mach number to be normally distributedMa∞∈N0.72,0.022. The following optimized mathematical model is built up for the multiobjective design problem:(18)objectives:minμCDp,σCDp,s.t.:CL=0.5,thicknessmax≥0.089,μCDv−μCDvoriginalμCDvoriginal≤0.1,where μCDp and σCDp are the mean and the variance of the pressure drag, respectively. μCDv and μCDvoriginal are the separate mean friction drags of the optimized and original ones.We focus on the marked individuals of the Pareto-optimal front named as optimization_2, as shown in Figure14. The deterministic optimized airfoil named as optimization_1 also plots in Figure 14. The optimization_2 has a smaller camber before the 40% chord compared with optimization_1. However, the camber distributions for these two airfoils are aligned with each other after 40% chord.Figure 14
(a) Pareto-optimal front; (b) optimized airfoils.
(a)(b)Table5 lists the aerodynamic force coefficients of the optimization ones at on-design point. It can be summarized that the total and pressure drag coefficients of optimization_2 are smaller than those of optimization_1. The lift-to-drag ratio of optimization_2 is better than that of optimization_1. Table 6 displays the geometric characteristics for two optimized airfoils. Compared with these of optimization_1, the positions of both the maximum thickness and camber of optimization_2 move forward slightly. And its leading edge radius becomes smaller demonstrated in Table 6.Table 5
Aerodynamic force coefficients for the optimized airfoils at on-design point.
CLCDK = CL/CDCDpCDvCmXtriOptimization_10.50.017428.70.01160.0058−0.05760.488Optimization_20.50.017129.20.01150.0056−0.05520.476Table 6
Geometric characteristics for the optimized airfoils.
Maximum thicknessLocation of maximum thicknessMaximum camberLocation of maximum camberLeading edge radiusOptimization_10.0890.3060.01680.6970.00351Optimization_20.0890.3040.01700.6960.00345Figure15 shows the comparisons of aerodynamic characteristics between deterministic and robust optimizations. When Ma < 0.72, optimization_2 has the larger pressure drag coefficients than the optimization_1. However, the contradictory trend is found for them, when Ma > 0.72. With a tradeoff of the pressure drag, the varying trend of pressure drag becomes gentler over the range of Mach number from 0.70 to 0.76. The friction drag coefficients of optimization_2 are smaller than the those for optimization_1. From the drag divergence curves of the optimized airfoils, we can see that the rapid change of drag calms down through the robust design. For the performance of optimization_2, although it has a trade-off of the pressure drag, its friction drag acquires a slighter drop. Therefore, the airfoil has a wide range of lower drag. The SP and RP varying with the Mach numbers are displayed in Figure 15(d). Two optimizations have the similar LSB length; however, the SP and RP of optimization_2 are both forward compared with those of optimization_1.Figure 15
The comparison of aerodynamic characteristics between the deterministic and robust optimizations: (a) pressure drag divergence curves; (b) friction drag divergence curves; (c) drag divergence curves; (d) SP and RP.
(a)(b)(c)(d)Figure16 shows the aerodynamic characteristics of optimized airfoils at on-design point. There is a stronger inverse pressure gradient occurring in 5% to 20% chord on optimization_2, and then, a weaker inverse pressure gradient can be found from 40% to 60% chord. Furthermore, the suction peak of the optimization_2 becomes larger than that of optimization_1. For clarity, the SPs and RPs are also marked in Figure 16. The SP and RP of optimization_2 would move forward at the positions of 11% and 58% chord, respectively. As mentioned before, the identical circumstance also happens here; that is, the strong inverse pressure gradient on the leading edge and the forward SP both lead to a forward TP. The position and length of LSB also can be seen from surface friction distributions, as shown in Figure 16. The separation bubble of optimization_2 on upper surface moves forward and becomes smaller than that of optimization_1. Figure 17 displays the aerodynamic characteristics of optimized airfoils at Re=1.1×105andCL=0.5. There is an obvious lower drag region of optimization_2 over the range of lift coefficient from 0.1 to 0.55. Additionally, with the Mach number increasing, the suction peak becomes smaller and the pressure recovery is postponed.Figure 16
The aerodynamic characteristics of optimized airfoils at on-design point: (a) surface pressure distributions; (b) surface friction distributions.
(a)(b)Figure 17
The aerodynamic characteristics of optimized airfoils atRe=1.1×105andCL=0.5: (a) lift-to-drag coefficient curves; (b) the surface pressure distribution of optimization_2 varying with Mach number.
(a)(b)The pitching moment curves of the original and optimized airfoils are shown in Figure18. It can be seen that the pitching moments of the optimized airfoils are smaller than that of the original one. Additionally, as shown in Figures 8, 14, and 16, both the maximum camber and after load of the optimized airfoils are also smaller than those of the original one, and thereby, they all will result in the decrease in the pitching moment.Figure 18
The comparison of moment coefficient between the original airfoil and optimized airfoils: (a) varying with Mach number; (b) varying with AOA.
(a)(b)Regarding to all above discussions, the mulit-objective robust design could produce an excellent airfoil. The optimized surface pressure distribution not only keeps the low drag at on-design point but also catches the robustness of low drag to overcome the fluctuation of flight condition.


## 4. Conclusions

In this paper, we carry out the aerodynamic design optimization of a propeller tip airfoil in a solar energy UAV. The flight condition of this airfoil is treated as transonic at LRN. Under such conditions, we develop an optimization methodology based on the PSO algorithm applied to the design of a transonic NLF airfoil at LRN. This new proposed methodology briefly contains two optimization design processes. They are traditional deterministic optimization at on-design point on the basis of MCPSO and multi-objective robust design based on MOPSO. The non-intrusive PCE method is adopted to get sufficiently accurate estimations of the mean and variance. It facilitates to the fact that the optimal deterministic solution and Pareto-optimal front of aerodynamic design can be captured in the space quickly and accurately. From the findings of deterministic optimizations and multi-objective robust design, we obtained some important conclusions for transonic NLF airfoil at LRN, as follows:(1)
We construct a new optimization model; that is, the minimized pressure drag is typically treated as the optimization objective, while the friction drag increment is a constraint condition. According to the deterministic optimization and robust design results, this proposed model is reliable and efficacious to design a high-performance airfoil of the propeller tip.(2)
The pressure drag is mostly affected by the size and position of LSB. The deterministic optimized pressure distribution can lead to a small and forward LSB. The SP and RP of the airfoil move forward by virtue of the deterministic optimized pressure distribution.(3)
The multi-objective robust design can weaken the sensitivity of airfoil performance to the flight fluctuation. The proper inverse pressure gradient at the leading edge and early pressure recover show the good abilities to control the LSB.(4)
The proposed optimization methodology is an efficient and reasonable approach for the aerodynamic design of transonic NLF airfoil at LRN.



---

*Source: 5812129-2020-07-16.xml*