Influence of Forward and Descent Flight on Quadrotor Dynamics

The focus of this chapter is an aircraft propelled with four rotors, called the quadrotor. Quadrotor was among the first rotorcrafts ever built. The first successful quadrotor flight was recorded in 1921, when De Bothezat Quadrotor remained airborne for two minutes and 45 seconds. Later he perfected his design, which was then powered by 180-horse power engine and was capable of carrying 3 passengers on limited altitudes. Quadrotor rotorcrafts actually preceded the more common helicopters, but were later replaced by them because of very sophisticated control requirements Gessow & Myers (1952). At the moment, quadrotors are mostly designed as small or micro aircrafts capable of carrying only surveillance equipment. In the future, however, some designs, like Bell Boeing Quad TiltRotor, are being planned for heavy lift operations Anderson (1981); Warwick (2007).


Introduction
The focus of this chapter is an aircraft propelled with four rotors, called the quadrotor. Quadrotor was among the first rotorcrafts ever built. The first successful quadrotor flight was recorded in 1921, when De Bothezat Quadrotor remained airborne for two minutes and 45 seconds. Later he perfected his design, which was then powered by 180-horse power engine and was capable of carrying 3 passengers on limited altitudes. Quadrotor rotorcrafts actually preceded the more common helicopters, but were later replaced by them because of very sophisticated control requirements Gessow & Myers (1952). At the moment, quadrotors are mostly designed as small or micro aircrafts capable of carrying only surveillance equipment. In the future, however, some designs, like Bell Boeing Quad TiltRotor, are being planned for heavy lift operations Anderson (1981); Warwick (2007).
In the last couple of years, quadrotor aircrafts have been a subject of extensive research in the field of autonomous control systems. This is mostly because of their small size, which prevents them to carry any passengers. Various control algorithms, both for stabilization and control, have been proposed. The authors in Bouabdallah et al. (2004) synthesized and compared PID and LQ controllers used for stabilization of a similar aircraft. They have concluded that classical PID controllers achieve more robust results. In Adigbli et al. (2007); Bouabdallah & Siegwart (2005) "Backstepping" and "Sliding-mode" control techniques are compared. The research presented in Adigbli et al. (2007) shows how PID controllers cannot be used as effective set point tracking controller. Fuzzy based controller is presented in Varga & Bogdan (2009). This controller exhibits good tracking results for simple, predefined trajectories. Each of these control algorithms proved to be successful and energy efficient for a single flying manoeuvre (hovering, liftoff, horizontal flight, etc.).
This chapter examines the behaviour of a quadrotor propulsion system focusing on its limitations (i.e. saturation and dynamic capabilities) and influence that the forward and descent flights have on this propulsion system. A lot of previous research failed to address this practical problem. However, in case of demanding flight trajectories, such as fast forward and descent flight manoeuvres, as well as in the presence of the In Ground Effect, these aerodynamic phenomena could significantly influence quadrotor's dynamics. Authors in Hoffmann et al. (2007) show how control performance can be diminished if aerodynamic effects are not considered. In these situations control signals could drive the propulsion system well within the region of saturation, thus causing undesired or unstable quadrotor behaviour. This effect is especially important in situations where the aircraft is operating at its limits (i.e. carrying heavy load, single engine breakdown, etc.).
The proposed analysis of propulsion system is based on the thin airfoil (blade element) theory combined with the momentum theory Bramwell et al. (2001). The analysis takes into account the important aerodynamic effects, specific to quadrotor construction. As a result, the chapter presents analytical expressions showing how thrust, produced by a small propeller used in quadrotor propulsion system, can be significantly influenced by airflow induced from certain manoeuvres.

Basic dynamic model
This section introduces the basic quadrotor dynamic modeling, which includes rigid body dynamics (i.e. Euler equations), kinematics and static nonlinear rotor thrust equation. This model, based on the first order approximation, has been successfully utilized in various quadrotor control designs so far. Nevertheless, recent shift in Unmanned Aerial Vehicle research community towards more payload oriented missions (i.e. pick and place or mobile manipulation missions) emphasized the need for a more complete dynamic model.

Kinematics
Quadrotor kinematics problem is, actually, a rigid-body attitude representation problem. Rigid-body attitude can be accurately described with a set of 3-by-3 orthogonal matrices. Additionally, the determinant of these matrices has to be one Chaturvedi et al. (2011). Since matrix representation cannot give a clear insight into the exact rigid body pose, attitude is often studied using parameterizations Shuster (1993). Regardless of the choice, every parameterization at some point fails to fully represent rigid body pose. Due to the gimbal lock, Euler angles cannot globally represent rigid body pose, whereas quaternions cannot define it uniquely. Chaturvedi et al. (2011) Although researchers proved the effectiveness of using quaternions in quadrotor control Stingu & Lewis (2009), Euler angles are still the most common way of representing rigid body pose. To uniquely describe quadrotor pose using Euler angles, a composition of 3 elemental rotations is chosen. Following X − Y − Z convention, a world reference coordinate system is first rotated Ψ degrees around X axis. After this, a Θ degree rotation around an intermediate Y axis is applied. Finally, a Φ degree rotation around a newly formed Z axis is applied to yield a transformation matrix from the world coordinate system W to the body frame B, as shown in figure 1. Equations 1 and 2 formalize this procedure: where cφ and sφ stand for cos(φ) and sin(φ), respectively. The same abbreviations are applied to other angles as well.

Dynamic motion equations
Forces and torques, produced from the propulsion system and the surroundings, move and turn the quadrotor. In this paragraph, the quadrotor is viewed as a rigid body with linear and circular momentum, − → L and − → M respectively. According to the 2nd Newtons law, the force applied to the body equals the change of linear momentum. Using the principal of the change of momentum used in Jazar (2010), the following equation maps the change of quadrotor's position with respect to the applied force: where m q represents quadrotor mass and − → v its velocity vector. Due to the fact that most unmanned quadrotors are electrically driven, it is safe to assume that quadrotor mass does not change over time, resulting in a simple equation 3.
Same analysis can be applied to angular momentum, having in mind, the angular momentum is produced from the quadrotor motion as well as from the rotors spinning to produce the desired thrust. There are four important variables concerning angular momentum: quadrotor angular speed vector -− → ω , rotor angular speed vector -− → Ω , quadrotor inertia tensor -I q and rotor inertia tensor -I r . Angular motion equations can be derived as follows: Quadrotors are normally constructed to be completely symmetric. Therefore, their tensor of inertia is a diagonal matrix 5. The same rule applies for rotors as well(otherwise they would be misbalanced and completely useless). Furthermore, rotors spin in one direction only, so the rotor angular speed vector − → Ω has only one component Ω z . Evaluating 3 yields a circular motion equation 6.
Equation 6 calculates rotation speeds in the body frame coordinate system. To transform these body frame angular velocities into world frame rotations, one needs a transformation matrix 8. This matrix is derived from successive elemental transformations 7 similarly as kinematics equation 2. Infinitesimal changes in Euler angles, affect the rotation vector in a way that the first Euler angle Ψ undergoes two additional rotations, the second angle Θ only one additional rotation, and the final Euler angle Φ no additional rotations Jazar (2010):

Rotor forces and torques
Four quadrotor blades are placed in a square shaped form. Blades that are next to each other spin in opposite directions, thus maintaining inherent stability of the aircraft. The same four blades that make the quadrotor hover enable it to move in the desired direction. Therefore, in order for quadrotor to move, it has to be pitched and rolled in the desired direction. To pitch and roll the quadrotor, some blades need to spin faster, while others spin slower. This produces the desired torques, which in term affect aircraft attitude and position Orsag et al. (2010).
Depending on the orientation of the blades, relative to the body coordinate system, there are two basic types of quadrotor configurations: cross and plus configuration shown in figure 2. In the plus configuration, a pair of blades spinning in the same direction, are placed on x and y coordinates of the body frame coordinate system. With this configuration it is easier to control the aircraft, because each move (i.e. x or y direction) requires a controller to disbalance only the speeds of two blades placed on the desired direction. The cross configuration, on the other hand, requires that the blades are placed in each quadrant of the body frame coordinate system. In such a configuration each move requires all four blades to vary their rotation speed. Although the control system seems to be more complex, there is one big advantage to the cross construction. Keeping in mind that the amount of torque needed to rotate the aircraft is very similar for both configurations, it takes less change per blade if all four blades change their speeds. Therefore, when the aircraft carries payload and operates near the point of saturation, it is wiser to use the cross configuration.
Changing the speed of each blade for a small amount, as opposed to changing only two blades but doubling the amount of speed change, will keep the engines safe from saturation point. Basic control sequences of cross configuration are shown in figure 3. First approximation of rotor dynamics implies that rotors produce only the vertical thrust force. As the rotors are displaced from the axis of rotation (i.e. x and y axis) they produce corresponding torques,

Aerodynamics
As the quadrotor research shifts to new research areas (i.e. Mobile manipulation, Aerobatic moves, etc.) Korpela et al. (2011);Mellinger et al. (2010), the need for an elaborate mathematical model arises. The model needs to incorporate a full spectrum of aerodynamic effects that act on the quadrotor during climb, descent and forward flight. To derive a more complete mathematical model of a quadrotor, one needs to start with basic concepts of momentum theory and blade elemental theory.

Combining momentum and blade elemental theory
The momentum theory of a rotor, also known as classical actuator disk theory, combines rotor thrust, induced velocity (i.e. airspeed produced in rotor) and aircraft speed into a single equation. On the other hand, blade elemental theory is used to calculate forces and torques acting on the rotor by studying a small rotor blade element modeled as an airplane wing so that the airfoil theory can be applied. Bramwell et al. (2001) A combination of these two views, macroscopic and microscopic, yields a base ground for a good approximative mathematical model.

Momentum theory
Basic momentum theory offers two solutions, one for each of the two operational states in which the defined rotor slipstream exists. The solutions refer to rotorcraft climb and descent, the so called helicopter and the windmill states. Quadrotor in a combined lateral and vertical move is shown in figure 4. The figure shows the most important airflows viewed in Momentum theory: V z and V xy that are induced by quadrotor's movement, together with the induced speed v i that is produced by the rotors.
Unfortunately, classic momentum theory implies no steady state transition between the helicopter and the windmill states. Experimental results, however, show that this transition exists. In order for momentum theory to comply with experimental results, the augmented momentum theory equation 10 is proposed Gessow & Myers (1952), where V 2 z 7.67 term is introduced to assure that the augmented momentum theory equation complies with experimental results, R stands for rotor radius and ρ is the air density. It is easy to show that in case of autorotation with no forward speed, thrust in equation 10 becomes

Blade element theory
Blade element theory observes a small rotor blade element ∆r 5. Figure 5 shows this infinitesimal part of quadrotor's blade together with elemental lift and drag forces it produces Bramwell et al. (2001). For better clarity angles are drawn larger than they actually are: where C L and C D are lift and drag coefficients, S is the surface of the element and V str the airflow around the blade element. The airflow is mostly produced from the rotor spin ΩR and therefore depends on the distance of each blade element to the center of blade rotation. Adding to this airflow is the total air stream coming from quadrotor's vertical and horizontal movement, V S = V xy + V z . Finally, blade rotation produces additional induced speed v i . The ideal airfoil lift coefficient C L can be calculated using equation 12 Gessow & Myers (1952).
where a is an aerodynamic coefficient, ideally equal to 2π. The effective angle of attack α ef , is the angle between the airflow and the blade. Its value changes with the change of airflow direction and due to the blade twist.
Standard rotor blades are twisted because the dominant airflow coming from blade rotation increases linearly towards the end of the blade. According to equation 11 this causes the increase of lift and drag forces. The difference in forces produced near and far from the center of rotation would cause the blade to twist, and ultimately brake. To avoid that, a linear twist,  (2009) is introduced to the blade design.
The effect of varying airflow can be calculated separating the vertical components V z + v i and horizontal ones V xy + Ωr. The airflow direction angle Φ can be easily calculated from the equation As lift and drag forces are not aligned with body frame of reference, horizontal and vertical projection forces need to be derived. Keeping in mind that Ωr ≫ {V z , v i , V x y} small angle approximations cos(Φ) ≈ 1 and sin(Φ) ≈ Φ can be used. Moreover, in a well balanced rotor blade, drag force should be negligible compared to the lift Gessow & Myers (1952). Applying this considerations to 11 and keeping in mind the relations from figure 5 enables the derivation of horizontal and vertical force equations 15.

Applying blade element theory to quadrotor construction
This section continues with the observation of a small rotor blade element ∆r from the previous section, placing it in real surroundings shown in figure 6. Since the blades rotate, the forces produced by blade elements tend to change both in size and direction. This is the reason why an average elemental thrust of all blade elements should be calculated. Figure 6 shows the relative position of one rotor as it is seen from quadrotor's body frame. This rotor is displaced from the body frame origin and forms an angle of 45 • with quadrotor's body frame x axis. Similar relations can be shown for other rotors. Accounting for the number of rotor blades N, the following equation for rotor vertical thrust force calculation is proposed Orsag & Bogdan (2009): The term inside the brackets of equation 17 is known as a thrust coefficient, and is given separately in 18.
Variables μ,λ i and λ c are speed coefficients V xy RΩ , V z RΩ and v i RΩ respectively. New constantc is the average cord length of the blade element shown in figure 5.
The same approach can be applied for the calculation of horizontal forces and torques produced within the quadrotor Orsag & Bogdan (2009). Calculated lateral force has x and y components, coming both from the drag and lift of the rotor, given in 19.
In case of torque equations the angles between the forces and directions are easily derived from basic geometric relations shown in figure 6, resulting in the elemental torque equations Orsag & Bogdan (2009): Using the same methods which were used for force calculation, the following momentum coefficients were calculated: It is important to notice that equations 20 have two solutions, since the rotors spin in different directions, as seen in figure 3. Different rotational directions have the opposite effect on torques. This is why the ± sign is used in torque equations. These differences, induced from the specific quadrotor construction, along with the augmented momentum equation provide an improved insight to quadrotor aerodynamics. Regardless of the flying state of the quadrotor, by using these equations one can effectively model its behavior.

Building a more realistic rotor model
Building a more realistic rotor model begins with redefining its widely accepted static thrust equation 22 with real experimental results. No matter how precise, static equation is valid only when quadrotor remains stationary (i.e. hover mode). In order for the equation to be valid during quadrotor maneuvers, aerodynamic effects from 3.1 need to be incorporated into the equation.

Experimental results
This section presents the experimental results of a static thrust equation for an example quadrotor. Most of researched quadrotors use DC motors to drive the rotors. Although new designs use brushless DC motors (BLDC), brushed motors are still used due to their lower cost. Some advantages of brushless over brushed DC motors include more torque per weight, more torque per watt (increased efficiency) and increased reliability Sanchez et al.
where Θ 3 4 is a mechanical angle at the 3 4 of the blade length R 13. Θ tw can later be assessed from the blade construction. Rearranging equation 23 yields an equation for solving the mechanical angle problem 24.
Using experimental data from table 1 it is easy to calculate rotor angle Θ 3 4 . For given set of data the average λ i = v i ΩR = 0.0766. Therefore the mechanical angle Θ 3 4 = 11.6291 o , which is well between the expected boundaries.
Obtained data is piecewise linearized, in order to clearly demonstrate the differences between various voltage ranges. From Fig. 7 it can be seen how thrust declines near the point of saturation. This is important to notice, when deriving valid algorithms for quadrotor stabilization and control. Linearizaton coefficients are given in table 2.

Applying aerodynamics to rotor dynamic model
To apply aerodynamic coefficient 18 to the static thrust experimental results, one needs to multiply experimental results with dynamic-to-static aerodynamic coefficient ratio 25. For the calculation of the aerodynamic coefficient C T it is crucial to know three airspeed coefficients μ, λ c and λ i . Two of them, μ, λ c , can easily be obtained from the available motion data V xy , V z and ΩR. λ i however, is very hard to know, because it is impossible to measure the induced velocity v i .
One way to solve this problem is to calculate the induced velocity coefficient λ i from the two aerodynamic principals, momentum and blade element theories. The macroscopic momentum equation 10 and the microscopic blade element equation 17 provide the same rotor thrust using different physical approach: When squared, equation 26 can be easily solved as a quadrinome: The results of solving this quadrinome can be shown in a 3D graph 8. Although equations 27 look straightforward to solve, it still requires a substantial amount of processor capacity. This is why an offline calculation is proposed. This way, the calculated data can be used during simulation without the need for online computation. By using calculated values of the induced velocity, it is easy to calculate the dynamic thrust coefficient from equation 18. The 3D representation of final results is shown in figure 9.
Due to an increase of airflow produced by quadrotor movement, the induced velocity decreases. This can be seen in figure 8. Although both movements tend to increase induced velocity, only the vertical movement decreases the thrust coefficient. As a result, during takeoff the quadrotor looses rotor thrust, but during horizontal movement that same thrust is increased and enables more aggressive maneuvers.

Quadrotor model
A complete quadrotor model, incorporating previously mentioned effects is shown in figure  10. A control input block feeds the voltage signals to calculate statics thrust, which is easily interpolated from the available experimental data, using an interpolation function as shown in figure 7.
Static rotor thrust is applied to equation 25 along with aerodynamic coefficient C T (μ, λ c , λ i ). Induced velocity and aerodynamic coefficient are calculated using inputs from the current flight data (i.e. λ c , μ). This data is supplied from the Quadrotor Dynamics block. The calculation can be done offline, so that a set of data points from figure 9 can be used to A combination of the results provided from these two blocks using equation 25 gives the true aerodynamic rotor thrust. The same procedure is used to calculate the induced speed from the data shown in figure 8. Once the exact induced speed is known it can be applied to horizontal coefficients 19 and torque coefficients 21. In this way, quadrotor dynamics block can calculate quadrotors angular and linear dynamics using equations 6 and 3.
Dynamics data is finally fed into the kinematics block, that calculates quadrotor motion in world coordinate system using transformation matrices 2 and 7.

Conclusion
As the unmanned aerial research community shifts its efforts towards more and more aggressive flying maneuvers as well as mobile manipulation, the need for a more complete aerodynamic quadrotor model, such as the one presented in this chapter arises.
The chapter introduces a nonlinear mathematical model that incorporates aerodynamic effects of forward and vertical flights. A clear insight on how to incorporate these effects to a basic quadrotor model is given. Experimental results of widely used brushed DC motors are presented. The results show negative saturation effects observed when using this type of DC motors, as well as the phenomenon of thrust variations during quadrotor's flight.
The proposed model incorporates aerodynamic effects using offline precalculated data, that can easily be added to existing basic quadrotor model. Furthermore, the model described in the paper can incorporate additional aerodynamic effects like the In Ground Effect.