Development of Response Surface Model of Endurance Time and Structural Parameter Optimization for a Tailsitter UAV

This study designed a vertical take-off and landing tailsitter unmanned aerial vehicle (UAV) with a long endurance time. Nine parameters of the tailsitter UAV were investigated. Using a 2k full factorial test, 512 experiments on the nine parameters were conducted at their maximum and minimum values. The time coefficient and air resistance were calculated using the computational fluid dynamics (CFD) method under different parameter combinations. The analysis of variance determined that the specific factors influencing the time coefficient and air resistance were the root chord, wingtip chord, wingspan, and sweep angle. By carrying out a central composite design (CCD) test, 25 sample points of the four particular factors were constructed. The time coefficient and air resistance were simulated under different structural parameter combinations using the CFD method. CFD simulation was verified by carrying out a wind tunnel test, and the results revealed that the aerodynamic coefficient error was less than 5%, while the air resistance error was less than 6%. The response surface methodology (RSM) for the time coefficient and air resistance was established using a genetic aggregation method. A multi-objective genetic algorithm (MOGA) was used to optimize the parameters with regard to the maximum time coefficient and minimum air resistance. The optimal structural parameters were wing root chord length at 315 mm, wingtip chord length at 182 mm, wingspan length at 1198 mm, and sweep angle at 16°. Compared with the original layout and size, the time coefficient of the new design of the tailsitter UAV improved by 19.5%, while the air resistance reduced by 34.78%. The results obtained by this study are significant for the design of tailsitter UAVs.


Introduction
The tailsitter unmanned aerial vehicle (UAV) has multi-rotor vertical take-off and landing (VTOL) characteristics and fixed-wing dynamic cruising characteristics. These characteristics solve the low efficiency and short-range of failures in multi-rotor UAVs, the lengthy fixed-wing UAV placement preparation, auxiliary equipment requirement, and various landing problems. Additionally, tailsitter UAVs can satisfy the requirements of increasingly complex agricultural low-altitude remote sensing tasks [1][2][3][4][5]. Compared with other types of VTOL vehicles, the tailsitter UAV does not need additional rotating control units and has the advantages of a compact structure, a light weight, high overall natural stability, and smooth operation.
Many studies used computational fluid dynamics (CFD) in tailsitter UAV design and aerodynamic optimization [6][7][8][9][10][11]. In 2014, Ang et al. designed a foldable drone U-lion using a CFD numerical flight space, wing load, and batteries to optimize the flight parameters of electric-driven tailsitter drones [31]. In 2015, Wang et al. used the battery discharge curve and flight parameters to estimate spatial distance [32]. The structure and layout of a tailsitter drone is an essential factor that affects its battery life. However, existing research cannot explain the relationship between battery life and structural parameters.
With regard to the design of a long endurance tailsitter UAV, this study investigated the relationship between structural parameters and endurance time using a genetic clustering algorithm, and used a multi-objective genetic algorithm to optimize the drone's structural parameters.

UAV Structure Layout
The tailsitter UAV adopts a double-pull forward layout with flying wings, as shown in Figure 1. It has a left-right symmetrical structure, including the wings, winglets, rudder surface, motor seats, motors, and propellers. When the drone takes off, it performs a vertical ascent in the multi-rotor flying mode, and the rudder surface performs attitude control and adjustment. After reaching a predetermined height, the drone enters into the transition flying mode, the motor speed increases, and the rudder surface provides a lowering torque for transformation into a fixed-wing flying mode. Then, the UAV collects the remote sensing images of the working area while cruising. Before landing, the motor speed decreases to reduce the flying velocity. When it reaches a certain speed, the rudder surface provides a head-up torque to pull the nose. During landing, the drone enters into the multi-rotor flying mode again. flight space, wing load, and batteries to optimize the flight parameters of electric-driven tailsitter drones [31]. In 2015, Wang et al. used the battery discharge curve and flight parameters to estimate spatial distance [32]. The structure and layout of a tailsitter drone is an essential factor that affects its battery life. However, existing research cannot explain the relationship between battery life and structural parameters. With regard to the design of a long endurance tailsitter UAV, this study investigated the relationship between structural parameters and endurance time using a genetic clustering algorithm, and used a multi-objective genetic algorithm to optimize the drone's structural parameters.

UAV Structure Layout
The tailsitter UAV adopts a double-pull forward layout with flying wings, as shown in Figure  1. It has a left-right symmetrical structure, including the wings, winglets, rudder surface, motor seats, motors, and propellers. When the drone takes off, it performs a vertical ascent in the multi-rotor flying mode, and the rudder surface performs attitude control and adjustment. After reaching a predetermined height, the drone enters into the transition flying mode, the motor speed increases, and the rudder surface provides a lowering torque for transformation into a fixed-wing flying mode. Then, the UAV collects the remote sensing images of the working area while cruising. Before landing, the motor speed decreases to reduce the flying velocity. When it reaches a certain speed, the rudder surface provides a head-up torque to pull the nose. During landing, the drone enters into the multirotor flying mode again. During the whole flight course, the fixed-wing flying mode accounts for more than 80% of the total battery life. In this study, the development of a response surface model of endurance time for the tailsitter drone was mainly during the fixed-wing flying mode.

Structural Parameter Range
The main structural parameters ( Figure 2) affecting the drone's endurance time are the wingspan b, root chord cr, wingtip chord ct, wing sweep angle Λw, winglet length lys, winglet span bv, winglet height lv, winglet swept angle Λv, winglet thickness h, and winglet foot length ljc (the distance between the two legs of the winglet). During the whole flight course, the fixed-wing flying mode accounts for more than 80% of the total battery life. In this study, the development of a response surface model of endurance time for the tailsitter drone was mainly during the fixed-wing flying mode.

Structural Parameter Range
The main structural parameters ( Figure 2) affecting the drone's endurance time are the wingspan b, root chord c r , wingtip chord c t , wing sweep angle Λ w , winglet length l ys , winglet span b v , winglet height l v , winglet swept angle Λ v , winglet thickness h, and winglet foot length l jc (the distance between the two legs of the winglet). Switzerland's senseFly and senseFly eBee, France's Parrot disco-pro, and China's skywalker X5 are all electric driven fixed-wing remote sensing UAVs currently in agriculture. The UAV parameters obtained from the manufacturer's official website are listed in Table 1. The range of b can be determined by the wingspan range of the existing UAV presented in Table  1. The scope of cr (240-500 mm) was determined by considering the embedded flight control installation, digital transmission system, and other equipment. The range of ct, Λw, Λv, and lv can be carried out according to the small flying wing drone's taper ratio (0.6) [33], sweepback angle (0-60°) [34], winglet sweep angle (30-60°), and winglet height (0.05b-0.3b) range [35,36]. The expanded polypropylene (EPP) material stiffness is considered to determine the h range. The range of bv, lys, and ljc was selected to determine the longitudinal stability of the UAV's VTOL timing body. The sensitivity analysis results showed that the influences of parameter bv on air resistance and the time coefficient are very small. The sensitivity coefficient of air resistance is 0.4% and 0.2% the time coefficient. Therefore, parameter bv was set to a fixed value (180 mm) to reduce computation load. The final structural parameter range is presented in Table 2.

Governing Equations and Numerical Method
The tailsitter UAV combines the advantages of a multi-rotor VTOL and fixed-wing efficient cruising but also has the disadvantages of short endurance time and poor attitude conversion stability. Therefore, improving the aerodynamic efficiency, extending the endurance time, reducing Switzerland's senseFly and senseFly eBee, France's Parrot disco-pro, and China's skywalker X5 are all electric driven fixed-wing remote sensing UAVs currently in agriculture. The UAV parameters obtained from the manufacturer's official website are listed in Table 1. The range of b can be determined by the wingspan range of the existing UAV presented in Table 1. The scope of c r (240-500 mm) was determined by considering the embedded flight control installation, digital transmission system, and other equipment. The range of c t , Λ w , Λ v , and l v can be carried out according to the small flying wing drone's taper ratio (0.6) [33], sweepback angle (0-60 • ) [34], winglet sweep angle (30-60 • ), and winglet height (0.05b-0.3b) range [35,36]. The expanded polypropylene (EPP) material stiffness is considered to determine the h range. The range of b v , l ys , and l jc was selected to determine the longitudinal stability of the UAV's VTOL timing body. The sensitivity analysis results showed that the influences of parameter b v on air resistance and the time coefficient are very small. The sensitivity coefficient of air resistance is 0.4% and 0.2% the time coefficient. Therefore, parameter b v was set to a fixed value (180 mm) to reduce computation load. The final structural parameter range is presented in Table 2.

Governing Equations and Numerical Method
The tailsitter UAV combines the advantages of a multi-rotor VTOL and fixed-wing efficient cruising but also has the disadvantages of short endurance time and poor attitude conversion stability. Therefore, improving the aerodynamic efficiency, extending the endurance time, reducing the air resistance, reducing the lateral slip, and enhancing the security of the attitude conversion are key Sensors 2020, 20, 1766 5 of 21 considerations for optimizing the UAV's structural parameters. Additionally, the CFD numerical simulation method can be used to simulate the external flow field of the UAV.
The equations governing the continuity, momentum, and energy in the computational domain can be expressed as follows: Continuity: Momentum: Energy: The k-ω-based shear stress transport (SST) model provides highly accurate predictions for the onset and amount of flow separation under adverse pressure gradients by considering the transport of turbulent shear stress. k-equation: ω-equation: In addition to the independent variables, the density ρ and velocity vector U are treated as known quantities in the Navier-Stokes method; P k is the production rate of turbulence and is calculated in the same manner as in the k − ε model. The model constants are as follows: β = 0.09, α = 5/9, and β = 0.075.
The proper transport behavior can be obtained by applying a limiter to the formulation of eddy-viscosity, as follows: where v t = µ t ρ .
Again, F 2 is a blending function restricting the limiter to the wall boundary layer because the underlying assumptions are not valid for free shear flows; S is an invariant measure of the strain rate.
The three-dimensional model of the UAV was built by computer aided tri-dimensional interface application (CATIA), which is widely used in aircraft design [37]. Owing to the symmetrical layout of the UAV, a three-dimensional model was built for the left side of the fuselage, and the structures of the propeller, motor, and motor seat were simplified. The UAV 3D model was imported into the ansys (ANSYS) geometry module to create the outflow field. The outflow field was constructed with the dimensions of 5 × 5 × 3 m, and the three-dimensional solid model of the outflow field was obtained by a Boolean operation [38]. The grid of the outflow field was generated with the ICEM software [39,40]. An unstructured grid was used for grid generation. Additionally, to capture the flow details of the Sensors 2020, 20, 1766 6 of 21 UAV close to the wall and outside the flow field, the grid around the UAV was locally encrypted, as shown in Figure 3a. The number of the final generated grid of the outflow field was 2.3 million nodes, and the grid close to the wall of the UAV is shown in Figure 3b.
of the UAV, a three-dimensional model was built for the left side of the fuselage, and the structures of the propeller, motor, and motor seat were simplified. The UAV 3D model was imported into the ansys (ANSYS) geometry module to create the outflow field. The outflow field was constructed with the dimensions of 5 × 5 × 3 m, and the three-dimensional solid model of the outflow field was obtained by a Boolean operation [38]. The grid of the outflow field was generated with the ICEM software [39,40]. An unstructured grid was used for grid generation. Additionally, to capture the flow details of the UAV close to the wall and outside the flow field, the grid around the UAV was locally encrypted, as shown in Figure 3a. The number of the final generated grid of the outflow field was 2.3 million nodes, and the grid close to the wall of the UAV is shown in Figure 3b. Fluent is a common CFD software that is widely used in the optimal design of aircraft and can accurately calculate aerodynamic parameters. The lift coefficient CL, drag coefficient CD, and air resistance of the UAV were calculated by Fluent with the external flow field [41][42][43][44].
The boundary conditions in Fluent software were set as follows: the inlet speed was 12 m/s, the included angle was 4-12°, and the outlet pressure was standard atmospheric pressure. The working fluid was air, and its physical properties were considered constant. For the boundary of the calculation area, no-slip conditions and the standard law for the wall were assumed. The semiimplicit method for pressure-linked equations (SIMPLE) method in Fluent software was adopted for the coupling of pressure and velocity. The convergence criterion was considered as normalized residuals being less than 1 × 10 -4 for the flow equations.

Theoretical Model of Endurance Time
The endurance time is mainly related to the battery power, take-off weight, and aerodynamic characteristics. According to Tan [45], the endurance time can be expressed as follows: Fluent is a common CFD software that is widely used in the optimal design of aircraft and can accurately calculate aerodynamic parameters. The lift coefficient C L , drag coefficient C D , and air resistance of the UAV were calculated by Fluent with the external flow field [41][42][43][44].
The boundary conditions in Fluent software were set as follows: the inlet speed was 12 m/s, the included angle was 4-12 • , and the outlet pressure was standard atmospheric pressure. The working fluid was air, and its physical properties were considered constant. For the boundary of the calculation area, no-slip conditions and the standard law for the wall were assumed. The semi-implicit method for pressure-linked equations (SIMPLE) method in Fluent software was adopted for the coupling of pressure and velocity. The convergence criterion was considered as normalized residuals being less than 1 × 10 −4 for the flow equations.

Theoretical Model of Endurance Time
The endurance time is mainly related to the battery power, take-off weight, and aerodynamic characteristics. According to Tan [45], the endurance time can be expressed as follows: In this study, The endurance time can be expressed as follows: The meaning of each variable: ρ-air density,1.185kg/m 3 ; s-effective projected area, m 2 ; G-UAV weight, N; C D -drag coefficient; C L -lift coefficient.

Feature Factor Extraction Method
According to Equation (8), the battery power Q and time coefficient determine the endurance time. During the design of the UAV, the take-off weight and battery power were generally fixed. Therefore, the endurance time was mainly related to the projection area s, and the lift coefficient and Sensors 2020, 20, 1766 7 of 21 drag coefficient were related. An UAV typically employs a tapered wing, and the projection area is primarily associated with the root chord length, wingtip chord length, and wingspan length.
The two levels of the factor in the 2 k design are usually represented as -1 (for the first level) and 1 (for the second level). Note that this representation is reversed from the coding used in general full factorial designs for the indicator variables that represent two level factors in analysis of variance (ANOVA) models [46]. For ANOVA models, the first level of the factor is represented using a value of 1 for the indicator variable, while the second level is represented using a value of -1. When k = 3, the design can be represented geometrically using a cube with the eight treatment combinations lying at the eight corners, as shown in Figure 4.

Feature Factor Extraction Method
According to Equation (8), the battery power Q and time coefficient determine the endurance time. During the design of the UAV, the take-off weight and battery power were generally fixed. Therefore, the endurance time was mainly related to the projection area s, and the lift coefficient and drag coefficient were related. An UAV typically employs a tapered wing, and the projection area is primarily associated with the root chord length, wingtip chord length, and wingspan length.
The two levels of the factor in the 2 k design are usually represented as -1 (for the first level) and 1 (for the second level). Note that this representation is reversed from the coding used in general full factorial designs for the indicator variables that represent two level factors in analysis of variance (ANOVA) models [46]. For ANOVA models, the first level of the factor is represented using a value of 1 for the indicator variable, while the second level is represented using a value of -1. When k = 3, the design can be represented geometrically using a cube with the eight treatment combinations lying at the eight corners, as shown in Figure 4. A two-level factorial design is a type of multi-factor cross-group experimental design [47]. In 1993, Moore and Dalva applied a 2 k test for the first time to determine the influence of temperature and groundwater level on the CO2 and CH4 of soil. Not only were the effects of the two variables on gas emissions investigated but also the effects of their interaction on emissions [48]. Fogue et al. used a statistical analysis based on the 2 k factorial methodology to determine the most representative A two-level factorial design is a type of multi-factor cross-group experimental design [47]. In 1993, Moore and Dalva applied a 2 k test for the first time to determine the influence of temperature and groundwater level on the CO 2 and CH 4 of soil. Not only were the effects of the two variables on gas emissions investigated but also the effects of their interaction on emissions [48]. Fogue et al. used a statistical analysis based on the 2 k factorial methodology to determine the most representative factors that affect traffic safety applications under real roadmaps [49]. The 2 k factorial design can test not only the differences between the levels of each factors but also check the interactions between elements [50][51][52].
In this study, a 2 k factorial design was used to conduct a full combination test for the maximum and minimum values of the nine structural parameters to obtain 512 sets of experimental data. A variance analysis was carried out for the main effect and second-order interaction effect of each factor to determine their impacts on the endurance time coefficient.
In the half normal plot, the effects with an absolute value exceeding the margin of error (ME) are labeled as significant. The Minitab software was used to calculate the ME. The margin of error is expressed as follows: where t is the (1-α/2) quantile of a t-distribution. Lenth's pseudo standard error (PSE) is based on the concept of sparse effects, which assumes that the variation in the smallest effects is caused by a random error. The PSE was calculated in the Minitab software as follows: 1.
To calculate the absolute value of the effects; 2.
To calculate ∆, which is 1.5× the median of the effects in step 1; 3.
To calculate the median of effects that are less than 2.5 × ∆; 4.
To calculate the PSE, which is 1.5× the median calculated in step 3.

Endurance Time Response Surface
Central composite design is a five-level fractional factorial design that is suitable for calibrating the quadratic response model [53]. There are three types of central composite designs (CCDs) that are commonly used in experiment designs: circumscribed, inscribed, and face-centered CCDs. The five-level coded values of each factor are represented by (−a, −1, 0, +1, +a), where (−1, +1) corresponds to the physical lower and upper limit of the explored factor space. It is obvious that (−a, +a) establishes new "extreme" physical lower and upper limits for all factors. Figure 5 is a geometrical representation of a circumscribed CCD of two factors: concept of sparse effects, which assumes that the variation in the smallest effects is caused by a random error. The PSE was calculated in the Minitab software as follows: 1. To calculate the absolute value of the effects; 2. To calculate Δ, which is 1.5× the median of the effects in step 1; 3. To calculate the median of effects that are less than 2.5 × Δ; 4. To calculate the PSE, which is 1.5× the median calculated in step 3.

Endurance Time Response Surface
Central composite design is a five-level fractional factorial design that is suitable for calibrating the quadratic response model [53]. There are three types of central composite designs (CCDs) that are commonly used in experiment designs: circumscribed, inscribed, and face-centered CCDs. The fivelevel coded values of each factor are represented by (-a, -1, 0, +1, +a), where (-1, +1) corresponds to the physical lower and upper limit of the explored factor space. It is obvious that (-a, +a) establishes new "extreme" physical lower and upper limits for all factors. Figure 5 is a geometrical representation of a circumscribed CCD of two factors: The central composite design can analyze the importance of factors and the interactive responses between factors. CCD is widely used in chemical reagent ratios and processing technology optimization [54,55]. In 2006, Idris et al. conducted a central composite design study and built a response surface methodology to identify the essential interfacial reaction factors that influence membrane performance [56]. In the CCD test, each variable has five levels that can be used to fit the response surface and optimize the relevant structural parameters. In 2016, Wei et al. combined the response surface methodology (RSM) and CCD to optimize the structural parameters of a multiobjective analysis and achieved dynamic performance improvement of the cross-spring compliant microdisplacement magnifying mechanism [57].
The four essential influence factors for the time coefficient are root chord, wingtip chord, wingspan, and sweep angle as discussed in Section 3.1. The CCD was used to design the four structural parameters and obtain 16 cubic points, eight axial points, and one center point, which The central composite design can analyze the importance of factors and the interactive responses between factors. CCD is widely used in chemical reagent ratios and processing technology optimization [54,55]. In 2006, Idris et al. conducted a central composite design study and built a response surface methodology to identify the essential interfacial reaction factors that influence membrane performance [56]. In the CCD test, each variable has five levels that can be used to fit the response surface and optimize the relevant structural parameters. In 2016, Wei et al. combined the response surface methodology (RSM) and CCD to optimize the structural parameters of a multi-objective analysis and achieved dynamic performance improvement of the cross-spring compliant microdisplacement magnifying mechanism [57].
The four essential influence factors for the time coefficient are root chord, wingtip chord, wingspan, and sweep angle as discussed in Section 3.1. The CCD was used to design the four structural parameters and obtain 16 cubic points, eight axial points, and one center point, which amounts to a total of 25 combined sample points. The CFD method was used to calculate the time coefficient and air resistance values of the 25 combined sample points, and the genetic aggregation response surface method was used to construct the response surface models of the time coefficient t FA and air resistance D.
Four types of metamodels (polynomial regression, Kriging regression, support vector regression, and moving least squares) with different settings were used as the initial metamodels to build the response surface model. To increase the chance of obtaining the most effective response surface, DesignXplorer generated a population of metamodels with different types and settings [58,59]. The next population was obtained by the cross-over and mutation of previous populations using the genetic algorithm run by DesignXplorer.
The genetic aggregation generated response surface can be expressed as an ensemble using the weighted average of different metamodels, as follows: whereŷ ens is the prediction of the ensemble,ŷ i is the prediction of the i-th response surface, N M is the number of metamodels used, and w i is the weight factor of the i-th response surface. The weight factors satisfy the following conditions:

Multi-Objective Genetic Algorithm (MOGA)
The tailsitter UAV collects remote sensing images of the target area in the fixed-wing cruising mode. During the design process, it should be ensured that the drone can achieve a long endurance time. During the attitude conversion process, the pitch angle gradually increases from 0 • to 90 • . To ensure the stability of the attitude conversion process, the air resistance of the drone should be as small as possible: where D is the air resistance of the UAV; x is a structural parameter; x d is the structural parameter range; and m, n, i, j, and d denote the number of invariant states. The MOGA used in this paper was NSGA-II (non-dominated sorted genetic algorithm-II) [60]. The parameters and operators of the NSGA-II were as follows: the population size was 100, the number of generations was 100, the mutation rate was 0.01, the crossover rate was 0.98, the type of crossover was one point, and the type of mutation was simple mutation.
The response surface and optimized MOGA were built with DesignXplorer, which described the relationship between the design variables and the performance of the product by using design of experiments (DOE), combined with the response surfaces [61]. In the DesingXplorer software, the initial population was constructed by the MOGA module, the data of the constructed population were transmitted to the RSM module to solve the endurance time and air resistance, and the solution result was transmitted to the MOGA module for convergence verification [62]. If the requirements were not met, the population crossover and mutation were carried out by the MOGA for a new iterative calculation, which was iterated until finding the optimal individual solution with the longest endurance time and least air resistance, which are the optimal structural parameters.

Test Materials and Equipment
A wind tunnel test was carried out to verify the CFD numerical simulation results. A sample point (wingspan 1200 mm, sweep angle 28 • , root chord length 500 mm, and wingtip chord length 280 mm) was randomly selected from 25 CCD combination samples. The aerodynamic coefficients C FA and D of the prototype were measured under different wind speeds. Additionally, 3D printing technology was used to process the sample (as shown in Figure 6); the wings of the sample were scaled to 600 mm according to the similar criteria used in wind tunnel tests.
A wind tunnel test was carried out to verify the CFD numerical simulation results. A sample point (wingspan 1200 mm, sweep angle 28°, root chord length 500 mm, and wingtip chord length 280 mm) was randomly selected from 25 CCD combination samples. The aerodynamic coefficients CFA and D of the prototype were measured under different wind speeds. Additionally, 3D printing technology was used to process the sample (as shown in Figure 6); the wings of the sample were scaled to 600 mm according to the similar criteria used in wind tunnel tests.

Test Conditions and Scheme
The wings and winglets were produced by 3D printing technology and assembled into a UAV sample. The grooves and gaps of the wings were skinned to ensure the surface smoothness. The wind tunnel test system includes aircraft support, six-component strain balance, fan, frequency converter, and data analysis system, as shown in Figure 7. A six-component strain balance was embedded into the UAV body to reduce the interference of airflow to the balance, and the UAV body was fixed with an extension strut to minimize the intervention of the support to the tail airflow. During the test, the UAV was kept at the center of the wind field (1.4 m above the ground) to obtain the best wind field data, and the inverter was used to control the fan speed. After the wind speed was stable, the lift coefficient and drag coefficient were determined by the six-component strain balance. The data analysis system analyzed 1000 data sets and finally obtained the aerodynamic coefficients and air resistance of the drone.

Test Conditions and Scheme
The wings and winglets were produced by 3D printing technology and assembled into a UAV sample. The grooves and gaps of the wings were skinned to ensure the surface smoothness. The wind tunnel test system includes aircraft support, six-component strain balance, fan, frequency converter, and data analysis system, as shown in Figure 7. A six-component strain balance was embedded into the UAV body to reduce the interference of airflow to the balance, and the UAV body was fixed with an extension strut to minimize the intervention of the support to the tail airflow. During the test, the UAV was kept at the center of the wind field (1.4 m above the ground) to obtain the best wind field data, and the inverter was used to control the fan speed. After the wind speed was stable, the lift coefficient and drag coefficient were determined by the six-component strain balance. The data analysis system analyzed 1000 data sets and finally obtained the aerodynamic coefficients and air resistance of the drone. The wind tunnel test was conducted at the Chinese National Key Laboratory of Science and Technology on Aerodynamic Design and Research. During testing, the pitch angle of the UAV was set at 8°, and the wind speed was stabilized between 12 and 20 m/s using a frequency converter at 2 m/s intervals. The lift-drag ratio and air resistance values of the UAV sample were determined at different wind speeds and the simulated values of the CFD was compared with the measured values in the wind tunnel test to verify the accuracy of CFD simulation methods.

Numerical Simulation Results Verified
Change in the UAV aerodynamic operating conditions is essentially a change in the Reynolds number. According to the Reynolds number similarity theory and geometric similarity, the same Reynolds number produces similar aerodynamic characteristics. The scatter diagrams for the aerodynamic coefficient CFA and D of the UAV under different Reynolds numbers were drawn and compared with the wind tunnel data as shown in Figure 8. The wind tunnel test was conducted at the Chinese National Key Laboratory of Science and Technology on Aerodynamic Design and Research. During testing, the pitch angle of the UAV was set at 8 • , and the wind speed was stabilized between 12 and 20 m/s using a frequency converter at 2 m/s intervals. The lift-drag ratio and air resistance values of the UAV sample were determined at different wind speeds and the simulated values of the CFD was compared with the measured values in the wind tunnel test to verify the accuracy of CFD simulation methods.

Numerical Simulation Results Verified
Change in the UAV aerodynamic operating conditions is essentially a change in the Reynolds number. According to the Reynolds number similarity theory and geometric similarity, the same Reynolds number produces similar aerodynamic characteristics. The scatter diagrams for the aerodynamic coefficient C FA and D of the UAV under different Reynolds numbers were drawn and compared with the wind tunnel data as shown in Figure 8. m/s intervals. The lift-drag ratio and air resistance values of the UAV sample were determined at different wind speeds and the simulated values of the CFD was compared with the measured values in the wind tunnel test to verify the accuracy of CFD simulation methods.

Numerical Simulation Results Verified
Change in the UAV aerodynamic operating conditions is essentially a change in the Reynolds number. According to the Reynolds number similarity theory and geometric similarity, the same Reynolds number produces similar aerodynamic characteristics. The scatter diagrams for the aerodynamic coefficient CFA and D of the UAV under different Reynolds numbers were drawn and compared with the wind tunnel data as shown in Figure 8. The average error of the aerodynamic coefficient CFA was less than 5%, while the air resistance D less than 6%. Since the UAV model was constructed using 3D printing, the surface roughness The average error of the aerodynamic coefficient C FA was less than 5%, while the air resistance D less than 6%. Since the UAV model was constructed using 3D printing, the surface roughness affected the flow field distribution close to the wings. The model's air resistance was small, and the model's jitter, signal interference, and other factors in the test process affected the accuracy of the air resistance measurements. The errors in the test results were less than 6% within the allowable range, which indicates that the CFD numerical simulation method is reliable and can simulate the outflow field of tailsitter UAVs. Figure 9 shows the half-normal distribution of the nine structural parameters with the time coefficient and air resistance.

Feature Factor Extraction Results
The half-normal distribution chart shows that the specific factors of the time coefficient and air resistance are the root chord, wingtip chord, wingspan, and sweep angle. In the local sensitivity of air resistance, the root chord, wingtip chord, and wingspan were positively affected, while the sweep angle was negatively affected, and the sensitivity was Cr > b > C t > Λ w . In time coefficient terms, the wing root chord length and wingspan length were positively affected, while the sweep angle and wing chord length were negatively affected, and the sensitivity was b > Λ w > C r > C t .
The mathematical expression of air resistance and time coefficient is Figure 9 shows the half-normal distribution of the nine structural parameters with the time coefficient and air resistance.

Feature Factor Extraction Results
It can be seen from the mathematical expressions and Figure 9 that there are no interactions between the four influence factors of air resistance, however, there is an interaction between the wingspan and sweep angle, as well as between the root chord length and sweep angle for time coefficient. The aerodynamic coefficient of the drone is mainly determined by the shape and the incoming angle. Moreover, the two factors interact in the distribution of the fuselage flow field. The sweep angle is the main factor affecting the incoming angle, and the wing root and wingspan are the main factors of the shape. The interactions between the wingspan, sweep angle, and root chord length should be considered when calculating endurance time.  It can be seen from the mathematical expressions and Figure 9 that there are no interactions between the four influence factors of air resistance, however, there is an interaction between the wingspan and sweep angle, as well as between the root chord length and sweep angle for time coefficient. The aerodynamic coefficient of the drone is mainly determined by the shape and the incoming angle. Moreover, the two factors interact in the distribution of the fuselage flow field. The sweep angle is the main factor affecting the incoming angle, and the wing root and wingspan are the main factors of the shape. The interactions between the wingspan, sweep angle, and root chord length should be considered when calculating endurance time. Figure 10 shows the relationship between the wingspan and the time coefficient. Under the specified attack angle, the time coefficient increased with an increase of wingspan b. The endurance time was determined by the mass coefficient and aerodynamic coefficient. The mass coefficient was determined by the projected area, which was positively correlated with wingspan b. The aerodynamic coefficient was related to the lift-drag ratio, which was determined by the aspect ratio. Therefore, the aspect ratio had a linear relationship with the wingspan length. Figure 10 shows the relationship between the wingspan and the time coefficient. Under the specified attack angle, the time coefficient increased with an increase of wingspan b. The endurance time was determined by the mass coefficient and aerodynamic coefficient. The mass coefficient was determined by the projected area, which was positively correlated with wingspan b. The aerodynamic coefficient was related to the lift-drag ratio, which was determined by the aspect ratio. Therefore, the aspect ratio had a linear relationship with the wingspan length.     Figure 11 shows the variation of the time coefficient with the wing root chord length under different attack angles. At the same angle of attack, the time coefficient increased with the wing root chord length. When the attack angle was 8 • , the time coefficient with a wing root chord length of 500 mm was 3% higher than that with a wing root chord length of 5 mm. specified attack angle, the time coefficient increased with an increase of wingspan b. The endurance time was determined by the mass coefficient and aerodynamic coefficient. The mass coefficient was determined by the projected area, which was positively correlated with wingspan b. The aerodynamic coefficient was related to the lift-drag ratio, which was determined by the aspect ratio. Therefore, the aspect ratio had a linear relationship with the wingspan length.  Figure 11 shows the variation of the time coefficient with the wing root chord length under different attack angles. At the same angle of attack, the time coefficient increased with the wing root chord length. When the attack angle was 8°, the time coefficient with a wing root chord length of 500 mm was 3% higher than that with a wing root chord length of 5 mm.   Figure 12 shows the variation of the time coefficient with the wingtip chord length under different attack angles. The results illustrate that the time coefficient first decreased and then increased as the wingtip chord length increased. Under the specified wingtip chord length, there existed a minimum time coefficient value. Moreover, the time coefficients with wingtip chord lengths of 180 mm and 300 mm were 2.6% and 0.8% higher than those with a wingtip chord length of 270 mm, respectively. As the size of the wingtip increased, the turbulence of the wing tip gradually weakened, while the projected area of the wing increased. Accordingly, the endurance time first exhibited a decreasing, and then an increasing, trend with an increase of the wingtip.

Effects of the Structural Parameters on the Time Coefficient
increased as the wingtip chord length increased. Under the specified wingtip chord length, there existed a minimum time coefficient value. Moreover, the time coefficients with wingtip chord lengths of 180 mm and 300 mm were 2.6% and 0.8% higher than those with a wingtip chord length of 270 mm, respectively. As the size of the wingtip increased, the turbulence of the wing tip gradually weakened, while the projected area of the wing increased. Accordingly, the endurance time first exhibited a decreasing, and then an increasing, trend with an increase of the wingtip.  Figure 13 shows the variation of the time coefficient with the wing sweep angle under different attack angles. When the sweep angle increased from 15° to 35°, the time coefficient of the UAV linearly decreased. The time coefficient with a sweep angle of 35° was 4% lower than that with a sweep angle of 15°. This result can be explained by the aerodynamic coefficient CFA when the sweep angle changed. As the sweep angle increased, the lift coefficient CL first increased and then decreased. The lift-drag ratio K first decreased and then increased, and the combination of the two exhibited a linear downward trend. Under the specified attack angle, the aerodynamic coefficient CFA increased with the wing sweep angle Λw, and the wing sweep angle had less of an impact on the mass coefficient GFA. Thus, the time coefficient increased with the wing sweep angle in a specified geometrical structure.  When the angle of attack changed from 4° to 12°, the time coefficient increased in the range of 4-8° and decreased in the range of 10-12°. The endurance time significantly declined with an attack angle of 6°. As the angle of attack increased, the laminar flow effect on the wing surface gradually weakened, the turbulence effect intensified, the aerodynamic coefficient efficiency of the UAV decreased, and the change rate of the endurance time decreased. Figure 14 shows the relationship between the wingspan and air resistance. Under the specified attack angle, the air resistance increased with wingspan b. Moreover, the rate of the air resistance growth became faster as the wingspan increased in a specified geometrical structure, which can be When the angle of attack changed from 4 • to 12 • , the time coefficient increased in the range of 4-8 • and decreased in the range of 10-12 • . The endurance time significantly declined with an attack angle of 6 • . As the angle of attack increased, the laminar flow effect on the wing surface gradually weakened, the turbulence effect intensified, the aerodynamic coefficient efficiency of the UAV decreased, and the change rate of the endurance time decreased. Figure 14 shows the relationship between the wingspan and air resistance. Under the specified attack angle, the air resistance increased with wingspan b. Moreover, the rate of the air resistance growth became faster as the wingspan increased in a specified geometrical structure, which can be explained by the differential pressure between the upper and lower wing surfaces. Owing to the effect of the air pressure difference on the wingtip, their induced resistance was present from the bottom-up. As the wingspan increased, the induced air resistance increased, and thus the air resistance increased.

Effects of the Configuration Parameters on Air Resistance Performance
4-8° and decreased in the range of 10-12°. The endurance time significantly declined with an attack angle of 6°. As the angle of attack increased, the laminar flow effect on the wing surface gradually weakened, the turbulence effect intensified, the aerodynamic coefficient efficiency of the UAV decreased, and the change rate of the endurance time decreased. Figure 14 shows the relationship between the wingspan and air resistance. Under the specified attack angle, the air resistance increased with wingspan b. Moreover, the rate of the air resistance growth became faster as the wingspan increased in a specified geometrical structure, which can be explained by the differential pressure between the upper and lower wing surfaces. Owing to the effect of the air pressure difference on the wingtip, their induced resistance was present from the bottom-up. As the wingspan increased, the induced air resistance increased, and thus the air resistance increased. The air resistance variation with the root chord length under different wingtip chord lengths is shown in Figure 15. As can be seen, the air resistance increased with the root chord length under the specified length of the wingtip chord. When the wingtip chord length changed from 180 to 240 mm, the air resistance increased under the specified length of the wing root chord length. Additionally, the growth rate of the air resistance increased at 180-240 mm and gradually decreased at 270-300 mm. The air resistance variation with the root chord length under different wingtip chord lengths is shown in Figure 15. As can be seen, the air resistance increased with the root chord length under the specified length of the wingtip chord. When the wingtip chord length changed from 180 to 240 mm, the air resistance increased under the specified length of the wing root chord length. Additionally, the growth rate of the air resistance increased at 180-240 mm and gradually decreased at 270-300 mm. This result can be explained by the wing root effect on the swept-wing. When the air flowed through the swept-wing, the streamline was slightly skewed toward the wingtip, which resulted in an s-shaped streamline distribution. On the upper surface of the wing root, the streamline exhibited an expanding trend at the front region, and the streamline exhibited a shrinking trend in the back This result can be explained by the wing root effect on the swept-wing. When the air flowed through the swept-wing, the streamline was slightly skewed toward the wingtip, which resulted in an s-shaped streamline distribution. On the upper surface of the wing root, the streamline exhibited an expanding trend at the front region, and the streamline exhibited a shrinking trend in the back area. The offset of the streamline increased with the wing root chord length and wingtip chord length, which increased the air resistance on the wing. The root-tip ratio decreased with an increase of the wingtip chord length, which reduced the induced resistance of the wingtip. The air resistance gradient decreased as the induced resistance decreased. In short, the air resistance increased with the wing root chord length and wingtip chord length.

Multi-Objective Optimization
Multi-objective optimization, including the maximum t FA and minimum D, was applied to obtain a set of optimal solutions (as shown in Figure 16). The initial population was constructed using the screening method, and the time coefficient and air resistance were calculated using the response surface model. Table 3 presents the three optimized configurations. Compared with the original prototype, the time coefficient t FA of optimal configuration 1 was enhanced by 19.5%, while the air resistance D was reduced by 34.78%. For all three optimal configurations, the time coefficient t FA increased by an average of 10.5%, while the air resistance D was reduced by an average of 40.58%. This study selected optimal configuration 1 as the solution with the best compromise. The parameters of the final tailsitter UAV structure are a wing root chord length of 315 mm, a wingtip chord length of 182 mm, a wingspan length of 1198 mm, and sweep angle of 16 • .  In this study, three optimal configurations were simulated and verified. The CFD simulation value was used as the actual value for comparison with the response surface prediction value, as presented in Table 3. The relative error of D and tFA increased with the sweep angle. For all three verification points, the relative error of the time coefficient tFA was less than 3%, while the relative error of the air resistance D was less than 5%.
Based on previous research and increasing demand for agricultural remote sensing, a VTOL UAV was designed with vertical take-off and landing characteristics, as well as dynamic cruise characteristics. These characteristics solve the low efficiency and short-range of failures for multirotor UAVs, the lengthy fixed-wing UAV take-off preparation, auxiliary equipment requirements, and various landing problems. Ang et al. designed a foldable tailsitter U-lion in 2014 [12]. The wing of this U-lion could be retracted through a link mechanism, which improved the maneuverability and portability of the drone. Aksugur and Inalhan designed a ducted tailsitter drone and adjusted the attitude in the vertical take-off and landing state by using a ducted steering gear [19]. In contrast, our tailsitter UAV used a dual-wing symmetrical layout, and there was no additional auxiliary mechanism for attitude adjustment. This UAV has the characteristics of a compact structure, light weight, and smooth operation.
The predecessors of this system mainly used CFD calculations to compare the lift-drag ratio and  In this study, three optimal configurations were simulated and verified. The CFD simulation value was used as the actual value for comparison with the response surface prediction value, as presented in Table 3. The relative error of D and t FA increased with the sweep angle. For all three verification points, the relative error of the time coefficient t FA was less than 3%, while the relative error of the air resistance D was less than 5%.
Based on previous research and increasing demand for agricultural remote sensing, a VTOL UAV was designed with vertical take-off and landing characteristics, as well as dynamic cruise characteristics. These characteristics solve the low efficiency and short-range of failures for multi-rotor UAVs, the lengthy fixed-wing UAV take-off preparation, auxiliary equipment requirements, and various landing problems. Ang et al. designed a foldable tailsitter U-lion in 2014 [12]. The wing of this U-lion could be retracted through a link mechanism, which improved the maneuverability and portability of the drone. Aksugur and Inalhan designed a ducted tailsitter drone and adjusted the attitude in the vertical take-off and landing state by using a ducted steering gear [19]. In contrast, our tailsitter UAV used a dual-wing symmetrical layout, and there was no additional auxiliary mechanism for attitude adjustment. This UAV has the characteristics of a compact structure, light weight, and smooth operation.
The predecessors of this system mainly used CFD calculations to compare the lift-drag ratio and air resistance of each sample in the discrete sample points to find the optimal sample point. Cetinsoy et al. used CFD simulations to calculate the winglet turbulence under three airfoils and two different angles of attack. The optimal winglet size was determined by the wingtip vortex and lift [16]. Rao et al. used CFD simulations to analyze the UAV's drag coefficient and the lift-drag ratio at ten wind speeds and nine angles and determined the aerodynamic characteristics of the UAV in transition mode [14]. In this study, a calculation model of endurance time was established, and a response surface model of endurance time and air resistance was constructed. A multi-objective genetic algorithm was used to determine the optimal structural parameters. Unlike previous studies, this study adopts battery life as a comprehensive performance index to evaluate the versatility of tailsitter UAVs. The response surface model has sufficient continuity to achieve the optimization of global structural parameters. Compared with the original prototype, the final solution structure with a good compromise was determined by using the multi-objective genetic algorithm to improve the time coefficient by 19.5% and reduce the air resistance by 34.78%. These results show that a multi-objective genetic algorithm based on the genetic aggregation response surface can be used to optimize tailsitter drones.

Conclusions
This study designed an agricultural VTOL tailsitter UAV with symmetrical winglets and wings. The time coefficient and resistance were calculated under different structural parameter combinations using the CFD numerical simulation method. Additionally, a response surface model of the time coefficient and air resistance with genetic aggregation was developed, and the structural parameters were optimized using a multi-objective genetic algorithm.
(1) The root chord, wingtip chord, wingspan, and sweep angle were the specific factors that influenced the time coefficient and air resistance. In terms of air resistance, the root chord, wingtip chord, and wingspan were positively affected, while the sweep angle was negatively affected, and the sensitivity was Cr > b > C t > Λ w . In time coefficient terms, the wing root chord length and wingspan length were positively affected, while the sweep angle and wing chord length were negatively affected, and the sensitivity was b > Λ w > C r > C t .
(2) The CFD simulation was verified by conducting a wind tunnel test, and the results revealed that the aerodynamic coefficient error was less than 5%, while the resistance error was less than 6%. This result indicates that the CFD numerical simulation method can be used to simulate the external flow field of a tailsitter UAV.
(3) The time coefficient increased with the wing root chord length and wingspan, decreased with an increase of the sweep angle, and first decreased and then increased as the wingtip chord length increased. The air resistance increased as the wing root chord length, wingtip chord length, and wingspan increased, and first decreased and then increased with an increase of the sweep angle.
(4) The solution parameters for the tailsitter UAV with the best compromise were a wing root chord length of 315 mm, a wingtip chord length of 182 mm, a wingspan length of 1198 mm, and a sweep angle of 16 • . Compared with the original shape, the time coefficient t FA was enhanced by 19.5% while the air resistance D declined by 34.78%.
Author Contributions: W.H. was the supervisor of the whole system and experiment. X.Y. conceived and designed this system and experiment and drafted the manuscript; X.Y. and W.L. discussed and drafted the manuscript; W.L., and G.L. conducted the collection of experimental data; G.L helped during laboratory experiments; Q.M. gave advice on the manuscript modification. All authors have read and agreed to the published version of the manuscript.