# Aerodynamic Optimization of a UAV Wing subject to Weight, Geometric, Root Bending Moment, and Performance Constraints

Authors: Durmuş Sinan Körpe; Öztürk Özdemir Kanat
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

In this study, the optimization of a low-speed wing with functional constraints is discussed. The aerodynamic analysis tool developed by the coupling of the numerical nonlinear lifting-line method to Xfoil is used to obtain lift and drag coefficients of the baseline wing. The outcomes are compared with the results of the solver based on the nonlinear lifting-line theory implemented into XLFR5 and the transition shear stress transport model implemented into ANSYS-Fluent. The agreement between the results at the low and moderate angle of attack values is observed. The sequential quadratic programming algorithm of the MATLAB optimization toolbox is used for the solution of the constrained optimization problems. Three different optimization problems are solved. In the first problem, the maximization ofCL3/2/CD is the objective function, while level flight condition at maximum CL3/2/CD is defined as a constraint. The functional constraints related to the wing weight, the wing planform area, and the root bending moment are added to the first optimization problem, and the second optimization problem is constructed. The third optimization problem is obtained by adding the level flight condition and the available power constraints at the maximum speed and the level flight condition at the minimum speed of the baseline unmanned air vehicle to the second problem. It is demonstrated that defining the root bending moment, the wing area, and the available power constraints in the aerodynamic optimization problems leads to more realistic wing planform and airfoil shapes.


---

## Body

## 1. Introduction

Recent advances in computer technology have advanced the conceptual design of an aircraft using advanced optimization tools rather than rough sketches and calculations on paper. In the early years of aviation, a small group of engineers, led by a design engineer who had experience in both design and manufacturing, were responsible for the conceptual design. However, at present, the design group has expanded and divided into subgroups because of the complexity of the aircraft, which is a result of the technological advances [1]. Simultaneously, faster solutions of the optimization problems have enabled the definition of optimization problems that have objective functions and constraints defined by different disciplines to satisfy the expected performance merits of the aircraft in the conceptual design. In addition to satisfying the performance merits, developing a mathematical proof of the optimized aircraft at the end of the design phase is crucial. To obtain accurate results in the optimization problems, the function evaluator should be a verified and widely used analytical or numerical analysis tool.Xfoil, developed by Mark Drela, is utilized as the airfoil analysis code in the numerous studies that are related to the design of fixed and rotary wings, wind turbine blades, and marine propellers [2]. Aerodynamic coefficients are obtained by the coupling of the panel method to the two-equation integral boundary-layer formulation that is the integral momentum and kinetic energy shape parameter equations. Falkner-Skan and Swafford velocity profile formulas are used for laminar closure and turbulent closure, respectively. The en method is implemented for transition point detection [3].Koreanschi et al. used Xfoil for the evaluation of the base airfoil performance in their study that presented numerical optimization and experimental wind tunnel testing of a morphing wing tip equipped with an adaptable upper surface, and a rigid aileron [4]. Gabor et al. used Xfoil in the high angle of attack α optimization of a morphing airfoil as a function evaluator [5]. Della Vecchia et al. investigated the effect of a morphing trailing edge to the performance of a high-altitude long-endurance (HALE) aircraft. The variation of the lift coefficient Cl and the profile drag coefficient Cd of the airfoils with trailing edge deflection was obtained with Xfoil [6]. Magrini and Benini investigated the aerodynamic optimization of a morphing leading-edge airfoil. The optimization problem was solved by using the transition SST model, Spalart-Allmaras model, and Xfoil [7]. Liu et al. studied the optimization of the airfoil at an ultra-low Reynolds number Re for nanorotor performance [8]. In the study, Xfoil results were compared with the results of a two-dimensional incompressible Navier-Stokes solver in which the artificial compressibility method was utilized to deal with incompressible flow. In this validation study, the National Advisory Committee for Aeronautics (NACA) 0006 airfoil was used and Re was 6000. Xfoil captured Cl and Cd values up to stall α correctly. Silvestre et al. coupled the blade element momentum theory with Xfoil and developed a new propeller design code named as JBLADE [9].The nonlinear numerical lifting-line method proposed by Anderson and Corda is an iterative solution to Prandtl’s lifting-line theory [10]. Gamboa et al. used the nonlinear lifting-line method as a function evaluator for the morphing wing optimization. In the study, airfoil aerodynamic data was provided by Xfoil [11]. Merz et al. used the numerical nonlinear lifting-line method to find the most effective α on a pitching rotor blade [12].


## 2. Method

First of all, this paper describes the coupling of the nonlinear numerical lifting-line method to Xfoil. Then, the results of the newly designed aerodynamic analysis tool are compared with the results of the solver based on the nonlinear lifting-line theory (LLT) implemented into XLFR5 and the transitionγ−Reθ shear stress transport (SST) model implemented into ANSYS-Fluent. After that, three different optimization problems in which CL3/2/CD maximization is the objective function are constructed. The optimization problems are defined by adding new functional constraints that are the wing root bending moment Mb, the wing weight Ww, the wing planform area S, the level flight condition and power available PA at the maximum speed Vmax, and the level flight condition at the stall speed Vstall of the baseline unmanned air vehicle UAV. The change of the objective function, the constraints, the wing planform, and the airfoil shapes is discussed. This study differs from the studies in the literature because the performance of the UAV at both Vmax and Vstall conditions at the level flight is taken into account with the functional constraints in the optimization problems that are discussed below.


### 2.1. The Nonlinear Numerical Lifting-Line Theory

In this part of the study, the methodology of the nonlinear numerical lifting-line theory and integration of Xfoil to it are discussed. The numerical solution starts with the division of the spanb of the finite wing into N number of spanwise stations in the y-direction as shown in Figure 1 [13]. In Figure 1, the distance between the equidistant stations is Δy. The method needs an initial circulation distribution Γy along the span. The elliptic circulation distribution in equation (1) is acceptable for this requirement.
(1)Γy=Γ01−2yb2.Figure 1
The positions of the stations.In the above equation,Γ0 is the circulation at the root chord cr. By using initial Γy, the induced angle of attack αi for the nth panel is obtained as follows:
(2)αiyn=14πV∞∫−b/2b/2dΓ/dydyyn−y.This integral is evaluated numerically by using Simpson’s rule [13]:
(3)αiyn=Δy12πV∞∑j=2,4,6N−1dΓ/dyj−1yn−yj−1+dΓ/dyjyn−yj+dΓ/dyj+1yn−yj+1,where V∞ is the freestream velocity. For the calculation of αi, dΓ/dy is needed. Although dΓ/dy is analytically available for the initial calculation [13], as shown in equation (4), it must be calculated by using a finite difference method for the other iterations as shown in equation (5).The forward difference method is selected for the process.(4)dΓdy=−4Γ0b2y1−4y2/b21/2,(5)dΓdyj=Γj+1−Γjyj+1−yj.In order to avoid singularity whenyn=yj−1, yj, or yj+1, the average of the neighbor stations’ results is taken as the result of the corresponding term [13]. After the calculation of αi, the effective angle of attack αeff is obtained by subtracting αi from the geometric angle of attack αgeo.
(6)αeffyn=αgeoyn−αiyn.At this point, it should be stated that the method is a numerical iterative solution to Prandtl’s lifting-line theory that states that the circulation is zero wheny=±b/2 [14]. That implies that Cl values of the airfoils at the tip chords ct of the wing are zero and αeff is equal to zero lift angle of attack value αCl=0 of the airfoils at the tips as shown in
(7)αeff±b2=αCl=0±b2.Onceαeff distribution along the span is known, Cl and Cdof the airfoils at the stations are calculated by using Xfoil [15]. An interface developed to run Xfoil in MATLAB is used [16]. It is modified for the MAC OSX operating system. Then, the new circulation distribution Γnew is calculated by using
(8)Γnewy=V∞cyCly2.cy in the above equation is the chord value at the station. The input circulation values are calculated by using a damping factor D as follows:
(9)Γy=Γy+DΓnewy−Γy.The iterative solution requires a strong damping factor. SettingD values as one leads to divergent iteration [10]. The iteration process stops if the Euclidean norm of the difference between the new circulation vector and the input circulation vector is less than a convergence criterion ε:
(10)Γnewy−Γy<ε.By using converged circulation distribution, the liftL [10], the induced drag Di [10], and the root bending moment Mb [17] of the wing are calculated as follows:
(11)L=ρ∞V∞∫−b/2b/2Γnydy,(12)Di=ρ∞V∞∫−b/2b/2Γnysinαinydy,(13)Mb=ρ∞V∞∫0b/2Γnyydy,where ρ∞ is the free stream density. These integrals are also evaluated numerically by using Simpson’s rule. The parasite drag of the wing Dpr is calculated by using
(14)Dpr=12ρ∞V∞2Δy4∑j=1N−1cnj+cnj+1Cdj+Cdj+1.After the calculations of the aerodynamic forces and moments, the lift coefficientCL, the drag coefficient CDw, and the root bending moment coefficient CMb of the wing are obtained as follows:
(15)CL=2Lρ∞V∞2S,(16)CDw=2Dpr+Diρ∞V∞2S,(17)CMb=4Mbρ∞V∞2Sb.As a result, Xfoil and the numerical nonlinear lifting-line method are coupled and the aerodynamic analysis tool is obtained.


### 2.2. XFLR5 LLT and ANSYS-Fluentγ−Reθ SST Model

The XFLR5 LLT and ANSYS-Fluentγ−Reθ SST model are introduced in this section. XFLR5 is an analysis tool for airfoils and wings. There are four different solvers based on the solution of the inviscid potential flow implemented into it. These are the solvers based on the horseshoe vortex lattice method, the ring vortex lattice method, the 3-D panel method, and the nonlinear lifting-line theory. They were used to calculate the drag polar of the different wings [18–20].The solver based on nonlinear lifting-line theory calculates theCL and CDw results by taking the drag polar data of the airfoil form Xfoil that is also implemented into XFLR5. The results are obtained by using the default values defined from the solver except for the iteration number. It is set to 1000.Another computational fluid dynamics (CFD) tool, ANSYS-Fluent, which solves the problems according to the finite volume method, is used to compare theCL and CDw results obtained by the aerodynamic analysis tool.Langtry and Menter developed a new correlation-based transition model to simulate laminar to turbulent transition [21]. The model is essentially based on two transport equations. One equation is for a transition onset criterion and the other is for intermittency. To validate this model, they gave examples of the cases which they studied such as 2-D three-element flap, 2-D airfoil, and a transonic wing. They stated that the transitional CFD simulation was in very good agreement with the experimental results. Aftab et al. compared the following turbulence models for the flow over the NACA4415 airfoil for Re of 120000 with the experimental results: Spalart-Allmaras, K−ω SST, intermittency SST, k−kl−ω, and γ−Reθ SST [22]. According to the results, only the γ−Reθ SST provided reliable results for low and high α values. Lanzafame et al. developed a horizontal axis wind turbine 3-D CFD model to predict wind turbine performance [23]. They compared K−ω SST and γ−Reθ SST with experimental data. According to the result, γ−Reθ SST captured the trend of the aerodynamic coefficients whereas K−ω SST overpredicted Cl values and underpredicted Cd values.According to these studies in the literature, the computations are performed withγ−Reθ SST in this study. The numerical studies are carried out according to the pressure-based method. There are some specific criteria that determine the accuracy of the results obtained with ANSYS-Fluent. The first is to determine the size of the boundary region where flow analysis is intended. For this study, the size of the boundary is formed approximately 30 times larger than the maximum wing thickness and wing chord length as shown in Figure 2 [24, 25]. Since the control volumes are 3-D, the volumetric adaptation of the inflation layers and the prism cells are used. The dimensionless distance to the wall y+ is selected as 1 [26]. Thereby, first-layer height is calculated as 0.015 mm for the selected y+. The growth rate and the number of inflation layers are selected as 1.1 and 61, respectively. These selections create a mesh with 2569785 elements and with a maximum skewness value of 0.95 as shown in Figure 3. The turbulence intensity is selected as 1%. For the convergence criterion, the continuity equation is evaluated as 2·10−4 for 0° α. In order to provide convergence, approximately 1600 iterations are performed. The convergence criterion for high α values increases up to 10 times. Therefore, the results are accepted as accurate in these calculations when fluctuations in CL and CDw lie within the range of ±0.001.Figure 2
The boundaries of the solution domain.Figure 3
Cells around the airfoil atcr.


### 2.3. Comparison of the Aerodynamic Analysis Tool with ANSYS-Fluent

The comparison of the aerodynamic analysis tool with XFLR5 6.47 LLT and the ANSYS-Fluent 16.2γ−Reθ SST model is made in this part of the paper. The flow solutions of the aerodynamic analysis tool are obtained on a 2.4 GHz Intel Core i7, 8 GB MacBook Pro. The computation of ANSYS-Fluent is performed on 2 processors of 2.4 GHz Intel Xeon E5, 32 GB HP Z820.The interface is modified so thatCl and Cd values are obtained for the airfoils at each station by parallel processing. When Xfoil fails to converge, the initial panel number, which is 200, is altered between 194 and 206. NACA 4412 airfoil with closed trading edge is used as the candidate airfoil. Equations (18) and (19) are used for the half thickness zt and the camber generation zc, respectively. In these equations, t is the maximum thickness in hundredths of the chord, p is the maximum camber location in tenths of the chord, and m is the maximum camber in hundredths of the chord. x coordinates are generated by the cosine distribution.Thez coordinates of the upper surface are obtained by summing equations (18) and (19), whereas equation (18) is subtracted from equation (19) for the lower surface as is done in the NACA airfoil generation subroutine of Xfoil [27].
(18)ztx=t0.20.2969x0.5−0.126x−0.3516x2+0.2843x3−1036x4,(19)zcx=mp22px−x2,if0≤x<p,zcx=m1−p21−2p+2px−x2,ifp≤x≤1.cr, ct, and b are 0.45m, 0.45m, and 4m, respectively. ρ∞ is 1.225kg/m3. The dynamic viscosity μ∞ is 1.7974·10−5Ns/m2 and V∞ is 22 m/s. In Xfoil, the default critical n value of the en transition prediction method is nine. This value is changed as 2.62, because it corresponds to the selected turbulent intensity in the ANSYS-Fluent analysis.In Figure4, the CL results for the different panel station numbers and the results of ANSYS-Fluent and XFLR5 are shown. According to the results in Figure 4, ANSYS-Fluent calculates the least CL at the same α. CL values of N=21, N=31, and N=41 cases overlap. The results show slight difference with the results of XFLR5 LLT and capture the trend line of ANSYS-Fluent up to the stall α. CDw values of ANSYS-Fluent are almost identical to N=11 case values when α is between 0° and 10° as depicted in Figure 5. But as α gets close to 15°, the ANSYS-Fluent, N=31 case, and N=41 case yield similar CDw. After 15°, ANSYS-Fluent results start to diverge from the results of the aerodynamic tool. In order to investigate the situation in detail, the flow field around the wing for 10° and 18° α values is presented in Figures 6–9.Figure 4
CL vs. α curves of the aerodynamic design tool, ANSYS-Fluent, and XFLR5.Figure 5
CDw vs. α curves of the aerodynamic design tool, ANSYS-Fluent, and XFLR5.Figure 6
Cp contours of the aerodynamic analysis tool and ANSYS-Fluent for 10°.Figure 7
x-wall shear stress contours of the aerodynamic analysis tool and ANSYS-Fluent for 10°.Figure 8
Cp contours of the aerodynamic analysis tool and ANSYS-Fluent for 18°.Figure 9
x-wall shear stress contours of the aerodynamic analysis tool and ANSYS-Fluent for 18°.In the aerodynamic analysis tool, the transition region can be determined by finding the region at which a sudden increase in thex-wall shear stress is observed. The turbulent region is the region at which the x-wall shear stress values get closer to zero. Contrary to this, the transition region and the turbulent separation region are identified when x-wall shear stress is less than zero in the results of ANSYS-Fluent [28].According to the results for 10°, apart from the slight difference around the leading edge, pressure contours are almost identical for both solvers as shown in Figure6. The transition prediction by the aerodynamic solver is understood from the color change from green to yellow in Figure 7. ANSYS-Fluent captures the transition location almost at the same position as the aerodynamic analysis tool does that is shown by a color change from green to blue. Both solvers capture the separation around the trailing edge. When the results for 18° are discussed, one can easily see the difference between the pressure contours especially on the aft of the wing as depicted in Figure 8. As is shown in Figure 9, both solvers predict the transition region almost at the same location. ANSYS-Fluent predicts the turbulent separation ahead. But it is very close to the turbulent separation location of the aerodynamic analysis tool.SinceN=31 and N=41 cases capture the CDw values of ANSYS-Fluent at high α values, they are selected as the candidate station numbers.In order to select the optimumN, D, and ε values, the total solution time ttot to obtain drag polar results between the α values of -6° and 20° is compared for the different N, D, and ε values in Table 1. The maximum deviations in CL∆CL and CDw∆CDw are obtained with respect to the results when the N, D, and ε values are equal to 41, 0.05, and 0.01, respectively. According to the results, choosing the N value as 31 saves the computation time approximately by 38% at the same D and ε values when it is compared with the result of N=41.Table 1
Effect ofN, D, and ε to the results and ttot.
N
D
ε
∆CL(%)
∆CDw(%)
ttot(s)
41
0.05
0.01
0
0
3007.7
31
0.05
0.01
0.4
1.1
2532.9
31
0.1
0.01
0.4
1.1
1096.2
31
0.2
0.01
0.4
1.1
519.3
31
0.2
0.05
5.9
1.1
418.5
31
0.2
0.1
11.2
1.0
414.3
31
0.2
0.2
15.2
0.9
347.6IncreasingD decreases the computation time by 86.9% without boosting ∆CL and ∆CDw. When D is selected as 0.4 and above, the convergence problem is observed at the high α values. Contrary to this, increasing ε yields drastic increase in ∆CL especially at the negative α and the high α values. Apart from this region, ∆CL is around 1% at the low α and moderate α values. In light of this information, N, D, and ε values are selected as 31, 0.2, and 0.01, respectively. ttot of ANSYS- Fluent results is approximately 120 hours.


### 2.4. Optimization Solver

In this part, the advantages of the selected optimization solver and its use are discussed. Vanderplaats compared different optimization methods by using a structural optimization problem [29]. In the problem, the objective function was the minimization of the structural volume of a cantilevered beam that was fixed at the right end and had a vertical load applied at the free left end. The cantilevered beam had five segments and thicknesses, and heights of the segments were the design variables. In the study, the allowable bending stress at the right end of each segment was defined as the constraints. In addition to this, the allowable left end deflection was also defined as the constraint. The ratio of the height to the thickness of each segment was limited by using constraints. In summary, the identified optimization problem had 10 design variables and 11 constraints. Genetic search, sequential linear programming, method of feasible directions, generalized reduced gradient method, modified feasible directions method, and sequential quadratic programming (SQP) were compared in the study. The first method is an evolutionary optimization method whereas the others are the gradient-based optimization methods. The optimum results obtained by the methods were very close to each other. The genetic search method obtained the optimum result after 106 function evaluations. Among the gradient-based optimization methods, SQP had the least number of iterations and the least number of gradient evaluations. The total number of function call is 125 for the SQP analysis. As discussed above, the main advantage of SQP is that it does not attempt to satisfy constraints at each iteration that results in less number of function calls [30]. Therefore, the SQP algorithm is used as given in the MATLAB optimization toolbox for the solution of the optimization problems that are discussed below.The procedure for the optimization process is represented in Figure10.Figure 10
The flowchart of the optimization process.The design variable array is used to calculate the objective function and the constraints for the current step. For this case, the aerodynamic design tool is called in parallel computations and the threads are used to calculateCl and Cd values for the airfoils at each station. After that, the design variable array is perturbed by an increment array that is 1% of each design variable array element. Parallel computing is applied for the calculation of the objective function and the constraints by using perturbed design variable array elements. The aerodynamic design tool runs in serial for this case. The forward difference method is used for the calculation of the gradient of the objective function and Jacobian matrix of the constraints.The process stops if the change of the design variable array (TolX) or the change in the objective function (TolFun) is less than 10−3. The violation of the constraints (TolCon) is limited to 10−2.


### 2.5. The Characteristics of the Baseline UAV

Before going into the details of the optimization problems, the geometric and performance characteristics of the baseline UAV and basic assumptions related to it are revealed in this part of the work. The baseline UAV is a propeller-driven aircraft andPA of the engine is 2000 W. The drag coefficient of the other components (CDoth) of the UAV except the wing is assumed as a varying parameter as follows:
(20)CDoth=0.02SbS,where Sb is the planform area of the wing of the baseline UAV. As a result, the total drag coefficient CD of the baseline UAV is the summation of CDoth and CDw. The total weight of the UAV except for the wing (Woth) is 250 N. The wing weight Ww is calculated as follows [31]:
(21)Ww=Sc¯tcmaxρmatKρARnult0.6λ0.04g,where c¯ denotes the wing mean aerodynamic chord, t/cmax the maximum thickness of the airfoil, AR the aspect ratio, and λ the taper ratio. ρmat is the density of the construction material that is graphite/epoxy. Kρ is the wing density factor. They are selected as 1575 kg/m3 and 0.0016, respectively, by calculating the average values in the appropriate reference tables for homebuilt aircraft [31]. The relation between ultimate load factor (nult) and maximum positive load factor (nmax) is defined as follows:
(22)nult=1.5nmax.Unlike a fighter aircraft, the produciblenmax by the engine power is less than the tolerable nmax by the aircraft structure for a general aviation aircraft [32]. For this reason, the equation of the producible nmax by the engine power for the propeller-driven aircraft that is in equation (23) is used in the Ww calculation.
(23)nmax=CL3CD2PA2ρ∞S21/31Woth+Ww.As seen in the above equations,Ww and nmax are implicit. The secant method is applied to find Ww.The geometric parameters of the wing of the baseline UAV are shown in Table2.Table 2
Geometric properties of the baseline UAV’s wing.
Airfoil
NACA4412
cr
0.45 m
ct
0.45 m
b
4
Sb
1.8 m2The wing of the baseline UAV is a rectangular and unswept wing. The airfoil profile is constant along the span. Table3 depicts the performance parameters of the baseline UAV.Table 3
Performance parameters of the baseline UAV.
maxCL3/2/CD
14.08
α at maxCL3/2/CD
8°
V at maxCL3/2/CD
15.52 m/s
Mb at maxCL3/2/CD
127.46 Nm
Ww
24.06 N
Vmax
39.45 m/s
α at Vmax
-2.32°
Vstall
13.25
α at Vstall
18°It has a maximumCL3/2/CD ratio of 14.08 when α is 8° and speed is 15.52 m/s. Mb is 127.46 Nm for this flight condition. Ww is 24.06 N. Vmax is 39.45 m/s when α is -2.32° and Vstall is 13.25 m/s when α is 18°.


## 2.1. The Nonlinear Numerical Lifting-Line Theory

In this part of the study, the methodology of the nonlinear numerical lifting-line theory and integration of Xfoil to it are discussed. The numerical solution starts with the division of the spanb of the finite wing into N number of spanwise stations in the y-direction as shown in Figure 1 [13]. In Figure 1, the distance between the equidistant stations is Δy. The method needs an initial circulation distribution Γy along the span. The elliptic circulation distribution in equation (1) is acceptable for this requirement.
(1)Γy=Γ01−2yb2.Figure 1
The positions of the stations.In the above equation,Γ0 is the circulation at the root chord cr. By using initial Γy, the induced angle of attack αi for the nth panel is obtained as follows:
(2)αiyn=14πV∞∫−b/2b/2dΓ/dydyyn−y.This integral is evaluated numerically by using Simpson’s rule [13]:
(3)αiyn=Δy12πV∞∑j=2,4,6N−1dΓ/dyj−1yn−yj−1+dΓ/dyjyn−yj+dΓ/dyj+1yn−yj+1,where V∞ is the freestream velocity. For the calculation of αi, dΓ/dy is needed. Although dΓ/dy is analytically available for the initial calculation [13], as shown in equation (4), it must be calculated by using a finite difference method for the other iterations as shown in equation (5).The forward difference method is selected for the process.(4)dΓdy=−4Γ0b2y1−4y2/b21/2,(5)dΓdyj=Γj+1−Γjyj+1−yj.In order to avoid singularity whenyn=yj−1, yj, or yj+1, the average of the neighbor stations’ results is taken as the result of the corresponding term [13]. After the calculation of αi, the effective angle of attack αeff is obtained by subtracting αi from the geometric angle of attack αgeo.
(6)αeffyn=αgeoyn−αiyn.At this point, it should be stated that the method is a numerical iterative solution to Prandtl’s lifting-line theory that states that the circulation is zero wheny=±b/2 [14]. That implies that Cl values of the airfoils at the tip chords ct of the wing are zero and αeff is equal to zero lift angle of attack value αCl=0 of the airfoils at the tips as shown in
(7)αeff±b2=αCl=0±b2.Onceαeff distribution along the span is known, Cl and Cdof the airfoils at the stations are calculated by using Xfoil [15]. An interface developed to run Xfoil in MATLAB is used [16]. It is modified for the MAC OSX operating system. Then, the new circulation distribution Γnew is calculated by using
(8)Γnewy=V∞cyCly2.cy in the above equation is the chord value at the station. The input circulation values are calculated by using a damping factor D as follows:
(9)Γy=Γy+DΓnewy−Γy.The iterative solution requires a strong damping factor. SettingD values as one leads to divergent iteration [10]. The iteration process stops if the Euclidean norm of the difference between the new circulation vector and the input circulation vector is less than a convergence criterion ε:
(10)Γnewy−Γy<ε.By using converged circulation distribution, the liftL [10], the induced drag Di [10], and the root bending moment Mb [17] of the wing are calculated as follows:
(11)L=ρ∞V∞∫−b/2b/2Γnydy,(12)Di=ρ∞V∞∫−b/2b/2Γnysinαinydy,(13)Mb=ρ∞V∞∫0b/2Γnyydy,where ρ∞ is the free stream density. These integrals are also evaluated numerically by using Simpson’s rule. The parasite drag of the wing Dpr is calculated by using
(14)Dpr=12ρ∞V∞2Δy4∑j=1N−1cnj+cnj+1Cdj+Cdj+1.After the calculations of the aerodynamic forces and moments, the lift coefficientCL, the drag coefficient CDw, and the root bending moment coefficient CMb of the wing are obtained as follows:
(15)CL=2Lρ∞V∞2S,(16)CDw=2Dpr+Diρ∞V∞2S,(17)CMb=4Mbρ∞V∞2Sb.As a result, Xfoil and the numerical nonlinear lifting-line method are coupled and the aerodynamic analysis tool is obtained.


## 2.2. XFLR5 LLT and ANSYS-Fluentγ−Reθ SST Model

The XFLR5 LLT and ANSYS-Fluentγ−Reθ SST model are introduced in this section. XFLR5 is an analysis tool for airfoils and wings. There are four different solvers based on the solution of the inviscid potential flow implemented into it. These are the solvers based on the horseshoe vortex lattice method, the ring vortex lattice method, the 3-D panel method, and the nonlinear lifting-line theory. They were used to calculate the drag polar of the different wings [18–20].The solver based on nonlinear lifting-line theory calculates theCL and CDw results by taking the drag polar data of the airfoil form Xfoil that is also implemented into XFLR5. The results are obtained by using the default values defined from the solver except for the iteration number. It is set to 1000.Another computational fluid dynamics (CFD) tool, ANSYS-Fluent, which solves the problems according to the finite volume method, is used to compare theCL and CDw results obtained by the aerodynamic analysis tool.Langtry and Menter developed a new correlation-based transition model to simulate laminar to turbulent transition [21]. The model is essentially based on two transport equations. One equation is for a transition onset criterion and the other is for intermittency. To validate this model, they gave examples of the cases which they studied such as 2-D three-element flap, 2-D airfoil, and a transonic wing. They stated that the transitional CFD simulation was in very good agreement with the experimental results. Aftab et al. compared the following turbulence models for the flow over the NACA4415 airfoil for Re of 120000 with the experimental results: Spalart-Allmaras, K−ω SST, intermittency SST, k−kl−ω, and γ−Reθ SST [22]. According to the results, only the γ−Reθ SST provided reliable results for low and high α values. Lanzafame et al. developed a horizontal axis wind turbine 3-D CFD model to predict wind turbine performance [23]. They compared K−ω SST and γ−Reθ SST with experimental data. According to the result, γ−Reθ SST captured the trend of the aerodynamic coefficients whereas K−ω SST overpredicted Cl values and underpredicted Cd values.According to these studies in the literature, the computations are performed withγ−Reθ SST in this study. The numerical studies are carried out according to the pressure-based method. There are some specific criteria that determine the accuracy of the results obtained with ANSYS-Fluent. The first is to determine the size of the boundary region where flow analysis is intended. For this study, the size of the boundary is formed approximately 30 times larger than the maximum wing thickness and wing chord length as shown in Figure 2 [24, 25]. Since the control volumes are 3-D, the volumetric adaptation of the inflation layers and the prism cells are used. The dimensionless distance to the wall y+ is selected as 1 [26]. Thereby, first-layer height is calculated as 0.015 mm for the selected y+. The growth rate and the number of inflation layers are selected as 1.1 and 61, respectively. These selections create a mesh with 2569785 elements and with a maximum skewness value of 0.95 as shown in Figure 3. The turbulence intensity is selected as 1%. For the convergence criterion, the continuity equation is evaluated as 2·10−4 for 0° α. In order to provide convergence, approximately 1600 iterations are performed. The convergence criterion for high α values increases up to 10 times. Therefore, the results are accepted as accurate in these calculations when fluctuations in CL and CDw lie within the range of ±0.001.Figure 2
The boundaries of the solution domain.Figure 3
Cells around the airfoil atcr.


## 2.3. Comparison of the Aerodynamic Analysis Tool with ANSYS-Fluent

The comparison of the aerodynamic analysis tool with XFLR5 6.47 LLT and the ANSYS-Fluent 16.2γ−Reθ SST model is made in this part of the paper. The flow solutions of the aerodynamic analysis tool are obtained on a 2.4 GHz Intel Core i7, 8 GB MacBook Pro. The computation of ANSYS-Fluent is performed on 2 processors of 2.4 GHz Intel Xeon E5, 32 GB HP Z820.The interface is modified so thatCl and Cd values are obtained for the airfoils at each station by parallel processing. When Xfoil fails to converge, the initial panel number, which is 200, is altered between 194 and 206. NACA 4412 airfoil with closed trading edge is used as the candidate airfoil. Equations (18) and (19) are used for the half thickness zt and the camber generation zc, respectively. In these equations, t is the maximum thickness in hundredths of the chord, p is the maximum camber location in tenths of the chord, and m is the maximum camber in hundredths of the chord. x coordinates are generated by the cosine distribution.Thez coordinates of the upper surface are obtained by summing equations (18) and (19), whereas equation (18) is subtracted from equation (19) for the lower surface as is done in the NACA airfoil generation subroutine of Xfoil [27].
(18)ztx=t0.20.2969x0.5−0.126x−0.3516x2+0.2843x3−1036x4,(19)zcx=mp22px−x2,if0≤x<p,zcx=m1−p21−2p+2px−x2,ifp≤x≤1.cr, ct, and b are 0.45m, 0.45m, and 4m, respectively. ρ∞ is 1.225kg/m3. The dynamic viscosity μ∞ is 1.7974·10−5Ns/m2 and V∞ is 22 m/s. In Xfoil, the default critical n value of the en transition prediction method is nine. This value is changed as 2.62, because it corresponds to the selected turbulent intensity in the ANSYS-Fluent analysis.In Figure4, the CL results for the different panel station numbers and the results of ANSYS-Fluent and XFLR5 are shown. According to the results in Figure 4, ANSYS-Fluent calculates the least CL at the same α. CL values of N=21, N=31, and N=41 cases overlap. The results show slight difference with the results of XFLR5 LLT and capture the trend line of ANSYS-Fluent up to the stall α. CDw values of ANSYS-Fluent are almost identical to N=11 case values when α is between 0° and 10° as depicted in Figure 5. But as α gets close to 15°, the ANSYS-Fluent, N=31 case, and N=41 case yield similar CDw. After 15°, ANSYS-Fluent results start to diverge from the results of the aerodynamic tool. In order to investigate the situation in detail, the flow field around the wing for 10° and 18° α values is presented in Figures 6–9.Figure 4
CL vs. α curves of the aerodynamic design tool, ANSYS-Fluent, and XFLR5.Figure 5
CDw vs. α curves of the aerodynamic design tool, ANSYS-Fluent, and XFLR5.Figure 6
Cp contours of the aerodynamic analysis tool and ANSYS-Fluent for 10°.Figure 7
x-wall shear stress contours of the aerodynamic analysis tool and ANSYS-Fluent for 10°.Figure 8
Cp contours of the aerodynamic analysis tool and ANSYS-Fluent for 18°.Figure 9
x-wall shear stress contours of the aerodynamic analysis tool and ANSYS-Fluent for 18°.In the aerodynamic analysis tool, the transition region can be determined by finding the region at which a sudden increase in thex-wall shear stress is observed. The turbulent region is the region at which the x-wall shear stress values get closer to zero. Contrary to this, the transition region and the turbulent separation region are identified when x-wall shear stress is less than zero in the results of ANSYS-Fluent [28].According to the results for 10°, apart from the slight difference around the leading edge, pressure contours are almost identical for both solvers as shown in Figure6. The transition prediction by the aerodynamic solver is understood from the color change from green to yellow in Figure 7. ANSYS-Fluent captures the transition location almost at the same position as the aerodynamic analysis tool does that is shown by a color change from green to blue. Both solvers capture the separation around the trailing edge. When the results for 18° are discussed, one can easily see the difference between the pressure contours especially on the aft of the wing as depicted in Figure 8. As is shown in Figure 9, both solvers predict the transition region almost at the same location. ANSYS-Fluent predicts the turbulent separation ahead. But it is very close to the turbulent separation location of the aerodynamic analysis tool.SinceN=31 and N=41 cases capture the CDw values of ANSYS-Fluent at high α values, they are selected as the candidate station numbers.In order to select the optimumN, D, and ε values, the total solution time ttot to obtain drag polar results between the α values of -6° and 20° is compared for the different N, D, and ε values in Table 1. The maximum deviations in CL∆CL and CDw∆CDw are obtained with respect to the results when the N, D, and ε values are equal to 41, 0.05, and 0.01, respectively. According to the results, choosing the N value as 31 saves the computation time approximately by 38% at the same D and ε values when it is compared with the result of N=41.Table 1
Effect ofN, D, and ε to the results and ttot.
N
D
ε
∆CL(%)
∆CDw(%)
ttot(s)
41
0.05
0.01
0
0
3007.7
31
0.05
0.01
0.4
1.1
2532.9
31
0.1
0.01
0.4
1.1
1096.2
31
0.2
0.01
0.4
1.1
519.3
31
0.2
0.05
5.9
1.1
418.5
31
0.2
0.1
11.2
1.0
414.3
31
0.2
0.2
15.2
0.9
347.6IncreasingD decreases the computation time by 86.9% without boosting ∆CL and ∆CDw. When D is selected as 0.4 and above, the convergence problem is observed at the high α values. Contrary to this, increasing ε yields drastic increase in ∆CL especially at the negative α and the high α values. Apart from this region, ∆CL is around 1% at the low α and moderate α values. In light of this information, N, D, and ε values are selected as 31, 0.2, and 0.01, respectively. ttot of ANSYS- Fluent results is approximately 120 hours.


## 2.4. Optimization Solver

In this part, the advantages of the selected optimization solver and its use are discussed. Vanderplaats compared different optimization methods by using a structural optimization problem [29]. In the problem, the objective function was the minimization of the structural volume of a cantilevered beam that was fixed at the right end and had a vertical load applied at the free left end. The cantilevered beam had five segments and thicknesses, and heights of the segments were the design variables. In the study, the allowable bending stress at the right end of each segment was defined as the constraints. In addition to this, the allowable left end deflection was also defined as the constraint. The ratio of the height to the thickness of each segment was limited by using constraints. In summary, the identified optimization problem had 10 design variables and 11 constraints. Genetic search, sequential linear programming, method of feasible directions, generalized reduced gradient method, modified feasible directions method, and sequential quadratic programming (SQP) were compared in the study. The first method is an evolutionary optimization method whereas the others are the gradient-based optimization methods. The optimum results obtained by the methods were very close to each other. The genetic search method obtained the optimum result after 106 function evaluations. Among the gradient-based optimization methods, SQP had the least number of iterations and the least number of gradient evaluations. The total number of function call is 125 for the SQP analysis. As discussed above, the main advantage of SQP is that it does not attempt to satisfy constraints at each iteration that results in less number of function calls [30]. Therefore, the SQP algorithm is used as given in the MATLAB optimization toolbox for the solution of the optimization problems that are discussed below.The procedure for the optimization process is represented in Figure10.Figure 10
The flowchart of the optimization process.The design variable array is used to calculate the objective function and the constraints for the current step. For this case, the aerodynamic design tool is called in parallel computations and the threads are used to calculateCl and Cd values for the airfoils at each station. After that, the design variable array is perturbed by an increment array that is 1% of each design variable array element. Parallel computing is applied for the calculation of the objective function and the constraints by using perturbed design variable array elements. The aerodynamic design tool runs in serial for this case. The forward difference method is used for the calculation of the gradient of the objective function and Jacobian matrix of the constraints.The process stops if the change of the design variable array (TolX) or the change in the objective function (TolFun) is less than 10−3. The violation of the constraints (TolCon) is limited to 10−2.


## 2.5. The Characteristics of the Baseline UAV

Before going into the details of the optimization problems, the geometric and performance characteristics of the baseline UAV and basic assumptions related to it are revealed in this part of the work. The baseline UAV is a propeller-driven aircraft andPA of the engine is 2000 W. The drag coefficient of the other components (CDoth) of the UAV except the wing is assumed as a varying parameter as follows:
(20)CDoth=0.02SbS,where Sb is the planform area of the wing of the baseline UAV. As a result, the total drag coefficient CD of the baseline UAV is the summation of CDoth and CDw. The total weight of the UAV except for the wing (Woth) is 250 N. The wing weight Ww is calculated as follows [31]:
(21)Ww=Sc¯tcmaxρmatKρARnult0.6λ0.04g,where c¯ denotes the wing mean aerodynamic chord, t/cmax the maximum thickness of the airfoil, AR the aspect ratio, and λ the taper ratio. ρmat is the density of the construction material that is graphite/epoxy. Kρ is the wing density factor. They are selected as 1575 kg/m3 and 0.0016, respectively, by calculating the average values in the appropriate reference tables for homebuilt aircraft [31]. The relation between ultimate load factor (nult) and maximum positive load factor (nmax) is defined as follows:
(22)nult=1.5nmax.Unlike a fighter aircraft, the produciblenmax by the engine power is less than the tolerable nmax by the aircraft structure for a general aviation aircraft [32]. For this reason, the equation of the producible nmax by the engine power for the propeller-driven aircraft that is in equation (23) is used in the Ww calculation.
(23)nmax=CL3CD2PA2ρ∞S21/31Woth+Ww.As seen in the above equations,Ww and nmax are implicit. The secant method is applied to find Ww.The geometric parameters of the wing of the baseline UAV are shown in Table2.Table 2
Geometric properties of the baseline UAV’s wing.
Airfoil
NACA4412
cr
0.45 m
ct
0.45 m
b
4
Sb
1.8 m2The wing of the baseline UAV is a rectangular and unswept wing. The airfoil profile is constant along the span. Table3 depicts the performance parameters of the baseline UAV.Table 3
Performance parameters of the baseline UAV.
maxCL3/2/CD
14.08
α at maxCL3/2/CD
8°
V at maxCL3/2/CD
15.52 m/s
Mb at maxCL3/2/CD
127.46 Nm
Ww
24.06 N
Vmax
39.45 m/s
α at Vmax
-2.32°
Vstall
13.25
α at Vstall
18°It has a maximumCL3/2/CD ratio of 14.08 when α is 8° and speed is 15.52 m/s. Mb is 127.46 Nm for this flight condition. Ww is 24.06 N. Vmax is 39.45 m/s when α is -2.32° and Vstall is 13.25 m/s when α is 18°.


## 3. Results

First of all, the details of the optimization problems are discussed in this section. Afterward, their results are analyzed.


### 3.1. The Optimization Problems

In the optimization problems,t, p, and m are defined as the design variables for airfoil generation. cr, ct, and b are the design variables for the generation of the wing planform. α and V for different level flight conditions are defined as the design variables. Three different optimization problems are solved. The first optimization problem is defined between equations (24) and (33) as follows:
(24)min−CL3/2/CDatα1andV1,(25)subjectto−6°≤α1≤25°,(26)5m/s≤V1≤60m/s,(27)Woth+Ww=Latα1andV1,(28)0.08≤t≤0.16,(29)0.15≤p≤0.6,(30)0.02≤m≤0.08,(31)0.2m≤cr≤0.6m,(32)0.2m≤ct≤0.6m,(33)4m≤b≤8m.The maximization ofCL3/2/CD is defined in equation (24). The upper limit of the α at the maximum CL3/2/CD condition α1 is 25° whereas the lower limit is -8°. The speed at the maximum CL3/2/CD condition V1 is also defined as a design variable so that the aerodynamic design tool calculates the aerodynamic forces for the appropriate Re. Equation (27) is defined to satisfy the level flight condition at the maximum CL3/2/CD condition. The equations between equations (28) and (33) define the upper and the lower limit of the airfoil shape and the wing planform geometry design variables. The second optimization problem is obtained by adding three functional inequality constraints to the first optimization problem as follows:
(34)Mb≤127.46NmatmaxCL3/2/CD,(35)Ww≤24.06N,(36)S≥1.8m.Equations (34) and (35) constrain the root bending moment at the maximum CL3/2/CD condition and the wing weight in the second optimization problem. In addition to these, equation (36) sets the minimum planform area as 1.8 m2.The equations between Equations (37) and (43) are added into the second optimization problem, and the third optimization problem is obtained.
(37)Woth+Ww=Latα2andV2,(38)39.45m/s≤V2≤60m/s,(39)PR≤2000Watα2andV2,(40)−6°≤α2≤25°,(41)Woth+Ww=Latα3andV3,(42)5m/s≤V3≤13.25m/s,(43)−6°≤α3≤25°.In the third optimization problem, the level flight condition for theV2 and α2 values are defined in equation (37). This corresponds to the level flight of the baseline UAV with Vmax. V2 has a lower limit of 39.45 m/s so that the optimized UAV shall have an equal or greater Vmax than the baseline UAV. In addition to this, PR at V2 and α2 is constrained with PA of the UAV in equation (39). Equation (41) describes the level flight of the baseline UAV with Vstall. In the same manner, the upper limit of V3 is set as 13.25 m/s. The baseline UAV is the initial starting point for all optimization problems. The objective function, the constraints, the design variables, and their upper and lower limits are scaled with their initial values. In summary, the first optimization problem has eight design variables and one equality constraint. The second optimization problem has eight design variables and one equality and three inequality constraints. Finally, the third optimization problem has 12 design variables and three equality and four inequality constraints.


### 3.2. Comparison of the Optimization Problem Solutions

According to the results, the highest maximumCL3/2/CD is obtained for the first optimization problem when α1 is 5.56° and V1 is 11.97 m/s as depicted in Table 4. This corresponds to the neck in the CL−α curve in Figure 11. Contrary to this, Mb at maximum CL3/2/CD nearly doubles and the increment in Ww is about 55%. According to the results in Figure 12, Vmax of the first optimization problem is 34.2 m/s that is nearly 5 m/s less than the Vmax of the baseline UAV. Vstall of the first optimization problem is 11.33 m/s that is almost the same as its V1. The main reason for this situation is CL of the maximum CL3/2/CD is very close to the maximum lift coefficient CLmax. When the airfoil shape and the wing planform details in Table 5 and Figures 13 and 14 are investigated, it is seen that the airfoil of the first optimization problem reaches the lower limit of the thickness and the upper limit of the camber. ct and b are set to the lower limit and upper limit, respectively. cr is decreased from 0.45 m to 0.386 m that alters the taper ratio λ from 1 to 0.52.Table 4
Performance parameters of baseline UAV and the optimum UAVs.
Baseline UAV
Case 1
Case 2
Case 3
maxCL3/2/CD
14.08
30.79
18.95
18.21
α1atmaxCL3/2/CD (°)
8
5.56
5.69
6.44
V1atmaxCL3/2/CD (m/s)
15.52
11.97
13.41
14.41
MbatmaxCL3/2/CD (Nm)
127.46
249.4
127.3
127.4
WwN
24.06
37.64
19.41
21.8
V2m/s
39.45
39.45
α2atV2 (°)
-2.32
-5.3
V3m/s
13.25
12.58
α3atV3 (°)
18
17.86Figure 11
CL vs. α curves for the baseline UAV and the optimum UAV wings.Figure 12
PR vs. V curves for the baseline UAV and the optimum UAVs.Table 5
Geometric properties of the wings of the baseline UAV and the optimum UAV wings.
Baseline UAV
Case 1
Case 2
Case 3
t
0.12
0.08
0.08
0.092
p
0.4
0.524
0.566
0.494
m
0.04
0.08
0.08
0.063
crm
0.45
0.386
0.586
0.596
ctm
0.45
0.2
0.2
0.2
bm
4
8
4.578
4.524
Sm2
1.8
2.34
1.8
1.8Figure 13
The airfoils of the baseline UAV and the optimum UAV wings.Figure 14
The wing planform views of the baseline UAV and the optimum UAV wings.The airfoil shape of the second optimization problem is almost identical to the airfoil shape of the first optimization problem. Contrary to this,cr and b have significant change due to the addition of the functional constraints related with Mb at the maximum CL3/2/CD and S. Mb reaches its upper limit whereas S is at its lower limit. But Ww is distant from its bound. As a result, the increment at the maximum CL3/2/CD is decreased from 220% to 34.6%. Vmax of the second optimization problem is 37.2 m/s according to the PR−V curve in Figure 12 whereas the Vstall is 12.3 m/s. Since S is decreased from 2.34 m2 to 1.8 m2, the increase in Vstall is observed.The wing planform shape of the third optimization problem is almost the same as the wing planform shape of the second optimization problem. As clearly seen in the results, the functional constraints related withMb at the maximum CL3/2/CD and S are the key reason for the similar outcomes from the different optimization problems. They reach their limits whenever they are defined as the functional constraints. Vmax of the third optimization problem is identical to V2 of the same problem according to Table 4 and Figure 12. This implies that the third optimization problem barely has the same Vmax of the baseline UAV. The airfoil of the third optimization problem is thicker and less cambered when it is compared with the airfoil of the second optimization problem. This change is caused by the addition of the level flight performance analysis at the V2 speed. Vstall of the third optimization problem is approximately 0.7 m/s less than Vstall of the baseline UAV. According to this result, it can be said that defining a functional constraint related with Vstall in CL3/2/CD maximization problems has no significant contribution, because all optimized wings satisfy the level flight condition at a speed less than Vstall of the baseline UAV.Table6 depicts the number of iterations I, the number of the function count fc, and the elapsed solution time te of the optimization problems. According to the result, the third optimization problem has the longest solution time because the aerodynamic design tool is called three times at each function count due to the analysis of three different level flight phases. te of the third problem is 6.64 hours. If ANSYS-Fluent was used as a function evaluator in the third optimization problem, te would be around 500 hours. This comparison shows the cost and time effectiveness of the aerodynamic analysis tool.Table 6
I, fc, and te of Case 1, Case 2, and Case 3.
Case 1
Case 2
Case 3
I
13
18
12
fc
32
26
38
tes
8520
10505
23896The evaluation of the design variables and the objective function with the iteration number for the optimization problems are presented in Figures15, 16, and 17.Figure 15
The variation of the objective function and the design variables in Case 1.Figure 16
The variation of the objective function and the design variables in Case 2.Figure 17
The variation of the objective function and the design variables in Case 3.


## 3.1. The Optimization Problems

In the optimization problems,t, p, and m are defined as the design variables for airfoil generation. cr, ct, and b are the design variables for the generation of the wing planform. α and V for different level flight conditions are defined as the design variables. Three different optimization problems are solved. The first optimization problem is defined between equations (24) and (33) as follows:
(24)min−CL3/2/CDatα1andV1,(25)subjectto−6°≤α1≤25°,(26)5m/s≤V1≤60m/s,(27)Woth+Ww=Latα1andV1,(28)0.08≤t≤0.16,(29)0.15≤p≤0.6,(30)0.02≤m≤0.08,(31)0.2m≤cr≤0.6m,(32)0.2m≤ct≤0.6m,(33)4m≤b≤8m.The maximization ofCL3/2/CD is defined in equation (24). The upper limit of the α at the maximum CL3/2/CD condition α1 is 25° whereas the lower limit is -8°. The speed at the maximum CL3/2/CD condition V1 is also defined as a design variable so that the aerodynamic design tool calculates the aerodynamic forces for the appropriate Re. Equation (27) is defined to satisfy the level flight condition at the maximum CL3/2/CD condition. The equations between equations (28) and (33) define the upper and the lower limit of the airfoil shape and the wing planform geometry design variables. The second optimization problem is obtained by adding three functional inequality constraints to the first optimization problem as follows:
(34)Mb≤127.46NmatmaxCL3/2/CD,(35)Ww≤24.06N,(36)S≥1.8m.Equations (34) and (35) constrain the root bending moment at the maximum CL3/2/CD condition and the wing weight in the second optimization problem. In addition to these, equation (36) sets the minimum planform area as 1.8 m2.The equations between Equations (37) and (43) are added into the second optimization problem, and the third optimization problem is obtained.
(37)Woth+Ww=Latα2andV2,(38)39.45m/s≤V2≤60m/s,(39)PR≤2000Watα2andV2,(40)−6°≤α2≤25°,(41)Woth+Ww=Latα3andV3,(42)5m/s≤V3≤13.25m/s,(43)−6°≤α3≤25°.In the third optimization problem, the level flight condition for theV2 and α2 values are defined in equation (37). This corresponds to the level flight of the baseline UAV with Vmax. V2 has a lower limit of 39.45 m/s so that the optimized UAV shall have an equal or greater Vmax than the baseline UAV. In addition to this, PR at V2 and α2 is constrained with PA of the UAV in equation (39). Equation (41) describes the level flight of the baseline UAV with Vstall. In the same manner, the upper limit of V3 is set as 13.25 m/s. The baseline UAV is the initial starting point for all optimization problems. The objective function, the constraints, the design variables, and their upper and lower limits are scaled with their initial values. In summary, the first optimization problem has eight design variables and one equality constraint. The second optimization problem has eight design variables and one equality and three inequality constraints. Finally, the third optimization problem has 12 design variables and three equality and four inequality constraints.


## 3.2. Comparison of the Optimization Problem Solutions

According to the results, the highest maximumCL3/2/CD is obtained for the first optimization problem when α1 is 5.56° and V1 is 11.97 m/s as depicted in Table 4. This corresponds to the neck in the CL−α curve in Figure 11. Contrary to this, Mb at maximum CL3/2/CD nearly doubles and the increment in Ww is about 55%. According to the results in Figure 12, Vmax of the first optimization problem is 34.2 m/s that is nearly 5 m/s less than the Vmax of the baseline UAV. Vstall of the first optimization problem is 11.33 m/s that is almost the same as its V1. The main reason for this situation is CL of the maximum CL3/2/CD is very close to the maximum lift coefficient CLmax. When the airfoil shape and the wing planform details in Table 5 and Figures 13 and 14 are investigated, it is seen that the airfoil of the first optimization problem reaches the lower limit of the thickness and the upper limit of the camber. ct and b are set to the lower limit and upper limit, respectively. cr is decreased from 0.45 m to 0.386 m that alters the taper ratio λ from 1 to 0.52.Table 4
Performance parameters of baseline UAV and the optimum UAVs.
Baseline UAV
Case 1
Case 2
Case 3
maxCL3/2/CD
14.08
30.79
18.95
18.21
α1atmaxCL3/2/CD (°)
8
5.56
5.69
6.44
V1atmaxCL3/2/CD (m/s)
15.52
11.97
13.41
14.41
MbatmaxCL3/2/CD (Nm)
127.46
249.4
127.3
127.4
WwN
24.06
37.64
19.41
21.8
V2m/s
39.45
39.45
α2atV2 (°)
-2.32
-5.3
V3m/s
13.25
12.58
α3atV3 (°)
18
17.86Figure 11
CL vs. α curves for the baseline UAV and the optimum UAV wings.Figure 12
PR vs. V curves for the baseline UAV and the optimum UAVs.Table 5
Geometric properties of the wings of the baseline UAV and the optimum UAV wings.
Baseline UAV
Case 1
Case 2
Case 3
t
0.12
0.08
0.08
0.092
p
0.4
0.524
0.566
0.494
m
0.04
0.08
0.08
0.063
crm
0.45
0.386
0.586
0.596
ctm
0.45
0.2
0.2
0.2
bm
4
8
4.578
4.524
Sm2
1.8
2.34
1.8
1.8Figure 13
The airfoils of the baseline UAV and the optimum UAV wings.Figure 14
The wing planform views of the baseline UAV and the optimum UAV wings.The airfoil shape of the second optimization problem is almost identical to the airfoil shape of the first optimization problem. Contrary to this,cr and b have significant change due to the addition of the functional constraints related with Mb at the maximum CL3/2/CD and S. Mb reaches its upper limit whereas S is at its lower limit. But Ww is distant from its bound. As a result, the increment at the maximum CL3/2/CD is decreased from 220% to 34.6%. Vmax of the second optimization problem is 37.2 m/s according to the PR−V curve in Figure 12 whereas the Vstall is 12.3 m/s. Since S is decreased from 2.34 m2 to 1.8 m2, the increase in Vstall is observed.The wing planform shape of the third optimization problem is almost the same as the wing planform shape of the second optimization problem. As clearly seen in the results, the functional constraints related withMb at the maximum CL3/2/CD and S are the key reason for the similar outcomes from the different optimization problems. They reach their limits whenever they are defined as the functional constraints. Vmax of the third optimization problem is identical to V2 of the same problem according to Table 4 and Figure 12. This implies that the third optimization problem barely has the same Vmax of the baseline UAV. The airfoil of the third optimization problem is thicker and less cambered when it is compared with the airfoil of the second optimization problem. This change is caused by the addition of the level flight performance analysis at the V2 speed. Vstall of the third optimization problem is approximately 0.7 m/s less than Vstall of the baseline UAV. According to this result, it can be said that defining a functional constraint related with Vstall in CL3/2/CD maximization problems has no significant contribution, because all optimized wings satisfy the level flight condition at a speed less than Vstall of the baseline UAV.Table6 depicts the number of iterations I, the number of the function count fc, and the elapsed solution time te of the optimization problems. According to the result, the third optimization problem has the longest solution time because the aerodynamic design tool is called three times at each function count due to the analysis of three different level flight phases. te of the third problem is 6.64 hours. If ANSYS-Fluent was used as a function evaluator in the third optimization problem, te would be around 500 hours. This comparison shows the cost and time effectiveness of the aerodynamic analysis tool.Table 6
I, fc, and te of Case 1, Case 2, and Case 3.
Case 1
Case 2
Case 3
I
13
18
12
fc
32
26
38
tes
8520
10505
23896The evaluation of the design variables and the objective function with the iteration number for the optimization problems are presented in Figures15, 16, and 17.Figure 15
The variation of the objective function and the design variables in Case 1.Figure 16
The variation of the objective function and the design variables in Case 2.Figure 17
The variation of the objective function and the design variables in Case 3.


## 4. Conclusion

This paper demonstrates a rapid aerodynamic design tool for the wing optimization by using the maximum thickness, the maximum camber, the maximum camber location, the root chord, the tip chord, the span, the free stream velocity, and the angle of attack as the design variables. The function evaluator of the aerodynamic design tool consists of the Xfoil and nonlinear numerical lifting-line theory. Sequential quadratic programming is the optimization solver in the design tool. The derivatives of the objective function and the constraints with respect to the design variables are obtained by the use of parallel programming. Three different optimization problems are solved. According to the results,(1)
Defining the functional constraint related to the stall speed performance of the UAV inCL3/2/CD maximization problems has no effect because all optimized UAVs have less stall speed than the baseline UAV. Moreover, the stall speed analysis increases the solution time
(2)
When the bending moment at the root, the wing weight, and the wing planform are limited at the same time, the bending moment reaches its upper limit whereas the wing area is at its lower limit. Contrary to this, the wing weight is distant from its upper limit. That means the bending moment at the root and the planform area constraints form the planform shape. But their effect on the airfoil shape change is insignificant
(3)
Limiting the maximum speed of the UAV yields thicker and less cambered airfoilIn summary, defining the constraints related to the bending moment at the root, the wing area, and the maximum speed of the UAV ensures the evaluation of the outcomes of the optimization problems to the more down-to-earth design results.



---

*Source: 3050824-2019-10-16.xml*