Modeling, System Measurements and Controller Investigation of a Small Battery-Powered Fixed-Wing UAV

: In this paper, a complete set of nonlinear modeling and controller design process for a small electric ﬁxed-wing unmanned aerial vehicle (UAV) is presented. The nonlinear mathematical model and aerodynamic model of the small ﬁxed-wing UAV are derived. The computational ﬂuid dynamics (CFD) method was used to obtain the aerodynamic coefﬁcients of the UAV, and the models of propulsion system components were established through experiments. Since the linearized and decoupled model of the ﬁxed-wing UAV has a large error, a nonlinear model is established based on Simulink, which is utilized to design and verify the control algorithms. Based on the established nonlinear model, a stability controller, path following controller and path management controller of the aircraft are set up. The results indicate that system parameters of the aircraft can be quickly acquired and an efﬁcient and practical model can be established by the methods. In addition, the controller designed and applied in this paper has good performance and small steady-state error, which can meet the basic ﬂight mission requirements, including stability of ﬂight attitude, path following and switching of different waypoints. These modeling and control methods can also be employed in other small battery-powered ﬁxed-wing UAV projects.


Introduction
UAVs have attracted increasing research attention in the last decades and gained many applications in both civilian and military operations with improvements in avionics devices. The UAVs have been employed for agricultural monitoring, reconnaissance, photography and remote sensing [1][2][3][4]. Due to the convenience and high flexibility of small UAVs, the demand for their development and application continues to increase.
Small UAVs refer to UAVs weighing less than 5 kg, including small fixed-wing aircraft, helicopters and multicopters [5]. UAVs are mostly powered by battery, and a few UAVs are powered by fuel [6]. Compared with rotary-wing UAVs, fixed-wing UAVs have longer endurance and better flight efficiency [7]. The aerodynamic layout of a fixed-wing aircraft can affect its lift-to-drag ratio, which in turn affects endurance [8]. Flying wing UAV is a non-canard and tailless aircraft with integrated fuselage and wing [9]. The unconventional aircraft has lower aerodynamic drag compared with the conventional configuration, while the special aerodynamic shape will add more challenges to the stability of the aircraft [10,11]. In order to develop a small electric UAV with flying wing configuration, it is vital to measure the system parameters, establish an accurate mathematical model and design a reliable controller.
Predicting aerodynamic forces and moments of the UAV is very difficult due to the external disturbances [12] and there are several approaches utilized to assess the aerodynamic derivatives and characteristics of the fixed-wing UAV. Wind tunnel testing [13], DATCOM software [14], empirical relations [15], flight testing [16] and CFD simulation [17] are employed to acquire aerodynamic characteristics of the UAVs in different studies. Wind tunnel testing is one of the most important tools because it can accurately control the environmental conditions and the measurement results are accurate; however, its equipment is expensive and the tests are complex. Flight testing is easily affected by external disturbances. The results obtained by DATCOM software and empirical relations usually have more errors than other methods. CFD method may not obtain accurate results, and the simulation results will be affected by parameters such as grid quality and boundary settings; however, it is very suitable for the preliminary design of aircraft.
For a small electric fixed-wing UAV, in addition to the aerodynamic coefficients, it is also important to obtain the parameters of the propulsion system. The thrust and torque data of a small UAV's propulsion system can be obtained by wind tunnel tests [18]. ZHAO proposed a novel airborne test method suitable for fixed-wing aircraft, which can replace the high-cost wind tunnel test [19]. Piljek proposed a method for characterization of the multirotor UAV propulsion system, and conducted experimental measurements to validate the method [20]. Virginio proposed a low-cost Thrust Benchmarking System for an electric UAV propulsion system, and three kinds of propellers were modeled by the instrument [21]; however, the propulsion system and actuator models required for UAV modeling have not been well summarized.
Since the fixed-wing UAV is a nonlinear and strongly coupled system, linearization of equations of motion using small perturbation theory is often used in literature. UAV equations of motion can be separated into longitudinal and lateral directional dynamics and decoupled between them. The perturbation with little difference from the reference motion is called the small perturbation motion. Such difference is related to the specific situation and has no absolute value. This assumption is no longer applicable to a UAV that needs to complete multiple flight missions, such as hovering and landing. A mini-UAV designed for agricultural application was separated into longitudinal and lateral state-space models, and five different controllers were designed for the UAV lateral dynamics [22]. An MPC strategy based on the lateral and longitudinal linear models was proposed to examine a fixed-wing UAV as it undergoes five flight scenarios [23]. Melkou adopted an uncertain linear model decoupled into longitudinal and lateral directional modes to represent the aircraft and added a PID outer loop to improve the accuracy [24]. The method of linearizing and decoupling a fixed-wing UAV is very common, and it is convenient for researchers to design controllers; however, the linearized model will produce larger errors in the process of simulation and flight test, and the use of nonlinear models can better predict the aerodynamics of the UAV [25].
The previous studies indicate the importance of small electric UAVs and how difficult it is to design and develop them. For aircraft conceptual design, there are some methodologies [11,26]. A methodology for sizing the electric propulsion subsystems of UAVs was developed by Gokcin [27]. Avanzini proposed a methodology to obtain best-endurance and best-range conditions by optimizing the aircraft size [28]; however, there are few studies on flying wing configuration UAVs with only two elevons, and most of the design and research on fixed-wing UAVs are based on linearized models. This work focuses on the mathematical modeling and controller design of the unconventional UAV configuration. The methodologies to acquire the parameters of the UAV system, including parameters of the aerodynamic structure and electric equipment, are well summarized, and the nonlinear model of the small electric flying wing UAV is quickly established. The stability controller of the UAV is established based on PID algorithm by combining the characteristics of the small electric flying wing UAV. The shortcomings of some previous path-following methods are summarized in this paper, and a comprehensive and optimized path-following algorithm is proposed. As an initial stage in this work, the nonlinear mathematical model, Machines 2021, 9,333 3 of 25 including the 6-DoF equations of motion and aerodynamic model, is derived. Aerodynamic derivatives are obtained using the CFD software Fluent and propulsion system parameters are acquired by corresponding tests. A complete, efficient and practical simulation model was quickly established to verify the effectiveness of the algorithm. In the final stage of the work, a stability controller, path following controller and path management controller were designed and verified in turn. A set of systematic development procedures of a small electric fixed-wing UAV is presented in this paper, which provides a reference for the design and development of other small electric fixed-wing UAVs. MATLAB/Simulink was utilized to establish the mathematical model, solve the nonlinear equations of motion numerically and validate the control method in our research. The remaining part of the paper is organized as follows. In Section 2, the nonlinear mathematical modeling of the fixed-wing UAV is indicated. In Section 3, we detail the parameters of the UAV system obtained by experiments and simulations. The controllers, including the stability controller and the upper controllers, are explained in Section 4. Discussion on system parameters measurements and control methods is presented in Section 5. Finally, the summary and conclusions are presented in Section 6.

Mathematical Model of the Fixed Wing UAV
The geometry model and aerodynamic shape data of the UAV investigated in our study are taken from our previous work [29], which is shown in Figure 1. The geometrical values are given in Table 1. In Table 1, AR is the aspect ratio, S is the area of the wing, b is the wingspan, c is the mean aerodynamic chord of the wing and m a is the mass of the UAV. The small fixed-wing UAV we designed only has two control surfaces, as known as elevon, and they have the same effects as both ailerons and elevators, and it has a brushless DC electric motor at the end side of the UAV to provide thrust. Based on this aerodynamic shape, we made a prototype for accurate parameter identification, but there are some differences with theoretical parameters. For example, the theoretical gross weight is 1.5 kg, but the actual mass is 0.9 kg.
Machines 2021, 9, x FOR PEER REVIEW 3 of 28 controller of the UAV is established based on PID algorithm by combining the characteristics of the small electric flying wing UAV. The shortcomings of some previous pathfollowing methods are summarized in this paper, and a comprehensive and optimized path-following algorithm is proposed. As an initial stage in this work, the nonlinear mathematical model, including the 6-DoF equations of motion and aerodynamic model, is derived. Aerodynamic derivatives are obtained using the CFD software Fluent and propulsion system parameters are acquired by corresponding tests. A complete, efficient and practical simulation model was quickly established to verify the effectiveness of the algorithm. In the final stage of the work, a stability controller, path following controller and path management controller were designed and verified in turn. A set of systematic development procedures of a small electric fixed-wing UAV is presented in this paper, which provides a reference for the design and development of other small electric fixedwing UAVs.MATLAB/Simulink was utilized to establish the mathematical model, solve the nonlinear equations of motion numerically and validate the control method in our research. The remaining part of the paper is organized as follows. In Section 2, the nonlinear mathematical modeling of the fixed-wing UAV is indicated. In Section 3, we detail the parameters of the UAV system obtained by experiments and simulations. The controllers, including the stability controller and the upper controllers, are explained in Section 4. Discussion on system parameters measurements and control methods is presented in Section 5. Finally, the summary and conclusions are presented in Section 6.

Mathematical Model of the Fixed Wing UAV
The geometry model and aerodynamic shape data of the UAV investigated in our study are taken from our previous work [29], which is shown in Figure 1. The geometrical values are given in Table 1. In Table 1, AR is the aspect ratio, S is the area of the wing, b is the wingspan, c is the mean aerodynamic chord of the wing and ma is the mass of the UAV. The small fixed-wing UAV we designed only has two control surfaces, as known as elevon, and they have the same effects as both ailerons and elevators, and it has a brushless DC electric motor at the end side of the UAV to provide thrust. Based on this aerodynamic shape, we made a prototype for accurate parameter identification, but there are some differences with theoretical parameters. For example, the theoretical gross weight is 1.5 kg, but the actual mass is 0.9 kg.

6-DOF Equations of Motion
In order to design and validate the UAV control system, it is very important to describe the details of the UAV's dynamics and kinematics. The six degrees of freedom nonlinear model of the UAV consists of force equations, moment equations and kinematic equations for position and navigation. These dynamic equations for the UAV are given as follows [30].
are velocities of the UAV in the x, y and z body axis, while ω = [p, q, r] are the angular velocities about the same axis, respectively. F b is the sum of all the forces acting on the UAV in the body axis. Moment equation: .
In Equation (2), I b is the 3 × 3 inertia matrix of the UAV and M b is the sum of external moments of the body axis.
where p e represents the position of the UAV in the Earth coordinate system; R e b is transformation matrix from body axis to Earth axis. Navigation equation: where Θ = [ϕ, θ, ψ] are Euler angles, W is transformation matrix from angular velocities to Euler rates. When the initial conditions of the UAV are determined, the attitude and position information of the UAV can be obtained by inputting the force and the moment into the system. For the UAV, the force can be split into thrust F t , gravity F g and aerodynamics forces F a , as seen in Equation (5). Since the gyroscopic moment generated by the rotation of the propeller and the motor is far less than other moments, it is approximately zero here. Similarly, the moment on the UAV, which is assumed to act around the center of gravity, can be decomposed into the counter torque M t exerted by the propeller [31] and the aerodynamic moments M a , as seen in Equation (6).
The thrust and the counter torque generated by the rotating propeller were obtained through experiments, detailed in subsequent sections. The aerodynamic moments are modeled in the body frame, but the aerodynamic forces are modeled in the wind frame, which need to transform to the body frame, as seen in Equation (7).
where F x , F y and F z indicate the aerodynamic forces along x, y, z axis of the body frame; D, Y and L denote the aerodynamic drag force, side force and lift force acting along the negative x, positive y and negative z wind axis, respectively, R b w is transformation matrix from wind axis to body axis. More details about the modeling of the UAV are shown in (A1)-(A3).

Aerodynamic Model
Since UAVs have different aerodynamic characteristics under different conditions, modeling the aerodynamic forces and moment is one of the most vital and complicated parts of designing a fixed-wing UAV. A commonly used model for aerodynamic forces and moments acting on a fixed-wing UAV is shown in Equations (8) and (9). where ρ is the air density, V a is the airspeed of the UAV. C D , C Y and C L are the nondimensional coefficients of drag, side and lift force, respectively. l, m and n denote the aerodynamic moments along x, y and z-axis of the body frame and can be modeled in terms of the dimensionless functions C l , C m and C n , respectively. All the aerodynamic forces are modeled in the wind frame. The aerodynamic coefficients of forces and moments depend on many variables with complicated nonlinear correlation, and it is a mathematically convenient way to build up the coefficients as a sum of contributing components. There are several kinds of expressions for the force and moment coefficients, and a simple but reasonable expression [32] is employed in this paper, as seen in Equations (10) and (11).
where α is the angle of attack and β is the sideslip angle, δ e , δ a , δ r denote the elevator, aileron and rudder deflection of the UAV, respectively. Due to the particular structure of the flying wing UAV, δ r is set to 0, δ e and δ a can be decoupled into expressions of left and right elevon deflections by Equation (12). All the other aerodynamic coefficients are presented in the following section and their values are given in Table 2.  δ er and δ el represent the angular deflections of the right and left elevon, respectivelyassuming that the aerodynamic effects of the left and right elevons are the same and there is no coupling between them. In Section 3, the aerodynamic coefficient values affected by δ e and δ a are obtained, which are twice the values affected by a single elevon. So aerodynamic forces and moments can be calculated with coefficients of δ er and δ el .

System Parameters Measurements
For establishing an accurate nonlinear simulation model of the UAV, some system parameters need to be measured to improve the UAV model, including aerodynamic derivatives, propulsion system parameters and basic parameters of the UAV. All the simulation and test data are fitted through the curve fitting toolbox in MATLAB, and the fitting results are presented in the following subsections.

Aerodynamic Derivatives Measurements
All the aerodynamic derivatives can be divided into static stability derivatives (influenced by α and β), damping derivatives (influenced by p, q, r) and control derivatives (influenced by δ e and δ a ) and their values determine the static stability, dynamic stability and control effects of the UAV, respectively. CFD method is investigated and used in this paper to obtain the aerodynamic characteristics of the UAV numerically through comprehensive consideration of accuracy and convenience.
The stability derivatives and the control derivatives can be obtained under steady-state simulation conditions, while the damping derivatives need dynamic simulation conditions. In this work, the stability derivatives and the control derivatives are collectively referred to as static derivatives, and the damping derivatives are referred to as dynamic derivatives. The aerodynamic analysis of the UAV and the calculation of all the aerodynamic derivatives were implemented with the Fluent software. The finite volume method is employed to solve the N-S equations to acquire aerodynamics [33], and the SST k-w turbulence model is utilized as the fluid model.

Static Derivatives
In the simulation of calculating the static derivatives, the air velocity is set to 15 m/s. The tests focus on the quality of the model around α = β = δ e = δ a = 0. The corresponding changes of aerodynamic coefficients can be obtained by changing the values of the parameters. The Reynold number for the UAV corresponding to cruise speed is about 400,000 when the air temperature is 288.15 K and the air density is 1.225 kg/m 3 .
In Figure 2, the lift force coefficient, drag force coefficient and pitch moment coefficient are presented as a function of the angle of attack. When the angle of attack is −10 • to 12 • , the lift force coefficient and pitch moment coefficient increase linearly. The airflow begins to separate and the nonlinear phenomena of the lift force and the pitch moment occur when the attack angle exceeds 12 • . The UAV is close to the stall at 20 • angle of attack. In this study, in order to facilitate modeling, they are approximated as linear models. Since the UAV designed in this paper rarely has a high attack angle during flight missions, the error caused by linearization is tolerable; however, the linearized model may not be applicable to complex missions that require a high attack angle. It should be noted that only the data with an attack angle less than 16 • is used for the fitting of the lift and pitch moment coefficients. This is to make the model more accurate with a small attack angle. As shown in Figure 2b, the drag force coefficient curve can be well fitted with a quadratic curve. In Figure 2, the lift force coefficient, drag force coefficient and pitch moment coefficient are presented as a function of the angle of attack. When the angle of attack is −10° to 12°, the lift force coefficient and pitch moment coefficient increase linearly. The airflow begins to separate and the nonlinear phenomena of the lift force and the pitch moment occur when the attack angle exceeds 12°. The UAV is close to the stall at 20° angle of attack. In this study, in order to facilitate modeling, they are approximated as linear models. Since the UAV designed in this paper rarely has a high attack angle during flight missions, the error caused by linearization is tolerable; however, the linearized model may not be applicable to complex missions that require a high attack angle. It should be noted that only the data with an attack angle less than 16° is used for the fitting of the lift and pitch moment coefficients. This is to make the model more accurate with a small attack angle. As shown in Figure 2b, the drag force coefficient curve can be well fitted with a quadratic curve. The sideslip angle is one of the parameters that have the greatest impact on the UAV due to external wind interference, and it can affect a number of aerodynamic parameters of the fixed-wing UAV, which we can see in Figure 3. The drag force, side force, roll moment and yaw moment coefficients are plotted versus the sideslip angle in Figure 3. In the range of sideslip angle of −20° to 20°, the drag force coefficient uses quadratic fitting, and the other coefficients use linear fitting. The sideslip angle is one of the parameters that have the greatest impact on the UAV due to external wind interference, and it can affect a number of aerodynamic parameters of the fixed-wing UAV, which we can see in Figure 3. The drag force, side force, roll moment The sideslip angle is one of the parameters that have the greatest impact on the UA due to external wind interference, and it can affect a number of aerodynamic paramete of the fixed-wing UAV, which we can see in Figure 3. The drag force, side force, roll m ment and yaw moment coefficients are plotted versus the sideslip angle in Figure 3. In t range of sideslip angle of −20° to 20°, the drag force coefficient uses quadratic fitting, a the other coefficients use linear fitting. In Figures 4 and 5, the aerodynamic force and moment coefficients are plotted functions of the elevator deflection and the aileron deflection. The values of the contr derivatives are equal to the slopes of these curves. As shown in Figure 4, the lift coefficie and pitching moment coefficient increase linearly with the increase in the elevator ang Due to the asymmetric shape of the elevator's upper deflection and lower deflection, t drag coefficient is not completely symmetric on both sides of δe = 0°. For the convenien of modeling, we only keep the quadratic term of the curve, which will result in a ma mum error of 12% when δe = 12°. Such error introduced by the hypothesis of keeping t drag symmetric is within the acceptable range, so we ignore it here. Due to the symm rical shape of the UAV, it is only necessary to calculate the aileron angle of 0° to 20° reduce the simulation time. As shown in Figure 5, when the aileron angle is below 1 with the increase in the aileron angle, the side force coefficient and roll moment coefficie increase approximately linearly, and the yaw moment coefficient decreases linearly. Sin the aileron has little effect on the side force coefficient and the yaw moment coefficie they are approximated as a linear relationship to simplify the model. For the same reas as the attack angle, the fitting curve does not include the data of side force and yaw m ment coefficients when δa is 16° and 20°  Figure 4, the lift coefficient and pitching moment coefficient increase linearly with the increase in the elevator angle. Due to the asymmetric shape of the elevator's upper deflection and lower deflection, the drag coefficient is not completely symmetric on both sides of δ e = 0 • . For the convenience of modeling, we only keep the quadratic term of the curve, which will result in a maximum error of 12% when δ e = 12 • . Such error introduced by the hypothesis of keeping the drag symmetric is within the acceptable range, so we ignore it here. Due to the symmetrical shape of the UAV, it is only necessary to calculate the aileron angle of 0 • to 20 • to reduce the simulation time. As shown in Figure 5, when the aileron angle is below 16 • , with the increase in the aileron angle, the side force coefficient and roll moment coefficient increase approximately linearly, and the yaw moment coefficient decreases linearly. Since the aileron has little effect on the side force coefficient and the yaw moment coefficient, they are approximated as a linear relationship to simplify the model. For the same reason as the attack angle, the fitting curve does not include the data of side force and yaw moment coefficients when δ a is 16 • and 20 • . reduce the simulation time. As shown in Figure 5, when the aileron angle is below 16°, with the increase in the aileron angle, the side force coefficient and roll moment coefficient increase approximately linearly, and the yaw moment coefficient decreases linearly. Since the aileron has little effect on the side force coefficient and the yaw moment coefficient, they are approximated as a linear relationship to simplify the model. For the same reason as the attack angle, the fitting curve does not include the data of side force and yaw moment coefficients when δa is 16° and 20°

Dynamic Derivatives
In addition to the settings of calculating static coefficients, the rotating reference frame is adopted to gain the dynamic derivatives by simulation of the steady movement [34]. The unsteady flow field required for aircraft simulation can be converted into a steady flow field by the rotating reference frame method. The SRF (Single Reference Frame) model is utilized to simulate the rotation with the constant angular velocity of the UAV about the x, y, z-axis of the body frame in Fluent software, respectively. Boundary settings and meshes of the far-field are shown in Figure 6. An unstructured mesh is used to save computing time. The first height of the boundary layer mesh is 2 × 10 −5 m and y+ is 1. The growth ratio is 1.2 and 10 boundary layers are generated. The total number of 3D computational grid is about 4 × 10 6 . As shown in Figure 6b, velocity-inlet with magnitude and direction is selected in the far-field, which can ensure the airflow flows vertically to the aircraft when the fluid domain rotates. The unit of angular velocity is expressed in terms of rad/s in this subsection. In Figure 7, the lift force and pitch moment coefficients are presented as functions of the pitch rate for various attack angles. Since the fluid domain needs to simulate the steady climb movement of the UAV when obtaining aerodynamic coefficients related to pitch

Dynamic Derivatives
In addition to the settings of calculating static coefficients, the rotating reference frame is adopted to gain the dynamic derivatives by simulation of the steady movement [34]. The unsteady flow field required for aircraft simulation can be converted into a steady flow field by the rotating reference frame method. The SRF (Single Reference Frame) model is utilized to simulate the rotation with the constant angular velocity of the UAV about the x, y, z-axis of the body frame in Fluent software, respectively. Boundary settings and meshes of the far-field are shown in Figure 6. An unstructured mesh is used to save computing time. The first height of the boundary layer mesh is 2 × 10 −5 m and y+ is 1. The growth ratio is 1.2 and 10 boundary layers are generated. The total number of 3D computational grid is about 4 × 10 6 . As shown in Figure 6b, velocity-inlet with magnitude and direction is selected in the far-field, which can ensure the airflow flows vertically to the aircraft when the fluid domain rotates. The unit of angular velocity is expressed in terms of rad/s in this subsection.

Dynamic Derivatives
In addition to the settings of calculating static coefficients, the rotating reference frame is adopted to gain the dynamic derivatives by simulation of the steady movement [34]. The unsteady flow field required for aircraft simulation can be converted into a steady flow field by the rotating reference frame method. The SRF (Single Reference Frame) model is utilized to simulate the rotation with the constant angular velocity of the UAV about the x, y, z-axis of the body frame in Fluent software, respectively. Boundary settings and meshes of the far-field are shown in Figure 6. An unstructured mesh is used to save computing time. The first height of the boundary layer mesh is 2 × 10 −5 m and y+ is 1. The growth ratio is 1.2 and 10 boundary layers are generated. The total number of 3D computational grid is about 4 × 10 6 . As shown in Figure 6b, velocity-inlet with magnitude and direction is selected in the far-field, which can ensure the airflow flows vertically to the aircraft when the fluid domain rotates. The unit of angular velocity is expressed in terms of rad/s in this subsection. In Figure 7, the lift force and pitch moment coefficients are presented as functions of the pitch rate for various attack angles. Since the fluid domain needs to simulate the steady climb movement of the UAV when obtaining aerodynamic coefficients related to pitch rate, only positive rate values are considered. The lift force coefficient and pitch moment coefficient increase linearly, with pitch rate increasing at different attack angles. In this In Figure 7, the lift force and pitch moment coefficients are presented as functions of the pitch rate for various attack angles. Since the fluid domain needs to simulate the steady climb movement of the UAV when obtaining aerodynamic coefficients related to pitch rate, only positive rate values are considered. The lift force coefficient and pitch moment coefficient increase linearly, with pitch rate increasing at different attack angles. In this study, because the slope of the drag force coefficient curve is approximately zero, the influence of the pitch rate on the drag force is ignored.
Machines 2021, 9, x FOR PEER REVIEW 10 study, because the slope of the drag force coefficient curve is approximately zero, th fluence of the pitch rate on the drag force is ignored. The side force and roll moment coefficients versus the roll rate are plotted in Fi 8. The side force and roll moment coefficients decrease linearly with the increase in roll rate. In this study, the roll rate has little effect on the yaw moment coefficient, so ignored here. The side force, roll moment and yaw moment coefficients versus yaw rate are sented in Figure 9. With the increase in yaw rate, the side force coefficient and yaw ment coefficient decrease linearly, while the roll moment coefficient increases linearl The side force and roll moment coefficients versus the roll rate are plotted in Figure 8. The side force and roll moment coefficients decrease linearly with the increase in the roll rate. In this study, the roll rate has little effect on the yaw moment coefficient, so it is ignored here.
Machines 2021, 9, x FOR PEER REVIEW study, because the slope of the drag force coefficient curve is approximately zero, t fluence of the pitch rate on the drag force is ignored. The side force and roll moment coefficients versus the roll rate are plotted in F 8. The side force and roll moment coefficients decrease linearly with the increase roll rate. In this study, the roll rate has little effect on the yaw moment coefficient, s ignored here. The side force, roll moment and yaw moment coefficients versus yaw rate ar sented in Figure 9. With the increase in yaw rate, the side force coefficient and yaw ment coefficient decrease linearly, while the roll moment coefficient increases linea The side force, roll moment and yaw moment coefficients versus yaw rate are presented in Figure 9. With the increase in yaw rate, the side force coefficient and yaw moment coefficient decrease linearly, while the roll moment coefficient increases linearly.
All the values of the aerodynamic derivatives calculated from Fluent software and the results with 95% confidence bounds fitted by MATLAB are summarized in Table 2. Angles and control surface deflections are given in radians, and angular rates are expressed in terms of radians per second when calculating the aerodynamic derivatives. Among these coefficients, C mα < 0, C nβ >0 and C lβ < 0 indicates the UAV has a longitudinal, yaw and lateral static stability. Machines 2021, 9, x FOR PEER REVIEW 11 of 28 All the values of the aerodynamic derivatives calculated from Fluent software and the results with 95% confidence bounds fitted by MATLAB are summarized in Table 2. Angles and control surface deflections are given in radians, and angular rates are expressed in terms of radians per second when calculating the aerodynamic derivatives. Among these coefficients, Cm α < 0, Cn β >0 and Cl β < 0 indicates the UAV has a longitudinal, yaw and lateral static stability ..

Propulsion System Measurements
The propulsion system components consist of a propeller, brushless DC (Direct Current) motor, Li-Po (Lithium Polymer) battery and ESC (Electronic Speed Control) unit. In addition to the propulsion system components, the actuator of the elevon is also tested in this section, as shown in Figure 10. The propulsion system components were selected according to the structure, size and weight of the UAV investigated in this study. A SunnySky 2216 KV2400 electric brushless motor with a 6040 propeller and a Grepow 2200 mA·h 4S (16.8V) Li-Po battery were adopted and tested in this research.

Propeller Model
Generally, the thrust generated by the propeller during the flight of the UAV depends not only on the rotating speed of the propeller but also on the airspeed [35]. Since the airspeed of the UAV used in this study is relatively small, only the influence of the rotating speed on the thrust is considered to simplify the model. The thrust and the counter torque exerted by the propeller [36] can be obtained by Equation (13).

Propulsion System Measurements
The propulsion system components consist of a propeller, brushless DC (Direct Current) motor, Li-Po (Lithium Polymer) battery and ESC (Electronic Speed Control) unit. In addition to the propulsion system components, the actuator of the elevon is also tested in this section, as shown in Figure 10. The propulsion system components were selected according to the structure, size and weight of the UAV investigated in this study. A SunnySky 2216 KV2400 electric brushless motor with a 6040 propeller and a Grepow 2200 mA·h 4S (16.8 V) Li-Po battery were adopted and tested in this research. All the values of the aerodynamic derivatives calculated from Fluent software and the results with 95% confidence bounds fitted by MATLAB are summarized in Table 2. Angles and control surface deflections are given in radians, and angular rates are expressed in terms of radians per second when calculating the aerodynamic derivatives. Among these coefficients, Cm α < 0, Cn β >0 and Cl β < 0 indicates the UAV has a longitudinal, yaw and lateral static stability ..

Propulsion System Measurements
The propulsion system components consist of a propeller, brushless DC (Direct Current) motor, Li-Po (Lithium Polymer) battery and ESC (Electronic Speed Control) unit. In addition to the propulsion system components, the actuator of the elevon is also tested in this section, as shown in Figure 10. The propulsion system components were selected according to the structure, size and weight of the UAV investigated in this study. A SunnySky 2216 KV2400 electric brushless motor with a 6040 propeller and a Grepow 2200 mA·h 4S (16.8V) Li-Po battery were adopted and tested in this research.

Propeller Model
Generally, the thrust generated by the propeller during the flight of the UAV depends not only on the rotating speed of the propeller but also on the airspeed [35]. Since the airspeed of the UAV used in this study is relatively small, only the influence of the rotating speed on the thrust is considered to simplify the model. The thrust and the counter torque exerted by the propeller [36] can be obtained by Equation (13).

Propeller Model
Generally, the thrust generated by the propeller during the flight of the UAV depends not only on the rotating speed of the propeller but also on the airspeed [35]. Since the airspeed of the UAV used in this study is relatively small, only the influence of the rotating speed on the thrust is considered to simplify the model. The thrust and the counter torque exerted by the propeller [36] can be obtained by Equation (13).
where D p denotes the diameter of the propeller (m), C T and C M are dimensionless coefficients representing thrust and torque, N denotes the rotating speed (RPM). Since D p and ρ are constants, T and M can be expressed as quadratic functions of N with coefficient K T and K M , respectively, as shown in Equation (14).
In order to acquire the aerodynamics of the propeller, an electric propulsion test bench made by the Wing-Flying company is investigated and utilized, as shown in Figure 11. Table 3 gives its nominal forces/moments and accuracy, along with the basic parameters of the test bench.
where Dp denotes the diameter of the propeller (m), CT and CM are dimensionless coefficients representing thrust and torque, N denotes the rotating speed (RPM). Since Dp and ρ are constants, T and M can be expressed as quadratic functions of N with coefficient KT and KM, respectively, as shown in Equation (14).
In order to acquire the aerodynamics of the propeller, an electric propulsion test bench made by the Wing-Flying company is investigated and utilized, as shown in Figure  11. Table 3 gives its nominal forces/moments and accuracy, along with the basic parameters of the test bench. Figure 11. Electric propulsion measurement system.  Figure 12 shows the test results of the propeller, where Figure 12a is the thrust curve and Figure 12b is the torque curve. The thrust force and torque increase with the increase of the rotating speed and the fitting results are shown in Table 4 with 95% confidence bounds of KT and KM.   Figure 12 shows the test results of the propeller, where Figure 12a is the thrust curve and Figure 12b is the torque curve. The thrust force and torque increase with the increase of the rotating speed and the fitting results are shown in Table 4 with 95% confidence bounds of K T and K M .

Motor Model
The parameters that need to be measured for the motor include the relationship tween throttle and rotating speed, the response curve of the motor and the relations between the rotating speed, torque, current and voltage. The propeller is installed on motor in all tests. The equipment used for the motor measurement is the same as the t bench utilized for the propeller.

Motor Model
The parameters that need to be measured for the motor include the relationship between throttle and rotating speed, the response curve of the motor and the relationship between the rotating speed, torque, current and voltage. The propeller is installed on the motor in all tests. The equipment used for the motor measurement is the same as the test bench utilized for the propeller. Figure 13 shows the relationship between the rotating speed of the motor and the throttle setting. Since the UAV usually works in the cruise status for a long time during the flight missions and the high throttle is rarely used, we limit the maximum throttle setting in the test to 0.8. The motor starts to work when the throttle increases to 0.09, and the throttle below 0.09 is the throttle dead zone. In order to simplify the modeling of the motor, the throttle-speed curve is linearly fitted.

Motor Model
The parameters that need to be measured for the motor include the relationship between throttle and rotating speed, the response curve of the motor and the relationship between the rotating speed, torque, current and voltage. The propeller is installed on the motor in all tests. The equipment used for the motor measurement is the same as the test bench utilized for the propeller. Figure 13 shows the relationship between the rotating speed of the motor and the throttle setting. Since the UAV usually works in the cruise status for a long time during the flight missions and the high throttle is rarely used, we limit the maximum throttle setting in the test to 0.8. The motor starts to work when the throttle increases to 0.09, and the throttle below 0.09 is the throttle dead zone. In order to simplify the modeling of the motor, the throttle-speed curve is linearly fitted. The motor is approximately a first-order system in this study and the step response curve of the motor is shown in Figure 14. It can be obtained from the response curve that the time constant Tm of the motor system is equal to 0.19 with 95% confidence bounds between 0.129 and 0.243. The motor is approximately a first-order system in this study and the step response curve of the motor is shown in Figure 14. It can be obtained from the response curve that the time constant T m of the motor system is equal to 0.19 with 95% confidence bounds between 0.129 and 0.243. Step response curve of the brushless motor.
In order to calculate the energy consumption of the propulsion system, it is vital to obtain the voltage and current of the motor. In this research, the motor and ESC are simplified into one model, and the equivalent voltage and equivalent current can be calculated by Equation (15) when the thrust force and torque are given [37]. Step response curve of the brushless motor.
In order to calculate the energy consumption of the propulsion system, it is vital to obtain the voltage and current of the motor. In this research, the motor and ESC are simplified into one model, and the equivalent voltage and equivalent current can be calculated by Equation (15) when the thrust force and torque are given [37].
where U m0 is the nominal no-load voltage (V), I m0 is the nominal no-load current (A), R m is the internal resistance (Ω), K V0 is the KV value of the motor (RPM/V), M m denotes motor load torque (N·m), N m denotes the rotating speed of the motor, which is equal to the rotating speed of the propeller (RPM). U m0 is 10 V, I m0 is 1.7A, R m is 72.5 mΩ and K V0 is 2400 RPM/V. The model of the motor equipped in this research can be obtained as shown in Equation (16).

Servo Motor Model
Two DC servo motors are equipped in the UAV to drive aileron deflection and elevator deflection. The response characteristics of the servo motors will directly affect the attitude control of the UAV, so the modeling of the servo motors is very important. The modeling of the servo motor is conducted by a second-order system approach [38] and the step response is used to identify the model.
A Pixhawk flight controller equipped with a transmitter and receiver is utilized to provide PWM signal to servo motors, and an MPU6050 sensor is adopted to obtain the rotation angle, as shown in Figure 15. It should be noted that the connecting cables for data acquisition might affect the test results due to the limitation of the hardware. The magnitude of the influence caused by the cables is quite small compared to the torque of the servo motor, and we ignore it here.    Y(s) X(s) = 9.774 2 s 2 +2 × 9.774 × 0.801s+9.774 2 (17) where the natural frequency and damping of the servo motor are 9.77 and 0.801, respectively. The servo motor is an under-damped plant, and it can track the command with small oscillation.  Figure 16 represents the transient state response of the rotation angle of the servo motor with a step input. The transfer function of the servo motor adopted in terms of the natural frequency and damping is shown in Equation (17). Y s X s s s (17) where the natural frequency and damping of the servo motor are 9.77 and 0.801, respectively. The servo motor is an under-damped plant, and it can track the command with small oscillation.

Battery Model
The battery is the core and foundation of the propulsion system. Acquiring the discharge dynamics of the battery is crucial for the entire flight process. Moreover, we can also predict the endurance of the UAV for reference. A Li-Po battery was equipped and investigated in this work due to its high energy density. The main equation [39] of the battery is shown in Equation (18) and assumes the following: the temperature and the internal resistance of the battery are constant; neither the aging effect nor the Peukert effect [40] were considered. Step response curve of the servo motor.

Battery Model
The battery is the core and foundation of the propulsion system. Acquiring the discharge dynamics of the battery is crucial for the entire flight process. Moreover, we can also predict the endurance of the UAV for reference. A Li-Po battery was equipped and investigated in this work due to its high energy density. The main equation [39] of the battery is shown in Equation (18) and assumes the following: the temperature and the internal resistance of the battery are constant; neither the aging effect nor the Peukert effect [40] were considered.
where V batt denotes the battery's actual voltage; E 0 is the battery constant voltage; K is the polarization constant; C denotes the maximum battery capacity; i denotes actual constant discharge current and it is the discharged capacity; A is the exponential voltage and B represents the exponential capacity; R is the internal resistance. In order to obtain the battery discharge characteristics, it is necessary to record the voltage of the battery during discharging. An Ultra Power UP400AC charger was adapted to discharge the battery at a constant current. Its maximum discharge power is 40 W, and the maximum discharge current can reach 8 A. The propulsion test bench was also utilized to record the real-time voltage of the battery. Figure 17 shows the layout of the test apparatus. In this test, the battery is discharged with a constant current of 2 A until the single cell of the battery drops to 3.5 V, which is approximately completely discharged. (18) where Vbatt denotes the battery's actual voltage; E0 is the battery constant voltage; K is the polarization constant; C denotes the maximum battery capacity; i denotes actual constant discharge current and it is the discharged capacity; A is the exponential voltage and B represents the exponential capacity; R is the internal resistance.
In order to obtain the battery discharge characteristics, it is necessary to record the voltage of the battery during discharging. An Ultra Power UP400AC charger was adapted to discharge the battery at a constant current. Its maximum discharge power is 40 W, and the maximum discharge current can reach 8 A. The propulsion test bench was also utilized to record the real-time voltage of the battery. Figure 17 shows the layout of the test apparatus. In this test, the battery is discharged with a constant current of 2 A until the single cell of the battery drops to 3.5 V, which is approximately completely discharged.  Figure 18a presents the experimental data and the fitted parametric curve. The parameters can be obtained by fitting the test data (time input and voltage output). Table 5 shows the optimal values of the discharge curve. The error between the predicted battery capacity and the nominal capacity is 0.41%, and the error between the predicted discharge  Figure 18a presents the experimental data and the fitted parametric curve. The parameters can be obtained by fitting the test data (time input and voltage output). Table 5 shows the optimal values of the discharge curve. The error between the predicted battery capacity and the nominal capacity is 0.41%, and the error between the predicted discharge current and the actual current is 3.95%. The obtained battery model has high accuracy. The discharge curves at different discharge rates can be simulated through the obtained battery model, as shown in Figure 18b.  Figure 18a presents the experimental data and the fitted parametric curve. Th rameters can be obtained by fitting the test data (time input and voltage output). Ta shows the optimal values of the discharge curve. The error between the predicted ba capacity and the nominal capacity is 0.41%, and the error between the predicted disch current and the actual current is 3.95%. The obtained battery model has high accu The discharge curves at different discharge rates can be simulated through the obta battery model, as shown in Figure 18b.

Mass Moment of Inertia Measurements
Solving the moment equation requires obtaining the UAV's mass moment of inertia. A simple approach [41] was adopted to measure the mass moment of inertia of the UAV in this work. The prototype was suspended from two points equidistant from its center of gravity. Then, the prototype was rotated by a small angle around its center of gravity and released. The moment of inertia can be calculated from the recorded oscillation period. We set a vertical rotation axis x v perpendicular to y body axis to make the angle with x body axis is α v and the angle with z body axis is 90 • −α v , and the value of I xz can be gained by Equation (19). Since the UAV is symmetrical about the x-z plane in the body frame, I xy and I yz can be neglected. Table 6 shows the measurement results of the UAV's moment of inertia.
where I xx , I yy and I x v are the moment of inertia of the x body axis, y body axis and x v axis, respectively.

Simulink Model
All the system parameters obtained from the simulations and experiments are imported into Simulink. The nonlinear model consists of six main blocks: the controller block, which contains stability and upper controllers; the propulsion system block, which contains the models of the propulsion system components; the environment block, which contains wind disturbances; the aerodynamic forces and moments calculation block, which calculates the aerodynamic forces and moments generated during flight; the 6-DOF equations of motion block, which contains the previously derived equations; the visualization block, which applies the co-simulation of FlightGear and MATLAB to exhibit the flight attitude of the UAV. Figure 19 shows the overall block diagram of the UAV in Simulink.

Control Methods of the Fixed Wing UAV
Based on the established simulation model, we can efficiently design and verify the control algorithm in Simulink software. In this section, three controllers with different functions are designed to control the movement of the UAV intelligently. The simulation step is set to 1 ms and the specific simulation results are presented in the following subsections.

Control Methods of the Fixed Wing UAV
Based on the established simulation model, we can efficiently design and verify the control algorithm in Simulink software. In this section, three controllers with different functions are designed to control the movement of the UAV intelligently. The simulation step is set to 1 ms and the specific simulation results are presented in the following subsections.

Stability Controller Design
The stability controller is the most fundamental and important controller of UAV. Its performance directly determines the safety and reliability of UAVs during flight. The stability controller designed in this research can be divided into three parts: longitudinal controller, lateral controller and speed controller. The longitudinal controller includes pitch rate control, pitch angle control and height control and the lateral controller includes roll rate control, roll angle control and course angle control, while the speed controller only controls the airspeed of the aircraft. PID control method is applied to the stability controller because it is efficient and easy to implement. Figure 20 shows the specific structure of the stability controller.
The Least Square Method based on MATLAB was employed to gain PID parameters (K p , K i , K d ). On the basis of the acquired PID parameters, the parameters are adjusted according to the expected performance by the Response Optimization toolbox in Simulink. The three parts of the stability controller are adjusted separately. In the longitudinal and lateral controllers, the PID parameters of the angular rate loop, angular loop and position loop are tuned successively. Figure 21 displays the step response of the pitch angle based on the Least Square Method and the toolbox. After tuning by Least Square Method, the response curve has a relatively small rise time, but a large overshoot and oscillation. The system's overshoot and oscillation can be reduced when optimized through the toolbox. Since the PID tuning of other controllers is similar to the tuning of the pitch angle, the results will not be given here and all the parameters of the stability controller are given in Table A1 of the Appendix B.

Stability Controller Design
The stability controller is the most fundamental and important controller of UAV. Its performance directly determines the safety and reliability of UAVs during flight. The stability controller designed in this research can be divided into three parts: longitudinal controller, lateral controller and speed controller. The longitudinal controller includes pitch rate control, pitch angle control and height control and the lateral controller includes roll rate control, roll angle control and course angle control, while the speed controller only controls the airspeed of the aircraft. PID control method is applied to the stability controller because it is efficient and easy to implement. Figure 20 shows the specific structure of the stability controller. The Least Square Method based on MATLAB was employed to gain PID parameters (Kp, Ki, Kd). On the basis of the acquired PID parameters, the parameters are adjusted according to the expected performance by the Response Optimization toolbox in Simulink. The three parts of the stability controller are adjusted separately. In the longitudinal and lateral controllers, the PID parameters of the angular rate loop, angular loop and position loop are tuned successively. Figure 21 displays the step response of the pitch angle based on the Least Square Method and the toolbox. After tuning by Least Square Method, the response curve has a relatively small rise time, but a large overshoot and oscillation. The system's overshoot and oscillation can be reduced when optimized through the toolbox. Since the PID tuning of other controllers is similar to the tuning of the pitch angle, the results will not be given here and all the parameters of the stability controller are given in Table A1 of the Appendix B. Step responses of pitch angle under different methods. Figure 22 shows responses of the airspeed, course angle and altitude of the UAV under respective desired commands. As can be seen in Figure 22a, the UAV can quickly follow the desired command for an increase in airspeed by increasing the motor speed.  The Least Square Method based on MATLAB was employed to gain PID parameters (Kp, Ki, Kd). On the basis of the acquired PID parameters, the parameters are adjusted according to the expected performance by the Response Optimization toolbox in Simulink. The three parts of the stability controller are adjusted separately. In the longitudinal and lateral controllers, the PID parameters of the angular rate loop, angular loop and position loop are tuned successively. Figure 21 displays the step response of the pitch angle based on the Least Square Method and the toolbox. After tuning by Least Square Method, the response curve has a relatively small rise time, but a large overshoot and oscillation. The system's overshoot and oscillation can be reduced when optimized through the toolbox. Since the PID tuning of other controllers is similar to the tuning of the pitch angle, the results will not be given here and all the parameters of the stability controller are given in Table A1 of the Appendix B. Step responses of pitch angle under different methods. Figure 22 shows responses of the airspeed, course angle and altitude of the UAV under respective desired commands. As can be seen in Figure 22a, the UAV can quickly follow the desired command for an increase in airspeed by increasing the motor speed. Step responses of pitch angle under different methods. Figure 22 shows responses of the airspeed, course angle and altitude of the UAV under respective desired commands. As can be seen in Figure 22a, the UAV can quickly follow the desired command for an increase in airspeed by increasing the motor speed. However, there will be a big delay when the expected airspeed drops; this is because the airspeed of the UAV can only be reduced by wind resistance. Since the UAV investigated in this study achieves a turn by controlling the roll of the body, the tracking speed of the course angle is relatively slow, as shown in Figure 22b. Detailed step response values such as rising time, steady-state error and peak overshoot are given in Table 7 to evaluate the performance of the stability controller. However, there will be a big delay when the expected airspeed drops; this is because the airspeed of the UAV can only be reduced by wind resistance. Since the UAV investigated in this study achieves a turn by controlling the roll of the body, the tracking speed of the course angle is relatively slow, as shown in Figure 22b. Detailed step response values such as rising time, steady-state error and peak overshoot are given in Table 7 to evaluate the performance of the stability controller.    The biggest interference encountered by the UAV during flight missions is atmospheric disturbance. Turbulence is one of the atmospheric disturbances and its effects on all parts of the fixed-wing UAV are continuous and uneven. This paper mainly studies the stability and robustness of the stability control system when the UAV encounters turbulence. For solving engineering issues, Dryden model and Von Karman model are usually implemented to simulate turbulence [42]. Dryden model is applied in this study to test and verify the controller because it is easy to implement in simulation. Dryden continuous filter can be expressed as following transfer function: where σ u , σ v and σ w denote turbulence intensities; L u , L v and L w denote scale length. The turbulence intensities are related to flying height and typical values for wind speed at a height of 20 feet; specific parameter settings can be found in [43]. Figure 23 displays the response of the UAV's airspeed, course angle and altitude under the effect. The airspeed is most affected by the effects of turbulence. The response of the airspeed is higher than the desired command with a maximum error of 3.5 m/s, which is because the turbulence can directly affect the airspeed and it is slow to track the descending airspeed command. The course angle and altitude will oscillate around the desired command under the influence of turbulence, and the course angle has a maximum error of 0.02 rad. Since the UAV is not in a steady state at the beginning of the simulation, it will continuously adjust to stabilize the attitude and enter the steady state at 5 s. The height error during this period is not only due to the influence of turbulence. The maximum error of height is 0.11 m at 31.2 s. The UAV system can run stably with the controller under turbulent conditions. where σu, σv and σw denote turbulence intensities; Lu, Lv and Lw denote scale length. The turbulence intensities are related to flying height and typical values for wind speed at a height of 20 feet; specific parameter settings can be found in [43]. Figure 23 displays the response of the UAV's airspeed, course angle and altitude under the effect. The airspeed is most affected by the effects of turbulence. The response of the airspeed is higher than the desired command with a maximum error of 3.5 m/s, which is because the turbulence can directly affect the airspeed and it is slow to track the descending airspeed command. The course angle and altitude will oscillate around the desired command under the influence of turbulence, and the course angle has a maximum error of 0.02 rad. Since the UAV is not in a steady state at the beginning of the simulation, it will continuously adjust to stabilize the attitude and enter the steady state at 5 s. The height error during this period is not only due to the influence of turbulence. The maximum error of height is 0.11 m at 31.2 s. The UAV system can run stably with the controller under turbulent conditions.

Path following Controller Design
Realizing the position control of aircraft is the ultimate objective of autonomous flight control. In this part, a path following controller was designed and verified for the fixedwing UAV based on the stability controller. The path following controller adds autonomy and intelligence to the UAV, which can make the UAV fly along a desired straight or circular path. For the UAV investigated in this study, the course angle is the basis and core of path following. The lateral distance between the UAV and the reference path can be reduced by controlling the course angle. The vector field guidance law proposed in reference [44] was investigated and implemented in this work.

Path following Controller Design
Realizing the position control of aircraft is the ultimate objective of autonomous flight control. In this part, a path following controller was designed and verified for the fixedwing UAV based on the stability controller. The path following controller adds autonomy and intelligence to the UAV, which can make the UAV fly along a desired straight or circular path. For the UAV investigated in this study, the course angle is the basis and core of path following. The lateral distance between the UAV and the reference path can be reduced by controlling the course angle. The vector field guidance law proposed in reference [44] was investigated and implemented in this work. Figure 24 shows the vector field constructed by the guidance strategy when following a straight line. The UAV needs to keep tracking the direction of the vectors to make the flight direction consistent with the vector direction, and the UAV can finally fly along the desired path. The expected course angle χ d of the aircraft can be expressed by Equation (21) in the vector field.
where χ ∞ denotes the reference course angle far from the waypoint path; k path denotes the turning coefficient of following a straight-line path; e y denotes the lateral distance error; χ q denotes the course angle of the straight path. The values of χ ∞ and k path can affect the performance of following a straight path. Relatively appropriate values are selected for the two parameters according to the specific situation of the aircraft and the influence of different values on the flight mission is not given here.   Figure 25a shows the desired path and tracking trajectory and Figure  25b shows the distance error. It can be seen that the distance error increases in the initial stage of the simulation which is because the initial speed of the aircraft will keep the aircraft away from the path. The response curve of the distance error is not only related to the parameters of the vector field, but also closely related to the flight speed of the aircraft. In the simulation, the distance error converges towards zero with 0.1 m steady-state error when the airspeed is 15 m/s. The vector field guidance law and its parameter settings can satisfy the flight mission of following a straight-line path in this research. Following an orbit path is more difficult than a straight-line path. The vector field of the circular path is similar to that of the straight-line path and is not shown here. The expected course angle can be expressed by Equation (22) when following an orbit path.
where φorbit is the angle of the vector from the aircraft to the center of the circle relative to the north direction; λ is the turn direction, and λ = 1 signifies clockwise direction and λ =   Figure 25a shows the desired path and tracking trajectory and Figure 25b shows the distance error. It can be seen that the distance error increases in the initial stage of the simulation which is because the initial speed of the aircraft will keep the aircraft away from the path. The response curve of the distance error is not only related to the parameters of the vector field, but also closely related to the flight speed of the aircraft. In the simulation, the distance error converges towards zero with 0.1 m steady-state error when the airspeed is 15 m/s. The vector field guidance law and its parameter settings can satisfy the flight mission of following a straight-line path in this research.   Figure 25a shows the desired path and tracking trajectory and Figure  25b shows the distance error. It can be seen that the distance error increases in the initial stage of the simulation which is because the initial speed of the aircraft will keep the aircraft away from the path. The response curve of the distance error is not only related to the parameters of the vector field, but also closely related to the flight speed of the aircraft. In the simulation, the distance error converges towards zero with 0.1 m steady-state error when the airspeed is 15 m/s. The vector field guidance law and its parameter settings can satisfy the flight mission of following a straight-line path in this research. Following an orbit path is more difficult than a straight-line path. The vector field of the circular path is similar to that of the straight-line path and is not shown here. The expected course angle can be expressed by Equation (22) when following an orbit path.
where φorbit is the angle of the vector from the aircraft to the center of the circle relative to Following an orbit path is more difficult than a straight-line path. The vector field of the circular path is similar to that of the straight-line path and is not shown here. The expected course angle can be expressed by Equation (22) when following an orbit path.
where ϕ orbit is the angle of the vector from the aircraft to the center of the circle relative to the north direction; λ is the turn direction, and λ = 1 signifies clockwise direction and λ = −1 signifies counter-clockwise direction; k orbit is turning coefficient of following an orbit path; d is the distance from the aircraft to the center of the orbit path; r is the radius, r is selected as 100 m according to the flight task requirements in this paper. However, the application of the vector field guidance law to follow the circular trajectory has a large steady-state error in the simulation. This is because the desired course angle will change with the position of the aircraft when the aircraft follows an orbit path, and the response of the UAV to the course angle command has a certain delay, resulting in the UAV not being able to track the desired course angle in time. L1 guidance law is implemented in this work to solve the problem, and specific details of L1 guidance law can be found in reference [45]. It should be noted that relatively suitable parameters have been selected through simulation tests for the two guidance strategies, and the influence of different parameters is not given here. Figure 26 shows the simulation results of implementing the two guidance laws to track a circular trajectory with a radius of 100 m, where Figure 26a represents the flight trajectory and Figure 26b represents the distance error. The UAV can follow the orbit path with both guidance laws. After 140 s of simulation, the distance error of the vector field guidance law is 16.83 m, while that of L1 guidance law is only 0.56 m.
Machines 2021, 9, x FOR PEER REVIEW 23 can be found in reference [45]. It should be noted that relatively suitable parameters h been selected through simulation tests for the two guidance strategies, and the influe of different parameters is not given here. Figure 26 shows the simulation results of implementing the two guidance law track a circular trajectory with a radius of 100 m, where Figure 26a represents the fl trajectory and Figure 26b represents the distance error. The UAV can follow the orbit p with both guidance laws. After 140 s of simulation, the distance error of the vector f guidance law is 16.83 m, while that of L1 guidance law is only 0.56 m. Although the application of the L1 guidance law has a smaller steady-state error, only suitable for the case that the aircraft is not far from the path due to some assumpt made in its derivation process. In this study, the two guidance strategies are combi The vector field guidance law is used when the distance error is greater than 40 m, the L1 guidance law is utilized in other cases. Figure 27 shows the flight result of the craft with a start point away from the orbit path, where Figure 27a shows the flight tra tory and Figure 27b shows the distance error. The steady-state error is 0.68 m during st flight, which can meet the requirements of the flight mission. Although the application of the L1 guidance law has a smaller steady-state error, it is only suitable for the case that the aircraft is not far from the path due to some assumptions made in its derivation process. In this study, the two guidance strategies are combined. The vector field guidance law is used when the distance error is greater than 40 m, and the L1 guidance law is utilized in other cases. Figure 27 shows the flight result of the aircraft with a start point away from the orbit path, where Figure 27a shows the flight trajectory and Figure 27b shows the distance error. The steady-state error is 0.68 m during stable flight, which can meet the requirements of the flight mission.
only suitable for the case that the aircraft is not far from the path due to some assumpt made in its derivation process. In this study, the two guidance strategies are combi The vector field guidance law is used when the distance error is greater than 40 m, the L1 guidance law is utilized in other cases. Figure 27 shows the flight result of the craft with a start point away from the orbit path, where Figure 27a shows the flight tr tory and Figure 27b shows the distance error. The steady-state error is 0.68 m during st flight, which can meet the requirements of the flight mission.

Path Management Controller Design
On the basis of the path following controller, a path management controller was signed in this part. The path management controller allows the UAV to fly between ferent waypoints, which can help the UAV complete some designated cruise tasks. F path management controller, the most important thing is to determine whether the U has reached the previous waypoint and whether it needs to switch to the next wayp

Path Management Controller Design
On the basis of the path following controller, a path management controller was designed in this part. The path management controller allows the UAV to fly between different waypoints, which can help the UAV complete some designated cruise tasks. For a path management controller, the most important thing is to determine whether the UAV has reached the previous waypoint and whether it needs to switch to the next waypoint. A simple algorithm was implemented for path management [46] in this work to accomplish the flying mission. Figure 28 shows the situation of UAV switching from one waypoint to the next waypoint.
Machines 2021, 9, x FOR PEER REVIEW 24 of 28 A simple algorithm was implemented for path management [46] in this work to accomplish the flying mission. Figure 28 shows the situation of UAV switching from one waypoint to the next waypoint.
where wi−1, wi, wi+1 represent the coordinates of the waypoints. H(wi, ni) is the plane perpendicular to ni through point wi. When the UAV flies over the plane H(wi, ni), which indicates that the UAV has reached waypoint wi at this time and it will fly towards the next waypoint wi+1. Figure 29 shows the trajectory of the aircraft and performance metrics when completing a flight mission. The aircraft starts from the starting point, climbs to the desired altitude, flies along the desired path, and finally returns to the initial waypoint. The switching between waypoints can be realized by the algorithm; however, the aircraft will fly a short distance to adjust the position when it reaches the waypoints, which is the shortcoming of the path management controller. The maximum distance error is 47.26 m when the UAV switches from one waypoint to another. Adjusted distance refers to the distance the At this time, the aircraft flies from point w i−1 to point w i , and vector n i can be represented by Equation (23).
where w i−1 , w i , w i+1 represent the coordinates of the waypoints. H(w i , n i ) is the plane perpendicular to n i through point w i . When the UAV flies over the plane H(w i , n i ), which indicates that the UAV has reached waypoint w i at this time and it will fly towards the next waypoint w i+1 . Figure 29 shows the trajectory of the aircraft and performance metrics when completing a flight mission. The aircraft starts from the starting point, climbs to the desired altitude, flies along the desired path, and finally returns to the initial waypoint. The switching between waypoints can be realized by the algorithm; however, the aircraft will fly a short distance to adjust the position when it reaches the waypoints, which is the shortcoming of the path management controller. The maximum distance error is 47.26 m when the UAV switches from one waypoint to another. Adjusted distance refers to the distance the UAV needs to fly to reduce the distance error to within 3 m.
indicates that the UAV has reached waypoint wi at this time and it will fly towards the next waypoint wi+1. Figure 29 shows the trajectory of the aircraft and performance metrics when completing a flight mission. The aircraft starts from the starting point, climbs to the desired altitude, flies along the desired path, and finally returns to the initial waypoint. The switching between waypoints can be realized by the algorithm; however, the aircraft will fly a short distance to adjust the position when it reaches the waypoints, which is the shortcoming of the path management controller. The maximum distance error is 47.26 m when the UAV switches from one waypoint to another. Adjusted distance refers to the distance the UAV needs to fly to reduce the distance error to within 3 m.

Discussion
Through simulation and test results, there are some issues that need to be discussed and improved.

1.
Some of the aerodynamic coefficient curves obtained by the CFD method in Section 3 were simplified in order to facilitate modeling, which would bring some errors. It is an effective method to verify the CFD simulation results through additional wind tunnel tests to increase the accuracy of the aerodynamic model.

2.
When obtaining the parameters of the propulsion system and the servo motor through corresponding tests, the layout and performance of the test equipment might affect the test results. For example, the dynamic response of the servo motor would be affected by connecting cables used for data acquisition.

3.
A simple logic was used in the path management controller to determine whether the drone has reached the waypoint so that the UAV can switch to the next waypoint; however, the aircraft would fly a short distance to adjust the position when it reached the waypoint, which was the shortcoming of the path management controller.
For future work, a functional prototype will be developed and assembled, and flight tests will be conducted to evaluate the control algorithms and measure the difference between the simulation model and the real aircraft.

Conclusions
In this study, a complete design process was proposed and implemented for a small electric fixed-wing UAV. Firstly, 6-DOF equations of motion and nonlinear aerodynamic model of the fixed-wing UAV investigated in this work were derived, which were the basis for aircraft analysis and design.
The CFD software Fluent was used to quickly analyze the aerodynamic characteristics of the UAV and obtain the aerodynamic coefficients based on the CFD method. Among them, the rotating reference frame was applied to gain the damping derivatives by simulation of the steady movement. The parameters in the simulation were set according to the actual flight conditions, such as the angle of attack and angles of control deflections. The aerodynamic model of the UAV was well established by obtaining aerodynamic coefficients through the CFD method. The propulsion system of the aircraft was tested through different experiments, and the mathematical model of each component was established by using the test data. Some simplifications were carried out for the convenience of modeling when obtaining the motor and servo motor models. Based on Simulink software, a complete simulation model for assisting design and verifying the algorithms was quickly established by importing all the test data.
The stability controller could be well applied to the nonlinear model of the unconventional UAV with only two elevons and ensure the UAV to fly stably under the turbulent condition. Vector field guidance law and L1 guidance law were combined to follow a circular trajectory and the distance error was 0.68 m during stable flight. It could solve the problems that vector field guidance law had a large steady-state error and overcome the limitations of the assumption of the L1 guidance law. The simulation results indicate that all the control algorithms are easy to implement with excellent performance, which can meet the requirements of basic flight missions. A complete nonlinear model of the small battery-powered flying wing UAV was well established in this paper to solve the fatal problem of the small application range of the linearized model, and stability controller, path following controller and path management controller were designed and optimized. These modeling and control methods of the small electric fixed-wing UAV can provide some reference for relevant designers. In addition, all the UAV parameters measured in our work are presented in this paper, which can be used to assist in the design of the controller and verify the effectiveness of the control algorithm.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Detailed Transformation Matrices of the UAV's Mathematical Model
The detailed expression of the transformation matrix from the body axis to Earth axis is shown in (A1). The detailed expression of transformation matrix from wind axis to body axis is shown in (A3).