C Gellida-Coutiño

V Dominguez-De La Cruz

A Sanchez-Orta

O Garcia-Salazar

P Castillo
email: pedro.castillo@hds.utc.fr

The tailsitter autogiro UAV: modeling, design, and CFD simulation

This paper focuses on the development of a novel UAV configuration, the convertible tailsitter autogiro, based on simulation and modeling of the most significant aerodynamic and dynamic effects that can affect the autogiro flight during the expected flight conditions. The aerodynamic forces acting on the rotor are simulated using a novel algorithm based on the blade element momentum (BEM) theory (to simulate wind turbines), and on semi-empirical models used for helicopter rotors in autorotation. The tail design is derived from a model of the velocity vector field in the slipstream of the propeller. Finally, aerodynamic control surfaces in the tail are evaluated using a transient k-epsilon Computational Fluid Dynamics (CFD) simulation, which validates the design computation.

I. INTRODUCTION

The autogiro is a rotary wing aircraft equipped with an engine-driven propeller, which provides forward thrust, and a passive rotor (unpowered) that rotates as a result of the inflow wind due to the aircraft motion, including external wind currents. It is not capable of performing hover flight, however, it is safer than both helicopters and fixed-wing air vehicles due to the absence of the stall velocity, which arises in airplanes at low speeds, and the dead man's curve, which occurs in helicopters. Additionally, the autogiro consumes less energy at moderate speeds [START_REF] Gordon | Development of the autogiro: A technical perspective[END_REF]. Nevertheless, autogiros have some disadvantages, including the inability to achieve the same speeds as fixed-wing aircraft (as it is the case with helicopters), moreover, autogiros need runways during takeoff and landing, which is shorter than those used by airplanes.

A few air vehicles with passive rotors were proposed to overcome the disadvantages of autogiros; the most welldocumented vehicles are the McDonnell XV-1 Convertiplane [START_REF] Franklin | An overview of autogyros and the mcdonnell xv-1 convertiplane[END_REF] and the Fairey Rotodyne [START_REF] Gibbings | The fairey rotodyne-technology before its time?[END_REF], both incorporate a tip jet engine, which ejects a jet of hot gases and thus enables vertical take off and landing. Despite the fact that both vehicles were capable of flying, they faced problems, especially noise and vibration.

Several researches have been conducted research on unmanned autogiros vehicles, which have a better performance than multirotors performing foward flight [START_REF] Petritoli | Unmanned autogyro for advanced sar tasks: a preliminary assessment[END_REF], [START_REF] Cai | Control and dynamics analysis for miniature autogyro and compound autogyro[END_REF], [START_REF] Petritoli | Unmanned autogyro for mars exploration: A preliminary study[END_REF]. Additionally, companies such as Carter Aviation Technologies (Carter A.T.) and PAL-V are actively developing V-TOL human scale vehicles based on the autogiro concept. For instance, in Carter A.T., a hybrid airplane-autogiro was offered [START_REF] Jr | Electric motor powered rotor drive for slowed rotor winged aircraft[END_REF], whilst PAL-V proposed an autogiro equipped with wheels, transmission, and collapsible propeller and rotor that could be used as a car [START_REF] Ju | Flying car related technology trends[END_REF].

To address the inability of fixed-wing aircraft to achieve the vertical takeoff and landing, convertible vehicles were proposed. These vehicles can be classified as tilt-rotors, tiltwing, and tailsitters [START_REF] Richard | Design and flight testing of a convertible quadcopter for maximum flight speed[END_REF], in which the first two vehicles have a high mechanical complexity and add the weight of additional motors. Convertible tailsitters are mechanically identical to fixed-wing aircraft, but they require a vertical-to-horizontal transition to switch from hover to horizontal flight and vice versa. The main problem with tailsitters is that their wings have a stall velocity, which is the lowest velocity at which they can generate considerable lift force, which means that the vehicle must fly fast enough in order to achieve a smooth transition.

To model the rotor aerodynamics of both helicopters and wind turbines, the blade element momentum (BEM) theory is a feasible option. Nevertheless, this theory fails to describe the autorotation. To overcome this limitation, a semi-empirical engineering correction was proposed in [START_REF] Gordon J Leishman | Principles of helicopter aerodynamics with CD extra[END_REF], however, the proposed strategy is solved through linearization and simplification of the transcendental equations in BEM theory and using only one control volume. In contrast, in [START_REF] Hansen | Aerodynamics of wind turbines[END_REF], [START_REF] Nørkaer | The general momentum theory[END_REF] a more sophisticated computational algorithm to solve the transcendental equations was presented. This algorithm allows the control volume to be divided into small control volumes and the blades into sections, however, it cannot directly be applied to describe the autorotation. The formulation of an algorithm to solve the BEM equations with small control volumes and the blade analyzed in sections, as well as the implementation of engineering corrections for autorotation is proposed as part of our contributions.

In this paper, our tailsitter autogiro (TA) is a novel concept that performs the transition from horizontal to vertical flight in a soft manner since the passive rotor can provide an important amount of thrust at low flight speeds. Moreover, the angular velocity of the rotor can be directly measured, allowing the estimation of the rotor thrust, which is useful information for the control algorithm. The main objective of this work is to propose the development of a small TA UAV by using its aerodynamic-dynamic model and simulations of the flight conditions. To the best of our knowledge, this proposed vehicle has not been published; thus, the proposed approach has not been explored before.

The paper is divided into the following sections: A description of the design algorithm employed to develop the TA is described in section II. The proposed algorithm to estimate the aerodynamic properties of the rotor is presented in section III. The dynamic model using the Newton-Euler formulation is derived in section IV. A novel simple model to describe the flapping motion of blades is presented in section V. The selection of blades and their information generated by applying the proposed novel algorithm to solve the BEM equations are described in section VI. The mass positioning problem is solved and exposed in section VII, while the dimensions and location of the wings in the tail are described in section VIII, which also discusses the fuselage design. A CFD simulation is performed to determine the effect of moment generation on tail and fuselage drag; the results are presented in section IX. Finally, a general conclusion is presented at the end of the paper.

II. DESIGN ALGORITHM

The design objectives include the weight restriction of the vehicle and its endurance, range and autonomy. The problem of design is how to select the design variables such that the design objectives are satisfied while the restrictions imposed by the closed-loops of relations are maintained. A design algorithm based on the spiral principles of design is proposed in [START_REF] Howe | Aircraft conceptual design synthesis[END_REF]. This algorithm yields a solution to the design problem.

The first step of the design algorithm is to select the blades based on the target weight (for instance the blades of a commercial RC autogiro). The second step is to select the motor-propeller so that it provides twice the target weight as thrust force. The third step is to select the battery based on the target time-flight for hover and the propeller-motor power consumption. The fourth step is to define the weight of the fuselage-tail subsystem as the subtraction of the weight of the selected components from the target weight. The fifth step is to propose a tail and fuselage geometry based on aerodynamic models. The sixth step is to propose materials and the structure for the fuselage-tail subsystem to obtain the strength and stiffness that allow flight. This step must be repeated in order to maintain the weight inside a target range. If the full weight exceeds the target weight plus an acceptable gap when the sixth step iterations reach the selected limit, then the rotor must be replaced and the design algorithm begins again, otherwise, a solution of the design is found.

III. ALGORITHM FOR SIMULATION OF ROTOR

AERODYNAMICS

An important problematic to analyze the TA rotor is that the momentum theory fails exactly in the region when autorotation appears; thus, the classical algorithms used to study helicopter rotors or wind turbines cannot be directly applied. In [START_REF] Gordon J Leishman | Principles of helicopter aerodynamics with CD extra[END_REF] a model was presented to correct the momentum theory in order to describe the autorotation; however, to the best of our knowledge, this model has not been combined with an algorithm to solve the rotor aerodynamics by segmenting the blades and the rotor disk. Therefore, the correction for autorotation combined with an algorithm is presented in this section.

The first step of this algorithm is to modify the hover velocity definition proposed in [START_REF] Gordon J Leishman | Principles of helicopter aerodynamics with CD extra[END_REF] in order to describe the flux on partial annular volume control

v h = ∆T 2ρ ∆A (1)
where ∆A = πr∆r is the area of an annular section of the rotor, ρ is the air density and ∆T the axial force generated by such annular section. The climbing ratio in annular control volumes can be then defined as follows

l = v c /v h (2)
where v c is the local climbing velocity for a given blade section, also known as the free wind velocity for an observer in the section blade frame. The induced velocity at the rotor, in agreement with [START_REF] Gordon J Leishman | Principles of helicopter aerodynamics with CD extra[END_REF], can be estimated as

v i v h =                          - l 2 + l 2 2 + 1 if l > 0 κ + κ 1 l + κ 2 l 2 + κ 3 l 3 + κ 4 l 4 if -2 ≤ l ≤ 0 - l 2 - l 2 2 + 1 if l < -2
(3) where the case -2 ≤ l ≤ 0 corresponds to the transition region between the wind mill and helicopter regime, κ = 1, κ 1 = -1.125, κ 2 = -1.372, κ 3 = -1.718, κ 4 = -0.655. From [START_REF] Hansen | Aerodynamics of wind turbines[END_REF], the tangential induced velocity, v t , can be computed as follows:

v t = ∆M 2ρ ∆A(v c -v i ) ( 4
)
where ∆M is the aerodynamic moment of the force acting on an annular section of the rotor. Therefore, both the axial and tangential induced velocities at a given annular section of the rotor are estimated. These quantities can be stacked in the induced velocity vector w = (v i , v t ) T , where v i and v t are the magnitude of the axial and tangential velocity, respectively. Thus, the following equations are used [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]:

∆T = B(L cos(ϕ) + D sin(ϕ)) 2πrF ∆r (5)
∆M = B(L sin(ϕ) -D cos(ϕ)) 2πrF ∆r ( 6
)
where F is the Prandtl tip lost factor, B is the number of blades, ∆r is the width of annulus (which is the same as the length of blade sections), r is the mean radio of annulus, ϕ is the flux angle given by

ϕ = arctan v c -v i ωr + v t ( 7
)
where ω is the angular velocity of the rotor, while L and D are the lift and drag forces produced by each blade section, respectively, and they can be estimated as follows:

L, D = 1 2 ρ (v c -v i ) 2 + (ωr + v t ) 2 v 2 rel c C i ∆r ( 8
)
where c is the length of chord, C i , for i ∈ {L, D}, stands for the lift and drag coefficients, and v rel is the magnitude of the relative wind at each blade section. Finally, the proposed algorithm is given by the pseudocode III.1:

Algorithm III.1: BEM SOLUTION(T , M = f (v c )) w ← (0, 0) T T ← 0 M ← 0 for r ← r min to r max do                                                            while e < T OL do                                                ϕ with (7) L and D with (8) ∆M with (6) ∆T with (5) v h with (1) l with (2) v i with (3) v t with (4) e ← (v i , v t ) T -w w ← w + k p e e ← ||e|| 2 i ← i + 1 if i > i max then break T ← T + ∆T M ← M + ∆M r max = R
and ideally r min should be 0, but the relative wing speed is almost axial in the vicinity of the rotor center, posing a convergence problem for the algorithm, thus, r min = 0.35R is proposed. i max is the maximum number of iterations (i max = 100 is suggested to achieve convergence). e is an error that feeds a pseudo proportional control, which leads to a soft approximation of the solution. The proportional constant, k p , must be tuned in such a way that convergence problems caused by the lack of continuity in (3) are overcome.

IV. DYNAMIC MODEL OF THE TAILSITTER AUTOGIRO

The TA is considered as a pair of rigid bodies; the main body and the rotor. A reference frame is assigned to each body, denoted by Σ a and Σ r , respectively. Σ a is placed exactly at the center of mass of the TA, and Σ r is situated at the middle of the rotor head. Using the Newton-Euler formulation, the equation of motion to describe the attitude dynamics is given by

m i:a = I a ωa + ω a × I a ω a + P R (9)
where Σm i:a is the summation of all the moments of force acting on the TA, I a is the inertia tensor of the aerial vehicle excluding the rotor. ω a is the angular velocity vector expressed in a reference frame whose origin coincides with the center of mass, and P R is a disturbance due to rotor dynamics, which is computed as follows

P R = R rot|a (I rot ωrot + ω rot × I rot ω rot ) (10)
where R rot|a is the rotation matrix that maps the vector coordinates from Σ r to Σ a , I rot is the inertia tensor of the rotor, and ω rot is the angular velocity vector of the rotor expressed as

ω rot = ω r ẑ + R rot|a ω a (11)
and

ωrot = R rot|a ( ωa -(ẑ • ωa )ẑ) (12)
where ω r is the magnitude of the rotor angular velocity in the frame of the rotor, ẑ is a unit vector in the direction of the z axis of Σ a (which is the same for Σ r ), and R rot|a is a rotation matrix that maps coordinate vectors from Σ r to Σ a . The angular acceleration of the rotor is neglected in [START_REF] Nørkaer | The general momentum theory[END_REF] since no torque is transmitted from the rotor to the main body. Moreover, the forces acting on the rotor are considered to be only parallel to ẑ, forces acting in other directions are absorbed by the blade flapping motion (which is described in section V). Additionally, the rotor dynamics is considered to be strongly determined by the aerodynamic forces and depends on the free wind. The moments acting on the autogiro are expressed as

m i:a =   d w d w d w -d t cos(30 o ) d t cos(30 o ) 0 d t sin(30 o ) d t sin(30 o ) d t     L 1 L 2 L 3   (13)
where d w is the distance between the vehicle center of mass and the point of intersection of the aerodynamic axis of the wings (PIAA). d t is the distance between the PIAA and the mean aerodynamic radius of the wings. L 1 , L 2 , and L 3 are the lift forces of wings. Finally the position dynamics of TA is described as follows where m a is the mass of the TA (including the rotor) and p a is the position of the center of mass with respect to the origin of the inertial frame, which is expressed in the coordinates of Σ a and derived in the inertial frame. f i is the summation of the forces acting on the autogiro that is given by

f i = m a pa (14)
f i = k p ω 2 p x + k R ω 2 r ẑ - ṗa | ṗa | D a ( ṗa • ṗa ) + m a g ( 15
)
where k p is the thrust constant of the propeller, ω p is the magnitude of the propeller angular velocity, x is a unit vector in the direction of the x axis of Σ a , D a is a drag constant of the fuselage and g = R 0|a (0 0 -9.81)

T is the gravitational acceleration vector. R 0|a is a rotation matrix that maps vectors from the inertial reference frame to Σ 0 .

V. BLADES FLAPPING MODEL

Flapping is the sinusoidal motion of the blades in the same direction as the rotor angular velocity vector, which is caused by an asymmetric moment that arises due to the different relative air velocity acting on each blade. This motion allows autogiros and helicopters to counteract the asymmetric moment by producing an aerodynamic moment with the same frequency and magnitude but in the opposite direction. Despite this benefit, it can cause a collision between the rotor and the autogiro tail, thus, the design process must take into account for the maximum amplitude of oscillations. To accomplish this, in this section, a model, which can be used in autogiros, helicopters, and other rotorcraft of any size, is proposed.

The main objective of this method is to evaluate the case where the blades are flapping but the autogiro is static. Under this condition, a static drag force can be estimated, which generates a static drag moment (M Drag ). The flapping motion amplitude can be determined from the static drag moment since it has similar magnitudes to the asymmetric relative wind speed but acts in opposite directions. To derive the proposed model, it is considered that the sum of two wind vector fields has the same aerodynamic forces as those of each vector field, which results from a quasi-linear dependence between lift and drag functions with the angle of attack for small angles, [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. In addition, the algorithm of the BEM equations presented in section VI is solved assuming small angles in order to obtain the flapping angle that cancels both the asymmetric moment and the static drag moment.

To apply the proposed model, the first step is to define the flapping angle, η, as the angle formed by the tip of the blade relative to itself when the blade is fully tense. Then, the blade drag constant, K D , is defined, which is used to estimate the moment about the "flap axis" caused by the drag force presented in (17) (7.7e -3kg • m 2 for the current blade) Fig. 2: The schematic diagram of the flapping motion.

K D = Bρ R 0 r c dr ( 16
)
which is a general expression, also valid for non rectangular blades. The flapping static moment due to drag force M Drag as a function of flapping angular velocities can be computed as

M Drag = -K D η2 sign( η) (17)
Function M Drag ( η) in ( 17) is evaluated within the range η ∈ [-0.3ω max , 0.3ω max ], where ω max is the rotor maximum angular velocity, defined as the value of the angular velocity at which the thrust equals twice the weight of the TA. The obtained curve is linearized using a curve fitting algorithm that corresponds to a first-degree polynomial. The slope of the obtained function is denoted by K bl . Finally, the maximum amplitude of the flapping angle is given by

η max = - V 0,max T max ω 2 max K bl (18)
where V 0,max is the maximum free wind velocity of the rotor and T max is the maximum thrust of the rotor.

VI. ESTIMATION OF ROTOR PROPERTIES VIA

SIMULATION

The selected blades have a constant chord of 4 cm and no torsion, which are classified as rectangular blades. To characterize the airfoil, a photograph of the cross section of a blade is taken and used to generate a point cloud, which is then compared with the airfoils in an airfoil tools database. The airfoil GOE593 is the closest to the point cloud.

The algorithm presented in section III is evaluated for different values of v c and ω using a function that considers a non-dimensional quantity called tip speed ratio (TSR) as an independent variable,

T SR = ωR v c (19)
and the thrust and moment coefficients, c T and c M , respectively, as dependent variables

c T = T 0.5ρv 2 c A (20)
c M = M R 0.5ρv 3 c A (21)
The algorithm is solved for different number of blades B = [2, 3], and using T OL = 10 -4, k p = 0.2, i max = 100 and T SR = [START_REF] Cai | Control and dynamics analysis for miniature autogyro and compound autogyro[END_REF][START_REF] Ss Houston | Identification of autogyro longitudinal stability and control characteristics[END_REF]. Fig. 4 and Fig. 3 show the results of the proposed algorithm, where the obtained points are interpolated using a spline interpolation. The algorithm is implemented for two rotors, one with two and the other with three blades. However, for the the three-bladed rotor case, no solution to the design problem has been obtained. Analyzing Fig. 4 for the three-bladed rotor case, a function that is similar to a straight line with a negative slope can be observed; however, this function must decrease in magnitude until it approaches the point (0, 0), otherwise the law of conservation of energy is violated; this is more evident for the two-bladed rotor case. It can be observed from Fig. 4 that a unique root is obtained at T SR = 10.5, for the case of the two-bladed rotor. Therefore, this value represents a stable equilibrium point of the rotor dynamics denoted as T SR eq .

The obtained point by evaluating T SR eq at the thrust coefficient function in Fig. 3 is denoted as the equilibrium thrust coefficient c T,eq . The equilibrium thrust coefficients of both rotors are nearly equal; in this case c T = 1.14.

In order to estimate how the rotor behaves when the flow is not axial, a new semi-empirical model is proposed. This model approximates the observed trend in the experimental data reported in [START_REF] Gordon | Development of the autogiro: A technical perspective[END_REF]:

c T (α) = c T,eq sin π 2 tanh α r π 2 (22)
Evaluating ( 22) over some range and using [START_REF] Ss Houston | Identification of autogyro longitudinal stability and control characteristics[END_REF], an estimation of how the rotor thrust evolves at different values of v c and α r is obtained and presented in Fig. 5. Using the curves in Fig. 5, it is possible to deduce that a valid range of weights can be lifted by the blades is from 450gf up to 2000gf; however, in order to avoid exceeding the mechanical resistance of the blade, the weight of the vehicle is considered as 0.55 Kg with a gap of ±0.1 Kg. Under this consideration, the TA is able to produce sufficient lift with the main rotor at 4 m/s. Using c T,eq and the data in Fig. 3, the angular velocity of the rotor is plotted against the rotor thrust, resulting the data presented in Fig. 6, where a quadratic relationship is observed. This demonstrates that the climb velocity has no significant effect on the quadratic tendency. A curve fitting allows determining a quadratic coefficient denoted as b rotor , which is b rotor = 8.6e -5 Ns/rad, for the case of the threebladed rotor, and b rotor = 5.8e -5 Ns/rad for the case of the two-bladed rotor. Knowing the gyroscopic moment generated by the main rotor is essential for designing the aerodynamic control surfaces; this variable depends on the angular velocity vector of both ω rot and ω a and can be expressed as

M giro = I r ω r ω a × ẑ ( 23
)
where I r is the moment of inertia about the axis of rotation of the main rotor, which coincides with the TA z axis.

There is no physical reason for limiting rotor gyroscopic moments; however, the control algorithm reduces the inflow angle when the rotor thrust increase in an undesirable manner, maintaining, indirectly, the gyroscopic moment stable. Using (23) and the quadratic curve in Fig. 6, the magnitude of gyroscopic moment is evaluated as a function of the rotor thrust and the vehicle angular velocity, see Fig. 7. From Fig. 7 and with the estimated maximum moment of the tail derived in section IX, a maximum vehicle angular velocity of 0.8 rad/s and a maximum rotor thrust of 700 grams are proposed. According to the design strategy, the motor-propeller, the battery, and other important components are selected; see table I.

The TA can be equipped with either a pusher or a tractor propeller. However, in agreement with recent evidence on small propellers, [START_REF] Khan | Development and validation of a propeller slipstream model for unmanned aerial vehicles[END_REF], the slipstream has a greater transverse area as it moves away from the propeller plane. The greater the area of aerodynamic control surface (ACS) submerged in the slipstream, the greater the magnitude of the moment of force provided by the tail. Thus, in order to generate the maximum possible moment, the ACS must be located sufficiently away from the propeller. As a result, the tractor propeller is more convenient for tailsitter applications, thus, this propeller is selected for the proposed UAV. Additionally, designing a tail with the propeller on the front side is easier.

To avoid the moments of force caused by unbalanced payloads, the center of mass of the TA, CM V , is established at the intersection of the propeller thrust line (PTL) and the rotor thrust line (RTL), which are generally orthogonal in every autogiro when the tilt angle is zero; otherwise, an undesired moment of force is produced, see Fig. 8. Using a Cartesian coordinate system with PTL and RTL as axes and the origin at their intersection, an equation for estimating the center of mass over RTL is deduced

CM v = L r m rotor -L b m batt -L e m emb + i∈E L i m emb (24)
where L b , L r , and L e denote the distance between PTL and the battery, rotor, and embedded system centers of mass, respectively. The mass of the battery, rotor, and embedded system are denoted by m batt , m rotor and m emb , respectively. E is a set containing the remainder elements, which are assumed to be distributed symmetrically around RTL, thus, the summation result is approximately zero. Even when it is desirable to locate the battery and embedded components as close as possible to the RTL, this is difficult if the condition CM v = (0, 0) is imposed, as there are some factors that prevent L r from being small enough. To propose a minimum L r (L r,min ), the following assumptions must be established: collisions between the blades and tail are avoided, and the rotor is kept out of the slipstream of the propeller.

The blade coning angle, θ, is the angle formed by the blade longitudinal axis and the plane perpendicular to the RTL. It accomplishes two objectives: to reduce the mechanical loads on the blades and to prevent collisions with the tail. The cone angle is generally equal to the angle formed by the thrust force generated by a blade and the centrifugal force acting on it. Because both the centrifugal and thrust forces depend quadratically on the rotor angular velocity, a constant expression for such an angle is found

θ ideal = arctan 2 b rotor m pala R B ( 25
)
where b rotor is the quadratic coefficient defined in section VI. The ideal cone angle should be 5 deg for three and twobladed rotors. Applying the proposed blade flapping model of section V, the maximum amplitude of the flapping motion is 0.04 rad or 2.29 degrees. Therefore, a cone angle of 5 deg guarantees a balance of the forces over the blades and overcomes the flapping angle.

VIII. FUSELAGE AND CONTROL SURFACES MODEL

DESIGN

The slipstream model proposed in [START_REF] Khan | Development and validation of a propeller slipstream model for unmanned aerial vehicles[END_REF] is codified. The results of evaluating this model are shown in Fig. 9 and Fig. 10, where it can be seen that the vector field approaches zero at radial distances greater than 20 cm, indicating that this is the minimum distance to place the rotor-head. The control surfaces in the tail will be referred as wings for practical purposes. The selected airfoil for these wings is the NACA0012. Since the vector field described in Fig. 9 and Fig. 10 is almost zero at a radius greater than 0.2 m, this value is chosen as the wing span. The movable part of wings begins from 75% of the chord to the trailing edge. Three wings are chosen to give stability when the TA is on the ground, each separated by 120 deg. The lift and drag coefficients of this airfoil, with the flap turning at 10 deg and 30 deg and at various angles of attack, are obtained through simulations that are conducted using XFOIL, which is a professional software to solve panel method equations to approximate the aerodynamic properties of airfoils [START_REF] Drela | Xfoil: An analysis and design system for low reynolds number airfoils[END_REF]. Using the model ( 13) and the blade element theory, the following equations to estimate the moment provided by the tail wings are obtained

M roll = 1 N w L 0 0.5 W act1 ρ c v 2 slip C L y dy (26) M pitch , M yaw = 1 N w L 0 0.5 W act2 ρ c v 2 slip C L x dy (27)
where W acti , for i ∈ [START_REF] Gordon | Development of the autogiro: A technical perspective[END_REF][START_REF] Franklin | An overview of autogyros and the mcdonnell xv-1 convertiplane[END_REF], denotes the number of wings capable of fully providing force to generate moment in that direction, which is three for i = 1 and two for i = 2, as a result of the selected wing angles. N w is the number of wings. Evaluating equations ( 26) and ( 27) for different distances between the propeller and tail, and chord length, Fig. 11 and Fig. 12 are obtained The fuselage is shaped based on the NACA0020 airfoil rotating around its chord. The length of the airfoil of the fuselage is 0.6m, which provides sufficient space for the wings to be attached at the calculated distance. The fuselage is constructed using planar pieces cut from balsa wood sheets. To facilitate simple replacement of the battery or internal components, the fuselage is divided into two equal halves, as it is illustrated in Fig. 13. Balsa wood is chosen as the primary material for the fuselage, which is covered with Monokote paper. Extruded polystyrene is used for the aerodynamic control surfaces.

IX. FUSELAGE AND CONTROL SURFACES CFD

SIMULATION

Transient CFD simulation is performed in order to characterize the forces and moments generated by the control surfaces in relation to the TA center of mass and to validate the proposed model for the tail. This characterization allows the TA hover flight behavior to be predicted. Numerical CFD simulations are carried out using ANSYS Fluent [START_REF]ANSYS Fluent Theory Guide[END_REF].

To conserve computational resources during the analysis, a simplified model is considered, see Fig. 15. The rotor is omitted from this analysis. Three distinct control volumes are used in the simulation: one fixed for the fuselage and wings, and two movable for each coaxial propeller. The mesh is generated using the Fluent meshing default setting tool, with body sizes of 1e-4 m for the propellers and 1e-3 m for the remaining bodies.

To solve second-order equations and incompressible flow, a simulation is developed using a pressure-based algorithm and the k -epsilon turbulence model. The system is initialized using hybrid initialization. Up to 20 iterations are performed, and 500 steps are defined with a time value of 0.001 s.

The propellers are turning at different angular speeds in order to avoid the reaction torque observed in previous simulations. It is observed that a relationship of 1.13 exists between the angular velocities, ensuring the least amount of reaction torque.

CFD Results

To visualize the flux around each wing, three planes were chosen, as illustrated in Fig. 16. Fig. 16: Independent planes generated over every control surface.

Pressure contours are generated for each generated plane, as illustrated in Fig. 16. Additionally, a velocity contour is generated to quantify the velocity homogeneity around three ailerons in Fig. 18. Finally, velocity contours around control surfaces are generated at a distance of 10 cm from the tailsitter center of mass. This enables visualization of the air stream impacting the control surfaces and generating the quantified moments. It is possible to observe that a higher velocity is allocated at the rear of the deflected aileron, implying that this zone receives less pressure, as can be seen in Fig. 1920. The combined thrust generated by the two propellers is 628 grams-force, simulating hover flight. The moment generated by the three wings is 0.096 Nm per wing on roll moment, which matches acceptably to the 0.089 Nm predicted in section VIII, indicating that the applied technique is valid.

X. CONCLUSIONS

The proposed model that describes the aerodynamic forces in hover flight was compared with CFD results, and an acceptable match was found. The proposed algorithm for rotor aerodynamics was implemented to analyze the selected blades, inferring important performance curves of the main rotor to guide the design process. Even if the models are valid for an operation region near to hover, the obtained results confirm that the main rotor can provide sufficient lift at 4 m/s, where the models on the tail are valid. Nevertheless, experimental studies must be performed in order to validate the proposed design methodology. The main goal of this research was the design of a novel tailsitter autogiro configuration based on models and simulation results, in which the tail overcomes the dynamic and aerodynamic disturbance of the passive rotor.

Fig. 1 :

 1 Fig. 1: The schematic diagram of the TA.

Fig. 3 :

 3 Fig. 3: Thrust coefficient vs tip speed ratio obtained by solving the proposed algorithm.

Fig. 4 :

 4 Fig. 4: Moment coefficient vs tip speed ratio obtained by solving the proposed algorithm.

Fig. 5 :

 5 Fig. 5: Estimated thrust force of the rotor at different angles of attack and flight velocity, for both 2 and 3 blades.

Fig. 6 :

 6 Fig.6: Thrust vs rotor angular velocity using the data obtained from the solution of the proposed algorithm

Fig. 7 :

 7 Fig. 7: Magnitude of the gyroscopic moment due to the main rotor thrust and the vehicle angular velocity.

Fig. 8 :

 8 Fig. 8: Main concepts to define the conceptual design.

Fig. 9 :

 9 Fig. 9: Magnitude of the velocity vector field of slipstream estimated at 4000 RPM.

Fig. 10 :

 10 Fig. 10: Magnitude of the velocity vector field of slipstream estimated at 6100 RPM.

Fig. 11 :Fig. 12 :Fig. 13 :Fig. 14 :

 11121314 Fig. 11: Estimation of the maximum roll moment generated by the aerodynamic control surfaces on the slipstream.

Fig. 15 :

 15 Fig. 15: Tailsitter autogiro simplified model for CFD analysis.

Fig. 17 :

 17 Fig. 17: Pressure contour for 240 deg aileron.

Fig. 18 :

 18 Fig. 18: Velocity contour for 240 deg aileron.

Fig. 19 :

 19 Fig. 19: Velocity contour for 120 deg aileron at 10 cm away from center.

Fig. 20 :

 20 Fig. 20: Velocity contour for 240 deg aileron at 10 cm away from center.

TABLE I :

 I Selected actuators and avionic components

	Type	Component	weight (g)
	Brushless motor	CRM2413A	170 g
	ESC	-	60 g
	Battery	LiPo 11.1V (3S)	180g
	Embedded	NodeMCU32-Gyro-RF	50g
	Servomotor	SG90	9g
	VII. CONCEPTUAL DESIGN AND CENTER OF MASS
		PROBLEMATIC


E U Rojo-Rodriguez

E G Rojo-Rodriguez

E J Ollervides-Vazquez

P Castillo
email: castillo@hds.utc.fr

O Garcia-Salazar

E U Rojo-Rodriguez

E G Rojo-Rodriguez

Obstacle-aware UAV navigation based on a leader-follower consensus protocol

This work aims to present a navigation scheme whose topmost layer further considers the presence of obstacles and steers the aircraft towards a safe path by altering its navigation references. The complete framework is divided into two control layers carrying out the orientation and navigation, while a third obstacle-aware layer, which is based on a Virtual Leader-Follower coordination problem, generates the positioning references for the navigation controller. This outline allows designing the inner control layers independently, without obstacle awareness, giving flexibility to different combinations of UAVs and controller laws. To validate the proposed structure, numerical simulations are performed where moving obstacles are considered. In addition, experimental testing is carried out by implementing the three-layer navigation in a quadrotor aircraft, performing an independent trajectory while in the presence of a static obstacle. Both tests show acceptable results, and the obstacle collision avoidance is reached.

I. INTRODUCTION

There are numerous applications and fields where an Unmanned Aerial Vehicle (UAV) could be used. In the past decade, many researchers have focused their goals toward the solution or study of different application problems with specific criteria, such as navigation, obstacle avoidance, artificial vision, load carrying, among similar, [START_REF] Zhilenkov | Intelligent autonomous navigation system for uav in randomly changing environmental conditions[END_REF], [START_REF] Radmanesh | Overview of path-planning and obstacle avoidance algorithms for uavs: a comparative study[END_REF], [START_REF] Guerrero-Sánchez | Nonlinear control strategies for a uav carrying a load with swing attenuation[END_REF].

Controlling the dynamics of the aircraft has proven to be the core problem to solve. Nonetheless, very successful developments have already proven the feasibility of their implementation and usage; different approaches ranging from traditional Proportional Integral Derivative (PID) controllers to more complex non-linear schemes have been tested and shown different advantages and benefits [START_REF] Babu | Designing of self tuning pid controller for ar drone quadrotor[END_REF], [START_REF] Gonzalez-Hernandez | Real-time altitude control for a quadrotor helicopter using a super-twisting controller based on high-order sliding mode observer[END_REF]. These advances in the field allow the development of more sophisticated algorithms and applications that were not fully covered. Even though it may seem simple enough to implement an already existing scheme, different scenarios require a specific combination of UAV and corresponding control strategies, preventing form standardization of higher-level algorithms.

This drawback is considerable enough when security concerns arise as a consequence of taking the implementation outside the controlled testing environment. In an actual UAV application, the complete mission must be performed, accounting for external disturbances, unexpected faults, and all possible unwanted interactions with the working environment, such as obstacles. The latter gains special attention as it represents a high percentage of UAV mission failures [START_REF] Hu | Risk assessment model for uav cost-effective path planning in urban environments[END_REF], so avoiding obstacle collisions is a navigation problem that has been heavily discussed in recent years. Having a different combination of aircraft types with even more diverse approaches of controller laws leads to designing obstacle avoidance strategies specifically for the implementation; as a consequence of this, it is challenging to come up with a general usage algorithm.

There are many schemes developed to get a suitable response against an obstacle. For instance, the most common approach is the anticipated path planning, where the trajectory followed by the UAV is designed with some degree of obstacle awareness. One example is the work presented in [START_REF] Krishnan | Path planning of an autonomous quadcopter based delivery system[END_REF], where an offline path planning is carried out by means of a combination of a rapidly-exploring random tree (RRT) and potential fields to generate the most efficient path, this for a delivery quadrotor; although the authors reached their proposed goal, the main disadvantage become evident as all obstacles locations must be known before the algorithm is calculated. A similar development was presented in [START_REF] Bhandari | Path-planning around obstacles for a quadrotor uav using the rrt algorithm for indoor environments[END_REF], in which an RRT algorithm is also used, in which the detection of static obstacles is also considered non-previously known obstacles, with the limitation of being static due to the nature of the path generation algorithm. This type of approach could also be implemented using different kinds of schemes, from evolutionary algorithms such as the work of [START_REF] Gutierrez-Martinez | Collision-free path planning based on a genetic algorithm for quadrotor uavs[END_REF], and [START_REF] Chen | Obstacle avoidance strategy for quadrotor uav based on improved particle swarm optimization algorithm[END_REF], or even direct interaction algorithms, more focused on reaction time, and avoiding the collision with simple commands, [START_REF] Falanga | Dynamic obstacle avoidance for quadrotors with event cameras[END_REF].

While all these developments present a suitable way of solving the avoidance problem, they must be specifically designed or at least adapted to the specific mission and aircraft used. This has been addressed by defining working constraints; for instance, the dynamic obstacle avoidance solves the problem of moving obstacles, and in many cases, the algorithms are designed as velocity references most aircraft can follow; always considering its type constraints, such as turning radius for fixed-wing UAVs. Applications such as the presented in [START_REF] Mac | Improved potential field method for unknown obstacle avoidance using uav in indoor environment[END_REF] or in [START_REF] Abdellatif | Artificial potential field for dynamic obstacle avoidance with mpc-based trajectory tracking for multiple quadrotors[END_REF], validate the usage of virtual fields to steer the UAV far from the obstacle, in which, if the navigation is independent, the field only needs to consider velocity constraints. On the other hand, the most recent development is focused on mechanisms of self-adaptation and learning as a possible solution against obstacles, which can be considered uncertainties. As an example, authors in [START_REF] Dai | Automatic obstacle avoidance of quadrotor uav via cnn-based learning[END_REF] developed a framework of avoiding dynamic obstacles by sensing and identifying the obstacle using a convolutional neural network, and with the relative positioning data, determine the necessary turning angle and relative velocity. More recent developments deal with the problem by applying the Reinforced Learning approach, as a way of altering the ideal reference the UAV follows, [START_REF] Rubí | Quadrotor path following and reactive obstacle avoidance with deep reinforcement learning[END_REF], [START_REF] Zheng | Policy-based monocular vision autonomous quadrotor obstacle avoidance method[END_REF].

The aim of this work is directed towards the proposal of a UAV navigation scheme that can be considered independent of the inner control layers of the aircraft while at the same time assuring the success of the navigation in an environment where static and dynamic obstacles could be present. A 3-layer control scheme is proposed, where the topmost layer corresponds to a reference generation, which in this sense can be separated from the remaining layers, so different configurations can be implemented. The algorithm is based on a multi-agent system (MAS) approach, where virtual coordination is defined between the ideal trajectory reference, considered as a virtual leader, and the actual navigation states of the UAV, while considering a virtual velocity vector in each obstacle, as added velocities in the coordination protocol.

The remaining document is divided as follows: in Section II necessary mathematical tools are presented and defined to further develop the reference generator. Section III deals with the navigation carried out by the reference generation, which is based on a coordinated MAS, while validations through numerical simulations and a real-time experiment are described in Section IV. Conclusions are finally presented in Section V.

II. MATHEMATICAL TOOLS

A. Multi-rotor UAV dynamics

As an starting procedure, it is necessary to define the dynamic equation governing the aircraft, which for this case a quadrotor configuration is used. Different reference frames are considered while defining the equations of motion; for instance, a ground-fixed inertial framed defined as I i ={x Ii , y Ii , z Ii } in which translational effects are modeled, a body frame attached to the center of gravity of the UAV as B i ={x Bi , y Bi , z Bi }, and an aerodynamic frame, W i ={x Wi , y Wi , z Wi }, that considers external forces such as wind gusts, changes in angle of attack and sideslip angle, and other aerodynamic effects, [START_REF] Stengel | Flight dynamics[END_REF].

Newton-Euler formulation is used to describe the equations of motion of the UAV as follows

ξi = V i (1) m i Vi = (-T Ti )R i e 3 + m i ge 3 + D ξi (2) Ṙi = R i Ωi (3) J i Ωi = -Ω i × J i Ω i + τ ai + D ηi (4)
with i = 1, . . . , N corresponding to an specific UAV if a formation of multiple aircrafts is defined, ξ i = (x i , y i , z i ) ∈ R 3 are the position coordinates relative to the inertial fixed frame and η i = (φ i , θ i , ψ i ) ∈ R 3 describes the rotation states of the i-th UAV, which are described by an orthogonal rotation matrix R i ∈ SO(3) : B i → I i parametrized by the Euler angles φ i roll, θ i pitch and ψ i yaw. Angular velocities in the body frame B i are defined as Ω i = (p i , q i , r i ) ∈ R 3 , while linear velocities in the inertial fixed frame I i as V i = ( ẋi , ẏi , żi ) ∈ R 3 . T Ti ∈ R >0 is the total thrust of all the UAV actuator, and τ ai ∈ R 3 are the moments relative due the body frame, and as a result of actuators combinations.

Vectors of canonical basis of R 3 are also considered, and are represented by e 1 ,e 2 , and e 3 . The term m i ∈ R defines the mass of the i -th UAV, while J i ∈ R 3×3 describes the moments of inertia. Ωi is the skew-symmetric matrix associated with the cross product (i.e. [START_REF] Lozano | Unmanned aerial vehicles: Embedded control[END_REF], [START_REF] Ramirez-Rodriguez | Robust backstepping control based on integral sliding modes for tracking of quadrotors[END_REF], [START_REF] Sanchez-Orta | Position-yaw tracking of quadrotors[END_REF].

âb = a × b ∀a, b ∈ R 3 ). D ξi = (d ξi 1 , d ξi 2 , d ξi 3 ) ∈ R 3 and D ηi = (d ηi 1 , d ηi 2 , d ηi 3 ) ∈ R 3 are disturbances vectors,

B. The General Consensus Protocol

Multi-agent systems are widely used in the coordinated flight of UAVs. Many approaches have been proposed; for instance, the Consensus Protocol scheme has been proven to be a reliable and robust approach to achieve coordinated motion.

On its most basic structure, a Consensus protocol gets the following form

u ζ d i = K Pi ( N j=1 a ij (ζ j -ζ i )) + K Di ( N j=1 a ij ( ζj -ζi )) (5)
where K Pi is a positive diagonal matrix whose elements guarantee the consensus on the position states, K Di represents a positive diagonal matrix containing terms for reaching consensus of velocity states. In addition, ζ i represents the states of the agent, while ζ j the state of its neighbor agent, with i = 1, 2, 3, . . . , N , j = 1, 2, 3, . . . , N , being N the total number of agents in the formation. It is worth to note that a ij sets the communication strength from agent j to i.

C. Virtual-leader formation flight

One of the most common approaches to achieving a coordinated flight by using a multi-agent system is the so-called Leader-Follower Consensus (LFC), whose main objective is to drive the i -th follower to the states of a time-varying reference, that can be a simulated agent or an ideal numerical reference. The advantages of this model rely on the simplicity of the protocol design due to the independence of the leader, which in this sense, partially solves some phenomena such as noise amplification through the multi-agent system communication topology.

From the point o view of the Consensus protocol, no changes are introduced in the basic structure, despite the fact that the neighbor agent is the virtual leader and thus is not subjected to the consensus. Nevertheless, it is important to guarantee that the actual agent is reaching a consensus with the virtual leader, meaning the following definition can be addressed.

Definition 1: An agent is said to reach consensus with the virtual leader if there exist a Consensus protocol [START_REF] Gonzalez-Hernandez | Real-time altitude control for a quadrotor helicopter using a super-twisting controller based on high-order sliding mode observer[END_REF], such that the closed-loop system satisfies

lim t→∞ ||(ζ j (t) -ζ i (t)) -ξ i (0)|| = 0 (6)
for all possible geometric formation constrains, considered as the initial condition vector ξ i (0).

For the case of this implementation, no formation distances are considered, so ξ i (0) = 0.

III. NAVIGATION BASED ON A CONSENSUS PROTOCOL

As an alternative, to drive the states of the UAV to a secure position, considering its reference vector and the position of obstacles, a 3-layer control scheme is proposed. Rewriting equations ( 1)-( 4), these are described as follows

ξi = V i (7) Vi = u pi + d ξi (8) Ṙi = R i Ωi (9) Ωi = u ai + d Ri ( 10
)
with u pi = ge 3 -

T T i mi (R i e 3 ), τ ai = Ju ai , d ξi = D ξ i mi and d Ri = J -1 i [-Ω i × J i Ω i + D ηi ]
. This leads to the definition of virtual control inputs u pi ∈ R 3 and u ai ∈ R 3 for the translational dynamics and rotational dynamics, respectively. As the most inner and fast layer u ai , a local closed-loop controller is designed, whose reference vector is the output of the following controller loop u pi , which deals with the local navigation of the aircraft. These first two layers are designed as saturated Proportional Integral Derivative (PID) controllers. For the in-detail development of these loops, see the previous work on [START_REF] Ollervides-Vazquez | Twolayer saturated pid controller for the trajectory tracking of a quadrotor uav[END_REF].

The actual references for the navigation loop come from the topmost layer, which in this case is proposed as a modified LFC protocol whose configuration guarantees coordination between the virtual leader and the states of the UAV, but at the reference level. In this sense, the protocol acts as a reference generator for the navigation layer by coordinating the tracking errors between the aircraft and the virtual leader; this can be defined as

E ζi = e ζj -e ζi Ėζi = ė ζj -ė ζi (11)
being e ζj the tracking error for the j agent, or in this case the virtual leader, e ζi , the tracking error for the UAV, ė ζj and ė ζi its corresponding time derivatives. It is worth to note that having a virtual leader configuration leaves the agent j errors at zero, being an ideal tracking.

For the actual LFC protocol acting as a reference generator, the following structure is proposed

u ζ d i = K S ζ i (t) + K CL e ζi + K ob êζi + K Pi ( N j=1 a ij E ζi ) + K Di ( N j=1 a ij Ėζi ) (12)
where the terms on the first line correspond to a closedloop designed to generate unit compatible references for the navigation layer; K S is a matrix of gains that guarantee that the output is kept when consensus is reached, K CL is a matrix of gains that drives the output to the actual desired states to reach, and K ob is a matrix of gains that assures the generated reference will avoid the present obstacles. The second line of terms is designed to assure the tracking of the time-varying reference of the virtual leader, where K Pi are the gains that assure position consensus while the K Di matrix the gains for the velocity consensus. The introduced term êζi has the task to coordinate the rate at which the protocol change the output taken as navigation reference; and at the same time assures that, in the case of a possible collision, it would drive the output references out of a detection zone. This term is defined as

êζi = ( ζ di (t) + V obi ) -ζi (t) (13)
where V obi is and induced velocity depending on the position and orientation relative to the detection zone, further defined as the contribution of all obstacles present in the field of sensing of the UAV as

V obi = M m=1 ||F m (x i , y i , z i )|| (14)
with m = 1, 2, 3 . . . , M having M as the total number of obstacles, and F m (x i , y i , z i ) as a vector space that defines the detection zone of each obstacle defending on its own euclidean position and the UAV state vector

||F m || = a x e - (x i -bx (t)) 2 2c 2 x a y e - (y i -by (t)) 2 2c 2 y a z e - (z i -bz (t)) 2 2c 2 z ( 15
)
having a x , a y , a z as the magnitude of added velocities, b x (t), b y (t), b z (t) are the coordinates of the obstacle relative to the inertial fixed frame, and c x , c y , c z the geometrical extension of the detection zone in each axis. As an example, Fig. 1, shows the plot for a velocity field with a coverage of 500mm in radius, and placed at the origin. To summarize the complete scheme, in Fig. 2, a graphical representation is shown. Note that the reference the UAV must perform is independent of the other algorithms, so it is possible to have collision scenarios. To deal with the previous situation, the actual references that the UAV must follow, are coordinated by the LFC protocol, which has the task of generate a compatible modified reference, considering the obstacles states, so the navigation controller assures its following, by means of the orientation controller.

IV. VALIDATION

For validation purposes, a set of simulations and an experimental test are carried out. Two scenarios are simulated, a hover flight with a moving obstacle and a circular trajectory tracking with a crossing obstacle. For the case of the experimental test, a circular trajectory with an obstacle in a collision position is tested.

A. Numerical simulations

Numerical simulations were performed in the C++ programming language, using the Euler integration method for solving the ODEs, and at 50Hz with a high priority in the OS scheduler. All algorithms ran in a computer with an Intel Core i7 9750H processor and 24 GB of RAM memory.

The simulated quadrotor is based on a Parrot AR. Drone 2.0, whose parameters and coefficients were proposed by [START_REF] Ollervides-Vazquez | A sectorial fuzzy consensus algorithm for the formation flight of multiple quadrotor unmanned aerial vehicles[END_REF].

1) Hover Flight: For the first numerical simulation scenario, the quadrotor aircraft is commanded with a fixed reference point at the origin, with a reference height of 1500mm. An obstacle with a detection zone of 750mm of radius follows a straight path starting at X = 1750mm and 1500mm, with a velocity of 100mm/s., which allows the obstacle to interact directly with the aircraft, as its detection zone fully covers the original position of the UAV at some point.

In Fig. 3 the resultant trajectory is displayed in the XY plane, as well as the obstacle path and two instants of interaction. Note that at 19s, the detection zone of the obstacle gets close to the UAV, but the last does not move until they start to interact. At 26s, the detection zone is fully covering the original position of the aircraft; however, it is moved to a secure position, as depicted by the blue trajectory; meaning that, while the obstacle is getting close, the quadrotor is gradually moving by the effect of the added velocities. Figure 4 shows the previous behavior in the form of the individual X -axis component. The same can be seen; this is, how the trajectory of the aircraft is disturbed from its reference in order to avoid the obstacle. In Fig. 5, the Y -axis component is presented. Similar behavior is displayed, with the only difference being the sign and magnitude of the movement due to the direction of approximation of the obstacle. Finally, the Z -axis component can be seen in Fig. 6, where no evident changes are presented, due to the fact that the obstacle is at the same height that the aircraft, so these do not interact in this axis. The core of the scheme is the actual references generated by the FLC protocol, which in essence have the task of dealing with external interactions; this is, the obstacles. As shown in Fig. 7, the behavior displayed by the aircraft is a consequence of the references generated by the protocol, whose outputs are the result of the avoidance of the obstacle passing by. Not less important, it is worth to display the trajectory of the obstacle, which for this case is graphed in Fig. 8, accounting for a constant height. As previously stated, both coordinates do not start at the same value; this is a way of preventing an ideal 45 • approximation and having a better evaluation scenario. 2) Trajectory tracking with obstacle intersection: As another test scenario, a circular trajectory is performed. The radius of the circumference trajectory is 1750mm, while the commanded speed is fixed at 0.03rad/s. The trajectory for the obstacle remained at the same speed as the previous scenario, but with equal starting coordinates at 2100mm, in order to assure two interactions with the actual aircraft. In Fig. 9, an XY plane representation of the simulation is shown, where two instants of the whole trajectory are extracted, the first at 11.4s and the second when the time is at 31.6s. It is worth to note that in the first instant, the detection zone is interacting with the ideal reference; thus, the LFC protocol commanded a change of trajectory to the actual aircraft, allowing it to avoid the first encounter. When the quadrotor is at about 3/4 of the cycle, the second interaction happened, and as seen in the Figure, now the LFC protocol commanded an avoidance getting ahead of the obstacle instead of retarding the ideal trajectory. Both interactions can also be seen in Fig. 10, and Fig. 11, for the individual X -axis and Y -axis, noting that at the two instants when the collision is occurred, a sudden deviation from the ideal reference happened, as a result of the LFC protocol changing the navigation reference. Figure 12, displays the behavior in the Z -axis, with null interaction as both the obstacle and the UAV are at equal heights. Note the references trajectories generated by the LFC protocol, which are displayed in Fig. 13. These signals are responsible for the avoidance behavior described by the UAV. The trajectory of the obstacle remains equal in both movement axis, so they overlap in Fig. 14, meaning it moved at the same rates in both directions.

X Position [mm]

B. Real-time experiment

A real-time experimental test is also carried out to further validate the behavior of the whole scheme. A circular trajectory is commanded to the quadrotor aircraft, while a static obstacle is placed directly in the ideal path.

1) Experimental platform: As experimental platform a quadrotor aircraft is used, as well as a motion capture system. The testing facilities are located at the Navigation Laboratory, of the Aeronautical Engineering Research and Innovation Center (CIIIA), of the Autonomous University of Nuevo Leon (UANL).

The motion capture system consists of 16 VICON T40 cameras, streaming position and orientation data to a centralized computer at 100Hz, with an accuracy of 0.1mm and 0.1 • . The quadrotor aircraft used is designed and built for testing purposes in the same facilities test, and it has basic characteristics shown in Table I. Fig. 15 shows the actual aircraft. The embedded avionics consists of a real-time dual core processor with up to 200 Mhz, where one core deals with communications and the remaining with the attitude estimation and controllers execution, an onboard 9-DOF inertial measurement unit used for attitude estimation, and 8 PWM outputs running at 50Hz of sample rate, of which 4 are used as motor outputs. FLC protocol and navigation control are processed in the same central computer, which deals with the position estimation, while in the onboard avionics, the attitude control is processed, using the estimated Euler angles. 2) Trajectory tracking with obstacle intersection: The experimental circumference followed is defined with a radius of 1500mm and an angular velocity of 0.02rad/s, with a fixed height at 1200mm. The obstacle is placed at X = -1425mm and Y = -20mm with a virtual height of 1200mm; this is, although the physical obstacle is short, a grater height is commanded as a security measure. The actual size of the obstacle is 350mm of radius, while the detection zone is established at 750mm.

In Fig. 16 the ideal references and actual trajectory performed by the quadrotor is shown. The obstacle and its detection zone are also displayed as a blue circle and red dotted circle, respectively. It is possible to note that in this case, the aircraft actually crosses through the detection zone; despite this, no collision occurred due to the fact that the detection zone is established with a section that is greater than the actual obstacle. This behavior is expected knowing the fact that the actual quadrotor used for the physical experiment has a slower response time in comparison with the simulated platform. Nonetheless, the performance exhibited is taken as acceptable, as the goal of avoiding the obstacle is met. It is important to note, that the velocity field not only affects in close proximity or within the detection zone; having an exponential decay of the field means that a small residue of velocity is going to be present, more notoriously in the experimental test. Individual axis signals are shown in Fig. 18 for the Xaxis, in Fig. 19 for the Y -axis and for the Z -axis in Fig. 20. Obstacle interaction occurred at approximately at 55s, which can be observed as the disturbance in the signal in comparison with the ideal reference. This experimental test was captured on video and is available on the following link https://youtu.be/ IiYC9jNboEQ.

V. CONCLUSIONS

In this work, a navigation scheme is presented, based on an implementation of the Leader-follower consensus, having the ideal reference as a virtual leader. In addition, the LFC is designed with awareness of obstacles, defining a virtual velocity field, so the consensus condition is changed during the fly. Numerical simulations are performed to validate the effectiveness of the approach, as well as an experimental scenario. Acceptable results are obtained, where the main goal is achieved, considering the limitations of the aircraft used and the necessary detection zones and time response.

Fig. 1 :

 1 Fig.1: Representation of the directions of a added velocities field with 500 radius, positioned at the origin.

Fig. 2 :

 2 Fig. 2: Graphical representation of the complete navigation scheme.

Fig. 3 :

 3 Fig. 3: Trajectories followed by the UAV and the obstacle in the hover flight scenario.

Fig. 4 :

 4 Fig.4: X-axis trajectory followed by the UAV in the hover flight scenario.

Fig. 5 :

 5 Fig.5: Y-axis trajectory followed by the UAV in the hover flight scenario.

Fig. 6 :

 6 Fig.6: Z-axis trajectory followed by the UAV in the hover flight scenario.

Fig. 7 :

 7 Fig. 7: Generated consensus references in the hover flight scenario.

Fig. 8 :

 8 Fig. 8: Trajectories of the obstacle in the hover flight scenario.

Fig. 9 :

 9 Fig. 9: Trajectories followed by the UAV and the obstacle in the trajectory tracking flight scenario.

Fig. 10 :Fig. 11 :

 1011 Fig. 10: X-axis trajectory followed by the UAV in the trajectory tracking flight scenario.

Fig. 12 :

 12 Fig.12: Z-axis trajectory followed by the UAV in the trajectory tracking flight scenario.

Fig. 13 :

 13 Fig. 13: Generated consensus references in the trajectory tracking flight scenario.

Fig. 14 :

 14 Fig. 14: Trajectories of the obstacle in the trajectory tracking flight scenario.

Fig. 15 :

 15 Fig. 15: Quadrotor used for the experimental tests.

Fig. 16 :

 16 Fig. 16: Trajectories followed by the UAV and the obstacle in the experimental trajectory tracking flight scenario.

Fig. 17

 17 Fig. 17: 3D trajectories followed by the UAV and the obstacle in the experimental trajectory tracking flight scenario.

  Fig. 17: 3D trajectories followed by the UAV and the obstacle in the experimental trajectory tracking flight scenario.

Fig. 18 :Fig. 19 :Fig. 20 :Fig. 21 :

 18192021 Fig. 18: X-axis trajectory followed by the UAV in the experimental trajectory tracking flight scenario.

TABLE I :

 I Quadrotor specifications and parameters.

ACKNOWLEDGMENT

This research work is supported by the Office of Naval Research Global through the grant number N62909-20-1-2030