The Role of Propeller Aerodynamics In the Model of a Quadrotor UAV

Proceedings of the European Control Conference 2009 • Budapest, Hungary, August 23–26, 2009 MoB3.2 The role Roleof of propeller Propeller Aerodynamics in the Model The aerodynamics in model of of aaQuadrotor quadrotor UAV UAV Pierre-Jean Bristeau, Philippe Martin, Erwan Salaün, Nicolas Petit Abstract— We study models of a quadrotor UAV. These models include various aerodynamic effects of the propellers and their interactions with the motion of the rigid body of the UAV. Two main assumptions are formulated: the twisting of the propellers is such that the local angle of attack is constant along the blades in stationary flight, and, secondly, that the local induced velocity is also constant along the blades. Under these assumptions (which correspond to an optimal hovering rotor), we conclude that the flexibility of the propellers plays an important role in the dynamics of the UAV (including the effects of the vertical location of the center of gravity), and that it is interesting to account for it when designing a closed loop controller. I. INTRODUCTION Lately, the quadrotor has emerged as a popular aerial platform for control design experiments among small rotarywing vehicles [1]. It has been the subject of several studies [2], [3], [4], [5], [6]. While numerous teams have reported successful automatic flights [2], [6], [7], [8], this aerial platform, which has a low cost and agile dynamics, is still a subject of a collection of theoretical questions. In particular, the role of several key parameters of its mechanical and aerodynamic design remains to be clearly determined. A prime example is the location of the center of gravity which can be chosen either in the horizontal plane of the rotors, slightly above, or slightly below by selecting the location of various masses (e.g. batteries, sensors, or on-board computer). As we will briefly discuss it, both positions are relevant either for wind gust disturbance or forward flight stability. Others questions focus on the observed aerodynamics drag induced by the rotors during a translational motion. Their analytic expressions play key roles in the quadrotor stability and observability. While the mechanical structure of a quadrotor is simple (four rotors with simple propellers and rigid frame), the dynamic behavior is surprisingly involved. In this paper, we explain how the aerodynamics of the propellers and the motion of the rigid frame can interfere to produce a coupled dynamical model. We formulate several modeling assumptions and expose the derivations of the corresponding models. Their properties are compared at the light of the previously mentioned control-oriented questions. Finally, we discuss the relevance of these models, based on simulations. The paper is organized as follows. This study is based on a simple modeling of the aerodynamics of a propeller, twisted or non-twisted, flexible or not in Section II. Analytic formulas for the generated forces and torques are determined and incorporate the motion of the quadrotor. In Section III, with these torques and forces, we obtain a 6 D.O.F. model with four inputs being the rotation speeds of the four propellers. By successively accounting for the coupling of the propellers dynamics, the rigid body dynamics, and the flexibility of the propellers, the study proposes three models of increasing complexity. In Section IV, we expose analysis results which permit to stress the merits and limitations of the developed models. II. AERODYNAMIC EFFECTS OF A ROTOR MOUNTED ON A RIGID BODY We consider a quadrotor UAV as depicted in Fig. 1. We note u,v,w its translational speed and p,q,r its angular rates in body axes (xb ,yb ,zb ). The body frame is obtained by a rotation of an inertial frame by the Euler angles (ψ,θ ,φ ). Considering a general wind speed which coordinates in the body frame are uw , vw , ww , and the induced velocity vi being a consequence of the rotation of the rotor blades, we use the following notations : ū = uw − u, v̄ = vw − v, w̄ = vi + ww − w. Each motor acts on the body depending on its rotation speed and the lever-arm L which is the length between the rotation axis and the center of gravity (Fig. 1). The mass of the quadrotor is around 1 kilogram. A. Calculus Method Our modeling is widely inspired by works on largesize helicopter rotors [9], [10], [11], and models proposed specifically for quadrotors [6], [7]. We transpose the largesize rotors modeling techniques to small-size rotors, taking into account angular rates, which are negligible at larger scales. By contrast with [12], [13], we do not neglect the forward flight speed (u,v). The aerodynamic effects applied to the rotor are evaluated by integrating, along each rotor blade, the aerodynamic resultant force per surface increment. xb P.J. CAS, Tech, Bristeau is PhD candidate in Mathematics and Control, Unité Mathématiques et Systèmes, MINES Paris60 bvd Saint-Michel, 75272 Paris, FRANCE pierre-jean.bristeau@mines-paristech.fr P. Martin and N. Petit are with the CAS, Unité Mathématiques et Systèmes, MINES ParisTech, 60 bvd Saint-Michel, 75272 Paris, FRANCE E. Salaün is with Decision and Control Laboratory, School of Aerospace Engineering, Georgia Institute of Technology, Atlanta GA 30332-0150, USA ISBN 978-963-311-369-1 © Copyright EUCA 2009 683 c.g. v, ɽ, q yb M1 ʘ1 < 0 zb xb u, ʔ, p zb Fig. 1. w ʗ, w, ʗ r M4 ʘ4 > 0 M3 3 ʘ3 < 0 Notations in the body frame. yb L M2 ʘ2 > 0 P.-J. Bristeau et al.: The Role of Propeller Aerodynamics in the Model of a Quadrotor UAV We assume that, at a current point y along the rotor blade, the resultant force decomposes into a lift force dL and a drag force dD 1 dL = ρU(y)2CLα α(y)cdyzb , (1) 2 1 dD = ρU(y)2 (CD0 +CDi α(y)2 )cdyeθ (2) 2 where ρ is the air density, U is the airspeed, c is the chord, cdy is the surface increment, CLα is the lift coefficient, CD0 is the parasitic drag coefficient, and CDi is the lift-induced drag coefficient. For each point y, we decompose the airspeed (Fig. 2) to calculate the local angle of attack α(y) as a function of the position y along the rotor blade, and of the rotation speed ω of the rotor. Importantly, we also consider the motion of the rigid body of the quadrotor in those variables. After integrating the aerodynamic resultant force with respect to the space variable y, the result is averaged over a period of rotation while assuming the overall motion of the quadrotor is smooth and slow during a period of rotation of the rotor (this is a quasi-stationarity hypothesis). Then, the obtained formulas are simplified considering the following orders of magnitude for the various variables: ω ≈ 4200 tr/min ≈ 440 rad.s−1 , |p|max ≈ |q|max ≈ |r|max ≈ 1 rad.s−1 , |ū|max ≈ |v̄|max ≈ 10 m.s−1 . In summary, using blade element theory, we obtain formulas for the aerodynamic forces and moments applied at the center of the rotor disc. We now present the results. B. Lift Force and Drag Moment By integrating the lift per surface increment, we obtain the lift force LF (or thrust)   α0 w̄ + L(ε1 q − ε2 p) 3 2 − LF = ρcR ω CLα zb (3) 3 2R|ω| where R is the radius of the rotor blade, α0 its pitch angle at rest. The coefficients ε1 and ε2 depend on the rotor under consideration (Mi stands for the ith rotor). In details, ε1 = {1 for M1, 0 for M2, − 1 for M3, 0 for M4}, ε2 = {0 for M1, 1 for M2, 0 for M3, − 1 for M4}. By integrating the drag per surface increment relative to the rotation axis and considering the drag due to lift, we find the drag moment DM  ⎛  ⎞ CD0 2 α0 − 2w̄ +C α Di 0 4 3R|ω|  4  ⎠ zb DM = −sgn(ω)ρcR4 ω 2 ⎝ α0 w̄ w̄ − −CLα R|ω| 3 2R|ω| (4) We have found two expressions. Both feature a quadratic term in the rotation speed, and a linear term with respect to this same variable. Actually, the distribution of the pitch ɲ0 ɲ Urelative wind eɽ UP UT zb Fig. 2. Decomposition of the airspeed to calculate the angle of attack. MoB3.2 angle α0 along the rotorblade is an hyperbolic function of the space variable y. In consequence, the angle of attack α (considering the induced velocity) is assumed constant along the rotorblade. We call αt this constant for all y. vi = αt α(y) = α0 (y) − (5) y|ω| In this paper, from now on, we consider a modeling around stationary flight, and we make the assumption that the induced velocity is uniform and does not vary more than ±1 m.s−1 around its hover value of 4.5 m.s−1 . This is a reasonable assumption as suggested in [10], [2]. Under this hypothesis, we admit that the parameter αt is also independent of the rotation speed since vi is proportional to |ω| around stationary flight. With this twisted rotor blade, the following expressions of the lift force and the drag moment hold   L αt − (ε1 q − ε2 p) zb , LF = ρcR3 ω 2CLα (6) 3 2R|ω|   CD0 +CDi αt2 CLα αt w̄ − DM = −sgn(ω)ρcR4 ω 2 zb (7) 4 3R|ω| C. Drag Force By integrating the drag per surface increment, the drag force DF can be decomposed into two parts. On the one hand, there is the parasitical drag DF0 1 DF0 = ρcR2CD0 (|ω| (ūxb + v̄yb ) − sgn(ω)w̄ (pxb + qyb )) 2 (8) On the other hand, there is a term coming from induced drag DFi   w̄ ω 3 (ūxb + v̄yb ) + (pxb + qyb ) (9) DFi = ρcR CDi αt R2 3 And finally, we have to consider the component of lift in the rotation plane. We call it "drag due to lift" and note it DFL Rαt w̄ (ūxb + v̄y   b) R3 αt ω 2 + sgn(ω)R w̄ − 3 (pxb + qyb ) (10) One can notice the existence of a term proportional to the speed u,v. Interestingly, unlike the usual quadratic drag of a body, this linear dependence can be used to estimate the translational speed using accelerometers measurements during slow forward flight. Conversely, if the interactions of the vehicle dynamics onto the aerodynamic effects are neglected, the drag force caused by the rotors is zero. 1 DFL = − ρcCLα 2 D. Lift Moment Lastly, by integrating the lift per surface increment relative to the rotation axis, we obtain the lift moment   |ω| w̄sgn(ω) 4 (pxb + qyb ) + (ūxb + v̄yb ) LM = ρcR CLα 8 4R2 (11) In this expression, a side effect on the vehicle can be observed (last term in (11)). In particular, a rotor facing a non-zero airspeed along xb is subjected to flip and roll. 684 Proceedings of the European Control Conference 2009 • Budapest, Hungary, August 23–26, 2009 MoB3.2 E. Consideration of flapping dynamics A. Modeling without interactions In fact, these expressions can be rendered more precise by adding a degree of freedom to the rotor dynamics. We take into account the flexibility of the rotor blade, and, implicitly, we allow it to rotate out of the plane normal to the rotation axis. The angle between the blade and this plane is called flapping angle β . The flapping dynamics of the rotor blade can be approximately determined using the Fourier expansion of β , and the conservation of angular momentum around the flapping axis. The flapping angle is determined by the equilibrium between aerodynamic moment, centrifugal moment and stiffness moment. Additionally, the flapping dynamics can be considered, but in this paper, we simply use the expressions of stabilized flapping angles (this is another quasi-stationarity assumption). In details, one has: First, we examine a model without interactions of the motion onto the aerodynamics effects. The linearization around stationary flight makes the aerodynamic effort on the body vanish. To help the interested reader, we put stars to show where these terms would appear. In details, one has ⎤ ⎡ ∗ 0 0 0 ∗ 0 0 −g 0 ⎢0 ∗ 0 ∗ 0 0 g 0 0  ⎢ ⎢0 0 ∗ 0 0 0 0 0 0  ⎢ ⎢0 ∗ 0 ∗ 0 0 0 0 0  ⎢  A=⎢ ⎢∗ 0 0 0 ∗ 0 0 0 0 ⎢0 0 0 0 0 ∗ 0 0 0  ⎢ ⎢0 0 0 1 0 0 0 0 0  ⎢ ⎣0 0 0 0 1 0 0 0 0⎦ 0 0 0 0 0 1 0 0 0 a = (Kr q + Krq p + Kv ū + Kvq v̄)/Kd b = ( − Kr p + Krq q + Kv v̄ − Kvq ū)/Kd (12) (13) ⎡ 0 ⎢ 0 ⎢ ⎢−Kthr ⎢ ⎢ 0 ⎢ B=⎢ ⎢ Kcmd ⎢ Kyaw ⎢ ⎢ 0 ⎢ ⎣ 0 0 with Kd = γω 2 + kβ2 /(γIβ2 ω 2 ), γ = ρc|CLα |R4 /(8Iβ ), (14) Kr = 2|ω| + kβ /(Iβ |ω|), (15) Krq = sgn(ω) γ|ω| − 2kβ /(γIβ |ω|) , (16) Kv = 2γ w̄/R2 , Kvq = sgn(ω)2kβ w̄/(Iβ R2 ω 2 ), (17) where a (respectively b) is the flapping angle along the xb axis (respectively the yb axis), γ is the Lock number of the blade (ratio between aerodynamics effects and inertial effects), kβ its stiffness and Iβ a moment of inertia on the flapping axis. Thus, we find new expressions of the lift effects on a rotor ⎡ ⎤  −a  αt L − (ε1 q − ε2 p) ⎣−b⎦ , (18) LF = ρcR3 ω 2CLα 3 2R|ω| 1  T LM = kβ b −a 0 (19) Naturally, the expressions (6,11) for the rigid rotor are the limit of the expressions (18,19) when the stiffness kβ approaches ∞. III. DYNAMIC COUPLING BETWEEN THE VEHICLE AND THE ROTORS In this section, we consider the whole vehicle with its four rotors. The UAV is assumed to fly out of ground effects. Applying Newton’s second law, the efforts created by the four rotors are incorporated into the rigid body dynamics. To easily observe the impact of each modeling assumption, the equation will be linearized around stationary flight, i.e. |ωi | = ω0 (nominal speed). Let X denote the vector state and U the command vector. The linearized state model writes under the form Ẋ = AX + BU with X = [u v w p q r φ θ ψ]T and U = [δ1 δ2 δ3 δ4 ]T where δi is a small variation of rotation speed of the motor Mi around ω0 (we could discuss the transfer between motor current and rotation speed). By convention, a positive value for δi means that the motor Mi thrusts more. 0 0 −Kthr −Kcmd 0 −Kyaw 0 0 0 0 0 −Kthr 0 −Kcmd Kyaw 0 0 0 ⎤ 0 0   −Kthr   Kcmd   0   −Kyaw   0   0 ⎦ 0 where Kthr , Kcmd , Kyaw are positive coefficients and g represents gravity (g = 9.81 m.s−1 ). This model emphasizes a general lack of observability of the motion of the UAV, since aerodynamic effects are quadratic and the matrix A is quite sparse. As we have assumed the inertia matrix to be diagonal, the dynamics of the three axis are decoupled. Therefore, one can study the dynamics on each axis separately. Precisely, we now focus on the pitch dynamics. The roll dynamics is similar to it up to a rotation, whereas the yaw and vertical dynamics are less critical. We use a reduced state vector X = [u q θ ]T and the control vector U is reduced to the pitch torque U = [δ1 − δ3 ]. B. Rigid modeling Now, we take in account the coupling between the rigid body motion and the aerodynamic effects. Here, rigid rotor blades are considered. Consequently, we find new terms in the (reduced) matrix A and the (reduced) vector B appear. ⎤ ⎡ ⎤ ⎡ 0 0 −g −Kd f 0⎦ (20) B = ⎣Kcmd ⎦ A = ⎣−Kdm −Kl f 0 0 1 0 The coefficient Kd f is positive and stems from the drag force (8-10) which is summed over the four motors. Velocity terms sum up, but angular rates terms cancel out at first order. The coefficient Kl f is positive. It represents the fact that each rotor counteract the vehicle rotation. This phenomenon also appears in the lift effects (6,11). Finally, the coefficient Kdm has the same sign as h, which is the height of the center of gravity relative to the rotors plane (h being positive means that the center of gravity is above the rotors). It stems from 685 P.-J. Bristeau et al.: The Role of Propeller Aerodynamics in the Model of a Quadrotor UAV the moment exerted at the center of gravity by the drag force (8-10). Typical numerical values for these parameters are Kd f ≈ 0.049, Kl f ≈ 4.1, Kdm ≈ 3.8h. Compared to Section III-A, one can notice that some coefficients of the diagonal of the matrix A have been enlarged. Interestingly, these negative terms have a stabilizing effect.The coefficient Kdm will be discussed later on, in Section IV-A. It can also be noticed that the advantage of the quadrotor, compared to the helicopter, is to cancel the side effect visible in (11) by associating contrarotative rotors. C. Flexible modeling At last, we incorporate the flexibility of the rotor blades about stationary flight. Consequently, the matrix A is updated as follows while the vector B remains identical ⎤ ⎡ −K1 K3 −g (21) A = ⎣ K4 −K2 0 ⎦ 0 1 0 We have obtained four new positive coefficients. First, K1 reinforces Kd f with a more important term coming from the tilt phenomenon of the rotor discs (18). This phenomenon takes place as follows. During a forward flight, the advancing blade experiments stronger lift than the retreating one. Under the rigid modeling assumption, one would obtain roll moments on each rotor which would cancel on the whole vehicle. By contrast, by considering the flexibility, stronger lifts do not cause any moment but induces flapping speeds. The flapping angle has the same extrema for all the rotation direction (Fig. 3). In fact, one has to notice that, as well as in the rigid model in Section III-B, on a quadrotor, axial effects are added while side effects cancel (12-17). In the same manner, K2 contributes to incrementing Kl f thanks to the same phenomenon (19), stemming this time from the angular rates. The coefficient K3 represents the inertia of rotor discs which leads to tilt the lift (18). Finally, the coefficient K4 comes from the stiffness of the blades which tends to keep rotor discs orthogonal to the rotation axis of the motors (19). Experimental measurements of the blade stiffness yield K1 ≈ 0.072, K2 ≈ 5.6, K3 ≈ 0.079, K4 ≈ 0.41 − 3.8h. Climbing blade xb Climb speed goes from positive to negative u yb Advancing blade ʘ<0 Retreating blade Climb speed goes from negative to positive MoB3.2 In summary, the flexibility of the propellers results into larger coefficients into the matrix A and stresses a stronger coupling phenomenon between speed and angular rates. IV. DISCUSSIONS A. The role of the location of the center of gravity A common thought, which refers to an analogy between a quadrotor and the celebrated inverted pendulum dynamical system, is that the location of the center of gravity must play an important role in the overall stability of the quadrotor. Due to physical limitations, actual prototypes (both commercial products and academic platforms) consider only center of gravity slightly distant from the "equatorial plane" defined by the rotors. Yet, it remains unclear whether this can improve stability. In fact, two properties must be distinguished [7]. The first one is forward flight stability and the second one is wind gust rejection. In forward flight, the induced wind generates (at the center of gravity) a pitch torque and an aerodynamic force. Both can be calculated from (8-10). Depending on the vertical location of the center of gravity, the torque (which is the dominant effect) can either tend to increase, or decrease, the translational speed of the UAV. This point is pictured in Fig. 4. In details, the torque which depends on the translational velocity u will drive the attitude dynamics θ . If the center of gravity is below the equatorial plane, the absolute value of the pitch angle will decrease and this will damp the system. Otherwise, the absolute value of the pitch angle will increase, and the system will go unstable. Interestingly, the reasoning is symmetric for wind gust disturbances. The situation is pictured in Fig. 5. While the generated effects are similar to the previous case, for stability, it is desired that the attitude of the vehicle counteracts the wind gust. In details, to avoid that the translational velocity follows the wind gust, the UAV should have a negative pitch angle. Having the center of gravity above the equatorial plane generates a contra rotating torque which contributes to this desired effect. It is certainly not a sufficient condition for stability, but the situation would be much worst if the torque would have the opposite sign, i.e. if the center of gravity would be below the equatorial plane. Fig. 4 and Fig. 5 give insight into the "static stability" (as defined in [9], [10], [11]). Further, the situation can also be seen on transfer functions. In fact, the location of the center of gravity plays a role on the poles of the transfer function between U and X and on the zeros of the transfer function between the disturbance ū and X. In particular, the sign of Unstable bl forward f d flight fl h Falling blade Lift forces Falling blade xb Climb speed goes from positive to negative u yb Retreating blade ʘ>0 Ad Advancing i blade bl d Climb speed goes from negative to positive c.g. Drag forces u>0 cg c.g. Stable forward flight Climbing blade Fig. 3. Tilt phenomenon in case of forward flight: counter-clockwise rotation (top) and clockwise rotation (bottom). Fig. 4. flight. 686 Impact of the location of the center of gravity during forward Proceedings of the European Control Conference 2009 • Budapest, Hungary, August 23–26, 2009 MoB3.2 Disturbance b rejection 2 Lift forces 1 Real part of the eigenvalues c.g. uw < 0 Drag forces c.g. Glide with the wind gust 0 -1 -2 Rigid Modeling Flexible Modeling -3 -4 -5 Fig. 5. Impact of the location of the center of gravity on disturbance rejection. K4 determines the sign of the zeros of the transfer between ū and the states. Under the rigid modeling assumption, a (large) positive value of h induces a non minimum phase transfer u K1 s2 + (K1 K2 − K3 K4 )s + gK4 (s) = ū s(s + K1 )(s + K2 ) + K4 (g − sK3 ) q −K4 s2 (s) = ū s(s + K1 )(s + K2 ) + K4 (g − sK3 ) θ −K4 s (s) = ū s(s + K1 )(s + K2 ) + K4 (g − sK3 ) (22) (23) (24) These discussions show a necessary trade-off between stability and disturbance rejection. It is interesting to know in which configuration the quadrotor is, depending on the location of the center of gravity. We saw in Section III that, in the rigid model, the motion speed exerts a moment which value is proportional to the height of the center of gravity. Considering the flexible modeling, the torque exerted by the tilt phenomenon adds an offset to this value. This implicitly appears in (19), and is summarized by the affine expression of the coefficient K4 . The pole-zero map of the two models (Fig. 6), obtained for varying locations of the center of gravity, supports these conclusions. With the rigid modeling of Section III-B, the system always has a stable node, while an unstable focus becomes a pair of stable/unstable nodes as the center of gravity goes from below the rotor plane to above it. With the flexible modeling of Section III-C, there is no remarkable change on the pole-zero map apart from a faster stable node. It is instructive to look at the real part of the eigenvalues, as function of the height of the center of gravity (Fig. 7). The flexible modeling changes the switching value of the center of gravity where the focus split in two stable/unstable nodes. In other words, it impacts on the bifurcation diagram in Figure 7 by sliding it to the right. Typically, the critical point is reached when the center Pole-zero map 1.5 0.84 3 2 0.91 0.84 0.68 0.4 Rigid Modeling 0.998 6 5 4 1 0.998 -0.5 0.99 -1 -1.5 0.91 0.99 0.5 0 0.955 0.976 1 0.976 0.955 -6 -5 -4 -3 0.68 -2 0.4 -1 0 1 2 Pole-zero map 1.5 0 0.91 0.84 3 2 0.91 0.84 0.68 0.4 Flexible Modeling 0.99 0.998 6 5 4 1 0.998 -0.5 0.99 -1 -1.5 0.955 0.976 1 0.5 0.976 -6 0.955 -5 -4 Fig. 6. -3 -2 0.68 -1 0.4 0 Pole-zero map. 1 2 -6 -20 -15 -10 -5 0 5 10 Height of the center of gravity (cm) 15 20 Fig. 7. Real part of the eigenvalues depending as a function of the location of the center of gravity. of gravity is located nearly 10 centimeters above the rotor plane. When the flexibility of the propellers is accounted for, we conclude that, when the center of gravity is close to the rotor plane (|h| < 5 cm), we face a slow unstable focus and the precise location of the center of gravity does not matter. The location of the center of gravity is not used to avoid an unstable node. By contrast, the rigid model features a bifurcation when h is near zero which is inconsistent with every reported experimental flights. We conclude here that, because of the flexibility of the propellers, the location of the center of gravity about the equatorial plane does not play any critical role on the dynamics of the quadrotor. The above reasoning is interestingly complemented by the following discussion. In [3], [4], a quadrotor has been proposed where flapping dynamics is taken in account, but where the rotors have the particularity to be connected to their motors through teetering hubs. In this case, the coefficient K4 can be neglected and only Kdm needs to be considered. Its sign can be determined by the location of the center of gravity. The particularity of teetering hubs is to generate stronger drag forces due to tilt phenomenon without undergoing the corresponding torques. In this case, the dynamics are indeed qualitatively depending on the location of the center of gravity. B. Comparison of the models From a closed-loop control perspective, it is very instructive to compare the frequency domain representation of the previously considered modeling assumptions. In Fig. 8, the Bode diagram of the model (21), which accounts for the rigid body and the flexibility of the propellers, is represented. The input signal is the pitch torque, δ3 − δ1 and the outputs are u, q and θ . To evaluate the impact of the above mentioned hypothesis, we compare the frequency response for various modelings derived from the flexible modeling. It appears that the interactions of the angular velocities on the flexible modeling are negligible: considering (0, Kl f ) or (K3 , K2 ) yields relatively similar results (see Fig. 9). Similarly, the drag force created by the rotor-disc tilt does not impact drastically the Bode diagram. The only coefficient which significantly varies between the rigid and the flexible model is K4 . Neglecting the flexibility moments yields an important error on the whole dynamics in the low-frequency domain. Therefore, it would be prejudicial to consider the rigid modeling for control design purposes since it induces a large 687 P.-J. Bristeau et al.: The Role of Propeller Aerodynamics in the Model of a Quadrotor UAV Flexible Modeling R EFERENCES 50 0 -50 Magnitude of q -100 -2 10 0 -20 -40 -60 -80 -2 10 Magnitude of T Magnitude of u Rigid Modeling 20 0 -20 -40 -60 -80 10 -2 10 -1 10 0 10 10 10 1 10 2 -1 10 0 10 1 10 2 -1 10 Frequency (rad.s-1) 0 10 1 10 2 Fig. 8. Bode diagram (rigid and flexible models). Center of gravity is located 1 cm below the equatorial plane. error for low-frequencies (below 1 rad.s−1 , which includes the typical frequency of the desired closed-loop behavior) and an offset in the location of the resonance frequency (see Fig. 8). Using the flexible model, slight displacements of the center of gravity induce only negligible modifications of the Bode diagram as can be observed in Fig. 9. V. C ONCLUSION The paper exposes the derivation of simple models accounting for the aerodynamics of the propellers on a quadrotor UAV. Using such models can be interesting, in a control perspective, to estimate states variables and to control more accurately such an UAV. Yet, it is certainly possible to control such a system without taking into account the presented involved phenomenons, as has been done for several years. In any case, a conclusion of this study is that, if one is interested in such aerodynamics, one should explicitly consider the flexibility of the propellers or wrong conclusions could be drawn. This is exemplified by the discussion of the location of the center of gravity. In the future, we will focus on UAV design. We will try to use the presented models to determine the structure of a different aerial platform which could fully take advantage of the flexibility of the propellers. Coordinated larger rotors and more flexible propellers could be considered. Acknowledgement The authors wish to thank Johann Forgeard and David Vissière for their technical support and fruitful discussions on the topics of this paper. h = +4, K3 = 0 h = +4 h = -1 50 0 -50 -100 -2 10 Magnitude of T Magnitude of q Magnitude of u MoB3.2 20 0 -20 -40 -60 -80 -2 10 10 -1 10 0 10 -1 10 10 -1 10 Frequency (rad.s-1) 10 1 10 2 0 10 1 10 2 0 10 1 10 2 20 0 -20 -40 -60 -80 10 -2 Fig. 9. Bode diagram (flexible model) for various locations of the center of gravity. [1] L.A.Young, E. Aiken, J. Johnson, R. Demblewski, J. Andrews, and J. Klem, “New concepts and perspectives on micro-rotorcraft and small autonomous rotary-wing vehicles,” in Proceedings of the AIAA Applied Aerodynamics Conference, St. Louis, MO, June 2002. [2] G. Hoffmann, H. Huang, S. Waslander, and C. Tomlin, “Quadrotor helicopter flight dynamics and control: Theory and experiment,” in Proceedings of the AIAA Guidance, Navigation, and Control Conference, Hilton Head, SC, August 2007. [3] P. Pounds, J. Gresham, R. Mahony, J. Robert, and P. Corke, “Towards dynamically favourable quad-rotor aerial robots,” in Proceedings of the Australasian Conference on Robotics and Automation, Canberra, Australia, December 2004. [4] P. Pounds, R. Mahony, and P. Corke, “Modelling and control of a quad-rotor robot,” in Proceedings of the Australasian Conference on Robotics and Automation, Auckland, New-Zealand, December 2006. [5] D. Gurdan, J. Stumpf, M. Achtelik, K.-M. Doth, G. Hirzinger, and D. Rus., “Energy-efficient autonomous four-rotor flying robot controlled at 1 khz,” IEEE International Conference on Robotics and Automation, April 2007. [6] H. Romero, S. Salazar, A. Sanchez, P. Castillo, and R. Lozano, “Modelling and real-time control stabilization of a new VTOL aircraft with eight rotors,” IEEE International Conference on Intelligent Robots and Systems, November 2007. [7] N. Guénard, “Optimisation et implémentation de lois de commande embarquées pour la téléopération intuitive de micro drones aériens x4-flyer,” Ph.D. dissertation, Université Nice Sophia Antipolis, 2007. [8] A. Tayebi and S. McGilvray, “Attitude stabilization of a VTOL quadrotor aircraft,” IEEE Transactions on Control Systems Technology, vol. 14, no. 3, pp. 562–571, 2006. [9] W. Johnson, Helicopter Theory. Princeton, NJ: Princeton University Press, 1980. [10] J. Leishman, Principles of Helicopter Aerodynamics. New York, NY: Cambridge University Press, 2000. [11] R. Prouty, Helicopter Performance, Stability, and Control. Malabar, FL: Krieger Publishing Company Inc., 2003. [12] B. Mettler, Identification Modeling and Characteristics of Miniature Rotorcraft. Boston, MA: Kluwer Academic Publishers, 2003. [13] D. Vissière, P.-J. Bristeau, A. Martin, and N. Petit, “Experimental autonomous flight of a small-scaled helicopter using accurate dynamics model and low-cost sensors,” in Proceedings of the 17th IFAC World Congress, Seoul, Korea, July 2008. A PPENDIX A. Rotor induced drag force The results of the section II can be used to evaluate the Hforce variable which is a standard representation of all linear drag effects. This variable stands for all aerodynamic forces having a non-vanishing first order translational velocity term in their development. Here, it gathers terms coming from the different drag forces and from the rotor-disc tilt. Around stationary flight, we have w̄ ∝ |ω| and Kv /Kd can be approximated by an affine function of |ω|. Thus, the Hforce can be expressed as a function of the rotation speed and the translational speed    |ω| |ω| |ω| H = −λ1 1 + λ2 1 − λ3 (uxb + vyb ) (25) ω0 ω0 ω0 where λ1 ≈ 0.049, λ2 ≈ 0.54 and λ2 ≈ 0.30. Interestingly, this H-force is composed at 70% of drag force and the other 30% come from the tilt of the propellers enabled by their flexibility. In details, this can be observed in (8-10, 18). This flexible effect can be evaluated from the stiffness of the propellers using the equation (12-17). This drag force has a stabilizing effect. It is certainly possible to take advantage of this well-conditioned equation to estimate the longitudinal velocity with a longitudinal accelerometer. 688