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 is the objective function, while level flight condition at maximum 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.
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 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 and the profile drag coefficient 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 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 was 6000. Xfoil captured and 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 shear stress transport (SST) model implemented into ANSYS-Fluent. After that, three different optimization problems in which maximization is the objective function are constructed. The optimization problems are defined by adding new functional constraints that are the wing root bending moment , the wing weight , the wing planform area , the level flight condition and power available at the maximum speed , and the level flight condition at the stall speed of the baseline unmanned air vehicle . 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 at both and 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 of the finite wing into number of spanwise stations in the y-direction as shown in Figure 1 [13]. In Figure 1, the distance between the equidistant stations is . The method needs an initial circulation distribution along the span. The elliptic circulation distribution in equation (1) is acceptable for this requirement.
In the above equation, is the circulation at the root chord . By using initial , the induced angle of attack for the panel is obtained as follows:
This integral is evaluated numerically by using Simpson’s rule [13]: where is the freestream velocity. For the calculation of , is needed. Although 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.
In order to avoid singularity when , , or , the average of the neighbor stations’ results is taken as the result of the corresponding term [13]. After the calculation of , the effective angle of attack is obtained by subtracting from the geometric angle of attack .
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 [14]. That implies that values of the airfoils at the tip chords of the wing are zero and is equal to zero lift angle of attack value of the airfoils at the tips as shown in
Once distribution along the span is known, and 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 is calculated by using
in the above equation is the chord value at the station. The input circulation values are calculated by using a damping factor as follows:
The iterative solution requires a strong damping factor. Setting 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 :
By using converged circulation distribution, the lift [10], the induced drag [10], and the root bending moment [17] of the wing are calculated as follows: where is the free stream density. These integrals are also evaluated numerically by using Simpson’s rule. The parasite drag of the wing is calculated by using
After the calculations of the aerodynamic forces and moments, the lift coefficient , the drag coefficient , and the root bending moment coefficient of the wing are obtained as follows:
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 SST Model
The XFLR5 LLT and ANSYS-Fluent 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 and 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 and 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 of 120000 with the experimental results: Spalart-Allmaras, SST, intermittency SST, , and SST [22]. According to the results, only the 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 SST and SST with experimental data. According to the result, SST captured the trend of the aerodynamic coefficients whereas SST overpredicted values and underpredicted values.
According to these studies in the literature, the computations are performed with 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 is selected as 1 [26]. Thereby, first-layer height is calculated as 0.015 mm for the selected . 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 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 and lie within the range of .
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 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 and 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 and the camber generation , respectively. In these equations, is the maximum thickness in hundredths of the chord, is the maximum camber location in tenths of the chord, and is the maximum camber in hundredths of the chord. coordinates are generated by the cosine distribution.
The 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].
, , and are , , and , respectively. is . The dynamic viscosity is and is 22 . In Xfoil, the default critical value of the 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 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 at the same . values of , , and 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 . values of ANSYS-Fluent are almost identical to case values when is between 0° and 10° as depicted in Figure 5. But as gets close to 15°, the ANSYS-Fluent, case, and case yield similar . 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.
In the aerodynamic analysis tool, the transition region can be determined by finding the region at which a sudden increase in the -wall shear stress is observed. The turbulent region is the region at which the -wall shear stress values get closer to zero. Contrary to this, the transition region and the turbulent separation region are identified when -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 and cases capture the values of ANSYS-Fluent at high values, they are selected as the candidate station numbers.
In order to select the optimum , , and values, the total solution time to obtain drag polar results between the values of -6° and 20° is compared for the different , , and values in Table 1. The maximum deviations in and are obtained with respect to the results when the , , and values are equal to 41, 0.05, and 0.01, respectively. According to the results, choosing the value as 31 saves the computation time approximately by 38% at the same and values when it is compared with the result of .
Increasing decreases the computation time by 86.9% without boosting and . When is selected as 0.4 and above, the convergence problem is observed at the high values. Contrary to this, increasing yields drastic increase in especially at the negative and the high values. Apart from this region, is around 1% at the low and moderate values. In light of this information, , , and values are selected as 31, 0.2, and 0.01, respectively. 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 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.
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 and 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 () or the change in the objective function () is less than . The violation of the constraints () is limited to .
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 of the engine is 2000 W. The drag coefficient of the other components () of the except the wing is assumed as a varying parameter as follows: where is the planform area of the wing of the baseline UAV. As a result, the total drag coefficient of the baseline UAV is the summation of and . The total weight of the UAV except for the wing () is 250 N. The wing weight is calculated as follows [31]: where denotes the wing mean aerodynamic chord, the maximum thickness of the airfoil, the aspect ratio, and the taper ratio. is the density of the construction material that is graphite/epoxy. 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 () and maximum positive load factor () is defined as follows:
Unlike a fighter aircraft, the producible by the engine power is less than the tolerable by the aircraft structure for a general aviation aircraft [32]. For this reason, the equation of the producible by the engine power for the propeller-driven aircraft that is in equation (23) is used in the calculation.
As seen in the above equations, and are implicit. The secant method is applied to find .
The geometric parameters of the wing of the baseline UAV are shown in Table 2.
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.
It has a maximum ratio of 14.08 when is 8° and speed is 15.52 m/s. is 127.46 Nm for this flight condition. is 24.06 N. is 39.45 m/s when is -2.32° and is 13.25 m/s when is .
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, , , and are defined as the design variables for airfoil generation. , , and are the design variables for the generation of the wing planform. and 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:
The maximization of is defined in equation (24). The upper limit of the at the maximum condition is 25° whereas the lower limit is -8°. The speed at the maximum condition is also defined as a design variable so that the aerodynamic design tool calculates the aerodynamic forces for the appropriate . Equation (27) is defined to satisfy the level flight condition at the maximum 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:
Equations (34) and (35) constrain the root bending moment at the maximum 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.
In the third optimization problem, the level flight condition for the and values are defined in equation (37). This corresponds to the level flight of the baseline UAV with . has a lower limit of 39.45 m/s so that the optimized UAV shall have an equal or greater than the baseline UAV. In addition to this, at and is constrained with of the UAV in equation (39). Equation (41) describes the level flight of the baseline UAV with . In the same manner, the upper limit of 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 is obtained for the first optimization problem when is 5.56° and is 11.97 m/s as depicted in Table 4. This corresponds to the neck in the curve in Figure 11. Contrary to this, at maximum nearly doubles and the increment in is about 55%. According to the results in Figure 12, of the first optimization problem is 34.2 m/s that is nearly 5 m/s less than the of the baseline UAV. of the first optimization problem is 11.33 m/s that is almost the same as its . The main reason for this situation is of the maximum is very close to the maximum lift coefficient . 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. and are set to the lower limit and upper limit, respectively. is decreased from 0.45 m to 0.386 m that alters the taper ratio from 1 to 0.52.
The airfoil shape of the second optimization problem is almost identical to the airfoil shape of the first optimization problem. Contrary to this, and have significant change due to the addition of the functional constraints related with at the maximum and . reaches its upper limit whereas is at its lower limit. But is distant from its bound. As a result, the increment at the maximum is decreased from 220% to 34.6%. of the second optimization problem is 37.2 m/s according to the curve in Figure 12 whereas the is 12.3 m/s. Since is decreased from 2.34 m2 to 1.8 m2, the increase in 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 at the maximum and 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. of the third optimization problem is identical to of the same problem according to Table 4 and Figure 12. This implies that the third optimization problem barely has the same 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 speed. of the third optimization problem is approximately 0.7 m/s less than of the baseline UAV. According to this result, it can be said that defining a functional constraint related with in maximization problems has no significant contribution, because all optimized wings satisfy the level flight condition at a speed less than of the baseline UAV.
Table 6 depicts the number of iterations , the number of the function count , and the elapsed solution time 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. of the third problem is 6.64 hours. If ANSYS-Fluent was used as a function evaluator in the third optimization problem, would be around 500 hours. This comparison shows the cost and time effectiveness of the aerodynamic analysis tool.
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.
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 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.