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

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 of $C_L^{3/2}/C_D$ is the objective function, while level flight condition at maximum $C_L^{3/2}/C_D$ 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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 $e^n$ 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 $\left(\alpha \right)$ 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 $\left(C_l\right)$ and the profile drag coefficient $\left(C_d\right)$ 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 $\left(\mathrm{Re}\right)$ 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 $\mathrm{Re}$ was 6000. Xfoil captured $C_l$ and $C_d$ values up to stall $\alpha$ 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 $\alpha$ 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 $\left(\gamma -{\mathrm{Re}}_{\theta }\right)$ shear stress transport (SST) model implemented into ANSYS-Fluent. After that, three different optimization problems in which $C_L^{3/2}/C_D$ maximization is the objective function are constructed. The optimization problems are defined by adding new functional constraints that are the wing root bending moment $\left(M_b\right)$, the wing weight $\left(W_w\right)$, the wing planform area $\left(S\right)$, the level flight condition and power available $\left(P_A\right)$ at the maximum speed $\left(V_{\max }\right)$, and the level flight condition at the stall speed $\left(V_{\text{stall}}\right)$ of the baseline unmanned air vehicle $\left(\text{UAV}\right)$. 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 $\text{UAV}$ at both $V_{\max }$ and $V_{\text{stall}}$ 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 span $\left(b\right)$ 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 $\Delta y$. The method needs an initial circulation distribution $\left(\Gamma \left(y\right)\right)$ along the span. The elliptic circulation distribution in equation (1) is acceptable for this requirement. $$\begin{array}{}\text{(1)}& \Gamma \left(y\right)={\Gamma }_0\sqrt{1-{\left(\frac{2y}{b}\right)}^2}.\end{array}$$

[figure omitted; refer to PDF]

In the above equation, ${\Gamma }_0$ is the circulation at the root chord $\left(c_r\right)$. By using initial $\Gamma \left(y\right)$, the induced angle of attack $\left({\alpha }_i\right)$ for the $n^{\text{th}}$ panel is obtained as follows: $$\begin{array}{}\text{(2)}& {\alpha }_i\left(y_n\right)=\frac{1}{4\pi V_{\infty }}{\int }_{-b/2}^{b/2}\frac{\left(d\Gamma /dy\right)dy}{y_n-y}.\end{array}$$

This integral is evaluated numerically by using Simpson’s rule [13]: $$\begin{array}{}\text{(3)}& {\alpha }_i\left(y_n\right)=\frac{\Delta y}{12\pi V_{\infty }}\sum _{j=2,4,6}^{N-1}\left(\frac{{\left(d\Gamma /dy\right)}_{j-1}}{y_n-y_{j-1}}\right.+\frac{{\left(d\Gamma /dy\right)}_j}{y_n-y_j}+\left.\frac{{\left(d\Gamma /dy\right)}_{j+1}}{y_n-y_{j+1}}\right),\end{array}$$where $V_{\infty }$ is the freestream velocity. For the calculation of ${\alpha }_i$, $d\Gamma /dy$ is needed. Although $d\Gamma /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. $$\begin{array}{}\text{(4)}& \frac{d\Gamma }{dy}=-\frac{4{\Gamma }_0}{b^2}\frac{y}{{\left(1-4y^2/b^2\right)}^{1/2}},\text{(5)}& {\left(\frac{d\Gamma }{dy}\right)}_j=\frac{{\Gamma }_{j+1}-{\Gamma }_j}{y_{j+1}-y_j}.\end{array}$$

In order to avoid singularity when $y_n=y_{j-1}$, $y_j$, or $y_{j+1}$, the average of the neighbor stations’ results is taken as the result of the corresponding term [13]. After the calculation of ${\alpha }_i$, the effective angle of attack $\left({\alpha }_{\text{eff}}\right)$ is obtained by subtracting ${\alpha }_i$ from the geometric angle of attack $\left({\alpha }_{\text{geo}}\right)$. $$\begin{array}{}\text{(6)}& {\alpha }_{\text{eff}}\left(y_n\right)={\alpha }_{\text{geo}}\left(y_n\right)-{\alpha }_i\left(y_n\right).\end{array}$$

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 when $y=\pm b/2$ [14]. That implies that $C_l$ values of the airfoils at the tip chords $\left(c_t\right)$ of the wing are zero and ${\alpha }_{\text{eff}}$ is equal to zero lift angle of attack value $\left({\alpha }_{C_l=0}\right)$ of the airfoils at the tips as shown in $$\begin{array}{}\text{(7)}& {\alpha }_{\text{eff}}\left(\pm \frac{b}{2}\right)={\alpha }_{C_l=0}\left(\pm \frac{b}{2}\right).\end{array}$$

Once ${\alpha }_{\text{eff}}$ distribution along the span is known, $C_l$ and $C_d$of 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 $\left({\Gamma }_{\text{new}}\right)$ is calculated by using $$\begin{array}{}\text{(8)}& {\Gamma }_{\text{new}}\left(y\right)=\frac{V_{\infty }c\left(y\right)C_l\left(y\right)}{2}.\end{array}$$

$c\left(y\right)$ in the above equation is the chord value at the station. The input circulation values are calculated by using a damping factor $\left(D\right)$ as follows: $$\begin{array}{}\text{(9)}& \Gamma \left(y\right)=\Gamma \left(y\right)+D\left({\Gamma }_{\text{new}}\left(y\right)-\Gamma \left(y\right)\right).\end{array}$$

The iterative solution requires a strong damping factor. Setting $D$ 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 $\left(\epsilon \right)$: $$\begin{array}{}\text{(10)}& ‖{\Gamma }_{\text{new}}\left(y\right)-\Gamma \left(y\right)‖<\epsilon .\end{array}$$

By using converged circulation distribution, the lift $\left(L\right)$ [10], the induced drag $\left(D_i\right)$ [10], and the root bending moment $\left(M_b\right)$ [17] of the wing are calculated as follows: $$\begin{array}{}\text{(11)}& L={\rho }_{\infty }{V}_{\infty }{\int }_{-b/2}^{b/2}{\Gamma }_n\left(y\right)dy,\text{(12)}& D_i={\rho }_{\infty }{V}_{\infty }{\int }_{-b/2}^{b/2}{\Gamma }_n\left(y\right)\sin \left({{\alpha }_i}_n\left(y\right)\right)dy,\text{(13)}& M_b={\rho }_{\infty }{V}_{\infty }{\int }_0^{b/2}{\Gamma }_n\left(y\right)ydy,\end{array}$$where ${\rho }_{\infty }$ is the free stream density. These integrals are also evaluated numerically by using Simpson’s rule. The parasite drag of the wing $\left(D_{\text{pr}}\right)$ is calculated by using $$\begin{array}{}\text{(14)}& D_{\text{pr}}=\frac{1}{2}{\rho }_{\infty }{V}_{\infty }^2\frac{\Delta y}{4}\sum _{j=1}^{N-1}\left[\left(c_n\left(j\right)+c_n\left(j+1\right)\right)\left(C_d\left(j\right)+C_d\left(j+1\right)\right)\right].\end{array}$$

After the calculations of the aerodynamic forces and moments, the lift coefficient $\left(C_L\right)$, the drag coefficient $\left(C_{\text{Dw}}\right)$, and the root bending moment coefficient $\left(C_{\text{Mb}}\right)$ of the wing are obtained as follows: $$\begin{array}{}\text{(15)}& C_L=\frac{2L}{{\rho }_{\infty }{V}_{\infty }^{2}S},\text{(16)}& C_{\text{Dw}}=\frac{2\left(D_{\text{pr}}+D_i\right)}{{\rho }_{\infty }{V}_{\infty }^{2}S},\text{(17)}& C_{\text{Mb}}=\frac{4M_b}{{\rho }_{\infty }{V}_{\infty }^{2}Sb}.\end{array}$$

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 $\gamma -{\mathrm{Re}}_{\theta }$ SST Model

The XFLR5 LLT and ANSYS-Fluent $\gamma -{\mathrm{Re}}_{\theta }$ 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 the $C_L$ and $C_{\text{Dw}}$ 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 the $C_L$ and $C_{\text{Dw}}$ 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 $\mathrm{Re}$ of 120000 with the experimental results: Spalart-Allmaras, $K-\omega$ SST, intermittency SST, $k-kl-\omega$, and $\gamma -{\mathrm{Re}}_{\theta }$ SST [22]. According to the results, only the $\gamma -{\mathrm{Re}}_{\theta }$ SST provided reliable results for low and high $\alpha$ values. Lanzafame et al. developed a horizontal axis wind turbine 3-D CFD model to predict wind turbine performance [23]. They compared $K-\omega$ SST and $\gamma -{\mathrm{Re}}_{\theta }$ SST with experimental data. According to the result, $\gamma -{\mathrm{Re}}_{\theta }$ SST captured the trend of the aerodynamic coefficients whereas $K-\omega$ SST overpredicted $C_l$ values and underpredicted $C_d$ values.

According to these studies in the literature, the computations are performed with $\gamma -{\mathrm{Re}}_{\theta }$ 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 $\left(y^{+}\right)$ 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° $\alpha$. In order to provide convergence, approximately 1600 iterations are performed. The convergence criterion for high $\alpha$ values increases up to 10 times. Therefore, the results are accepted as accurate in these calculations when fluctuations in $C_L$ and $C_{\text{Dw}}$ lie within the range of $\pm 0.001$.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

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 $\gamma -{\mathrm{Re}}_{\theta }$ 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 that $C_l$ and $C_d$ 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 $\left(z_t\right)$ and the camber generation $\left(z_c\right)$, 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.

The $z$ 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]. $$\begin{array}{c}\text{(18)}& z_t\left(x\right)=\frac{t}{0.2}\left(0.2969x^{0.5}-0.126x-0.3516x^2+0.2843x^3-1036x^4\right),\text{(19)}& z_c\left(x\right)=\frac{m}{p^2}\left(2px-x^2\right),\hspace{1em}\text{if}0\le x<p,z_c\left(x\right)=\frac{m}{\left(1-p^2\right)}\left(1-2p+2px-x^2\right),\hspace{1em}\mathit{\text{if}}p\le x\le 1.\end{array}$$

$c_r$, $c_t$, and $b$ are $0.45\text{m}$, $0.45\text{m}$, and $4\text{m}$, respectively. ${\rho }_{\infty }$ is $1.225\text{kg}/{\text{m}}^3$. The dynamic viscosity $\left({\mu }_{\infty }\right)$ is $1.7974·{10}^{-5}\text{Ns}/{\text{m}}^2$ and $V_{\infty }$ is 22 $\text{m}/\text{s}$. In Xfoil, the default critical $n$ value of the $e^n$ 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 Figure 4, the $C_L$ 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 $C_L$ at the same $\alpha$. $C_L$ 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 $\alpha$. $C_{\text{Dw}}$ values of ANSYS-Fluent are almost identical to $N=11$ case values when $\alpha$ is between 0° and 10° as depicted in Figure 5. But as $\alpha$ gets close to 15°, the ANSYS-Fluent, $N=31$ case, and $N=41$ case yield similar $C_{\text{Dw}}$. 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° $\alpha$ values is presented in Figures 6–9.

[figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF]

In the aerodynamic analysis tool, the transition region can be determined by finding the region at which a sudden increase in the $x$-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 Figure 6. 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.

Since $N=31$ and $N=41$ cases capture the $C_{\text{Dw}}$ values of ANSYS-Fluent at high $\alpha$ values, they are selected as the candidate station numbers.

In order to select the optimum $N$, $D$, and $\epsilon$ values, the total solution time $\left(t_{\text{tot}}\right)$ to obtain drag polar results between the $\alpha$ values of -6° and 20° is compared for the different $N$, $D$, and $\epsilon$ values in Table 1. The maximum deviations in $C_L$$\left(∆C_L\right)$ and ${\text{C}}_{\text{Dw}}$$\left(∆C_{\text{Dw}}\right)$ are obtained with respect to the results when the $N$, $D$, and $\epsilon$ 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 $\epsilon$ values when it is compared with the result of $N=41$.

Table 1

Effect of $N$, $D$, and $\epsilon$ to the results and $t_{tot}$.

$N$	$D$	$\epsilon$	$∆C_L$(%)	$∆C_{\text{Dw}}$(%)	$t_{\text{tot}}$(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.6

Increasing $D$ decreases the computation time by 86.9% without boosting $∆C_L$ and $∆C_{\text{Dw}}$. When $D$ is selected as 0.4 and above, the convergence problem is observed at the high $\alpha$ values. Contrary to this, increasing $\epsilon$ yields drastic increase in $∆C_L$ especially at the negative $\alpha$ and the high $\alpha$ values. Apart from this region, $∆C_L$ is around 1% at the low $\alpha$ and moderate $\alpha$ values. In light of this information, $N$, $D$, and $\epsilon$ values are selected as 31, 0.2, and 0.01, respectively. $t_{\text{tot}}$ 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 ${10}^6$ 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 Figure 10.

[figure omitted; refer to PDF]

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 calculate $C_l$ and $C_d$ 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 ($\text{TolX}$) or the change in the objective function ($\text{TolFun}$) is less than ${10}^{-3}$. The violation of the constraints ($\text{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 and $P_A$ of the engine is 2000 W. The drag coefficient of the other components ($C_{\text{Doth}}$) of the $\text{UAV}$ except the wing is assumed as a varying parameter as follows: $$\begin{array}{}\text{(20)}& C_{\text{Doth}}=0.02\frac{S_b}{S},\end{array}$$where $S_b$ is the planform area of the wing of the baseline UAV. As a result, the total drag coefficient $\left(C_D\right)$ of the baseline UAV is the summation of $C_{\text{Doth}}$ and $C_{\text{Dw}}$. The total weight of the UAV except for the wing ($W_{\text{oth}}$) is 250 N. The wing weight $\left(W_w\right)$ is calculated as follows [31]: $$\begin{array}{}\text{(21)}& W_w=S\overline{c}{\left(\frac{t}{c}\right)}_{\max }{\rho }_{\text{mat}}{K}_{\rho }{\left(\text{AR}{n}_{\text{ult}}\right)}^{0.6}{\lambda }^{0.04}g,\end{array}$$where $\overline{c}$ denotes the wing mean aerodynamic chord, ${\left(\text{t}/c\right)}_{\max }$ the maximum thickness of the airfoil, $\text{AR}$ the aspect ratio, and $\lambda$ the taper ratio. ${\rho }_{\text{mat}}$ is the density of the construction material that is graphite/epoxy. $K_{\rho }$ is the wing density factor. They are selected as 1575 kg/m³ and 0.0016, respectively, by calculating the average values in the appropriate reference tables for homebuilt aircraft [31]. The relation between ultimate load factor ($n_{\text{ult}}$) and maximum positive load factor ($n_{\max }$) is defined as follows: $$\begin{array}{}\text{(22)}& n_{\text{ult}}=1.5n_{\max }.\end{array}$$

Unlike a fighter aircraft, the producible $n_{\max }$ by the engine power is less than the tolerable $n_{\max }$ by the aircraft structure for a general aviation aircraft [32]. For this reason, the equation of the producible $n_{\max }$ by the engine power for the propeller-driven aircraft that is in equation (23) is used in the $W_w$ calculation. $$\begin{array}{}\text{(23)}& n_{\max }={\left(\frac{C_L^3}{C_D^2}\frac{P_A^2{\rho }_{\infty }S}{2}\right)}^{1/3}\frac{1}{\left(W_{\text{oth}}+W_w\right)}.\end{array}$$

As seen in the above equations, $W_w$ and $n_{\max }$ are implicit. The secant method is applied to find $W_w$.

The geometric parameters of the wing of the baseline UAV are shown in Table 2.

Table 2

Geometric properties of the baseline UAV’s wing.

$\text{Airfoil}$	$\text{NACA}4412$
$$c_r$$	0.45 m
$$c_t$$	0.45 m
$$b$$	4
$$S_b$$	1.8 m2

The wing of the baseline UAV is a rectangular and unswept wing. The airfoil profile is constant along the span. Table 3 depicts the performance parameters of the baseline UAV.

Table 3

Performance parameters of the baseline UAV.

$\max C_L^{3/2}/C_D$	14.08
$\alpha$ at $\max C_L^{3/2}/C_D$	8°
$V$ at $\max C_L^{3/2}/C_D$	15.52 m/s
$M_b$ at $\max C_L^{3/2}/C_D$	127.46 Nm
$$W_w$$	24.06 N
$$V_{\max }$$	39.45 m/s
$\alpha$ at $V_{\max }$	-2.32°
$$V_{\text{stall}}$$	13.25
$\alpha$ at $V_{\text{stall}}$	18°

It has a maximum $C_L^{3/2}/C_D$ ratio of 14.08 when $\alpha$ is 8° and speed is 15.52 m/s. $M_b$ is 127.46 Nm for this flight condition. $W_w$ is 24.06 N. $V_{\max }$ is 39.45 m/s when $\alpha$ is -2.32° and $V_{\text{stall}}$ is 13.25 m/s when $\alpha$ 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. $c_r$, $c_t$, and $b$ are the design variables for the generation of the wing planform. $\alpha$ 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: $$\begin{array}{}\text{(24)}& \min \left(-C_L^{3/2}/C_D\right)\text{at}{\alpha }_1\text{and}{V}_1,\text{(25)}& \text{subject}\text{to}\hspace{1em}-6^{°}\le {\alpha }_1\le {25}^{°},\text{(26)}& 5\text{m}/\text{s}\le V_1\le 60\text{m}/\text{s},\text{(27)}& W_{\text{oth}}+W_w=L\text{at}{\alpha }_1\text{and}{V}_1,\text{(28)}& 0.08\le t\le 0.16,\text{(29)}& 0.15\le p\le 0.6,\text{(30)}& 0.02\le m\le 0.08,\text{(31)}& 0.2\text{m}\le c_r\le 0.6\text{m},\text{(32)}& 0.2\text{m}\le c_t\le 0.6\text{m},\text{(33)}& 4\text{m}\le b\le 8\text{m}.\end{array}$$

The maximization of $C_L^{3/2}/C_D$ is defined in equation (24). The upper limit of the $\alpha$ at the maximum $C_L^{3/2}/C_D$ condition $\left({\alpha }_1\right)$ is 25° whereas the lower limit is -8°. The speed at the maximum $C_L^{3/2}/C_D$ condition $\left(V_1\right)$ is also defined as a design variable so that the aerodynamic design tool calculates the aerodynamic forces for the appropriate $\mathrm{Re}$. Equation (27) is defined to satisfy the level flight condition at the maximum $C_L^{3/2}/C_D$ 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: $$\begin{array}{}\text{(34)}& M_b\le 127.46\text{Nm}\text{atmax}{C}_L^{3/2}/C_D,\text{(35)}& W_w\le 24.06\text{N},\text{(36)}& S\ge 1.8\text{m}.\end{array}$$

Equations (34) and (35) constrain the root bending moment at the maximum $C_L^{3/2}/C_D$ condition and the wing weight in the second optimization problem. In addition to these, equation (36) sets the minimum planform area as 1.8 m².

The equations between Equations (37) and (43) are added into the second optimization problem, and the third optimization problem is obtained. $$\begin{array}{}\text{(37)}& W_{\text{oth}}+W_w=L\text{at}{\alpha }_2\text{and}{V}_2,\text{(38)}& 39.45\text{m}/\text{s}\le V_2\le 60\text{m}/\text{s},\text{(39)}& P_R\le 2000\text{W}\text{at}{\alpha }_2\text{and}{V}_2,\text{(40)}& -6^{°}\le {\alpha }_2\le {25}^{°},\text{(41)}& W_{\text{oth}}+W_w=L\text{at}{\alpha }_3\text{and}{V}_3,\text{(42)}& 5\text{m}/\text{s}\le V_3\le 13.25\text{m}/\text{s},\text{(43)}& -6^{°}\le {\alpha }_3\le {25}^{°}.\end{array}$$

In the third optimization problem, the level flight condition for the $V_2$ and ${\alpha }_2$ values are defined in equation (37). This corresponds to the level flight of the baseline UAV with $V_{\max }$. $V_2$ has a lower limit of 39.45 m/s so that the optimized UAV shall have an equal or greater $V_{\max }$ than the baseline UAV. In addition to this, $P_R$ at $V_2$ and ${\alpha }_2$ is constrained with $P_A$ of the UAV in equation (39). Equation (41) describes the level flight of the baseline UAV with $V_{\text{stall}}$. In the same manner, the upper limit of $V_3$ 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 maximum $C_L^{3/2}/C_D$ is obtained for the first optimization problem when ${\alpha }_1$ is 5.56° and $V_1$ is 11.97 m/s as depicted in Table 4. This corresponds to the neck in the $C_L-\alpha$ curve in Figure 11. Contrary to this, $M_b$ at maximum $C_L^{3/2}/C_D$ nearly doubles and the increment in $W_w$ is about 55%. According to the results in Figure 12, $V_{\max }$ of the first optimization problem is 34.2 m/s that is nearly 5 m/s less than the $V_{\max }$ of the baseline UAV. $V_{\text{stall}}$ of the first optimization problem is 11.33 m/s that is almost the same as its $V_1$. The main reason for this situation is $C_L$ of the maximum $C_L^{3/2}/C_D$ is very close to the maximum lift coefficient $\left({C_L}_{\max }\right)$. 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. $c_t$ and $b$ are set to the lower limit and upper limit, respectively. $c_r$ is decreased from 0.45 m to 0.386 m that alters the taper ratio $\left(\lambda \right)$ from 1 to 0.52.

Table 4

Performance parameters of baseline UAV and the optimum UAVs.

	Baseline UAV	Case 1	Case 2	Case 3
$$\max C_L^{3/2}/C_D$$	14.08	30.79	18.95	18.21
${\alpha }_1\text{atmax}{C}_L^{3/2}/C_D$ (°)	8	5.56	5.69	6.44
$V_1\text{at}\max C_L^{3/2}/C_D$ (m/s)	15.52	11.97	13.41	14.41
$M_b\text{at}\max C_L^{3/2}/C_D$ (Nm)	127.46	249.4	127.3	127.4
$$W_w\left(\text{N}\right)$$	24.06	37.64	19.41	21.8
$$V_2\left(\text{m}/\text{s}\right)$$	39.45			39.45
${\alpha }_2\text{at}{V}_2$ (°)	-2.32			-5.3
$$V_3\left(\text{m}/\text{s}\right)$$	13.25			12.58
${\alpha }_3\text{at}{V}_3$ (°)	18			17.86

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

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
$$c_r\left(\text{m}\right)$$	0.45	0.386	0.586	0.596
$$c_t\left(\text{m}\right)$$	0.45	0.2	0.2	0.2
$$b\left(\text{m}\right)$$	4	8	4.578	4.524
$$S\left({\text{m}}^2\right)$$	1.8	2.34	1.8	1.8

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

The airfoil shape of the second optimization problem is almost identical to the airfoil shape of the first optimization problem. Contrary to this, $c_r$ and $b$ have significant change due to the addition of the functional constraints related with $M_b$ at the maximum $C_L^{3/2}/C_D$ and $S$. $M_b$ reaches its upper limit whereas $S$ is at its lower limit. But $W_w$ is distant from its bound. As a result, the increment at the maximum $C_L^{3/2}/C_D$ is decreased from 220% to 34.6%. $V_{\max }$ of the second optimization problem is 37.2 m/s according to the $P_R-V$ curve in Figure 12 whereas the $V_{\text{stall}}$ is 12.3 m/s. Since $S$ is decreased from 2.34 m² to 1.8 m², the increase in $V_{\text{stall}}$ 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 with $M_b$ at the maximum $C_L^{3/2}/C_D$ and $\text{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. $V_{\max }$ of the third optimization problem is identical to $V_2$ of the same problem according to Table 4 and Figure 12. This implies that the third optimization problem barely has the same $V_{\max }$ 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 $V_2$ speed. $V_{\text{stall}}$ of the third optimization problem is approximately 0.7 m/s less than $V_{\text{stall}}$ of the baseline UAV. According to this result, it can be said that defining a functional constraint related with $V_{\text{stall}}$ in $C_L^{3/2}/C_D$ maximization problems has no significant contribution, because all optimized wings satisfy the level flight condition at a speed less than $V_{\text{stall}}$ of the baseline UAV.

Table 6 depicts the number of iterations $\left(I\right)$, the number of the function count $\left(\text{fc}\right)$, and the elapsed solution time $\left(t_e\right)$ 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. $t_e$ of the third problem is 6.64 hours. If ANSYS-Fluent was used as a function evaluator in the third optimization problem, $t_e$ would be around 500 hours. This comparison shows the cost and time effectiveness of the aerodynamic analysis tool.

Table 6

$I$, $fc$, and $t_e$ of Case 1, Case 2, and Case 3.

	Case 1	Case 2	Case 3
$$I$$	13	18	12
$$\text{fc}$$	32	26	38
$$t_e\left(\text{s}\right)$$	8520	10505	23896

The evaluation of the design variables and the objective function with the iteration number for the optimization problems are presented in Figures 15, 16, and 17.

[figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF]

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 in $C_L^{3/2}/C_D$ 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 airfoil

In 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.

References

[1] I. M. Kroo, S. Altus, R. Braun, P. Gage, I. Sobieski, "Multidisciplinary optimization methods for aircraft preliminary design," 5th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization,DOI: 10.2514/6.1994-4325, .

[2] M. Drela, "XFOIL: an analysis and design system for low Reynolds number airfoils." Low Reynolds number aerodynamics, 1989.

[3] M. Drela, M. B. Giles, "Viscous-inviscid analysis of transonic and low Reynolds number airfoils," AIAA Journal, vol. 25 no. 10, pp. 1347-1355, DOI: 10.2514/3.9789, 1987.

[4] A. Koreanschi, S. G. Oliviu, T. Ayrault, R. M. Botez, M. Mamou, Y. Mebarki, "Numerical optimization and experimental testing of a morphing wing with aileron system," 24th AIAA/AHS Adaptive Structures Conference,DOI: 10.2514/6.2016-1083, .

[5] O. Ş. Gabor, A. Simon, A. Koreanschi, R. M. Botez, "Improving the UAS-S4 Éhecal airfoil high angles-of-attack performance characteristics using a morphing wing approach," Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 230 no. 1, pp. 118-131, DOI: 10.1177/0954410015587725, 2016.

[6] P. Della Vecchia, S. Corcione, R. Pecora, F. Nicolosi, I. Dimino, A. Concilio, "Design and integration sensitivity of a morphing trailing edge on a reference airfoil: the effect on high-altitude long-endurance aircraft performance," Journal of Intelligent Material Systems and Structures, vol. 28 no. 20, pp. 2933-2946, DOI: 10.1177/1045389x17704521, 2017.

[7] A. Magrini, E. Benini, "Aerodynamic optimization of a morphing leading edge airfoil with a constant arc length parameterization," Journal of Aerospace Engineering, vol. 31 no. 2, article 04017093,DOI: 10.1061/(ASCE)AS.1943-5525.0000812, 2017.

[8] Z. Liu, L. Dong, J. M. Moschetta, J. Zhao, G. Yan, "Optimization of nano-rotor blade airfoil using controlled elitist NSGA-II," International Journal of Micro Air Vehicles, vol. 6 no. 1, pp. 29-42, DOI: 10.1260/1756-8293.6.1.29, 2014.

[9] M. A. Silvestre, J. P. Morgado, J. Pascoa, "JBLADE: a propeller design and analysis code," 2013 International Powered Lift Conference,DOI: 10.2514/6.2013-4220, .

[10] J. D. Anderson, S. Corda, "Numerical lifting line theory applied to drooped leading-edge wings below and above stall," Journal of Aircraft, vol. 17 no. 12, pp. 898-904, DOI: 10.2514/3.44690, 1980.

[11] P. Gamboa, J. Vale, F. J. P. Lau, A. Suleman, "Optimization of a morphing wing based on coupled aerodynamic and structural constraints," AIAA Journal, vol. 47 no. 9, pp. 2087-2104, DOI: 10.2514/1.39016, 2009.

[12] C. B. Merz, C. C. Wolf, K. Richter, K. Kaufmann, A. Mielke, M. Raffel, "Spanwise differences in static and dynamic stall on a pitching rotor blade tip model," Journal of the American Helicopter Society, vol. 62 no. 1,DOI: 10.4050/jahs.62.012002, 2017.

[13] J. D. Anderson, Fundamentals of aerodynamics, 2010.

[14] J. Lee, K. B. Tan, P. C. Wang, "Parametric investigation of box wing configuration in viscous flow regime," 2018 Applied Aerodynamics Conference,DOI: 10.2514/6.2018-3814, .

[15] https://web.mit.edu/drela/Public/web/xfoil/

[16] https://www.mathworks.com/matlabcentral/fileexchange/30446-xfoil-interface

[17] D. J. Pate, B. J. German, "Lift distributions for minimum induced drag with generalized bending moment constraints," Journal of Aircraft, vol. 50 no. 3, pp. 936-946, DOI: 10.2514/1.c032074, 2013.

[18] B. Sterligov, S. Cherkasov, "Reducing magnetic noise of an unmanned aerial vehicle for high-quality magnetic surveys," International Journal of Geophysics, vol. 2016,DOI: 10.1155/2016/4098275, 2016.

[19] M. Hassanalian, G. Throneberry, M. Ali, S. B. Ayed, A. Abdelkefi, "Role of wing color and seasonal changes in ambient temperature and solar irradiation on predicted flight efficiency of the albatross," Journal of Thermal Biology, vol. 71, pp. 112-122, DOI: 10.1016/j.jtherbio.2017.11.002, 2018.

[20] D. Communier, M. F. Salinas, O. Carranza Moyao, R. M. Botez, "Aero structural modeling of a wing using CATIA V5 and XFLR5 software and experimental validation using the Price-Païdoussis wing tunnel," AIAA Atmospheric Flight Mechanics Conference, .

[21] R. B. Langtry, F. R. Menter, "Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes," AIAA journal, vol. 47 no. 12, pp. 2894-2906, DOI: 10.2514/1.42362, 2009.

[22] S. M. A. Aftab, A. M. Rafie, N. A. Razak, K. A. Ahmad, "Turbulence model selection for low Reynolds number flows," PloS one, vol. 11 no. 4, article e0153755,DOI: 10.1371/journal.pone.0153755, 2016.

[23] R. Lanzafame, S. Mauro, M. Messina, "Wind turbine CFD modeling using a correlation-based transitional model," Renewable Energy, vol. 52, pp. 31-39, DOI: 10.1016/j.renene.2012.10.007, 2013.

[24] ANSYS, Fluent Theory Guide, 2013.

[25] J. L. Thomas, M. D. Salas, "Far-field boundary conditions for transonic lifting solutions to the Euler equations," AIAA Journal, vol. 24 no. 7, pp. 1074-1080, DOI: 10.2514/3.9394, 1986.

[26] F. R. Menter, R. B. Langtry, S. R. Likki, Y. B. Suzen, P. G. Huang, S. Völker, "A correlation-based transition model using local variables—part I: model formulation," Journal of Turbomachinery, vol. 128 no. 3, pp. 413-422, DOI: 10.1115/1.2184352, 2006.

[27] https://web.mit.edu/drela/Public/web/xfoil/xfoil6.99.tgz

[28] E. R. Prytz, Ø. Huuse, B. Müller, J. Bartl, L. R. Sætran, "Numerical simulation of flow around the NREL S826 airfoil at moderate Reynolds number using delayed detached Eddy simulation (DDES)," AIP Conference Proceedings, vol. 1863 no. 1, article 560086,DOI: 10.1063/1.4992769, 2017.

[29] G. N. Vanderplaats, Numerical optimization techniques for engineering design, vol. 1999, 1984.

[30] L. S. Lasdon, A. D. Waren, "Large scale nonlinear programming," Computers & Chemical Engineering, vol. 7 no. 5, pp. 595-604, DOI: 10.1016/0098-1354(83)80005-2, 1983.

[31] M. H. Sadraey, Aircraft design: a systems engineering approach, 2012.

[32] M. Sadraey, Aircraft performance: analysis, 2009.