Model-based Design

Development and Control of a

Wind Resistant Multirotor UAV

Christian Mnsson
Daniel Stenberg

Department of Automatic Control

MSc Thesis
ISRN LUTFD2/TFRT--5947--SE
ISSN 02805316

Department of Automatic Control

Lund University
Box 118
SE-221 00 LUND
Sweden

2014 by Christian Mnsson and Daniel Stenberg. All rights reserved.

Printed in Sweden by Media-Tryck.
Lund 2014

Abstract
Multirotor UAVs have in recent years become a trend among academics, engineers
and hobbyists alike due to their mechanical simplicity and availability. Commercial
uses range from surveillance to recreational flight with plenty of research being
conducted in regards to design and control.
With applications towards search and rescue missions in mind, the main objective of this thesis work is the development of a mechanical design and control algorithm aimed at maximizing wind resistance. To these ends, an advanced multirotor
simulator, based on helicopter theory, has been developed to give an accurate description of the flight dynamics. Controllers are then designed and tuned to stabilize
the attitude and position of the UAV followed by a discussion regarding disturbance
attenuation.
In order to study the impact of different design setups, the UAV model is constructed so that physical properties can be scaled. Parameter influence is then investigated for a specified wind test using a Design of Experiments methodology. These
results are combined with a concept generation process and evaluated with a control engineering approach. It was concluded that the proposed final design should
incorporate a compact three-armed airframe with six rotors configured coaxially.

Acknowledgments
The authors of this thesis would like to thank their supervisors, Simon Yngve at
Combine Control Systems and Professor Rolf Johansson at the Department of Automatic Control LTH, for their guidance and support.

Contents
List of Figures

List of Symbols and Abbreviations

1.

2.

3.

12

Introduction
1.1 Background . . . . . . . . . .
1.2 Related work . . . . . . . . . .
1.3 Aims of the thesis . . . . . . .
1.4 Autonomous flight regulations

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

18
18
19
19
21

Modeling
2.1 Multirotor dynamics . . . . . . .
Basic flight control maneuvers . .
6-DoF coordinate systems . . . .
Rigid body equations . . . . . . .
2.2 DC motors . . . . . . . . . . . .
2.3 Rotor dynamics . . . . . . . . .
Rotor motion . . . . . . . . . . .
Simplified rotor model . . . . . .
Aerodynamic forces and moments
Induced inflow model . . . . . . .
Rotor model validation . . . . . .
2.4 Airframe dynamics . . . . . . .
2.5 Atmospheric model . . . . . . .
Wind model . . . . . . . . . . . .
Air density . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

22
22
22
24
25
26
29
29
30
32
35
38
40
44
44
47

Control
3.1 Overview . . . . . .
3.2 PID control . . . . .
Cascaded control . .
Attitude control . . .
Position control . . .
3.3 Control allocation .
3.4 Disturbance analysis

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

48
48
50
51
52
57
61
62

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

Contents
4.

Multifactor analysis
4.1 Design of experiments . . . . . . . . . . . . . . . . . . . . . . .
Screening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.

Design development
5.1 Rotor and airframe assessment
Rotor configurations . . . . . .
Airframe configurations . . . . .
5.2 Concept evaluation . . . . . . .
Design criteria . . . . . . . . . .
Results . . . . . . . . . . . . . .
5.3 Proposed design . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

65
66
68
76
84
84
85
100
107
107
114
114

6.

Performance comparison

117

7.

Summary and conclusions

120

8.

Future work

123

Bibliography

124

A.

6-DoF Kinematics

128

B.

Geometric and inertia calculations

B.1 Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 Moment of inertia . . . . . . . . . . . . . . . . . . . . . . . . .
B.3 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

130
131
131
132

C.

Aerodynamic theory
C.1 Momentum theory . .
C.2 Blade element theory
C.3 Blade flapping . . . .
C.4 Coaxial rotors . . . .

134
134
140
144
147

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

List of Figures
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12
2.13
2.14
2.15
2.16
2.17
2.18
2.19
2.20
2.21

Thrust actuation. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Roll actuation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Pitch actuation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Yaw actuation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Roll, pitch and yaw orientation angles. . . . . . . . . . . .
Rotor blade motion. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Incident velocity distribution in hover and forward flight. . . . . . . .
Rotor types generally used for multirotors. . . . . . . . . . . . . . . .
Thrust, horizontal force and drag moment acting on a flapping rotor. .
Definition of the sideslip angle. . . . . . . . . . . . . . . . . . . . . .
Aerodynamic forces and moments expressed in Pi , xw , yw , zw . . . . . .
Schematic overview of the rotor modeling, implemented in Simulink.
Induced velocity solution. . . . . . . . . . . . . . . . . . . . . . . . .
APC 10x4.7 Slowflyer propeller. . . . . . . . . . . . . . . . . . . . .
Wind tunnel test setup. . . . . . . . . . . . . . . . . . . . . . . . . .
Chord length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Blade twist. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Rotor model validation for CT and CQ. . . . . . . . . . . . . . . . .
Body drag force. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Centroid to center of gravity distance. . . . . . . . . . . . . . . . . .
Wind gust. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23
24
24
24
30
31
32
33
33
34
35
37
38
39
39
39
41
43
43
45

3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10

The closed-loop multirotor system. . . . . . . . . .

A generalized control problem. . . . . . . . . . . .
The cascaded PID control structure. . . . . . . . .
Roll/pitch reference and disturbance step responses.

Gcl Bode plot. . . . . . . . . . . . . . . . . . . .

Su Bode plot. . . . . . . . . . . . . . . . . . . . .
Yaw reference and disturbance step response. . . .

Gcl Bode plot. . . . . . . . . . . . . . . . . . . . .

Su Bode plot. . . . . . . . . . . . . . . . . . . . .
Position reference and disturbance step response. .

48
49
51
53
54
54
56
56
57
59

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

List of Figures
3.11 Gxcl Bode plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.12 Sux Bode plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.13 Dryden wind PSD. . . . . . . . . . . . . . . . . . . . . . . . . . . .

59
60
63

4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
4.18
4.19
4.20
4.21
4.22
4.23

COST approach optimization. . . . . . . . . . .

Design of experiments optimization. . . . . . . .
Wind test design. . . . . . . . . . . . . . . . . .
Full factorial design. . . . . . . . . . . . . . . .
Fractional factorial design. . . . . . . . . . . . .
Qualitative center points. . . . . . . . . . . . . .
Linear regression coefficients for emax . . . . . . .
Linear regression coefficients for Pmean . . . . . .
Linear regression validation for emax . . . . . . . .
Linear regression validation for Pmean . . . . . . .
Student chart for emax . . . . . . . . . . . . . . .
Student chart for Pmean . . . . . . . . . . . . . . .
Workspace in modeFRONTIER. . . . . . . . . .
Central composite designs. . . . . . . . . . . . .
Quadratic regression coefficients for emax . . . . .
Quadratic regression coefficients for Pmean . . . . .
Quadratic regression validation for emax . . . . . .
Quadratic regression validation for Pmean . . . . .
Student chart for emax . . . . . . . . . . . . . . .
Student chart for Pmean . . . . . . . . . . . . . . .
Pareto frontier from the optimization run. . . . .
The best designs with respect to wind resistance.
The best designs with respect to efficiency. . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

65
66
68
71
71
72
73
73
74
74
74
75
77
78
79
79
79
79
80
80
81
81
81

5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
5.12
5.13
5.14
5.15

Equilibrium angle explanation. . . . . . . . . . . . . . . . .

Equilibrium angle and power loading for the APC 10x4.7 sf.
Design of experiments for the propeller simulations. . . . . .
Design of experiments result: eq . . . . . . . . . . . . . . .
Design of experiments result: PL . . . . . . . . . . . . . . .
Three-bladed propeller comparison. . . . . . . . . . . . . .
Scaled three-bladed propeller comparison. . . . . . . . . . .
Implementation of a variable pitch propeller. . . . . . . . . .
Coaxial rotor comparison. . . . . . . . . . . . . . . . . . .
Scaled coaxial rotor comparison. . . . . . . . . . . . . . . .
Ducted fan implementation. . . . . . . . . . . . . . . . . . .
Quadrotor: I4 config. . . . . . . . . . . . . . . . . . . . . .
Quadrotor: X4 config. . . . . . . . . . . . . . . . . . . . . .
Quadrotor: H4 config. . . . . . . . . . . . . . . . . . . . . .
Quadrotor: V-Tail config. . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

86
87
87
88
89
90
91
92
94
95
96
100
100
100
101

10

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

List of Figures
5.16
5.17
5.18
5.19
5.20
5.21
5.22
5.23
5.24
5.25
5.26

Tricopter: Y3 config. . . . . . . . . . . . . . . . . . . . . . . . . .
Tricopter: T3 config. . . . . . . . . . . . . . . . . . . . . . . . . .
Hexacopter: X6 config. . . . . . . . . . . . . . . . . . . . . . . . .
Octorotor: I8 config. . . . . . . . . . . . . . . . . . . . . . . . . .
Octorotor: V8 config. . . . . . . . . . . . . . . . . . . . . . . . . .
X4 concept. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Y6 concept. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
X6 concept. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
X4S concept. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Y6 and X4S fusion concept. . . . . . . . . . . . . . . . . . . . . .
Y6 configuration shown in the graphical interface of SimMechanics.

.
.
.
.
.
.
.
.
.
.
.

102
103
104
105
106
109
110
111
112
115
116

6.1
6.2

Wind performance test. . . . . . . . . . . . . . . . . . . . . . . . . .

Wind performance. . . . . . . . . . . . . . . . . . . . . . . . . . . .

118
119

C.1
C.2
C.3
C.4
C.5
C.6

Momentum theory for hover. . . . . . . . . . . . . .

Momentum theory for ascent flight. . . . . . . . . .
Momentum theory for descent flight. . . . . . . . . .
Momentum theory for forward flight [18]. . . . . . .
Aerodynamic environment at a typical blade element.
Momentum theory model for a coaxial rotor in hover.

135
137
138
139
141
147

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

11

List of Symbols and

Abbreviations
Symbols
A
A
A
a0
a1s
Ax
Axy
Ay
Az
b
b1s
c
C
Cd
Cd
CH
Cl
Cl
CP
CQ
CT
d
dCH
dCQ
dCT
dD

12

Rotor disk area [m2 ]

Airframe reference area [m2 ]
Vena contracta area [m2 ]
Coning angle [rad]
Longitudinal flapping angle [rad]
Projected surface area perpendicular to the x-axis [m2 ]
Projected surface area perpendicular to the x- and y-axes [m2 ]
Projected surface area perpendicular to the y-axis [m2 ]
Projected surface area perpendicular to the z-axis [m2 ]
Thrust coefficient [Ns2 ]
Lateral flapping angle [rad]
Chord length [m]
Centroid [m]
Airframe drag coefficient []
Local drag coefficient []
Nondimensional horizontal force coefficient []
Local lift coefficient []
2-D lift-curve-slope of the airfoil section [rad1 ]
Nondimensional power coefficient []
Nondimensional drag moment coefficient []
Nondimensional thrust coefficient []
Drag moment coefficient [Nms2 ]
Infinitesimal horizontal force coefficient []
Infinitesimal drag moment coefficient []
Infinitesimal thrust coefficient []
Infinitesimal drag force [N]

List of Symbols and Abbreviations

dFx
dFz
dH
dL
dm
dQ
dT
e
e
e(t)
emax
Faero
Fd
Fdx , Fdy , Fdz
Fr (s)
fxi , fyi , fzi
Fxi , Fyi , Fzi
Fxy
Fy (s)
g
G
Gcl (s)
h
H
ha
HIB
hc
hm
hs
Hu (s) , Hv (s), Hw (s)
i, j, k
I
I
iw , jw , kw
IXX
Iyb
IYY
IZZ
J
Jp

Infinitesimal force parallel to the rotor disk [N]

Infinitesimal force perpendicular to the rotor disk [N]
Infinitesimal horizontal force [N]
Infinitesimal lift force [N]
Gust length [m]
Infinitesimal drag moment [Nm]
Infinitesimal thrust force [N]
Flapping hinge offset [m]
Distance error [m]
Control error signal
Maximum position error [m]
Input force disturbance [N]
Drag force [N]
Drag force body frame components [N]
Reference to control signal transfer function
Force components expressed in Pi , xw , yw , zw [N]
Force components expressed in Pi , x, y, z [N]
Horizontal force component expressed in Pi , x, y, z [N]
Measurement to control signal transfer function
Gravitational acceleration [m/s2 ]
Center of gravity [m]
Reference closed-loop transfer function
Flight altitude [m]
Horizontal force [N]
Height of the arm [m]
Direction Cosine Matrix
Height of the center [m]
Height of the motor [m]
Height of the shaft [m]
Dryden linear filters
unit vectors for Pi , x, y, z
Norm of the inertia tensor [kgm2 ]
Current [A]
unit vectors for Pi , xw , yw , zw
x-axis moment of inertia [kgm2 ]
Blade moment of inertia around the flapping hinge [kgm2 ]
y-axis moment of inertia [kgm2 ]
z-axis moment of inertia [kgm2 ]
Rotational inertia experienced by the motor [kgm2 ]
Rotational inertia of the propeller [kgm2 ]

13

List of Symbols and Abbreviations

Jr
Js
Kc
Ke
KP , KI , KD
Ks
Kt
k
l
L
la
LIB
Lu , Lv , Lw
m
ma
mc
Md
mp
ms
mtot
Mxi , Myi , Mzi
mxi , myi , mzi
N
N
na
Nb
nr
p, q, r
p, q, r
Pi , x, y, z
Pi , xw , yw , zw
PL
Pmean
pw , qw
q
Q
R
R
r(t)
r0

14

Rotational inertia of the rotor [kgm2 ]

Rotational inertia of the shaft [kgm2 ]
Induced inflow distribution slope []
Back-emf constant [V/(rad/s)]
Proportional, integral and derivative control parameters []
Geometrical scaling factor []
Motor torque constant [Nm/A]
Torsional spring constant [Nm/rad]
Arm length [m]
Armature inductance [H]
Arm length [m]
Rotation Cosine Matrix
Turbulence scale lengths [m]
Mass flow rate [kg/s]
Arm mass [kg]
Center mass [kg]
Drag moment [Nm]
Propeller mass [kg]
Shaft mass [kg]
Total mass [kg]
Moment components expressed in Pi , x, y, z [Nm]
Moment components expressed in Pi , xw , yw , zw [Nm]
Gear box reduction ratio []
Control-loop tuning factor []
Number of arms []
Number of blades []
Number of rotors []
Body-frame rotational velocities [rad/s]
Rotational wind turbulence components [rad/s]
Local coordinate system for rotor i
Local wind coordinate system for rotor i
Power loading [N/W]
Mean power [W]
Airframe angular velocities expressed in Pi , xw , yw , zw [rad/s]
Dynamic pressure [N/m2 ]
Drag moment [Nm]
Rotor radius [m]
Armature resistance []
Reference signal
Distance between the blade root and the rotational axis [m]

List of Symbols and Abbreviations

rc
rd
rdx , rdy , rdz
rGb
rm
rs
S(s)
Su (s)
T
T (s)
t0
tf
Tmax
u, v, w
u, v, w
U
u(t)
U1
U2
U3
U4
UP
UT
v
V
V
Vc
vh
vi
vl
vM
vu
vw
vwg
vwm
vwt
V
w
W
W20

Radius of the center [m]

Centroid to CoG distance vector [m]
Centroids in the x, y and z directions [m]
Flapping hinge and blade CoG distance [m]
Radius of the motor [m]
Radius of the shaft [m]
Sensitivity transfer function
Input sensitivity transfer function
Thrust [N]
Complementary sensitivity transfer function
Wind disturbance start time [s]
Wind disturbance final time [s]
Maximum thrust [N]
Body frame translational velocities [m/s]
Translational wind turbulence components [m/s]
Resultant air flow speed [m/s]
Control signal
Thrust force, control input [N]
Roll moment, control input [Nm]
Pitch moment, control input [Nm]
Yaw moment, control input [Nm]
Vertical air flow speed component [m/s]
Horizontal air flow speed component [m/s]
Fluid onset velocity [m/s]
Voltage [V]
Mean-wind relative airspeed [m/s]
Climb velocity [m/s]
Induced velocity in hover [m/s]
Induced velocity [m/s]
Lower rotor induced velocity in a coaxial configuration [m/s]
Wind gust magnitude [m/s]
Upper rotor induced velocity in a coaxial configuration [m/s]
Wind speed [m/s]
Wind gust speed [m/s]
Wind mean speed [m/s]
Wind turbulence speed [m/s]
Wind magnitude [m/s]
Induced vena contracta velocity [m/s]
Weight force [N]
Specific wind speed 20 ft above the ground [m/s]

15

List of Symbols and Abbreviations

wa
x, y, z
y(t)
z
za
zCG
zm
zPG
zs

SS

0
eq
tw

c
h
i
(i )local

u , v , w

aero
load

, ,
u (),
w ()

16

v (),

Width of the arm [m]

Inertial frame position [m]
Measurement signal
Vertical deflection [m]
Vertical distance for the CoG of the arm [m]
Vertical distance between the CoG and the centroid [m]
Vertical distance for the CoG of the motor [m]
Vertical distance between the CoG and the rotor plane [m]
Vertical distance for the CoG of the shaft [m]
Rotor disk angle of attack [rad]
Blade flapping angle [rad]
Thrust margin []
Sideslip angle [rad]
Convergence error []
Blade twist [rad]
Propeller blade twist at the root of the blade [rad]
Equilibrium angle [rad]
Linear slope of the propeller blade twist [rad]
Measured induced power factor []
Motor damping coefficient [Nm/(rad/s)]
Rotor inflow ratio [rad1 ]
Climb inflow ratio [rad1 ]
Inflow ratio during hover [rad1 ]
Induced inflow ratio [rad1 ]
Locally induced inflow ratio [rad1 ]
Rotor advance ratio [rad1 ]
Air density [kg/m3 ]
Rotor solidity []
Turbulence intensities []
Time constant [s]
Input moment disturbance [Nm]
External motor load torque [Nm]
Induced angle [rad]
Roll, pitch and yaw Euler angles [rad]
Dryden power spectral density functions
Azimuth angle [rad]
Rotor rotational speed [rad/s]

List of Symbols and Abbreviations

Abbreviations
BET
CAD
CCC
CCF
CCI
CFD
COESA
CoG
DC
DCM
DOE
DoF
EASA
EPFL
ETH
GPS
HP
LTH
MOGA
MT
PID
PSD
RSM
TPP
TWS
UAV
UIUC
VRS
VTOL
WBS

Blade Element Theory

Computer-Aided Design
Central Composite Circumscribed
Central Composite Face-centered
Central Composite Inscribed
Computational Fluid Dynamics
U.S. Committee On Extension to the Standard Atmosphere
Center of Gravity
Direct Current
Direction Cosine Matrix
Design Of Experiments
Degree of Freedom
European Aviation Safety Agency
cole Polytechnique Fdrale de Lausanne
Eidgenssische Technische Hochschule
Global Positioning System
Hub Plane
Lunds Tekniska Hgskola
Multi Objective Genetic Algorithm
Momentum Theory
Proportional Integral Derivative
Power Spectral Density
Response Surface Methodology
Tip Path Plane
Turbulent Wake State
Unmanned Aerial Vehicle
University of Illinois at Urbana-Champaign
Vortex Ring State
Vertical Take-Off and Landing
Windmill Brake State

17

1
Introduction
1.1

Background

Small scale unmanned aerial vehicles (UAV) with vertical take-off and landing
(VTOL) capabilities have produced increased interest in recent years. One of the
most common concepts is the multirotor, which, as its name suggests, is an aircraft
based on utilizing several propellers. The multirotor provides a simple platform
to achieve flight control while delivering high levels of maneuverability without
the complex mechanics of a conventional helicopter. Not only are their properties
studied in the commercial and academical spheres, but they have also intrigued the
public with growing communities of hobbyists. There are a number of designs and
configurations available on the market, and whereas some have distinct qualities,
many are simply designed for entertainment.
Normally, a store-bought multirotor will include a radio controller requiring a pilot
on the ground to control the flight path. For more advanced applications and possibly more accurate flight, the aircraft is often made to act in a more autonomous
fashion. It is common to make the UAV adhere to preplanned flight trajectories and
even make logical decisions based on sensor input. Such flight puts extra requirements on the control setup and algorithms as well as the onboard equipment to
ensure safe and steady flight in varying conditions.
Perhaps the largest areas of application of the multirotor is within the photography
and surveillance industries. Possible tasks include surveying of buildings, border
patrol or aerial filming of events, each with specific requirements on design specifications and control performance. An application of great interest to this thesis is
the use of multirotors to assist in search and rescue missions. Outfitting a multirotor with a thermal camera may provide significant aid in locating missing people
and thus becoming a valuable tool for search parties. In fact, Swedish law enforcement are currently investigating the possibility of using UAVs for surveillance and
search applications, indicating a genuine interest for applying UAV technology to
the aforementioned ends [28, 29]. To ensure reliable performance of such multirotor
UAVs, providing a design able to cope with challenging wind conditions is of major
interest and will be the main focus of this thesis.
18

1.2

1.2

Related work

Related work

A number of masters thesis projects concerning multirotors have previously been

conducted at Combine Control Systems AB. The first project was conducted in cooperation with Chalmers University of Technology where Niklas Ohlsson and Martin Sthl presented an autonomous multirotor prototype with the ability to navigate
with computer vision [33]. Since then, Johan Fogelberg has performed a sensor fusion project for GPS-deprived environments [13] and Sara Gustafzelius based her
thesis on navigation with an onboard camera. Another thesis is currently running in
parallel at Combines Gothenburg office where optimal search algorithms are being
developed for multirotors acting as a flock.
A great number of articles and reports regarding multirotors is easily accessible
thanks to the huge interest during recent years. Inspiration has been found in various work from control and aerospace departments at institutions such as EPFL in
Lausanne, ETH in Zrich, Cranfield University, Australian National University and
Lund University. In addition, publications from organizations such as the Institute
of Electrical and Electronics Engineers (IEEE), the American Institute of Aeronautics and Astronautics (AIAA), the American Helicopter Society (AHS) and the National Aeronautics and Space Administration (NASA) have been reliable sources
of information. Dedicated robotics research groups have shown a lot of interest
towards multirotors, contributing with innovative research resulting in numerous
publications. Two of the most prominent are the DAndrea Group, lead by Professor DAndrea at ETH Zrich, and the Kumar Lab Group, lead by Professor Vijay
Kumar at the University of Pennsylvania. Both of these have been sources of fascination and inspiration to the authors.
The extensive modeling in this thesis project has mainly been inspired by the works
of Martnez [20] and Pounds [35], reinforced by helicopter theory based on Leishmans Principles of Helicopter Aerodynamics [18] and Proutys Helicopter Performance, Stability and Control [37]. Control analysis and design has been performed
using methods and conventions presented in Glad and Ljungs Reglerteori - Flervariabla och olinjra metoder. Inspiration on different control strategies has primarily been found in Bouabdallahs extensive work [3]. In the proceeding design work,
the multifactor analysis has been performed in accordance to data analysis company
Umetrics book Design of Experiments: Principles and Applications [11].

1.3

Aims of the thesis

In control analysis and design, it is often the case that a physical process is predetermined with very little ability to influence its core functionality and characteristics.
Any adjustments to the actual plant, with the intent to increase control capability,
are therefore generally limited. Here, however, modeling and control based design
19

Chapter 1. Introduction
development is put to practice with the ability to alter and test physical designs in
a simulated environment. Control and mechanical design can thus be developed in
tandem, aided by simulation and optimization software along with design development methodology.
With application towards search and rescue missions in mind, the main goal of this
thesis work is to create a favorable multirotor design for flight in windy conditions.
This is seen as a major requirement in order to guarantee mission performance and
reliability in a larger range of external conditions. In order to reach this goal, several
milestone objectives must be performed to create an environment where designs can
be evaluated with sufficient accuracy. These objectives represent the main contributions of this thesis work and are summarized in the list below.
Multirotor modeling: The development of a realistic model describing multirotor flight conditions and external phenomena with accuracy is highly important in enabling design evaluation. This includes research of motor dynamics as well as rotor and airframe aerodynamics. Hence, a lot of attention
is paid to this, whereas in other multirotor control projects much simpler models are often sufficient. Implementation of the deduced model in a simulator
environment, where it is possible to configure various dimension and mass
properties, enables analysis of parameter influence. A realistic simulator can
also provide a more realistic testing platform, on which further work concerning control and design can be conducted.
Wind modeling: Modeling of typical wind characteristics to include in the
model simulator is important to add realism to the magnitude and variations
of forces caused by the onset of wind. Using this, suitable wind test cases can
be created to enable evaluation of design and control performance.
Control algorithm development: Developing a suitable control algorithm
to control the multirotor during flight is essential. Since it is presumed that a
search and rescue mission will be performed by trajectory planning, possibly
following a search algorithm, the control architecture must provide positional
control. Also, a control theory analysis is performed on how external disturbances enter the system and what effects they might have.
Multifactor analysis: An extensive evaluation of a generic multirotor is performed by studying the influence of parameters such as geometrical dimensions and mass properties in a multifactor analysis methodology. This is performed by running the simulator with different parameter combinations while
the performance is evaluated.
Design development: The design development concerns the strategy of how
to reach a suitable multirotor design. This is achieved by evaluating current
20

1.4

Autonomous flight regulations

multirotor configurations and assessing their properties and uses for the aims
of this thesis work. The research is broken down into study of airframes and
rotor configurations separately. To reach a final design, a set of performance
criteria is assembled and given a relative weighting according to their rated
importance. Utilizing the conclusions drawn from the multifactor analysis
and design assessments, a few candidate concepts are put together and given
a score for each criteria according to their anticipated performance.
Considering all of the topics of research mentioned above makes the task of multirotor design a multi-faceted one, requiring knowledge and research in a variety of
engineering branches. It should be noted that this thesis work is done at a conceptual level and, although the physical model is aimed at being as realistic as possible, many other realism effects and implementation aspects are neglected (such as
real-time constraints and sensor readings). The final design results and performance
evaluations may be seen as a fore-study to a continuation of multirotor UAV control
and implementation for search and rescue missions.

1.4

Autonomous flight regulations

For search and rescue applications, it is desirable to have multirotors act autonomously according to predetermined flight paths. In order to operate an aircraft
autonomously in Sweden for scientific, commercial or non-recreational purposes,
permission must be granted by the Swedish Transport Agency (Transportstyrelsen)
[43]. Different types of UAVs are categorized into different classes depending on
their take-off weight, maximum kinetic energy and operational premises. Classes
1A (0-1.5 kg) and 1B (1.5-7 kg) are for lighter UAVs, which will most likely be the
weight region of the UAV designed in this thesis. These classes, however, only concern UAVs operated within the line of sight of the operator. Class 3 classification
concerns those UAVs which may operate beyond the line of sight which is likely for
search and rescue missions.
Class 3 UAVs need to be certified by the European Aviation Safety Agency (EASA)
before flight permission may be granted by the Swedish Transport Agency. It is also
possible to apply for special permission for occasional flight or flight in specific
areas. Among a list of things, the operator must be able to override the automatic
control at all times and in the event of communication failure the aircraft must
be disabled safely. A full report on UAV regulations from the Swedish Transport
Agency can be downloaded from their website [43]. In a practical implementation of
UAVs operating in those application areas mentioned in this thesis, these regulations
will have to be considered and adhered to.

21

2
Modeling
The importance of accurate modeling of a multirotor for a variety of flight conditions is seen as critical for a successful design development process. Hence, in
this chapter an advanced multirotor model is presented considering motors, rotors,
airframe and wind dynamics. The modeled multirotor system has six degrees of
freedom, where each type of movement is associated with a number of forces and
moments. Derivation of the aerodynamic forces and moments has been made with
well regarded helicopter rotor theory primarily based on the works of Leishman
[18] and Prouty [37]. An adaption to the small scale multirotor case is inspired by
the extensive thesis work made by Martnez [20].
The model presented in this chapter has been implemented in MATLABs simulation environment to create a multirotor simulator, enabling analysis and evaluation
of control and design. MATLAB is a numerical computation software developed
by Mathworks, providing a block-programming simulation environment known as
Simulink. Within Simulink [23], the Simscape toolboxes SimMechanics and SimElectronics have been used extensively. These environments provide useful tools to
implement mechanical and electronic systems. A graphical 3D environment is also
provided by SimMechanics, where the simulation results can be visualized.

2.1

Multirotor dynamics

Basic flight control maneuvers

In a multirotor UAV, there are four principal means of control actuation to achieve
desired movement thrust, roll, pitch and yaw. Since any aircraft operates in a
six degrees of freedom (6-DoF) space, the system is underactuated. This basically
means that the vehicle will not be able to fully control all of the translational and rotational movements at the same time. It is however, possible to indirectly control the
position by determining the thrust and its direction so that a desired flight trajectory
can be achieved.
Below, the four means of actuation, thrust force and moments, are listed and described. The accompanying figures illustrate the resulting movements for a simple
22

2.1

Multirotor dynamics

quadrotor structure with coordinates fixed to the vehicles body. The circles beside
each i , i = 1, .., 4 indicate the propeller speeds of each rotor, where a larger circle
means higher speed. The force created at each rotor is indicated by a vertical arrow.

Thrust force - U1 .
The thrust force is generated by increasing or decreasing the speed of the
multirotor propellers. The force from each rotor contributes to a total thrust
force directed upwards in an body coordinate frame. For a common quadrotor
configuration, Figure 2.1 shows the applied thrust due to the rotational speed
of the rotors.
Roll moment - U2 .
A rolling moment is achieved by altering the rotational speeds of the rotors
around the body x-axis. Rotors to the left of the x-axis contribute to a positive
moment. The generation of a roll maneuver is shown in Figure 2.2.

Figure 2.1: Thrust actuation.

Figure 2.2: Roll actuation.

Pitch moment - U3
A pitch moment is achieved by altering the rotational speeds of the rotors
around the body y-axis. Rotors to the left of the y-axis contribute to a positive
moment. Figure 2.3 shows how a pitch maneuver is performed.
Yaw moment - U4
The rotation of the rotors encounter a drag force when rotating in the air. This
results in a reaction moment on the airframe directed in the opposite direction of a rotors rotational direction. Figure 2.4 shows how yaw movement is
achieved by altering the propellers individual speeds depending on if they rotate in a clockwise or counter-clockwise direction. For a multirotor to control
its yaw orientation it is important that half of its propellers rotate clockwise
and the other half counter clockwise. Thus, a large part of any unwanted yaw
moment is eliminated for most flight cases.
23

Chapter 2. Modeling

Figure 2.3: Pitch actuation.

Figure 2.4: Yaw actuation.

6-DoF coordinate systems

In flight applications, an aircrafts movement is typically expressed by six basic
types of movement or degrees of freedom (6-DoF). In an inertial (fixed) coordinate
system, these are the three translational movements along the x y z directions
as well as the rotational angles , commonly known as roll, pitch and yaw
angles. These three angles represent the attitude of the aircraft and are defined as
the aircrafts rotated angles around each of the inertial axes. The flight dynamics
conventions presented in this section follow the theory of Robert F. Stengel [41].
It is often convenient to describe an aircrafts dynamics in a body-fixed coordinate
system with the origin placed at the center of gravity. Good reasons for this includes
that the actuation of forces and moments are decoupled, and thus more easily described, in such a coordinate system. Some of the sensor readings are also typically
expressed in the body frame. Also, the moment of inertia tensor does not change
with the body frame orientation since the body coordinates follow the rotational
motion of the aircraft. The orientation of the body frame with respect to the inertial
frame is illustrated in Figure 2.5 below. According to the DIN 9300 standard [1],
the vertical z-axis should be pointing in the downwards direction in aerospace
applications. This convention will be followed in this thesis work.

Figure 2.5: Roll, pitch and yaw orientation angles.

24

2.1

Multirotor dynamics

To express translational and rotational movement in both coordinate systems, transformation matrices between them are used. For translational motion the matrix HIB
in (2.1), often called the direction cosine matrix, is used to transform velocity and
acceleration from the body to the inertial frame. Likewise, angular velocity and
angular acceleration in the body frame can be transformed to equivalences in the
inertial frame using LIB (2.2). For instance, velocity measured in the body frame
vB can be represented in the inertial frame by the transformation vI = HIB vB . In
order to represent a quantity measured in the inertial frame in the body frame, the
inverses of HIB and LIB are used. Reusing the example for the linear velocity, but
instead transforming from the inertial frame, the expression vB = (HIB )1 vI yields
the corresponding velocity in the body frame.

cos cos
HIB = cos sin
sin

cos sin cos + sin sin

cos sin sin sin cos
cos cos
(2.1)

cos tan
sin
(2.2)
cos / cos

sin sin cos sin cos

sin sin sin + cos cos
sin cos

1
LIB = 0
0

sin tan
cos
sin / cos

The direction cosine matrix (2.1) is orthogonal, which makes its inverse equal to its
transpose, (HIB )1 = (HIB )T . It is, however, important to note that (LIB )1 6= (LIB )T .

Rigid body equations

The rigid body kinematics for the multirotor dynamical system will here be described using Newton-Euler formalism in the body coordinate system. For a 6-DoF
rigid body, such as the multirotor aircraft, this takes the generalized form of (2.3) below. In the system description I3x3 is a 3 3 identity matrix, 03x3 a 3 3 containing
only zeros, I is the 3 3 inertia tensor matrix, vB is the body frame translational velocity and B is the body frame rotational velocity. Finally, FB and B are the forces
and moments acting on the body described in the body frame. The derivations made
in this section follow the methods presented by Stengel [41] and Bresciani [6].


I3x3
03x3

03x3
I



 
  
vB
B (mvB )
F
+
= B
B (IB )
B
B

(2.3)

As previously mentioned, it has been assumed that the body coordinate system has
its origin at the center of gravity and that the axes coincide with the principal inertia
25

Chapter 2. Modeling
axes (so that I is diagonal). This results in an easier description of the forces and
moments applied to the system as well as that the inertia tensor I is a diagonal
matrix, which further simplifies the system description. Henceforth, the notations
for the body frame velocities and inertia tensor components are shown in Equations
(2.4) and (2.5).

vB = u

T

IXX
I= 0
0


B = p

0
IYY
0

T
0
0

T

(2.4)

(2.5)

IZZ

Using these notations, the dynamics in (2.3) can be written as a system of equations, which can be seen in (2.6). It can be seen clearly that the system is nonlinear
with coupled terms between the translational and rotational velocities. The expressions for the forces and moments are kept in a generalized form for now, but more
nonlinearities enter with them.

2.2

u =

v =

w =

vr wq

wp ur

p =

q =

r =

IYY IZZ
qr +
IXX
IZZ IXX
pr +
IYY
IXX IYY
pq +
IZZ

uq vp +

1
Fx
m
1
Fy
m
1
Fz
m
1
x

(2.6)

IXX
1
y
IYY
1
z
IZZ

DC motors

The DC motors are the actuators of the multirotor propellers. Rotational motion is
achieved by running a current through the coils where the resulting magnetic field
is cyclically attracted and repelled by the stator magnet. Resulting rotational speeds
are transmitted from the motors to the propellers to generate an upward thrust. There
are a number of constants that define a common DC motors characteristics, which
are listed in Table 2.1 below. It is common to equip DC motors with some kind
26

2.2

DC motors

of transmission that adjusts the torque and the rotational speed of the load. For
multirotors, this can be needed to find a suitable operational range for the specific
propeller.
Parameter
Armature inductance
Armature resistance
Rotational inertia
Damping coefficient
External load torque
Back-emf constant
Motor torque constant
Gear reduction ratio

Notation
L
R
J

load
Ke
Kt
N

Unit
H

kgm2
Nm/(rad/s)
Nm
V/(rad/s)
Nm/A
-

Table 2.1: List of motor parameters.

Motor torque and back-emf constants Ke and Kt associate voltage and current
with output torque and speed of the motor. In SI-units , they take the same numerical
values although they are associated with different phenomena. Henceforth, K will
be used for both (K = Ke = Kt ).
Gear reduction ratio The gear reduction ratio N determines the ratio between
motor and propeller torques and rotational speeds M = NP . It is chosen such
that the propeller operates in a desirable and safe speed range.
Rotational inertia The rotor is connected to the propeller and its shaft through
the motor gear. The rotor rotational inertia Jr around the propeller shaft axis is described by (2.7), where the contributions of the shaft Js and propeller J p are added.
Assuming an idealized motor, with zero internal rotational inertia, the the total inertia J experienced by the motor is described by the right-hand equation in (2.7).

Jr = Js + J p

J=

1
Jr
N

(2.7)

Damping The damping of the DC motor is often very hard to estimate and is
generally very small. Therefore, the coefficient will be approximated to zero in
this thesis work.

Motor dynamics
The mathematical model of a DC motor is obtained using Kirchoffs voltage law
and Newtons second law. The system equations for a general DC motor are shown
27

Chapter 2. Modeling
below in (2.8). This section has largely been inspired by the DC motor theory presented by Javier R. Movellan [31] but include some alterations.

V=
J

dI
+ RI + KM
dt

dM
= KI M load
dt

(2.8)

Considering the expressions derived above, as well as neglecting the very fast transient effects due to the inductance as well as the damping, the resulting equation for
the DC motor is presented in Equation (2.9). Note that the external inertia as well
as the external load torque are divided by N. The load torque for a multirotor is the
drag moment caused by the rotation of the propeller. It has a quadratic dependence
on the rotational speed of the propeller, load = d2P .
V = RI + KM
d
K
K2
1
Jr M = M 3 M + V
N
R
N
R

(2.9)

Rewriting this equation for the propeller speed results in the dynamic equation below (2.10).
KN
K2N2
P d2P +
V
(2.10)
Jr P =
R
R
Steady state analysis When the motors have been running with a constant input
voltage for a short amount of time, the rotational speed of the motor has reached a
steady state. In general, the time constants of the motors are very fast. This is due to
the low electrical inductance and rather low rotational inertia in the kind of motors
and propellers used in multirotor applications. This means that the rotational accelerations M will be zero in the motor system equations (2.8) almost instantaneously
for a given input voltage V . The steady state equations for the DC motors take the
form of (2.11), after some rearranging.
KN 3
K2N3
M
V =0
Rd
Rd
K2N2
KN
2P +
P
V =0
Rd
Rd

2M +

(2.11)

Conversely, the required voltage needed to reach a specific speed can easily be obtained from the polynomials above. Since the motor dynamics are fast, these ex28

2.3

Rotor dynamics

pressions will be used to calculate the input voltage to the motor from a desired
propeller speed when actuating control output in the simulator. The expression used
to obtain a certain steady state propeller speed is given below in (2.12).
V=

Rd 2
+ KNP
KN P

(2.12)

Transient behavior Since the dynamic equation in (2.10) is non-linear, due to the
load torques quadratic dependence on the rotational speed. The time constant of
the motor can therefore not be deduced directly from this equation. For zero load
torque (d = 0) on the motor, the time constant becomes = JR/(K 2 N 2 ), which for
the modeled motor and propeller configuration becomes 0.03 s. Since d will
be relatively small, no significant increase in the time constant will be imminent.
It should be noted that the time constant identified by Fogelberg [13] for a similar
motor as the one modeled here was much larger. It is, however, believed that there
are some unwanted system delays in the previously used hexacopter at Combine.

2.3

Rotor dynamics

The aim of this section is to develop a realistic model of a conventional rotor. Different rotor phenomena will be discussed, and the aerodynamic forces and moments
acting on the rotor will be derived. Rotors can be modeled in many ways and the
choice of method is determined by a compromise between accuracy and effort, depending on which is more important for the task at hand. It was decided to use a
theoretical model for this project because of the possibility to provide information
about how different propeller parameters affect the behavior of the rotor. Additionally, an empirical model would be limited to the specific rotor that would be used
during necessary identification tests. During this project, a combination of blade
element theory (BET) and momentum theory (MT) have been used. These methods
have frequently been used for helicopter studies [18, 37] , and more recently for
multirotor modeling where satisfying results have been achieved by Martnez and
others [20, 2, 3]. There are more complex computational methods available such
as solving the Navier-Stokes equations or using different wake methods, but this
requires a tremendous amount of effort and doesnt fit into the time frame of this
thesis work.

Rotor motion
Aerodynamic forces and moments are generated when a rotor is propelled through
the surrounding air. Thrust is the lift force which is produced by rotating the rotor
blades at a certain velocity . This rotation will also cause a drag moment around
the rotational axis which will oppose the DC motor. In addition to the rotation, the
blades will also experience three different types of motion which are called feathering, blade flapping and lagging, which can be seen in Figure 2.6. These motions
29

Chapter 2. Modeling
affect the flight characteristics of the multirotor and are all dependent on the stiffness of the propeller.

Figure 2.6: Rotor blade motion [18].

Feathering is the angular motion around the longitudinal axis of the blade and is
caused by the aerodynamic moment acting on the blade around this axis. This motion will affect the blades aerodynamic angle of attack which is important for the
flight dynamics of the multirotor. The feathering motion was not deemed important
enough for the objectives of this project, considering the time it would take to model
this accurately.
Blade flapping motion is caused by an uneven incident velocity distribution along
the length of the blade which can be seen to the left in Figure 2.7. This will, in turn,
cause an uneven lift force which is responsible for the flapping motion. To the right
in Figure 2.7, it can be seen that during forward flight, the incident velocity normal
to the propeller blade are no longer constant during a revolution. This means that
the advancing and retreating blades will have non equal flapping angles which will
slightly tilt the direction of the aerodynamic forces and moments.
Lagging motion is caused by the Coriolis force since the distance to the rotational
axis for each blade element is changed when the the blade is flapping and rotating at
the same time. Because of this motion, the center of gravity of the rotor can deviate
from the rotational axis resulting in vibrations. The lagging motion doesnt affect
the control of the multirotor according to Martnez [20] and will therefore not be
modeled in this report.

Simplified rotor model

There are several types of rotors that are designed to accommodate the described
blade motions. For helicopters there are generally four types of rotors [18]:
30

2.3

Rotor dynamics

Figure 2.7: Incident velocity distribution in hover and forward flight [18].
Teetering rotor: The teetering rotor uses a single hinge which is placed along
the rotational axis, to accommodate the flapping and lagging motion. In addition, feathering bearings are used to allow cyclic and collective pitch capability which is needed for the roll and pitch actuation in helicopters.
Articulated rotor: Articulated rotors use separate flapping and lagging
hinges for each blade in addition to the feathering bearings.
Hingeless rotor: Hingeless rotors uses flextures to accommodate the flapping
and lagging in addition to the feathering bearings.
Bearingless rotor: The motion of the blades in a bearingless rotor is solely
controlled by bending, flexing and twisting the hub structure.
The roll and pitch actuation for multirotors are not controlled by cyclic and collective pitch which eliminates the need for feathering bearings. Multirotor rotors are
also much smaller which reduces the structural demands. Multirotors generally use
one of the three rotor types shown in Figure 2.8 where the flexural blade is the most
commonly used option.
However, simulating a flexible blade is complicated and too time consuming for
a project such as this. A simplified model is therefore needed. It was decided to
use rigid rotor blades with spring loaded hinge offsets. The spring stiffness and the
hinge offset are then determined so that the characteristics match that of the real
flexible blade. There have also been a number of other simplifications during this
section such as neglecting the reverse flow region, which can be seen in Figure
2.7, the potential stalling of the propeller blades and the so called tip losses [18].
Additional simplifications will be described in Appendix C.2.
31

Chapter 2. Modeling

Figure 2.8: Rotor types generally used for multirotors [35].

Aerodynamic forces and moments

A rotor mainly generates three types of forces and moments while rotating; the
thrust force, the horizontal force and the drag moment which can be seen in Figure
2.9. The thrust force T is acting perpendicularly to the tip path plane (TPP) whose
boundary is described by the blade tips. The expressions for T can be seen in (2.13).
T = CT AR2 2

(2.13)

The horizontal force H is the rotor drag force which is parallel to the TPP. The
horizontal force is nonexistent during hover and axial flight. The expressions for H
can be seen in (2.14).
H = CH AR2 2

(2.14)

The drag moment Q is also acting perpendicularly to the TPP, but the direction is
determined by the rotation of the rotor. The expressions for Q can be seen in (2.15).
Q = CQ AR3 2

(2.15)

The nondimensional coefficients CT , CH and CQ depends on a number of different variables such as the inflow and advance ratios of the rotor, several parameters
that describe the geometry of the propeller, the flapping angles as well as the induced inflow velocity which will be discussed in Section 2.3. The expressions for
the nondimensional coefficients have originally been developed for helicopters by
Prouty [37] and was later modified by Martnez [20] to make them more suitable
for multirotors. This was accomplished with the blade element theory which is explained in more detail in Appendix C.2.
32

2.3

Rotor dynamics

Figure 2.9: Thrust, horizontal force and drag moment acting on a flapping rotor.
The thrust, horizontal force and the drag moment are either parallel or perpendicular
to the TPP, but ultimately we want to express the forces and moments produced
by the rotor in a coordinate system which is aligned with the body frame of the
multirotor, but has its origin in the resulting center of gravity of the rotor blades.
This coordinate system will be called Pi , x, y, z. It is assumed that the CoG of the
rotor blades remain at a fixed distance from the CoG of the entire multirotor.
To achieve the necessary force and moment transformations it is convenient to introduce yet another coordinate system. This coordinate system will also have its origin
at the center of gravity for the rotor blades, but instead the x-axis will be aligned
with the air velocity component contained in the hub plane (HP) which is perpendicular to the rotational axis. This coordinate system will be called Pi , xw , yw , zw . The
relation between these two newly introduced coordinate systems are defined by the
sideslip angle SS , which can be seen in Figure 2.10.

Figure 2.10: Definition of the sideslip angle.

As defined by Martnez [20], the relationship between the forces and moments,
33

Chapter 2. Modeling
expressed in the different coordinate systems, are defined by the expressions in
(2.16).
Fi = Fxi i + Fyi j + Fzi k = fxi iw + fyi jw fzi kw
Mi = Mxi i + Myi j + Mzi k = mxi iw + myi jw + mzi kw

(2.16)

The F and M variables corresponds to forces and moments expressed in Pi , x, y, z

while f and m, which can be seen in Figure 2.11, corresponds to forces and moments
expressed in Pi , xw , yw , zw .

Figure 2.11: Aerodynamic forces and moments expressed in the coordinate system
which is defined by the air velocity component contained in the HP.
The expressions in (2.17) are then obtained by combining the definition of the
sideslip angle SS with (2.16)

Fxi = fxi cos SSi fyi sin SSi

Mxi = mxi cos SSi + myi sin SSi
Fyi = fxi sin SSi + fyi cos SSi Myi = mxi sin SSi + myi cos SSi

Mzi = mzi
Fzi = fzi

(2.17)

The final step is to transform the thrust, horizontal force and drag moment to the
Pi , xw , yw , zw coordinate system which has been done by Martnez [20]. The relations can be seen in (2.18), where T , H and Q are the absolute values of the thrust,
horizontal force and drag moment. The sign is determined by the rotation of the
rotor:
= + if the rotor is turning counter clockwise.
= if the rotor is turning clockwise.
34

2.3

T tan a1s

fx = H cos a1s + p

1 + tan a1s + tan2 (b1s )

T tan (b1s )
fy = p
2 a + tan2 (b )

1
+
tan
1s
1s

fz = H sin a1s + p
1 + tan2 a1s + tan2 (b1s )

Rotor dynamics

Q tan a1s

mx = p

1 + tan2 a1s + tan2 (b1s )

Q tan (b1s )
my = p
2 a + tan2 (b )

1
+
tan
1s
1s

mz = p
1 + tan2 a1s + tan2 (b1s )
(2.18)

An overview of the rotor modeling can be seen in Figure 2.12, where the theory
from this section has been implemented in Simulink.

Figure 2.12: Schematic overview of the rotor modeling, implemented in Simulink.

Induced inflow model

The aerodynamic forces and moments that are acting on the rotor blades have been
derived in the previous section, but there is still one variable which cannot be determined by the blade element theory alone. During the derivation of the expressions
in (2.13-2.15), Prouty [37] used Glauerts inflow model. This suggests a linear solution for the local induced inflow ratio, which can be seen in (2.19). In this equation,
Kc is a constant that determines the slope of the inflow distribution and i is the average induced velocity ratio, which is typically given by the uniform inflow solution
from the momentum theory [18].


r
(i )local = i 1 + Kc cos
R

(2.19)

The proposed strategy is to combine the two expressions for the thrust coefficient
in (2.20), which have been derived in the Appendixes C.1 and C.2 for MT and BET
respectively.
35

Chapter 2. Modeling




e
2
1
Cl 

2
2

1
+ 0 + (1 + )tw

CT = 4
R
3
2
p

CT = 2( tan ) 2 + 2

= c + i

(2.20)

By doing so, a new equation (2.21) is obtained where the inflow ratio is the only
unknown variable.

=

p


2
8( tan 2 + 2 )
1
2
2
+ 0 + (1 + )tw
3
2
Cl (1 Re )

(2.21)

Equations for the inflow ratio are normally solved numerically and the approach
that was used during this project is inspired by Leishman [18]. It was decided to use
the Newton-Raphson method because of the relatively fast convergence time, despite the necessity of a fairly accurate initial condition. The expression in Equation
(2.22) is solved repeatedly until the error < 0.0005 for which it is assumed that
convergence has been reached.

n+1 = n

f ( )
f 0 ( )


n



 n+1 n 

= 
n+1 

(2.22)

The expressions for f ( ) and f 0 ( ) are derived from (2.21) and can be viewed in
(2.23).

f ( ) =
f 0 ( ) = 1 +

p


1
8( tan 2 + 2 )
2
2
2
+ 0 + (1 + )tw +
3
2
Cl (1 Re )

8( 2 + 2 2 tan )
p
Cl (1 Re ) 2 + 2

(2.23)

The proposed method will determine the induced inflow ratio, but there are unfortunately some major limitations. The reason being that MT is based on the assumption
that there is a well-defined airflow direction. This is not the case when there is an
upward velocity component 2vi V sin 0 through the rotor disk, which typically happens during descent flight. According to Leishman there are four different
types of working states during axial flight [18]:
Normal working state.
36

2.3

Rotor dynamics

Vortex ring state (VRS).

Turbulent wake state (TWS).
Windmill brake state (WBS).

MT is only valid in the normal working state where the air flow is passing down
through the rotor in a well defined slipstream, and in the Windmill Brake State
where the air is flowing in the opposite direction. The air flow is very chaotic in
between these two states and, hence, it is not possible to predict the induced inflow analytically during the vortex ring state or the turbulent wake state. Instead, a
continuous approximation in the form of a quartic curve can be used [18]. This approximation, shown in Equation (2.24), is based on a large number of experiments
from flight tests and wind tunnel measurements, and it is valid for the entire range
of 2vi V sin 0. In this equation, is the measured induced power factor
and k1 = 1.125, k2 = 1.372, k3 = 1.718, k4 = 0.655.
 
 2
 3
 4
vi
Vc
Vc
Vc
Vc
= + k1
+ k2
+ k3
+ k4
vh
vh
vh
vh
vh

(2.24)

The approximated curve can be seen in Figure 2.13 together with the solution produced by the momentum theory. Note that there is a nonphysical solution to the MT
which predicts the induced velocity fairly well until a certain limit where the quartic
function drops of to meet the MT solution for descent flight.

Figure 2.13: Solution for the induced velocity across the four working states [18].
37

Chapter 2. Modeling
In summary, the current working state has to be detected during the simulation. A
suitable method for determining the inflow ratio is then decided i.e. the quartic curve
is used in the range of 2vi V sin 0, while the MT and BET combination is
used for the rest of the flight envelope. Momentum theory is used to test whether
V sin 2vi even if the solution for vi isnt quite theoretically sound in the
specified range.

Rotor model validation

Although the structure of the rotor model is presented in previous section, there
are still a few constants missing whose values depend on the design of the rotor.
To completely validate the proposed model a large number of experiments would
be required, both in axial and forward flight. Because of the time constraints, it
was decided to limit the scope of this thesis to only one type of propeller, the APC
10x4.7 Slowflyer [24] which can be seen in Figure 2.14.

Figure 2.14: Rendered picture of an APC 10x4.7 Slowflyer propeller, modeled in

Solidworks [27].
This propeller has been used for the previous master thesis projects conducted at
Combine Control Systems [33, 13] and its dimensions are suitable for this project
as well. It was also decided to use experimental data for the nondimensional coefficients of thrust CT and power CP CQ to further reduce the effort needed
for this part of the project. This information is available at the UIUC Propeller
Database from the University of Illinois [4]. Unfortunately this means that the rotor model will only be validated for axial flight. The experiments have been conducted for different angular velocities ranging from 420 680 rad/s while the propeller was exposed to winds of varying magnitudes. The corresponding inflow ratios
c = V /(R) were increased until zero thrust was generated, which for the APC
10x4.7 sf translated into roughly 0.12 rad1 . The experimental setup can be seen in
Figure 2.15.
The UIUC Propeller Database [4] also provides geometrical data for the tested propellers. In Figure 2.16, the chord length is plotted as a function of the position along
the length of the blade. Remember that the rotor has been modeled as having a constant chord length. The green line represents the true average of the chord length,
while the red line is the chord length at 70% of the blade length which was used
38

2.3

Rotor dynamics

Figure 2.15: Experimental setup for the wind tunnel test conducted by [5].
during the simulations. According to Martnez [20] it is common practice to use
this value since it will be fairly representative for the whole blade, and thus leads
to the most accurate results. The twist variation along the length of the blade can
be seen in Figure 2.17. The green line represents the linearly fitted function while
the red line represents the modified function which was tuned to get a better overall
model accuracy.

Figure 2.16: Chord length.

Figure 2.17: Blade twist.

There are several more propeller parameters that are needed for the proposed rotor
model. All of the parameters are listed in Table 2.2 together with the values that
were used during the simulations.
ly b and rGb were determined from a CAD model, originally made by Nathan Kau
[17], which is based on the geometrical dimensions provided by [4] as well as the
weight of the propeller. The sectional drag coefficient was tuned based on the experimental data for CQ . The torsional spring constant for the propeller k was deter39

Chapter 2. Modeling
Parameter

Notation

Value

Unit

Chord length
Linear slope of the propeller blade twist
Propeller blade twist at the root of the blade
Rotor radius
Propeller mass
Blade moment of inertia around the flapping hinge
Flapping hinge and blade center of gravity distance
2-D lift-curve-slope of the airfoil section
Sectional drag coefficient
Torsional spring constant

c
tw
0
R
mp
lyb
rGb
Cl
Cd
k

0.0267
0.4
0.55
0.1270
0.0116
2.6275 105
0.0629
2
0.025
1.14

m
rad
rad
m
kg
kgm2
m
rad1
Nm/rad

Table 2.2: Propeller parameters

mined by applying a vertical load F at the tip of a real APC 10x4.7Sf propeller blade
while measuring the vertical deflection z. The spring constant was then calculated
according to 2.25.

k =



FR
Nm
z
rad
arctan
R

(2.25)

Interpolated values based on the experimental data have been divided by the values
calculated by the proposed model to obtain compliance ratios for CT and CQ . These
ratios can be viewed in Figure 2.18.
Conclusions The model has shown relatively satisfying results for the majority of
the tested flight envelope. Some of the propeller parameters such as the linear twist
coefficients had to be slightly tuned to improve the overall accuracy. The highly
transient notches at around 525 and 420 rad/s are caused by irregularities in the
interpolated data as well as a potential decrease in model accuracy for higher inflow
ratios. The model accuracy is also seemingly decreased for rotational velocities
above 650 rad/s which, considering the maximum recommended velocity of 680
rad/s, might be caused by structural phenomena such as increased vibrations. The
validation in this section doesnt cover descent or forward flight, but it is assumed
that the model is accurate enough to base further conclusions upon.

2.4

Airframe dynamics

Any rigid body that moves through a field of viscous fluid, such as gases or liquids,
is affected by a fluid friction often known as drag force. This force is proportional
to the squared velocity v of the onsetting fluid. Calculating the drag force acting
on an arbitrary body is, however, a complicated task. The fluid exerts a pressure on
40

2.4

Airframe dynamics

Figure 2.18: Rotor model validation for CT and CQ.

the exposed reference surface A, which results in a net force on the body scaled by
the drag coefficient Cd . In this thesis, the reference surface is seen as the projected
surface area. To determine the drag force to full accuracy, complete knowledge of
a bodys projected surface area and drag coefficient in all directions must therefore
be known. The general form of the drag equation is shown in (2.26).
1
Fd = Cd Av2
2

(2.26)

The onset of airflow resulting from a multirotors movement through the air as well
as from applied wind are the sources of aerodynamic drag. Realistically, rotational
movement of the aircraft will also cause drag moments, but these will not be modeled here because of the calculation complexity. A relative airspeed V in the body
frame, fully described by its magnitude together with its sideslip angle SS and angle of attack , will create a drag force according to (2.27). In (2.27), Cd is the drag
coefficient and A is the projected surface area in the airspeed direction, which both
depend on the direction of V , i.e. SS and . Since the wind speed is described
in the inertial frame, this velocity must first be transformed using the DCM matrix
(2.1) before being added to the body frame velocity to create the relative airspeed
V .

41

Chapter 2. Modeling
1
Fd = Cd (SS , )A(SS , )V2
2

(2.27)

Using the SS and angles, the force Fd in (2.27) can be described along the three
body axes of the body frame as (2.28).

Fdx = Fd cos SS cos

Fd = Fdy = Fd sin SS cos

Fdz = Fd sin

(2.28)

A multirotor airframe is typically quite complex, which makes it difficult to derive

a general expression for the projected surface area A and the drag coefficient Cd for
an arbitrary direction of V . Multiple arms and other parts partly obscuring each
other depending on different SS and yields heavy geometrical calculations with
many special cases. Considering this, some simplifications to how the projected
surface area A is calculated from an arbitrary point of view is presented in (2.29),
where the airframe may be seen as a rectangular cuboid with three different areas,
Ax , Ay and Az . These areas are calculated in Appendix B and represent the surface
areas in the principal body frame directions. This way, the areas may not be exactly
correct, but at least an approximation is provided where the horizontal and vertical
dimensions of the vehicle are given some representative relevance in regards to the
wind direction. Anyhow, as will soon be discussed, Cd is highly uncertain, which
decreases the importance of an absolutely accurate calculation of the surface area.
A Ax cos SS cos + Ay sin SS cos + Az sin

(2.29)

The drag coefficient Cd is difficult to estimate even in the principal body frame directions x, y and z. For it to be calculated accurately, CFD analysis or wind-tunnel
testing is needed which are tools not available to this project. Even if they were, the
situation is worsened by the fact that the coefficient is greatly affected by the spinning rotors wake turbulence, possibly increasing the drag coefficient. It should also
be noted that the airframe disrupts the inflow to the rotors, which could be viewed
as additional drag. Accurately describing Cd is therefore not deemed possible in this
thesis. However, testing the design for a range of reasonable values [42] gives an
indication of the importance of creating a streamlined airframe. In the simulator, Cd
is assumed to be constant for every airflow direction.
The point of attack of Fd will be in the centroid or the center of geometry C from the
relative airspeeds direction. It is a point calculated in similar fashion as the center
of gravity, but using the projected two-dimensional areas in the airspeed direction
42

2.4

Airframe dynamics

instead. Since this point isnt necessarily located in the center of gravity G, the
force Fd will cause a rotation moment. The situation is described by Figure 2.19,
which shows an arbitrary body from the airflows point of view. The location of C
will be on the surface of the body, whereas G will be inside. The angle d between
Fd and the distance vector rd will yield a distance component perpendicular to
the surface, resulting in a drag moment. Viewing the problem from a cross section
plane containing rd as in Figure 2.20, it is clear that there also is outward distance
component that will not contribute toward any moment.

Figure 2.19: An arbitrary body seen

from the airflows point of view. The
resulting drag force will act in the
centroid C, creating translational and
rotational acceleration.

Figure 2.20: The distance of the centroid in relation to the center of gravity creates a lever for the drag force
Fd resulting in a moment Md .

Since the location of the centroid varies with the airspeed direction, rd = C G will
vary in magnitude and direction depending on SS and . A general description of
the drag moment Md is presented in (2.30).

Md = rd (SS , ) Fd

Md = rd (SS , ) sin d Fd

(2.30)

As previously explained, depending on the airspeed direction the location of the

centroid will differ so C = C(SS , ). To simplify this matter, a similar approximation as for the surface area is performed here for the distance vector rd . This is
then used to determined the moments Md in the simulator. The approximation is
presented in (2.31), where rdx , rdy and rdz are the distance vectors in each of the
body frame principal directions. How the centroid in each direction and the center
of gravity G are calculated is skipped here but is included in Appendix B.

rd rdx cos SS cos + rdy sin SS cos + rdz sin

(2.31)
43

Chapter 2. Modeling
To summarize, a number of simplifications have been made to the drag forces acting
on the airframe. It should also be mentioned that possible lift force of the airframe
has not been considered in the model, but is likely to be very small. The validity
of the airframe model must be rated as relatively poor, especially in regards to the
drag constant Cd . However, the main characteristics of wind will translate to forces
and moments with relevant characteristics on the airframe. Even if the model has
its flaws, it was still decided to include the drag force effect to obtain these characteristics in the simulator. Thus, the benefit of creating a streamlined and compact
airframe may be shown.

2.5

Atmospheric model

For accurate description of external phenomena affecting multirotor flight, atmospheric conditions were modeled. This primarily concerned wind, but also variation
in air density were considered. Since the thesis main aim is to provide a multirotor design able to withstand windy conditions, a lot of attention was put into wind
modeling. Using theory for wind characteristics presented in this section, specific
wind scenarios were created to test different designs ability to resist wind.

Wind model
Wind is modeled as a stationary mean wind from a given direction to which gusts
and turbulence are added according to (2.32). In the equation, the total wind is
denoted vw whereas its mean, gust and turbulence components are vwm , vwg and vwt .
The wind modeling has largely been inspired by the work of Praceus thesis work
[36] and Mathworks content on aerospace modeling [22].

vw = vwm + vwg + vwt

(2.32)

The resulting wind from (2.32) is added to the ground speed of the multirotor to
form the vehicles relative airspeed in the inertial frame. Transforming it to the
body frame yields V , which is used to obtain the aerodynamic forces and moments
affecting the multirotor system.
Wind gust model A wind gust is a temporary increase in the mean wind speed. Its
rate of increase is governed by the present wind speed, the vehicle airspeed and the
increase length of the gust. A gust is commonly modeled to have a 1-cosine shape
according the the US Department of Defense Military Specification MIL-F-8785C
[22]. The equation describing the gust is given by (2.33) where x is the distance into
the gust, dm the gust length and vm the magnitude.
44

2.5

 
v 
x
m
vwg =
1 cos

2
dm

vm

Atmospheric model

x<0
0 x dm

(2.33)

x > dm

Note that Equation (2.33) is not time-dependent and only describes the physical
length of the gust. The duration of which an aircraft is exposed to the relative onset
of a gust depends on the aircrafts current airspeed, where a faster airspeed will
result in a quicker onset of the wind gust. Figure 2.21 shows a typical wind gust
modeled with (2.33).

Figure 2.21: Wind gust with dm = 80 m and vm = 10 m/s.

Dryden wind turbulence The Dryden wind turbulence model is used to predict
wind velocity turbulence, following the Military Specification MIL-F-8785C from
the US Department of Defense [22]. The model takes white noise as an input and
may be seen as a wide-sense stationary stochastic process, i.e. it is time-invariant.
White noise is a completely uncorrelated stochastic process with equal PSD magnitude for the entire frequency range. The power spectral density (PSD) functions
for the Dryden model (2.34) define the contributions towards the turbulence across
the frequency range. The functions describe the spectral density for a given height
h, mean-wind relative airspeed V as well as a specific wind speed 6 m (20 ft) above
the ground W20 .
45

Chapter 2. Modeling
2u2 Lu
1
 
V 1 + Lu
V
 2
2 Lv 1 + 3 Lv V
v () = v 
 2 2
V
1 + Lv
V
 2
2 Lw 1 + 3 Lw V
w () = w 
 2 2
V
1 + Lw
V
u () =

(2.34)

There are also PSD functions describing rotational wind velocity disturbances p, q
and r, resulting from an uneven wind distribution over the vehicles body. Only the
linear components in the u (longitudinal), v (lateral) and w (vertical) directions will,
however, be modeled in this project. This simplification is made because the turbulence magnitudes of these are small and a small scale multirotor will be exposed to
approximately the same wind over its airframe.
The parameters i and Li , where i = u, v, w, in (2.34) are the turbulence intensities and scale lengths. According to the MIL-F-8785C specifications, these can be
calculated as (2.35) for altitudes below 300 m, which will most likely cover the
operating range of the UAVs studied and designed in this thesis.

Lw = h
w = 0.1W20

h
(0.177 + 0.000823)1.2
w
u = v =
(0.177 + 0.000823)0.4
Lu = Lv =

(2.35)

To generate the turbulence, the transfer functions below 2.36 are fed with zero mean
unit variance white noise. The turbulence obtained from the Dryden transfer functions are then added to the constant wind components to generate the total wind.
Transfer functions H(i) are related to their corresponding PSD functions ()
through () = |H(i)|2 . This property has been used to obtain (2.36).

46

2.5

1
Lu
1+ s
V

3Lv
r
Lv 1 + V s
Hv (s) = v


V
Lu 2
1+ s
V

3Lw
r
Lw 1 + V s
Hw (s) = w


V
Lw 2
1+ s
V
Hu (s) = u

Atmospheric model

2Lu
V

(2.36)

Air density
Air density for a given altitude is calculated using the COESA 1976 atmospheric
model included in the Simulink Aerospace toolbox. It has no greater effect on the
model since the quantity varies very slowly, especially when flight is conducted
at a small range of altitudes. More information regarding atmospheric conditions
according to the COESA 1976 standard can be found in [45].

47

3
Control
In order to achieve the desired positional control, a suitable control algorithm must
be implemented in the simulator. The development of this will be the main focus of
this chapter, resulting in a proposed control structure along with a linear analysis of
the closed-loop system. In addition, the controllers ability to reject external disturbances will be discussed. The design and analysis performed in this chapter follows
the conventions and methods presented by Glad and Ljung [15].

3.1

Overview

The closed-loop problem of controlling the multirotor system through feedback can
be visualized as in Figure 3.1. It is clear that there are no direct means of actuating
the horizontal positions x and y directly. Since positional control is desired, the
angular subsystem must be controlled by a positional controller in such a way that
a desired position can be reached.

Figure 3.1: The closed-loop multirotor system.

The state vector x of the multirotor system is given below in (3.1). These states
correspond to translation and rotation described in the inertial frame. The control
48

3.1

Overview

signal vector U to influence these states is shown in (3.2). Its contents correspond to
the thrust force U1 as well as actuated roll, pitch and yaw moments, U2 , U3 and U4 .

X= x

x y

y z z


U = U1

U2 U3

U4

T

(3.1)

T

(3.2)

Using the vectors in (3.1) and (3.2), the multirotor system in the inertial frame can
be described compactly on the form (3.3). The full state-space equations for the
system was previously described in Section 2.1, of which some simplifications will
be made for the controller design.
= f (X, U)
X

(3.3)

A generalized control problem can be schematically viewed as in Figure 3.2, with a

generalized control law according to (3.4). Apart from the reference signal r, there
are other input signals to the system representing different types of disturbances
input disturbance du and measurement disturbance n. External forces and moments
due to aerodynamic flight effects enter the system as du and will be of central interest
for this thesis. For satisfactory flight, it is important that these forces are sufficiently
attenuated by the controller.
u(t) = Fr r(t) Fy y(t)

(3.4)

du
r

Fr

n
y

Fy
Figure 3.2: A generalized control problem.
In this generalized control problem, the process output z is defined by a series of
transfer functions that map each input to z. The relation is shown in (3.5) with
49

Chapter 3. Control
relevant closed-loop transfer functions defined in (3.6). The transfer function Gcl
describes the response from r to the z, whereas Su , known as the input sensitivity
function, describes the response from du to z. For effective attenuation of input disturbance du , Su should be made as small as possible across all frequencies in the
Laplace domain.
z(t) = Gcl r(t) + Su du (t) T n(t)

Gcl =

Fy
PFr
=T
1 + PFy
Fr

Su =

(3.5)

P
= PS
1 + PFy
(3.6)

S=

1
1 + PFy

T=

PFy
1 + PFy

In (3.6), the S is the sensitivity function and T the complementary sensitivity function. The measurement disturbance n is mapped to z by T , whereas S maps n to the
measured output y. Obviously, Su = S = T = 0 would be desirable, but will not be
possible due to physical constraints and the limitation S + T = 1. Generally, S is
kept low at low frequencies and T kept low at high frequencies to achieve desirable
disturbance rejection and robustness to model uncertainty. More information on this
topic can be read in Glad and Ljung [15].
In the proceeding control analysis, only the influence of the input disturbance du
will be analyzed. This thesis work is performed at a conceptual level and therefore doesnt concern the real-world implementation aspect of sensor measurement
noise. Thus, modification of Gcl and Su to achieve satisfactory reference tracking
and rejection of input disturbance are the main aims of the controller design.

3.2

PID control

The abbreviation PID stems from that the controller type employs proportional (P),
integral (I) and derivative (D) control action. It is one of the most conventional and
proven control strategies around and often proves sufficient in many applications.
The controller aims to minimize the control error e(t) presented in (3.7). The signal r(t) constitutes the desired process output, also known as the reference signal,
whereas y(t) is commonly known as the measurement signal and is the measured
output from the system.
e(t) = r(t) y(t)
50

(3.7)

3.2

PID control

The standard PID controller can be represented as shown below in (3.8). The constants KP , KI , KD are chosen, or tuned, to obtain satisfactory closed-loop response
and performance.
Z

u(t) = KP e(t) + KI

e(t) dt + KD

d
e(t)
dt

(3.8)

Cascaded control
The multirotor system can be seen as two subsystems, one rotational (inner) and one
translational (outer), and the cascaded control structure proposed in this section is
therefore well-suited. PD controllers are designed to control the attitude ( and ),
whereas PID controllers are employed to control the position (x, y and z) as well as
the yaw angle (). The proposed structure of the cascaded controller can be viewed
in Figure 3.3.

Figure 3.3: The cascaded PID control structure.

Viewing the system as an inner and outer control loop simplifies linear control analysis. When adjusting a cascaded PID controller one starts with the inner loop, which
in this case is the rotational subsystem. It is tuned to desired performance. Moving
on to the outer loop, which here is the translational control loop, its speed (time
constant or bandwidth) will be limited by the inner loop. It is common practice to
make it a factor N times slower than the inner loop, which makes it possible to view
the inner loop as static from the outer loops perspective. It was decided that the
factor N will be at least four for the multirotor controller case.

51

Chapter 3. Control

Attitude control
When designing the controller for the rotational loop and analyzing the system response, simplifications are made to the dynamical system in order to decouple the
control signals and to remove nonlinear terms. Thus, linear analysis can be applied
instead of more complex analysis. The rotational velocity in the body and inertial
T [p q r]T , which has most vaframes are approximated as equal, i.e. [ ]
lidity at low rotational velocities. This is most common in near-hover flight. The
simplified control model is presented in (3.9), in which U2 , U3 and U4 constitute the
moment control signals, and aero are the corresponding moment disturbances due
to aerodynamic effects.

1
1

=
U2 +
aero

I
I

XX
XX

1
1
(3.9)
=
U3 +

IYY
IYY aero

1
1

=
U4 +

IZZ
IZZ aero
Roll and pitch control The attitude control laws designed to control the and
angles of the rotational subsystem are expressed in (3.10), where the error e(t) is the
angle error. The derivative part has been modified, compared to a standard PD controller, to act as pure damping of the system. The control output u(t) is seen as an
angular acceleration and is therefore multiplied by the moment of inertia around the
appropriate axis to get the associated control moment. Since the positional PID controllers generate angle references for the roll and pitch controllers, integral action
was not deemed necessary to control these angles.

u(t) = KP e(t) KD

d
y(t)
dt

(3.10)

In the Laplace domain, the transfer functions from reference and process disturbance to the output angle is shown for in (3.11). These are analogous to the
transfer functions for . It is clear that two poles are present, which must be have
negative real parts for the system to be stable.

(s) = Gcl (s)re f (s) + Su (s)aero

(s)

KP

Gcl (s) =

s2 + KD s + KP

52

(3.11)
Su (s) =

IXX s2 + K s + K
D
P

3.2

PID control

The parameters in the proposed control laws in (3.10) are tuned such that the two
system poles p1 and p2 in (3.11) are placed on the real axis and in the same point in
the complex plane, so p = p1 = p2 . To achieve this pole placement, the PD parame

ter expressions KP = p2 and KD = 2p are used. In Figure 3.4, step responses from
reference and disturbance signals can be viewed for controllers tuned for different
performance. The blue line shows the response for poles placed in p = 10, green
for p = 15 and red for p = 20. The black dashed line shows the reference.

Figure 3.4: Roll/pitch reference and process disturbance step responses for three
different PD control settings. After 0.5 seconds, a reference step of 5 degrees was

made. After two seconds, a moment input disturbance step of aero = 0.1 Nm was
2
made during one second on a multirotor with IXX = 0.04 kgm .
As can be seen from Figure 3.4, there will be a steady state error due to the disturbance that a PD controller can not perfectly compensate for. The controller tuned
with highest gain will have a faster reference response as well as significantly
smaller error due to the disturbance. Another way of viewing the closed-loop per

formance is by looking at the Bode diagram of Gcl and Su . The Bode diagrams
can be viewed in Figures 3.5 and 3.6, where it can be seen that the faster controller

yields a higher cut-off frequency in Gcl and lower steady state gain in Su .
Although a faster control loop shows beneficial results, increasing it too much is illadvised since it is limited by the speed of the motors. Because of this, it is concluded
to tune the controller such as the closed-loop response of the rotational loop is at
least a factor N times slower than the motor response, using similar reasoning as
discussed for cascaded control loops. Thus, from the rotational loops point of view,
the motor dynamics may be seen as a static gain achieving a desired rotational speed
instantaneously.
53

Chapter 3. Control

Figure 3.5: Bode plot of the Gcl transfer function for roll/pitch with different controller settings.

Figure 3.6: Bode plot of the Su transfer function for roll/pitch with different controller settings.

54

3.2

PID control

Yaw control In the proposed control structure, the position is controlled without
use of the yaw angle . It must, however, be stabilized to ensure flight performance.
To achieve this, integral action must be employed since for this angle there is no

outer controller generating references. Unwanted moment disturbances aero will,

hence, create a steady state error. Such moments may come from the fact that in a
real setting the motors and propellers may not be perfectly identical or they may
experience local variations in conditions. For a coaxial rotor configuration the propellers above and below will produce different counter-moments on the air frame,
which will create a situation where integral action is certainly needed. Thus, the
yaw angle will be controlled by a PID controller (3.12), with U4 = IZZ u(t).
Z

u(t) = KP e(t) + KI

e(t) dt KD

d
y(t)
dt

(3.12)

Using this controller, the resulting transfer functions for the yaw system are presented in (3.13). The introduction of integral action introduces an additional pole,

but since also a control constant KI is added, pole-placement techniques may again
be used to tune the system.

(s) = Gcl (s)re f (s) + Su (s)aero

(s)

Gcl (s) =

KP s + KI

3
s + KD s2 + KP s + KI

Su (s) =

s
1
IZZ s3 + KD s2 + KP s + KI
(3.13)

Arbitrarily placing the three poles yields the control parameters KP = p1 p2 + p1 p3 +

p2 p3 , KI = p1 p2 p3 and KD = p1 p2 p3 . The pole locations are chosen according to a dominant pole placement technique [46], where two poles are placed
closer to the imaginary axis in the complex plane such that p2 = p1 , p3 = 5p1 .
The performance will thus mainly be defined by the two dominant poles p1 and
p2 , somewhat approximating the system to second order with reduced overshoot.
Figure 3.7 shows the reference and disturbance step responses for the yaw system
with controllers designed for poles in p1 = 1 (blue line), p1 = 2 (green line) and
p1 = 3 (red line). The dashed black line is the reference value.

55

Chapter 3. Control

Figure 3.7: Yaw reference and disturbance step response for three different PID
control settings. After 0.5 seconds, a reference step of 5 degrees is made. After 10

seconds, a moment input disturbance step of aero = 0.1 Nm is made on a multirotor

2
with IZZ = 0.08 kgm .

The Bode plots for the Gcl and Su transfer functions can be viewed in Figures 3.8
and 3.9. Relating these to the step response plot (Figure 3.7), it is clear that the best
performing controller has higher cut-off frequency (and hence faster response) as
well as lower closed-loop gain for disturbances across the frequency range.

Figure 3.8: Bode plot of the Gcl transfer function for yaw with different controller
settings.
56

3.2

PID control

Figure 3.9: Bode plot of the Su transfer function for yaw with different controller
settings.

As for roll and pitch, the yaw controller can not be made arbitrarily fast. It must
adhere to the constraints set by the motor actuators. The yaw reference will be zero
in the proposed control setup since only the roll and pitch angles will be used for
position control. Hence, the yaw controller will be tuned to provide at least a two
times slower response than for the prioritized control of the roll and pitch angles.

Position control
When considering the outer translational loop, similar model simplifications as for
attitude are made to obtain the linear system (3.14). These simplifications yield a
model that doesnt accurately describe all flight conditions, but may still be used for
control analysis to give an indication of the closed-loop response characteristics.
Simplifications include small angle approximations as well as neglecting the coupling between angular and translational velocity. This is done in order to decouple
the control of the angular references.

1 x
1

x = U1 + Faero

m
m

1
1 y
y = U1 + Faero

m
m

z = g + 1 U + 1 F z
1
m
m aero

(3.14)

The PID controllers designed to control the translational subsystem have the famil57

Chapter 3. Control
iar expression in (3.15), where the error e(t) is the position error. To ensure that
no steady state position errors will be present, integral action of the controller is
essential. Again, the derivative part has been modified to act as pure damping of the
system.
Z

u(t) = KP e(t) + KI

e(t) dt KD

d
y(t)
dt

(3.15)

Horizontal control In horizontal control, reference tracking of desired x and y

positions are considered. From the outer control loops perspective the inner loop
is static, achieving a desired attitude instantaneously. The controller output of the
horizontal flight provides references for the pitch and roll controllers. The analysis
in this section will be for control of the x position, but may be applied analogously
to y.
Since there is no way to actuate the horizontal position directly, the control is dependent on the available thrust U1 and current attitude. To simplify analysis, the
thrust will be approximated for near-hover flight as U1 mg. Since the horizontal
position control sends angle references to the roll and pitch controllers, the input to
the inner subsystem is = u(t)/g.

x
X(s) = Gxcl (s)Xre f (s) + Sux (s)Faero
(s)

KPx s + KIx
Gxcl (s) = 3
s + KDx s2 + KPx s + KIx

Sux (s) =

1
s
m s3 + KDx s2 + KPx s + KIx

(3.16)

The PID parameters are chosen such that all poles are placed on the real axis to
prevent oscillation. For a given set of poles, the parameters are KP = p1 p2 + p1 p3 +
p2 p3 , KI = p1 p2 p3 and KD = p1 p2 p3 . Two poles p1 and p2 are placed in
the same location and are the ones that have greatest impact on the response time.
The last pole is placed a factor 5 further into the left half of the complex plane,
which reduces the overshoot and makes the closed-loop system behave somewhat
like a second order system. This dominant pole placement approach is similar to
the one adopted for yaw control. In Figure 3.10, the reference and disturbance step
responses are shown for three control settings. These controllers correspond to p1 =
0.5 (blue line), p1 = 1 (green line) and p1 = 2 (red line). The black dotted line
is the reference signal.

58

3.2

PID control

Figure 3.10: Position reference and disturbance step response for three different
PID control settings. After 0.5 seconds, a reference step of 1 m is made. After 10
x = 1 N is made on a multirotor with
seconds, a force input disturbance step of Faero
m = 1.2 kg.

Analyzing the Bode diagrams in Figure 3.11 and 3.12 for these controller settings
clearly shows differences in bandwidth and disturbance attenuation evident in the
step responses from Figure 3.10.

Figure 3.11: Bode plot of the Gxy

cl transfer function for x position with different
controller settings.
59

Chapter 3. Control

Figure 3.12: Bode plot of the Sux transfer function for x position with different controller settings.
Altitude control Starting with the linear expression for the altitude z in (3.14),
system analysis is performed in similar fashion as for horizontal control. The
PID controller is designed as before (3.15) and related to the system input U1 by
U1 = m(u(t) g). This way, the same transfer function structure as for the horizontal control is obtained for altitude (3.17). Note that the gravitational acceleration g is subtracted from the control signal since z points downward according to
aerospace standards [1].

z
(s)
Z(s) = Gzcl (s)Zre f (s) + Suz (s)Faero

Gzcl (s) =

KPz s + KIz
3
s + KDz s2 + KPz s + KIz

Suz (s) =

1
s
m s3 + KDz s2 + KPz s + KIz

(3.17)

The system is tuned according to the same principles as for the horizontal system,
making use of pole-placement technique, and will have an identical response in this
linear approximation. Thus, step responses and Bode diagrams may be omitted here
and Figures 3.10, 3.11 and 3.12 for the horizontal system may be applied to this
altitude analysis as well.
When the controller is implemented in the multirotor system, an anti-windup
scheme is employed to handle the problem of integrator windup due to saturation
in attainable thrust. The limits are [0, Tmax ], where the zero is obvious since rotors
60

3.3

Control allocation

with fixed blades can not be controlled to provide any downward force. Tmax is not
as straightforward since the available thrust will depend on the flight regime. The
used anti-windup scheme is known as back calculation [40], which subtracts the
difference between the desired and actual control signal from the integral part of the
control. Undesired effects when re-entering the non-saturated area are reduced or
avoided that would otherwise cause large overshoots.

3.3

Control allocation

To actuate the system, the four physical control signals U1 , U2 , U3 and U4 , representing desired thrust force and moments (roll, pitch and yaw) respectively, must be
recalculated to corresponding rotor speeds and voltages. These signals are stored in
the control signal vector U and are related to the squared rotor speeds 2 according to (3.18). The matrix A is here presented for the quadrotor case of two rotors
(i = 1, 3) spinning clockwise and two counter-clockwise (i = 2, 4) with individual
rotor speeds i . In A, the coefficients b = CT AR2 and d = CQ AR3 are the thrust
and drag coefficients (which vary with the current flight conditions, see Section
2.3), whereas l is the arm length. Here, it is important to remember that the forces
and moments are defined according to a coordinate system where the vertical axis
points downwards [1].

U1
U2
2

U=
U3 = A
U4

b
bl
A=
bl
d

b
bl
bl
d

b
bl
bl
d

b
bl

bl
d

(3.18)

The equation system (3.18) can be inverted in the quadrotor case to obtain the corresponding squared propeller speeds 2i needed to actuate the desired forces and
moments. The resulting solution is shown in (3.19) below.
2
1
1
4b
2
1
1
2
2

= 2 = A U = 4b
1
3
4b
2
1
4
4b

1
2bl
1
2bl

1
2bl
1
2bl

1
4d
1
4d

1 U
4d
1
4d

(3.19)

For other multirotor configurations, using additional rotors, an equivalent matrix A

will result in an over-determined equation system, yielding an infinite amount of
possible solutions. A way to achieve a solution with an even control distribution on
each rotor is to make use of the Moore-Penrose pseudo-inverse. The pseudo-inverse
solution is given by 2 = (AT A)1 AU [10]. Since the motor dynamics are very fast
61

Chapter 3. Control
compared to the vehicles motion, the calculated rotor speeds are used to calculate
the steady state motor voltages (2.12), which are then applied to the corresponding
DC motors.

3.4

Disturbance analysis

Wind disturbances enter the system as forces and moments and will thereby cause
translational and rotational acceleration. The drag force acts in the centroid, whereas
the rotor forces act in each rotor plane. However, the resulting levers to the center
of gravity will be quite small and thus rolling and pitching disturbance moments
will not be of great significance. Hence, PD controllers were assumed sufficient for
control of the roll and pitch angles, especially since these were connected to PID
controllers in the outer loop. Since the drag moment Q affects the yaw angle and is
susceptible to wind, it will benefit from integral action.
The forces have different dependence on the wind speed. Airframe drag, that has
its point of attack in the airframe centroid, is proportional to the squared airspeed v
according to Fd = 12 Cd Av2 . The squared velocity dependence makes the airframe
drag more prominent at high speeds and also yields faster attenuation of high frequency content. Rotor forces T and H, drag-moment Q and the tilting of rotor plane
due to blade flapping have more intricate wind dependence. The advance ratio
and inflow ratio are both affected by the airspeed and contribute to variations
of the proportionality variables CT , CH and CQ , for which expressions are given in
Appendix C.2. These expressions depend linearly and quadratically on and as
well as the flapping angles. It has, however, been observed that the horizontal rotor
forces have an approximately linear dependence on the wind speed.

Wind variation and turbulence

The Dryden wind turbulence model derived in Section 2.5 is used here to analyze
the frequency content of wind disturbances. Note that the Power Spectral Densities
(PSD) (2.34) and corresponding transfer functions (2.36) bandwidth and magnitude changes with the relative airspeed V and the flight altitude h. Both the bandwidth and magnitude quantities decrease with the altitude and increase with the
velocity.
Plotting the PSD for wind in each linear direction yields the diagrams in Figure
3.13, where the red line is for the longitudinal direction, green the lateral direction
and red the vertical direction. The diagrams were created for an altitude of h = 10
m, a mean relative airspeed of V = 20 m/s and wind speed W20 = 15 m/s at 6 m
above the ground. The diagrams have been normalized relative to their respective
maximum PSD.

62

3.4

Disturbance analysis

Figure 3.13: Longitudinal (red), lateral (green) and vertical (blue) Dryden wind
turbulence normalized power spectral density diagrams.
As can be seen, the lateral and longitudinal variations are attenuated at a much
lower frequency (close to 1 rad/s for a factor 10 in the figure) than for the vertical
turbulence.
With the previous discussion in mind regarding force and moment disturbance due
to wind, it will be important that the controllers are designed in such a way that they
respond well to these. To cope with the variations shown in the PSD diagrams in
Figure 3.13 effectively, the positional controllers must be tuned such that the magnitude of Su is small enough (max(Su ) < 1) and bandwidth of Gcl is high enough for
these frequency ranges. As previously discussed, arbitrarily high control bandwidth
is not possible to obtain due to physical limits and may not be desirable due to influence and amplification of measurement noise present in a real-life implementation.
Thus, a trade-off is typically performed to achieve satisfactory performance.

Wind magnitude and gusts

An increase of wind magnitude will not directly require high-performance control
in terms of fast response (unless the onset of the wind is rapid). However, high
enough wind speeds will challenge the physical limits of the control algorithms and
actuators. For starters, enough thrust force must be available to counteract a wind in
steady state in order for the UAV to hold its position. Saturation of motors running
at their maximum voltage will be a problem. Also, if the UAV gets too far from its
reference position (as a result of temporary onset of strong wind or saturation), the
position error may become too large and yield undesirably aggressive control action
which may topple the multirotor. To amend for this, built-in saturation on control
signals or reference set-pointing is appropriate.
63

Chapter 3. Control

Lift and drag variations

Depending on the flight mode, the lift and drag coefficients b and d will vary. They
are directly proportional to CT and CQ which depend on the advance ratio and
inflow ratio , increasing with and decreasing with . This is explained in further
detail in Appendix C.2. Not accounting for the variations of lift and drag coefficients
when controlling the system has undesired effects to wind-robustness, making the
multirotor unsteady in the vertical direction. Knowledge of these coefficients will
help with both the control allocation as well as predicting the current motor load,
which is needed to calculate appropriate input voltages. In this thesis, the mean
value of these are assumed to be known at each instance, which isnt necessarily the
case. For robust forward flight and wind-rejection, using a wind-estimating algorithm with an airflow sensor may significantly improve the performance, although
responsive control may somewhat address the issue.

64

4
Multifactor analysis

The aim of this chapter is to investigate how the performance of a multirotor is affected by different parameters. To achieve this, a large number of simulations have
been conducted. A traditional and often intuitive approach is to test each parameter one at a time until the optimal performance is obtained for each parameter.
This method is called the COST approach [11]. Unfortunately there are some major drawbacks with this method. The solution depends on the starting values and it
doesnt guarantee that the found solution is the optimal design. The COST approach
is also inefficient, especially when a lot of parameters are considered. In Figure 4.1,
the COST approach is performed with two factors, X1 and X2, that are increased
separately until the performance is no longer improved and an optimum is believed
to be found. From the plot, however, it is evident that this is not the case. It was
therefore concluded that an alternative method is needed.

Figure 4.1: COST approach optimization.

65

Chapter 4. Multifactor analysis

4.1

Design of experiments

In Design Of Experiments (DOE), a carefully chosen set of experiments are conducted where all factors are varied simultaneously. The DOE approach is a reliable
and efficient tool that is used in a lot of different industrial sectors such as the chemical industry, the car manufacturing industry and the food industry [11]. An example
of this method can be seen in Figure 4.2, where it is applied to the same case as for
the COST approach mentioned above. The design space is in this case made of a
rectangular pattern, which means that a trend in the results can be identified. This
will eventually lead to an optimal solution.

Figure 4.2: Design of experiments optimization.

Design of experiments can be used for a number of different phases during the
development of a product. Familiarization, screening, optimization and robustness
testing are examples of such phases. This report will only cover the screening and
the optimization, but if the reader is interested, each of the experimental objectives
are covered in Design of Experiments: Principles and Applications [11].
For a multirotor, there are a vast number of different parameters which can potentially have an important effect on how the vehicle performs. Everything from
dimensions and masses of all the individual parts to the number of propellers and
their arrangements. It was concluded that it would not be possible to design a scalable multirotor model that can assume all possible design variations during the time
frame of this thesis work. Therefore a number of limitations have been introduced
in this chapter:
The DOE evaluation will only deal with the quadrotor design which is the
most common and basic multirotor concept.
66

4.1

Design of experiments

Only one type of propeller will be used.

Only the cascaded PID controller will be used during the DOE evaluation.
However, there will be a factor which represents either of two different PID
parameters setups. This will give an indication of how much the responses
are affected by software and physical characteristics respectively.
It should be noted that the multifactor analysis in this report is mainly aimed at
investigating the influence of different parameters, and not so much at finding the
optimal parameter combination which is usually the goal. The reason for this is that
the exact values of each parameter are not likely to be valid for different airframe
or propeller types, while the identified influence trends should be. The different
steps of the analysis will nevertheless be called screening and optimization since
the methods in each step will be the same.

DOE software
It is possible to conduct a DOE without the aid of specialized software, but a lot of
time can be saved by using such programs . There are several programs dedicated to
design of experiments on the market today, but the alternatives for this report were
narrowed down to the following two candidates:
1. MODDE by Umetrics [30].
2. modeFRONTIER by ESTECO [12].
Most of the theory in this section is based on the book developed by Umetrics [11].
However a number of different advantages spoke in the favor of modeFRONTIER
which will be used during this project. The most notable advantages are listed below:
The software is available at LTH.
The authors of this report have some previous experience with the software.
The software is compatible with MATLAB which means that each design in
the DOE can be run sequentially without any manual interaction.
It should be noted that MATLAB has a set of functions based on design of experiments that could have been used for this project, but it was preferred to use
modeFRONTIER which offers a more complete package.

67

Chapter 4. Multifactor analysis

Experiment setup
A specified flight scenario is needed to be able to draw any conclusions about the
performance of the quadrotor. It was decided to set the quadrotor in position hold
while an external wind acts as a disturbance. A typical setup for the wind experiment includes choosing a mean wind speed and direction angle. Additional settings
include turbulence and gust settings, which enhances the realism of the simulated
experiment. The wind test that was used during the DOE simulations is shown in
Figure 4.3. An initial mean wind of 8 m/s is used while a wind gust as well as
variations predicted by the Dryden wind turbulence model is included. Section 2.5
gives a theoretical description of how these phenomena are modeled.

Figure 4.3: Typical wind test design which has been used for evaluation of the multirotor performance. The plots show wind speeds in the inertial longitudinal, lateral
and vertical directions.

Screening
During a screening process, the main goals are to detect the factors that have the
most influence on the chosen responses, as well as finding the factor intervals in
which the optimal design is most likely found. According to the Pareto principle, 20
% of the factors are responsible for roughly 80 % of the influence on the responses
[11]. Therefore the most time efficient way of approaching an optimization problem
is to exclude the factors which doesnt have any major impact on the final result as
early in the process as possible. By doing so, more resources can be spent on tuning
the factors that greatly affect the responses. It is also important to make sure that the
ranges for the chosen factors includes the optimal design, or else that design wont
be found during the optimization phase. The information which is gained from the
screening is relatively easy to obtain, so the number of experiments needed per
factor is quite low.
68

4.1

Design of experiments

Factors Nine factors have been identified and used during the screening phase.
These were chosen in the most general way possible so that the knowledge gained
from the results could easily be applicable to different multirotor design concepts.
Axy [0.02, 0.04] [m2 ]: The projected surface area that is perpendicular to the
x- and y-axis. Because of the axisymmetric design, these areas are equal for
the quadrotor.
Az [0.04, 0.1] [m2 ]: The projected surface area that is perpendicular to the
z-axis.
Cd [0.2, 0.4]:The drag coefficient of the airframe. This is normally a function
of the relative airspeed direction. However, as mentioned above in Section
2.4, it was determined to use the same drag coefficient along the entire airframe .
I [0.048, 0.08] [kgm2 ]: The norm of the inertia tensor is defined as in (4.1),
where Ixx, Iyy and Izz are the moments of inertia around the x-, y- and z-axes.:
p
I = Ixx2 + Iyy2 + Izz2
(4.1)
la [0.3, 0.4] [m]: The arm length which is defined as the distance between the
CoG of the quadrotor and the rotational axis for the propellers.
mtot [1, 1.5] [kg]: The total mass of the quadrotor.
{0.5825, 0.8737} [s]: The time constant of the positional control loop.
This factor depends on the chosen control parameters, and it is therefore not
straightforward to change this factor in a continuous way. With this in mind, it
was decided to use the time constant as a qualitative factor with two different
settings.
zCG [0.1, 0.1] [m]: The vertical distance between the CoG of the quadrotor
and the centroid of the airframe.
zPG [0.1, 0.1] [m]: The vertical distance between the CoG of the quadrotor
and the CoG of the propellers.
The factor ranges were mainly determined by rough estimations for what could
be realistic limits with respect to the multirotor type and the chosen propellers.
The most obvious example being that the total mass of the quadrotor is limited by
the desired thrust margin. The quadrotor model consists of multiple body objects
from the SimMechanics environment. There is no easy way of manipulating the
characteristics of an entire multi-body system in SimMechanics, which means that
69

Chapter 4. Multifactor analysis

only masses and dimensions for the individual parts can be used as inputs to the
simulation. An algorithm was therefore needed to translate the previously listed
factors to compatible inputs. To simplify this algorithm it was decided to keep some
of the geometrical dimensions constant in the mechanical sub-model. By doing so
it was possible to determine the three different masses for the center, motors and
arms based on the total mass, the inertia and the arm length. The projected surface
areas were therefore only changed in the aerodynamic submodel, where the drag
forces acting on the airframe were calculated. The correlation between I, mtot and
la means that the intervals had to be carefully chosen so that the mass distribution
would be realistic.
Responses This report is mainly focused on wind resistance, but it is also interesting to study how the flight time is affected by changing the different factors. The
following responses were therefore chosen:

Maximum position error, emax = max(e) [m]: The maximum distance that the
quadrotor deviates from the desired position during the applied wind disturbance. The distance error e is defined as:
e=

q
(x xre f )2 + (y yre f )2 + (z zre f )2

(4.2)

Mean power, Pmean [W]: The total mean electrical power during the applied
wind disturbance is calculated as:
4

1
Pmean =
t

t0
f
i=1

Zt f

Vi Ii dt

(4.3)

t0

It is worth noting that a decreased power isnt necessarily equivalent to a

longer flight time. A quadrotor with an increased power, but higher thrust
margin, could potentially equip more batteries which could result in a longer
flight time. However, since accurate modeling of the battery was considered
to be outside the scope of this report, it was concluded that the mean power
is a strong indication of how the possible flight times are affected.
Selection of regression model and DOE design When the factors and responses
have been chosen, it is time to choose a DOE design. There are a number of different designs of varying complexity and they are all more or less suited to specific
problems. The choice of DOE design depends on which type of regressions model
that is desired, but also on the number of factors that are being considered. The three
most commonly used regression models are [11]:
70

4.1

Design of experiments

Linear: y = 0 + 1 x1 + 2 x2 + +
Interaction: y = 0 + 1 x1 + 2 x2 + 12 x1 x2 + +
Quadratic: y = 0 + 1 x1 + 2 x2 + 11 x12 + 22 x22 + 12 x1 x2 + +

The listed models are sorted after simplicity in descending order. If a more complex model is chosen, the potential to capture the behavior of the real process is
improved, but such a model will also increase the number of experiments that needs
to be conducted. A linear or interaction model is normally proposed for a screening
process [11], depending on the number of factors. A linear model can be accomplished by a fractional factorial design while an interaction model needs a full factorial. In addition to these designs, there are also a number of other DOE designs
that are suited for the screening process, such as Plackett-Burman, Rechtschaffner,
L-design and D-optimal design. These are, however, considered to be outside the
scope of this project.
The full factorial design is the most basic DOE design and therefore the most easily
explained. It is common to define a standard experiment which serves as a reference.
If the rest of the experiments are then symmetrically distributed around this center
point, a full factorial design is obtained. A two-level full factorial design for three
factors can be seen together with a reference point, also known as center point, in
Figure 4.4. This means that the three factors are tested for two values each. In this
report there have been nine different factors identified that may or may not have a
significant impact on the responses. The full factorial design is mostly just relevant
for up to five factors, after which the number of experiments are increased rapidly,
see Table 4.1 [11]. It was therefore decided to use the fractional factorial which
translates to 32 simulation runs, plus any number of center points that are being
used. A two-level fractional factorial design for three factors can be seen together
with a center point in Figure 4.5.

Figure 4.4: Full factorial design.

Figure 4.5: Fractional factorial design.

71

Chapter 4. Multifactor analysis

Number of
investigated factors
2
3
4
5
6
7
8
9
10

Number of runs,
Full factorial
4
8
16
32
64
128
256
512
1024

Number of runs,
Fractional factorial

4
8
16
16
16
16
32
32

Table 4.1: Number of experiments needed for full and fractional factorials for a
given number of investigated factors.

There are mainly two benefits of using center points in DOE [32]:
1. To provide a measure of process stability and inherent variability.
2. To check for curvature.
The experiments conducted in this project utilizes a simulated model with no variability, and hence there is no need for multiple simulations with the same center
point. However, one center point can be useful to detect whether the relationship
between a factor and a response is linear or not. Note that the time constant factor,
, is a qualitative factor which means that it is not possible to place a center point in
the middle of the hypercube. Instead two different pseudo center points have to be
used, one for each value of the qualitative factor. An example of how this can look
for three factors is depicted in Figure 4.6.

Figure 4.6: Center points for a full factorial with one qualitative factor.
72

4.1

Design of experiments

Results When all of the experiments from the DOE has been run it is time to interpret and present the results from the simulations. There are numerous tools available
in modeFRONTIER to help with this task. The RSM-tool automatically creates a
regression analysis based on the raw data from the simulations. The coefficient of
the linear regression models can be seen in Figure 4.7 and 4.8. Large coefficient
values means that the response is largely influenced by the corresponding factor.
Keep in mind that the coefficients have been normalized with respect to the factor
intervals.

Figure 4.7: Linear regression coefficients for emax .

Figure 4.8: Linear regression coefficients for Pmean .

The difference between the simulated results and the predicted values from the linear regression model is presented in Figure 4.9 and 4.10. The blue dots are the values
predicted by the regression model and the orange dots are the simulated values. The
linear regression models are relatively poor, but this was to be expected since the
quadrotor behavior is fairly complex. The average relative error is 35.3% for emax
and 5.53% for Pmean . Despite the questionable regression model, the collected data
was still considered good enough for the screening phase, where the aim is only to
roughly estimate the influence of the different factors.
An alternative method of presenting the information is to use a so-called student
73

Chapter 4. Multifactor analysis

Figure 4.9: Linear regression validation

for emax .

Figure 4.10: Linear regression validation for Pmean .

chart [12]. The influence of each factor is shown in Figure 4.11 and 4.12 for the
maximum position error and the mean power respectively. The height of the bars
is called effect size and shows the strength of the relationship between the output
and the input values. It is calculated as the difference between M and M + which
are the mean values of the responses for the lower and upper ranges of the factor
intervals. Low values indicate that there is not any detected relationship between the
input and output variables. The transparency of the bars is a measure of a variables
so-called significance. High transparency means that the effect size of that variable
isnt very reliable. The red bars have a direct effect on the response. This means
that an increased factor value will increase the response value. The blue bars have an
inverse effect which consequently means that an increased factor value will decrease
the response. One of the advantages of the student chart is that it is possible to see
the certainty of the influence of each factor individually. The information is also
presented in a more easily interpreted way, especially when higher order regression
models are needed.

Figure 4.11: Student chart for emax .

74

4.1

Design of experiments

Figure 4.12: Student chart for Pmean .

Conclusions The aim of the screening phase is normally to eliminate any potential
factors that dont have any influence over the responses. It is possible to distinguish
a couple of factors which dont appear to have any influence over the responses
by looking at the regression coefficients and the student charts. la , I and zPG dont
affect neither the position error or the power in any significant way. The mean power
isnt affected by or zCG either while the position error is somewhat influenced by
zCG ,even if that information isnt completely reliable which is suggested by the
transparent bar in Figure 4.11.
The arm length la determines the maximum attitude moments for a given actuator.
The desired angular velocities of the quadrotor are determined by the control algorithm which takes the arm length into account. Hence, for a given set of control
parameters, the length of the arm shouldnt have any impact on the position error.
For a longer set of arms it would theoretically be possible to tune the attitude controller more aggressively which can be beneficial for the wind resistance potential.
In reality, there is a strong correlation between the length of the arms and the projected surface areas which are preferably minimized according to 4.11. It is worth
noting that because of the potential aerodynamic interference between different rotors, there is a lower limit to the length of the arms which is determined by the
size of the propellers and the number of arms. If manual flight is considered, longer
arms will generally make multirotors easier to control in a stable manner. This is
caused by the fact that an increased arm length will make a certain angle error more
apparent to the human eye. la will not be further studied in the optimization, but it
would be interesting to study how the the total drag and performance of the quadrotor is affected by having the rotors close to each other and closer to the rest of the
equipment which is often situated at the center. This would however require more
advanced modeling techniques than those covered in this report.
The inertia of the quadrotor determines the angular acceleration caused by a given
75

Chapter 4. Multifactor analysis

moment. This inertia is also taken into account by the controller which means that
the inertia doesnt affect the maneuverability for a given set of control parameters.
However, the tuning of the attitude control is limited by the inertia which means
that a smaller and lighter multirotor generally allows more agile flight. Intuitively
an increased inertia should be beneficial for rejecting moment disturbances, and this
is also confirmed by the input sensitivity transfer function in (3.11). It is therefore
surprising that the results doesnt show any correlation between the inertia and the
maximum position error. One explanation for this could be that the translational
velocities are constant along the quadrotor body which could potentially mean that
the disturbance moments are relatively small compared to the translational forces.
However, this explanation is contradicted by the influence of the vertical distance
zCG which causes a disturbance moment. After some additional testing it was concluded that the inertia does have an affect on emax , but its influence is either direct
or inverse depending on the signs of zCG and zPG . It was therefore decided to further
study these factors during the optimization.
The total drag of the airframe, which is determined by Axy , Az and CD , greatly affects
the position error as well as the power for the chosen factor intervals. This doesnt
come as a surprise, especially since the drag force on the airframe is proportional
to the squared wind magnitude. It is desirable to minimize these three factors for
both of the responses and they will therefore not be further investigated during the
optimization since the optimal design can safely be assumed to be found when they
are minimized.
The same logic can be applied to the time constant . This factor has the biggest
influence over the position error, but the performance will only be increased for
a decreased time constant as long as the closed-loop system is stable. It was nevertheless interesting to see how much the software can improve the performance
compared to the physical characteristics of the vehicle.
The total mass mtot is also very significant for the performance of the multirotor.
It is already clear that it will have a conflicting influence on the maximum position
error and the mean power. It will therefore be further studied in the optimization.

Optimization
During an optimization the aim is normally to find the optimal combination of factors which will give the best performance. Therefore one or more objectives have to
be defined so that the solver knows if a response should be maximized or minimized.
Often these objectives are in conflict with each other so that a compromise has to be
made. The optimization requires a lot of runs per factor in contrast to the screening
process where just a few simulations have to be conducted for a lot of factors.The
aim of the optimization in this report isnt primarily about finding the exact value
for each factor, but rather to further investigate the effect of the individual factors
76

4.1

Design of experiments

and the relationship between the two response objectives.

Factors The influence of the different factors have been discussed during the
screening phase and it was concluded that mtot , zCG , zPG and I would be further
investigated during this optimization. The optimization was run with a time constant setting of 0.8737 s, which is the slower setting. This setting is probably closer
to what the time constant would be in a real implementation. Axy , Az and CD were
locked at the minimum values which will give the optimal design combination
within the specified factor intervals. The arm length la was locked to the center
value of the interval, but this choice shouldnt have any significant effect on the
results.
Responses The same responses as in the screening will be used during the optimization. It is desirable to minimize both of these responses, which means that two
objectives have to be added to the responses during the optimization. The workspace
in modeFRONTIER can be seen in Figure 4.13, where the optimization phase is
modeled.

Figure 4.13: Workspace in modeFRONTIER.

Selection of regression model and DOE design It is very important to have a
good model of the process during the optimization process so that the solver can
converge to the real optimum. Therefore a quadratic regression model is commonly
used, but sometimes even more complex models are necessary. In this report it was
decided to use a central composite design since it is an intuitive extension of the
factorial designs which also supports quadratic regression models. There are three
basic central composite designs; circumscribed (CCC), face-centered (CCF) and
inscribed (CCI) which are pictured in Figure 4.14 for two factors.
Only the CCI and the CCF designs are available in modeFRONTIER. The CCF
77

Chapter 4. Multifactor analysis

Figure 4.14: The three basic central composite designs visualized for two factors.
design covers all of the test range while the CCI design tests each factor at five
different levels which potentially gives the CCI design slightly better possibilities
to make an accurate quadratic model. In this case the authors decided to use the CCF
design since some of the factors will ultimately tend towards either the maximum
or the minimum setting during the optimization.
Optimizer It is necessary to choose an optimizer which will work towards the
specified objectives. There are a number of different optimizers available in modeFRONTIER and they are suited for different tasks depending on things such as
the number of objectives, the type of the objectives and the number of factors and
responses. In this project it was decided to use MOGA-II which is one of the basic
optimizer choices in modeFRONTIER. MOGA stands for multi objective genetic
algorithm and it is designed for fast Pareto convergence according to the software
guide. Pareto designs can be described by the following definitions [21]:

A design point dominates another if it is superior with respect to one of the

objectives while not being inferior with respect to any of the other objectives.
Two design points are indifferent with respect to each other if neither point is
dominating the other.
A design point is a Pareto optimal design if it is not dominated by any other
design points.

The so called Pareto frontier is a curve that is defined by all the Pareto optimal
design points. The optimal design lies somewhere along this curve, depending on
the relative importance of the different objectives.
Results During the optimization phase, a total number of roughly 550 simulations
were run including the center point and the ones that were generated by the central composite face-centered design. The RSM-tool was used to create a quadratic
regression model whose coefficients can be seen in Figure 4.15 and 4.16.
78

4.1

Design of experiments

Figure 4.15: Quadratic regression coefficients for emax .

Figure 4.16: Quadratic regression coefficients for Pmean .

Figure 4.17: Quadratic regression validation for emax .

Figure 4.18: Quadratic regression validation for Pmean .

79

Chapter 4. Multifactor analysis

The regression model values are compared to the simulation results in Figure 4.17
and 4.18. This regression model is much better than the linear version that was used
in the screening section. The average relative error was reduced to 0.409% for emax
and 0.0818% for Pmean .
The influence of each factor can also be seen in the student charts pictured in Figure
4.19 and 4.20.

Figure 4.19: Student chart for emax .

Figure 4.20: Student chart for Pmean .

The objectives for the optimizer were to minimize the mean power and the maximum position error. These objectives oppose each other for some of the parameters,
so that there is not a single best design solution. The Pareto frontier can be seen
in the scatter plot in Figure 4.21 which contains all of the simulated designs. The
highlighted dots are the Pareto optimal designs.
The best factor combinations for emax and Pmean can be seen in Figure 4.22 and 4.23
respectively.
80

4.1

Design of experiments

Figure 4.21: Pareto frontier from the optimization run.

Figure 4.22: The best designs with respect to wind resistance.

Figure 4.23: The best designs with respect to efficiency.

Conclusions The regression model that was developed during the optimization
provides a much better fit for the simulated data. Conclusions can therefore safely
be made from the regression coefficients as well as from the student charts, where
the transparency of the bars are now generally much lower.
The total mass mtot should be maximized to achieve the highest possible wind resistance. If agile flight is not a necessary criteria, the thrust margin should therefore be
kept relatively low. The rotational speeds of the rotors will generally decrease when
81

Chapter 4. Multifactor analysis

a wind disturbance is present. Saturation of the actuators should therefore not be a
problem up until a certain limit for the rotor disk angle of attack where the thrust
will start to decrease again. However, it is of utmost importance to minimize the
mass if the power is going to be minimized. This is the main reason why there has
to be a trade-off between power and the maximum position error, which is also confirmed by the Pareto frontier in Figure 4.21. It is clear that an eventual final design
choice would depend on the relative importance of the objectives.
A common belief among hobbyists is that the rotor plane should be placed above the
overall center of gravity to achieve the best stability. This belief originates from the
intuitive comparison to an inverted pendulum. However, this parable is faulty since
the thrust of a multirotor will always point through the CoG if properly balanced.
Pounds et al. [34] show that to achieve the best stability for the X-4 Flyer, the rotor
plane should optimally be aligned with the CoG. It is also noted that a small error of
this placement will lead to significant changes for the system dynamics, and hence
the rotor plane is placed with a small offset relative to the CoG. The results in 4.19
shows that zPG should be minimized i.e. the rotor plane should be placed under
the center of gravity. One explanation for this could be that the moment caused by
the horizontal force will support the actuated pitch/roll moment that is necessary to
tilt the thrust force so that it can counter the translation drag force. The desirable
minimization of zCG can be explained similarly.
The inertia of the quadrotor does seem to have a relatively large influence on the
position error, as opposed to the screening results. It is apparently desirable to minimize the moment of inertia with respect to the wind resistance, which is slightly
surprising by keeping (3.11) in mind. However, note that this input sensitivity function is for the attitude while the response during this DOE has been focused on the
position error. Also remember that the optimizer has minimized the vertical distances zPG and zCG , which means that the majority of the design combinations have
a negative value for these factors. Therefore, as explained above, the disturbance
moments caused by the wind will have a positive effect on the positional stability.
The angular acceleration towards the wind is increased by minimizing the moment
of inertia. If zPG and zCG would have been positive, the disturbance moments caused
by the wind would instead have a negative effect on the positional stability. This is
the reason why the inertia didnt seem to have any effect on emax during the screening.
To summarize the multifactor analysis of the scalable quadrotor, the wind resistance
and the efficiency are conflicting objectives for the quadrotor. This is mostly caused
by the mass of the vehicle which should be maximized with respect to emax and
minimized with respect to Pmean . Axy , Az and Cd , which determines the total drag
of the quadrotor airframe, should intuitively be minimized which is also supported
by the results from the screening. The rest of the factors might not be quite as
82

4.1

Design of experiments

influential, but if the airframe is designed in an aerodynamically favorable way, the

tuning of these factors will have an increased influence on the performance. The
vertical distances zPG and zCG should ideally be negative, but this can be hard to
implement in a conventional multirotor design, especially if camera equipment is
mounted below the body.

83

5
Design development
The work presented in this chapter aims to propose a final multirotor concept suitable for windy conditions. Considering the limited time frame of this project, it was
not deemed possible to construct fully functional models of every potential multirotor concept. Instead, it was decided to base the final concept choice on research,
simplified simulations and engineering intuition. The chapter is mainly divided into
two sections, listed below.
Rotor and airframe assessment.
Concept evaluation.
Concepts for rotor and airframe configurations will be assessed separately in regards
to wind resistance and efficiency properties. Knowledge and conclusions drawn
from this will then be used, together with the results from the multifactor analysis in Chapter 4, during the proceeding concept evaluation. During this evaluation,
a number of potentially favorable multirotor concepts are generated and compared
in a concept scoring process. The scoring is based on a number of specified criteria, which will determine the final design choice. Concept scoring is a frequently
used method for comparing different concepts during a design development process. More information on design development can be found in Product Design and
Development by Ulrich and Eppinger [44].

5.1

Rotor and airframe assessment

During this section, different rotor and airframe configurations are studied separately. By doing so, it was possible to exclude less advantageous solutions before combining rotor and airframe configurations into complete multirotor concepts. Thus, a decrease in the number of possible multirotor concept candidates
is achieved.

84

5.1

Rotor and airframe assessment

Rotor configurations
There are several different rotor types more or less suited to a specific task and/or
design restrictions. The aim of this section is to determine which configurations that
are most suited for a wind resistant multirotor, while also considering the influence
they might have on power efficiency. The conclusions from this chapter are mostly
based on available literature and logical reasoning, using experience gained when
constructing the model as well as the results from the multifactor analysis in Section
4. In addition, a number of simplified rotor simulations have been conducted to
verify some of the reasoning.
During these simulations, different rotor configurations were compared at certain
rotational speeds, influenced by wind with varying magnitude V and inflow angle
. The interval for V was decided based on the assumed limitations of the model.
Also, only positive were evaluated since multirotors generally pitch down towards
the incoming wind in order to compensate for the force applied by the wind. The
upper limit of was decided based on that the pitch angle of the modeled quadrotor
rarely exceeded this angle during high wind disturbances. It is important to keep in
mind that these simulations only make use of the rotor models and doesnt take the
dynamic movement of the multirotor into account. The responses that were used to
evaluate the performance of the rotor configurations are listed below.

Dynamic power loading PLd =

Equilibrium angle eq = arctan

Fz
P

Fxy
Fz

where Fxy =

q
Fx 2 + Fy 2

Power loading is an efficiency parameter frequently used for different types of rotorcraft. It is defined as the ratio between the thrust and power required to hover
[18].

PL =

W
P

(5.1)

However, the authors of this report were also interested to see how the thrust to
power ratio changes for different flight conditions. This ratio was decided to be
called dynamic power loading, but will henceforth simply be referred to as power
loading PL.
In Figure 5.1, it can be seen that the expression in (5.2) holds when a force equilibrium along the horizontal plane is achieved.
85

Chapter 5. Design development

tan eq =

Fxy
Fz

(5.2)

Note that the equilibrium angle is only an indication of how much the multirotor
would have to pitch down to be able to compensate for the wind at a specific moment. In reality, or when running the full multirotor simulator, the angle of attack
for the rotor disk will continuously change as the multirotor starts to pitch down
towards the wind. There will also be a delay until the multirotor has actuated the
pitch maneuver, during which the force will cause an acceleration that depends on
the total mass of the vehicle. However, this delay is reduced by minimizing the
equilibrium angle.

Figure 5.1: Equilibrium angle explanation.

Propeller
There are a lot of different factors to consider when evaluating the performance
of a rotor, even when only considering regular propellers. The propeller can have
varying dimensions and operate at different speeds depending on the payload. Additionally, more complex solutions such as variable pitch propellers can be used. In
this section, some of these factors will be investigated further to give a better understanding of how they may affect the performance of the multirotor. The two-bladed
propeller is without a doubt the most commonly used rotor setup. Therefore it will
be used as a reference throughout this section.
Simulation The APC 10x4.7 Slowflyer propeller, used during the quadrotor simulations, was used as the reference for the rotor simulations. The equilibrium angle
eq and the power loading PL are plotted as functions of V and in Figure 5.2.
These results can partially be predicted by the expressions for the nondimensional
coefficients in C.23, and they are sources for control disturbances which are normally not taken into account.
86

5.1

Rotor and airframe assessment

Figure 5.2: Equilibrium angle and power loading for the reference propeller, APC
10x4.7 sf, running at 425 rad/s.
Propeller size and operating speed It would be very time consuming to investigate all of the possible parameters that determine the performance of the propeller.
For this reason, it was decided to focus on the influence of the propeller size and the
rotational speed at which it operates. To this end, a factor Ks was introduced to scale
the radius and the average chord length of the propeller, with mass and rotational
inertia scaled accordingly.
Simulation Evaluating both the scaling factor Ks and the rotational speed simultaneously spoke in the favor of using another DOE design. Factor intervals for these
two quantities are shown in Figure 5.3, where a two-level full factorial is depicted.
Note that the center point is the mentioned reference whose results were presented
in Figure 5.2.

Figure 5.3: Design of experiments for the propeller simulations.

87

Chapter 5. Design development

The factor Ks was kept quite close to one since even small changes to the radius and
the average chord length have a relatively large impact on performance. Also, the
authors didnt want to deviate too much from the original propeller dimensions, for
which the model has been validated with real experimental results. The rotational
speeds were determined with the aim of covering the most frequented operating
ranges. Results from the simulations are shown in Figure 5.4 and 5.5 respectively.
For an easier interpretation of the results, eq and PL have been normalized with the
corresponding results from the center point.

Figure 5.4: Design of experiments result: eq

It is evident from the simulation results that, within the investigated intervals, it is
desirable to use small and fast spinning propellers for minimum wind disturbance.
In contrast, it is shown to be more efficient to use small and slow spinning propellers. These results are assumed to apply for all rotor configurations covered in
this report.
88

5.1

Rotor and airframe assessment

Figure 5.5: Design of experiments result: PL

Conclusions The results for the power loading can easily be predicted by considering the expressions derived by MT (C.5), where T (R)2 and P (R)3 . This
confirms that small and slow spinning propellers a favorable. However, such propellers will not generate as much thrust. When designing helicopters, the choice is
usually limited to either a small and faster spinning propeller, or a larger and slower
spinning propeller. Choosing between these two options, the general opinion is that
larger propellers are more energy efficient which is also confirmed by the results
in Figure 5.5. Multirotors are not limited to a single rotor and could therefore potentially use multiple small and slow spinning propellers. However, propellers need
to have some spacing between each other to avoid aerodynamic interference [16],
which makes it difficult to implement a large number of propellers in practice. Furthermore, if a maximum footprint diameter is given as in [9], the conclusion will be
that a large main rotor is favorable from an efficiency perspective.
The results for the equilibrium angle can not be predicted in a similar fashion. According to the results in Figure 5.4, a small and fast spinning propeller is favorable
89

Chapter 5. Design development

for wind resistance. It is also worth noting that a smaller propeller can spin faster
since the maximum allowed rotational speed is limited by the blade tip speed.
Larger propellers will have higher blade tip speeds, for a given rotational speed,
which might cause compressibility effects and noise. In addition, smaller propellers
have a smaller moment of inertia which will reduce the time constants of the motors.
Multi-bladed propeller Sometimes it is preferable to use more than two blades on
a propeller. Usually, this is caused by limitations concerning available space, structural properties or actuator performance. There is, however, a limit to the number of
blades that can be used since they will interfere with each other if placed too close
together. This limit depends on a number of factors such as rotational speed and
propeller dimensions, but the theory behind this will not be covered in this report.
Simulation It was decided to compare a three-bladed propeller with the reference propeller since this is the most common alternative. It has been assumed that
any benefits or drawback obtained by using three blades will also apply to other
multi-bladed propellers, only magnified. A comparison between propellers with a
different number of blades might also depend on the size of the rotors and their rotational speed, but these effects have been neglected. The simulation for the threebladed propeller has been conducted at the center point from Figure 5.3 to get the
most representative comparison. Once again, it was decided to normalize eq and
PL with the results for the reference rotor. The ratios for the three-bladed propeller
can be viewed in Figure 5.6, where the color scale used will be utilized for all of the
the proceeding rotor comparison graphs.

Figure 5.6: eq and PL ratios for an equally sized three-bladed propeller.

90

5.1

Rotor and airframe assessment

Using an equally sized propeller with three blades at the same rotational speed as
the reference propeller will result in an increase of the aerodynamic forces and
moments, including the thrust. A multirotor that uses propellers with three blades
can therefore carry more weight, which is favorable for wind resistance and cannot
be investigated in this simulation. Instead, it has to be investigated with the complete
multirotor simulator.
Alternatively, a scaling factor Ks can be used so that both propellers have equal
thrust margin = Tmax /(mg). The maximum thrust was calculated at 6500 RPM
680 rad/s for each propeller. During this simulation, it was also decided to estimate
the time constant for the motor. This was accomplished by using a linearization of
the DC motor model presented in Section 2.2. Since the reference propellers rotational speed was used for this comparison as well, the motor model was linearized
around 425 rad/s. The resulting eq and PL ratios from the simulation can be seen
in Figure 5.7. The time constant proved to be almost exclusively dependent on the
rotational inertia of the propeller. The time constant for the motor with the scaled
three-bladed propeller was therefore approximately 25 % longer than for the reference.

Figure 5.7: eq and PL ratios for a scaled three-bladed propeller. A scaling factor
Ks = 0.952 was used to obtain the same as for the reference propeller.

Conclusions The results from the two simulations are fairly equal, where the
scaled propeller showed slightly better efficiency and wind resistance performance.
This was expected from results shown in Figures 5.4-5.5. The power loading is
clearly lower than for the reference propeller, whereas a uniform result for the equi91

Chapter 5. Design development

librium angle was not evident. The increased rotational inertia, caused by the extra
blade, decreases the bandwidth of the motor, which is undesirable for the multirotors ability to deal with rapid wind variations.
Despite the disadvantages, there are certain implementation aspects which might
swerve the choice in favor of the three bladed propeller. The thrust for a given propeller radius is higher for the three bladed propeller, enabling an increased aircraft
payload without increasing the overall size of the vehicle. The possibility to carry
additional mass will have a positive effect on the multirotors ability to counteract
wind disturbance. Also, a three-bladed propeller might be favorable from a structural perspective. The reduced tip speeds and more evenly distributed inertia around
the rotational axis could potentially reduce the vibrations translated to the multirotor.
To summarize, pros and cons of the multi-bladed propeller are listed in Table 5.1.
PROs
The maximum possible payload is
increased for a given rotor size.

CONs
Increases the motor time constant.
Inefficient.

Allows a more compact airframe

design.
Potentially reduces vibrations.
Table 5.1: List of pros and cons of the three-bladed propeller.
Variable pitch propellers Using propellers that enable control of the blade pitch
presents a more complex rotor solution. By using such rotors, perhaps the largest
advantage with multirotors, mechanical simplicity, is somewhat neutralized. There
are, however, some major advantages which can compensate for the increased complexity. An implemented variable pitch propeller can be seen in Figure 5.8.

Figure 5.8: Implementation of a variable pitch propeller for the Stingray 500
quadrotor [25].
92

5.1

Rotor and airframe assessment

Conclusions Variable pitch propellers is not a rotor concept that will be modeled
during this project due to time constraints. However, it provides a decreased actuator time constant [8] which primarily increases the agility and responsiveness of
the multirotor. It can potentially also improve wind resistance if a suitable control
algorithm and allocation is applied. It is, however, hard to predict how a multirotor equipped with variable pitch propellers will be affected by wind disturbances
because of the added degree of freedom. Therefore, a complete model has to be
developed and controlled before any easily interpreted results can be achieved.
The added complexity of this rotor concept places the aircraft in a grey zone between helicopter and regular multirotor. Hence, a performance comparison with
helicopters of similar size could also be relevant if variable pitch propellers are considered.
To summarize, pros and cons of the variable pitch rotor concept are listed in Table
5.2.
PROs
Decreases the actuator time constant.
Acceleration during descent flight
can be greater than gravity.

CONs
Increase complexity, both in terms
of control and mechanical development.

More advanced acrobatic maneuvers are possible, including inverted

flight.
Table 5.2: List of pros and cons of the variable blade pitch propellers.

Coaxial rotor
The coaxial rotor setup was introduced as an alternative to the conventional main
and tail rotor setup on helicopters. It utilizes two contra-rotating propellers, with
one placed above the other. Using a coaxial rotor makes it possible to control the
yaw moment while still directing all thrust capacity upwards. The downside with
coaxial rotors is that the lower propeller is working in the down-wash of the upper
propeller, which results in an efficiency loss.
Simulation The two comparison alternatives that were used for the three-bladed
propeller will also be used for the coaxial rotor configuration. To make this possible,
some modifications to the rotor modeling had to be made to account for the downwash of the upper rotor. A more detailed explanation of these modifications are
included in Appendix C.4. During the first simulation, the propellers in the coaxial
93

Chapter 5. Design development

rotor configuration were of the same size as the reference propeller. As a result of
the down-wash from the upper propeller, the lower propeller usually operates at a
higher rotational speed so that the drag moments from the propellers are balanced.
Therefore, the rotational speeds for the coaxial propellers have been chosen so that
they deviate an equal amount from the reference speed. The ratios from the first
simulation can be seen in Figure 5.9

Figure 5.9: eq and PL ratio for a coaxial rotor, where the two propellers are of the
same size as the reference propeller. The upper propeller is running at 410 rad/s
and the lower at 440 rad/s.
For the second simulation, the scaling factor Ks was used to give the coaxial rotor the
same thrust margin as the reference rotor. Since it is usually the lower propeller
that saturates first in a coaxial rotor, it was decided to calculate Tmax while the lower
rotor was running at the maximum rotational speed and the upper rotor at a speed
of 640 rad/s. Thus, cancellation of the yaw moments was ensured. The equilibrium
angle and the power loading from the simulation can be seen in Figure 5.10, where
the propellers ran at the same rotational speeds as in the first coaxial simulation. The
time constant of the coaxial rotor was approximately 50% smaller than the reference
propeller.

94

5.1

Rotor and airframe assessment

Figure 5.10: eq and PL ratio for a coaxial rotor. A scaling factor Ks = 0.906 was
used to obtain the same as for the reference rotor.

Conclusions It was predicted that the coaxial rotor would be less efficient than a
normal propeller, which was also confirmed by the results in Figure 5.9. However, if
the rotor is scaled so that the thrust margin is equal to that of the reference propeller,
the power loading is almost as good in hover and even surpasses the reference for
high wind magnitudes acting at larger disk angle of attacks. There is not a uniform
result for eq , but it can be seen that the coaxial rotor is less wind resistant for small
wind magnitudes while it depends on for large magnitudes.
The coaxial rotor option offers a space-efficient way to generate more payload capacity which if utilized, is an advantage during challenging winds. If the radius of
the coaxial rotor is instead scaled for an equal to that of the reference rotor, the
time constant of the motor is significantly decreased. This is a great benefit for a
multirotor tasked with handling highly turbulent wind.
It is possible to improve the performance of a coaxial rotor by using different settings for top and bottom propellers [19]. However, this is a complicated task and
will not be performed in this thesis work.
Coaxial rotors will make the design process a bit more challenging. The lower propellers will make it harder to integrate the landing gear and the camera equipment.
They will also more easily interfere with the field of view of the camera.
A tandem propeller configuration is also a possible alternative where partly overlapping propellers will not be as inefficient as a coaxial rotor setup, but on the other
95

Chapter 5. Design development

hand will require more space. Tandem propellers can be seen as a compromise between single and coaxial configurations.
To summarize, pros and cons of the coaxial rotor concept are listed in Table 5.3.
PROs
The maximum possible payload is
increased for a given rotor diameter.
Or, the motor time constant is decreased for a given thrust margin.
Allows a more compact airframe
design.
Doesnt necessarily need another
rotor to counter the reaction moment on the airframe.

CONs
Inefficient, especially in hover.
Difficult to model accurately if new
control algorithms are to be tested
virtually.
Difficult to optimize coaxial rotors.
Requires a stiffer frame to handle
the extra force and weight.
More complex integration of landing gear and camera equipment.

Table 5.3: List of pros and cons of the coaxial rotor configuration.

Ducted fan
A ducted fan is a special kind of propeller that is protected by a shroud. The ducted
fans have been used since the 1950s when Alexander Lippisch designed the first
prototype for the VTOL aircraft Aerodyne [38]. However, the principle has been
used in marine applications for even longer. A quadrotor that utilizes ducted fans
can be seen in Figure 5.11.

Figure 5.11: The cyberQuad developed by Cyber technology [26].

Conclusions The ducted fan rotor concept offers an interesting approach that may
provide a steadier air stream through the rotor disk. However, BET predicts that
airflow parallel to the rotor disk is beneficial to rotor efficiency. This effect may be
reduced by using ducted fan rotors. On the other hand, the achievable thrust may not
fluctuate as much, allowing steadier flight control. The shroud will most likely be a
96

5.1

Rotor and airframe assessment

significant drag source at higher wind magnitudes, which will decrease the ability
to resist wind.
To summarize, pros and cons of the ducted fan rotor concept are listed in Table 5.4
[38].
PROs
The maximum possible payload can
be increased for a given rotor diameter under certain circumstances.

CONs
The duct will be a source of aerodynamic drag, which will be significant at higher wind magnitudes.

The shroud acts a structural barrier that protects the blades from the
surroundings and vice versa.
Table 5.4: Ducted fan rotor list of pros and cons.

Inverted propeller placement

A slightly alternative configuration of a normal propeller is to mount the propeller
and its motor beneath the multirotor arm. By doing so, the down-wash of the rotor
will no longer be interfered by the DC motor and the arm, which it is mounted on.
Conclusions The airframe and the motors disrupt the induced air flow passing
down through the rotor. This effect has not been modeled in this thesis work, but
realistically it should be beneficial from an efficiency perspective to place the propeller beneath the arm. The reasoning behind this assumption is that it would be
advantageous to interfere with the rotor slipstream before it has passed through the
rotor, where the airflow speed isnt as high or concentrated.
Also, by placing the propeller below the arm and the motor, the rotor plane is lowered relative to the CoG of the aircraft. According to the multifactor analysis in
Chapter 4, this would be advantageous from a wind resistance perspective. A somewhat marginal side effect of placing a propeller below the motor is that the motor
cooling, in the form of the induced airflow, is slightly reduced. Additionally, making
an optimal design of the multirotor will be more challenging, considering landing
gear and camera equipment integration.
To summarize, pros and cons of an inverted propeller placement are listed in Table
5.5.

97

Chapter 5. Design development

PROs
Potentially more efficient.

CONs
Less efficient motor cooling.

Easier to place the rotor plane lower

relative to the CoG.

More complex integration of landing gear and camera equipment.

Table 5.5: List of pros and cons of the inverted propeller rotor configuration.

Tiltable rotor
Another possible rotor configuration is to implement a servo motor on the multirotor
arm, enabling control of the rotor tilt angle. This will add additional degrees of
freedom to the multirotor system. Commonly, this is used for one of the rotors on
a tricopter since three normal rotors isnt enough to fully control the required four
degrees of freedom. Also, it is possible to implement an actuator system so that the
rotor can tilt around two axes, which would further increase the control possibilities.
However, this would require a very complex system design.
Conclusions If a multirotor were to use tiltable rotors for all of its arms, forward
flight and wind compensation can be achieved without pitching or rolling the body.
This has been implmented on a quadrotor prototype by Markus Ryll et al. [39],
which lead to promising results. Not necessarily having to tilt the multirotor body
is generally a very good thing concerning wind resistance, since the airframe drag
can be minimized in the horizontal directions. It would most likely also reduce the
fluctuation of the drag force acting on the airframe since the exposed area and its
drag coefficient will be more or less constant. In addition, it should also be possible
to increase the responsiveness of the multirotor while trying to compensate for an
applied drag force, depending on the time constants of the servo motors. Tilting the
rotors instead of the entire body should also provide a steadier platform for aerial
video and photography.
To summarize, pros and cons of the tiltable rotor configuration are listed in Table
5.6.

98

5.1

Rotor and airframe assessment

PROs
Adds additional degrees of freedom.

CONs
Servo motors add to the cost and
complexity.

Increased maneuverability.

Multirotors with tilting rotors needs

more complex control algorithms.

Forward flight possible without rotation of the aircraft body.

Difficult to model accurately if new

control algorithms are to be tested
virtually.

Table 5.6: List of pros and cons of the tilting rotor configuration.

Varying rotor sizes

Using propellers with varying dimensions is an unconventional approach to multirotor design. It has been successfully implemented by [9], where the aim was to
increase the flight time for operations close to hover.
Conclusions Using differently sized rotors may prove beneficial depending on the
purpose of the design. If theory predicts an improved performance with a large main
rotor, as was the case in [9], then it can certainly be worth the added implementation
difficulties. However, considering the aim of this thesis and the previously mentioned conclusion that smaller and faster spinning propellers are desirable, there is
no evident reason to use anything but small propellers of the same size.
To summarize, pros and cons of utilizing differently sized rotors are listed in Table
5.7.
PROs
Greater design flexibility.
Can improve efficiency under certain circumstances.

CONs
Harder to balance yaw moment and
more unknown constants are introduced in the control algorithm.
Varying motor time constants further complicates the control algorithm implementation.

Table 5.7: Varying sizes of rotor list of pros and cons.

99

Chapter 5. Design development

Airframe configurations
During this section, different airframe configurations will be discussed and compared. For simplicity, only one type of rotor will be considered during this study.
The airframes will be sorted after their conventional number of rotors.
Quadrotor
This part will cover multirotor airframes that normally utilizes four rotors.
I4/X4 config. The I4/X4 quadrotor configuration, shown in see Figures 5.12 and
5.13, is the most common airframe option for multirotors. It will here serve as an
overall design reference to which other designs concepts will be compared. The
quadrotor can be used either with an I configuration or a X configuration depending
on how the arms are oriented relative to the front of the multirotor. The different
configurations have slightly different flight characteristics for a given flight direction, but will not have any significant impact on the overall performance of the
multirotor. The X configuration is better suited for aerial photography.

Figure 5.12: Quadrotor: I4 config.

Figure 5.13: Quadrotor: X4 config.

H4 config. The H4 quadrotor is an alternative frame configuration and is shown

in Figure 5.14. It can be seen that the arms do not coincide in a single point as for
the standard quadrotor.

Figure 5.14: Quadrotor: H4 config.

100

5.1

Rotor and airframe assessment

Conclusions
The H4 provides some benefits for specific applications such as
aerial photography. However, it will be slightly harder to optimize the CoG placement which can reduce the stability of the vehicle. The H4 configuration will also
have slightly different flight characteristics, but it is hard to predict how much of
an impact this will have on the overall performance. Drag forces that are caused
by the airframe in the x- and y-directions of the body frame will differ due to the
rectangular design. Therefore, it might be worth implementing a yaw control which
forces the front of the quadrotor to face the incoming wind in certain situations.
Pros and cons of the H4 configuration, compared to the standard quadrotor configuration, are listed in Table 5.8.
Pros
Better for aerial photography.
More space to place battery and
equipment.

Cons
Slightly more difficult to balance
the center of gravity.
Behaves differently depending on
the direction of the wind.

Table 5.8: H4 list of pros and cons.

V-Tail config. The V-Tail is an unconventional quadrotor design with only three
arms, as shown in Figure 5.15, which can potentially reduce the airframe drag. It
utilizes a rear rotor configuration with two tilted rotors to achieve yaw control.

Figure 5.15: Quadrotor: V-Tail config.

Conclusions The V-Tail is a slightly more complex configuration than the regular quadrotor configurations. It is a good choice if more responsive yaw control is
required, which is accomplished by the two rotors that are tilted with a fixed angle.
This means that yaw control moment isnt limited to the reaction moments from the
DC motors. The drawback of this solution is that it reduces the potential upward
thrust of the multirotor. This issue may also be worsened by the interference of the
rear rotors.
101

Chapter 5. Design development

Pros and cons of the V-tail configuration, compared to the standard quadrotor configuration, are listed in Table 5.9.
Pros
Improved yaw control.
Potentially reduced airframe drag.

Cons
Potential thrust is limited by the
tilted rear rotors.
Rear rotors will most likely interfere with each other which results
in an efficiency loss.

Table 5.9: V-tail list of pros and cons.

Trirotor
This part will cover multirotor airframes that normally utilizes three rotors. However, this is not enough to provide the necessary means of control. For stable flight,
at least four DoF is required which is commonly achieved by implementing a servo
motor on one of the arms.
Y3 config. The Y3 tricopter is the standard configuration for trirotors. It can be
viewed in Figure 5.10, where the arrow around the rear arm indicate the use of a
servo that ensures yaw control.

Figure 5.16: Tricopter: Y3 config.

Conclusions With only three arms there is a possibility to make the design more
compact because of the increased spacing between the rotors. The drag force acting
on the airframe can therefore be reduced. It is also predicted that the drag coefficient will be reduced as a result of fewer arms. Naturally, the trirotor cant handle
much payload because of the limited amount of rotors. However it can be a suitable
platform for small-scale camera equipment since the unobstructed field of view is
increased. The trirotor is popular among hobbyists since it has a different feel
while flying manually.
102

5.1

Rotor and airframe assessment

Pros and cons of the Y3 configuration, compared to the standard quadrotor configuration, are listed in Table 5.10.
Pros
Reduced airframe drag.
Less likelihood of rotors being affected by one another and by the
airframe itself.
Only three three arms leaves space
for a wide camera angle.

Cons
Either servo or additional rotors are
necessary to obtain full controllability.
Less stabilization authority.
Poor payload capabilities.

Table 5.10: Y3 list of pros and cons.

T3 config. The T3 tricopter is an alternative trirotor configuration, which is depicted in Figure 5.17. The arms doesnt coincide in the CoG as in the case of the Y3
configuration.

Figure 5.17: Tricopter: T3 config.

Conclusions Just as with the H4 configuration, this concept is mainly beneficial

for special applications. The T3 configuration is for example advantageous if servo
motors are implemented on all of the arms, since no rotors will counteract each
other. In theory, however, this concept will be considered as slightly inferior to the
standard Y configuration because of its unsymmetrical design.
Additional pros and cons to those of the standard trirotor Y3 configuration are listed
in Table 5.11.

103

Chapter 5. Design development

Pros
Even wider camera angle during
aerial photography.

Cons
Slightly more difficult to balance
the center of gravity.
Behaves differently depending on
the direction of the wind.

Table 5.11: T3 list of pros and cons.

Hexacopter
This part will cover multirotor airframes that normally utilizes six rotors.
I6/X6 config. As for the quadrotor, the hexacopter can also be designed with either
an I or a X arm configuration, see Figure 5.18. Similarly, a slight difference in flight
characteristics is present due to the arms placement relative to the front of the
aircraft. This difference is even more marginal for the hexacopter since the total
thrust is distributed more evenly around the center of the multirotor.

Figure 5.18: Hexacopter: X6 config.

Conclusions
Hexacopters are generally used when quadrotors arent strong
enough or when some limited redundancy is required because of expensive equipment. The extra payload capacity comes at the price of increased drag forces on the
airframe. The airframe can not be made as compact as those with fewer arms due
to rotor spacing. The projected surface area will therefore increase and the effective drag coefficient will most likely also increase because of the increased number
of arms that the wind has to bypass. It is, however, hard to analytically determine
exactly how much the overall drag will increase.The additional motors and arms
will increase the overall mass and inertia of the multirotor which will increase the
translational and rotational stability, but somewhat limit the agility of the aircraft.
Pros and cons of the I6/X6 configuration, compared to the standard quadrotor configuration, are listed in Table 5.12.
104

5.1

Rotor and airframe assessment

Pros
Limited redundancy.

Cons
Increased airframe drag.

Larger payload capacity.

Generally not as agile.

Generally more stable with respect

to rotation and translation.
Table 5.12: I6/X6 list of pros and cons.

Octorotor
This part will cover multirotor airframes that normally utilize eight rotors.
I8 config. The I8 octorotor is the standard configuration for airframes with eight
rotors and can be viewed in Figure 5.19. The octorotor is generally used when the
requirement on payload capacity is highly demanding. The chances of a safe return
in case of motor failure is also significantly increased. Pros and cons of the I8 configuration, compared to the standard quadrotor configuration, are basically the same
as those for the hexacopter, only magnified. Hence, no new list and conclusion is
included.

Figure 5.19: Octorotor: I8 config.

105

Chapter 5. Design development

V8 config. The V8 configuration is an unconventional octorotor design where the
rotors are distributed evenly on either side of the center, see Figure 5.20.

Figure 5.20: Octorotor: V8 config.

Conclusions The V8 is a slightly more complex design with some advantages for
special applications. It is beneficial for larger camera equipment, thanks to the payload capabilities and the wider unobstructed field of view. As mentioned previously,
the unsymmetrical design will make it slightly harder to balance the CoG properly.
The vehicle will also react differently to wind from various directions.
Additional pros and cons to those of the standard octorotor configuration are listed
in Table 5.13.
Pros
Wide camera angle.

Cons
Slightly more difficult to balance
the center of gravity.
Behaves differently depending on
the direction of the wind.

Table 5.13: V8 list of pros and cons.

106

5.2

5.2

Concept evaluation

Concept evaluation

This section aims at presenting a favorable multirotor concept able to withstand

wind. This is achieved by combining the previously presented rotor and airframe
assessments, after which a concept generation and scoring process is performed. To
this end, a number of chosen design criteria are given a relative weighting which
are used to rate a number of concepts generated by the authors.

Design criteria
Before any concepts are presented, a set of criteria is established for which each
concept will be rated on a scale of 1 to 10. Each criteria is accompanied by a relative weight to signify its importance to the aims of this thesis. Apart from wind
resistance, flight time was also considered an important factor for mission performance. This is reflected in the moderately high weighting of efficiency.

Cost and complexity [0.05]

Rates the implementation simplicity and availability of the airframe design.
Issues such as CoG balancing are considered as well as the required amount
of material and components such as arms, motors and cabling. Along with
these considerations, the modeling and control complexity is also taken into
account. A low rating on this criteria indicates high levels of cost and complexity. Since it is assumed that any possible production of a proposed design
will be small-scale, both cost and complexity will be given roughly equal
significance.
Wind resistance [0.45]
Considers how favorable a design is regarding its ability to withstand and
minimize forces and moments caused by the onset of wind. Typically, the
projected surface area, aerodynamic drag coefficient, number of rotors and
payload capacity will have an influence on this criteria.
Maneuverability [0.10]
The agility of the aircraft along the six degrees of freedom is rated by this
criteria. Increased mass and moment of inertia will limit the ability to perform aggressive maneuvers. However, additional rotors around the CoG will
increase the attainable thrust and actuated moments, improving maneuverability. Adding tilting rotor servos will open up further control possibilities
not otherwise achievable.
Efficiency [0.25]
Increased efficiency of the aircraft is pivotal to maximizing the flight time.
Larger rotors operating at low speeds are generally more efficient, whereas
107

Chapter 5. Design development

coaxial rotors will suffer from some efficiency loss. Other factors that influence this criteria are the unloaded aircraft weight and the drag exposure,
including number of rotors and airframe design.
Redundancy [0.10]
The ability to maneuver the aircraft in the case of motor or rotor malfunction is considered by this criteria. Redundancy may be achieved by use of
additional rotors as well as placing them appropriately around the airframe.
Designs with several rotors will not suffer as much thrust margin loss as those
with fewer in case of rotor malfunction. Additionally, designs with fewer rotors may even suffer loss of control caused by an insufficient number of actuated degrees of freedom. It should, however, be kept in mind that an increase
in rotors will also increase the likelihood of rotor failure.
Aerial photography [0.05]
For photography and video capturing purposes, it is important that the propellers and arms do not interfere with the camera view. This criteria rates
the ease of a satisfactory camera system implementation on the multirotor, as
well as its ability to carry additional payload introduced by such systems.
To make the following concept comparison as fair as possible, a number of constraints have been introduced. These are listed below.
All of the concepts will have the same thrust margin = Tmax /(mg).
The minimum horizontal distance between the tips of any two propellers are
limited to R to avoid unwanted rotor wake interaction [16].
All multirotor concepts, that use equally sized propellers, will have the same
propeller dimensions.
All concepts will use the same control algorithm. The response time of the
control algorithm is determined by the size of the rotor so that a tuning factor
N is achieved relative to the actuator response time (See Chapter 3).
With design criteria defined, the following sections will present the generated concepts. Each concept is described thoroughly and given a rating for each of the design
criteria.

108

5.2

Concept evaluation

X4 quadrotor [X4]
The quadrotor is the traditional multirotor configuration concept. It is mechanically
simple and proven, and will serve as a reference concept during the concept scoring.

Figure 5.21: The X4 concept is a standard quadrotor configuration, used as a reference concept.
Cost and complexity [10] The quadrotor has had great success thanks to its relatively low cost and complexity. Four is the lowest amount of rotors needed for sufficient controllability, without adding additional complexity in the shape of servo
motors. Since the quadrotor is the most common platform, it is easy to obtain various information regarding design considerations and control.
Wind resistance [5] The airframe of the quadrotor has four arms, yielding a compact design with moderate airframe drag. Four rotors provide sufficient thrust for a
smaller payload, granting some stability and resistance to the onset of wind.
Maneuverability [5] The conventional quadrotor can be quite agile due to its
relatively low inertia, while only using four rotors.
Efficiency [7] The quadrotor utilizes regular rotors, which provide good efficiency. It will have a rather low unloaded weight and unnecessary drag exposure
thanks to few components such as arms and motors.
Redundancy [1] If a rotor is lost, crash is virtually unavoidable. A fourth of the
total lifting force as well as control of both attitude and position is lost, highlighting the severity of such a situation. Therefore, the quadrotor is not considered a
redundant system.
Aerial photography [5] The quadrotor is quite suitable for aerial video and photography. However, integrating heavier equipment in the total mass may require
more payload capability.
109

Chapter 5. Design development

Y6 three-armed frame with coaxial rotors [Y6]

The Y6 concept features a three-armed airframe with coaxial rotors attached to each
arm and is known among hobbyists for its wind robustness. The coaxial rotors provide means to achieve yaw control, which eliminates the need for a servo controlled
rotor or tail rotor.

Figure 5.22: The Y6 concept has a three-armed airframe and coaxial rotors.
Cost and complexity [5] The coaxial rotors of this concept has some drawbacks
when it comes to both cost and complexity. Since the upper and lower rotors provide
different thrusts and drag moments, it may be slightly more difficult to optimize the
control and design of the rotors. In addition, coaxial rotors are more difficult to
model accurately in a simulation environment. The three-armed airframe requires
less material, whereas the coaxial rotors will introduce additional motors and propellers. It should be noted that the airframe arms will have to be made slightly
stiffer than otherwise and that the utilization of coaxial rotors will complicate the
integration of landing gear.
Wind resistance [9] The airframe used in this concept will be of some benefit
compared to the X4 reference. Due to fewer and shorter arms, it is possible to make
it more compact and less susceptible to airframe drag. Utilizing coaxial rotors will
enable additional payload capacity useful in wind, but it should be kept in mind that
each added rotor is a source of additional drag.
Maneuverability [5] The Y6 suffers from a slight increase in inertia due to added
motors, but this is somewhat compensated by a reduced arm length. Its ability to
produce rolling and pitching moments is, however, increased by the presence of
additional rotors. Thus, it is hard to arbitrarily determine how the maneuverability
of the Y6 compares to the X4. Hence, it was concluded that this concept would
receive the same rating as the X4 reference.
Efficiency [3] The unloaded weight of the Y6 will likely be slightly increased
compared to the X4 reference, although the Y6 has one arm less. This is due to the
110

5.2

Concept evaluation

added motors and rotors. However, the main efficiency loss is caused by the coaxial
rotors. Since the lower rotors operate in the wakes of the upper, more power is used
per generated thrust. Additionally, the airframe is likely to be exposed by less drag,
yielding some benefits in forward flight.
Redundancy [7] Losing a rotor will essentially result in a significant decrease in
attitude control performance. Depending on the predetermined thrust margin and
flight conditions, the aircraft may still be landed safely either manually or by an
emergency control allocation routine. The lost thrust due to rotor failure will be
limited to less than a sixth of the total thrust since three of the rotors will always
operate without downwash disturbance.
Aerial photography [6] The three-armed airframe leaves a lot of space, which allows a wide field of view for the camera without interfering propellers. The coaxial
rotors, although providing some additional payload capacity for heavier equipment,
may cause some interference of the view depending on where the camera is directed.

X6 hexacopter [X6]
This concept employs a six-armed frame configuration with a single rotor placed on
each arm. It is similar to the design that Combine Control Systems currently use.
The six rotors provide additional payload capacity, but placing them with sufficient
spacing yields a larger airframe.

Figure 5.23: The X6 concept has a six-armed airframe with rotors placed beneath
the arms.

Cost and complexity [7] The X6 suffers from high material cost. It has a large
airframe due to additional arms, which have to be long enough to provide necessary
rotor spacing. Additionally, the six rotors introduce an increased number of motors
and propellers.
111

Chapter 5. Design development

Wind resistance [6] Six rotors provide a larger payload capacity, which stabilizes
flight in windy conditions. However, the airframe of the X6 will be significantly
larger than those of the previous concepts, which in combination with added rotors,
results in an increased overall drag.
Maneuverability [5] Similarly to the discussion on the Y6 maneuverability, the
X6 will receive the same rating as the X4 reference. The concept does provide
means to produce larger attitude moments, but also has a larger inertia an increased
number of longer arms.
Efficiency [6] Since the regular single rotor configuration is employed, the X6
will be rather efficient. However, the unloaded weight of the X6 is increased due
to the six-armed frame and the addition of extra components such as arms and motors. The large airframe of this concept will likely be a source of significant drag,
requiring extra thrust to compensate for during forward flight.
Redundancy [6] Similarly to the Y6 concept, rotor loss will essentially result in a
significant decrease in attitude control performance, but with some additional thrust
loss.
Aerial photography [4] The presence of additional arms, which will generally
be longer, and propellers will cause some obstruction of the camera view in forward flight. The added payload capacity will, however, enable integration of heavier
equipment.

X4 quadrotor with servo-tiltable rotors [X4S]

The last proposed concept has a simple quadrotor airframe, but with the added functionality of servo-tiltable rotors. The servos gives additional actuation options for a
control algorithm and enables simultaneous control of all six degrees of freedom.
An important aspect to consider for this concept is the increased complexity.

Figure 5.24: The X4S with tiltable rotors concept has a standard quadrotor frame
but with rotors that can be tilted by DC servos.
112

5.2

Concept evaluation

Cost and complexity [3] Although the airframe itself is of the simple quadrotor design, the introduction of servo-tiltable rotors yields a considerable increase
in both cost and complexity. The added servo motors themselves add to the material cost, whereas mechanical implementation add to the complexity. In addition,
extending the multirotor model to include the servo functionality and developing a
suitable control algorithm may prove difficult. With these considerations in mind,
the X4S concept is rated rather low for this criteria.
Wind resistance [8] As for the regular quadrotor, this concepts airframe will
result in moderate airframe drag during the onset of wind. Since servo-tiltable rotors are used, the airframe will mainly be exposed to wind in its horizontal plane.
Utilizing this, the horizontal properties of the airframe can be aerodynamically optimized to decrease the overall drag. In addition, the tilting of the rotors can be
actuated faster than rotation of the entire airframe, provided the server motors are
sufficiently powerful. This will therefore yield a faster response to rapidly varying
wind conditions.
Maneuverability [8] Servo-tiltable rotors allow simultaneous control of all six
degrees of freedom, enabling maneuvers not possible for regular multirotors. Also,
as previously discussed, tilting just the rotors to achieve forward flight may be done
faster than rotating the entire airframe like fixed-rotor concepts do. With these considerations in mind, the X4S concept is considered highly maneuverable.
Efficiency [5] Despite the presence of added unloaded weight from the servo motors, the X4S may still provide average efficiency thanks to its ability to tilt its rotors.
Positive effects resulting from this ability are most evident in forward flight, where
the airframe drag may be minimized and the need to rotate the entire airframe is
eliminated. However, it should be kept in mind that the added servo actuators continually require power to control the attitude.
Redundancy [1] In the even of rotor failure, attitude control will likely be lost as
well as a significant amount of the total thrust. Thus, the X4S is not considered to
be redundant with crash being an inevitable outcome. In case of servo motor failure,
the aircraft may still be operated safely although some emergency control routine
reallocating servo actuation to rotor speed would be needed.
Aerial photography [6] The quadrotor frame is quite suitable for aerial video
and photography, where the X4S draws further benefit by using its servo tiltable
rotors. These will enable forward flight without rotation of the airframe, thus reducing obstruction from the arms as well as the required effort by a camera stabilizer.
However, integrating heavier equipment in the total mass may require more payload
capacity than the X4S can provide.
113

Chapter 5. Design development

Results
The results from the concept scoring are presented in Table 5.14. The criteria ratings
for each generated concept are denoted as R, which are multiplied by a weighting
factor to yield an overall weighted score WS.

Criteria

Weight

X4
R

WS

Y6
R WS

X6
R WS

X4S
R WS

Cost & complexity

Wind resistance
Maneuverability
Efficiency
Redundancy
Aerial photography

0.05
0.45
0.1
0.25
0.1
0.05

10
5
5
7
1
5

0.50
2.25
0.50
1.75
0.10
0.25

5
9
5
3
7
6

7
6
5
6
6
4

3
8
8
5
1
6

Total score
Rank

5.35
4

0.25
4.05
0.50
0.75
0.70
0.30

6.55
1

0.35
2.70
0.50
1.50
0.60
0.20

5.85
3

0.15
3.60
0.80
1.25
0.10
0.30

6.20
2

Table 5.14: Concept scoring for the four generated design concepts.

5.3

Proposed design

From the concept scoring results in Table 5.14, it is clear that the Y6 is a preferred
concept along with the X4S. It was concluded that the proposed design concept of
this thesis will be a fusion of the Y6 and X4S concepts, combining the strengths
of each concept to provide a preferable design. Such a fusion design may look
like in Figure 5.25. Considering that the UAV will constantly be in forward flight
for the suggested applications, adding a wing-like airframe structure granting additional lift might be a consideration to improve flight time. Since it was deemed too
time-consuming to model and control the key features of both the Y6 and the X4S
accurately it was, however, decided to move forward with just the Y6 in this thesis
work.
The Y6 concept employs a coaxial rotor configuration. Although proving to be less
power-efficient, such rotors enable a smaller airframe with fewer arms while also
providing about 20 % extra payload capacity compared to a standard quadrotor.
Since the three-armed airframe can be also made smaller than for the X4 and X6
concepts, a larger portion of the mass can be placed in the center. It can be argued
that the center mass is the most useful mass, since it is typically where equipment
such as onboard electronics and battery is carried.
114

5.3

Proposed design

Figure 5.25: Combining the beneficial features of the Y6 and X4S concepts may
yield a superior design.
The DOE multifactor analysis in Chapter 4 indicated that a reduced airframe size
will have positive effects on both wind resistance capabilities and energy consumption. This suggests that a compact design is beneficial. However, the multirotor can
not be made arbitrarily small since it needs a reasonable arm length to actuate moments for attitude control. Placing the rotors in too close proximity to each other
will also cause interaction between the rotor slipstreams [16]. To stay clear of this
unwanted effect, it was concluded that the arms must be long enough so that the
distance between the blade tips of any two rotors is at least one rotor radius. For
different frame configurations, this will lead to different lower limits of the arm
lengths. Because of this, the three-armed frame of the final Y6 will have shorter
arms than otherwise possible for the X4 or X6 designs, yielding reduced projected
surface areas resulting in lower exposure to wind drag force.
Also, the DOE showed that the total mass had great significance to withstanding
wind. The downside of a high mass was shown in regards to the power consumption
and, obviously, the mass can not be made too large anyway since the thrust from the
rotors must be enough to lift the aircraft and perform maneuvers while maintaining
altitude. It was concluded that the total mass of the aircraft mt ot must respect a
minimum thrust margin of , so Tmax mtot g.
Lastly, another finding of the DOE analysis showed that the points at which the
airframe drag and the rotor forces act are significant to minimizing the position
error. Attempting to place them as close to the center of gravity as possible will
minimize any moments, since the levers of these are reduced. But even more benefit
can be achieved if these points are situated slightly below the center of gravity. The
forces would then act to tip the aircraft slightly into the wind, directing a portion of
the thrust against the wind forces without any controller action. To take advantage
of these effect, the arms are attached as low as possible on the center.
To summarize, combining the final Y6 design concept with the findings from the
DOE analysis yields a favorable design. The desired reduction of the airframe size
115

Chapter 5. Design development

is clearly accommodated by the Y6s three-armed frame due to fewer and shorter
arms. The six rotors, configured coaxially, enable some extra total mass although the
extra motors and rotors will be sources of additional airframe and rotor drag forces.
A possible final product may look like in Figure 5.26, which has been designed with
CAD and imported into SimMechanics.

Figure 5.26: Y6 configuration shown in the graphical interface of SimMechanics.

116

6
Performance comparison
In an attempt to demonstrate the advantages of the proposed Y6 design from Section
5.3, a Y6 and two X4 configurations have been assembled in the model simulator in
order to compare their respective performance. In Table 6.1, these designs are listed
with their most relevant specifications. The two X4 designs represent the parameter
combinations that provided the best wind rejection (X4W) and power consumption
(X4P) in the multifactor optimization from Section 4.1. They are compared to the
Y6 design, which is configured to have the same thrust margin as the wind resistant X4W. Along with this, a dimensioning principle has been applied to ensure
equal rotor spacing. Thus, the Y6 and X4 will have different airframe projected
surface areas where the same drag coefficient C has been assumed.
Specifications

Y6

X4W

X4P

Total mass

1.80 kg

1.50 kg

1.0 kg

Projected surface areas

Ax , Ay , Az

0.02, 0.02
0.03 m2

0.02, 0.02
0.04 m2

0.02, 0.02
0.04 m2

Thrust margin

1.83

1.83

2.75

Table 6.1: Design specifications for Y6, X4W and X4P.

The designs listed in Table 6.1 are run in the simulator with a wind scenario designed to test each designs ability to handle shorter and longer wind gusts, rotation
of the horizontal wind angle and turbulence. The wind test is shown if Figure 6.1,
where each line shows the wind in one of the inertial directions. When running
the designs through the simulator, the cascaded PID controller has been tuned such
that the positional control has a bandwidth BW 1 rad/s, enabling basic levels of
turbulence and wind variation robustness.
117

Chapter 6. Performance comparison

Figure 6.1: Wind test designed to test the performance of each multirotor design
listed in Table 6.1. Note the legend showing which wind direction corresponds to
which line color.
Resulting simulation performances for each design are summarized in Table 6.2 and
in Figure 6.2. The simulations shows that the Y6 has an approximate performance
increase of 10 % compared to the X4 optimized for wind resistance in regards to
holding its position. This is put down to a number of factors. Firstly, its airframe
is reduced, resulting in a smaller projected surface areas. With the same drag coefficients C assumed for the designs, this is clearly a favorable attribute. Secondly,
the mass is increased, decreasing the acceleration due to the onset of forces. Due to
the increased number of rotors of the Y6, there will be additional rotor disturbance
forces, but this is what enables the added mass. Lastly, the effective levers to the
centroid and resultant lever of the rotors is reduced greatly because of the construction with motors and rotors attached beneath the arms. Note that all of these effects
were showed to be influential in the DOE analysis in Chapter 4 and discussed at
length in Section 5.3.
Also worth noting from the performance comparison in Table 6.2 is the differences
in both power consumption and mass lifting efficiency. The X4P design, which
was shown to be optimal in regards to power-consumption, has superior efficiency
performance. Thus, it is clear that resistance comes at a price of efficiency, where a
suitable trade-off must be made.
Since the simulator doesnt accurately predict drag forces, as discussed in Section
2.4, and the fact that it would be difficult to achieve the same levels of control and
sensor feedback as in this idealized simulator, the results will likely not be the same
as if tested in reality. Especially the uncertainty that arises when comparing different
airframes is a major drawback of the comparison performed in this section, as the
118

Chapter 6. Performance comparison

Performance

Y6

X4W

X4P

Max pos. error emax

Mean pos. error emean
Mean power consumption Pmean
Power per mass Pmean /mtot

8.14 cm
3.34 cm
102 W
57 W/kg

9.30 cm
3.79 cm
76 W
51 W/kg

14.81 cm
5.61 cm
37 W
37 W/kg

Table 6.2: Design performance comparison.

Figure 6.2: The position error norms for the proposed Y6, X4W and X4P designs
during the wind performance test.
real drag coefficients were not possible to predict. Thus it is unclear how much of
an advantage the smaller Y6 airframe has in regard to drag forces and how much
the rotors downwash affect this. Although it is dangerous to jump to conclusions
on this issue, it is believed by the authors that that a smaller airframe with fewer
arms will result in less drag force.
In order to achieve a more accurate comparison, tools not available to and methods
out of the scope of this project must be applied. These could possibly include windtunnel testing and CFD analysis of the airframe especially, but also the rotors for
horizontal airflow. Applying methods like these, further conclusions may be made
about the advantages of each designs. Also, modeling of sensors, real-time implementation aspects and other non-ideal effects not investigated by this thesis work,
will likely increase the realism of the results.

119

7
Summary and conclusions
As discussed when formulating the aims of the thesis work, it has been attempted to
provide a favorable multirotor UAV design for autonomous missions in challenging
wind condtions. The goal was to be achieved using a model- and control-based design approach. The design of such an UAV is of great interest to the overall desire
to design multirotors suitable for aiding in search and rescue missions. Attaching
a thermal camera to the UAV would provide a powerful tool for such applications.
During missions when, for instance, searching for a missing person in the woods, it
is considered unacceptable that worsening weather conditions jeopardizes the ability to use the multirotor.
To these ends, an advanced aerodynamic multirotor model, based on theory normally used for modeling of helicopters [18, 37], has been assembled and implemented in a simulator environment. The model attempts to predict thrust and drag
moment variations in different flight states (hovering, vertical and forward flight)
by combining an air inflow model based on momentum theory with blade element
theory, where thrust, horizontal force and drag moment is derived and integrated
for each blade section. In addition, the blade flapping phenomena was modeled to
account for the tilting of the rotor in both the longitudinal and lateral directions,
deflecting the rotor forces.
Using data from the UIUC propeller database [4] some validation could be made
for the axial and hover flight states, which, after some parameter tweaking, showed
promising agreement. To validate the rotor model further, it is suggested that a series
of wind tunnel tests are performed for forward flight. Without further validation, the
model must was still considered sufficient since it is based on well regarded theory.
It is anticipated by the authors that it at least captures the main flight dynamics of
the rotor. The airframe model, used to apply drag forces, is, however, of low validity
since it was not possible to decide the drag coefficient Cd for a complex multirotor
body. Compared to the rotor forces, an aerodynamically favorable airframe will
likely have a limited contribution to the total drag acting on the aircraft.
Along with the actual aerodynamical multirotor model, modeling of wind characteristics have been of central interest. Wind gusts and variations have been modeled
120

Chapter 7. Summary and conclusions

and implemented using the Dryden wind turbulence model and standardized mathematical gust representations. Using this to create suitable wind test scenarios, it
was decided to not make the wind too strong since the rotor model validation indicated less data compliance for high rotor inflow speeds. Considering this as well
as the low validity of the airframe model, it was decided to not conduct a maximum wind test to test the upper limits of a multirotor design. Instead, evaluating
designs for moderately high winds was decided to be sufficient to draw conclusions
regarding performance at even larger wind velocities. When dimensioning a final
multirotor, the required wind tolerance must be decided, which includes choosing a
suitable thrust margin and arm length to avoid motor saturation. To further prevent
saturation, a more refined control algorithm may be implemented that, for example,
temporarily allows some positional deviation during the most extreme wind gusts.
Finalizing the simulator, a PID cascaded controller was designed and tuned to control the position of the multirotor. The overall control architecture was designed
such that an outer controller controls the position, sending references to the inner
attitude controller. A control allocation scheme was then used to actuate the actual
motors in the simulator. Analysis of the control performance was achieved by using
a simplified and linearized model for control design. Disturbance analysis has then
been made where the closed-loop responses of the linearized model are analyzed.
The Bode plot bandwidths and response times of these responses were discussed in
relation to the modeled wind turbulence, predicted by the Dryden wind turbulence
model.
The modeled simulator was made such that dimensions and masses of various multirotor parts are scalable. Enabling this, an optimization software was used to access
the simulator externally, running it with a specified wind test for different parameter
setups. After a large number of simulations (using different setups determined by a
design of experiments), each parameters influence on position error and power consumption were determined by linear and quadratic regression models respectively.
The results indicated a trade-off between these criteria.
Finally, a design development methodology has been applied to generate a suitable design concept. A wide variety of airframe and rotor configurations have been
assessed from which a small set of multirotor concepts were generated with a predetermined set of criteria in mind. It was concluded that the wind resistance criteria
was of utmost importance, followed by power efficiency, when generating these
concept. The generated concepts were compared through a concept scoring process were the final score was based on the results from the multifactor analysis, the
rotor and airframe assessments, logical reasoning and inspiration from on-going research and community discussions. It was concluded that a three-armed multirotor
with six propellers configured coaxially (Y6) would be the most beneficial, possibly
with the modification of added servo-tiltable rotors. In the conducted simulations,
121

Chapter 7. Summary and conclusions

the different concepts modeled showed limited impact on performance and may in
reality play a larger role, especially in regards to airframe drag. It will invariably
be the physical properties of a multirotor that defines its flight behavior, whereas
different concepts merely provide means to manipulate these properties. Choosing
a specific concept introduces both possibilities and boundaries to how much the
physical properties may be altered.
Because of the uncertainties in the model, the ability provide accurate performance
data is limited. An attempt was however made, showing the superiority of the final
design within the simulated environment. It should be stated that this project was
performed on a conceptual level, not considering implementation aspects such as
sensor readings and real-time constraints, which will likely be of significance in
reality. It is, however, assumed that the evaluated designs in this report are affected
similarly by these aspects, thus not affecting their relative performance.

122

8
Future work
This project has touched on a number of different engineering principles such
as control, mechanics, electronics, aerodynamics, design development and multiobjective optimization. Hence, there are also a lot of possible directions in which a
succeeding project might proceed.
The multirotor model can be further enhanced by including additional effects such
as feathering and tip losses. Perhaps the largest model improvement would be
gained by further studies of how the total drag on the multirotor is affected by the
interaction between the rotors and the airframe. This should realistically be made
after a design concept has been chosen. Alternatively if a better estimation of the
possible flight time is desired, a more detailed model of the battery would be needed.
Only the cascaded PID controller, which is probably not the most suitable for rejection of demanding disturbances and nonlinearities, has been covered in thesis work.
Therefore it is recommended that an in-depth study of different control algorithms
such as backstepping, L1 -adaptive control, quaternion based control or feedback
linearization is conducted. The proposed multirotor model would provide a suitable test platform for such an evaluation. It would also be interesting to investigate
the possibilities of an online estimation of the drag and thrust coefficients, and the
effects this would have on the wind resistance.
Another approach would be to further investigate the influence of propeller parameters such as blade twist and chord distribution. Optimizing the geometry of a propeller is no easy task. It gets even more complex when structural properties such as
vibrations, which has been ignored during this report despite the evident effects it
has on multirotor performance, are considered.
Moving forward with a multirotor design suitable for aiding in search and rescue
missions, there are several implementation aspects worth considering for further
project work. For instance, it would be interesting to model and evaluate the prosed
Y6 multirotor with the added feature of servo-tiltable rotors. Also of high interest toward these applications would be the integration and control of an on-board
thermal camera, attached to a motorized gimbal joint.
123

Bibliography
[1]

Aerospace; concepts, quantities and symbols for flight dynamics; aircraft motion relative to the air; ISO 1151-1:1988 modified. Deutsches Institut fr
Normung, Berlin, Germany, 1990.

[2]

M. Bangura and R. Mahoney. Nonlinear Dynamic Modeling for High Performance Control of a Quadrotor. In: Australasian Conference on Robotics
and Automation, 3-5 December 2012. Victoria University of Wellington,
New Zealand.

[3]

S. Bouabdallah. Design and Control of Quadrotors with Application to Autonomous Flying. Ecole Polytechnique Federale de Lausanne, Lausanne,
Switzerland, 2007.

[4]

J. Brandt and M. Selig. UIUC Propeller Database.

http://aerospace.illinois.edu/m-selig/props/propDB.html.
Accessed: 2014-04-06.

[5]

J. B. Brandt and M. S. Selig. Propeller Performance Data at Low Reynolds

Numbers. In: 49th AIAA Aerospace Sciences Meeting, 4-7 January 2011.
Orlando, USA.

[6]

T. Bresciani. Modelling, Identification and Control of a Quadrotor Helicopter. Department of Automatic Control, Lund University, Lund, Sweden,
2008.

[7]

C. P. Coleman. NASA Technical Paper 3675: A Survey of Theoretical and

Experimental Coaxial Rotor Aerodynamic Research. National Aeronautics
and Space Administration (NASA), Moffett Field, USA, 1997.

[8]

M. J. Cutler. Design and Control of an Autonomous Variable-Pitch Quadrotor Helicopter. Department of Aeronautics and Astronautics, Massachusetts
Institute of Technology, Boston, USA, 2012.

[9]

S. Driessens and P. E. I. Pounds. Towards A More Efficient Quadrotor

Configuration. In: 2013 IEEE/RSJ International Conference on Intelligent
Robots and Systems (IROS), 3-7 November 2013. Tokyo, Japan.

124

Bibliography
[10]

G. J. J. Ducard and M.-D. Hua. Discussion and Practical Aspects on Control Allocation for a Multi-Rotor Helicopter. In: Conference on Unmanned
Aerial Vehicle in Geomatics (UAV-g 2011), 14-16 September 2011. ETH,
Zrich, Switzerland.

[11]

L. Eriksson, E. Johansson, N. Kettaneh-Wold, C. Wikstrm, and S. Wold.

Design of Experiments: Principles and Applications. 3rd ed. Umetrics
Academy, Ume, Sweden, 2008.

[12]

ESTECO. Modefrontier.
http://www.esteco.com/modefrontier.
Accessed: 2014-06-16.

[13]

J. Fogelberg. Navigation and Autonomous Control of a Hexacopter in Indoor

Environments. Department of Automatic Control, Lund University, Lund,
Sweden, 2013.

[14]

E. Fresk and G. Nikolakopoulos. Full Quaternion Based Attitude Control

for a Quadrotor. In: European Control Conference (ECC), 17-19 July 2006.
Zrich, Switzerland.

[15]

T. Glad and L. Ljung. Reglerteori - Flervariabla och olinjra metoder.

2nd ed. Studentlitteratur AB, Lund, Sweden, 2011.
D. A. Griffiths and J. G. Leishman. A Study of Dual-Rotor Interference and
Ground Effect Using a Free-Vortex Wake Model. In: 58th Annual Forum
and Technology Display of the American Helicopter Society International,
11-13 June 2002. Montral, Canada.

[16]

[17]

N. Kau. APC 10x4.7 Propeller CAD model.

http://whenrobotsfly.tumblr.com/.
Accessed: 2014-04-21.

[18]

J. G. Leishman. Principles of Helicopter Aerodynamics. 2nd ed. Cambridge

University Press, New York City, USA, 2002.

[19]

J. G. Leishman and S. Ananthan. Aerodynamic Optimization of a Coaxial

Prorotor. In: 62th Annual Forum and Technology Display of the American
Helicopter Society International, 9-11 May 2006. Phoenix, USA.

[20]

V. Martnez. Modelling of the Flight Dynamics of a Quadrotor Helicopter.

Department of Aerospace Sciences, Cranfield University, Cranfield, United
Kingdom, 2007.

[21]

P. Marwedel. Pareto Points - Lecture slides.

http://ls12-www.cs.tu-dortmund.de/daes/media/documents/teaching/courses/
ws1213/es/lecture/es-marw-5.1-evaluation.pdf.
Accessed: 2014-06-05.
Mathworks. Aerospace Wind Model.
http://www.mathworks.se/help/aeroblks/wind.html.
Accessed: 2014-05-15.

[22]

125

Bibliography
[23]

Mathworks. MATLAB Simulink.

http://www.mathworks.se/products/simulink.
Accessed: 2014-06-16.

[24]

Advanced precision composites. Model Airplane Propellers.

http://www.apcprop.com/v/index.html.
Accessed: 2014-05-23.

[25]

CJ Youngblood Enterprises Inc. Stingray 500.

http://curtisyoungbloodhobbies.curtisyoungblood.com/Stingray-500Airframe-ND-YS5-K4001-9712.htm.
Accessed: 2014-06-20.

[26]

Cyber Technology. CyberQuad.

http://www.cybertechuav.com.au/.
Accessed: 2014-06-20.

[27]

Dassault Systmes. Solidworks.

http://www.solidworks.com/.
Accessed: 2014-06-16.

[28]

Svenska Dagbladet. Polisen vill anvnda drnare.

http://www.svd.se/nyheter/inrikes/polisen-vill-anvanda-dronare_3571454.svd.
Accessed: 2014-06-25.

[29]

Sydsvenska Dagbladet. Polisen vill ha spionplan.

http://www.sydsvenskan.se/skane/polisen-vill-ha-spionplan/.
Accessed: 2014-06-25.

[30]

Modde. Modde - Design of Experiments.

http://www.umetrics.com/products/modde.
Accessed: 2014-06-16.

[31]

J. R. Movellan. DC Motors. Machine Perception Laboratory, University of

California, San Diego, USA, 2010.

[32]

NIST/SEMATECH. E-Handbook of Statistical Methods.

http://www.itl.nist.gov/div898/handbook.
Accessed: 2014-06-02.

[33]

N. Ohlsson and M. Sthl. Model-Based Approach to Computer Vision and

Automatic Control using Matlab Simulink for an Autonomous Indoor Multirotor System. Department of Signals and Systems, Chalmers University of
Technology, Gothenburg, Sweden, 2013.

[34]

P. Pounds, R. Mahoney, and P. Corke. Modelling and Control of a QuadRotor Robot. In: Australasian Conference on Robotics and Automation, 6-8
December 2006. Auckland, New Zealand.

[35]

P. E. I. Pounds. Design, Construction and Control of a Large Quadrotor Micro Air Vehicle. Australian National University, Canberra, Australia, 2007.

126

Bibliography
[36]

J. Praceus. Modellering av vind och turbulens fr en High-Fidelity flygsimuleringsmodell. Department of Automatic Control, Lund University, Lund,
Sweden, 2013.

[37]

R. W. Prouty. Helicopter Performance, Stability and Control. Reprint edition

with corrections. Krieger, Malabar, USA, 1990.

[38]

M. Rutkowski and W. Kruszi. Design and analysis of ducted fan for multirotor VTOL UAV. Polytechnika Gdansk, 2013.

[39]

M. Ryll, H. H. Blthoff, and P. R. Giordano. First Flight Tests for a Quadrotor UAV with Tilting Propellers. Max Planck Institute for Biological Cybernetics, Tbingen, Germany, 2013.

[40]

K.-E. rzn. Real-Time Control Systems. 2012 edition. Department of Automatic Control, Lund University, Lund, Sweden, 2012.

[41]

R. F. Stengel. Aircraft Flight Dynamics MAE331 - Lecture slides.

http://www.princeton.edu/stengel/MAE331Lectures.html.
Accessed: 2014-02-25.

[42]

The Engineering Toolbox. Drag Coefficient.

http://www.engineeringtoolbox.com/drag-coefficient-d_627.html.
Accessed: 2014-06-02.

[43]

Transportstyrelsens freskrifter om obemannade luftfartyg (UAS) TSFS

2009:88. Transportstyrelsen, Norrkping, Sweden, 2009.

[44]

K. Ulrich and S. Eppinger. Product Design and Development. 5th ed.

McGraw-Hill/Irwin, New York City, USA, 2012.

[45]

U.S. Standard Atmosphere 1976. National Aeronautics and Space Administration (NASA), Washington D.C., USA, 1976.

[46]

Q.-G. Wang, Z. Zhang, K. J. strm, Y. Zhang, and Y. Zhang. Guaranteed

Dominant Pole Placement with PID Controllers. In: 17th World Congress
of the International Federation of Automatic Control, 6-11 July 2008. Seoul,
South Korea.

127

A
6-DoF Kinematics
A rigid body mechanical systems, such as the flying multirotor, is typically described by motion with six degrees of freedom. These include translational motion
along each coordinate axis as well as rotation around each of these axes. To simplify
the system description, a body fixed coordinate system is often used where the coordinate axes follow the orientation of the body. Usage of such a coordinate frame
is convenient in flight applications due to that applied forces and moments as well
as sensor readings often are defined in relation to the vehicle orientation. Also, with
the body frame origin placed in the center of gravity, the moment of inertia does not
change with the orientation of the body.
A convenient way of transforming between body and inertial system descriptions
is the application of rotation transformation matrices. In this section, these will be
derived according to flight dynamics theory according to Robert F. Stengel [41]. In
the simulator, the MATLAB toolbox SimMechanics has been employed to handle
these transformations, but will nonetheless be described here.
Rotation angles around a set of axes are known as Euler angles. When studying
flight dynamics, the angles , and are often referred to as roll, pitch and yaw
angles, which represent the angle rotated around the x, y and z axes respectively.
These angles are used to transform between the inertial (I) and body frame (B) descriptions through a series of transformations, visiting two intermediate coordinate
systems, (1) and (2). Step by step transformation for each rotation angle is described
by the matrices in (A.1).

cos
H1I () = sin
0

sin
cos
0

0
0
1

cos
H21 ( ) = 0
sin

0
1
0

sin
0
cos
(A.1)

1
HB2 ( ) = 0
0
128

0
cos
sin

0
sin
cos

Appendix A. 6-DoF Kinematics

Combining each rotation, (A.1) yields the direction cosine matrix (DCM) (A.2),
which transforms translational velocities and accelerations from the inertial frame
to the body frame. The inverse of the DCM describes the transformation in the
opposite direction. Because of the orthogonality of the DCM, (HBI )1 = (HBI )T .
HBI ( , , ) = H1I H21 HB2 =

cos cos
cos sin
sin

sin sin cos sin cos

sin sin sin + cos cos
sin cos

(A.2)

cos sin cos + sin sin

cos sin sin sin cos
cos cos

To obtain a similar transformation for the rotational velocities and accelerations, the
matrix formula in (A.3) is used, where p, q and r are the rotational velocities around
each of the body frame axes.




0
p
0

q = I3x3 0 + HB2 + HB2 H21 0

r

0
0

(A.3)

Rewriting (A.3), the transformation matrix LIB (A.4) is obtained to transform rotational motion from the body to the inertial frame. Again, the inverse can be used to
reverse the transformation, going from the inertial frame to the body frame coordinate systems. However, note that LIB is not orthogonal.

1
LIB ( , , ) = 0
0

sin tan
cos
sin / cos

cos tan
sin
cos / cos

(A.4)

As seen in (A.4), the matrix LIB becomes singular for angles = /2. At these
pitch angles, the transformation becomes undefined. Making use of a quaternion angle representation provides a remedy for this. However, further research into quaternions were not deemed necessary by the authors since such aggressive flight maneuvers are neither desired nor anticipated. Modeling and control with quaternions has,
however, been performed by Fresk [14] with promising results.

129

B
Geometric and inertia
calculations
Calculation of geometric dimensions and mass properties of the modeled multirotor
UAV is the main topic of this appendix chapter. This is necessary to evaluate the
relative influence of quantities such as total mass, moment of inertia and projected
surface area of different multirotor setups in the simulator.
The main components of the simulated multirotor, summarized in Table B.1, are
well described by two fundamental shapes for the airframe along with the propeller.
These are the cylinder and the rectangular cuboid, for which standardized formulas
exist to obtain the inertia and geometric properties. Using a rectangular cuboid and
a cylinder, the arm geometry can be approximated. A pure cylinder shape represents
the center, motors and propeller shafts. In addition to this, propellers modeled with
CAD have been imported along with their geometric and inertia properties.
Cylinder

Center, motor
shaft shapes.

and

Arm

Propeller

Approximated to a rectangular cuboid and a cylinder

in calculations.

Propeller properties
are imported from a
CAD model.

Table B.1: List of components.

130

B.1

B.1

Mass

Mass

The combined mass of each part of the UAV add up to the total mass. The total mass
mtot of a multirotor configuration with na arms and nr rotors is calculated by using
Equation (B.1) below.

mtot = mc + na ma + nr mm + nr ms + nr m p

(B.1)

The symmetrical design of the multirotor aircrafts in this thesis means that the location of the center of gravity G will only vary vertically in the body frame. For
arms attached to the middle of the cylindrical center part, the vertical location of the
center of gravity is given by expression (B.2). Note that the distances to each part is
measured from the middle of the center part to the parts local center of gravity.

G=

B.2



 

zs
1
za
za zm
nr mm
+
+ ms
+ zm +
mtot
2
2
2
2


za
+m p
+ zm + zs
2

(B.2)

Moment of inertia

A bodys moment of inertia describes the relation between angular acceleration and
applied moment around a given axis. This relationship is given by (B.3), where
is the applied moment, the angular acceleration and I the moment of inertia.
Hence, moment of inertia plays the same role for rotational motion as mass does for
translational.

= I

(B.3)

However, calculating the moment of inertia is generally more cumbersome than

mass. For an arbitrary body rotating around a given axis, the moment of inertia
I is given by (B.4). In this equation, is the bodys density, r a vector from the
rotational axis to a point in the body and V the bodys volume. Note that the used
here is not the air density.
Z

(r)r2 dV

(B.4)

131

Appendix B. Geometric and inertia calculations

For a three-dimensional body defined in a coordinate system with x-, y- and z-axes,
the moment of inertia takes the form of a matrix called the Inertia tensor defined as
(B.5). The tensor fully describes the rotational properties of a body rotating around
any axis defined in the coordinate system. It is a symmetric matrix, which means
that opposing terms mirrored in the diagonal are equal to each other (for instance
Ixy = Iyx ). Using the inertia tensor, the relationship between a three-dimensional

applied moment and angular acceleration is given by = I.

Ixx
I = Iyx
Izx

Ixy
Iyy
Izy

Ixz
Iyz
Izz

(B.5)

Since the multirotors studied in this thesis are of symmetric design, their arms are
only placed around the vertical axis in the body frame. To align a locally defined
inertia tensor, rotated an angle with respect to the global body frame axes, the
transformation in (B.6) is performed. In the equation, IL is the aligned local inertia
tensor, whereas ILR is the rotated local inertia tensor.

IL = TILR TT

cos
T = sin
0

sin
cos
0

0
0
1

(B.6)

The local moment of inertia tensors may be calculated for each component using the
formulas in Table B.2. To calculate the total moment of inertia tensor for the global
body frame, the local moment of inertia tensors for each part must first be rotated
accordingly so that their axes are aligned with those of the global body frame. Then,
by making use of the parallel axis theorem (B.7), the total moment of inertia tensor
for the entire body can be obtained. In the theorem, d is the distance between the
rotation axes and m the component mass.
IO = IL + md 2

B.3

(B.7)

Geometry

The projected surface area is the area which is seen from the perspective of the
incoming airflow. This can easily be calculated for its body frame x, y and z directions using knowledge of its setup and component dimensions. Proceeding with the
quadrotor case, the following expressions in (B.8) gives the projected surface areas
132

B.3
Cylinder

Cuboid

1
m(3r2 + h2 )
12
1
= m(3r2 + h2 )
12
1
= mr2
2

1 2
(h + w2 )
12
1
= (h2 + l 2 )
12
1
= (w2 + l 2 )
12

IXX =

IXX =

IYY

IYY

IZZ

Geometry

IZZ

Table B.2: Shapes local moment of inertia.

in the three main directions. The exposed surface area will shift depending on the
direction of the incoming airflow. A discussion of this is given in Section 2.4.
Ax = Ay = 2rc hc + 2(la rc )ha + 4rm hm + 4rs hs
2
Az = rc2 + 4(la rc )wa + 4rm

(B.8)

Airframe drag force will, as discussed in Section 2.4, act in the airframe centroid.
The centroid can be viewed as a center of geometry with area spread out evenly
around it. For an arbitrary polygon body surface, it is calculated using its corner
points. For a polygon with k corner points xi , i = 1, .., k, the centroid C is calculated
according to (B.9).

C=

1 k
xi
k
i

(B.9)

133

C
Aerodynamic theory
In this appendix chapter, derivation of aerodynamic rotor theory is given. Since the
model discussed in Section 2.3 makes use of both momentum theory and blade element theory, they are both described at length in the following sections. Momentum
theory provides a basis for rotor performance evaluation and is used in the deduce
the rotor inflow velocity in the proposed model. Blade element theory gives a finer
description of rotor forces and moment for each blade section. Applying integration
along each blade and for a full rotor revolution, a much more accurate rotor model
is obtained.

C.1

Momentum theory

Momentum theory provides a coarse prediction of rotor flight performance. It treats

the flow of air through a rotor disk as a tube of incompressible fluid, without any
external interaction or variations across each rotor segment. The airflow itself is
considered to be one-dimensional, i.e. its properties only vary perpendicular to the
rotor plane. The rotor disk is assumed infinitesimally thin with an infinite amount of
propeller blades, over which a pressure difference exists. The rotor work, introduced
by an applied moment around the rotational axis, leads to added kinetic energy
in the rotor slipstream. Due to the resulting airflow, an upward thrust is achieved.
The momentum theory derivations presented in this section are made according to
Leishman [18] and leads to equations for thrust force, applied moment and power
in hover, axial and forward flight.

Hover flight
The momentum theory model for hovering flight is described by Figure C.1. In the
figure, the numbers 1 and 2 indicate cross section planes just above and below the
rotor disk with area A. Between these planes, a pressure difference exists. The speed
vi is the airflow velocity induced by the rotor. Sufficiently far upstream, the cross
section plane indicated by a 0 is regarded as a boundary where the rotor no longer
has an effect on the air velocity. In hover, this will mean that the air velocity is zero,
V0 = 0, in this region. The sign indicates the far wake plane, also known as the
134

C.1

Momentum theory

vena contracta. This is where it is assumed that the flow is fully developed with
velocity V = w. Note that since the rotor tube contracts below the rotor and that
mass is conserved, the air velocity is increased below the rotor.

Figure C.1: Momentum theory for hovering flight [18].

Applying the principle of conservation of mass, the mass flow rate m must be constant within the rotor wake. This yields Equation (C.1).
m = A w = Avi

(C.1)

By using the principle of conservation of fluid momentum between the 0 and

planes, the resulting rotor thrust force T can be expressed as (C.2).
T = mw

(C.2)

Proceeding by applying conservation of energy between the 0 and planes, the

power P, i.e. work done per unit of time, is given by (C.3).
1
2
P = T vi = mw
2

(C.3)

Combining the expressions (C.2) and (C.3), the relation w = 2vi can be equated.
Using this along with the mass flow (C.1), the expression for the rotor thrust and
power can be re-written as T = 2Av2i and P = 2Av3i . Also, the contraction of the
135

Appendix C. Aerodynamic theory

rotor slipstream can be described by the relation A = 21 A between the rotor and the
vena contracta areas.
The induced inflow velocity vi , also denoted as vh for hover, can be written as (C.4),
where R is the rotor radius, the rotor rotational speed and h is called the induced
inflow ratio in hover. Again, using (C.4), the thrust, power and the applied rotor
axis moment Q can be re-written as functions of h and the rotor rotational speed
. These expressions are shown in (C.5).

s
vh vi = h R =

T = 2Ah2 R2 2

T
2A

Q = 2Ah3 R3 2

(C.4)

P = 2Ah3 R3 3

(C.5)

Introducing nondimensional coefficients to describe rotor thrust, applied moment

and power is common practice in helicopter engineering and are defined according
to (C.6). Note that since the applied moment Q and power are related as P = Q,
their corresponding coefficients will have the same numerical values.
Using the
p
thrust coefficient CT , the induced inflow ratio is described by h = CT /2 in hover.

CT =

T
AR2 2

CQ =

Q
AR3 2

CP =

P
AR3 3

(C.6)

Axial flight
When modeling axial flight with MT, there are two subcases to consider ascent
and descent where the net air flow has opposite directions. Starting with the
ascent case, Figure C.2 describes the MT model. The situation is similar to that of
hover, but with an added downward airflow speed resulting from the climb rate. As
can be seen in the figure, the far upstream plane now has a flow speed relative to
that of the ascending aircraft Vc , so V0 = Vc . The flow speed at the rotor is vi + Vc
and at the vena contracta the slipstream speed is V = Vc + w
Using the same procedure as for hover flight, the mass flow and resulting thrust
and power can be calculated as (C.7). Note that the added terms for Vc cancel each
other out for the thrust T at this stage, yielding the same thrust equation as for hover
flight.
136

C.1

Momentum theory

Figure C.2: Momentum theory for axial ascent flight [18].

m = A (Vc + w) = A(Vc + vi )
T = m(V
c + w) mV
c = mw

1
1
1
P = T (Vc + vi ) = m(V
c + w)2 mV
c2 = mw(2V

c + w)
2
2
2

(C.7)

It is clear from (C.7) that w can again be calculated as w = 2vi . This leads to that
the thrust can be re-written as (C.8) for ascent flight.

T = A(Vc + vi )w = 2A(Vc + vi )vi

(C.8)

Recalling the equation for vh (C.4), inserting the ascent flight thrust equation (C.8)
into (C.8) yields the quadratic expression v2h = (Vc + vi )vi for the induced inflow
velocity vi . Hence, with prior knowledge of h , MT can predict i and vi for ascent
flight, making it possible to evaluate thrust, applied moment and power. Note that
the positive solution of the quadratic function must be chosen.
For the descent axial flight case, the air is assumed to have a net flow in the upward
direction. In Figure C.3, this is reflected by that the locations of the planes 0, , 1
and 2 have been reversed.
137

Appendix C. Aerodynamic theory

Figure C.3: Momentum theory for axial descent flight [18].

The mass flow m is of the same magnitude as for ascent flight, but with a negative
sign. In descent flight, the thrust T = mw
will still be directed upwards, whereas
the power P will be negative. This is the case since the rotor is now extracting
power from the air, braking its speed as it passes the rotor disc. For this reason, the
described flight condition is known as the windmill brake state.
The same method used to calculate the induced inflow ratio i as for ascent flight
can be applied for descent. The quadratic function v2h = (Vc + vi )vi is used, where
the correct solution is the one that does not violate vi /vh > 1. The descent flight
assumption that the net flow is directed upwards limits the MT model to a great
extent. In the range 2vh Vc 0, the flow can be either upward or downward,
making the MT model invalid in this range of axial speed Vc .

Forward flight
To achieve forward flight velocity, the rotor is tilted forward creating an angle of
attack to the airflow. The flow model is described by Figure C.4, where it can be
seen that the airflow has components both perpendicular and parallel to the rotor
disk.
The air mass flow through the rotor is given by (C.9), where U is the resultant air
flow speed at the rotor disk.
138

C.1

Momentum theory

Figure C.4: Momentum theory for forward flight [18].

m = AU

U=

(V cos )2 + (V sin + vi )2

(C.9)

Applying the principle of conservation of momentum perpendicular to the rotor

disk, the thrust T can be calculated as (C.10) for forward flight. Again, the thrust is
described on the same form T = mw,
although with a different expression for m.

T = m(w
+V sin ) mV
sin = mw

(C.10)

By yet again applying conservation of energy, the power can be expressed as C.11.
As before, the power and thrust equations are used to obtain the relation w = 2vi ,
which is the same for hover and axial flight.
1
2
P = T (vi +V sin ) = m(2V

w sin + w )
2

(C.11)

Proceeding with these findings, the thrust can finally be written as (C.12).

T = 2Avi

q
(V cos )2 + (V sin + vi )2

(C.12)
139

Appendix C. Aerodynamic theory

As in the axial flight case, the expression for the hover induced velocity vh (C.4) is
used to determine vi for forward flight. Inserting the forward flight thrust equation
(C.12) into (C.4), vi can be expressed as (C.13).

vi = p

v2h
(V cos )2 + (V sin + vi )2

(C.13)

The advance ratio is definied as the velocity parallel to the rotor disk, =
V cos /(R). Using this, the total inflow ratio can be expressed as (C.14),
which includes the inflow induced by both the rotor and from the aircrafts relative airspeed.

V sin + vi
= tan + i
R

(C.14)

p
Using (C.13) and recalling h = CT /2, the induced inflow ratio i and the inflow
ratio can be expressed as in (C.15).

i = p

h2
2 + 2

CT
= p
2 2 + 2

CT
= tan + p
2 2 + 2

(C.15)

The equations for i and in (C.15) are applicable for all the flight modes where
momentum theory is valid. They allow calculation of the MT predictions of thrust
force, power and applied moment.

C.2

Blade element theory

Blade element theory (BET) is a widely used method for calculating the aerodynamic forces and moments acting on a blade. BET assumes that each section of the
blade is a quasi-2-D airfoil that produces aerodynamic forces and moments [18],
which can be seen in Figure C.5. These infinitesimal forces and moments are then
integrated over the length of the blade and averaged over one rotor revolution.
The infinitesimal lift force dL and the infinitesimal drag force dD is perpendicular
and parallel to the local velocity resultant U. The expressions for these forces can be
seen in Equation (C.16), where q = 12 U 2 is the dynamic pressure. In the equation,
Cl and Cd are the local lift and drag coefficients while cdr is the reference area
where c is the chord length.
140

C.2

Blade element theory

Figure C.5: Aerodynamic forces and moments acting at an infinitesimal blade element [18].

dL = qCl cdr
dD = qCd cdr

(C.16)

The determination of the moment dM around the aerodynamic center is neglected

in this report because of the complexity, and since it was decided not to model the
feathering motion.
A number of simplifications, listed below, are introduced to make the derivation of
the aerodynamic forces and moments significantly easier.
The local lift coefficient can be written as Cl = Cl ( ) [18] where Cl
can be approximated to 2 according to the thin-airfoil theory.
The local drag coefficient Cd is assumed to be constant and has been approximated in Section 2.3.
The chord length is assumed to be constant.
The pitch angle of the blade section is assumed to vary linearly along the
length of the blade = 0 + tw r/R
The velocity component that is parallel to the reference plane UT is assumed
to be much bigger than the perpendicular component UP . This leads to the
following two simplifications:
q
1. The resulting velocity U = UT2 +UP2 UT
 
UP
UP
2. The induced angle = arctan U
U
T
T
141

Appendix C. Aerodynamic theory

By applying these simplifications, the following expressions (C.17) for the lift and
drag forces are obtained.
1
r
dL = cCl ((0 + 1 )UT2 UPUT )dr
2
R
1
dD = UT2 cCd dr
2

(C.17)

It is convenient to express the forces in the directions that are parallel and perpendicular to the hub plane. By once again using the simplification that the induced
angle is small and that the drag force is at least one magnitude smaller than the
lift force, the following expressions (C.18) are obtained.
dFz = dL cos dD sin dL
UP
dL + dD
dFx = dL sin dD cos
UT

(C.18)

Multiplying the expressions in C.18 with the number of blades Nb yields the infinitesimal thrust, horizontal force and drag moment of the rotor (C.19).

dT = Nb dFz
dH = Nb dFx

(C.19)

dQ = Nb ydFx
In blade element theory, it is standard practice to introduce nondimensional coefficients according to (C.20), where the velocity components UT and UP are expressed
in (C.21) [20].

Nb 12 cCl ((0 + 1 Rr )UT2 UPUT )

dT
=
dr
R4 2
R4 2


UP 1
r
2 U U ) + 1 U 2 cC
N
cC
((
+

)U
P
T
0
1
b
l
d
T
T
UT 2
R
2
dH
dCH =
=
dr
R4 2
R4 2


UP 1
r
2 U U ) + 1 U 2 cC
N
y
cC
((
+

)U
P
T
0
1
b
l
d
T
T
U
2
R
2
dQ
T
dCQ =
=
dr
R5 2
R5 2
(C.20)
dCT =

142

C.2

Blade element theory

r+e
+ sin
R

 r
r
UP = R c i 1 + Kc cos (a1s sin b1s cos
R
R
!

UT = R

(C.21)

(a0 a1s cos b1 s sin )cos + (r + e)(qw cos + pw sin )

The expressions in (C.21) include several variables which are listed below.

= c + i =
=

V sin + vi
is the inflow ratio.
R

V cos
is the advance ratio.
R

is the angle of attack of the rotor disk. Note that this is not the same angle
of attack as the ones that each individual blade experiences.
V is the total magnitude of the free flow stream.
e is the hinge offset.
a0 , a1s and b1s are angles that describe the blade flapping motion. The expressions for these angles are derived in Appendix C.3.
The aerodynamic coefficients are determined by integrating the expressions in
Equation (C.20) over the length of the blade, i.e. from the root r0 to R. This integral is a function of the azimuth angle and the final result is calculated as the
average value over a revolution. It is obtained by integrating over the azimuth angle
from 0 to 2 and dividing by the length of interval, resulting in Equation (C.22).

1
CT =
2

Z2ZR
0 r0

1
dCT d CH =
2

Z2ZR
0 r0

1
dCH d CQ =
2

Z2ZR

dCQ d (C.22)
0 r0

Note that Prouty [37] assumes that the blade root r0 coincides with the hinge offset
e. The same assumption has therefore been made in this report to avoid reformulation of the expressions. It has also been used that e << R to simplify the expressions
143

Appendix C. Aerodynamic theory

somewhat [20]. It is assumed that the rotational speed of the rotor is constant over
the integrated revolution. This is true for conventional helicopters, but generally not
for multirotors. However, testing has shown that the rotational speed change is typically less than 2 % during one revolution and the assumption is therefore considered
to be reasonable. The final expressions for CT , CH and CQ can be seen in (C.23),
where is the rotor solidity which is the total blade area divided by the disk area
= Nb cR/R2 = Nb c/R.



e
2
Cl 
1
1
+ 2 0 + (1 + 2 )tw
CT =
4
R
3
2




Cla Cd
+ a1s
1 3
1
3
CH =

+ 2 0 +
1 + 2 tw +
3
4
Cl
3 2
2
2
1 + 2
2
 2


!
a0 1 2

1
1 2
a1s +
+
+ a0 i + i
1
2 9
2
3
8
1 + 2
2



Cla 1 Cd
+ a1s 1
1
2
CQ =
(1 + 2 )
(2 2 )0 +
1
tw +
3
4
2 Cl
3
2
2
1 + 2
2


 2
!
a0 1 2
2
1
1 2
+
+ a0 i + i
( + a1s ) ( + a1s )
1
2 9
2
3
8
1 + 2
2
(C.23)

C.3

Blade flapping

In forward flight, the rotor disk experiences an uneven lift distribution as the advancing and retreating blades experience different relative airspeeds. The advancing blade will generate additional lift whereas the retreating blade will generate less,
resulting in a slight tilt of the rotor plane. This is known as blade flapping motion,
where a blades tilted angle will vary depending on its location around the propeller shaft. Even when not in forward flight, the blades undergoes a constant tilting
known as coning. This effect is due to the lift force, which in axial flight is constant
over a revolution, and hence doesnt tilt the rotor plane as a whole. The theory and
mathematical modeling of this phenomena is well-described in the works of Leishman [18], Pounds [35] and Martnez [20], which have all inspired the blade flapping
modeling performed in this thesis work.
As discussed in Section 2.3, a mechanical rotor model with an effective flapping
hinge offset e and a torsional spring strength k is assumed. It is around this hinge
offset that the balancing of forces and moments is calculated to deduce the flapping
angle. The flapping angle of a blade rotating around the shaft is calculated using
a partial differential equation. The first-order solution to this equation, which is
144

C.3

Blade flapping

sufficient to model the flapping behavior accurately, can be expressed as a function

of the azimuthal angle , i.e. the blade position relative to the airflow, as (C.24).
The variables a1s and b1s are the longitudinal and lateral flapping angles whereas
a0 is the coning angle. These will be the angles of main interest when deducing the
rotor tip path plane.
= a0 a1s cos b1s sin

(C.24)

The three flapping angles, a0 , a1s and b1s , are determined by balancing the lift force,
centripetal force, gyroscopic moments and the weight of the propeller blades around
the hinge offset e. Performing the balancing around the flapping hinge, a0 , a1s and
b1s may be calculated from (C.25). M, Mpq , F and Fpq are matrices containing
expression for the various forces moments applied around the hinge listed in (C.27).
The added gyroscopic terms Mpq and Fpq are added to account for the bodys pitch
and roll rotation.

a0
(M + Mpq ) a1s = F + Fpq
(C.25)
b1s
The moments and forces matrices in (C.27) contain several variables that will be
briefly explained here.
e is the effective blade hinge offset from the rotor shaft.
k is spring constant of the blade.
Mb is the moment due to weight force of the propeller.
Iyb blade moment of inertia around the flapping hinge.
=

Cl cR4
is known as the Lock number.
Iyb

pw and qw rotational velocities of the rotor around the rotor wind axes.
The effective hinge offset e was deduced in similar fashion to [20] and is shown in
(C.26). The parameter r0 is the blade offset from the shaft joint, which was measured
for the modeled propeller.

e=R

1 1
Mb R
gIyb

1 = 1.19 + 1.57

r0
R

(C.26)

145

146



Mb
2 Iyb + e
k
g




e 2
2
1
Iyb 2 1
1
8
R
2
Mb
2 e
+ k
g


2

Mb
2 e
+ k
g



1
e 2
Iyb 2 1
1
8
R
2


1
e 2 e
Iyb 2 1
2
R
4R



1
e 2 1 2 1 e
1 e2
2
Iyb 1

2
2
R
4
8
6 R 12 R
Mb
2
e
+ k
g

0

1

e 2 pw
2

+ 2qw Iyb
Fpq = 8 Iyb 1 R


e 2 qw
2
Iyb 1
+ 2pw Iyb
8
R

+ k





2
2
2 1
1
e

1
e
1
e
Iyb 2 1

2
2
R
4
8
6 R 12 R
Mb
2 e

 







e 2 1
2 e 1  e 2
1 2 
e 1 e
1  e 2
1  e 3
1 1e
1
2
1 + 2 +
+
0 +
+
1

tw
+
2 Iyb 1 R

4
3R 3 R
5
6
R
10 R 15 R
30 R
3 6R

 









2
2
1
e
1 1e
1 1e
1 e

F=

2
Iyb 1
2
+
0 + 2

tw

2
R
3 6R
4 6 R 12 R
2

Mpq = 0

0
M=

2  1 1 e 

1
e
Iyb 2 1

+
2
R
3 6R

(C.27)

Appendix C. Aerodynamic theory

C.4

C.4

Coaxial rotors

Coaxial rotors

Compared to a conventional single propeller rotor, there are some additional difficulties which has to be considered when a coaxial rotor is modeled. The lower rotor
is affected by the induced inflow velocity generated by the upper rotor. This downwash will have a quite severe impact on the overall performance of the multirotor
and it is therefore deemed necessary to include this phenomenon when modeling a
coaxial rotor. The lower rotor will also affect the upper rotor, but this effect doesnt
influence the overall performance as much [7] and will therefore be neglected in
this report. The distance between the two propellers will not be considered either,
since this distance only affects the distribution of the thrust between the two rotors
according to Coleman [7].
The MT model for a coaxial rotor can be seen in Figure C.6. The slipstream of
the upper propeller flows down through the inner part of the lower propeller. The
slipstream typically contracts quickly (within 0.25R) according to Leishman [18]
and Ananthan [19] and it is therefore assumed that the slipstream of the upper rotor
is fully developed when it reaches the lower.

Figure C.6: Momentum theory model for a coaxial rotor in hover [18].
147

Appendix C. Aerodynamic theory

The contracted slipstream results in an uneven distribution of the inflow velocity
across the lower rotor disk. To still be able to use the expressions in Equation (C.23),
it was decided to use the ratio between the rotor disk area and the vena contracta
area A = 12 A, which was described in C.1. By doing so it was possible to calculate
an average inflow velocity according to Equation (C.28), where vu and vl are the
induced velocities for the upper and lower rotor.

m =

A
A
(2vu + vl +Vc ) + (vl +Vc ) = A (vu + vl +Vc )
2
2

(C.28)

The induced velocity of the upper rotor is therefore added to the inflow of the lower
rotor while calculating the aerodynamic forces and moments. It should be noted that
this is a simplification which had to be done because of the time constraints. The
simplification will have some affect on the overall results, but the amount is hard
to estimate. It was assumed that the error introduced by this simplification isnt any
larger than those introduced by not modeling for example tip losses and feathering.
It was also decided that an optimization of the performance of the coaxial rotor,
where twist and chord length along the upper and lower blades are determined,
would not be conducted.

148

Lund University
Department of Automatic Control
Box 118
SE-221 00 Lund Sweden

Document name

MASTERS THESIS
Date of issue

June 2014
Document Number

ISRN LUTFD2/TFRT--5947--SE
Author(s)

Supervisor

Christian Mnsson
Daniel Stenberg

Simon Yngve, Combine AB

Rolf Johansson, Dept. of Automatic Control, Lund
University, Sweden
Anders Rantzer, Dept. of Automatic Control, Lund
University, Sweden (examiner)
Sponsoring organization

Title and subtitle

Model-based Design Development and Control of a Wind Resistant Multirotor UAV

Abstract

Multirotor UAVs have in recent years become a trend among academics, engineers and hobbyists
alike due to their mechanical simplicity and availability. Commercial uses range from surveillance to
recreational flight with plenty of research being conducted in regards to design and control.
With applications towards search and rescue missions in mind, the main objective of this thesis work
is the development of a mechanical design and control algorithm aimed at maximizing wind
resistance. To these ends, an advanced multirotor simulator, based on helicopter theory, has been
developed to give an accurate description of the flight dynamics. Controllers are then designed and
tuned to stabilize the attitude and position of the UAV followed by a discussion regarding disturbance
attenuation.
In order to study the impact of different design setups, the UAV model is constructed so that physical
properties can be scaled. Parameter influence is then investigated for a specified wind test using a
Design of Experiments methodology. These results are combined with a concept generation process
and evaluated with a control engineering approach. It was concluded that the proposed final design
should incorporate a compact three-armed airframe with six rotors configured coaxially.

Keywords

Classification system and/or index terms (if any)

Supplementary bibliographical information

ISSN and key title

ISBN

0280-5316
Language

Number of pages

English

1-148

Security classification

http://www.control.lth.se/publications/

Recipients notes