ABSTRACT
Title of dissertation: DISTURBANCE REJECTION
FOR U.A.S. AIRCRAFT USING
BIO-INSPIRED STRAIN SENSING
Lina Castano, Doctor of Philosophy, 2015
Dissertation directed by: Professor Sean J. Humbert
Department of Aerospace Engineering
A bio inspired gust rejection mechanism based on structural inputs is pro-
posed. Insect wings possess a wealth of sensor systems which typically consist of
fast reflexive neuronal paths. Stretch and strain sensors on insect wings are used for
flight control and can be found across many species. These are used for monitoring
of bending and torsion during flight. The fast reflexive and proprioceptive mecha-
nisms based on strain sensing found in nature are the inspiration for this work. A
strain feedback controller allows for anticipation of the onset of rigid body dynam-
ics due to gust perturbations. This anticipation stems from sensing of higher order
states and the possibility of reacting before lower order states are reached. High
bandwidth inner loop compensation is therefore enabled. Forces and moments are
proportional to wing strain patterns and can be used in fast reaction inner loops.
Strain sensors are used for providing an indirect estimation of the differential forces
applied to the aircraft wing and therefore to the aircraft rigid body. These sensors
can be distributed over the surface of the aircraft wing to encode multiple degree
of freedom disturbances. Sensor locations for disturbance rejection are determined
based on metrics associated to the observability Grammian. The locations are pre-
selected based on modal energy analyses and are chosen according to wide field
integration patterns. A model for wide field integrated strain based on mass partic-
ipation factors is proposed as well as one which is based on the physics of the forces
and moments acting on the wing producing strain patterns which can be used for
disturbance rejection. Models of the differential forces via strains on the wings are
proposed. Strain feedback was implemented in four platforms under different types
of disturbances. The platforms consisted of a glider, a quadrotor, a wing section for
wind tunnel testing and an RC airplane with a full span wing. The disturbances
included discrete gusts as well as turbulence. The results of using strain feedback
showed not only to be faster than IMU estimations but also to be better when
compared to a classical attitude controller implementation.
DISTURBANCE REJECTION FOR U.A.S. AIRCRAFT
USING BIO-INSPIRED STRAIN SENSING
by
Lina Castano
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2015
Advisory Committee:
Professor Sean J. Humbert, Chair/Advisor
Professor Inderjit Chopra
Professor Alison Flatau
Professor James Baeder
Professor Amr Baz, Dean’s Representative
c© Copyright by
Lina Maria Castano Salcedo
2015
Acknowledgments
I would like to acknowledge a number of people without whom this degree
would not have been possible, or at least extremely unfeasible. Thank you God for
allowing me to be part of this wonderful research institution and to learn every day
from its wonderful professors, peers and friends. Physics for undergrad gave me the
tools for understanding engineering problems, however engineering is a challenging
field on its own and I am grateful to all the people who contributed to this work
one way or the other in many different areas.
My gratitude goes out to my parents Flor and Fernando, and uncle Alvaro
for their love and support during these years of being away and for encouraging me
every single day, no matter what. Being an only child, their constant encouragement
allowed me to have resilience among many life and academic circumstances. Thanks
as well to friends Ps. Mario and Leonor, for their support and encouragement as
well as aunt Cielo, who was also there during these years. I would like to thank a
couple of special friends, Rosalia and Zohaib. Rosalia, thank you for the opportunity
of joining Maryland, I love this place! and you are an amazing person. Thanks for
your friendship and mentor-ship along the way. Zohaib thank you for your help
along these years and for being a solid friend, thank you for your friendship. You
are an amazing academic and person.
To my sponsors, Aurora Flight Sciences for supporting my degree, and Tony
Thompson, from the USAF, whose vision gave birth to the inspiration of this won-
derful project. Thanks also to Riley, program manager for the PWings project for
ii
being great to work with. Thanks also to Simone and Terrence for the work with
the glider while visiting UMD.
To Dr. Humbert, my advisor, thank you for giving me the opportunity to
grow and become productive in the Autonomous Vehicle Lab (AVL). I appreciate
the exciting project I worked on for this degree and which I will continue to develop.
Thank you for sharing your controls expertise throughout your courses, I think the
world of controls now!
Thanks also to Dr. Chopra, Dr. Flatau, Dr. Baeder and Dr. Baz, the dean’s
representative, for being part of my committee. I appreciate your time and feedback.
I would like to thank the AVL folks as well. Greg, thanks for being a colleague
in the lab with whom it was always great to work with, thanks for your help and work
collaboration. Thanks Hector, Badri and Alex for your help at different instances
in time. Julia, thanks for being a great office mate and for helping out as well.
I would also like to thank Dr. Winkelmann and his lab, for supporting my
projects throughout the years on many fronts: hardware, supplies, lab space, etc.
The time as Dr. Winkelmann’s TA in my second semester was definitely very
formative for me in my career at UMD. I would like to thank Dr. Pines as well, for
having taken me as his TA when I started out with Masters. I learned a lot that
first semester! To Dr. Flatau, for her Comsol and hardware support. The work
with Dr. Flatau for Masters layed out the sensors foundation for the present work.
Would like to thank Anand, and Marc Stora, from the Alfred Gessow Rotor-
craft center on campus, as well as its director, Dr. Chopra. Thank you Anand and
Marc for the work with wing properties testing in your lab. Thanks also to Mike
iii
Sytsma, who helped us out at REEF in Florida with our wind tunnel tests as well
as Andrew Kehlenbeck for helping out during the first round of tests. Thanks Max
for wing fabrication and specs.
Would like to send a special thanks to Derrick Yeo, who not only piloted the
aircraft for the flight tests but also did such an amazing job with the plane daq code
and daq hardware. Thanks a lot Derrick!
Thanks also to James and Camli for starting out to collaborate on the CFD
structures simulations. I hope that we can resume that work at some point.
Thanks to Ahmad Kasaee at the wind tunnel for the brief test of our setup
and Evandro for the override box and hardware checkups.
I would like to thank also folks at the aerospace office; Lavita, Laura, Mike,
Otto and Tom. Without y’all we couldn’t possibly do our work!
iv
Contents
List of Figures viii
List of Abbreviations xiv
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Bio-Inspiration . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Strain feedback for gust alleviation . . . . . . . . . . . . . . . 5
1.1.3 Proprioceptive platforms . . . . . . . . . . . . . . . . . . . . . 7
1.1.4 Contributions in the perspective of the state of the art . . . . 9
1.1.5 Organization of the dissertation . . . . . . . . . . . . . . . . . 11
2 Theoretical Framework 13
2.1 Bio-inspired approach and distributed sensing . . . . . . . . . . . . . 14
2.1.1 Proprioception in nature . . . . . . . . . . . . . . . . . . . . . 15
2.1.1.1 Proprioception in humans . . . . . . . . . . . . . . . 15
2.1.1.2 Proprioception in birds . . . . . . . . . . . . . . . . . 19
2.1.1.3 Proprioception and strain sensing in insects . . . . . 21
2.1.2 Location and distribution of proprioceptors . . . . . . . . . . . 26
2.1.3 Proprioception in insect wings and flight control . . . . . . . . 28
2.1.3.1 Location and distribution of sensilla on wings . . . . 31
2.1.4 Distributed sensing in natural fliers . . . . . . . . . . . . . . . 33
2.1.5 Wave propagation . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.1 Atmospheric disturbances . . . . . . . . . . . . . . . . . . . . 37
2.2.2 Gust alleviation . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.2.3 Static and dynamic loading of wing . . . . . . . . . . . . . . 43
2.2.4 Aeroelastic models . . . . . . . . . . . . . . . . . . . . . . . . 47
2.2.5 Strain sensing . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.3 Strain feedback and controller framework . . . . . . . . . . . . . . . . 51
2.3.1 Strain feedback . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.2 Robust controller framework . . . . . . . . . . . . . . . . . . . 55
3 Theoretical considerations, Sensor placement and simulations 62
3.1 Theoretical approach to strains vs disturbance force differentials . . . 63
3.1.1 Roll Moment . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
v
3.1.2 Pitch Moment . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.1.3 Vertical Force . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.1.4 6-DOF model . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.2 Strain based wide field integration (WFI) . . . . . . . . . . . . . . . . 78
3.2.1 1D strain WFI . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.2.2 2D strain WFI . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.3 Sensor placement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.3.1 Sensor placement techniques . . . . . . . . . . . . . . . . . . . 83
3.3.2 Theoretical Framework EFI and Grammian . . . . . . . . . . 87
3.4 1D Sensor Placement Analysis . . . . . . . . . . . . . . . . . . . . . . 90
3.4.1 Wing FEM model . . . . . . . . . . . . . . . . . . . . . . . . . 90
3.4.2 Sensor placement using MPF/gradient and EFI/Grammian . . 91
3.5 Gust simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
3.5.1 ASWING aeroelastic model simulation . . . . . . . . . . . . . 99
3.5.1.1 ASWING nonlinear simulation -3D gust . . . . . . . 100
3.5.1.2 ASWING linear simulation -1D gust . . . . . . . . . 101
3.5.2 Fluid structure interaction simulations- Static Aeroelastic (Com-
sol) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
3.6 2D sensor placement analyses (Comsol) . . . . . . . . . . . . . . . . . 111
3.6.1 PWing test section . . . . . . . . . . . . . . . . . . . . . . . . 111
3.6.1.1 Wing FE model and sensor placement analysis . . . 111
3.6.2 Flight test PWing . . . . . . . . . . . . . . . . . . . . . . . . . 123
3.6.2.1 Wing FE model and sensor placement analysis . . . 125
4 Experimental section and analysis 130
4.1 Glider experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
4.1.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . 130
4.1.2 Controls Setup . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.1.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . 136
4.1.4 Force/Strain Feedback using Proportional Gains . . . . . . . . 136
4.1.5 Aileron System Identification . . . . . . . . . . . . . . . . . . 144
4.1.6 Static Proportional Model . . . . . . . . . . . . . . . . . . . . 144
4.1.7 Dynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . 148
4.2 Quadrotor experiments . . . . . . . . . . . . . . . . . . . . . . . . . 152
4.2.1 Strain Sensing . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
4.2.2 Sensor design for quadrotor arm . . . . . . . . . . . . . . . . . 158
4.2.3 Heave disturbance rejection . . . . . . . . . . . . . . . . . . . 163
4.3 Wing structural properties testing . . . . . . . . . . . . . . . . . . . . 166
4.3.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . 166
4.3.2 Bending and torsion properties . . . . . . . . . . . . . . . . . 169
4.4 Wind tunnel tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
4.4.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . 174
4.4.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . . 181
4.4.3 System identification of output matrix . . . . . . . . . . . . . 191
4.4.4 Discrete gust tests . . . . . . . . . . . . . . . . . . . . . . . . 195
vi
4.5 Flight testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
5 Summary and Conclusions 208
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
5.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214
5.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217
Appendix A Strain sensors for aircraft wing applications 219
Appendix B ATG characteristics 224
Bibliography 228
vii
List of Figures
1.1 Commercial and military applications of UAS: a. Amazon prime air,
a proposed multiple-rotor delivery system, b. Google project wing, a
first generation R&D fixed wing delivery prototype, c. Aurora Orion,
a long endurance unmanned aircraft system for military applications. 2
1.2 High bandwidth inner loop strain feedback . . . . . . . . . . . . . . . 7
1.3 Pwings test articles for wind tunnel testing: UMD PWing (OTS foil
gauges), Stanford PWing (Interconnected sensor network), Aurora
PWing (fiber optic sensors) . . . . . . . . . . . . . . . . . . . . . . . 8
2.1 Spinal reflex order and components [23] . . . . . . . . . . . . . . . . . 16
2.2 Spinal reflex in the human body, in particular with stretch proprio-
ceptive receptors and corresponding afferent/sensory neurons [23] . . 17
2.3 Stretch proprioceptive sensors in humans: a. muscle spindles and b.
Golgi tendons [24] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Avian muscles for flight control [26]. . . . . . . . . . . . . . . . . . . . 20
2.5 Different types of cuticular sensilla. a. sensillum trichoideum (hairs),
b. sensillum basiconium, c.sensillum campaniformium, d.sensillum
placodeum, e. sensillum coeloconium, f. sensillum ampullaceum, [35]. 25
2.6 Proprioceptors of the foreleg of periplaneta. The numbers in brackets
show the number of sensilla in each group. Ellipses show orientations
of campaniform sensilla [36]. . . . . . . . . . . . . . . . . . . . . . . . 26
2.7 Distribution of camaniform sensilla in the stick insect’s body [37]. . . 27
2.8 Campaniform sensilla on grasshopper wings . . . . . . . . . . . . . . 30
2.9 Distribution of camaniform sensilla on the wingblases of flies of dif-
ferent families of Diptera [39]. . . . . . . . . . . . . . . . . . . . . . . 32
2.10 Wave propagation measured with strain gauges: bimaterial beam
(left) and experimentally recorded strain histories (right) [46] . . . . . 35
2.11 Order of timescales in disturbance sensing . . . . . . . . . . . . . . . 36
2.12 Turbulence production in the atmosphere [49] . . . . . . . . . . . . . 37
2.13 Typical 1-cos gust shapes [50] . . . . . . . . . . . . . . . . . . . . . . 38
2.14 Ratio of Von Karman and Dryden PSD functions [52] . . . . . . . . . 39
2.15 State space equations for rigid body motion and structural dynamics 48
2.16 6-DOF coupled aeroelastic equations . . . . . . . . . . . . . . . . . . 49
2.17 Augmented attitude controller: load feedback to optimize the re-
sponsiveness of the attitude controller and to reject disturbances and
errors in actuated body torque (excerpt from Thompson et al. [16]). . 52
viii
2.18 Strain Feedback implementation. . . . . . . . . . . . . . . . . . . . . 55
2.19 Fast inner loop for disturbance rejection . . . . . . . . . . . . . . . . 56
2.20 MIMO general case with input and output disturbances and noise . . 59
2.21 High bandwidth inner loop strain feedback . . . . . . . . . . . . . . . 61
3.1 Calculation of wide field integrated constants . . . . . . . . . . . . . . 66
3.2 Balance of bending moments for calculation of disturbance roll moment 67
3.3 Pitching moment at root of aircraft as a function of shear stresses
along the wings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.4 Different pitching moments at root as a consequence of maneuver or
external inputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.5 Vertical force on aircraft and shear forces on the wings . . . . . . . . 73
3.6 6-DOF coupled model for force/strain feedback . . . . . . . . . . . . 75
3.7 Formulation for the strain feedback MIMO controller based on SISO
strain balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.8 6-DOF formulation for the strain feedback MIMO controller based on
system identified output matrices . . . . . . . . . . . . . . . . . . . . 76
3.9 6-DOF formulation for the strain feedback MIMO controller with
uncertainty in the Cmatrix . . . . . . . . . . . . . . . . . . . . . . . . 77
3.10 Strain wide field integration . . . . . . . . . . . . . . . . . . . . . . . 82
3.11 1D dimensional aluminum beam eigenstrain superposition (left) and
2D aluminum beam eigenstrain superposition (right). . . . . . . . . . 86
3.12 First and second finite element eigenmodes of UAS wing . . . . . . . 91
3.13 UAS normalized eigenstrains from finite element simulation . . . . . . 92
3.14 UAS normalized weighted eigenstrains using equal weights and MPF 93
3.15 Linearly weighted and MPF eigenstrain gradients . . . . . . . . . . . 94
3.16 UAS wing MPF based sensor placement using weighted eigenstrains
and gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.17 UAS wing MPF based sensor placement using weighted eigenstrains
and gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.18 UAS wing EFI/Grammian sensor placement . . . . . . . . . . . . . . 96
3.19 Comparison of EFI/Grammian and MPF/Gradient sensor placement
methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.20 ASWING definitions (excerpt from reference [72]) . . . . . . . . . . . 100
3.21 Top view of aircraft model for ASWING simulation . . . . . . . . . . 101
3.22 Airplane vertical velocity and bank angle with and without gust, non-
linear ASWING simulation . . . . . . . . . . . . . . . . . . . . . . . 102
3.23 Strain feedback used in ASWING reduced order simulation . . . . . . 104
3.24 Roll disturbance on aircraft states . . . . . . . . . . . . . . . . . . . . 104
3.25 Impact of roll disturbance on sensors . . . . . . . . . . . . . . . . . . 105
3.26 Roll angle at different strain feedback proportional grains K, with
derivative D=1 and integral I=2 gains. . . . . . . . . . . . . . . . . . 105
3.27 Wind tunnel simulation of fixed wing with impinging vertical gust . . 108
3.28 Wind tunnel simulation of fixed wing with impinging vertical gust,
strains are produced as a result of the impinging vertical gust . . . . 109
ix
3.29 Wind tunnel simulation of fixed wing with impinging rotational gust:
Fluid pressure profiles (left) and wing strains (right) . . . . . . . . . . 110
3.30 Principal strains for different impinging gusts: vertical gust (left) and
torsional gust (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
3.31 Airfoil geometric characteristics . . . . . . . . . . . . . . . . . . . . . 112
3.32 Normal element size mesh of carbon layup wing . . . . . . . . . . . . 113
3.33 First six eigenstresses of carbon fiber wing . . . . . . . . . . . . . . . 115
3.34 2D view of first six eigenstresses of carbon fiber wing . . . . . . . . . 116
3.35 Superposition of first six eigenstress modes (top view) . . . . . . . . . 118
3.36 Strain gauge stretchable network in spanwise orientation (Structures
and Composites Lab SACL, Stanford University) [8] . . . . . . . . . . 119
3.37 Sensor distributions based on SACL lab stretchable sensor network . 120
3.38 First six principal strain directions of carbon fiber wing . . . . . . . . 122
3.39 Pwing Airfoil and full wingspan geometry . . . . . . . . . . . . . . . 123
3.40 Flight Pwings for assembly with Apprentice Aircraft . . . . . . . . . 124
3.41 FE model of flight test wing and meshed geometry . . . . . . . . . . 125
3.42 First few eigenstrain modes of flight wing . . . . . . . . . . . . . . . . 126
3.43 Eigenstrain modes of flight wing: linear weighting (left) and MPF-v
weighted (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
3.44 MPF-v weighted eigenstrain map for 8 sensors . . . . . . . . . . . . . 127
4.1 Experimental setup of glider and launching mechanism . . . . . . . . 131
4.2 Glider plate and strain gauges . . . . . . . . . . . . . . . . . . . . . . 132
4.3 Experimental setup of glider and launching mechanism . . . . . . . . 132
4.4 VICON Setup: Vehicle marker setup(left) and VICON camera setup
and vehicle (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
4.5 Gust encountered by glider. . . . . . . . . . . . . . . . . . . . . . . . 135
4.6 Controls used in strain feedback tests. . . . . . . . . . . . . . . . . . . 135
4.7 Classical controller used for controller types comparison. . . . . . . . 136
4.8 Glider baseline roll angle curves: with no gust perturbation (top),
and under a vertical asymmetric gust (bottom). . . . . . . . . . . . . 137
4.9 Glider mean roll for free flight and with vertical gust perturbation. . . 138
4.10 Glider mean roll for proportional gains K=3, K=4, K=5, and K=6 . 139
4.11 Glider roll angle and strain feedback at gains K=3, 4, 5 ,6 and baseline
flight. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.12 D2 metric between the perturbed data with no applied control and
the one with applied strain feedback. . . . . . . . . . . . . . . . . . . 140
4.13 Controller comparison for roll gust rejection. TF= torque feedback,
P=proportional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.14 Glider trajectories: in free flight, under the influence of a gust and
using strain feedback. . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.15 Impact of weight on glider roll angle: Roll angle for glider weight and
a payload of 4.4 grams (top) and Roll angle for glider weight and a
payload of 10.9 grams (bottom) . . . . . . . . . . . . . . . . . . . . . 143
4.16 Vicon camera setup for system identification . . . . . . . . . . . . . . 144
x
4.17 Step input of wing deflections and corresponding output of aileron
deflections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
4.18 Aileron deflection vs wing deflection for K=1. Each point is the result
of averaging local points that correspond to the duration of the step
input. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
4.19 Filtered data, where aileron output has been drawn around its mean
value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
4.20 DJI Flamewheel quadrotor with strain sensing arm . . . . . . . . . . 153
4.21 (a) Quad arm with amplifier circuit, (b) Static testing of RPM quadro-
tor arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
4.22 (a)Hysteresis at a loading and unloading cycle, (b) Creep behavior of
RPM quad arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
4.23 (a) Photopolymer prototyped arm sections for MTS testing and (b)
CAD model of the quadrotor arm without optimized cross-section . . 161
4.24 Finite element simulations of stress under maximum expected quadro-
tor load: (a) Solid cross section, (b) Square hollow cross section and
(c) Smoothed square hollow section, (d) striated cross section, (e)
load cell cross section, (f)bi-rectangular cross section . . . . . . . . . 162
4.25 Strain sensing implementation: (a) Quad arms with different cross
sections and (b)Quad arm vehicle implementation . . . . . . . . . . . 163
4.26 Fast inner loop for disturbance rejection . . . . . . . . . . . . . . . . 164
4.27 Strain feedback performance metric . . . . . . . . . . . . . . . . . . . 164
4.28 Heave Velocity and Altitude Time History for Impulse Heave Distur-
bance for Varying Feedback Gain from Strain Gauge Sensing . . . . . 165
4.29 Height gauge setup for calculation of bending stiffness . . . . . . . . . 167
4.30 VICON setup for calculation of bending stiffness . . . . . . . . . . . . 168
4.31 Clamps used for static testing and Vicon wing setup . . . . . . . . . . 169
4.32 Load vs. tip displacement curves for estimation of EI with height gauge170
4.33 VICON processed data . . . . . . . . . . . . . . . . . . . . . . . . . . 171
4.34 Strain gauge top half of full bridge . . . . . . . . . . . . . . . . . . . 175
4.35 Wing with sensor network . . . . . . . . . . . . . . . . . . . . . . . . 175
4.36 Strain gauge amplification electronics . . . . . . . . . . . . . . . . . . 176
4.37 REEF wind tunnel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
4.38 Wind tunnel wing setup . . . . . . . . . . . . . . . . . . . . . . . . . 178
4.39 Load cell axes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
4.40 Test setup with ATG . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
4.41 6-axis sample data from load cell during steady test, mechanical
pulses were used to synchronize data . . . . . . . . . . . . . . . . . . 181
4.42 Calculated forces and moments vs wind tunnel speed for angle of
attack) degrees and angle of attack 5 degrees . . . . . . . . . . . . . . 183
4.43 Calculated lift and drag coefficients from wind tunnel data . . . . . . 184
4.44 AG35-r Drela aifoil, from [109], rotated airfoil and lift/drag coeffi-
cients vs angle of attack. . . . . . . . . . . . . . . . . . . . . . . . . . 185
4.45 Ramp test to simulate sustained gusts (20 kts, 25 kts, 15 kts, 20 kts),
AoA =5 deg: load cell forces and moments and strain sensor data. . 186
xi
4.46 Representative of interpolated strain rms measurements for steady
wind velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
4.47 Interpolated strain rms measurements for steady wind velocity at
AoA 0 deg and AoA 5 deg . . . . . . . . . . . . . . . . . . . . . . . . 188
4.48 Comparison of max. rms strains during steady runs . . . . . . . . . . 189
4.49 Comparison of strains produced with wing tip loadings of 100gr and
200gr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
4.50 ATG produced Lift force with and without strain feedback . . . . . . 190
4.51 Setup for Static system identification in the lift direction . . . . . . . 194
4.52 Stanford wing wind tunnel gust setup . . . . . . . . . . . . . . . . . . 195
4.53 Static manual forcing with coupled components . . . . . . . . . . . . 196
4.54 Static manual decoupled forces and moments Fx, forces and moments
(left), sensor data (right) . . . . . . . . . . . . . . . . . . . . . . . . . 196
4.55 Static manual decoupled forces and moments Fy, forces and moments
(left), sensor data (right) . . . . . . . . . . . . . . . . . . . . . . . . . 197
4.56 Static manual decoupled forces and moments Fz, forces and moments
(left), sensor data (right) . . . . . . . . . . . . . . . . . . . . . . . . . 197
4.57 Cmatrices for manual Fx, Fy, Fz calculated with least squares esti-
mator and colored noise . . . . . . . . . . . . . . . . . . . . . . . . . 198
4.58 Cmatrices for design points: AoA 9 deg, V=15kts (left) and AoA 3
deg, V=20 kts (right). . . . . . . . . . . . . . . . . . . . . . . . . . . 198
4.59 AoA=3 deg, V=20 kts, gusts, open loop . . . . . . . . . . . . . . . . 199
4.60 AoA=3 deg, V=20 kts, gusts, closed loop feedback S2 . . . . . . . . . 200
4.61 AoA=3 deg, V=20 kts, gusts, closed loop feedback C average . . . . 201
4.62 AoA=3 deg, V=20 kts, gusts, closed loop feedback C matrix (sysID) . 202
4.63 Filtered AoA 3 deg, V=20 Fy data to measure performance of feed-
back loop, d represents a decrease in the shift for the first gust which
is proportional to the shift observed in the other gusts . . . . . . . . 203
4.64 Apprentice aircraft with modified instrumented wing . . . . . . . . . 203
4.65 Strain signals and IMU recording and feedback for flight tests. . . . . 204
4.66 Apprentice aircraft and ground station during flight test preparation . 204
4.67 Apprentice flight tests . . . . . . . . . . . . . . . . . . . . . . . . . . 205
4.68 Open loop flight data, eight bridges, 32 sensors . . . . . . . . . . . . . 205
4.69 Flight data open loop, differential strain for sensors 4 and 5 and
aircraft roll rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
4.70 Flight data closed loop, differential strain for sensors 4 and 5 and
aircraft roll rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
4.71 Flight data closed loop, differential strain for sensors 4 and 5 and
aircraft roll rate for flight section. . . . . . . . . . . . . . . . . . . . . 207
A.1 Strain sensors for aircraft wings . . . . . . . . . . . . . . . . . . . . . 221
A.2 Strain sensors characteristics comparison chart . . . . . . . . . . . . . 222
A.3 Strain sensors advantages and disadvantages chart . . . . . . . . . . . 223
B.1 ATG turbulence grid . . . . . . . . . . . . . . . . . . . . . . . . . . . 224
xii
B.2 ATG turbulence statistics with wind tunnel throttle at 43.5 % and
at x/M=12.5, from [1] . . . . . . . . . . . . . . . . . . . . . . . . . . 225
B.3 ATG turbulence statistics mean core velocities: a. Mean core velocity
at 7.21 m/s, b. Mean core velocity at TIx=3.9%, c. Mean core
velocity at TIx=17.9%, d. Mean core velocity at TIx=20.7% . . . . . 227
xiii
List of Abbreviations
P Force acting on one wing
Fv Total vertical force
EI Bending stiffness
εR Strain on the right wing
εL Strain on the left wing
∆ε Strain differential
cw Wide field integration constant
x Span wise coordinate
y Normal coordinate
z Thickness
H(s) Aileron transfer function
Kp Proportional constant
Tp First order constant
Td Time delay
ys Sensor outputs
Φfs Matrix of target modes
q Response vector of target modes
w Measurement noise
Q Fisher information matrix
Yo Observability Grammian
R Intensity matrix
Ψ Matrix of eigenvectors of Q
Γ Matrix of eigenvalues of Q
EDi Effective Independence
εi Strain at location i
M Bending Moment
qw Distributed force acting on one wing
qi Discrete force acting at location i
Cw Wide field integration matrix
a Local distance to fuselage
xi Spanwise location of sensor i
L Roll moment
ML Left wing moment
MR Right wing moment
Lw Wing length
RPM Rapid prototype material
IMU Inertial measurement unit
U.A.S. Unmanned aerial system
xiv
Chapter 1: Introduction
1.1 Motivation
Unmanned aerial vehicles, systems and platforms are a growing field with in-
creasing number of applications in areas where human life needs to be protected,
areas which are of difficult access, unexplored goals or simply to perform tasks at a
reduced cost compared to a human operator. Therefore some of the most promis-
ing and useful applications include: reconnaissance, aerial photography, search and
rescue, transportation of supplies within nearby locations, civilian infrastructure
maintenance and disaster prevention and mitigation.These tasks need to be carried
in a way that is safe for the general population and the environment.
All these advantages however come at the price of increased instability, power
constraints, payload limitations, reduced flight endurance, among other challenges.
Autonomous flight is even more challenging, as ideally the aircraft will know how to
avoid obstacles and stabilize itself amid a number of disturbances of the environment.
In urban environments, avoidance of obstacles may be problematic, not only because
of gusts, but because of ground-effect and related flow field phenomena when flying
in the vicinity of walls and buildings. Light weight and slow flight speed of small
air systems creates profound complications in aerodynamic efficiency and agility [1].
1
(a) Amazon Prime Air (b) Google Project Wing
(c) Aurora Orion
Figure 1.1: Commercial and military applications of UAS: a. Amazon prime air, a
proposed multiple-rotor delivery system, b. Google project wing, a first
generation R&D fixed wing delivery prototype, c. Aurora Orion, a long
endurance unmanned aircraft system for military applications.
2
Low wing loading, which is necessary for aerodynamic efficiency and reasonable
flight endurance, exacerbates the problem of gust response. These platforms are
also prone to aeroelastic or inertial coupling with the rigid body dynamics and the
flexible structures, since many of these platforms are much more flexible than their
large scale counterparts. Coupling may also occur with vortex shedding and other
flow-separation phenomena [2].
Due to the promising advantages, both academia and industry have explored
ways to address the challenges that these vehicles pose. The reduced scale for flight
dynamics makes these smaller platforms more unstable and more susceptible to
atmospheric disturbances.
1.1.1 Bio-Inspiration
This work is inspired in the fast reflexive mechanisms found in nature as ap-
plied to stabilization of small UAS platforms. In particular, strain sensing has not
been fully explored in the sensing features presented herein. Furthermore, proprio-
ceptive strain sensing mechanisms have been found in numerous examples in nature.
Starting from the mamalian muscle spindles, golgi tendons and ruffini corpuscles,
then going to herbst corpuscles in birds, and sensilla of insects, nature is rich in
fast and reflexive mechanisms which allow for bioinspired controls. Insects not only
have strain sensors but also possess a number of visual sensors, tactile sensors, and
airspeed sensors for physiological and survival tasks. In particular, we are focus-
ing on the strain sensing features of flying insects, as their distributed sensing is
3
constructed to aid in flight and navigation.
Campaniform sensilla is the closest mechanism to that of strain sensing in
engineering. It’s structure is such that in-plane deformations get translated into out
of plane campaniform displacements, which activate the corresponding neuron. The
distribution of these organs coincide for the most part with the location of insect
wing veins, or the beams and trusses of insects. Therefore indicating that these
organs benefit from bending and stresses produced during flight. It has also been
observed that insects have these organs distributed along their wings and not just
at the base of joints, which is also another popular location for sensilla along the
insect body. These sensors allow the animal to monitor rate of flapping by firing
neurons on each wing flap. These signals are combined with other sensing organs
for flight stabilization. It has been observed as well that they use the campaniform
close to the tip of their wings to monitor the passage of torsional waves during the
flapping cycle. It is this precisely what inspires this work. The goal is to anticipate
the onset of disturbances by measuring the passage of the disturbance by means of
proprioceptive mechanisms such as those found in nature. The early detection of the
disturbance is what would allow for faster reaction to it. The reaction mechanisms
are proprioceptive in the sense that they are based on relative measurements to
other parts of the insect body and/or wings.
Structural inputs are much faster than rigid body inputs to any controls system
located at the fuselage of any aircraft. The reason being is that the frequencies are
much higher than rigid body ones. However, for certain aircraft these aeroelastic
modes become coincidental and therefore lead either to static or dynamic aeroelastic
4
divergence. If aeroelastic divergence is assumed not to be problematic for the system
at hand, then there is an opportunity to use these waves for controls purposes.
1.1.2 Strain feedback for gust alleviation
Strain signals had been used in the past for load characterization studies, mode
suppression, vibration suppresion and others. The approach of this work however,
intends to utilize strain signals not to provide precise rigid body state estimates
but rather use them to access higher order states and correct for flight disturbances
before they reach the aircraft plant dynamics. That is to say that the controller is
not typical in the sense that full states are estimated but rather state differentials
by means of strain measurements. Forces and moments at the center of gravity are
unambiguously measured by inertial measurement units. Strain sensing in the focus
of this work does not intend to replace these units but rather complement attitude
control and provide fast inner loop stability augmentation reflected in more stable
inner loop states.
In this sense, strain feedback can provide this stability, albeit studies of modal
impact will most likely have to be undertaken to guarantee that the strain controller
will be able to not be affected by system vibrations. To this end, filtering could be
necessary as well as tuning of the filter to the particular expected modal frequencies.
Gust mitigation literature review [3] suggests that the trend for more struc-
turally efficient aircraft yields both lighter and more flexible aircraft which brings
challenges in sensitivity to atmospheric disturbances. These lighter more flexible
5
aircraft face two significant challenges: reduced separation between rigid-body and
flexible modes, and increased sensitivity to gust encounters due to decreased wing
loading and improved lift-to-drag ratios.
Disturbances can be discrete or continuous, and are generally modeled as be-
ing isotropic in the real atmosphere. Discrete gusts are deterministic, continuous
disturbances or turbulence is represented as a random variation of gust velocity with
frequency spectra. This work mainly focuses on addressing disturbance rejection for
discrete gusts based on structural inputs. The framework that will be presented
however can be extended to frequency domain or time domain disturbances with
some modifications.
Wing loading due to gusts can be idealized as one dimensional phenomena
however it is actually a fully three dimensional phenomena. In this sense, the
bioinspired mechanisms for gust rejection that will be presented take the three di-
mensional nature of gusts into account in that load balancing of control surfaces by
means of strain associated signals can yield a wealth of information for stability aug-
mentation. Albeit we are not measuring absolute forces, the strain measurements
are a result of a combination of forces and moments on the aircraft.
1D and 2D theoretical approaches are proposed for stabilization by means
of strain differentials or calculated force differentials, for implementation in fast
reactive loops. Sensor placement techniques are applied to 1D and 2D cases as
well and simulations are used in closed and open loop to test for these strain based
mechanisms [4].
A high bandwidth inner loop was tested in hardware [5] [6] to use strain
6
signals for fast reaction to disturbances. The controls loop corresponding to the
strain feedback is shown in Figure 2.21.
Figure 1.2: High bandwidth inner loop strain feedback
This feedback loop can be implemented in SISO and MIMO systems, and in
the case of the latter one a robust framework can be devised in order to take into
account uncertainties in the placement of strain sensors. The robust framework can
also help with design of performance characteristics.
1.1.3 Proprioceptive platforms
Four platforms with strain feedback were implemented and/or tested. The
glider hardware and avionics was provided by Aurora Flight Sciences [7], the quadro-
tor which is an AVL testbed, the UMD Pwings and the full scale Pwing for flight
testing with the Apprentice RC aircraft, were instrumented with strain gauges.
There were also two more Pwings test section platforms. One with a stretch-
7
able sensor network built by the SACL lab at Stanford University [8] and one which
was instrumented with fiber optic sensors by Aurora Flight Sciences [7]. These wings
are shown in Figure 1.3.
Figure 1.3: Pwings test articles for wind tunnel testing: UMD PWing (OTS foil
gauges), Stanford PWing (Interconnected sensor network), Aurora PWing
(fiber optic sensors)
The strain sensing implementation is presented in Chapter 3 with all the exper-
imental results from open loop and closed loop tests in the presence of atmospheric
disturbances.
8
1.1.4 Contributions in the perspective of the state of the art
• A bioinspired gust mitigation reactive mechanism was analyzed and imple-
mented. This is a different methodology altogether from the conventional
approach where gust alleviation is based on estimation of gust states or rigid
body state disturbances. In the proprioceptive perspective of this work, the
focus is on the relative strain measurements within the body of the aircraft
and how they can be used to aid in gust alleviation.
• In the theory and analysis aspect, this work contributes a differential strain/force
balance theoretical approach for pitch, roll and vertical force, where stability
is achieved by monitoring the balance of forces on the corresponding controls
surfaces.
• It also proposes another theoretical approach where if the forces and moments
are empirically correlated to the strain sensing patterns i.e. by least squares
regression methods, then stability can be achieved by monitoring the 6-DOF
force differentials at the root of each wing.
• Another theoretical development is the wide field integration for strain as
applied to disturbance rejection. A wide field integration for one dimensional
wing was implemented in simulation producing successful gust rejection.
• Eigenstrains can be weighted according to the expected directionality of the
gusts according to calculated mass participation factors or MPFs. These direc-
tional 2D weighting patterns can be considered to be basis functions for strain
9
measurements where a set of gains can be calculated to extract the different
types of gusts. In this work, it was found that the MPF weighted patterns
are strongly correlated to the observability Grammian, with the added direc-
tionality feature. Therefore the 2D energy weighting pattern which encodes
six different types of gusts can be considered to be a directional Grammian.
The basis functions for the strain wide field integration can be considered to
be directional Grammians.
• In terms of hardware, bioinspired strain sensor networks were implemented
on different platforms. These include a quadrotor, a UAV carbon fiber wing
section and a full carbon fiber wing for a small UAV for flight testing. Bioin-
spired strain sensing had not been previously applied to UAS platforms for
reflexive disturbance rejection in prior art. Previous developments had used
strain sensors mainly for monitoring of loads for health monitoring, vibration
suppression, wing shaping, among others.
• A MIMO multidegree of freedom framework was formulated with the proposed
strain dependent output matrices. The uncertainty in the output matrix and
estimation matrix can be modeled as either a structured uncertainty or an
unstructured uncertainty by tuning the weighting functions related to the C
matrix multiplicative uncertainty. This is such that the system will be nomi-
nally stable, have good performance characteristics such as good tracking and
good disturbance rejection, and can be robustly stable and have robust per-
formance amid the uncertainties. This is typically done through minimization
10
of the maximum singular value of the corresponding transfer function.
1.1.5 Organization of the dissertation
This work is organized as follows. An introduction is provided in the first
chapter as an overview of the entire dissertation.
In chapter 2, the theoretical framework is presented. It begins with the pro-
prioceptive bioinspired approach that was applied to the force disturbance sensing
and rejection mechanisms. Then it presents some mechanisms found in insects when
they use their wing strain sensors, or sensilla, to sense cuticle strain for many dif-
ferent purposes, including flight control. Then the most relevant gust models are
presented as well as disturbance rejection mechanisms that have been used in the
literature. After this is presented, a theoretical framework of the state of the art is
introduced and after that, the theory and analytical model for the force and moment
disturbance rejection. It presents the general intended framework in which all the
strain feedback concepts are applied and it provides an overview of the general equa-
tions employed in the aeroelastic simulations used later on in chapter 2. The robust
controller framework is also introduced for output matrices with uncertainties.
Also in chapter 3, the theoretical aspect of strain feedback is presented. These
equations are based on the bioinspired approach presented in the previous chapter.
A wide field integration for roll is presented for a 1D case. After this, the sensor
placement analysis and results are presented for a 1D case. Gust simulations in open
loop for the 2D CFD/structures simulation is presented as well as the closed inner
11
loop full 6-dof linear aeroelastic gust simulation obtained using ASWING. For a 2D
case, the results are presented and applied to the PWings (Proprioceptive wings)
test section and flight wing section.
In chapter 4, the experimental section and analysis is presented. A small
UAS glider, a mini quadrotor, a wing section and a full wing were outfitted with
strain gauges, results and analysis are provided. These platforms were tested under
sharp disturbances including an actuator disturbance, a 1-cos disturbance, a turbu-
lent wind disturbance, and a discrete finite gust disturbance. All these platforms
which had strain sensing reactive mechanisms showed good disturbance rejection
properties.
In chapter 5, a summary of all the tests and results is presented and some
conclusions, as well as future work related to the developments of this work are
envisioned.
12
Chapter 2: Theoretical Framework
This section introduces the bioinspired motivation for strain sensing as applied
to flight control of fixed wing aircraft.
Bioinspired sensory mechanisms have shown to have multiple applications in
ground navigation [9] as well as for UAS aircraft [10]. Feedback laws motivated by
biological approaches can be used to reject gusts. Biological systems have natural
strain gauges to aid navigation and stability [11]. These are triggered by stimuli
caused by locomotion, changes in posture and external events. Examples of sens-
ing encoders found in nature include the campaniform sensilla of insects and the
slit sense organs of arachnids [12]. Insect wings have cuticular sensors which help
them respond to environmental inputs. Campaniform sensilla, the load/strain sen-
sors of insects, can induce a response in muscles that is an order of magnitude
faster than that achieved by the vision system. These stress/strain sensors are
directly connected to dendrites of the primary afferent neurons allowing a quick
reflexive response [13]. The possibility of a low latency response is the initial moti-
vation to use strain feedback for gust rejection. We will refer to strain feedback and
force/torque/moment feedback interchangeably throughout this work.
In this work we not only explore the conventional approach to sensor place-
13
ment but also a wide field integrated approach [9] based on structural inputs. This
approach is motivated by the mechanisms observed in nature. Strain, stress and load
sensing mechanoreceptors are found in all species for navigation and motor control.
Examples of sensing encoders found in nature include the cuticular sensors of insects
and the slit sense organs of arachnids [12]. In insects, the stress/strain sensors are
directly connected to dendrites of the primary afferent neurons allowing them to
have a quick reflexive response [13]. The cuticle of the spider has a large number of
slit sense organs, which respond to forces causing strain in the cuticle. These slits
are connected to bipolar neurons which reflexively transmit the chemical synaptic
potentials directly to the central nervous system [14]. A low latency response is
possible by feedback of structural quantities which are directly correlated to higher
order states. A distributed sensing approach can reduce the error covariance in the
measurements while reducing feedback latency [9].
Strain sensing found in nature can be imitated by constructing sensing mech-
anisms of maximal desired output. A method of sensing strains, that may encode
force information, based on the campaniform sensilla was developed by [?] and [15]
to be embedded in space structures. Controllers based in these strain mechanisms
have also been proposed and simulated for attitude control augmentation [16].
2.1 Bio-inspired approach and distributed sensing
The present work aims to harness the bio-inspired mechanisms of propriocep-
tion both in function and in functional morphology. Proprioception in biological
14
entities is a constantly running natural function which determines the body’s posi-
tion and the way it is moving through space. This mechanism allows for the brain
to determine the body’s state of motion and position. For a better understanding of
the concept of proprioception, this overview will go through the corresponding mech-
anisms found in humans and will then translate and compare to the mechanisms
found in insects and small fliers.
2.1.1 Proprioception in nature
2.1.1.1 Proprioception in humans
Proprioception, from the latin propius, is the perception of ones own motion.
It is the sense of the relative position of body segments in relation to other body
segments. Proprioceptive organs provide the central nervous system with informa-
tion about the relative position and movement of body parts [17]. Proprioception
allows for indication of whether the body is moving with the appropriate effort and
where the various segments of the body are located in relation to each other [18].
It allows for the understanding of where each joint is at that moment in time. This
relative position sensing is different from exteroception which entails sensitivity to
stimuli originating outside of the body and it is different from interoception which
is sensitivity to stimuli originating inside of the body. Conscious proprioception or
kinesthesia allows the brain to bring to a conscious or knowing level where the body
is in 3D space. Unconscious proprioception allows the brain and more specifically
the cerebellum to know unconsciously where the body is in 3D space [19]. Much
15
of the afferent information generated by self-initiated movements does not reach
consciousness.
Kinesthesia usually refers to the perception of one’s own motion whereas pro-
prioception not only entails the perception of movement but also awareness of the
body’s orientation in space [20]. This cognitive ability is channeled from the dif-
ferent proprioceptive receptors through neuronal afferents to the central nervous
system, as shown in Figure 2.1. Spinal reflexes only travel to the spinal chord and
back and are therefore subconscious. Reflexes are different from reactions, which are
voluntary responses to stimuli from the environment. The reflex time, or the time
between the onset of a stimulus and action of the effector, is determined mainly by
the conduction time in afferent and efferent pathways. The conduction velocities of
human nerve fibers are given by the speed at which a given electrochemical impulse
propagates down a neural pathway [21]. In humans, the reaction times vary depend-
ing on which one of the senses is participating in the reflex and on a number of other
biological factors. On average, the reflex time for visual inputs is 0.2 seconds [22].
The reflex arc is mainly composed of the elements shown in Figure 2.1.
Figure 2.1: Spinal reflex order and components [23]
16
Figure 2.2: Spinal reflex in the human body, in particular with stretch proprioceptive
receptors and corresponding afferent/sensory neurons [23]
Figure 2.2 shows the components of a reflex arc through muscle stretch pro-
prioceptive receptors. The receptor is the organ that responds to the stimulus by
converting it into a neural or electrical signal. Then the afferent neuron carries the
signal to the central nervous system and encounters either association neurons (in-
terneurons) or efferent neurons. The outcome of the signal transfer depends on the
complexity of the reflex mechanism [23]. The output neurons or efferent neurons
then carry the electric signal to the organ that is to respond to the stimulus. Re-
ceptors convert all the sensory information into action potentials or impulses which
are ultimately relayed to the cerebellum.
In humans, the proprioceptive organs entail the muscles, joint capsules, liga-
ments, tendons, skin and the internal ear. These organs provide the central nervous
system with information about the position and movement of body parts. The two
17
types of muscle receptors are the muscle spindles and the Golgi tendon organs, as
shown in Figure 2.3. Muscle spindles are the principal proprioceptors [20], these are
receptors which provide feedback information on the degree of muscle stretch for the
muscle reflex mechanism. They can sense muscle length (limb position) and rate
of length change (velocity or sense of movement) by means of the primary endings
or afferents. Golgi tendon organs are receptors located at the junction of muscle
and tendon and provide mechanosensory information to the central nervous system
about muscle tension. Joint kinesthetic receptors are located within and around
the articular capsules of synovial joints. They provide feedback information on the
degree and rate of angulation of a joint. They act as limit detectors in that they op-
erate mainly at the extremes of the normal range of biologically possible joint angles.
Skin receptors have a similar function to that of muscle spindles. Skin mechanore-
ceptors are activated when the skin suffers rotation or compression. Proprioceptors
in the internal ear function for equilibrium.
Figure 2.3: Stretch proprioceptive sensors in humans: a. muscle spindles and b. Golgi
tendons [24]
18
Information about limb position and movement however is not generated by
individual sensing receptors but by populations of afferents [20]. Afferent signals
are generated during motion and are processed to encode limb endpoint positions.
The afferent inputs are then referred to a central body map to determine the spatial
location of the limbs. Afferent nerves are peripheral nerves that transmit to the
spinal cord and ultimately the brain. These nerves therefore carry out a sensory
function. Injuries to these may decrease proprioception as well as diseases of the
neurologic system.
Proprioception also entails the ability to reposition a joint to a predetermined
position. It gives an estimation of the relative position of parts of the body as well
as the effort being employed in movement [17]. Proprioceptive senses include the
senses of position and movement of limb the sense of effort, sense of tension or force
and the sense of heaviness and the sense of balance [20]. These are always associated
to motor commands. In the normal limb tendon organs and possibly also muscle
spindles contribute to the senses of tension or force and heaviness. In particular, we
would like to apply the mechanisms for the sense of force and balance in the present
work.
2.1.1.2 Proprioception in birds
Birds also possess a number of mechanoreceptors, some of which have propri-
oceptive functions. These are located near the follicles of feathers allowing birds to
detect turbulence [25], as well as in muscles for flight control. In avian flight, birds
19
and bats, possess modified forelimbs which contain muscles that control wing shape
by moving the bones of the arm and hand (Figure 2.4).
Figure 2.4: Avian muscles for flight control [26].
Wing receptor types in birds can be free nerve endings or encapsulated endings.
Free endings sense temperature and pain, while encapsulated endings are mechanore-
ceptors or sensory corpuscles. Kestrels for instance, can withstand elevated levels
of turbulence while hovering and maintaining a spatial reference posture. There are
a number of corpuscles, out of which the Herbst and the Ruffini corpuscles have a
flight control function in addition to other functions. Herbst corpuscles are widely
distributed in the skin and are associated with feather follicles and with the muscles
of the feathers. Ruffini corpuscles are axon endings in close contact with collagen
fibers , probable serving as stretch receptors and numerous in-joint capsules [27].
Herbst corpuscles sense vibrations and are directionally sensitive. These are
clustered in the leading edge of birds wings and the alula, or the freely moving first
20
digit. These organs are tuned to certain wind frequencies, which is also observed in
insects and their analogous proprioceptive mechanoreceptors. It has been suggested
that this frequency is a way to measure airspeed. Flight steering muscles of birds are
controlled not only indirectly via neck reflexes but also in parallel, directly by the
central nervous system. This combination of direct and indirect control of wing and
tail muscles can be coordinated or switched off [28]. Stretch receptors in limb muscles
are also potential load sensors and could be auxiliary to the cutaneous receptors in
determining gust induced displacements of the wing. Feather displacements are also
detected by additional sets of mechanoreceptors.
2.1.1.3 Proprioception and strain sensing in insects
This work is inspired in the proprioceptive mechanisms of humans and other
vertebrates as mentioned in the previous sections, but finds most of its inspiration in
the fast reflexive mechanisms observed in insects. In particular, those mechanisms
which use wing mechanosensory information for flight stabilization. Insects are not
limited by the cellular biophysics muscle contraction limits. Insects have simpler
neurological mechanisms as compared to vertebrates and also have been observed
to have faster reflexes. For instance, the visual startle reflex of a Condylostylus fly
has been measured in the order of 5 ms [29].
Insects receive and respond to a wide variety of mechanical stimuli, they are
sensitive to physical contact with solid surfaces, they detect air movements, includ-
ing sound waves and they have a gravitation sense [30]. Some of these mechanisms
21
are exteroceptive, interoceptive and some proprioceptive. Proprioception in insects
has a different structure altogether as compared to humans and birds, mainly due
to the different physics associated to smaller scales. This information is gathered
by a spectrum of mechanosensilla. These are sense organs that are able to respond
continuously to deformations and stresses in the body, providing the insect with pos-
ture and position information. They send information to the insect central nervous
system about mechanical strains in the cuticular exoskeleton [31].
Mechanoreceptors can be found in many places on the surface of an insects
body Figure 2.6. These receptors are innervated by one or more sensory neurons that
fire in response to stretching, bending, compression, vibration or other mechanical
disturbance [32]. Arthropod mechanoreceptors are divided into two morphological
groups: type 1 or cuticular and type 2 or multipolar. Type 1 are ciliated receptors,
associated with the cuticle and have their nerve cell bodies in the periphery close
to the sensory endings, they can be subdivided into three major groups: Hair-like
receptors, campaniform sensilla, and chordotonal organs. Type 2 are nonciliated
neurons whose central cell bodies have many fine dendritic endings. Other types
of sensing mechanisms are also present in insects, including chemosensors. Sensilla
can be both exteroceptive as well as proprioceptive.
Five types of proprioceptors occur in insects: hair plates, campaniform sensilla,
chordotonal organs, stretch receptors, and nerve nets [30]:
(a) Hair plates:
Sensory hairs (sensilla trichoidea) occur on all parts of the insect body, however
are found in greatest concentration on parts of the body which come into
contact with the targeted sensing object. In its simplest form a sensillum
comprises a rigid, poreless hair set in a membranous socket and four associated
22
sensing cells. Hairs differ in sensitivity, long hairs respond to small forces and
pressure changes, whereas shorter thicker hairs (sensilla chaetica) respond to
larger stimuli. Hair plates on the legs contribute to the insects gravitational
sense.
(b) Campaniform sensilla:
A campaniform sensillum is similar to a tactile hair but has a dome shaped
plate of thin cuticle instead of the hair shaft. The plate may be slightly raised
above the surrounding cuticle, flush with it, or recessed but in all cases it
contacts the neuron and serves as a stretch or compression sensor. This plate
can be round or elliptical depending on the species. The elliptical plates have
a stiffening rod of cuticle along the longitudinal direction on the ventral side of
the plate. The plate has directional sensitivity to stress that is perpendicular
to the longitudinal axis of the rod. Sensilla are typically arranged in groups
with directional arrangements, i.e. such that rods are perpendicular.
(c) Chordotonal organs:
also known as scholophorous sensilla or scolopidia, are another widely dis-
tributed form of proprioceptor in insects. These differ from campaniforms in
that there is no specialized exocuticular component. They are associated with
the body wall, internal skeletal structures, tracheae and structures in which
pressure changes occur. They exist as strands of tissue that stretch between
two points. A change in length will produce stretching or bending of the den-
tritic membrane and will stimulate the sensing cell. Chordotonal organs of
arthropods mediate resistance reflexes that are similar to stretch reflexes of
vertebrates [33].
(d) Stretch receptors:
A multipolar neuron whose dendrites terminate in a strand of connective tissue
or a modified muscle cell, which is then connected to a body wall, interseg-
mental membranes and muscles. If the terminals of the neuron are stretched,
the receptor is stimulated. These provide information about rhythmically oc-
curring events within the insect, like breathing.
(e) Nerve nets:
Interconnected nerve fibers and cells are connected together to form a nerve
net. A nerve net provides a simple level of coordination. Most insect’s nerve
cells are clustered into groups which are called ganglia. Ganglia are intercon-
nected by nerve cords and receive impulses from sense organs and direct them
to other parts of the body, particularly muscles. Insects have multiple nerve
chords, this contrasts with humans and other vertebrates, as these only have
a single tubular one [34].
Amongst these types of proprioceptors, which involve insect organs which are
superficial and non-superficial, insects have a number of cuticular sensilla Figure 2.5.
23
Table 2.1: Some insect mechanoreceptors
Mechanoreceptor Description Function Location
Trichoform
sensilla
Hairs Tactile, motion
Behind head
on legs, near joints
,
Campaniform
sensilla
Flattened
oval discs
Flex receptors
Exoskeleton, legs,
base of wings
Stretch
receptors
Multipolar neurons Stretch receptors
Interseg. membranes,
muscular walls
Pressure
receptors
Hair like structures Pressure detectors Tracheal system
Chordotonal
organs
Subjenual, Tympanal,
Johnstons
Vibration,acoustic,
resonance
Legs, thorax,
pedicel of antenna
The cuticle is mostly superficial, however it’s also layered into endo-, meso-, and
exo- cuticle, with the tougher layer being the most external layer. The exocuticle
is very strong under compressive forces but comparatively weak under tension. The
endocuticle is flexible and is able to resist tensile forces. The cuticle is usually
characterized by surface patterns which reflect the arrangement of cells beneath
them. Sensilla found in the cuticle are the basic functional and structural units of
cuticular mechanoreceptors and chemoreceptors.
The different types of cuticular sensilla are found inFigure 2.5. These illus-
trate the mechanoreceptor contraptions of insects which use these organs to sense
tensile and compressive forces for survival and navigation. An example of the vari-
ous functions of these types of insect sensors is found in Figure 2.6, where different
24
types of sensilla serve different propriosensing purposes by functional location on
the insect’s body and by orientation within the locations. There are groups of cam-
paniform sensilla on the trochanter, a group proximally on the femur and another
on the tibia, and a small number on the dorsal surface. These are typically accom-
panied by other types of proprioceptors, as discussed above. Probably the major
function of campaniform sensillum is to register stresses produced by the weight of
the insect’s body on the limbs, the forces caused by the muscles on the cuticle, and
external stimuli [31].
Figure 2.5: Different types of cuticular sensilla. a. sensillum trichoideum (hairs),
b. sensillum basiconium, c.sensillum campaniformium, d.sensillum pla-
codeum, e. sensillum coeloconium, f. sensillum ampullaceum, [35].
25
Figure 2.6: Proprioceptors of the foreleg of periplaneta. The numbers in brackets
show the number of sensilla in each group. Ellipses show orientations of
campaniform sensilla [36].
2.1.2 Location and distribution of proprioceptors
Insect platforms have different performance requirements. Proprioceptive or-
gans in insects are located in many different places along the insects body, grouped
around regions which will produce maximal outputs, i.e. grouping around a joint.
Stretch proprioceptive organs are typically found near regions of maximal strain, as
in the case of the stick insect trochanter where the campaniform sensilla are located
at the leg joint Figure 2.7. In particular, stretch and stress sensors are found along
insect bodies and wings.
Insects distribution of proprioceptive sensors is different from that of humans.
In humans, the stretch sensors lie on skeletal muscle fiber and tendons. These
move with the limbs to generate motion and adapt to different situations. For
natural fliers, the muscle function however gets to be modified in actuating wings
26
Figure 2.7: Distribution of camaniform sensilla in the stick insect’s body [37].
for flapping and gliding.
The muscles of insect wings are however less complex than those of birds and
do not extend beyond the axilla, which is the point at which thoracic muscles are
attached to the insect wing and therefore control of shape must be exerted remotely.
That is, they transmit forces via rigid wing components such as veins and areas of
thickened membrane which act as levers. However, active muscular forces cannot
entirely control then insect wing shape in flight. The muscular forces interact dy-
namically with the aerodynamic and inertial forces that the wings experience during
flight and they also interact with the wings structural properties. The effects of these
interactions are determined by the architecture of the wing itself, its planform and
relief, the distribution and local mechanical properties of the veins, the local thick-
ness and properties of the membrane as well as the position and form of the lines of
flexion [38].
27
2.1.3 Proprioception in insect wings and flight control
Insect wings are muscularly driven aerodynamic surface that generate lift and
propel the insect through the air. However, the wing is also a sensory surface,
with a variety of sensory organs. Many of these sense mechanical and chemical
inputs. They can detect relative movement between wing and body and mechanical
deformation of the wing cuticle [39]. Similarly, as the sensilla distributed on the body
(Figure 2.6), the sensilla on the wings is also composed of hairs, campaniform sensilla
and chordotonal organs. The hairs are sensitive to external forces, the campaniform
and chordotonal organs are sensitive to self motion of the insect. Campaniform
sensilla are sensitive to cuticular deformation, and the chordotonal sensilla act as
stretch sensors, therefore providing proprioceptive information.
Wings also possess a variety of surface structures mainly for sensing purposes
[38]. These include tactile and proprioceptive sensilla, scales, microtrichia, a variety
of spines and papillae, among others. Sensilla on insect wings have been reported
previously in various species . Another insect organ that is well populated with
campaniform sensilla is the haltere. This is an organ that serves as a gyroscope.
Any departure from level flight will produce torques on the bases of halteres at
twice the vibration frequency [31]. The unique flight control equilibrium sense is
provided by the sensilla at the base of the halteres. This sense is a robust equilibrium
reflex in which angular rotations of the body elicit compensatory changes in both
the amplitude and stroke frequency of the wings [40]. Halteres encode the angular
velocity of the body by detecting the Coriolis forces that result from the linear
28
motion of the haltere within the rotating frame of reference of the fly’s thorax.
On wings, arrays of campaniform sensillae monitor the rate and extent of local
bending. These sensory receptors are found in arrays which are concentrated near
the front of the wing and at the wing base where stiffness and bending moments are
highest. The high surface density of campaniform sensilla and their direct neural
connections with central locations of flight control suggest that these sensors monitor
cyclical wing motions. Their spanwise topographical arrangement suggest a role in
monitoring progression of torsional waves during wing rotation [41]. On fly wings,
the campaniform are tuned to the wing beat frequency. Campaniform also respond
to deflections if camber is actively imposed to the wing, for instance during the
upstrokes and downstrokes. A correlation between number of wing sensillae and
flight maneuverability has been proposed.
Location of sensilla is species dependent. Grasshoppers, in particular the
Melanoplus sanguinipes, have groups of campaniform sensilla on the proximal area
of the subcostal vein of both fore and hind wings. These monitor the degree of
wing twisting. Figure 2.8 shows the veins and distributed sensilla along the wings
of a grasshopper. These figures show the number and distribution of sensilla on a
representative fore wing (left) and hind wing (right) of a grasshopper (melanoplus
sanguinipes). The primary veins are: Sc, subcostal, R, R1, Rs, radius; M, Ma, Mp,
media, Cu, cubitous [42]. In each wing, the main nerves follow the primary veins and
their branches [42]. The intervenal nerve nets are formed by axons from surrounding
sensilla, as these join to form the nerve branches of the main nerves. Both wings
of the grasshopper have trichoid and campaniform sensilla, which is identified as
29
chordotonal. These groups of campaniform sensilla are responsible for reflex control
of wing twisting, which is required for normal flight as well as body orientation and
stability during flight [31].
Campaniform are also present in the wings of a number of other insect species
such as Drosophila and Calliphora, as shown by Dickinson [43]. In Drosophila, the
distribution of campaniform is somewhat along specific veins near the leading edge of
the wing, this is shown in Figure 2.9. In the Dipteran wings shown, the campaniform
sensilla lies exclusively on two branches of the radial vein R1 and R4+5. These are
used for flight control and wing coordination.
Figure 2.8: Campaniform sensilla on grasshopper wings
30
2.1.3.1 Location and distribution of sensilla on wings
In insect wings, sensilla are typically located near veins, as these are the struc-
tural support of the wing and therefore encode bending and torsion.
Veins are typically cuticular tubes which contain hemolymph and often tra-
cheae and nerves. Hemolymph is the circulatory fluid of arthropods and is analogous
to blood and lymph in vertebrates. There is a distinction between main or longitu-
dinal veins which radiate from the base and cross-veins, which link the longitudinal
veins. These main veins can also contain trachea, nerves, and can also conduct oxy-
gen and sensory information. Aside from vital functionality they act as supporting
members, preventing rips from spreading across the membrane. Many longitudinal
veins are essentially beams, principally experiencing bending and twisting forces.
Like typical man-made cantilevered beams, their diameter tends to decrease from
the base to the tip, reflecting the fact that the bending moments on them do the
same. This span wise tapering has two other functions; it reduces the moment of
inertia, and hence the energy required as well as stresses generated during oscilla-
tion. It also makes the wing-tip more flexible and more easily deflected by excessive
inertial forces and unpredictable impacts. Junctions between veins vary depending
on the level of rigidity required for the function of the joint.
The location and distribution of sensilla is species dependent. It was shown
that campaniform sensilla is placed at specific locations on the wing of Drosophila
for the mainenance of control and phase relationships between the two wings. This
allows the insect to steer and turn as well [43]. The presence of sensilla on the wings
31
Figure 2.9: Distribution of camaniform sensilla on the wingblases of flies of different
families of Diptera [39].
to control wing motion is common across all insect species. For the grasshopper
wings, sensilla is located near the proximal end of both fore and hind wings. These
are used for monitoring the degree of wing twist.
Campaniform sensilla occur in veins and in intervenal membranes. Those in
membranes are surmounted by dome shaped plaques, those in veins have a less pro-
nounced dome. Both types are innervated by a single neuron with a large soma
and a short dentrite. The dendrite is inserted centrally inside the dome. Chor-
dotonal sensilla were noted at the base of both wings [42] of the grasshopper. In
both wings, the sensilla are concentrated along the longitudinal veins, only a few
are on transverse veins and in membranous areas. On the fore wing, the sensilla
are generally evenly distributed along the length of the wing, except for the wing
tip. Heavier concentrations occur along the trailing margin of the wing. Hair and
campaniform sensilla were observed mostly on the top surfaces of both wings of the
32
grasshopperFigure 2.8. Similarly for Diptera, the sensilla were mainly located near
principal longitudinal veins.
Campaniform sensilla are situated primarily in the membranes of the wings,
but always close to veins. Campaniform sensilla undoubtedly monitor the degree of
cuticular stress with changes in air pressure at various locations on the wings. They
are suitably located for this function on or near major veins and in clusters at the
base of the hind wing. These sensilla in the hind wing are essential in regulating
fore wing twisting during locust flight.
Venation patterns have an effect on the flexural stiffness of wings [44]. It has
been found that spanwise flexural stiffness is 1 to 2 orders of magnitude larger than
chordwise flexural stiffness. It has also been observed that there is an increased
density of veins towards the leading edge, which causes anisotropy. This anisotropy
would serve to strengthten the wing in bending while allowing chordwise bending to
generate camber. This could also facilitate spanwise torsion, which is seen in many
species during supination. Some insects have wings that are stiffer in both directions,
for instance dragonflies, hawkmoths, flies and bumblebees have wings that are stiffer
than expected for their size, whereas damselflies, craneflies and lacewings have more
flexible wings than expected.
2.1.4 Distributed sensing in natural fliers
It is in the perspective of strain sensing and it’s bioinspired functionality that
this dissertation is constructed upon. Similar to the Golgi tendons and joint recep-
33
tors in humans, as well as the sensilla organs, in particular those which sense in
plane tension. Strain sensors were installed in the proprioceptive perspective such
that they are able to sense relative motions of the aircraft. The location of the sen-
sors and their distributed nature resembles that of insect wings, in morphology and
in function. It has been shown that rapid reaction times can be achieved by means
of fast integration of distributed sensors, such as in the fly’s compound eye [9].
In the insect wings, natural stretch sensors are placed along wing venations,
or the structural components of insect wings. They are typically distributed along
the leading edge of the wings and concentrated around certain wing veins.
2.1.5 Wave propagation
This work is inspired not only in the proprioceptive mechanisms of nature but
also in the physics of propagation of mechanical disturbances. As discussed earlier,
some insects use the campaniform sensilla to sense the propagation of torsional
waves through their wings for flight control [41].
Mechanical wave propagation has not been explored for disturbance rejection
in the UAS aeronautical field, as typically not much information is obtained from the
aircraft structure itself while on flight. Studies by [45], have measured and modeled
the propagation of mechanical waves on a free free beam structure [46] by means of
accelerometers.
Strain sensing was used for impact testing [46] and measurement of wave prop-
agation related quantities Figure 2.10. The bimaterial beam shown is made of strips
34
of aluminum and epoxy. The impactor consisted of a steel ball which was dropped
from a short height on one end of the bimaterial beam.
Figure 2.10: Wave propagation measured with strain gauges: bimaterial beam (left)
and experimentally recorded strain histories (right) [46]
The propagation times depend on the material properties of the beam or wing.
The time difference from the maximum peak on strain gauge A and the maximum
peak of strain gauge B is the time of travel between the two points. In this case,
the propagation times are in the order of tens of microseconds as the material
for the beam is aluminum, which is very dense, rigid and therefore would yield
fast propagation speeds. Computational models have been built to study force
propagation , mainly for health monitoring purposes [45].
This work is based on anticipation of the rigid body motion by early detection
of disturbance propagation, used in a proprioceptive way. The following diagram
depicts the order of response to a disturbance in a fixed wing aircraft. First, discrete
35
or continuous disturbances will hit the wings, second, the mechanical disturbance
will propagate throughout the structure until it finally reaches the fuselage, where
IMUs will typically be recording rigid body motion. With strain feedback, antici-
pation is gained in the mitigation of the disturbance as reaction can occur before
causing the fuselage to be affected, or at least completely affected. Anticipation
times were measured for the glider experiments of Section (4.1.1).
Figure 2.11: Order of timescales in disturbance sensing
2.2 State of the art
Relevant state of the art is presented in the context of strain sensing as applied
to atmospheric disturbance rejection.
36
2.2.1 Atmospheric disturbances
Atmospheric phenomena is of many different origins Figure 2.12. Turbulence
producing phenomena can be due to sharp thermal changes, convection changes,
clear air turbulence, low level terrain induced turbulence, mountain wave turbu-
lence, cloud induced or convective induced turbulence, cumulus clouds, and in-cloud
turbulence [47], [48]. Thunderstorms produce the most severe type of turbulence.
Tip vortices from other aircraft may also produce air turbulence in the vicinity of
their flight path.
Figure 2.12: Turbulence production in the atmosphere [49]
Atmospheric wind can be considered a superposition of three types of flow:
mean wind, waves and turbulence [3]. The variation time of mean winds is a few
hours, for waves is tens of minutes and the timescales of turbulence range from
seconds to minutes.
37
Disturbances can be categorized into two ideal categories: discrete and contin-
uous, and are generally modeled as being isotropic in the real atmosphere. Discrete
gusts are deterministic and have defined forms. Three forms of discrete gusts are of
interest: step gusts or sharp-edged gusts, ramped gusts and 1-cosine gusts, which
capture the form of a solitary gust. Two typical 1-cosine gust shapes are shown in
Figure 2.13.
Figure 2.13: Typical 1-cos gust shapes [50]
Continuous disturbances or turbulence is represented as a random variation
of gust velocity with frequency spectra due to the irregular and anisotropic nature
of these wind variations. Dryden and von Karman distributions are two common
continuous turbulence distributions. Their spectral representation can be used to
generate turbulence by filtering band limited white noise with an appropriate shap-
ing filter [51]. Figure 2.14 shows a ratio of the psd functions of the Von Karman
and dryden continuous gust models. Both of these distributions are defined by
characteristic scale wavelength which is a function of altitude in the atmosphere
and the root mean square velocity. The von Karman distribution is often regarded
38
as a more accurate description of measured atmospheric distributions however the
Dryden distribution is simpler to implement in the time domain.
Figure 2.14: Ratio of Von Karman and Dryden PSD functions [52]
If uniform gusts are considered, the lift changes as a function of the angle
of attack induced by the gust for lateral and vertical uniform gusts. For head on
gusts, the change is in the dynamic pressure. These changes in lift are shown in
Equation (3.43), from [52].
∆L = ρ2V∞SCLαWg gustvertical
∆L = ρ2V∞SCLαVg gustlateral
∆L = ρV∞SCLUg gustlongitudinal
(2.1)
39
Where ρ is the air density, V∞ is the freestream velocity, S is the planform
area, and Ug, Vg, Wg are the longitudinal, lateral and vertical gust velocities. CLα
is the lift curve slope and CL is the lift coefficient.
Wing loading due to gusts can be idealized as one dimensional phenomena.
However, it has been shown that important effects arise from 3D modeling of gusts.
For instance, a comparison of 3D and 1D gust load modeling on the C-5A Galaxy
aircraft (Lockheed Martin, Bethesda, Maryland) showed that 3D gust modeling
reduces wing loads and increases tail loads increased during turbulence. This was
verified through flight test experiments [3]. In addition to this, uniform spanwise
gust distributions can predict larger loads than non-uniform gust distributions.
In nature gust rejection is effectuated by means of visuals [53], tactile feedback,
and other types of mechanosensory feedback. Halteres for instance keep the angular
momentum to act as a gyro [?]. In visual disturbance rejection, principles of optic
flow and neuronal integration produce patterns of motion which can be matched
to the different visual cues. This highly reactive mechanism uses low bandwidth
sensing organs in flying insects to quickly map out the nature of the disturbance.
2.2.2 Gust alleviation
Many different types of gust mitigation mechanisms have been devised and
investigated. These mechanisms include, amongst others, passive and active mech-
anisms. Passive gust mitigation mechanisms are incorporated at the design phase.
For instance adding wing dihedral adds roll stability and therefore mitigates gust
40
disturbances indirectly. Wings which have been designed to dampen out oscillations
from gusts are another type of passive solution to the gust problem. Other mech-
anisms have modified wings which passively react to gusts, decreasing the load on
the overall structure. An example of a passive mechanism includes having rotating
mid-sections which can let vertical gusts pass by [54].
Other mechanisms that have been proposed consist of cutting out sections
of the wing [55], enhancing spoilers, and varying the dihedral angle through wing
geometric changes. Articulated wings by means of passive wing deflection [56] and
inclusion of additional control surfaces such as winglets [57].
Active controls for gust mitigation can be grouped into three distinct groups:
operational aircraft, flight test experiments, and wind tunnel experiments [58].
Examples of operational aircraft active controls include: Active lift distri-
bution control (ALDCS), which uses wing tip accelerometers to determine wing
bending and torsional motion. Similarly the Active control system (ACS) [59] used
wingtip and fuselage forward and aft vertical accelerometers as well as fuselage pitch
gyroscopes, for load alleviation.
Structural mode control has been widely explored. Structural vibration sup-
pression can be found in many embodiments. One of them, the one implemented on
the B-1 Lancer, which actively suppresses the unfavorable motion at the pilot station
using dedicated active canard-like control surfaces and a co-located accelerometer.
Flutter suppresion has also widely been explored [60].
Wing morphing is another approach to vibration suppresion [61] as well as
adaptive controllers for recursive estimation of vibration modes [62]. Vibration
41
suppression can also be achieved by geometry changes of the airfoil [63]. Optimal
and adaptive controllers for gust alleviation have also been developed to achieve
gust mitigation [64].Another method entails the measurement of the angle of attack
computed as the difference between the measured angle of attack and the inertially
derived angle of attack [65], [66].
Gust load alleviation surfaces can be implemented in out of phase configura-
tions to damp out unwanted modes. Other methods actively control the rudder and
elevators to increase the fuselage damping response. Trailing edge flaps can also be
used to actively suppress the effect of gusts [67].
Optics mechanisms have also been proposed, such as LIDAR gust anticipations
by measuring the flow properties of the wind ahead of the aircraft [68] as well as
LIDAR variants such as the ultraviolet Doppler LIDAR for higher altitudes.
Wind tunnel experiments for gust alleviation have been similar to some of the
flight testing methods.
Most of the methods in the literature use inertial loading measurements for
gust alleviation, ride quality, and fatigue life extension of the aircraft.
The loading due to gusts is an important airworthiness requirement and is
regulated by the FAA, FAR [3], [52]. Requirementes on discrete, sharp gusts as well
as continuous turbulence can be found in [58].
42
2.2.3 Static and dynamic loading of wing
Aircraft are subjected to different forces while on flight, causing the airframe
to undergo structural deformations and internal loads. The total force is a sum-
mation of the inertial loads, the aerodynamic forces and the forces introduced by
atmospheric disturbances such as gusts and turbulence. In the case of gust induced
structural response, its flexibility will depend on either the amount of energy trans-
ferred from the gust disturbance to the structural bending modes or if energy is
absorbed from the gust, the dissipation of that energy by some form of damping.
One dimensional wing static loading models can be developed as a function of
Lift based on the available literature. Under the right assumptions, a wing can be
modeled as a cantilever beam and the strains generated by Lift and other forces can
be incorporated into a static force model.
A correlation of the strains produced by gusts and the gust characteristics
will allow the necessary control surface corrections to account for the lift changes
produced by the gust. Therefore gust disturbances can be aeroeleastically correlated
to strain profiles in the wing structure. In the simple case of a vertical gust Wg.
If ε = f (Wg) and Wg = g (L) →  = F ( Laero, Pinertial) (2.2)
Where Laero = Lg + Lm is the total lift or total aerodynamic loading, cor-
responding to the gust and the motion generated lift, respectively. Coriolis and
Centrifugal forces and their effect on the wing strain, are neglected for this initial
43
analysis.
If the wing is sufficiently slender, then strip theory can be applied and the
force per unit span due to the gust velocity encountered by the airfoil’s leading edge
at the instant t=sb/U can be calculated as follows [69], [70]:
Lg = 2piρUb
∫ s
0
Wg (σ)ψ
′
(s− σ)dσ (2.3)
Where U is freestream velocity. The distance traveled in semichord lengths is
defined as s = U ∗ t/br = U ∗ t/b, with br a reference chord value which equals b for
a uniform distribution. This expression assumes an equal spanwise chord length. If
the gust distribution is made constant instead of arbitrary, then the lift produced
by a sharp gust disturbance is an integral function of the Kussner’s function ψ.
Lg = 2piρUbWg
s
∫
0
ψ′ (s− σ) dσ (2.4)
The second part of the aerodynamic loading, Lm can be calculated as a function
of normal coordinates which represent the vibratory bending modes of the wing and
the Wagner functions for indicial lift φ [69]. Therefore,
Lm = f (εi, φi) , i = 1, 2, .., n (2.5)
Inertial terms can also be calculated in terms of these two parameters, the
44
integrand corresponding to inertial loads in a vibrating wing corresponds to:
Iloads =
(
n∑
i=1
φi (η) ε
′′
i (s)
)
m (η) (2.6)
Using Euler-bernoulli assumptions for the wing approximated in its mechanical
behavior to that of a beam, the calculation of strain at a vertical distance z from
the neutral axis and at a distance x from the root of the wing or the center line can
be calculated as follows:
εx =
z
EIy
M =
z
EIy
l∫
x
[Lg + Lm − Iloads] (η − x) dη (2.7)
Where E is the Youngs modulus of the wing, Iy the moment of area in the y
direction at the centerline and l the length of the wing. This static output equation is
based on the nature of the loads that the wing experiences and can be implemented
in a feedback loop to account for the gust disturbances. If dynamic effects are
considered, the state vector for longitudinal motion now depends on the truncated
bending modes. It can be defined as follows, i.e. for the first two modes [71]:
x = [α, q, ε1, ε˙1, ε1, ε˙1], (2.8)
And the output is a function of the modified state vector and the measured
strains. This state vector would then be formulated in the setting of an aeroelastic
system of equations, as shown in Section (2.2.4). Another option for this calculation
45
consists of a quasistatic force analysis for the duration of the gust.
The forces on a beam produce stresses and strains. Assuming that the wing
can be modeled as a simple cantilever, the lift per unit span would equal the shearing
force per unit span on the beam. The lift per unit span on a finite wing of wing
length L and a total wingspan of 2L is:
L′ = ρ∞V∞Γ0
√
1−
(x
L
)2
(2.9)
The lift is elliptical due to the effects of the downwash created by finite behavior
of the flow over a finite wing. V∞ is the freestream velocity and ρ∞ is the density
of the freestream. Γ0 is the circulation at the root of the wing and is equal to:
Γ0 =
V∞SCL
pi
(2.10)
Integrating the lift per unit span or shearing force per unit span q to get the
total shearing force as a function of wing span x:
∂Fs
∂x
= q (2.11)
∂Fs
∂x
= L′ (2.12)
46
Integrating to get the shearing force caused by the elliptical lift distribution:
Fs = ρ∞V∞Γ0
∫ (
1−
x2
L2
) 1
2
dx (2.13)
Fs = ρ∞V∞Γ0




x
2
√
1−
x2
L2
+
log
(
x
√
−1
L2 +
√
1− x
2
L2
)
2
√
−1
L2



 (2.14)
Using this expression of the shearing force, the bending moment can be found.
The shearing force is the rate of change of bending moment:
Fs =
dM
dx
(2.15)
where the bending moment is:
M(x) =
∫
Fsdx (2.16)
and the bending strain is calculated using:
x = −z
M(x)
EI
(2.17)
2.2.4 Aeroelastic models
The coupling of structural and aerodynamic forces acting on can be described
by aeroelastic models. These can be non linear, linear, dynamic, static or servoelastic
47
[70].
Aeroelastic models are fundamented in a system of equations which consist of
rigid body dynamics and structural dynamics. This system of equations is shown in
Figure 2.15. In this equation AR ,BR, CR ,DR are the linearized rigid body equations
of the aircraft and AE ,BE are the linearized structural dynamics equations. All
these matrices contain stability derivatives which can be identified with analytical
or numerical methods. The elastic degrees of freedom correspond to the modes
associated to each eigenfrequency.
Figure 2.15: State space equations for rigid body motion and structural dynamics
The fully coupled 6-DOF aeroelastic equations are presented in Figure 2.16,
in these equations, the elastic to rigid body aerodynamic coupling matrix is AER
and the rigid body to elastic aerodynamic coupling matrix is ARE.
ASWING is an aeroelastic simulation environment developed by [72]. In this
48
Figure 2.16: 6-DOF coupled aeroelastic equations
package, A fully nonlinear Bernoulli-Euler beam representation is used for all the
surface and fuselage structures, and an enhanced lifting-line representation is used
to model the aerodynamic surface characteristics. The lifting-line model employs
wind-aligned trailing vorticity, a Prandtl-Glauert compressibility transformation,
and local-stall lift coefficient limiting to economically predict extreme flight situa-
tions with reasonable accuracy. ASWNG also allows for linear models to be con-
structed about reference flight condition points generated from a list of possible
reference points.
2.2.5 Strain sensing
Strain gauges have been used for load calibrations [73]. These calibrations
permit the measurement in flight of the shear or lift, the bending moment, and
the torque or pitching moment on the principal lifting or control surfaces, landing-
49
gear structures, and relatively complex built-up wing and empennage assemblies.
This investigation did not use these signals for control, but rather to monitor loads
during flight with the objective of monitoring structural integrity demonstrations,
and aiding developmental flight testing.
Many applications have strain bridges have been used to measure bending and
torsional loads on the wing [74]. These however are only monitoring the bending
and torsion of the wing and not used in feedback gust alleviation, an accelerometer
was used in the active control of the full scale model instead.
Several wind tunnel models have been fitted with strain gauges for load mon-
itoring. A feedback implementation is realized by [75] where a wind tunnel airfoil
model was fitted with strain gauges at the root of an airfoil. Then a PID controller
used measuring strain as a control input and the flap acceleration as the controller
output .
Wind turbine blades have also been fitted with strain gauges to study fatigue
load reduction potential when applying trailing edge flaps. Strain gauges are used
as input for the flap controller such that local small fluctuations in the aerodynamic
forces can be alleviated by deformation of the airfoil flap [76].
Wing shape control is another field where strain gauges are used. It was shown
by [77] in a numerical simulation study, that the shape of the wings can be changed
by integrated strain actuators. This allows for tasks such as gust load alleviation
(i.e., disturbance rejection),flutter and vibration suppression (i.e., plant regulation),
and maneuver enhancement (i.e., command following) and redistribution of maneu-
ver loads. This approach albeit requires special piezoelectric actuation. Other strain
50
actuation applications include induced strain actuators, such as piezo ceramics and
electrostritives [78], which directly deform the structure through electromechanical
coupling terms that appear in the actuator constitutive relations. Other aeroelastic
wings have used strain gauges in closed loop to activate flaps [79] in wind tunnel
settings for gust load alleviation.
Aeroelastic controls have been extensively explored [80]. Control laws applied
to aeroelastic gust alleviation and flutter suppression include both MIMO [79], SISO
laws [75] , as well as PID [81], LQG [79], LQR [82], H∞ controllers [83], among
others.
2.3 Strain feedback and controller framework
Previous work related to the measurement of structural parameters of flexible
aircraft [84] has had a variety of goals and applications, i.e. enhance aircraft per-
formance, decrease loads and reduce vibrations [85]. In particular, the estimation
of flexible structural motion for gust load alleviation has been implemented in the
form of adaptive and non-adaptive control technologies [62]. Control laws associ-
ated to these systems relate the motion of control devices to measurements taken on
the vehicle [86]. The various schemes found in the literature utilize different struc-
tural motion measurements to activate control devices through the proper control
laws [87]. Most of these techniques focus on the aeroelastic and aeroservoelastic in-
teractions of aircraft and typically aim to suppress vibrations and prevent the onset
of flutter [62]. In UAV’s these effects are much more critical due to the relative
51
closeness of the rigid body modes and the structural modes [88]; flexible UAV’s
have a much more agile flight potential as compared to large aircraft [84]. However,
most of these approaches address vibration suppression. In this work we do not
aim to suppress structural vibrations but rather use structural signals to generate a
reflexive reaction to gusts. Reflexive strain feedback has not been widely explored.
Strain measurements can influence the attitude control design [16]. For instance,
regulation of the actuation commands sent to the control surfaces through augmen-
tation of the attitude controller. A second way in which strain measurements can
be implemented is to directly use them to detect angular acceleration to potentially
obtain a more optimal tracking response [16].
Figure 2.17 shows a load feedback loop to optimize the attitude controller and
reject disturbances.
Figure 2.17: Augmented attitude controller: load feedback to optimize the respon-
siveness of the attitude controller and to reject disturbances and errors
in actuated body torque (excerpt from Thompson et al. [16]).
52
2.3.1 Strain feedback
In this work, strain is measured in order to provide an estimation of the forces
and moments acting on the aircraft. Strain measurements are taken with metal film
strain sensors. These measurements are then projected or wide field integrated [9]
to obtain the actual force state estimations. Given our main focus is the extraction
of the force states, we will assume a quasi-static approach. To illustrate the state
extraction procedure a 1D cantilever beam model is used. For a simple cantilever
beam the second order ode shown in Equation (4.10) describes the behavior of the
displacement with respect to span x.
d2
dx2
(EI
d2y
dx2
) = Q (2.18)
where EI is the bending rigidity, y is the vertical displacement and Q is the
applied force. If Q is a uniformly applied force, then the strain as a function of span
location x is given by Equation (4.12):
ε(x) =
z
2EI
∗ f(x) ∗Q (2.19)
And z is the distance from the neutral axis, L is the wing semi-span. At the
root, of the right wing for instance, the strain simply becomes Equation (4.13):
εR =
zL2
2EI
∗Q (2.20)
53
The uniform vertical force applied to the wing can be calculated as Equa-
tion (4.11):
Q =
2EI
zL2
∗ εR = C ∗ εR (2.21)
Therefore, the vertical force Fv acting at the root of the wing is calculated
by simply scaling the strain sensor output. In the case of a fixed wing aircraft
and assuming no inertial considerations, the total force on the aircraft rigid body
becomes:
Fv = 2Q = 2CεR (2.22)
Fvlevelflight = 2CεR = 2CεL (2.23)
Fv = 2C(εR − εL) = Cw ∗∆ε (2.24)
Where 2C equals the integrated output or wide-field integrated constant Cw.
For level flight the left and right total strain output is approximately the same.
An imbalance in forces produced by a gust disturbance is then proportional to
the difference between the corresponding strains. The force/torque calculation is
proportional to the strain differential. This basic calculation can be extended to
encompass a chord wise structural response as well as to describe the extraction of
54
the remaining forces and moments. An implementation of the strain feedback is
shown in Figure 2.18.
Figure 2.18: Strain Feedback implementation.
2.3.2 Robust controller framework
A robust control scheme is proposed in Figure 4.26. In this proposed multi-
variable control scheme, the input disturbance on the plant output can be mitigated.
Therefore, the disturbances are absorbed by the inner loop before they reach the
plant.This controller was demonstrated for accelerometer and strain on a quadrotor
by [6] in correcting roll and heave perturbations.
This feedback improves tracking of requested controls in the presence of dis-
turbances by regulating errors between desired torques and actual torques. Thus the
system is simultaneously capable of rejecting external disturbances, such as gusts,
and internal disturbances, such as actuator variations.
In order to quantify the performance of this controller, robust analysis is used.
The performance is quantified by means of input and output transfer functions
55
Figure 2.19: Fast inner loop for disturbance rejection
between the different signals of interest. Typical signals of interest include input
and output disturbances, measurement noise, control commands, and states.
In the robust control framework, most signals live in bounded energy spaces,
such as the Lebesgue space, which is the space of signals f(t) with finite or bounded
energy, as shown in Equation (2.25).
Ln2 (−∞,∞) = {f : (−∞,∞)→ R
n,
∫ ∞
∞
||f(t)||2dt <∞} (2.25)
If a topological relationship is given to the Lebesgue space a number of spaces
can be defined, among these, the Hilbert space. This is a space which has an inner
product defined and is complete, meaning that the Cauchy sequences all converge
within the space. The inner product has an induced norm, which gives a measure
of the size of elements and in this case, of the gain of a particular transfer function.
The Lebesgue space induced norm is defined as follows from Equation (2.26), where
56
the inner product is defined in Equation (2.27).
< f(t), g(t) >=
∫ −∞
∞
||f(t)||2dt <∞ (2.26)
||f(t)||2 = (< f(t), f(t) >)
1/2 (2.27)
Now applying the previous structures to a space of transfer functions we
will obtain a transfer function space with bounded energy. This space is the Lˆ∞
space, which is a space of transfer functions that maps signals from L2(−∞,∞) to
L2(−∞,∞).
Lˆ∞ = {G : ||G||∞ <∞} (2.28)
For a stable real rational transfer function G, with input w and output z =
Gw, we may represent the induced operator norm in Ln2 of G as:
||G||L2→L2 = sup
w 6=0
||z||L2
||w||L2
= sup
w 6=0
||Gw||L2
||w||L2
(2.29)
This represents the input energy to output energy gain. This norm is equiva-
lent to the infinity norm, which is defined as the supremum of the maximum singular
value over all frequencies ω. The induced norm in Ln2 relates to the H∞ norm || · ||∞
57
as
||G||L2→L2 = ||G||∞ = sup
ω
σ¯[G(jω)] (2.30)
where σ¯[G(jω)] is the maximum singular value of G as a function of frequency
ω.
The H∞ space is a subset of the Lˆ∞ space, in that the transfer functions need
to be real, and therefore stable, as well as rational. The norm of the H∞ space is
then defined as:
||G||∞ = sup
ω
σ¯[G(jω)] (2.31)
where σ¯ is the maximum singular value .
Thus the input signal energy to output signal energy gain may be represented
via the maximum singular value σ¯, and we may further note that the worst case
scenario input energy to output energy gain is given by the H∞ norm, or equivalently
sup
ω
σ¯.
In a MIMO setup, the effect of the input and output disturbances on the
plant output can be quantified by means of different transfer functions. These
transfer functions can be defined from block diagrams which describe the systems
being modeled. In a general scenario, the block diagram for a system with input
disturbance, output disturbance, and measurement noise is described by the diagram
shown in Figure 2.20.
58
Figure 2.20: MIMO general case with input and output disturbances and noise
From this diagram the following equations for the output y and controller
input u can be obtained:
y = Tor + SoGdI + Sodo − Ton
u = KSor + SIdI −KSodo −KSon (2.32)
where the open loop transfer functions depend on where the loop is broken,
as in MIMO systems the transfer functions do not commute (Equation (2.33)).
Lo = GK So = (I + Lo)−1 To = (I + Lo)−1 breakoutput
LI = KG SI = (I + LI)−1 TI = (I + LI)−1 breakinput (2.33)
When quantifying the performance of a MIMO system in terms of the qualities
of the transfer functions in Equation (2.34) the norm of these transfer functions is
59
considered for desirable qualities. The maximum singular value σ¯ and the minimum
singular value act as collective gains of the MIMO system at every frequency ω and
can be used for controller design. The singular values are also interpreted as a ratio
of the input signal energy over the output signal energy, for each corresponding case.
The desired range of values for the transfer functions associated to good track-
ing and disturbance rejection are described as follows in Equation (2.34):
y = SoGdI → σ¯(SoG) small
y = Sodo → σ¯(So) small
y = −Ton → σ¯(To) small
u = KSor → σ¯(KSo) small
e = Sor → σ¯(So) small
(2.34)
The first expression in the collection of partial dependencies as shown in Equa-
tion (2.34) requires that the maximum singular value of SoG be small for the fre-
quencies of interest in order to have attenuation of input disturbance over the output
y. The second expression needs the maximum singular value of So be small for at-
tenuation of output disturbances on the output. Similarly, for the third expression,
the max. singular values of the complimentary output sensitivity function need to
be small in order for the measurement noise to have less of an impact at desired
frequencies. This analysis applies as well to other signals, such as the last expression
which allows for avoidance of large control signals due to reference commands.
In this framework, we would like to introduce strain feedback in a high band-
60
width inner loop, as show in Figure 2.21. In this fast loop, which is emphasized
with a red dotted box, we use the structural response scales which are faster than
the rigid body scales. Strain patterns encode instantaneous forces and moments
being applied to the vehicle. Spatial weightings improve the signal to noise ratio by
means of least squares estimated higher order states. The latency is reduced and
no dynamic filtering is required. This feedback loop ensures that what is applied at
the vehicle is what the controller asks for in the presence of disturbances.
Figure 2.21: High bandwidth inner loop strain feedback
This controller is different from conventional ones in that the fast inner loop
precedes the plant and absorbs input disturbances to the plant, preventing the plant
from receiving highly perturbed inputs. This approach however also reduces com-
manded inputs, albeit one way to correct for this is to feedforward the commanded
inputs to ensure that they will be properly fed to the plant. This approach was
proposed by [89] for implementation in a quadrotor platform.
61
Chapter 3: Theoretical considerations, Sensor placement and simu-
lations
A quasistatic one dimensional approach based on strain sensing is an approach
to gust rejection, detecting imbalances to correct for disturbances.
The forces on the airplane’s fuselage are a result of the applied forces on the
aerodynamic and control surfaces. These forces propagate along the material out
of which the wings and empennage are made. The speed of propagation of these
mechanical inputs is given by the composition of the aircraft. A mechanical wave
traveling on a softer material will take more time than one traveling in a material
with high stiffness. Wave propagation is a natural phenomena which some insects
use for flight control, like the diptera. Wing twist propagation is detected by clusters
of campaniform sensilla, as it was shown in the previous chapter.
The quasistatic approach basically needs to obtain the distributed sensing
profiles and match filter the wide field integrated constants to a particular flight
mode. Disturbances will then be rejected when the interrelated strain parameters
go beyond certain thresholds. This is a proprioceptive approach in that differences
among different sensors allow for a disturbance rejection mechanism.
62
3.1 Theoretical approach to strains vs disturbance force differentials
3.1.1 Roll Moment
Assuming we can model the wings of the aircraft as a one dimensional can-
tilever beam, we have that the strain at any point x along the wingspan can be
calculated as:
ε(x) =
z
EI
M(x) (3.1)
If we have a set of strain measurements distributed along the wingspan we can
calculate the discrete bending moments at each point:
Mi =
EI
z
εi (3.2)
Using the static Euler beam equation we can relate the moment and shearing
force applied locally at every sensor position:
d2
dx2
(EI
d2w
dx2
) = q(x) (3.3)
And since the moment can also be calculated as follows and considering only
changes in the x-direction:
M(x) = EI
d2w
dx2
(3.4)
63
we have that:
d2
dx2
(EI
d2w
dx2
) =
d2
dx2
(M(x)) =
EI
z
d2ε(x)
dx2
= q(x) (3.5)
And therefore calculate the force per unit span that generates each strain:
qi =
d2Mi
dx2
=
EI
z
d2εi
dx2
(3.6)
Each load can be calculated as follows by using the strain at three points
by using the following numerical scheme, which is a result of applying the three
point central difference formula for the approximation of the second derivative of a
function at a given point [90].
d2εi
dx2
=
εi(x+ h)− 2εi(x) + εi(x− h)
h2
(3.7)
Therefore the load at each point i can be calculated as:
qi(x) =
EI
z
εi(x+ h)− 2εi(x) + εi(x− h)
h2
(3.8)
Where the load at x is calculated by measuring the strain at three different
locations on the beam centered at x. The value of h depends on the hardware
constraints in terms of connections and space required to attach each sensor to the
wing surface. By calculating the bending moments on each sensor location and
calculating the shear force per unit span, a distribution of the shear force per unit
64
span can be obtained as shown in Figure 3.1. In order to calculate the moment
produced by the calculated point loads qi and considering a discrete set of n sensors:
MT =
n∑
i=1
qixi (3.9)
Where xi indicates the position of the sensors along the wing span. This
distance has an offset equal to half the width of the fuselage if the moments are
calculated with respect to the center of gravity of the aircraft. If the fuselage has
width 2a, the moment arm to the center of gravity becomes:
xi′ = xi + a (3.10)
and the moment becomes:
MT =
n∑
i=1
qi(xi + a) (3.11)
In the case of having four strain sensors, the moment on one wing can be
calculated as:
MT =
4∑
i=1
qi(xi + a) = q1x1 + q2x2 + q3x3 + q4x4 + a ∗ (q1 + q2 + q3 + q4) (3.12)
In the case of 6 sensors per wing then this equation becomes:
65
Figure 3.1: Calculation of wide field integrated constants
ML,R =




ML
MR



 =




CLq1 C
L
q2 C
L
q3 C
L
q4 C
L
q5 C
L
q6
CRq1 C
R
q2 C
R
q3 C
R
q4 C
R
q5 C
R
q6
























q1
q2
q3
q4
q5
q6




















(3.13)
Where the wide field integration matrix (Section (3.2)) becomes:
Cw =




CLq1 C
L
q2 C
L
q3 C
L
q4 C
L
q5 C
L
q6
CRq1 C
R
q2 C
R
q3 C
R
q4 C
R
q5 C
R
q6



 (3.14)
and its components can be calculated depending on the chosen locations along
66
the wingspan of each wing xi and the local distance to the fuselage ai:
Cw =




a1 + x1 a2 + x2 a3 + x3 a4 + x4 a5 + x5 a6 + x6
a7 + x7 a8 + x8 a9 + x9 a10 + x10 a11 + x11 a12 + x12



 (3.15)
In this case we have considered that the first sensor location is at the wingtip
of the left wing. The application of the output matrix allows for calculation of the
resulting moment on each wing and therefore the differential rolling moment of the
aircraft or disturbance rolling moment is given by:
Ld = ML −MR (3.16)
An equivalent C matrix can be obtained to extract the rolling moment differ-
ential from the measured strains. Figure 3.2
Figure 3.2: Balance of bending moments for calculation of disturbance roll moment
The sensors will then be placed in a strain feedback loop using proportional
gains and weighted using the wide field integrated constants that were derived. A
67
more elaborate controller using H-infinity analysis can also be implemented for roll.
3.1.2 Pitch Moment
Assuming we can model the torsional properties of an aircraft fixed wing by
means of the shear strain on the surface of the wing, we can calculate the total
pitching moment applied to the center of mass of the aircraft. For simplicity pur-
poses we will assume that the cross section of the wing is uniform and that it has
a uniform second moment of area along the span wise direction. The wing is also
assumed to have no sweep angle or dihedral. In this case the torque at any point of
the wing (beam) can be calculated as shown in Equation (3.17).
T = GIp
dθ
dx
(3.17)
Where G is the shear modulus of the beam at a particular location x = L and
θ is the shear deflection angle at every location x. This equation is related to the
shear strain at that particular location by means of Equation (3.18).
T = GIp
dθ
dx
= GIp
(γ
r
)
(3.18)
Where γ the shear strain and r is is an approximate torsion radius about a
torsion axis. In the case of a wing, this axis is also known as the moment elastic axis.
There are more detailed formulations of this mechanism; however, we will focus on
the simpler approach in this discussion. The pitching moment induced by aircraft
68
controls has an immediate effect on the structural strains of the aircraft wings. This
effect is time dependent but it can be approached in the way of a quasi-static torque
balance. This pitching moment which starts at the center of mass of the aircraft,
generates a torque pattern on the aircraft wings. The quasi-static solution entails
the summation of the torques generated along the wingspan.
Conversely, if these torques are measured along the wingspan, then the varia-
tion in the pitching moment can be calculated using shear strain sensors at specific
locations. In the case of a pure pitching motion, the total torque will be equivalent
in both wings. However, this balance may not apply for combined motions such as
pitch-yaw motions or pitch-roll motions.
Considering one wing, the total pitching moment applied corresponds to Equa-
tion (3.19):
M = T1 + T2 + T3 + TN (3.19)
In the case of N shear sensors distributed at locations xi along the wingspan.
Where each torque is described by Equation (3.18).
If these torques are measured by shear strain sensors then the total pitching
moment will be given by Equation (3.20).
M = T1 + T2 + T3 + . . . TN =
GIp
req
∗ (γ1 + γ2 + γ3 + . . . γN) (3.20)
It is worth mentioning that shear strain sensors need to be setup in pairs
69
Figure 3.3: Pitching moment at root of aircraft as a function of shear stresses along
the wings
of strain sensors at angles of 45 degrees. This setup however will depend on the
structural and geometric properties of the wing. The location of the bending and
torsion elastic axis will determine the exact location of these sensors.
In the case of a pure pitching motion, the total pitch moment on each wing
will be identical and ML = MR.
However in the case the aircraft is undergoing combined motions, this may
not apply. For instance, if the aircraft is undergoing a right turn while pitching:
ML > MR , and if it is undergoing a left turn ML < MR , due to the increased
loading on the inner turning wing.
The time propagation of the torque is described by the wave equation, therefore
the temporal displacements will follow harmonic oscillations. The time dependent
study however needs to take into account the appropriate initial conditions, which
may vary during flight.
In general the total pitching moment applied to the center of mass of the
70
Figure 3.4: Different pitching moments at root as a consequence of maneuver or ex-
ternal inputs
aircraft will then be:
MLorR =
n∑
i
GIpi
ri
∗ γi = τLorR (3.21)
In case of uniform cross section, composition and straight wings for pure pitch-
ing motion:
ML=R =
GIp
req
n∑
i
γi = Md (3.22)
Where Md is the pitch moment disturbance. An alternative approach to cal-
culating the difference in pitching moment of the aircraft by means of strain mea-
surements entails the use of bending strain sensors on the wings as well as on the
aircraft empennage horizontal stabilizer.
If we have that the vertical force on the wings in a one dimensional simplified
71
approximation, if the wings can be modeled by a beam with appropriate cross-
sectional characteristics:
Zwings = 2 ∗ Fs = 2 ∗
n∑
i=1
qi = 2 ∗
n∑
i=1
qi = 2 ∗
n∑
i=1
εi (x+ h)− 2εi (x) + εi (x− h)
h2
(3.23)
Where ε corresponds to the bending strains at sensor locations i spaced by
a distance h. In a similar way, the total load can be calculated on the horizontal
stabilizer of the aircraft:
Zhoriz = 2Fsh (3.24)
Such that the combined strain measurements on the wings and empennage
will yield a qualitative nose up or nose down pitch moment estimate, and will allow
for rapid difference calculations. After a certain threshold a compensator should
engage to reject the disturbance.
Zwings > Zhoriz, pitch up motion
Zwings < Zhoriz, pitch down motion
(3.25)
And it can also yield the total disturbance pitching moment:
Md = Zwingsdw − Zhorizdh (3.26)
72
Where dw is the distance between the wings effective center of force to the
aircraft center of mass and dh is the distance between the horizontal stabilizers
effective center of force to the aircraft center of mass. Ideally these will be co-planar
on the plane defined by the center of mass of the aircraft, which divides the aircraft
in symmetric halves along the nose axis.
3.1.3 Vertical Force
Similarly, based on Equation (3.23), a vertical force disturbance can be calcu-
lated as:
Zd = 2 ∗ Fs −mg (3.27)
Figure 3.5: Vertical force on aircraft and shear forces on the wings
The disturbance vertical force calculation may assume that the shearing force
on each wing is balanced for steady wings level flight Figure 3.5. However this can
73
be modified for a set tolerance for different maneuvers, i.e. steady turns, where the
shearing force on the left wing to the right may have a defined ratio.
The disturbance rejection mechanism based on these static strain comparisons
can be posed in different ways. One dimensional wide field integration constants
can be obtained for each sensor location and maneuver.
3.1.4 6-DOF model
Using results from the Roll, Pitch, and vertical roll formulations, a controller
which regulates the balance of forces by means of strain measurements can be for-
mulated in the fast inner loop discussed in Section (2.3.2) .
Ld = ML −MR = δ, Roll disturbance with bending strain
Md = τL − τR = δ, Pitch disturbance with shear strain
Md = Zwingsdw − Zhorizdh = δ, Pitch disturbance with bending strain
Zd = 2 ∗ Fs −mg = δ, Vertical force disturbance
(3.28)
Where δ is either zero or some calculated or measured value for steady wings
level flight or steady turns. A 6-DOF model can be obtained by obtaining the
remaining balancing equations, in rejection of larger values of δ. This is a somewhat
de-coupled method.
A coupled 6-DOF model can be obtained by computing least squares calcula-
tions on system identified matrices. It is proposed that the forces and moments at
the root of each wing need to balance out with respect to the center of gravity. By
74
parallel theorem this is only a constant away for rotational inertias of interest.
In this way, if a set of sensors has been system identified with respect to their
corresponding wing forces and moments, then a simple force balance of the forces
calculated by the regression matrices can be computed fast, if the matrices are
previously obtained by experimental load tests. In this sense a MIMO controller
can be formulated and built to regulate the following relationship:
Fˆrootleftwing − Fˆrootrightwing = ∆ (3.29)
A diagram of the method is illustrated in Figure 3.6.
Figure 3.6: 6-DOF coupled model for force/strain feedback
A formulation of the MIMO controller based on results from the theoretical
considerations is found in Figure 3.7 and Figure 3.8.
75
Figure 3.7: Formulation for the strain feedback MIMO controller based on SISO strain
balance
Figure 3.8: 6-DOF formulation for the strain feedback MIMO controller based on
system identified output matrices
76
Multiplicative uncertainty in the proposed output C matrices can be formu-
lated as shown in Figure 3.9. In this diagram the output matrices and the force
estimation matrices have been clustered into a single block. Techniques from robust
control can be then applied to synthesize the model based controller. A generalized
plant needs to be derived based on generalized inputs w, generalized outputs z. For
the uncertainty case a similar linear fractional transformation can be calculated us-
ing block diagram manipulations. Uncertainty in the Cmatrix will yield a family of
uncertain plants which will make the closed loop stable depending on the stability
margins found for the MIMO system. Weighting transfer function matrices can be
tuned for better closed loop response.
Figure 3.9: 6-DOF formulation for the strain feedback MIMO controller with uncer-
tainty in the Cmatrix
77
3.2 Strain based wide field integration (WFI)
3.2.1 1D strain WFI
Considering the static Euler Bernoulli beam models for the wing, which were
presented in Section (2.2.3), a one dimensional wide field integration can be formu-
lated based on the expected strains as a function of wingspan.
Alternatively, using the trigonometric substitution x = Lcosθ to get a simpler
form of the shearing force for an elliptic lift distribution from Equation (2.13):
Fs = ρ∞V∞Γ0
∫ (
1−
x2
L2
) 1
2
dx = −Lρ∞V∞Γ0
∫
(
sin2 θ
)
dθ (3.30)
Where L is the length of the wing and θ = arccos xL . Carrying out another
implicit integration, back substituting and expanding in Taylor series up to the third
degree:
Fs = −Lρ∞V∞Γ0
∫
(
sin2 θ
)
dθ = Lρ∞V∞Γ0 ·
(6L2x− 3piL3 + x3)3
648L9
(3.31)
Truncating this expression to the fourth term, the shear force becomes:
Fs = −Lρ∞V∞Γ0 ·
(
pi3
24
−
pi2x
4L
+
pix2
2L2
−
x3
3L3
)
(3.32)
78
Using this expression of the shearing force, the bending moment can be found.
The shearing force is the rate of change of bending moment:
Fs =
dM
dx
(3.33)
and the moment is:
M(x) =
∫
Fsdx (3.34)
Using the approximation found for the shear force, the moment is:
M(x) = −Lρ∞V∞Γ0
∫ (
pi3
24
−
pi2x
4L
+
pix2
2L2
−
x3
3L3
)
dx (3.35)
= Lρ∞V∞Γ0
(
x4
12L3
−
pi3x
24
+
pi2x2
8L
−
pix3
6L2
)
(3.36)
and the corresponding strain is calculated using:
εx = −z
M(x)
EI
(3.37)
using the approximation found for the moment, the strain is:
εx =
−zLρ∞V∞Γ0
EI
·
(
x4
12L3
−
pi3x
24
+
pi2x2
8L
−
pix3
6L2
)
(3.38)
79
which is the strain produced by an elliptical span wise aerodynamic lift. This
contribution of strain needs to be added to the contribution produced by the other
flight dynamics variables for a complete calculation.
εx = εxnominal + εxgusts (3.39)
where the nominal strain for an elliptic lift distribution corresponds to Equa-
tion (3.38):
εnominal =
−zLρ∞V∞Γ0
EI
·
(
x4
12L3
−
pi3x
24
+
pi2x2
8L
−
pix3
6L2
)
(3.40)
Assuming that the young’s modulus and second moment of area are constant
or can be reduced to a single equivalent value for the beam model of the wing. The
strain due to uniformly distributed gusts, is a summation of the strains produced
by
With the strain defined at each point of the cantilever, a set of coefficients
for each linearly independent term can be calculated when passed through an inner
product of orthogonal polynomials. Some of the possible polynomials which are
orthogonal in the -1 to 1 value range for the spanwise coordinate x are: Legen-
dre polynomials (from 0 to 1), Chebyshev polynomials, Laguerre polynomials. For
instance,the Laguerre polynomials:
L0 = 1, L1 = x, L2 =
1
2
(3x2−1), L3 =
1
2
(5x3−3x), L4 =
1
8
(35x4−30x2 + 3) (3.41)
80
where x goes from -1 to 1. Therefore, we have to normalize the span by
dividing by L. In this case, if we set y = x/L the normalized strain equation which
can be integrated in an inner product is:
εx = εnominal +
z
2EI
ρ∞V∞S
(
UgCL +
Uw
2
CLα +
Uv
2
CLα
)
(3.42)
Where the added expression is calculated under the assumption that the gusts
are impacting the aircraft uniformly and therefore the lift differential can be ex-
pressed according to Equation (2.1).
Lift differentials can be calculated based on a nominal expression for lift such
that any differences are due to gust perturbations. Coefficients can be calculated
via inner product due to the orthogonality of the polynomials in the range from
x=-1 to x=1. This approach can also be regarded as a one dimensional wide field
integration by means of basis functions. In a later chapter the WFI constants will
be derived based on the physics of the formulation.
Similarly, basis functions can be chosen based on the eigenstrains obtained
from the finite element calculations. Based on the type of gust affecting the struc-
ture, the patterns change according to mass participation factors.
3.2.2 2D strain WFI
Eigenstrains on the surface of the wing can be used as basis functions to gen-
erate a 2D eigenstrain pattern as well as weighting matrices. A general formulation
of strain field integration is shown in Figure 3.10.
81
Figure 3.10: Strain wide field integration
In particular, if the MPF weighted directional eigenstrains are used as basis
functions:
∫
εΦMPFv , Vgust
∫
εΦMPFu , Ugust
∫
εΦMPFw ,Wgust
(3.43)
and if the sensors are a discrete set, then they are appropriately weighted by:
∫
εdiscreteΦi (3.44)
A weighting function which would capture the three types of gusts would then
be given by:
∫
ΦMPFvΦMPFuΦMPFw = M (3.45)
82
3.3 Sensor placement
3.3.1 Sensor placement techniques
Optimal location of sensors can be a determining factor in both state estima-
tion and disturbance rejection. These methods are usually limited to energy and
modal analyses. In this work we not only explore the conventional approach to
sensor placement but also a wide field integrated approach [9] based on structural
inputs.
Sensor placement methods have been explored for a wide variety of purposes in-
cluding vibration suppression, actuator effectiveness, and active control [91]. These
methods can be deterministic or heuristic. Heuristic methods typically require a
cost function to be minimized or a performance function to be maximized. Deter-
ministic methods use more stringent conditions to find the optimal set of locations.
A mapping between the design variables and the physical sensor locations is required
for both. Some of these methods are appropriate for large undefined sets whereas
other methods require a small set of possible locations around local minima or max-
ima. For instance, some of the methods are based on neighbor search of locations
while others focus on selecting a location from widely scattered possible locations.
Genetic algorithms (GA) are effective at finding optimal solutions that are widely
scattered throughout the design space. The algorithm is initialized by generating
a large number of random subsets which are combined and mutated to search for
improved subsets.
83
Some of the most common methods will be described in what follows. Ex-
haustive search looks at all the possible sets of sensor locations to determine the
best; this can turn out to be computationally expensive. Local neighbor search or
Tabu Search (TS) uses less computational resources by keeping memory structures
that describe the visited solutions of a provided set of rules. Another method which
uses pre-selected locations is the Neural Network approach where sensor location are
discarded based on each sensor contribution to the optimization criteria. Similarly,
Simulated Annealing entails a nearest neighbor search based on optimization of
subsets which are selected upon improvement of the cost function. Energy methods
include: maximization of the modal kinetic energy (MKE) of a sensor configuration,
selection of a set with the highest entropy and maximization of the modal strain
energy (MSE) [92].
Modal techniques for sensor placement are yet another subcategory of tech-
niques. Driving point residue method (DPR) [93] [94] is a modal participation factor
methods where the sensors are placed based on the highest average response over
all of the target modes. The eigenvalue product method (EVP) forms an indicator
consisting in the absolute value of the product of every eigenvector. The best loca-
tions are those with the highest product value. Singular value decomposition (SVD)
of a given candidate set of sensor locations can also be used to find an optimal set.
Effective independence (EFI) [95] is a heuristic ranking scheme that attempts
to maximize the determinant of the Fisher information matrix (FIM) [96] or equiv-
alently, maximize its minimum eigenvalue. The FIM is directly correlated to the
error covariance, which gets minimized in the process. The Fisher information ma-
84
trix is also the inverse of the observability Grammian of sensor locations [97]. The
observability Grammian is used as an assessment of the contribution of the individ-
ual modes. Each contribution should be high in order to warrant high observability
in a specific set of sensor locations. In this work we will be focusing on three ap-
proaches: mass participation factor weighting, effective independence /observability
Grammian and wide field integrated strains.
As a first approach to identifying the best sensor solutions using the EFI
methodology, we set up an Aluminum 1D beam with free-free boundary conditions
and carried an eigenfrequency analysis. The first 6 modeshapes that resulted from
the finite element simulation were extracted. Figure 3.11 shows the results of the
modal displacements of such beam and the corresponding eigenstrains. These are
symmetric due to the boundary conditions. The same technique applies to cantilever
beams and wings. The results of the Comsol structural mechanics eigenvalue solver
yielded the numerical solutions shown in Figure 3.11. This technique will be used
in following analyses of 1D and 2D wings.
The optimization of sensor locations techniques can be applied to any struc-
ture, and in particular to the aircraft wings. By customizing the cost function to be
a performance metric, one can tailor the locations to best observe the quantities of
interest, i.e. Forces and Moments, gust forces, etc... The performance metrics may
also include observability and controllability metrics.
85
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
truss span (m)
t
o
t
a
l 
d
is
p
la
c
e
m
e
n
t
 
(
m
)
 
 
440.04 Hz
880.16 Hz
1320.79 Hz
1762.96 Hz
2208.43 Hz
2659.79 Hz
(a) First 6 modal displacements of free-
free 1D beam (truss)
0
1
2
3
4
5
6
7
8
9
0
0.5
1
1.5
2
2.5
-5
-4
-3
-2
-1
0
1
2
3
4
5
 
6 modal displacements
beam span (m)
 mode shape displacements
 
2.62 Hz
16.38 Hz
20.36 Hz
24.66 Hz
45.88 Hz
62.93 Hz
(b) First 6 modal displacements of can-
tilevered 2D beam)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
0
500
1000
1500
2000
2500
3000
3500
4000
4500
truss span (m)
m
o
d
a
l 
s
t
r
a
in
Truss modal strain -6 modes
 
 
440.04 Hz
880.16 Hz
1320.79 Hz
1762.96 Hz
2208.43 Hz
2659.79 Hz
(c) First 6 modal strains of free-free 1D
beam (truss)
0
1
2
3
4
5
6
7
8
9
0
0.5
1
1.5
2
2.5
-5
-4
-3
-2
-1
0
1
2
3
4
5
 
beam span (m)
6 modal strains
 strain mode shapes
 
2.62 Hz
16.38 Hz
20.36 Hz
24.66 Hz
45.88 Hz
62.93 Hz
(d) First 6 modal strains of cantilevered
2D beam
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
11000
truss span
s
u
p
e
r
im
p
o
s
e
d
 
m
o
d
a
l 
s
t
r
a
in
Superimposed modal strain
(e) Eigenstrain superposition
0
2
4
6
8
10
0
0.5
1
1.5
2
2.5
-10
-5
0
5
10
15
beam span (m)
strain mode shape superposition
 strain mode shape superposition 
(f) Eigenstrain superposition
Figure 3.11: 1D dimensional aluminum beam eigenstrain superposition (left) and 2D
aluminum beam eigenstrain superposition (right).
86
3.3.2 Theoretical Framework EFI and Grammian
A variety of methods can be applied; however, we have chosen to combine
the most representative ones based on the nature of our application. A wide field
integrated approach will also be pursued; it includes physics that are not normally
considered for sensor placement. In this sense, the conventional approaches can
define the local subsets of maxima based on energy and structural considerations.
Then a physical approach can determine sensor output weightings for an estimated
maximal output. Sensor locations can be first identified by using mass participation
factors. These coefficients provide a method for judging the significance of a vibra-
tion mode. That is, modes with higher factors can be readily excited whereas ones
with lower factors cannot. This is an appropriate analysis for an aircraft wing in
that it undergoes different types of loads during flight which may produce deforma-
tions and strains. It allows a study of the more significant impacts of the loads on
the wing, which produce strains. There is a participation factor for every available
degree of freedom.
In addition to mass and energy considerations, an effective independence ap-
proach is useful in that it allows for the error covariance to be minimized, or equiv-
alently for the observability Grammian to be maximized over a range of possible
sensor locations. It is also useful in the sense that a modified cost function can
be applied. This function entails not only the modal participation factors but can
also contain information about the flight dynamics in a closed loop strain feedback
87
configuration. In EFI, the output equation can be expressed as:
ys = Φfsq + w (3.46)
In this equation, ys is a vector of sensor outputs at a large candidate set of sen-
sor locations. Φfs is a matrix of FEM target modes which have been Φf partitioned
to the corresponding sensor location degrees of freedom. This partitioning refers to
the removal of the targeted sensor location to evaluate its contribution. The target
mode response vector is q in normal coordinates. The noise is represented by w,
and it is assumed to be a stationary random observation disturbance possessing zero
mean and covariance intensity matrix satisfying the expectation operator. Estimate
error covariance matrix for q is found to be the inverse of the Fisher information
matrix Q [95] and equivalent to the observability Grammian Yo :
E[(q − qˆ)(q − qˆ)T ] = [(Φfs)
TRΦfs]
−1 (3.47)
Maximizing determinant of the FIM will maximize the spatial independence
of the target modal partitions as desired and maximize the signal strength of the
target modal responses in the sensor output. This results in the best estimate of the
modal response. Sensors should be placed such that Q is maximized or equivalently
such that Yo is maximized. The Effective Independence Distribution vector can be
88
used as the cost function.
ED = [Φ¯fsΨ]
2(Λ)−11k (3.48)
Where Ψ is the matrix of eigenvectors of Q and Λ is the corresponding matrix
of eigenvalues. 1k is a k-dimensional column vector of ones .Each term within vector
ED represents the contribution of the corresponding candidate sensor location to
the rank of the information matrix and thus the linear independence of the target
modes. Another equivalent calculation of the effective independence vector [98] is
found from:
ED = diag[Φ¯fs(Φfs
TΦfs)
−1Φfs
T ] (3.49)
A ranking of candidate sensor locations can be calculated by:
EDi =
det(Q)− det(QTi)
det(Q)
(3.50)
and 0 < E(Di) < 1.0 is the range of possible values of the ranking of candidate
sensor locations. A zero value indicates that the corresponding candidate sensor
location will contribute nothing to the identification of the target modes. Additional
cost functions can be implemented in a closed loop entailing the entire set of sensor
measurements, weighted by wide field integrated constants.
89
3.4 1D Sensor Placement Analysis
3.4.1 Wing FEM model
A typical UAS plane was used for the aeroelastic simulations and sensor place-
ment analysis. The wing types used in these one dimensional linear structural models
is a unidirectional carbon fiber with typical Young’s modulus and material charac-
teristics. The wing was modeled using Comsol Multiphysics Figure 3.12 with slender
beam elements and had an aspect ratio of approx. 7. The mesh consisted of nearly
1000 elements in two dimensions. An eigenfrequency analysis was performed with-
out considering geometric nonlinearities. A parallel sparse direct solver was used
to solve the system of equations. The wing root was modeled as a fixed constraint
whereas the wing tip was modeled as a free constraint. The eigenvalue analysis
shows that the first six bending eigenfrequencies are: 0.69 Hz, 4.37 Hz, 12.44 Hz,
24.95 Hz, 42.5 Hz, and 65.95 Hz. These frequencies were obtained with the assump-
tion of a uniform cross-section and may vary depending on the internal structure
refinements such as ribs, spars, etc.
The following sensor placement analysis will be based on the modal results
rendered by the finite element model. Two approaches will be presented. The first
approach considers a modal energy approach complimented by the gradient, and the
second approach entails a measure of the information matrix as provided by each
sensor location. This second approach is equivalent to computing the observability
matrix, which is also a measure of the effective independence of the targeted modal
90
strains.
Figure 3.12: First and second finite element eigenmodes of UAS wing
3.4.2 Sensor placement using MPF/gradient and EFI/Grammian
Each eigenmode has an associated eigenstrain shape. From Figure 3.12 the
associated stresses and strains can be seen around the contour of the eigenmode.
The eigenstrains were normalized with respect to wing length Lw and with respect
to maximum strain, as shown in Figure 3.13. They directly correspond to the
eigenfrequencies in the same corresponding order. A complete set of finite element
eigenfrequencies and mass participation factors is shown in table 3.1. Participation
factors in the normal direction are called MPFv and in the axial direction are named
MPFu. In particular we are interested in the influence of structural and energy
characteristics in the vertical or normal direction since this is the direction most
impacted by gusts. The different magnitude scales of mass participation factors
can be seen from 3.1 which shows that lateral loading doesn’t have much of an
91
impact on bending modes and therefore chordwise bending moment which is directly
proportional to strain.
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1 first mode
normalized span
nor
ma
lized
 stra
in
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1 second mode
normalized span
nor
ma
lized
 stra
in
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1 third mode
normalized span
nor
ma
lized
 stra
in
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1 fourth mode
normalized span
nor
ma
lized
 stra
in
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1 fifth mode
normalized span
nor
ma
lized
 stra
in
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1 sixth mode
normalized span
nor
ma
lized
 stra
in
Figure 3.13: UAS normalized eigenstrains from finite element simulation
Table 3.1: Mass participation factor results for 2D model of wing
Eigenfrequency (Hz) MPF u MPF v Total Strain Energy (J)
0.69 -7.35E-13 1.35 4.72
4.37 -4.30E-12 0.74 189.17
12.44 -9.17E-12 0.43 1529.18
24.95 -1.60E-11 0.30 6144.57
42.50 -2.30E-11 0.23 17828.18
65.95 6.15E-10 0.18 42927.74
92
The mass participation factors in the vertical or normal direction are used to
weigh each of the six modes. This will give a relative indication of how well a partic-
ular mode is excited from a given degree of freedom. A comparison between linearly
weighted eigenstrains and MPF weighted strains can be seen from Figure 3.14. It
is noticeable that the mass participation factors dramatically modify the expected
areas of maximal sensor outputs, or areas where the absolute value of strain is ex-
pected to be larger. A comparison of the calculated gradients for both the linearly
weighted and the MPF weighted strain distributions can be seen from Figure 3.15.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
1
2
3
4
5
6 Linearly weighted and MPF weighted strain
normalized span
stra
in g
rad
ien
t
 
 
linearly weighted
MPF weighted
Figure 3.14: UAS normalized weighted eigenstrains using equal weights and MPF
The MPF weighted eigenstrain curve can be combined with its corresponding
gradient to avoid locations of fast strain rate. This normally depends on the sensor
93
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−70
−60
−50
−40
−30
−20
−10
0
10
20 Strain gradients: linearly weighted and MPF weighted
normalized span
stra
in g
radi
ent
 
 Linearly weighted gradientMPF weighted gradient
Figure 3.15: Linearly weighted and MPF eigenstrain gradients
location targeted application but typically areas of fastly changing strain need to
be avoided due to increased risk of measurement nonlinearities. These areas of
sharply changing strain can induce noise and increase possible debonding, in the
case of conventional strain gauges. The MPF weighted curve and its corresponding
gradient are shown in figure 5. Sensor placement that results from this analysis
is depicted with shaded columns. A stronger colored column indicates a higher
expected overall strain value whereas a fading one indicates the strain value is not
expected to be as high. The columns intersect areas of higher strain and constantly
increasing gradients, avoiding gradient peaks.
A similar analysis was performed using EFI/Grammian analysis based on the
finite element results for eigenshapes. The fisher information matrix was calcu-
94
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−60
−50
−40
−30
−20
−10
0
10
normalized span
MP
F e
ige
nsh
ape
MPF Weighted Eigenshape and Gradient
 
 
MPF weighted eigenshapeGradient
Figure 3.16: UAS wing MPF based sensor placement using weighted eigenstrains and
gradient
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-60
-50
-40
-30
-20
-10
0
10
normalized span
MPF 
e
i
ge
n
s
ha
p
e
MPF Weighted Eigenshape and Gradient
 
 
MP F weighted eigenshape
Gradient
Figure 3.17: UAS wing MPF based sensor placement using weighted eigenstrains and
gradient
95
lated using six eigenstrains. Then the effective independence distribution vector
was calculated by obtaining the diagonal of the matrix in (3.49) . Each term in the
vector represents the contribution of the corresponding candidate sensor locations
to the rank of the Fisher matrix and therefore the linear independence of the target
eigenstrains. The resulting independence vector ED was calculated as a function of
normalized span, as shown in Figure 3.18 .
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
normalized span
E D
EFI sensor placement
 
 EFI contributions
Figure 3.18: UAS wing EFI/Grammian sensor placement
The first ten locations with higher contributions to the observability Gram-
mian are found in table 3.2. These values were picked such that they would corre-
spond to rather distributed areas as opposed to being in close proximity. The table
is arranged from higher to lower values of the effective independence vector.
A comparison of results obtained between the energy method of mass partici-
96
Table 3.2: Distributed sensor locations with higher EFI contribution
Normalized span EFI contribution
0.4935 0.1532
0.2208 0.1236
0.06 0.1088
0.72 0.1065
0.2727 0.1035
0.1299 0.1032
0.05 0.0992
0.6494 0.0966
0.1429 0.0965
0.2857 0.0931
pation factors with gradient considerations and the effective independence/ observ-
ability Grammian method are presented in Figure 3.19. It can be observed from
this figure that the observability metric allows for more potential sensor locations
while the energy weighted modal strain approach provides an emphasis of the areas
of higher energy for the chosen direction of the participation factor. In this case,
this MPF corresponds to the direction normal to the span/chord wing plane, which
is the direction in which gusts would have a greater impact. Other MPF factors
can be calculated and applied to the mass weighting method based on the intended
sensing application.
97
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
normalized span
E D
 
an
d M
PF
/Gr
adi
ent
EFI and MPF sensor placement comparison
 
 
EFI
MPF/gradient
Figure 3.19: Comparison of EFI/Grammian and MPF/Gradient sensor placement
methods
98
3.5 Gust simulations
3.5.1 ASWING aeroelastic model simulation
An aeroelastic model of the UAS was implemented in ASWING [72] to evaluate
the strain wide field integration approach presented in section 3.3.2 as well as to look
for the impact of sensor placement considerations from section 3.4. ASWING allows
for fully non-linear aeroservoelastic simulations where lifting line theory is used for
aerodynamic computations and is coupled with non-linear Euler Bernoulli beam
elements for structural computations (Figure 3.20). Control surfaces can be defined
within the model geometry as well as full 3D gust fields.
The aircraft model used in the simulations is consistent with the wing beam
model used for finite element analisys of section 3.43.4.1. It consisted of four beams,
fuselage and wing planform, vertical and horizontal stabilizer, as can be seen in
Figure 3.21. Each beam had structural properties as well as aerodynamic properties
such as lift and drag coefficients. The propeller can be modeled as a point mass
with thrust vector properties. Sensors can be placed at different locations on the
wing to measure both aerodynamic (i.e. angle of attack) and structural properties
(i.e. strains). A large number of coupled aeroelastic modes as well as typical rigid
body dynamics can be obtained from the fully non-linear simulation. These can be
reduced if the system is linearized about a trim condition. A reduced order model
was obtained about a specific flight condition, namely steady wings level flight at
a constant speed of approximately 15 m/s. This linear time invariant model allows
99
for extraction of rigid body states as well as structural dynamics states associated
to bending of the wing. A coupled dynamics matrix was obtained and used for
extraction of the strain outputs at particular locations on the wing. Sensors can only
be placed at 12 locations at a time in ASWING, therefore a different geometry needs
to be executed for every change in sensor location. The linear reduced order model
implemented around the steady wings reference flight condition had 6 displacement
sensors on each wing, for a total of 12 sensors. The strains were also calculated
along the span of the wing.
Figure 3.20: ASWING definitions (excerpt from reference [72])
3.5.1.1 ASWING nonlinear simulation -3D gust
A heuristic non-linear simulation was implemented in ASWING using the UAS
model described in the previous section. The aircraft was set to interact with a 3-D
gust field of two different gust intensities. The simulation was set such that the
100
propeller 
Fuselage 
Horizontal  
Stabilizer Vertical  
Stabilizer 
Wing sensors 
Figure 3.21: Top view of aircraft model for ASWING simulation
aircraft would follow a straight path and then encounter a vertical gust field. This
produced a disturbance of the rolling angle and other rigid body states, as can be
observed from the time histories of the bank angle. Initial gust simulation results
are shown in Figure 3.22, where the is a measurable impact of gust vertical velocity
on the vehicle’s flight states .
3.5.1.2 ASWING linear simulation -1D gust
In the reduced order model, a similar roll gust effect was obtained by injecting
a disturbance in the left most aileron. The impact of this disturbance on the linear
model can be seen in Figure 3.24. This artificial roll disturbance has a similar effect
on the roll angle as that observed in the fully non-linear case, where the aircraft was
101
Figure 3.22: Airplane vertical velocity and bank angle with and without gust, nonlin-
ear ASWING simulation
set to fly over a simulated gust field.
Strain wide field integration was implemented in closed loop feedback by ob-
taining the strains from each sensor location (Figure 3.21) and weighting them using
the output Cw matrix. The calculated actuator command was then fed back to the
symmetric aileron on the right wing, to account for the effect of the roll distur-
bance generated by the aileron on the left wing. The strain feedback loop however
was only implemented for roll disturbance rejection in the way of an inner loop
(Figure 3.23). There were no additional outer loop controllers in place for further
attitude stabilization, therefore the simulations only evaluated the performance of
strain feedback alone. The simulated disturbance consisted on rapidly changing the
left aileron angle by a few degrees for one second after one second of flight. The
effect of the simulated roll disturbance on the vehicle states can be observed from
Figure 3.24. In this figure, it can be observed from the roll angle time history that
102
the roll disturbance causes the roll angle to increase over time. Other states also
become unstable as a consequence of the disturbance. The roll disturbance has an
impact on the wing strains, this can be observed from Figure 3.25. This information
is used in the strain feedback loop.
A combined set of sensor locations from EFI and MPF calculations were used
and the impact on the closed loop in the presence of a roll disturbance was ob-
served. The Cw matrix was calculated using these locations and used in closed loop
to weigh the values of strain coming from each sensor. These EFI/MPF sensor loca-
tions were then compared to the performance of random sensor locations to mitigate
the roll disturbance. The random set produced instabilities in the rigid body states,
therefore the EFI/MPF sensor locations were conserved. A PID controller was im-
plemented in the strain feedback loop using the EFI/MPF locations and at different
gains, higher gains enable a more effective strain feedback control, as can be seen
from Figure 3.26.
Inner loop stabilization by means of bioinspired mechanisms was demonstrated
for a one degree of freedom roll disturbance case. Distributed strain sensing showed
to be effective in counteracting roll disturbances. The weighting patterns in the
output matrix are crucial to the bioinspired reactive mechanism which addresses the
disturbance at higher order states before they propagate onto lower order states. A
physics based weighting pattern was implemented and it showed to be useful in roll
disturbance mitigation. This pattern was obtained by means of force and moment
considerations.
Sensor placement is also a fundamental part of distributed sensing arrays.
103
G(s)
K
Aircraft  
Plant 
Controller 
+ L + +  Aero 
Strain Sensors 
- 
Wing wC
Spatial 
Weights 
rL distL
Figure 3.23: Strain feedback used in ASWING reduced order simulation
0 2 4 6 8 1014
15
16
17
Time [s]
U [m/
s]
0 2 4 6 8 10−0.5
0
0.5
1
Time [s]
V [m/
s]
0 2 4 6 8 100.2
0.3
0.4
0.5
Time [s]
W [m
/s]
0 2 4 6 8 100
10
20
30
40
Time [s]
φ [de
g]
0 2 4 6 8 10−5
0
5
Time [s]
θ [de
g]
0 2 4 6 8 100
200
400
600
Time [s]
ψ [de
g]
0 2 4 6 8 100
10
20
30
Time [s]
P [de
g/s]
0 2 4 6 8 10−2
0
2
4
Time [s]
Q [de
g/s]
0 2 4 6 8 10−20
0
20
40
Time [s]
R [de
g/s]
Figure 3.24: Roll disturbance on aircraft states
104
0 500 1000 1500
2000 25000
2
4
6
8
10
12
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1 x 10
−3 scaled displacement gradients at sensor locations for flight duration
time(samples)=10 seconds
Sensor locations along wingspan
sca
led 
stra
ins
Figure 3.25: Impact of roll disturbance on sensors
0 1 2 3 4 5 6 7 8 9 100
5
10
15
20
25
30
35
40
Time [s]
φ [de
g]
Impact of strain feedback on roll
 
 
No controlK=1K=5K=10K=20
Figure 3.26: Roll angle at different strain feedback proportional grains K, with deriva-
tive D=1 and integral I=2 gains.
105
Two methods were explored, one based on system observability (EFI/Grammian)
approach and the other one based on modal energy superposition with gust direc-
tionality (MPF-v). A combined approach using these two methods yielded sensor
locations which were used in conjunction with the calculated output matrix. These
provided disturbance mitigation results which were compared to those of a random
sensor location set. It was observed that the randomly placed sensors did not pro-
duce a reduction in the roll disturbance. It was also observed that the optimal sensor
locations in the closed loop case may have been close to the ones calculated in the
open loop case. Symmetric sensor locations seemed to have consistently shown some
improvement in roll disturbance rejection. A PID controller was implemented in the
strain feedback loop. It was observed that roll disturbance rejection increased at
higher proportional gains and relatively small derivative and integral terms.
Robust controllers based on strain feedback can be implemented for higher
degree of freedom systems by means of higher order reactive loops. A closed loop
sensor location optimization can be used to compare to locations predicted by energy
and observability calculations.
3.5.2 Fluid structure interaction simulations- Static Aeroelastic (Com-
sol)
A static aeroelastic simulation was built using the FSI or fluid structure inter-
action module in Comsol Multiphysics. The module combines fluid flow with solid
mechanics to capture the interaction between the fluid and the solid structure. The
106
fluid structure interaction couplings appear on the boundaries between the fluid
and the solid. The fluid structure interaction interface uses an arbitrary Lagrangian
Eulerian or ALE method to combine the fluid flow formulated using an Eulerian
description and a spatial frame with solid mechanics which has been formulated
using a Lagrangian description in a material reference frame.
The fluid flow is described by the Navier Stokes equations which yield a solu-
tion for the fluid velocity field. The fluid exerts a force on the boundary of the solid,
which is given by the negative of the reaction force on the fluid [99]. This reaction
force is given by:
f = n.
(
−pI +
(
µ(∇u+ (∇u)T )−
2
3
µ(∇u)I
))
(3.51)
where p is the fluid pressure, n is the outward normal to the boundary, and
I the identity matrix. Since the Navier stokes equations are solved in the spatial
and therefore deformed frame and the solid mechanics interfaces are defined in the
undeformed frame, a transformation of the fluid force on the solid is necessary. This
transformation scales the mesh elements in both frames to come up with a mesh
ratio. This ratio is defined by the spatial element factor dv and the material element
dV . f is the force on the fluid, defined in the spatial frame and F is the force defined
in the material frame.
F = f
dv
dV
(3.52)
107
The solid mechanics counterpart entails the structural velocities to act as a
moving wall for the fluid domain. These bidirectional couplings allow for the fluid
structure interaction simulation. The mesh is free to move inside the fluid domains
and adjusts to the motion of the solid walls. The geometric change of the fluid
domain is automatically accounted for in Comsol.
The static aeroelastic simulations were set in the way of a wind tunnel test
scenario Figure 3.27. This approach therefore consists of a simulated wind tunnel
environment. There is an inlet upstream of the leading edge of the wing and an outlet
downstream of the trailing edge of the wing and four non-slip surfaces. This scenario
was chosen due to the computationally expensive nature of CFD computations and
due to the purpose of the simulations which is to observe the interactions of the
fluid and the wing, and the produced strains and stresses on the surface of the wing.
Figure 3.27: Wind tunnel simulation of fixed wing with impinging vertical gust
108
An acrylic wing was used for these simulations. The mesh elements used in
the simulation were tetrahedral elements in the order of 450000, as well as triangular
surface elements in the order of 30000 and about 9000 edge elements. The length of
the box was made shorter than typical virtual tunnel CFD simulations due to the
computational requirements, as it allowed to observe the interactions of the wing
material and the fluid.
A vertical gust was included in the simulation with the purpose of observing the
effects on the wing strains. Results from these simulations are shown in Figure 3.28.
The speed of the tunnel was varied from 0 to 30 m/s and the gust velocities varied
from 0 to 15 m/s, which are a bit unrealistic but were mainly for heuristic purposes.
Correlations of gust velocities, wind speeds and surface strain maps can be drawn
from different simulation cases.
Figure 3.28: Wind tunnel simulation of fixed wing with impinging vertical gust, strains
are produced as a result of the impinging vertical gust
Another interesting simulation setup which was implemented consists of a
109
rotational gust. This gust was simulated by placing gusts on opposite sides of the
wing. The fluid pressure profiles can be observed from Figure 3.29. The resulting
strain profile is quite different from the one observed for an impinging unidirectional
gust from Figure 3.28.
Figure 3.29: Wind tunnel simulation of fixed wing with impinging rotational gust:
Fluid pressure profiles (left) and wing strains (right)
From these simulations another important result was observed. The principal
strains and stresses are a function of the gust directions. These quantities vary de-
pending on the direction in which the fluid is interacting with the wing structure.
Further simulations will be conducted in other works to study and model this in-
teraction in depth. These directions can be observed from Figure 3.30. In these
figures, the stress surface is superimposed and the arrows show the first principal
stress in the wing due to the gust.
Time dependent simulations can also be undertaken to study dynamic aeroe-
lasticity effects, however, these are extremely computationally expensive with the
110
Figure 3.30: Principal strains for different impinging gusts: vertical gust (left) and
torsional gust (right)
current computational platform and will be left for further study. A single one
second simulation took approximately 8 days to successfully run. Another impor-
tant study that will be left for further study is the FSI validation. Comsol has
several CFD validation models and structural models, however these have not been
validated jointly in the FSI module. These need to be verified as part of the well
understood CFD and simulation requirements [100], [101].
3.6 2D sensor placement analyses (Comsol)
3.6.1 PWing test section
3.6.1.1 Wing FE model and sensor placement analysis
A wing section was used to test the concepts from the energy based sensor
placement approach. This approach was shown to have similar results to those
provided by the Grammian analysis for a vertical gust. The technique is now applied
111
to the surface of a three dimensional wing to find the best sensor locations to give
best outputs amidst vertical gusts.
The software used for the finite element simulations was COMSOL Multi-
physics. The wing was an unswept wing with airfoil profile AG-35 manufactured by
Borealis Composites [102]. It was constructed such that the torsional compliance
is much larger than the bending compliance. The layup configuration consisted of
foam core, 1k carbon fiber weave placed at +- 45 degrees with respect to the length-
wise axis and one layer of +-90 degree fiberglass cloth as the outtermost layer. The
wing geometric characteristics are shown in Figure 3.31 and Table 3.5. The aero-
dynamic properites are calculated based on the approximate characteristics of the
NACA 2408 airfoil.
Figure 3.31: Airfoil geometric characteristics
The wing also had two ribs, one at the root and one at the tip made out of
wood. These however are not carrying any significant loads. There are no ribs or
spars, therefore the entirety of the loads is carried by the carbon fiber shell. For a
simplified model, a carbon fiber at 0 degrees shell with extruded styrofoam core was
considered. The behavior of the composite is mostly unidirectional, therefore this
112
Table 3.3: Wing geometric characteristics
Wing characteristics Value Wing characteristics Value
Chord length 22.98 cm Approx. NACA 2408
Quarter chord 5.75 cm Lower flatness1 93.0%
Max. thickness 2 cm Max. CL1 0.953
Max. Camber 4 mm Max.CL angle1 13.5
Wing length 62 cm Max. L/D 1 41.453
Shell thickness 0.0528 cm Stall angle1 2.5
Second moment of area 5.9748e-08 mˆ4 Zero Lift angle1 -2.0
Dihedral 0 deg. Ailerons 2
approximation was made.
Figure 3.32: Normal element size mesh of carbon layup wing
113
The solvers used to solve for the eigenmodes is the MUMPS (or Multifrontal
Massively Parallel sparse direct Solver) [103]. This solver is a variant of Gauss elim-
ination which builds a LU or Cholesky matrix decomposition of the sparse matrix
and then solves it one subset, or front, at a time. The multifrontal approach refers
to the usage of parallel cores to process different frontals simultaneously. This solver
uses several preordering algorithms to permute the columns for prearrangement.
The eigenmodes obtained are in-vacuo and with cantilevered boundary condi-
tions. The finite element results are shown in Figure 3.33. The frequencies of these
eigenstresses are, respectively, 53.682 Hz, 166.65 Hz, 184.66 Hz, 235.80 Hz, 286.05
Hz, 325.02 Hz.
A top view of the 3D wings gives an idea of the 2D eigenstresses in the same
modes. This is shown in Figure 3.34.
Using these first six modes, we proceed now to weighting them and superim-
posing them by means of mass participation factor coefficients. These coefficients
provide a method for judging the significance of a vibration mode. That is, modes
with higher factors can be readily excited whereas ones with lower factors cannot.
This is an appropriate analysis for an aircraft wing in that it undergoes different
types of loads during flight which may produce deformations and strains. This al-
lows a study of the more significant impacts of the loads on the wing, which produce
strains.
The following mass participation factors were also calculated in the finite ele-
ment solver and are used for weighting the modes and combine them into a linearly
superimposed eigenstrain profile. The resulting profile is used for the sensor place-
114
(a) First eigenstress mode (b) Second eigenstress mode
(c) Third eigenstress mode (d) Fourth eigenstress mode
(e) Fifth eigenstress mode (f) Sixth eigenstress mode
Figure 3.33: First six eigenstresses of carbon fiber wing
115
(a) First eigenstress mode (b) Second eigenstress mode
(c) Third eigenstress mode (d) Fourth eigenstress mode
(e) Fifth eigenstress mode (f) Sixth eigenstress mode
Figure 3.34: 2D view of first six eigenstresses of carbon fiber wing
116
ment study. The resulting mass participation factors are shown in Table 3.6.
Table 3.4: Mass participation factors of carbon fiber wing
Eigenfreq. MPF-u MPF-v MPF-w
53.68259 -0.00297 0.34194 -7.62E-05
166.6564 -0.011898 0.07544 -5.04E-05
184.6618 0.00843 0.11441 1.95E-04
235.8009 -0.00412 0.05275 -1.30E-04
286.0515 -0.0026 0.03837 1.07E-05
325.0233 0.00675 -0.0038 2.35E-04
The eigenstrain profiles are calculated for two cases of linear superposition
of modes. One case with equal weighting coefficients and the other one with the
mass participation factor MPF weighted eigenstrains in the v-direction, which for
purposes of the simulation corresponds to an out of plane displacement, or displace-
ments which are produced by loads that are normal to the wing, i.e. lift. The
resulting eigenstrains are shown in Figure 3.35.
In this figure, the root of the wing is at zero wing length and the trailing edge
is at the zero chord line. In this figure, the areas of maximal stress are in red and the
areas of minimal stress are in dark blue. The equal weighting mode superposition
gives a very different strain field from that of the MPF-v weighted eigenstrains. It
can be observed from the graphs that the areas of maximal stress are rather close
to the trailing edge, this is probably due to a higher bending moment in thinner
sections of the wing, as the moment of area goes down in these regions. The areas of
117
Figure 3.35: Superposition of first six eigenstress modes (top view)
118
maximal strain which are suitable for placing strain gauges are those which have not
so high gradients. This combination will yield better results at the implementation
phase because the areas of larger strain gradients can cause the strain gauges to
detach and or produce rapidly changing signals.
The possible number of sensors and sensor arrangements depends on the avail-
able sensor platforms. These platforms will be presented in the following chapter.
One of the sensor platforms that were used was provided by the structures and com-
posites laboratory at Stanford University [8]. This platform consists of a stretchable
sensor network, where the strain gauges and interconnections are microfabricated
into a single stretchable mesh. The constraint of this network entails a maximum
of 15 sensors with 5 x 9 possible sensor locations in a grid which has an adjustable
length. A maximum of three sensors per lengthwise row could be chosen at a time.
Possible scenarios and distributions based on these constraints and the analysis
shown in Figure 3.35 are shown in Figure 3.36.
Figure 3.36: Strain gauge stretchable network in spanwise orientation (Structures and
Composites Lab SACL, Stanford University) [8]
119
Figure 3.37: Sensor distributions based on SACL lab stretchable sensor network
120
The three scenarios presented in Figure 3.37 were chosen based on the intended
measurements. The first scenario is a somewhat uniform vertical distance modal
capture distribution, as it covers most of the areas of maximal stress and lower
gradients. The second scenario also intends to capture most of the modal strain,
however with an adjusted vertical distance to look for more convenient areas for
sensor locations. The fixed vertical distance scenario is a distribution that may also
allow for possibly capturing other phenomena or serving to validate assumptions
during the modeling process.
Another important analysis feature is given by the sensor orientations. In
nature, as it was seen in the bio inspired section of this work, strain sensors of ani-
mals are oriented depending on what their intended functionality is. Campaniform
sensilla are elliptical and the major axes are oriented in the maximal strain output
directions. In strain analysis, the orientation of the sensors is given by the principal
stress directions. These are directions which are calculated to be the directions in
which strain is varying the most, according to geometry and full 3D material prop-
erties. These orientations are also calculated in Comsol, and can also be weighted to
find energy weighted patterns. These strain fields can also be superimposed accord-
ing to the mass participation factors. However this procedure entails 3D vectorial
summation to generate the superimposed vector field of principal stress/strain. Fig-
ure 3.38 are the principal directions for the first three modes of the wing.
121
1st  
2nd  
3rd  
4th  
5th  
6th  
L E  
TE  
Roo
t
 
Figure 3.38: First six principal strain directions of carbon fiber wing
122
3.6.2 Flight test PWing
A similar approach was followed for the actual test flight wing, which is to be
assembled with the Apprentice aircraft airplane.
The full wing for flight tests with the apprentice had a similar layup as com-
pared to the first Pwings test section. It has dihedral and a similar carbon fiber and
fiber glass layup. The wing characteristics are described as shown in the following
table. The airfoil is the same profile AG-35 which can be approximately described
by the NACA 2408. Figure 3.39 shows the geometric characteristics as designed by
the manufacturer [102] and by Aurora Flight Sciences [7]. Figure 3.40 shows a set of
two manufactured PWings without strain gauge instrumentation or servo motors.
Figure 3.39: Pwing Airfoil and full wingspan geometry
123
Table 3.5: Wing geometric characteristics
Wing characteristics Value Wing characteristics Value
Chord length 9 in Approx. NACA 2408
Quarter chord 2.25 in Lower flatness1 93.0%
Max. thickness 2 cm Max. CL1 0.953
Max. Camber 4 mm Max.CL angle1 13.5
Wing length 62 cm Max. L/D 1 41.453
Shell thickness 0.0528 cm Stall angle1 2.5
Second moment of area 5.9748e-08 mˆ4 Zero Lift angle1 -2.0
Dihedral 6 deg. Ailerons 4
Ribs 2 deg. Small joining spar 1
Figure 3.40: Flight Pwings for assembly with Apprentice Aircraft
124
(a) Comsol model (b) FE mesh
Figure 3.41: FE model of flight test wing and meshed geometry
3.6.2.1 Wing FE model and sensor placement analysis
The finite element model was generated with a sweep mesh. The boundary
conditions for the model were an elastic foundation at the section where the wing is
attached to the fuselage.
In this model, a spring foundation was used to model the actual attachment
of the wing to the apprentice fuselage. The rubber bands provide a firm attachment
to the wing, however the elastic foundation influences the eigenstrain modes found
as a result of the finite element calculation. Some of the first few eigenstrains are
shown in Figure 3.42.
Based on the first six modes, the equally weighted eigenstrain and the MPF-v
weighted eigenstrain superpositions are shown in Figure 3.43.
From this mapping, which uses the first low frequency damped modes, some
125
(a) First eigenstrain mode (b) Second eigenstrain mode
(c) Third eigenstrain mode (d) Fourth eigenstrain mode
(e) Fifth eigenstrain mode (f) Sixth eigenstrain mode
Figure 3.42: First few eigenstrain modes of flight wing
126
Figure 3.43: Eigenstrain modes of flight wing: linear weighting (left) and MPF-v
weighted (right)
Figure 3.44: MPF-v weighted eigenstrain map for 8 sensors
127
asymmetry can be ignored, as this is an artifact of the damping mode structure of
the FE solver. The sensor placement diagram based on these low frequency, damped
modes, is shown in Figure 3.44. Other diagrams may be obtained depending on the
number of modes that are going to be considered. Typically if 80 % of the model
mass is accounted for in the mass of the modes used for the superposition then is
a reasonable approximation. If however more mass needs to be accounted for, then
possibly more energy considerations on the modes need to be made. This model also
depends on the approximation of the spring foundation. It is worth noting that the
dihedral reduces the observability of the vertical gusts on the surface of the wing.
Another factor which reduces observability of the gusts is the spring foundation, as
now the impact of the vertical gust disturbance is partially absorbed and damped by
the rubber bands which hold the wing to the fuselage. The flight scenario is therefore
more complex and more factors need to be taken into consideration for a reasonable
approximation to be made. Figure 3.44 shows the sensor placement of 8 sensors, four
on each wing, on the areas of maximum expected strain. A symmetric extension was
made based on approximations made above. In this case, the placement is solely
made on observability considerations. If wave propagation were going to be studied
with the sensor setup, then most probably some sensors would need to go towards
the tip of the wing.
128
Table 3.6: Mass participation factor v of carbon fiber wing
Eigenfrequency MPF-v
0.01436i 2.11168e-4
0.02834i 1.02403e-5
0.01454 -2.63717e-5
0.44147 4.48709e-5
2.9061 0.35438
4.01844 0.49673
72.96186 -0.00189
159.43176 -9.06982e-9
168.91207 -1.16206e-8
200.00018 1.48484e-4
215.37157 5.81196e-5
226.2739 -5.74395e-9
257.33485 8.38917e-5
278.86208 -8.44229e-9
292.43694 1.53964e-9
300.08143 3.85717e-5
324.85174 -1.41942e-9
326.8127 -1.47191e-5
331.12235 -1.0998e-9
340.1933 -3.427e-9
129
Chapter 4: Experimental section and analysis
4.1 Glider experiments
4.1.1 Experimental Setup
A 1D strain feedback Glider experiment was conducted to test the feasibility
of using strain feedback for gust alleviation. The one dimensional model used is
found in Section (2.3.1).
The experimental setup intended to measure a glider’s response to a vertical
discrete gust while on flight. The setup was composed of a glider, rail launcher,
a fan blower and a landing net (Figure 4.1). The rail launcher entails a sliding
launching cradle pulled by a bungee cord. The launcher is compact and can be
placed at different launching angles simply by adjusting the height at one end. An
inelastic cord is used as a stopper for the cradle as well as for keeping the bungee
cord under tension. The cradle is released from the receiver upon a transmitter
command (Figure 4.1). The bungee cord then transfers its elastic potential energy
to the launching cradle which then stops near the end of the railing giving its kinetic
energy to the glider.
All the experiments were synchronized such that the Labview acquisition would
130
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
1.5
2
2.5
3
3.5
4
4.5
5
Gust profile encountered by Glider
width (m)
 Gu
s
t 
ve
lo
c
it
y m
/s
Glider with aileron 
On right wing  Vicon marker tracking (rigid body motion: attitude and trajectory) 
Gust profile  Net-end of trajectory 
Max.~ 5 m/s 
Figure 4.1: Experimental setup of glider and launching mechanism
start recording the glider trajectory upon launch. A MOTE transmitter was con-
nected to Labview in order to have repeatable launches. The glider was consistently
launched at a velocity of 0.6 m/s in the x-y plane (floor plane) in all the exper-
iments as a result of this mechanism. The strain gauges were mounted onto two
plate spars which were manufactured using a fast prototyping machine using Garo-
lite (Figure 4.2). These spars were fitted with full bridge strain gauges. The wings
were mounted onto a spar on each side of the central plate which also carried the
processing unit and XBee transmitter. The strain gauges installed were of 1000
Ohms in full bridge per each wing, 2 in tension and 2 in compression.
The microcontroller system on board performs all of the A/D conversion,
PWM generation and computation necessary for in flight feedback from up to 4 sep-
arate strain gauge channels to up to 4 separate PWM servo channels. The aileron
131
Strain  
Gauges 
To 
Flapperon 
Stiffener 
Plate 
Figure 4.2: Glider plate and strain gauges
servo was also connected to the microcontroller board. The connections were placed
inside the foam to help maintain laminar characteristics around the airfoil.
Figure 4.3: Experimental setup of glider and launching mechanism
The dynamics acquisition was implemented in Labview using VICON marker
132
data. All experimental results were obtained using the VICON motion capture sys-
tem of the Autonomous Vehicle Lab at the University of Maryland College Park
campus. The VICON motion capture system records position, attitude and orienta-
tion by triangulating marker points Figure 4.4. In the case of the glider, the markers
were distributed along the fuselage to avoid any error due to local deformation. This
would have been the case had the markers been placed on the flexible wings. A cal-
ibration of the VICON system was required before collecting the data and needed
to be recalibrated often to avoid errors due to marker displacements produced dur-
ing vehicle landings. The calibration included defining the vehicle’s marker setup
as well as defining the axes in a north, east down system such that the positive
x-direction corresponded to the direction in which the glider was launched and the
positive y-direction entailed the right hand side of the testing area. Adjustments
due to this nomenclature were accounted for in the data presented in this work.
Figure 4.4: VICON Setup: Vehicle marker setup(left) and VICON camera setup and
vehicle (right)
133
Roll dynamics were measured mainly for three sets of flight cases; free flight,
flight under gusting conditions and flight under gusting conditions using strain feed-
back. In particular the experiments focused on the effect of the application of strain
feedback to regulate roll. Different proportional gains were applied to investigate the
effects of strain feedback on the rolling motion of the glider under gust conditions.
The strain feedback results were also compared to classical controllers which were
implemented off board using the VICON system/Labview setup to send updated
aileron commands.
In order to achieve a good gust setup a couple of possibilities were evaluated.
First, a box fan was set as the gust source and the glider projected path directly over
it. This however doesn’t provide a sharp edged gust which is desired to generate a
strong rolling motion. Then a plate was incorporated to shape the gust. In this case
the flow was not strong enough and the air stream near the plate was turbulent.
Finally a leaf blower fan placed next to the plate was chosen for its more laminar
qualities and strength. Tufts were used to heuristically determine the orientation of
the flow field. The gust and plate were placed at 40 inches from the launch point.
Figure 4.5, shows the gust contour utilized for the tests. It resembles a 1-cos gust,
which is a commonly used gust model type. The measurements of the flow field
were obtained using a hand anemometer. The gust had an approximate peak value
of 5 m/s. After passing throught the sharp gust, the glider would roll and land on
a net set ahead of the gusting region. The VICON flight dynamics and trajectories
were recorded for each flight.
134
Figure 4.5: Gust encountered by glider.
4.1.2 Controls Setup
The controls implemented for the strain feedback closed loop tests is shown
in Figure 4.6. It entails a combination of off-board calculations and an on-board
calculation using strain sensing.
Vehicle  
States 
VICON  
data  
Ground Control  
Strain 
Gauges Calculation of  
Roll command  
Torque  
Calculation  
Calculation of  
Roll torque  
XBee  
transmission  
Strain  
Sensing 
Aileron  
Command  
Calculation of  
Aileron command  
Aileron  
Servo  
Roll  
Cmd  On Board Control  
Aileron  
Figure 4.6: Controls used in strain feedback tests.
Figure 4.6 shows the controls used in strain feedback tests.A reference roll
command is defined in the ground control which uses data from VICON. This in-
135
formation is used in the calculation of a desired roll torque which is sent via Xbee
to the on board controller. Once the desired torque has been transmitted to the
on-board controller, a calculation for the aileron command ensues. Finally this value
reaches the servo and changes the aileron angle, producing changes in the vehicle’s
states, which are being recorded in VICON at all times.
In order to compare the results obtained using strain feedback, a couple of
flights using classical controllers were undertaken. The classical controllers are clas-
sical in the sense that they use the recorded rigid body dynamics to calculate an
aileron response. For these tests, the controls scheme utilized is shown in Figure 4.7.
A comparison is shown later on in the experimental results section.
Vehicle  
States VICON  
data  
Calculation of  
Roll command  
Torque  
Calculation  
Calculation of  
Roll torque  
XBee  
transmission  
Aileron  
Command  
Calculation of  
Aileron command  
Aileron  
Servo  
Roll  
Cmd  
Aileron  
Figure 4.7: Classical controller used for controller types comparison.
4.1.3 Experimental Results
4.1.4 Force/Strain Feedback using Proportional Gains
The experimental results are presented in this section. Ten roll curves were
recorded in VICON for each case. Figure 4.8 a shows the free flight history of the
Glider. In this graph the four distinct flight phases are shown. First, the glider
136
remains on the launcher cradle (smooth horizontal region) after which it is launched
and the roll angle increases as the trajectory advances. The short flight ends when
the glider hits the net, generating scattered data points. These points however are
discarded for calculations. The free flight data set is used later on in this work as a
baseline for comparing the effectiveness of the structural feedback controller.
0 20 40 60 80 100 120 140 160−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Time (1 unit = 10 ms)
Roll A
ngle (R
adians)
Baseline Roll angle tests with no gust and no control
 
 Trial 1Trial 2Trial 3Trial 4Trial 5Trial 6Trial 7Trial 8Trial 9Trial 10
0 20 40 60 80 100 120 140 160−2
−1.5
−1
−0.5
0
0.5
1
1.5
2
Time (1 unit = 10 ms)
Roll A
ngle (R
adians)
Baseline Roll angle tests with gust disturbance and no control
 
 
Trial 1Trial 2Trial 3Trial 4Trial 5Trial 6Trial 7Trial 8Trial 9Trial 10
Figure 4.8: Glider baseline roll angle curves: with no gust perturbation (top), and
under a vertical asymmetric gust (bottom).
137
A similar test was carried out for the case where no feedback control was
applied in the presence of the sharp gust. The gust produces a maximum roll angle
in the order of 45 degrees. These tests are shown in Figure 4.8. The effect of the gust
can be clearly seen from Figure 4.8 (bottom). In Figure 4.9, each curve represents
an average of the 10 flights for each case; free flight and with gust perturbation.
0 20 40 60 80 100 120 140 160−0.4
−0.2
0
0.2
0.4
0.6
0.8
Time (1 unit=1ms)
Mean
 roll a
ngle (ra
dians)
Mean roll angle with and without gust perturbation
 
 
Free flightPerturbedNo−control
Figure 4.9: Glider mean roll for free flight and with vertical gust perturbation.
After measuring the impact of the gust on the glider’s roll angle a structural
feedback controller was implemented. The controller intends to bring the roll angle
back to an unperturbed baseline behavior. The major impact of the structural
feedback correction of the roll angle can be seen towards the end of the graphs due
to the shortness of flight duration. Figure 4.10 shows the effect of the controller
at different proportional gains. Each one of the curves shown in these figures is an
average of 10 to 20 flights. The effect of the structural feedback controller is shown
to be higher at higher gains.
138
0 50 100 150−0.4
−0.20
0.20.4
0.60.8
Time (1 unit=1ms)
Roll ang
le (radians
)
Gust disturbance and strain feedback K=3
 
 PerturbedNo−controlK=3
0 50 100 150−0.4
−0.20
0.20.4
0.60.8
Time (1 unit=1ms)
Roll ang
le (radians
)
Gust disturbance and strain feedback K=4
 
 PerturbedNo−controlK=4
0 50 100 150−0.4
−0.20
0.20.4
0.60.8
Time (1 unit=1ms)
Roll ang
le (radians
)
Gust disturbance and strain feedback K=5
 
 PerturbedNo−controlK=5
0 50 100 150−0.4
−0.20
0.20.4
0.60.8
Time (1 unit=1ms)
Roll ang
le (radians
)
Gust disturbance and strain feedback K=6
 
 PerturbedNo−controlK=6
Figure 4.10: Glider mean roll for proportional gains K=3, K=4, K=5, and K=6
Figure 4.11 shows a summary of the applied gains and compares it to the
baseline curves. This figure shows that when the gain is increased, the maximum
perturbed roll angle is reduced and the roll angle decreases faster.
The effectiveness of the structural feedback can be quantified using different
metrics. We have chosen the D2 metric for this purpose. The D2 metric provides a
summation of the differences in the trajectories of the glider for the perturbed case
with no applied control and the case where the controller is applied. Figure 4.12
shows that the D2 metric increases with the applied structural gain. This implies
that the glider roll angle will increasingly stay away from the roll angle corresponding
to the perturbed case.
As a way to compare strain feedback with classical controllers, we conducted
139
0 50 100 150−0.4
−0.2
0
0.2
0.4
0.6
0.8
Time (1 unit=1ms)
Roll
 ang
le (ra
dians
)
Strain feedback correction of roll angle under gust disturbance
 
 
Free flight
PerturbedNo−control
K=4
K=5
K=6
Figure 4.11: Glider roll angle and strain feedback at gains K=3, 4, 5 ,6 and baseline
flight.
1 2 3 4 5 6 70
0.20.4
0.60.8
11.2
1.41.6
1.8
Gain
d2met
ric
D2 metric between median curves as a function of controller gain
Figure 4.12: D2 metric between the perturbed data with no applied control and the
one with applied strain feedback.
140
four tests. The first two tests included a classical controller where VICON provided
the feedback and the following two tests included on board strain feedback. From
Figure 4.13 we can see that the two controllers with strain feedback are more effective
than the ones based on classical state feedback.
The effect of applying strain feedback can also be observed in the planar
trajectory followed by the glider (Figure 4.14). This graph shows a representative
curve of the x-y plane directional velocity evolution of the glider for three cases:
under no gusting conditions, under gusting conditions and under gusting conditions
while using strain feedback.
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
−20
0
20
40
60
80
Time(/10s)
Roll A
ngle (D
egrees)
Gust rejection controller types comparison 
 
 
No Controller (Trimmed State)Classical P Controller (Outer Loop, No TF)TF P Controller (Inner Loop Only)TF P Controller and VICON P Controller (Outer +Inner Loop)
Figure 4.13: Controller comparison for roll gust rejection. TF= torque feedback,
P=proportional.
Weights were attached to the fuselage of the glider. They were placed near
the center of mass and in a way that they would not affect the aileron mechanism or
produce any additional imbalance. The roll angle behavior was observed by testing
141
−3 −2 −1 0 1 2 3−2
−1.5
−1
−0.50
0.51
1.5
2
2.5
x−displacement
y−dis
place
ment
Glider trajectories−top view
 
 
free−flight trajectorywith gust and no controlwith gust and control (K=5)
Figure 4.14: Glider trajectories: in free flight, under the influence of a gust and using
strain feedback.
the glider without the gust disturbance, with the disturbance and with both the
disturbance and strain feedback to measure controller effectiveness. The following
graphs correspond to the three vehicle weights.
Each curve represents the average of ten runs of the glider. The two following
graphs amount to forty glider flights (Figure 4.15). These results show that strain
feedback is robust to small payload alterations. As an observation, the weights add a
different roll time history behavior, which is consistent within each weight. Another
trend observed is that at higher loads, the effect of the proportional controller seems
to be decreasing.
142
0 50 100 150−0.2
0
0.2
0.4
0.6
0.8
Roll angle for Weight 1
Time (1 unit = 10 ms)
Roll 
Angle
 (Radi
ans)
 
 No gust, no controlGust disturbance, no control
 with applied control, K=5
0 20 40 60 80 100 120 140−0.2
0
0.2
0.4
0.6
0.8
Roll angle for Weight 2
Time (1 unit = 10 ms)
Roll 
Angle
 (Radi
ans)
 
 
No gust, no controlGust disturbance, no control
 with applied control, K=5
Figure 4.15: Impact of weight on glider roll angle: Roll angle for glider weight and a
payload of 4.4 grams (top) and Roll angle for glider weight and a payload
of 10.9 grams (bottom)
143
4.1.5 Aileron System Identification
4.1.6 Static Proportional Model
A static proportional model for the aileron can be generated by providing step
inputs to the wing and recording the corresponding aileron deflection outputs. The
glider was placed on a solid surface and two of the vicon cameras were brought
up closer for increased resolution of the motion of wing and aileron relative to the
fuselage. For this purpose a three rigid body marker setup was devised in order to
record the corresponding angles of each (Figure 4.16).
fuselage 
Wing  
Aileron  
Vicon  
Camera  
Setup  
For  
System  
Identification  
Figure 4.16: Vicon camera setup for system identification
144
Figure 4.17 shows the raw VICON data for a strain feedback gain of K=1. It
depicts a series of step inputs of wing deflections and their corresponding output of
aileron angular deflections. A static model was derived from these tests by taking
an average of values for each step of aileron deflection and correlating it to the wing
angular deflection.
15 20 25 30 35 40
−0.3
−0.25
−0.2
−0.15
−0.1
−0.05
0
0.05
0.1
0.15 Ouput/input in time domain around mean
Time (s)
Angle 
(rads)
 
 Output−Aileron angle (rads) around meanInput−wing deflection (rads)
Figure 4.17: Step input of wing deflections and corresponding output of aileron de-
flections.
The static model is shown in Equation (4.16). It provides a direct correlation
between the input wing deflection and the output aileron deflection by means of
proportional constants. There is a different model for every strain feedback gain K.
A gain of K=1 had less error produced by oscillations in the servo than those seen
on higher gains such as K = 3, 4, 5. The model results for gains K=1 and K=3 are
145
shown in Table 4.4 and Table 4.2.
∆Wing = C ∗∆Ail + Co
∆Ail =
1
C
∆Wing +Do = D ∗∆Wing +Do
C = 5.7429, Co = 0.4090, D = 0.1734, Do = −0.0711
(4.1)
A fit for the raw data and the comparison with the extracted data points is
shown in Figure 4.18 for K=1. These values however depend on the method used
for data reduction.
−0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1−0.35
−0.3
−0.25
−0.2
−0.15
−0.1
−0.05
0
0.05
0.1 Aileron deflection θ vs Input wing deflection φ
Input wing deflection (rads)
Outpu
t ailero
n defle
ction (ra
ds)
Figure 4.18: Aileron deflection vs wing deflection for K=1. Each point is the result of
averaging local points that correspond to the duration of the step input.
Different methods were applied to obtain the model parameters. An average,
filtered and smoothed roll angle history was used in order to obtain the four con-
146
stants. The filter applied is a second order Butterworth with a sampling frequency of
33 Hz which is the VICON data rate and a cut-off frequency of 10 Hz. The smooth-
ing algorithm applied was a local regression using weighted linear least squares and
a 1st degree polynomial model. A comparison of the model fit parameters of inter-
est using different data processing methods is shown in Table 4.4 for K=1. Similar
results and method comparison are shown in Table 4.2 for K=3.
Table 4.1: Comparison of linear fit with different processing methods for static model
and K=1.
Parameter θ θ mean θ filtered θ filtered mean θ smooth θ smooth mean
C 1.9237 1.9237 1.9231 1.9231 1.9459 1.9459
Co 0.1139 -0.1291 0.1139 -0.1290 0.1148 -0.1308
D 0.5171 0.5171 0.5172 0.5172 0.5106 0.5106
Do -0.0589 0.0667 -0.0589 0.0667 -0.0586 0.0668
Table 4.2: Comparison of linear fit with different processing methods for static model
and K=3.
Parameter θ θ mean θ filtered θ filtered mean θ smooth θ smooth mean
C 5.7429 5.7429 5.7531 5.7531 5.8142 5.8142
Co 0.4090 0.2022 0.4093 0.2025 0.4039 0.1974
D 0.1734 0.1734 0.1730 0.1730 0.1690 0.1690
Do -0.0711 -0.0352 -0.0710 -0.0352 -0.0689 -0.0340
147
4.1.7 Dynamic Model
A dynamic model was obtained for the input wing deflection as a function
of the output aileron deflection. Inputs of different frequencies were provided to
the wing in order to characterize the aileron in its frequency response. Figure 4.19
shows the aileron angle and corresponding wing deflection time histories measured
in VICON.
0 5 10 15 20 25 30 35 40−0.4
−0.2
0
0.2
0.4
0.6 Ouput/input dynamic aileron data in time domain
Time (s)
Angl
e (rad
s)
 
 Output−Aileron angle (rads) around meanInput−Wing deflection (rads)
Figure 4.19: Filtered data, where aileron output has been drawn around its mean
value
A system identification of the aileron dynamics was obtained using the follow-
ing model which includes a delay term and one pole.
H(s) =
Kp
1 + Tps
∗ e−Tds (4.2)
Where H(s) is the transfer function that relates wing deflection to aileron
148
deflection:
H(s) =
Output
input
=
∆Aileron
∆Wing
(4.3)
For the filtered data, the following coefficients were obtained for the model
shown in Equation (4.17)
Kp = 1.5564
Tp = 0.0224
Td = 0.0869
Lossfunction = 0.01306
FPE = 0.01309
(4.4)
For completeness and to enable simulation capabilities of the current system,
a state space model was generated using a Pade approximation of order n=1 for the
delay term:
e−Tds =
1− Tps
1 + Tps
(4.5)
Thus, the following model can be used for the aileron transfer function by
using the approximation in Equation (4.5) and into Equation (4.17):
H(s) =
Kp
1 + Tps
∗
1− Tps
1 + Tps
(4.6)
149
For this transfer function a time domain model can be obtained by applying
an inverse Laplace transform, where the variables of interest are Φ, or the wing
deflection angle and Θ, or the aileron deflection angle. The differential equation that
relates the aileron and wing deflection can be written as shown in Equation (4.6).
Kp ∗ Φ(t)−
Td
2
∗
dΦ(t)
dt
= Θ(t) +
(
Tp −
Td
2
)
∗
dΘ(t)
dt
−
TdTp
2
∗
d2Θ
dt2
(4.7)
A state space representation can also be obtained for the model with delay:
A =




− 2TdTp ∗
(
Td
2 + Tp
)
− 2TdTp
1 0



 , B =




1
0



 , C =
[
−KpTp
2Kp
TdTp
]
, D = [0]
(4.8)
Or numerically for the last case of filtered data:
A =




−2.5403 −0.7678
1 0



 , B =




1
0



 , C =
[
−0.0201 0.0449
]
, D = [0]
(4.9)
The delay terms show that the aileron is fast enough to be used in conjunction
with the fast strain based loop. The delay constant is of the same order of magnitude
as the observed time required to turn the rigid body of the glider after a gust input.
A faster servo could decrease the overall system reaction time if the choice of strain
sensors and substrate conserve their high bandwidth.
The time scale of the structural response is much less than that of the rigid
150
body response. This time scale was measured using high speed videos of the glider
passing over the gust. The fuselage starts to turn roughly 0.03 seconds after the gust
hits the left wing and it levels out completely at 0.06 seconds. Thus, demonstrating
the advantage of using structural feedback as a way to achieve reflexive control. This
enables the implementation of a faster inner loop configuration. This loop can be
implemented before the plant dynamics, in a way that it absorbs unknown inputs
and disturbances. An inner loop force adaptive feedback approach.
Strain force feedback is demonstrated using a small UAV under a sharp gust
condition. The vehicle was instrumented with strain gauges which were sensitive to
vertical gust disturbances. The sensor output was used in a fast inner loop for reflex-
ive reaction to gusts. Linear strain measurement weightings or wide field integrated
strain can provide force state estimations which can be used in these fast inner loops.
This feedback is much more effective at higher proportional gains in the case of a
strain/force proportional controller. Less deviation in the roll trajectories was ob-
served by means of the D2 metric as the gain was increased. Structural feedback is
also robust to small payload additions, affecting the roll rate but conserving the fast
reaction characteristics. Structural feedback also demonstrated to be more effective
than classical rigid body dynamics controllers. The classical controllers were based
on VICON state feedback of the rigid body states, which would be equivalent to
on board IMU’s attached to the core of the fuselage. A system identification of the
aileron dynamics showed that the time constant is of the same order of magnitude
as the time that it takes for the disturbance input to propagate to the rigid body. A
faster servo would probably show a faster response to correct for disturbances, albeit
151
it also depends on the magnitude of the disturbance. Aileron dynamics were fast
enough for the employed combination of gust speed and roll correction. The aileron
was sufficiently fast to respond to the structural feedback inner loop. Structural
feedback using wide field integrated constants is a promising mechanism for UAV’s
to reject gust disturbances as they allow a faster reaction than conventional control
mechanisms.
4.2 Quadrotor experiments
A similar bio-inspired approach to strain sensing-based gust rejection can be
applied to quadrotor vehicles. In this case, the sensors and the underlying structure,
or quadrotor arm, need to be designed for such purpose. In quadrotors there are
added challenges such as additional and varying motor vibration as well as aerody-
namic phenomena present during flight and as a consequence of the rotor induced
aerodynamics. A quad strain sensing arm was developed to be used in gust mitiga-
tion, as shown in Figure 4.20.
In the case of strain sensing for a quadrotor, design and characterization of
substrate/sensor pairs is required in order to obtain optimal measurements. This
refers to the best locations, surface dimensions and geometries that can be used
for installation of strain gauges. A natural selection of these locations is found in
nature, these locations correspond to areas of high density of surface stresses or areas
of higher stimuli inputs. Given the small time scales of the structural responses as
compared to rigid body responses, a structural feedback scheme can be proposed
152
to reject gust disturbances. Moderately flexible quad arms are needed while also
preserving the control authority of the motors. This flexibility allows for encoding
of the forces and moments being exerted on the quad arms in the form of strains
and stresses.
Figure 4.20: DJI Flamewheel quadrotor with strain sensing arm
The forces and moments anticipate and precede the rigid body motion of the
quadrotor. Therefore, the mapping between forces and strains enables the imple-
mentation of a fast inner loop which can be used for disturbance rejection. The
resulting strain measurements are functions of the position and orientation of the
sensors, the structural modes of the airframe, and the forces and torques, which may
be externally applied or commanded. Assuming the structural loads are primarily on
the quadrotor arms, a network of strain sensors can be used to anticipate the quadro-
tor rigid body rotations and displacements. In this section, a heave perturbation to
153
the quadrotor is encoded in a set of strain gauges. The strain information associ-
ated to encoded forces is used in a fast inner loop which corrects for the disturbance
before it reaches the quadrotor rigid motion states. That is, it allows for a correc-
tion of the disturbance before it propagates to lower order states. Strain sensing
has not been explored for quadrotor UAVs, partially due to size, weight and power
constraints but also due to the complex structural effects it endures during flight.
The quadrotor arms are subjected to a number of forces, both aerodynamic [104]
and inherently mechanic. The four motors introduce not only mechanical vibra-
tions but also generate turbulent air patterns [104] which are constantly hitting the
quadrotor arms. The effects also vary depending on whether the quadrotor is close
to the ground, near a wall, or if it has forward velocity. Other effects inherent to
rotor motion such as feathering, flapping and lead-lag can also affect the strain time
history. Atmospheric disturbances will also introduce additional force components
which will also be seen in the strain measurements.
4.2.1 Strain Sensing
Structural deformations are gradients of three dimensional displacements. For
linear regimes shear deformations can be neglected allowing for a simpler formu-
lation. A quadrotor arm can be thought of as an Euler Bernoulli beam type of
structure cantilevered at the center plate. The fixed condition however is not en-
tirely met due to the flexible nature of the center plate, designed mainly to robustly
host electronics and batteries. As would be expected, the locations of higher strains
154
for an isolated beam are close to the root, or the clamped end. These higher strain
regions are mainly due to the nature of the predominant influence of thrust vectors
produced by the rotors, which are typically perpendicular to the plane formed by
the quadrotor arms. A one dimensional beam formulation can approximately model
the strains produced by the thrust vector forces. In a quasi-static formulation, the
displacement w at any location x along the arm spanwise direction produced by
applied forces P is described by:
d2
dx2
(EI
d2w
dx2
) = P (x) (4.10)
And the bending moment is given by:
d2
dx2
(M(x)) = P (x) (4.11)
This uniaxial approximation is made for simplicity, as the actual strains are
produced in all directions within the material. The larger strains however are along
this axial or spanwise direction. This formulation also becomes complex when con-
sidering beams with varying cross-sections and stiffenesses. The neutral axis be-
comes a function of the axial coordinate and a piecewise formulation is then more
appropriate.
In particular, for a beam of uniform cross-section and a force applied at a
location x on the quad arm of length L, the force and strain are related to each
other by means of
ε(x) =
zl
EI
∗ P (x) (4.12)
155
where the force is a function of lengthwise location of the force and zl refers
to the distance from the neutral axis of the surface at which the strain is being
recorded. If the strain is to be evaluated at a location Lo from the fixed end of the
arm, then the strain at this point becomes
ε(L0) =
zl
EI
∗ P ∗ (L0 − L) (4.13)
where P is the magnitude of the applied force. Similarly, the associated bend-
ing moment at the location L0 can be calculated as
M(L0) =
EI
zl
∗ ε(L0) (4.14)
Biased values of strain are present due to different factors which include:
weight, flight maneuvers or certain flight conditions. These however can be sub-
tracted so that the measurement vector, or vector of measured strains can be cal-
culated as in Equation (3.1) and Equation (3.2), as presented in [6]:
yε = Cεxε (4.15)
The dimension of the output y is related to the number of independent mea-
surements of strain. The absolute state vector xε is the force and moment associated
to all the forces and moments being exerted on the quadrotor at any given time:
156
xεT ≡
[
Px Mx
]T
≡
[
P (L0) M(L0)
]T
(4.16)
If we apply the bias, then the state vector directly correlates to the heave
disturbance, for this particular implementation.
xε ≡
[
∆P ∆M
]T
(4.17)
and the output matrix Cε is a function of the location L0 of the strain gauges:
Cε =




EI
zl
1
L0−L
EI
zl



 (4.18)
Such that the states related to the heave disturbance can be estimated as:
xˆε = Myε = (CTε Cε)
−1CTε yε (4.19)
or equivalently:




∆ˆP
∆ˆM



 = Cεxˆε, (4.20)
These states can be used in force feedback to attenuate the effect of distur-
157
bances and return the system to the unperturbed states. Distributing sensors on
more locations along the quadrotor arm may yield better estimations of the states
associated with the disturbance.
4.2.2 Sensor design for quadrotor arm
Sensor design and associated substrate optimization need to be considered in
order to achieve higher sensitivities to the targeted disturbances. In addition to this,
the overall bending stiffness of the quadrotor needs to be kept close to the one of the
original DJI quadrotor arm. This trade of between strain sensitivity and stiffness
was observed during design iterations. First, the possible loads on the quadrotor
arm were determined. This amounts to the loads due to gravity and the loads due
to flight maneuvering. Five maneuver tests, ranging from sharp turning to normal
hover, were recorded in VICON motion capture system. Analysis of the rigid body
states showed that the forces do not exceed approximately 0.2 N per arm. This force
range was utilized in the design of the quadrotor arms to maximize load sensitivity.
Characteristics of the DJI quadrotor arm are presented in table 4.3. The
bending stiffness characteristics were measured using an MTS machine and a three
point bending setup. In this procedure the quadrotor arm was placed evenly on
top of two supports. The loading point is found half way between the two supports.
Force and displacement data are recorded by means of a load cell and a displacement
gauge, respectively. The test was carried in a destructive way such that the yield
stress will also be measured. In the linear region of the graph, a direct relationship
158
is found between the force and displacement and this value is proportional to the
bending stiffness if the cross-section characteristics are known. It was interesting to
observe that by design, the lateral bending stiffness of the quadrotor is larger than
the longitudinal bending stiffness.
Table 4.3: Comparison of quad arm characteristics of Original DJI arm and proto-
typed arm.
Parameter Original DJI Arm RP Quad Arm
Composition Thermoplastic Photopolymer
Young’s Modulus 15000 MPa 2000-3000 MPa
Tensile strength 160 MPa 50-65 MPa
Density 1370 kg/m3 1170 kg/m3
Bending Stiffness (EI) 1.61 Nm
2.48 Nm (solid), 1.43 Nm (I-beam),
1.32 Nm (hollow)
Once the bending stiffness was measured, quadrotor arms of different cross
sections were built using rapid prototyping (Figure 4.23). The material used was
a photopolymer with medium stiffness properties. Static properties of the rapid
prototyped arm and sensors were obtained by means of a tip loading and unloading
set of tests at the Rotorcraft Center [105]. The setup is shown in Figure ??. Hys-
teresis was observed in the measurements. The creep behavior of the photopolymer
of the RPM material was also observed. The creep scale however was larger than
the expected time scales for our experiments. Results from these tests are found in
Figure 4.22.
159
Figure 4.21: (a) Quad arm with amplifier circuit, (b) Static testing of RPM quadrotor
arm
A series of tests with different cross-sections were performed until a good
candidate was found for maintaining good weight to strength ratio. Initially a solid
cross-section was used, this however did not yield good sensor sensitivity to loading.
Therefore a second part of the substrate optimization encompassed a local cross-
section optimization.
This local optimization was carried by means of finite element simulations of
the rapid prototyping material with different geometrical parameters. In order to
increase the strain output a local increase in strain needed to be generated without
reaching the material yield stress. This was achieved by matching high strain levels
and not so large stress surface distributions. It was observed that the best char-
acteristics were achieved by using an almost elliptical hollow section. Three of the
cross section designs are shown in figure 4.24. The best sensor locations were shown
160
Figure 4.22: (a)Hysteresis at a loading and unloading cycle, (b) Creep behavior of
RPM quad arm
a b  
Figure 4.23: (a) Photopolymer prototyped arm sections for MTS testing and (b) CAD
model of the quadrotor arm without optimized cross-section
to remain close to the highest curvature point of the hollow section near the root.
After determining the proper locations for the sensors a proper grid size and
backing material were chosen. The sensors used in this implementation were Omega
1000 Ohms dual strain gauges. These sensors are typical foil gauges with a gauge
factor of 2.02 and a tolerance of 0.3 %. These linear strain gauges were implemented
in a wheatstone bridge configuration in tension/compression pairs, and bonded to
the polymer substrate at the optimal locations. A surface treatment was necessary
to match the foil backing material to the photopolymer surface. The photopolymer
161
Figure 4.24: Finite element simulations of stress under maximum expected quadrotor
load: (a) Solid cross section, (b) Square hollow cross section and (c)
Smoothed square hollow section, (d) striated cross section, (e) load cell
cross section, (f)bi-rectangular cross section
was susceptible to creep typical of plastics but was accounted for by using updated
measurements. The bridge implementation is shown in figure 4.25. The bridge
output was amplified at a biased voltage so that the measured values of strain
can oscillate around a nominal value. Ground effects were observed in the strain
measurements as well as high altitude effects, where the rotor downwash had an
influence over the measurements.
It is worth mentioning that in this implementation there was a need to have
landing gear to protect the prototyped arm from undergoing sharp and concentrated
loads at the optimized sensor locations. The photopolymer looses strength over
162
a b  
Strain  
Sensing  
Quad arm  
Figure 4.25: Strain sensing implementation: (a) Quad arms with different cross sec-
tions and (b)Quad arm vehicle implementation
time, increasing sensor output but also making the prototype prone to breakage.
Therefore other material substrates can be considered to avoid these challenges, i.e.
carbon fiber, fiber glass, aluminum shells, etc...
4.2.3 Heave disturbance rejection
A robust control scheme is proposed in Figure 4.26. This controller was used
by [6] to demonstrate disturbance rejection with accelerometer data in correcting
roll perturbations. In this proposed scheme, the disturbance gets to be mitigated
way before it reaches the plant.
This feedback improves tracking of requested controls in the presence of dis-
turbances by regulating errors between desired torques and actual torques. Thus the
system is simultaneously capable of rejecting external disturbances, such as gusts,
and internal disturbances, such as actuator variations.
An artificial heave gust was implemented as a collective motor perturbation,
163
Figure 4.26: Fast inner loop for disturbance rejection
inducing the vehicle to lose altitude as might be experienced in a downdraft. An L1
norm was calculated for a sample of the perturbation intervals. The corresponding
strain output and the nominal average strain were used for this calculation, which
is shown in Figure 4.27.
Figure 4.27: Strain feedback performance metric
A control augmentation feedback control from the strain gauge network affixed
to the vehicle arm was implemented to return the vehicle to a previous strain state
corresponding to hover and thereby reject the disturbance. The improvement in
164
maintaining altitude can also be observed from the time histories for the altitude
and vertical velocity states for increasing feedback gains (Figure 4.28).
Figure 4.28: Heave Velocity and Altitude Time History for Impulse Heave Disturbance
for Varying Feedback Gain from Strain Gauge Sensing
Using novel sensor implementations of strain and acceleration measurements,
estimates of the forces and torques applied to a quadrotor sUAS were estimated.
Furthermore, these estimates were provided through static linear combination of the
raw measurements, yielding a computationally simple and analog implementable es-
timation architecture. The response of the systems was characterized in the presence
of force and torque stimuli. Feedback from these estimates was used to mitigate the
effect of motor perturbations in the roll and heave degrees of freedom. Two charac-
teristics were observed from these feedback implementations. Undesirable exogenous
disturbances, motor perturbations in this case, are attenuated in the output state
of the vehicle. However, desired actuation forces are also attenuated, resulting in
a slower reference tracking response. Future implementations of this feedback con-
165
trol will account for this undesired attenuation by removing contributions from the
intended actuator forces and torques in the control adaptation.
4.3 Wing structural properties testing
The carbon fiber wing described in Section (3.6.1.1) was manufactured by [102]
with the Aurora PWings design [7], which follows an AG-35 airfoil and has the
geometric characteristics shown in Table 3.5. The wing entails a foam core, no spars
or ribs, and a layup of carbon fiber and fiber glass.
The wing was tested in its static structural properties. Bending and torsion
tests yielded stiffness parameters for comparison with the finite element model.
The Vicon tests for structural wing property measurements were conducted at the
Rotorcraft center of the UMD, College Park [105]. The static load bending and
torsion tests were conducted at Aerolabs, UMD, College Park [106].
4.3.1 Experimental setup
In order to measure the bending stiffness of the carbon layup wing, two main
setups were used to measure vertical deflection and angular torsion. For bending,
static tip deflection and section wise vicon deflections were used. For torsion, two
methods were used, the mirror method and the angular displacement method.
The first setup consists of a height gage and a number of loads, as shown in
Figure 4.29. The loads are placed at the wingtip and the vertical tip displacement
is measured with the height gage. The maximum bending load was close to 1 kg,
166
therefore the weights used were only up to about 700 gr. It is worth mentioning
that these measurements were subject to errors stemming from clamping method,
loading speed, loading intervals, material creep from the wing and reading errors.
The reading error is less than 5 % and the other errors can be considered systematic.
The second setup is shown inFigure 4.30 and it consists of a Vicon system,
reflective markers and a number of loads . For the vicon marker placement the wing
had to be divided into equal length spanwise sections. Three markers were placed
at each spanwise location. The camera setup for the vicon test consists of four
cameras that are rather close to the wing as the volume needed to be minimized.
The surface of the wing had to be covered in non-reflective material, as the epoxy
finish of the wing induced errors in the marker position estimation of the Vicon
setup Figure 4.31.
Figure 4.29: Height gauge setup for calculation of bending stiffness
Both setups used reinforced RPM conformal clamps. These are clamps which
167
Figure 4.30: VICON setup for calculation of bending stiffness
have the exact airfoil shape to distribute the load evenly and provide appropriate
boundary conditions. The conformal clamps were placed at the root with aluminum
reinforcements, as the RPM material wasn’t stable enough for the loads applied and
produced artificially soft modulus results due to slippage. After the reinforcements
were included, results were much better. Another conformal clamp with reinforce-
ments was placed at the tip of the wing for torsion tests. The top reinforcement
of the tip conformal clamp was much longer to be able to produce a long moment
arm and therefore larger torsion. These clamps are shown in Figure 4.31. The an-
gle measurements for the torsion tests were measured with an electronic protractor
which was placed at the wingtip.
168
(a) Clamps used for static testing (b) VICON wing setup
Figure 4.31: Clamps used for static testing and Vicon wing setup
4.3.2 Bending and torsion properties
The testing procedure for the first setup consisted of recording the vertical
displacement of the tip of the airfoil under increasing load. For this purpose, the
following formula was applied and a least squares estimation of the slope was gen-
erated to calculate EI. The results can be observed from Figure 4.32.
After averaging the estimated least squares slopes, the bending stiffness EI for
each can be calculated using Equation (4.21).
EI =
PL3
3wt
(4.21)
Where P is the applied load, L is the effective length of the cantilever and wt
169
Figure 4.32: Load vs. tip displacement curves for estimation of EI with height gauge
is the tip deflection. The average calculated bending stiffness from the three curves
is then calculated to be EI = 25.4949Nm2 .
The second setup entails a similar test in that loads are placed at the tip
of the wing and the deflection evaluated. The deflection is now measured using
Vicon nexus for deflection tracking [105]. Four sets of three markers are placed
on the quarter chord along the length of the wing. These three points are used
in a Lagrange polynomial to calculate the corresponding deflection. The bending
stiffness is now calculated using position on the wing length and deflection slope.
This can be calculated using Equation (4.22) which entails the deflection slope of
the wing.
w′ =
P
EI
∗
(
Lx−
x2
2
)
(4.22)
170
By using this formula on every three markers, we can get an approximation
of the stiffness at the four different sections that the markers span. The thirteen
marker positions are depicted in Figure 4.33. Each deflection value is an average of
deflections taken by vicon at a 60Hz rate. The last three sets of markers had a lot
of noise. This could mean that these last markers which were closer to the wing tip
were probably not placed smoothly over the surface of the wing.
Figure 4.33: VICON processed data
Taking into consideration the first three averaged bending stiffness values we
can calculate an average of these, yielding EI = 31.10949Nm2 .
This value is very close to the one measured by the height gauge. The accuracy
of the height gage might be better than the vicon measurements in that the markers
used were not the optimal ones at that time. More accurate results can be obtained
with smaller markers. The other possible sources of error besides marker visibility
171
are numerical approximations and small deflections near the root of the wing. This
method has been vetted by comparing known material properties and actual results
[105], therefore the approximations are reasonable. Given the aspect ratio of the
wing, there might also be some error associated to the assumptions of the beam
equations.
Alternatively one can construct a CLPT model to try and estimate the bending
stiffness. Using a symmetric layup of perfectly bonded layers one can calculate
reduced stiffness matrices based on the assumption of planar stress. These reduced
stiffness matrices can be calculated for each layer, considering their orientation with
respect to the axial or longitudinal axis. The symmetric layering assumption can be
applied to chord wise sections of the airfoil, or panels. We can identify an effective
set of geometric properties that produce the required bending stiffness. For an
orthotropic laminate the reduced stiffness matrices in the material axes are shown
in Equation (4.23).
Qmat =








E1
1−υ12υ21
υ12E1
1−υ12υ21
0
υ12E1
1−υ12υ21
E2
1−υ12υ21
0
0 0 G12








(4.23)
The following layup was used in the calculation of the rotated stiffness matri-
ces: Layup = [−90,+90,+45,−45, 0, 0,−45,+45,+90,−90] and the rotation matrix
172
used is shown in Equation (4.24).
T =








cos2θ sin2θ −2cosθsinθ
sin2θ cos2θ 2cosθsinθ
cosθsinθ −cosθsinθ cos2θ − sin2θ








(4.24)
Using the available information on the layer thickness of each component ma-
terial, as well as the material properties specified for unidirectional fibers, placed at
the corresponding angles. Once these stiffness matrices have been constructed, one
can evaluate the generalized stiffness A, cross-coupling B and bending stiffness D
matrices. The bending stiffness can be calculated from Equation (4.25).
EI = c ∗D (1, 1) (4.25)
Where c is the airfoil chord and D(1,1) is the first entry of the resulting
bending stiffness matrix. The calculated bending stiffness for the airfoil is then
EI = 25.4949Nm2 .
This value is calculated based on a number of assumptions, however, the layup
and the resultant effective thickness provides a calculation which is close to the
measured results. Sources of error in this calculation are volume fractions of fibers
to epoxy, possible out of plane orthotropic properties and airfoil shape approximate
panel partitions. A better analytic calculation can be obtained from other methods.
173
4.4 Wind tunnel tests
The wing section described in Section (3.6.1.1) was instrumented with strain
gauges and tested at the REEF wind tunnel of the University of Florida [107].
The sensor locations chosen were according to Section (3.6) for sensing vertical
gusts. The wing sensor placement considerations were mainly for bending, however
with the goal of extracting forces and moments by means of system identification
tests done apriori.
4.4.1 Experimental setup
Fifteen locations were chosen for installing full bridges. The sensitivity of a full
bridge is the best as compared to half a bridge or quarter bridge. Each Wheatstone
bridge was set in bending configuration, that is a pair in tension and a pair in
compression. The sensor network was amplified using a 5V voltage source and a
National Instruments USB DAQ mass terminal was used to acquire the signals. The
sensors had a nominal resistance of 350 Ohms each and have a temperature range
of -100F to 350. The gauge factor as typical of constantan foil gauges on polyimide
film was of 2.160. The wiring consisted of small gauge wire with low -weight cable
jackets. A greater resolution is found by increasing the supply voltage to 10V,
however increasing the risk for overheating and gauge delamination.
The wing had to be prepared for sensor installation by means of surfactants
and abrassive removal of the superficial epoxy layer. Better adhesion and bond
endurance is achieved when the gauges are bonded to textured surfaces. Soldering
174
Figure 4.34: Strain gauge top half of full bridge
pads were used to add stability to the installations and prevent detachment of the
gauges. The cables were carefully placed to try and avoid tripping the boundary
layer. Although this may not be entirely the case, the streamlined configuration
used is much better than random orientation of the connections. Servo motors were
installed in the two ailerons.
Figure 4.35: Wing with sensor network
An amplification PCB board was designed to reduce the cable volume Fig-
175
(a) PCB setup (b) Prototype board setup
Figure 4.36: Strain gauge amplification electronics
ure 4.36. Initial tests however were performed with the prototyped amplification
board. The hardware preparation consisted of tuning each sensor bridge with two
trimming resistors. These produced a desired bridge output in the available voltage
output range so that it would be evenly split for positive and negative directions of
wing bending. Each bridge output was then connected to the PCB board with gain
resistors which were tuned to yield a gain of approximately 100 times the incoming
signals.
The data was then collected in labview for the open loop tests as well as the
closed loop tests. For the closed loop tests the weighting matrices were coded in fast
loops achieved with a compromise of sampling rate and samples sent to the ailerons.
After tuning this loop the latency was minimized, albeit non comparable to an on-
board processor with hard coded algorithms. Therefore the bandwidth seen with
the labview DAQ setup for data sensor acquisition and actuator signal generation
can be improved. For the closed loop both actuators were used simultaneously. The
176
tests entailed steady state as well as turbulence tests by means of an ATG or active
turbulence grid, as developed by [108].
The wind tunnel is a low speed wind tunnel of open section type. Velocities
can go up to 20 m/s. The overall length of the wind tunnel is approximately 45
feet while the inlet is a 10 by 10 ft area. The enclosure can hold large aerodynamic
models. The commanded throttle is controlled via an analog voltage 0-10V and the
airspeed inside the tunnel is monitored by a pitot probe.
Figure 4.37: REEF wind tunnel
The wing was set about 68 inches downstream from the wind tunnel inlet. It
was mounted on a six axis load cell Figure 4.39 and had a plywood floor which went
from the inlet up until the location of the load cell. The mounting fixture consisted
of an assembly of aluminum plates and c-channels. A load cell plate was adapted to
the load cell bolts directly. Attached to this plate a wing adapter plate was fixed by
means of four bolts. The wing adapter plate had two c-channels for wing positioning
over the plate. Initially one c-channel was used, however it was observed that there
177
were some oscillations induced by the lack of stability on the opposite end of the
fixture. Therefore the remaining tests were done using the additional c-channel
which was mounted on the opposite end of the first one. On top of the c-channel
there was the conformal reinforced clamps which clamped the wing. These clamps
too were bolted together to firmly grasp the wing, however without exerting too
much force and cause damage to the wing. The wing setup with the fixtures and
the load cell is shown in Figure 4.38.
The instrumentation used consisted of a voltage supply, amplifying electronics,
a DAQ board and a DAQ computer. This equipment was set close to the wing. The
ATG was rolled over the inlet for the turbulence tests, as shown in Figure 4.40.
Figure 4.38: Wind tunnel wing setup
178
Figure 4.39: Load cell axes
For each of the tests described in Table, Figure ?? two sets of data were
collected: first the data coming in from the wind tunnel and the data being recorded
from the strain gauge sensors installed on the wing.
Table 4.4: Overview of initial wind tunnel tests
Steady State: (2 Angles of Attack and 3 Speeds):
AoA 0 deg: 5m/s, 7.716 m/s, 10.288 m/s, 12.861 m/s (4x2)
AoA 5 deg: 5m/s, 7.716 m/s, 10.288 m/s, 12.861 m/s (4x2)
AoA 0 deg, Ramp: 10.288 m/s, 12.861 m/s, 7.716 m/s, 10.288 m/s (x2)
AoA 5 deg, Ramp: 10.288 m/s, 12.861 m/s, 7.716 m/s, 10.288 m/s (x2)
Turbulence Grid (AoA: 5 deg)
Turbulence Grid (ATG): 4.73 m/s and 7.11 m/s (x2)
Turbulence Grid (ATG) at 7.11 m/s:
Open loop (x1)
Closed Loop using S1.S15, One Aileron (x15)
Closed Loop using S1S15, Two Ailerons (x15)
Open Vanes (x2)
179
Figure 4.40: Test setup with ATG
The load cell data was collected at 2000 Hz and the strain gauge data was
collected at 10kHz. For preliminary testing two DAQs were used, the load cell data
was recorded using a NI-PXI system available at REEF whereas for the sensor data,
a USB mass terminal DAQ was used. In order to synchronize the measurements
coming out from the two DAQs a mechanical pulse was used, this is shown in
Figure 4.42.
The voltage data was passed through a calibration matrix M for the JR3 load
cell which was used so that the forces and moments produced were calculated in
Newtons (Forces )and Newtons.meters (Moments), as F = MV . V are the available
voltage values per channel. Each of the tests were approximately 2 minutes long,
with a steady velocity approximate duration of 30 seconds. This area of steady
flow for the steady tests can be seen from the flat portion of the signals. . In
180
Figure 4.41: 6-axis sample data from load cell during steady test, mechanical pulses
were used to synchronize data
addition to the conversion to forces and moments, each voltage channel had an
offset which translated to an offset in forces and moments. These offsets were due
to gravity effects in the balance of the heavy aluminum fixture assembly. The JR3
load cell however had sufficient load capacity. These offsets were subtracted from
the measurements to obtain the actual force and moment variations.
4.4.2 Experimental results
The frequency content of the wind tunnel data in both the steady and the
turbulent tests can be analyzed in the perspective of vibration modes coming both
from the wing structure, fixture structure and aerodynamic effects. It was observed
that the oscillations of aerodynamic behavior did not produce divergent aeroelastic
effects. It could also be observed that the Lift channel had the strongest impact
of aeroelastic interactions of the wing. In order to gauge for proper airfoil trends,
181
i.e. lift increases with speed and angle of attack, the load cell data vs. steady wind
tunnel speed is presented in Figure 4.42.
The lift coefficient was calculated for each case using the four speeds: 5 m/s
(9.72 kts), 7.716 m/s (15 kts), 10.288 m/s (20 kts), 12.861 m/s (25 kts). The values
in Figure 4.43 were calculated using a wing length of 0.555 m and air density of 1.2
Kg/m3 at 20 Celsius (68 F). The Reynolds number was calculated using the chord
as characteristic length and a dynamic viscosity of 1.8205E-5 Kg/m.s .
These results can be compared to the AG35 polar plots found in Figure 4.44
and the AG35 rotated airfoil, found in [109].
Even though the choice for a rotated airfoil to compare the results for CL and
CD needs to be further verified, the rotation by 1.563 degrees to orient the chord
line horizontally was very similar to what the UMD Pwings two point mount did.
In addition to this, the values for CL and CD are in close proximity with those
predicted by the low speed airfoil data of the rotated airfoil.
From the Lift coefficient graph in Figure 4.44, a zero degree angle of attack
will produce a lift coefficient close to 0.2 and for a 5 degree AoA a lift coefficient
close to 0.54. These values are in close agreement with experimental results. A
similar comparison can be drawn for the drag coefficient, albeit the measured values
seem to be in the same order of magnitude but higher than the results of the polar
curve in [109].
A sustained gust test was performed to test how well the strain gauge sensors
tracked the wind speed changes. Figure 4.45 shows the forces and moments as well
as the strain signals coming from each bridge.
182
5 6 7 8 9 10 11 12 13
0
0.5
1
1.5
2
2.5
3
3.5
Wind tunnel speed(m/s)
Fo
r
c
e
s
 
a
n
d
 
Mo
m
e
n
t
s
Forces and Moments vs Wind Tunnel Speed
 
 
channel 1 (Lift)
channel 2 (Drag)
channel 3 (Sideforce)
channel 4 (Yawing Moment N)
channel 5 (Rolling Moment L)
channel 6 (Pitching Moment M)
(a) Forces and moments AoA=0 deg
5 6 7 8 9 10 11 12 13
0
1
2
3
4
5
6
7
Wind tunnel speed(m/s)
Fo
r
c
e
s
 
a
n
d
 
Mo
m
e
n
t
s
Forces and Moments vs Wind Tunnel Speed
 
 
channel 1 (Lift)
channel 2 (Drag)
channel 3 (Sideforce)
channel 4 (Yawing Moment N)
channel 5 (Rolling Moment L)
channel 6 (Pitching Moment M)
(b) Forces and moments AoA=5 deg
Figure 4.42: Calculated forces and moments vs wind tunnel speed for angle of attack)
degrees and angle of attack 5 degrees
183
Figure 4.43: Calculated lift and drag coefficients from wind tunnel data
The rms level data collected from each bridge for each of the steady velocity
tests was plotted in a pattern that conserves the distances between each bridge
pair along the wing span. Interpolation of the bridge data renders a 3D surface of
measured rms strains on the surface of the wing. Figure 4.46 shows a representative
case where the higher strains are close to the root as expected. This map can be
better built if more sensors are available. During the preliminary tests four sensors
did not work properly and therefore had to be taken out of the set. If access is
possible to more sensor locations, then the following figure would have more bridge
points (hollow circles), as shown in Figure 4.46.
If a top view of these RMS calculated strains is taken, a similar pattern is
shown for each case of increased angle of attack and increased wind velocity. These
patterns are interpolations based on experimental measurements. The strains were
calculated by means of the gauge factor and the bridge voltage. The 0 deg and 5
deg angle of attack cases are shown in Figure 4.47.
The maximum rms strains also increase with wind speed, as can be observed
from a comparison of 3D interpolated rms strains at a zero angle of attack and low
184
(a) Rotated airfoil
(b) Lift/Drag coefficients vs angle of attack
Figure 4.44: AG35-r Drela aifoil, from [109], rotated airfoil and lift/drag coefficients
vs angle of attack.
185
(a) Forces and moments
(b) Strain sensor data
Figure 4.45: Ramp test to simulate sustained gusts (20 kts, 25 kts, 15 kts, 20 kts),
AoA =5 deg: load cell forces and moments and strain sensor data.
186
Figure 4.46: Representative of interpolated strain rms measurements for steady wind
velocity
speed, with a 5 degree angle of attack and fastest speed. This comparison is seen
from Figure 4.48, where the maximum strain goes from about 4.2 e−4 to about 6
e−4.
Wing tip static loading tests were also performed with tip loadings of 100 gr
to 800 gr. Two of the tests are shown in the following graphs of Figure 4.50. It
is interesting to observe that although strains are higher at the root some of the
sensors contribute more in the steady wind dynamic case. This pattern is consistent
with what is observed from the mass weighted strain simulation for the full wing.
Turbulence tests were carried by means of an active turbulence grid or ATG.
It consists of two orthogonal grids of rotating vanes. Each vane consists of a shaft
with attached winglets which are driven to rotate by a computer controllerd motor.
This device is placed in the wind tunnel at the upstream end of the test section
and the action of the spinning vanes causes the air flowing throught the device to
187
(a) Interpolated strains at AoA 0 deg
(b) Interpolated strains at AoA 5 deg
Figure 4.47: Interpolated strain rms measurements for steady wind velocity at AoA
0 deg and AoA 5 deg
188
Figure 4.48: Comparison of max. rms strains during steady runs
Figure 4.49: Comparison of strains produced with wing tip loadings of 100gr and
200gr
become highly turbulent.
Tests with the ATG grid on were carried using a core velocity of approx. 7.1
m/s. The random processes associated to this velocity have been characterized
with a hot-wire probe in [1]. The percentage variation from a mean core velocity
are attributed to external airflow events and do not impact the quality of the flow.
Overall the baseline flow is very uniform. The ATG was shown capable of generating
free stream turbulence intensities from 15 % to 35 %, with length scales on the
189
order of 0.6 m and with spectra approximating the von Karman model that showed
a possibility for tonal or smooth behavior depending on run mode. The statistical
characteristics of the ATG are given by the tables corresponding to x/M=12.5 and
43.5 % throttle. The statistical tables for the ATG are found in excerpts from [1] in
Appendix B.
The turbulent data can only be analyzed with statistical metrics, given that
this is a random process. The following graph shows the calculated Lift from load
cell data without strain feedback and with strain feedback. The strain feedback gain
used was K=1, both with one aileron and with two ailerons for feedback. Higher
gains could not be tested due to experimental constraints in that the wind tunnel
was available for limited amount of time.
Figure 4.50: ATG produced Lift force with and without strain feedback
190
One way to quantify improvements of using strain feedback is to calculate the
rms level of the measured lift force. Table 4.5 shows results of using strain feedback
with a gain of K=1, feeding back from different bridges individually, where bridge
1 is S1, bridge 8 is S8, etc All the ATG tests were about 2 minutes on average,
sampling at 2000 Hz.
Table 4.5: RMS lift using strain feedback at different strain gauge bridges
Test Case RMS level Lift Force (Newtons)
No control, turbulent Lift 12.9535
Strain feedback using one aileron SFB at S1 12.6901
Strain feedback using two ailerons at S1 12.0564
Strain feedback using two ailerons at S8 11.2332
Strain feedback using two ailerons at S9 11.1297
Strain feedback using two ailerons at S12 11.0696
Control Lift RMS two ailerons at S14 11.2330
4.4.3 System identification of output matrix
In the perspective of a bio inspired approach and to test the distributed sensing
concept, we used system identification techniques to relate forces and moments to
sensor measurements. If the set of measurements over a number of distributed
sensors are weighted in the least squares sense they can be used for feedback control.
The matrices can anticipate a force/moment measurement and enable the controller
191
to calculate an actuator command to anticipate impending rigid body motion. It
has been shown in previous work by the UMD AVL lab that the weighting algorithm
provides an improvement in bandwidth and allows for noisy sensor signals, if the
noise is non-correlated [9]. Other techniques can be used to overcome this type of
sensor noise.
For the system identification we used three methods to generate the loads and
calculate output matrices. The output matrices correlate the forces and moments
at the root of the wing -and therefore at the fuselage if there was one attached to
the wing- with data from sensors which have been distributed over the surface of
the wing. The following tests were carried for the system identification of output
matrices.
The first test consisted of a static loading in the Lift direction. The wing
was loaded at the wingtip and chordwise at the shear center by means of a pulley,
string, and static weights. This was done at only zero angle of attack. The second
test consisted of wind tunnel loading at different speeds and angles of attack. In
this case the wing was uniformly loaded by the air stream at different velocities.
Ten velocities and four angles of attack were used. In the third approach, static
manual loading was exerted on the wing to produce coupled and decoupled forces
and moments. The wing was randomly loaded to try and produce coupled and de-
coupled loads. Some loading directions were more challenging to obtain as can be
expected, in particular the yawing direction and all the moments.
A second setup consisted of discrete gusts to test for strain feedback proper-
ties and flight envelope increase under near stall conditions. The main UMD goal
192
however was to test the effectiveness of the strain feedback controller in mitigating
the gust disturbance.
The discrete gust test was designed to provide the wing with additional lift for
a different duration of time. The gust generator (Figure 4.51) was designed and built
by Aurora Flight Sciences [7] and it consisted of tubing, a fan and a latch mechanism
to control the duration of the gusts injected into the air stream. The gust generator
was placed upstream of the test section and it was adjusted in height such that it
would produce an increase in lift that would make the wing stall. This additional
objective had the purpose of demonstrating that the strain feedback was able to
increase the flight envelope of the wing. This is an added benefit and application
of the proposed sensing mechanism. For this setup the plywood tunnel floor was
removed and the amplification electronics were placed underneath the wing.
Three wings were tested in this setup as well. The fiber optics wing, the
stanford wing with a stretchable sensor network and the off the shelf UMD wing.
The stanford wing setup and amplification board are shown in Figure 4.52.
depicts a sample static forcing case to find the output matrices with coupled
components. The coupled forces should yield matrices which are non-sparse and
with coupled entries. The initial portion of the graph has small variations which
are also encoded in the strain signals, although the graph doesn’t show the scale of
these variations.
Another set of tests entailed forcing the wing with static directional forces with
the purpose of decoupling the forces and moments. Figure 4.54 shows the case for
Fx, where the the lift direction or the x direction as well as the rolling moment are
193
Figure 4.51: Setup for Static system identification in the lift direction
very pronounced. Figure ?? shows results for the data in Figure 4.54, Figure 4.55,
and Figure 4.56.
The least squares estimator used a linear regression to calculate the dependen-
cies of the overal strain field measured by the 15 sensing bridges and the forces and
moments at the root of the wing. The Cmatrices for manual Fx, Fy, Fz calculated
with least squares estimator and colored noise. These are shown in Figure 4.57.
For the gust tests, the following matrices in Figure 4.58 were calculated using
least squares for the two particular design points which were described earlier. These
matrices may be slightly different from those of the static case, due to possible
excitation of modes that were not excited in the static tests.
194
(a) Stanford wing with stretchable
sensor network
(b) Stanford wing amplification board
Figure 4.52: Stanford wing wind tunnel gust setup
4.4.4 Discrete gust tests
In order to measure the impact of gusts over the wing, open loop and closed
loop tests were carried at two specific design points. These design points were
originally intended to drive the wing close to a stall at the moment the gusts were
applied. However for the UMD wing the stall was not easy to set up or observe
therefore the focus will be in an overall reduction of the impact of gusts.
For the case of 3 deg angle of attack and a wind speed of 20kts, the following
impact on the wing forces and moments was observed, as shown in Figure 4.59.
As it can be observed from Figure 4.59, the impact of the upstream gusts on
the forces is mainly observed with a shift from the steady state forces, along some
oscillations due to aerodynamic and structural effects. Feedback was implemented
195
Figure 4.53: Static manual forcing with coupled components
Figure 4.54: Static manual decoupled forces and moments Fx, forces and moments
(left), sensor data (right)
by means of three methods: 1.Using only one sensor location (S2) for feedback,2.
Using an average Cmatrix, which kept an overall average of all the sensors at all
times, 3.Using system identified output matrices. The results of the closed loop
tests are shown in Figure 4.60, Figure 4.61 and Figure 4.62. For the feedback based
on a single sensor location (Figure 4.60), a sensor near the root of the wing and far
from the aileron was selected. This method shows that it produces some correction
to the overall force shift due to the gusts. However it also shows that is prone to
more oscillations due to sensor noise. The feedback which uses an average Cmatrix
196
Figure 4.55: Static manual decoupled forces and moments Fy, forces and moments
(left), sensor data (right)
Figure 4.56: Static manual decoupled forces and moments Fz, forces and moments
(left), sensor data (right)
to produce aileron commands, shows that it is able to correct for the overall shift
and reduce the amplitude of the oscillations observed when the gusts are applied.
The response shows that there is a good level of oscillation. Finally, with the sys
Id matrix the feedback is able to reduce the shift as well as reduce the amplitude of
the oscillations while the gust is acting on the wing.
A moving average was used for generating Figure 4.63 and allow for evaluation
of the improvement in counteracting the gust forces. From Figure 4.63 it can be
197
Figure 4.57: Cmatrices for manual Fx, Fy, Fz calculated with least squares estimator
and colored noise
Figure 4.58: Cmatrices for design points: AoA 9 deg, V=15kts (left) and AoA 3 deg,
V=20 kts (right).
observed that the overall shift from steady flight is reduced with strain feedback,
which was applied using three methods. The red curve represents the impact of gusts
on Fy without any feedback. The blue curve represents the correction achieved using
only S2 for feedback, which yields a good correction, although with some oscillation.
The magenta curve was obtained using a C average matrix, again the oscillation is
observed mostly during the application of the gusts. The green curve was obtained
by using the sys id matrix. It can be seen that the sys id matrix improves both
the shift and the amplitude of oscillations during the application of the gust. An
increased bandwidth can reduce the shift even further. Other methods can be
explored in order to continue to measure the improvement achieved when using
198
Figure 4.59: AoA=3 deg, V=20 kts, gusts, open loop
strain feedback and in particular the sensor to force and moment output matrices.
4.5 Flight testing
A similar wing to that of the wind tunnel tests was instrumented with strain
gauges. The locations of the gauges were distributed based on MPF-v calculations
from Figure 3.43. The outermost locations were chosen farther than the predicted
areas of larger strains for observation of strain phenomena near the wing tips. The
full wing was then assembled on top of the fuselage, with minor modifications to
the original Apprentice. The modified aircraft with the strain gauge instrumented
wing is shown in Figure 4.66.
Analog amplifiers were connected to the strain gauge bridges for signal output
to the xbee transmitter. Signals were transmitted and sent off to a ground station
199
Figure 4.60: AoA=3 deg, V=20 kts, gusts, closed loop feedback S2
where they were also recorded. Signals from the strain gauges were also used in
conjunction with the calculated Cmatrix for roll stabilization. This process is shown
in Figure 4.67. During the flight trials, an autopilot was also tested, the APM 2.6.
This setup will however require a multiplexer to input all the strain signals into the
APM for processing and output to the ailerons and other controls surfaces. Two
lipo batteries were used. One to power the ESC to the propeller throttle, or the
main power, and another one to power the strain gauge amplification board and
processor on board.
The voltage signals were obtained by the ground station and seemed to follow
the stages of the flight in open loop. This is shown in Figure 4.68, where the ground
data is clearly distinct from the flight data. In the flight data portion, several
observations can be made. First, the level of strain from each sensor bridge is
200
Figure 4.61: AoA=3 deg, V=20 kts, gusts, closed loop feedback C average
expected. The bridges closer to the wing tips, 1 and 8 have less strain response than
those closer to the root of the wing, 3,4,5,6. Strain gauges 2 and 7 have less of a
response as well although higher than than 1 and 8, also expected. Some channels
were filtered and others were not, this allows for observation of the actual strain
data being collected and the processing required.
After a few open loop flights, strain feedback was implemented with the out-
put matrix from the physics based approach. This matrix had proportionality of
expected forces and moments as well as strain signal intensity. The open and closed
loop sections chosen correspond to those where the strain differentials have similar
values. The open loop data and corresponding roll rate is shown in Figure 4.69 and
the closed loop data is shown in Figure 4.70.
Overall, it is noticeable that the roll rate decreases in magnitude when the
201
Figure 4.62: AoA=3 deg, V=20 kts, gusts, closed loop feedback C matrix (sysID)
closed loop strain feedback is enabled. A comparison of these for a single roll angle
oscillation is found in Figure 4.71.
202
Figure 4.63: Filtered AoA 3 deg, V=20 Fy data to measure performance of feedback
loop, d represents a decrease in the shift for the first gust which is pro-
portional to the shift observed in the other gusts
Figure 4.64: Apprentice aircraft with modified instrumented wing
203
Figure 4.65: Strain signals and IMU recording and feedback for flight tests.
Figure 4.66: Apprentice aircraft and ground station during flight test preparation
204
Figure 4.67: Apprentice flight tests
Figure 4.68: Open loop flight data, eight bridges, 32 sensors
205
Figure 4.69: Flight data open loop, differential strain for sensors 4 and 5 and aircraft
roll rate
Figure 4.70: Flight data closed loop, differential strain for sensors 4 and 5 and aircraft
roll rate
206
Figure 4.71: Flight data closed loop, differential strain for sensors 4 and 5 and aircraft
roll rate for flight section.
207
Chapter 5: Summary and Conclusions
5.1 Summary
Bioinspired mechanisms of force/strain feedback were investigated and applied
to small unmanned aircraft platforms. Regulation of the differential disturbance
force was used in the high bandwidth inner loop for gust mitigation.
Using a small glider UAS, it was verified that the structural time reaction scale
is much smaller in magnitude than that of rigid body motion, as it was observed
in footage of the glider’s flight. The rigid body motion lagged by roughly 0.03
seconds but was enough time for the strain sensors to activate the servo motor and
react to the sensed gust disturbance before it actually turned the fuselage. The
glider was instrumented with strain gauges [7] at the root of each wing and was
successfully flown in closed strain feedback loop. Therefore, strain/force feedback
was demonstrated using a small UAV under a sharp gust condition. The gust which
was used for the experimental tests resembled a 1-cosine gust. Strain sensing output
was successfully used in a fast inner loop for reflexive reaction to gusts. Aileron
dynamics were fast enough for the employed combination of gust speed and roll
correction.
Linear strain measurement weightings or wide field integrated strain can pro-
208
vide disturbance force state estimations which can be used in these fast inner loops.
It was observed that this feedback is much more effective at higher proportional
gains in the case of a strain/force proportional controller. Less deviation in the roll
trajectories was observed by means of the D2 metric as the gain was increased. The
structural feedback controller also demonstrated to be more effective than classical
rigid body dynamics controllers, as it was compared using vicon and xbee transmis-
sion as the IMU.
A quadrotor platform was also used to demonstrate strain feedback. One
of the arms of the quadrotor was instrumented with strain gauges, and was used
in closed loop to correct for a heave gust condition. The heave disturbance was
imparted by artificially injecting a motor disturbance and after the disturbance had
been generated, the strain feedback compensator was activated. Material creep
and therefore sensor creep was accounted for by carrying a parallel strain signal
to which the controller resorted to when the disturbance was activated. In this
mechanism, there can be sensor drift to a certain extent, as the loop is constantly
updating the measurement of level flight. This feedback is much more effective at
higher proportional gains in the case of a strain/force proportional controller. Less
deviation in the roll trajectories was observed by means of the L1 metric as the gain
was increased. Motor dynamics were fast enough for the employed combination
of artificial gust and heave correction. A robust controller was implemented using
strain feedback on the quadrotor, albeit in the case of the strain controller it was
only for one degree of freedom [6].
In the perspective of all the aeroelastic models, static and dynamic, three
209
methods for strain feedback compensation were developed. These methods corre-
spond to what would be more reasonable given the available degrees of freedom for
sensor implementation.
First, a 1D static model is proposed for computing spatial weightings to trans-
late strains into forces and moments. In this approach, the strain feedback provides
a measure of the force imbalance on the lifting surfaces of the aircraft, and therefore
proprioceptive. This approach can be implemented either with a single sensor on
each control surface or it can be implemented using several sensors. One way of in-
tegrating the multiple sensor outputs is through wide field integration, in the sense
that each sensor is weighted according to its expected physical contribution to the
system in stabilizing a particular force. In this way, a one dimensional roll controller
based on outputs weighted by physical constants was developed and was tested in
a full 6-DOF dynamic linear aeroelastic simulation (ASWING). The controller roll
tracking response was good, although the speed of response can be improved by
adding outter loop stabilization. The principle however was demonstrated by clos-
ing the inner loop with strain feedback. The distributed 1DOF mechanism can
also be used to detect disturbance propagation as insect wings such as the ones in
drosophila can achieve by detecting torsional waves from the moment they begin at
the wing tip onward.
Second, a wings only approach is presented as well, where a distributed sensor
network is empirically correlated to the forces and moments at the root of the
airfoil. This procedure was demonstrated by instrumenting a carbon fiber wing layup
manufactured by [102] and implementing force feedback by means of a pseudoinverse
210
matrix, which had been calculated by a least squares regression. A closed loop
improvement was observed for cases of discrete gusts of different durations as well
as an improvement in the case of the turbulent wind tunnel tests.
Third, with a system identified strain map vs forces and moments at the
root of the wing, estimations of the differential six degree of freedom forces can be
made. This approach basically means that once the wing has been identified, all
that needs to be done for control is to regulate the differences between the forces
at the root of the wings. This approach would only need the main lifting wings as
sensory surfaces for disturbance rejection. It is also a proprioceptive approach in
that relative measures of detected forces and moments on each wing are constantly
being compared and when any difference is detected it runs in the high bandwidth
inner loop which may be robust as the calculated forces and moments are 6-DOF.
A Full aero elastic model can be developed by measuring aeroelastic stability
derivatives through system identification or by aeroelastic simulation packages such
as NASTRAN or COMSOL, with dedicated computational resources.
A framework for sensor placement using an FEM model of the wing was de-
veloped for open loop sensor placement. The observability map, which correlates
well with the MPF-v weighted map that was obtained, relies only on the material
properties of the wing. Closed loop sensor placement mechanisms can be used when
a particular application is being targeted. Actual wing properties can be used to
tune the model. Experimental measurements of static and dynamic wing properties
were taken and it was found that the torsional stiffness of the wing is much larger
than the bending stiffness. This behavior however is expected from the +- 45 degree
211
carbon fiber layout. Some error in the vicon estimation is made due to corrections
in the Euler Bernoulli equations for a plate. A 2D wide field integration approach is
presented for mapping of the impact of gusts on wing strain by means of the MPF
weighted eigenstrains. A set of gains can be retrieved from the corresponding MPF
weighted patterns and can be applied to any number of sets of sensors that are
physically realizable. An integration of all the basis functions yield sensor locations
which are sensitive to all the different types of gusts, that is, in the corresponding
degree of freedom.
A CFD/structures gust framework is proposed to study the effect of gusts on
wing structures and help in validation of sensor locations. The length of the control
volume was not as long as it should be due to computational time constraints,
however, it was shown that this is a useful tool to study the static aeroelastic
interactions of fluid flow and wing structures. It was observed that the principal
strains change during the interaction with the gusts.
With the wind tunnel tests it was also shown that the observability Grammian
results for 2D sensor placement was coherent with the finite element simulation re-
sults. The rms strains showed to have energy patterns consistent with the mass
participation factors. As extrapolated from the 1D analysis where the effective
independence metric is proportional to the observability Grammian, these results
produced similar to those produced by eigenstrain weighting using directional mass
participation factors. This weighting produced best sensor location regions for best
observation of vertical forces on the wing and therefore vertical gusts and perhaps
lateral gusts if combined with another metric. Three estimation matrices were con-
212
structed with the linear regressors, they all produced improvements in the closed
loop measured lift in the presence of gust disturbances.
It was also observed that wind tunnel load cell data contains frequencies cor-
responding to structural frequencies, aerodynamic frequencies, and noise. Sensor
data also contains frequencies, although much less than the load cell data. Sensor
data was shown to be consistent with open loop sensor placement. ATG data shows
that strain feedback allows for a reduction in the average values of Lift force, at a
proportional feedback gain of 1. ATG data shows that using two ailerons for struc-
tural feedback produce better results in reducing the average values of Lift force.
The controller was only tested at K=1 but it produced a 10 % improvement over
the rms lift disturbance produced by the rotating vanes upstream of the wind tunnel
test section.
A flight test served to show that these principles work in actual flight, with the
uncertainty caveats, as small UAS aircraft materials undergo thermal and structural
changes. The sensor placement results depend on boundary conditions imposed on
the sparse equation solver. Therefore changes in these boundary conditions will
somewhat change the estimates for best observability. The extent to which these
changes dramatically affect the performance of the closed loop controller during
flight is a quantity that can be further studied.
213
5.2 Conclusions
From experimental work with four platforms as well as simulations of proposed
mathematical mechanisms, it was shown that strain sensors can be used in fast
inner loops for reflexive reaction to gusts. The mechanisms chosen in this work were
bioinspired in proprioceptive qualities as well as distributed sensing features.
Structural feedback provides faster reaction times than IMUs alone, meaning
that there indeed is an improvement by using signals coming from material defor-
mations of the aircraft wings.
Strain feedback can be implemented in different types of small unmanned aerial
platforms: glider, quadrotor, wind tunnel wing, flight test wing. It was effective in
rejecting discrete roll disturbances of the glider, lift disturbances of the wind tunnel
wing, roll disturbances of the flight wing, and heave disturbances of a quad rotor.
In the case of the quadrotor, these were injected at the rotors to simulate the effect
of the aerodynamic disturbances. For all the platforms, there is a certain degree of
flexibility to obtain an adequate signal to noise ratio. In particular for the quadrotor,
quad arms had to be designed to obtain a good signal to noise ratio without loosing
thrust power.
Effectiveness of the strain feedback mechanism heavily depends on the servo
motor response, i.e. time constant. An ideal case would have the servo response
be as fast as the strain sensor response, although typically, for off the shelf strain
sensors this response is much higher by at least one order of magnitude.
Sensor noise is mitigated by using the least squares estimator. This estimator
214
considers that the noise has a zero mean and is uncorrelated. When the pseudo
inverse is calculated it will correct for the noise in the sensor readings to the extent
of the number of sensors. Sensor noise however still affects the measurements,
especially if a reduced number of sensors are used. Analog filters and other faster
response filters can be implemented to provide improvements.
Mass participation factors for sensor placement lead to an estimation of the
observability Grammian in open loop, in the case of vertical gusts. This Grammian
is directly related to the effective independence index-EFI, which is an engineering
approach to mode targeting. It was observed that the MPF method produced more
restricted sensor placement areas than the EFI approach. This analysis was carried
in 1D and 2D. Wind tunnel experiments showed agreement of the strain sensor data
with the areas predicted by the MPF map to be of higher rms output.
Strain sensor implementation at the selected locations based on placement
analyses had good stability characteristics and was temperature compensated in
a full bridge bending configuration. Gauge wiring seemed to be strong yet unob-
trusive to the strain measurements. For the flight wing additional encapsulation
options were pursued to avoid strong shear over the sensors and prevent damage
to the sensors and wiring. Amplification and connections were customized for fast
troubleshooting.
Different types of C-matrices were calculated. One proposed method entails
an output matrix based on the physics of the different aircraft dynamics, i.e. roll.
Another one entails output matrices that were system identified according to partic-
ular flight conditions, both static and dynamic. MPF maps can be used to construct
215
spatial maps which yield areas of more or less influence of the different types of gusts.
Theoretical considerations produced a 6-DOF model for disturbance mitigation on
the assumption that the output matrix which relates the forces and moments at
the root of the wing has been identified. A 6-DOF decoupled model is also possi-
ble if considering differential balancing of strains on corresponding controls surfaces.
These differentials form output matrices to be used in SISO and MIMO strain based
controllers.
Closed loop linear aeroelastic simulations using both strain feedback and bioin-
spired Cmatrix based on the physics of the rolling motion, showed correction of the
roll angle as well as of the overall trajectory.
Strain feedback was effective in reducing the peak rms lift disturbance pro-
duced by discrete gusts of different durations in wind tunnel tests. It also showed
a promising correction to turbulent air with a core velocity of 7.1 m/s. Different C
matrices were used to estimate the lift force at the root of the wing and feed the
correction to the wing ailerons. All the matrices produced disturbance mitigation.
During flight tests, all the installed sensors showed a good response to different
flight characteristics. A roll rate reduction during closed loop was observed as well.
This would mean that the strain feedback was acting on the mitigation of roll
disturbances. Further tests will be pursued where vertical force rms and lateral
force rms will be measured for open and closed strain feedback loop. These are
similar tests to those taken in the wind tunnel with the wing section.
A robust controller framework was formulated for output matrices with multi-
plicative uncertainty. The least squares estimator will give structured uncertainties,
216
however a worse case scenario can be chosen for the unstructured uncertainty for-
mulation. The multi-degree of freedom case can be formulated in the forces and
moments of the aircraft be it at the center of gravity or be it at the root of each
wing for differential comparisons. In either event, the force estimations using strain
sensors can be of a distributed nature or a discrete nature. Distributed systems
benefit from least squares estimators as they provide a way to reduce noise in the
individual components to enable fast inner loops as no further measurements need
to be carried to estimate the states, in this case, force/strain differentials. The
bandwidth limit is imposed by the strain sensor bandwidth which is very high for
foil gauges, as compared to typical rigid body IMUs.
5.3 Future work
This work can be continued in many different directions as the strain platforms
can be tested with many different controllers and under many different types of
inputs.
Other modalities of strain sensing will also be explored in future work.
The CFD-structures or fluid structure interaction simulation can be extended
to a dynamic case by inputting gusts with frequency content matching von Karman
or Dryden profiles for a given altitude.
Closed loop sensor placement using minimization of the maximum singular
value of the system and comparison with results from open loop placements. This
procedure will allow for the closed loop to be optimized based on system require-
217
ments.
Implement the MIMO framework for the fixed wing system in the different
theoretical approaches that were presented in the theoretical section. Exploring
model uncertainties in the different parameters that were introduced is an important
theoretical direction. Physical implementation of these controls algorithms is also
an entire experimental area to be explored.
218
Appendix A: Strain sensors for aircraft wing applications
A number of off the shelf strain sensors can be used for aircraft applications.
The choice of sensor depends on the performance requirements and the application
goals. Different types of strain sensors can be applied to aircraft surfaces. These
main types, among others are: piezoelectric, piezoresistive, DEAP and resistive,
these are shown in Figure A.1.
Piezoelectric sensors use generated voltage which is produced on applied stress
and is mainly for dynamic measurements due to electric current leakage. They
are composed of piezoelectric materials such as Barium Titanate, cadmium sulfide,
Quartz, PZT, PVDF, PVC, PMN-PT. A common type of material for piezoelec-
tric applications is PVDF piezoelectric film due to its flexibility, high mechanical
strength and dimensional stability. It also has a larger area than resistive gauges,
which may make a difference in particular applications.
Resistive sensors or wire/metal strain gauges are typically made of constan-
tan foil, although there are many other alloys used as well. These may or may
not be temperature compensated but have a broad temperature range. They are
typically based on a polyimide carrier and have a gauge factor of about 2. The
maximum strain is not too large as it is mainly taylored to detect strains in the
219
micro-range. These need surface conditioning and are good for both dynamic and
static measurements.
Semiconductor strain gauges are typically made of Silicon and Germanium.
They have a higher unit resistance and higher strain sensitivity. The temperature
effects are greater too and typically are non-linear. These can be customized using
polymeric matrices and semiconductor fillers. These are also good for static and
dynamic measurements. While the higher unit resistance and sensitivity of semi-
conductor wafer sensors are definite advantages, their greater sensitivity to temper-
ature variations and tendency to drift are disadvantages in comparison to metallic
foil sensors.
Electroactive polymers are another type of strain sensors which can be applied
to an aircraft wing. They can be ionic or dielectric and possess greater flexibility.
These also have actuator features and are versatile materials. The fine film allows for
shape and size adaptation. The lenght is unrestricted which is a great advantage for
some applications. The gauge factor typically depends on substrate stiffness. The
maximum strain is greater than 100 %. These also allow for dynamic measurements
although the bandwidth may vary with the different types of materials.
In general, the surface where the strain gauges need to be bonded require sur-
face conditioning. This is to ensure a good physical and chemical bonding between
the two surfaces. There will be surfaces resistant to any bonding agent. For this, a
coating or abrasive method can be utilized to enable bonding. Solvent in bonding
agent needs to be non-reactive with the wing surface material. A different type of
glue may be used if this is the case. The impedance needs to match the moduli of
220
Piezoelectric film sensor  Resistive sensor  
Piezoresistive  sensor Electroactive polymer sensor  
Figure A.1: Strain sensors for aircraft wings
both materials as well as the shear modulus. The strain percentage on the surface
of the wing needs to match the maximum sensor strain or be below the maximum
value.
The installation requirements are such that for an airplane wing, the thickness
of sensors and connections needs to be non-disruptive to the airflow. The signal
conditioning and amplification is also very particular to each type of sensor instal-
lation, therefore light systems need to be kept in mind, especially for fiber optics
strain sensing platforms. The ultimate sensor characteristics will depend on aircraft
wing structural characteristics. These characteristics are compared using nominal
representative values for each in Figure A.2. Other additional factors for a success-
ful installation of strain sensors on aircraft wings include: effective encapsulation
221
for protection from humidity, electrical shielding and grounding, appropriate con-
nectors and mechanical stability. If the wires connecting the strain gauges to the
conditioning electronics are not protected against humidity a parasitic resistance
created between the wires and the substrate to which the strain gauge is glued or
between such wires. Acrylics adhesives, synthetic rubber resins, epoxies and cyano-
acrylates are all frequently employed in lamination and assembly. An overview of the
main advantages and disadvantages of these types of sensors is found in Figure A.3.
Figure A.2: Strain sensors characteristics comparison chart
222
Figure A.3: Strain sensors advantages and disadvantages chart
223
Appendix B: ATG characteristics
The active turbulence grid at the university of Florida was developed by [1]. It
consists of orthogonal rotating vanes, rotated by motors which rotate at a maximum
frequency of 18Hz.
Figure B.1: ATG turbulence grid
The statistical characteristics of the settings used for the wing strain tests are
found in Figure B.2. The ratio x/M is a key non dimensional parameter which
224
characterizes wind tunnel statistics and distance at which turbulence is developed
downstream from the inlet. It is calculated based on the distance of the wing to
the wind tunnel inlet upsteam. In this case the distance was approx. 68 inch which
divided by M, or the mesh size, gives a ratio of about 12.9, the closest chart that
describes the characteristics of the flow is the one corresponding to a ratio of 12.5.
Figure B.2: ATG turbulence statistics with wind tunnel throttle at 43.5 % and at
x/M=12.5, from [1]
Measured percent variations from core for freestream velocity and different
ATG run modes are shown in Figure B.3. This figure shows measured percent
variation from the mean core velocity close to 7.1 m/s for Uinf , run mode Ω =0
and ω = 3, as well as the measured percent variation from core for TIx at different
percentages with the baseline flow [1]. TIx is a relative turbulence parameter which
is calculated by dividing the square root of the flow velocity in the direction of the
225
freestream into the freestream velocity.
TIx =
rms(u)
Uinf
(B.1)
226
Figure B.3: ATG turbulence statistics mean core velocities: a. Mean core velocity at
7.21 m/s, b. Mean core velocity at TIx=3.9%, c. Mean core velocity at
TIx=17.9%, d. Mean core velocity at TIx=20.7%
227
Bibliography
[1] Michael J Sytsma. Effects of turbulence on fixed wing small unmanned aerial
systems. University of Florida, 2013.
[2] Michael Ol, Gregory Parker, Gregg Abate, and Johnny Evers. Flight controls
and performance challenges for mavs in complex environments. AIAA Paper,
6508, 2008.
[3] Christopher D Regan and Christine V Jutte. Survey of applications of active
control technology for gust alleviation and new challenges for lighter-weight
aircraft. Technical report, NASA/TM–2012 216008, Dryden Flight Research
Center, Edwards, California, April: 1–15, 2012.
[4] Lina M Castano, Simone Airoldi, Terrence McKenna, and J Sean Humbert.
Wing sensor placement for gust disturbance rejection.
[5] Lina M Castano, Simone Airoldi, Terrence McKenna, and J Sean Humbert.
Gust rejection using force adaptive feedback for roll. In AIAA Aviation Tech-
nology, Integration, and Operations Conference. AIAA, 2014.
[6] Gregory M Gremillion, Lina M Castano, and J Sean Humbert. Disturbance
rejection with distributed acceleration and strain sensing. In AIAA SciTech
Conference, 2015.
[7] Aurora Flight Sciences. Research and development,90 broadway suite 11 cam-
bridge, ma 02142.
http://www.aurora.aero/. Last accessed: June 22 2015.
[8] Fu-Kuo Chang. Structures and composites lab sacl.
http://structure.stanford.edu/. Last accessed: June 22 2015.
[9] James Sean Humbert and Andrew Maxwell Hyslop. Bioinspired visuomotor
convergence. Robotics, IEEE Transactions on, 26(1):121–130, 2010.
[10] Nicolas Franceschini, Franck Ruffier, and Julien Serres. A bio-inspired flying
robot sheds light on insect piloting abilities. Current Biology, 17(4):329–335,
2007.
228
[11] Graham K Taylor and Holger G Krapp. Sensory systems and flight stability:
what do insects measure and why? Advances in insect physiology, 34:231–316,
2007.
[12] A Skordos, PH Chan, JFV Vincent, and G Jeronimidis. A novel strain sensor
based on the campaniform sensillum of insects. Philosophical Transactions
of the Royal Society of London A: Mathematical, Physical and Engineering
Sciences, 360(1791):239–253, 2002.
[13] Thomas Werner, Shigeyuki Koshikawa, Thomas M Williams, and Sean B Car-
roll. Generation of a novel wing colour pattern by the wingless morphogen.
Nature, 464(7292):1143–1148, 2010.
[14] Mikko Juusola and Andrew S French. The efficiency of sensory information
coding by mechanoreceptor neurons. Neuron, 18(6):959–968, 1997.
[15] Carlo Menon, Rebecca Brodie, Sally Clift, and Julian FV Vincent. Concept
design of strain sensors inspired by campaniform sensilla. Acta Astronautica,
64(2):176–182, 2009.
[16] Rhoe A Thompson, Johnny H Evers, and Kelly C Stewart. Attitude control
augmentation using wing load sensing-a biologically motivated strategy. In
AIAA Atmospheric Flight Mechanics Conference, 2010.
[17] Joel A DeLisa, Bruce M Gans, and Nicholas E Walsh. Physical medicine
and rehabilitation: principles and practice, volume 1. Lippincott Williams &
Wilkins, 2005.
[18] Physiopedia. Proprioception.
http://www.physio-pedia.com/Proprioception. Last accessed: June 22 2015.
[19] Brian Coles. Essentials of avian medicine and surgery. John Wiley & Sons,
2008.
[20] Uwe Proske and Simon C Gandevia. The proprioceptive senses: their roles in
signaling body shape, body position and movement, and muscle force. Physi-
ological reviews, 92(4):1651–1697, 2012.
[21] Robert F Schmidt and Gerhard Thews. Human physiology. Springer, 1983.
[22] Robert J Kosinski. A literature review on reaction time. Clemson University,
10, 2008.
[23] Sharon A Plowman and Denise L Smith. Exercise physiology for health fitness
and performance. Lippincott Williams & Wilkins, 2013.
[24] Scott Kline Powers and Edward T Howley. Exercise physiology: Theory and
application to fitness and performance. 2007.
229
[25] A Mohamed, S Watkins, R Clothier, M Abdulrahim, K Massey, and R Saba-
tini. Fixed-wing mav attitude stability in atmospheric turbulencepart 2: Inves-
tigating biologically-inspired sensors. Progress in Aerospace Sciences, 71:1–13,
2014.
[26] Paul Noll. Avian muscles.
http://www.paulnoll.com/Oregon/Birds/Avian-Muscle.html. Last accessed:
June 22 2015.
[27] John J Videler. Avian flight. Oxford University Press, 2006.
[28] Mark NO Davies and Patrick Green. Perception and motor control in birds:
an ecological approach. Springer Science & Business Media, 2012.
[29] Andrei Sourakov. Faster than a flash: the fastest visual startle reflex response
is found in a long-legged fly, condylostylus sp.(dolichopodidae). Florida En-
tomologist, 94(2):367–369, 2011.
[30] Cedric Gillott. Entomology. Springer Science & Business Media, 2005.
[31] Robert B Northrop. Introduction to dynamic modeling of neuro-sensory sys-
tems. CRC Press, 2000.
[32] Vincent H Resh and Ring T Carde´. Encyclopedia of insects. Academic Press,
2009.
[33] SASHA N Zill. Plasticity and proprioception in insects. ii. modes of reflex
action of the locust metathoracic femoral chordotonal organ. Journal of ex-
perimental biology, 116(1):463–480, 1985.
[34] Brian S Beckett. Beginning Science: Biology. Oxford University Press, 1983.
[35] Rolf G Beutel, Frank Friedrich, Xing-Ke Yang, and Si-Qin Ge. Insect Mor-
phology and Phylogeny: A Textbook for Students of Entomology. Walter de
Gruyter, 2013.
[36] Reginald Frederick Chapman. The insects: structure and function. Cambridge
university press, 1998.
[37] Zill Lab Sasha N. Zill. Mechanoreceptors in stick insects.
http://www.marshall.edu/bms/2012/08/30/sasha-n-zill-ph-d/. Last accessed:
June 22 2015.
[38] Robin J Wootton. Functional morphology of insect wings. Annual review of
entomology, 37(1):113–140, 1992.
[39] MH Dickinson, S Hannaford, and J Palka. The evolution of insect wings and
their sensory apparatus. Brain, behavior and evolution, 50(1):13–24, 1997.
230
[40] Michael H Dickinson. Haltere–mediated equilibrium reflexes of the fruit fly,
drosophila melanogaster. Philosophical Transactions of the Royal Society of
London B: Biological Sciences, 354(1385):903–916, 1999.
[41] Robert Dudley. The biomechanics of insect flight: form, function, evolution.
Princeton University Press, 2002.
[42] PJ Albert, RY Zacharuk, and L Wong. Structure, innervation, and distribu-
tion of sensilla on the wings of a grasshopper. Canadian Journal of Zoology,
54(9):1542–1553, 1976.
[43] Michael H Dickinson, Fritz-Olaf Lehmann, and KG Gotz. The active control of
wing rotation by drosophila. The Journal of experimental biology, 182(1):173–
189, 1993.
[44] SA Combes and TL Daniel. Flexural stiffness in insect wings i. scaling and
the influence of wing venation. Journal of experimental biology, 206(17):2979–
2987, 2003.
[45] Pooya Ghaderi, Andrew J Dick, Jason R Foley, and Gregory Falbo. Spectral
domain force identification of impulsive loading in beam structures. In Topics
in Nonlinear Dynamics, Volume 3, pages 157–165. Springer, 2012.
[46] James F Doyle. Force identification from dynamic responses of a bimaterial
beam. Experimental Mechanics, 33(1):64–69, 1993.
[47] Peter F Lester. Turbulence: A new perspective for pilots. Jeppesen Sanderson,
1994.
[48] Francois Vandenberghe Robert Sherman. Turbulence forecast for manned and
unmanned aerial vehicle, 2010.
[49] Aviation Explorer Magazine. Atmospheric turbulence.
http://www.aviationexplorer.com/. Last accessed: June 22 2015.
[50] Michael V Cook. Flight dynamics principles: a linear systems approach to
aircraft stability and control. Butterworth-Heinemann, 2012.
[51] Stacey Gage. Creating a unified graphical wind turbulence model from multi-
ple specifications. In AIAA Modeling and simulation technologies conference
and exhibit, page 5529, 2003.
[52] Frederic M Hoblit. Gust loads on aircraft: concepts and applications. Aiaa,
1988.
[53] JT Vance, I Faruque, and JS Humbert. Kinematic strategies for mitigat-
ing gust perturbations in insects. Bioinspiration & biomimetics, 8(1):016004,
2013.
231
[54] Hugh J Schmittle. Ultra-light aircraft with freely rotating rigid wing, Febru-
ary 4 1986. US Patent 4,568,043.
[55] Vikramjit Singh, Logan Warren, Nathan Putnam, B Walther, P Becker,
A Danielson, Babar Koraishy, W Wood, Dan Jensen, and Andy Szmerekovsky.
A novel exploration into gust resistant operation of mavs/uavs through trans-
formation. In Second US-Euro MAV Conference, Destin, FL, 2006.
[56] Adetunji Oduyela and Nathan Slegers. Gust mitigation of micro air vehicles
using passive articulated wings. The Scientific World Journal, 2014, 2014.
[57] Matt Shields and Kamran Mohseni. Passive mitigation of roll stall for low
aspect ratio wings. Advanced Robotics, 27(9):667–681, 2013.
[58] Ted L. Lomax. Structural loads analysis for commercial transport aircraft:
theory and practice. Aiaa, 1996.
[59] WF Grosser, WW Hollenbeck, and DC Eckholdt. The c-5 a active lift distribu-
tion control system. AGARD Impact of Active Control Technol. on Airplane
Design 18 p(SEE N 75-30027 21-01), 1971.
[60] Mordechay Karpel. Design for active flutter suppression and gust alleviation
using state-space aeroelastic modeling. Journal of Aircraft, 19(3):221–227,
1982.
[61] Richard W Wlezien, Garnett C Horner, Anna-Maria R McGowan, Sharon L
Padula, Michael A Scott, Richard J Silcox, and Joycelyn S Harrison. Air-
craft morphing program. In 5th Annual International Symposium on Smart
Structures and Materials, pages 176–187. International Society for Optics and
Photonics, 1998.
[62] Jie Zeng, Boris Moulin, Raymond De Callafon, and Martin Brenner. Adaptive
feedforward control for gust load alleviation. Journal of Guidance, control, and
dynamics, 33(3):862–872, 2010.
[63] Santiago Basualdo. Load alleviation on wind turbine blades using variable
airfoil geometry. Wind Engineering, 29(2):169–182, 2005.
[64] D McLean. Gust-alleviation control systems for aircraft. In Proceedings of the
Institution of Electrical Engineers, volume 125, pages 675–685. IET, 1978.
[65] Volk J. A. Dreim D. R. & Applewhite K. A. Britt, R. T. aeroservoelastic
characteristics of the b-2 bomber and implications for future large aircraft,
structural aspects of flexible aircraft. 1999.
[66] H Bo¨hret, B Krag, and J Skudridakis. Olga: An open loop gust alleviation
system. AGARD Active Control Systems: Rev, Evaluation and Projections 16
p(SEE N 85-27883 17-08), 1985.
232
[67] RM Rennie and EJ Jumper. Gust alleviation using trailing-edge flaps. In 37th
AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 1999.
[68] David C Soreide, Rodney K Bogue, LJ Ehernberger, and Harold R Bagley.
Coherent lidar turbulence measurement for gust load alleviation. In SPIE’s
1996 International Symposium on Optical Science, Engineering, and Instru-
mentation, pages 61–75. International Society for Optics and Photonics, 1996.
[69] Raymond L Bisplinghoff and Holt Ashley. Half man, robert l.: Aeroelasticity,
1955.
[70] Raymond L Bisplinghoff, Holt Ashley, and Robert L Halfman. Aeroelasticity.
Courier Corporation, 2013.
[71] Donald McLean. Automatic flight control systems(book). Englewood Cliffs,
NJ, Prentice Hall, 1990, 606, 1990.
[72] Mark Drela. Aswing 5.81 technical description, steady formulation. Technical
report, Technical report, MIT, 2008.
[73] TH Skopinski, William S Aiken, and Wilber B Huston. Calibration of strain-
gage installations in aircraft structures for the measurement of flight loads. US
Government Printing Office, 1954.
[74] Boyd Perry III. Qualitative comparison of calculated turbulence responses
with with an active control system. 1981.
[75] M Frederick, EC Kerrigan, and JMR Graham. Gust alleviation using rapidly
deployed trailing-edge flaps. Journal of Wind Engineering and Industrial Aero-
dynamics, 98(12):712–723, 2010.
[76] Peter Bjoern Andersen, Lars Henriksen, Mac Gaunaa, Christian Bak, and
Thomas Buhl. Deformable trailing edge flaps for modern megawatt wind
turbine controllers using strain gauge sensors. Wind Energy, 13(2-3):193–206,
2010.
[77] Eric Lee Brown. Integrated strain actuation in aircraft with highly flexible
composite wings. PhD thesis, Massachusetts Institute of Technology, 2003.
[78] Kenneth B Lazarus, Edward F Crawley, and Charrissa Y Lin. Multivariable
active lifting surface control using strain actuation: analytical and experimen-
tal results. Journal of Aircraft, 34(3):313–321, 1997.
[79] Charrissa Y Lin, Edward F Crawley, and Jennifer Heeg. Open-and closed-
loop results of a strain-actuated active aeroelastic wing. Journal of Aircraft,
33(5):987–994, 1996.
[80] Nhan Nguyen and Ilhan Tuzcu. Flight dynamics of flexible aircraft with aeroe-
lastic and inertial force interactions. In AIAA Atmospheric Flight Mechanics
Conference, AIAA Paper, volume 6045, 2009.
233
[81] A De Gaspari, Sergio Ricci, Luca Riccobene, and Alessandro Scotti. Active
aeroelastic control over a multisurface wing: Modeling and wind-tunnel test-
ing. AIAA journal, 47(9):1995–2010, 2009.
[82] Bing Feng Ng, Rafael Palacios, J Michael R Graham, and Eric C Kerrigan.
Robust control synthesis for gust load alleviation from large aeroelastic models
with relaxation of spatial discretisation. Proceedings of the EWEA, 2012.
[83] Peter M Suh, Alexander W Chin, and Dimitri N Mavris. Virtual deformation
control of the x-56a model with simulated fiber optic sensors. In Proceedings
AIAA Atmospheric, Flight Mechanics Conference, 2014.
[84] Animesh Chakravarthy, Katie A Evans, and Johnny Evers. Sensitivities &
functional gains for a flexible aircraft-inspired model. In American Control
Conference (ACC), 2010, pages 4893–4898. IEEE, 2010.
[85] Stephen Chiu, Sujeet Chand, Doug Moore, and Ashwani Chaudhary. Fuzzy
logic for control of roll and moment for a flexible wing aircraft. Control Sys-
tems, IEEE, 11(4):42–48, 1991.
[86] M Cicala, VR Baraniello, and A Sollazzo. A load alleviation application of an
elastic motion estimator for a flexible uav. In Control & Automation (MED),
2011 19th Mediterranean Conference on, pages 1229–1234. IEEE, 2011.
[87] Y Matsuzaki, T Ueda, Y Miyazawa, and H Matsushita. Gust load alleviation
of a transport-type wing-test and analysis. Journal of aircraft, 26(4):322–327,
1989.
[88] Mayuresh J Patil and Dewey H Hodges. Flight dynamics of highly flexible
flying wings. Journal of Aircraft, 43(6):1790–1799, 2006.
[89] Gregory Gremillion. Bio-inspired disturbance rejection with ocellar and dis-
tributed acceleration sensing for small unmanned aircraft systems. 2015.
[90] Rizwan Butt. Introduction to Numerical Analysis Using MATLAB R©. Jones
& Bartlett Learning, 2009.
[91] Padula Sharon L and Kincaid Rex K. Optimization strategies for sensor and
actuator placement. Technical report, NASA Langley Technical Report Server,
1999.
[92] Firdaus E Udwadia. Methodology for optimum sensor locations for param-
eter identification in dynamic systems. Journal of Engineering Mechanics,
120(2):368–390, 1994.
[93] Giuseppe Quaranta, Giuseppe Carlo Marano, Francesco Trentadue, and Gior-
gio Monti. Numerical study on the optimal sensor placement for historic
swing bridge dynamic monitoring. Structure and Infrastructure Engineering,
10(1):57–68, 2014.
234
[94] Giuseppe C Marano, Giuseppe Quaranta, Rita Greco, and Giorgio Monti. Sen-
sor network design for monitoring a historic swing bridge. In Proceedings of
the International Symposium on Engineering under Uncertainty: Safety As-
sessment and Management (ISEUSAM-2012), pages 493–502. Springer, 2013.
[95] Daniel C Kammer. Sensor placement for on-orbit modal identification and
correlation of large space structures. Journal of Guidance, Control, and Dy-
namics, 14(2):251–259, 1991.
[96] Masoud Sanayei and Chitra N Javdekar. Sensor placement for parameter
estimation of structures using fisher information matrix. Proc., Application of
Advanced Technology in Transportation, pages 386–393, 2002.
[97] Wodek K Gawronski. Actuator and sensor placement. In Dynamics and
Control of Structures, pages 100–128. Springer, 1998.
[98] Daniel C Kammer. Optimal sensor placement for modal identification using
system-realization methods. Journal of Guidance, Control, and Dynamics,
19(3):729–731, 1996.
[99] COMSOL Multiphysics. Structural mechanics module–users guide, 2013.
[100] William L Oberkampf, MM Sindir, and AT Conlisk. Guide for the verification
and validation of computational fluid dynamics simulations. Am. Institute of
Aeronautics and Astronautics, 1998.
[101] NASA NPARC. Computational fluid dynamics (cfd) verification and valida-
tion.
http://www.grc.nasa.gov/WWW/wind/valid/. Last accessed: June 22 2015.
[102] Max Chtangeev. Borealis composites.
max.chtangeev@borealiscomposites.com. Last accessed: June 22 2015.
[103] Patrick R Amestoy, Ian S Duff, and J-Y L’Excellent. Mumps multifrontal
massively parallel solver version 2.0. 1998.
[104] Mark F Costello. Theory for the analysis of rotorcraft operating in atmospheric
turbulence. In Annual Forum Proceedings- American Helicopter Society, vol-
ume 2, pages 1003–1015, 1990.
[105] Professor Inderjit Chopra. Alfred gessow rotorcraft center, university of mary-
land, college park.
http://www.agrc.umd.edu//. Last accessed: June 22 2015.
[106] Professor Allen Winkelmann. Aerospace engineering laboratories -aerolab,
university of maryland, college park.
http://www.aero.umd.edu/research-labs. Last accessed: June 22 2015.
235
[107] University of Florida. Reef aerodynamic characterization facility wind tunnel.
www.reef.ufl.edu. Last accessed: June 22 2015.
[108] Michael Sytsma and Lawrence Ukeiley. Wind tunnel generated turbulence.
AIAA Paper, 1159, 2011.
[109] GA Williamson, BD McGranahan, BA Broughton, RW Deters, JB Brandt,
and MS Selig. Summary of low-speed airfoil data. University of Illinois Low
Speed Airfoil Tests, 2012.
236