ABSTRACT Title of dissertation: CONTROL-ORIENTED REDUCED ORDER MODELING OF DIPTERAN FLAPPING FLIGHT Imraan Faruque, Doctor of Philosophy, 2011 Dissertation directed by: Professor J Sean Humbert Department of Aerospace Engineering Flying insects achieve ight stabilization and control in a manner that re- quires only small, specialized neural structures to perform the essential components of sensing and feedback, achieving unparalleled levels of robust aerobatic ight on limited computational resources. An engineering mechanism to replicate these con- trol strategies could provide a dramatic increase in the mobility of small scale aerial robotics, but a formal investigation has not yet yielded tools that both quantita- tively and intuitively explain apping wing ight as an \input-output" relationship. This work uses experimental and simulated measurements of insect ight to cre- ate reduced order ight dynamics models. The framework presented here creates models that are relevant for the study of control properties. The work begins with automated measurement of insect wing motions in free ight, which are then used to calculate ight forces via an empirically-derived aerodynamics model. When paired with rigid body dynamics and experimentally measured state feedback, both the bare airframe and closed loop systems may be analyzed using frequency domain system identi cation. Flight dynamics models describing maneuvering about hover and cruise conditions are presented for example fruit ies (Drosophila melanogaster) and blow ies (Calliphorids). The results show that biologically measured feedback paths are appropriate for ight stabilization and sexual dimorphism is only a mi- nor factor in ight dynamics. A method of ranking kinematic control inputs to maximize maneuverability is also presented, showing that the volume of reachable con gurations in state space can be dramatically increased due to appropriate choice of kinematic inputs. CONTROL-ORIENTED REDUCED ORDER MODELING OF DIPTERAN FLAPPING FLIGHT by Imraan Faruque Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial ful llment of the requirements for the degree of Doctor of Philosophy 2011 Advisory Committee: Professor J Sean Humbert, Chair Professor Sarah Bergbreiter Professor Roberto Celi Professor Inderjit Chopra Professor Gerald Wilkinson ? Copyright by Imraan A Faruque 2011 Acknowledgments I owe my gratitude to all the people who made these discoveries possible. I thank my advisor Dr Sean Humbert, for the engaging problems and resources he has provided me, and for his apparently limitless enthusiasm. I?m surprised I ended up studying bugs for this long, and I am sure part of that is due to his curiosity, excitement, and ideas. I am quite grateful for my committee members? contributions as well. Dr Celi?s expertise with system identi cation and all things rotorcraft was invaluable, even if he refuses to y on them. Dr Sanner?s wealth of control theoretic knowledge allowed him issue an un-matched academic challenge and present incredibly focused questions. I?m thankful for Dr Chopra, who became an aerodynamics reference, and Dr Bergbreiter, who introduced me to a microscopic world of robots, as well as Dr Wilkinson, who introduced me to antlike bugs with eyes on the ends of stalks like aliens. Many in the AVL contributed to this research, most particularly through Dr Jason Vance?s biological expertise, Nick Kostreski?s experimental work with insects, and Mac MacFarlane?s computational work. I thank the rest of the AVL for their patience with my late work hours, and persistently inquisitive nature, both academic and philosophical. Several agencies contributed nancially, particularly the Army, O ce of Naval Research, and Air Force Research Lab. I?m especially thankful to the Army Research Lab for hosting me during a portion of my graduate tenure. ii I am thankful for the insects who proved that controlling ight at this scale was possible, and to the God who made them. Most of all, I?m grateful for my family, who made this possible. For my father, who challenged me to think and who inspired all of my academic pursuits, especially this one. I?m perpetually grateful for your sacri ce and, as ridiculously long as this document is, it?s not long enough to express your e ect on me. For my mother, who humored and even encouraged my investigative hobby that eventually grew into a full time job. For Ruel, who always believed I was smart enough to do anything. For Anya, who always believed I could work hard enough to achieve anything. And for Orion, who always believed I was brave enough to throw myself into anything. iii Table of Contents List of Tables viii List of Figures ix List of Nomenclature and Abbreviations xiii 1 Introduction 1 1.1 Motivation and Technical Challenges . . . . . . . . . . . . . . . . . . 1 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Previous Modeling Work . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5 Dissertation Contributions . . . . . . . . . . . . . . . . . . . . . . . . 9 1.6 Dissertation Organization . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Wing Kinematics Measurement 13 2.1 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Insect Subjects and Experimental Apparatus . . . . . . . . . . . . . . 15 2.3 Kinematics Digitization . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Coordinate De nitions . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Digitization Algorithm . . . . . . . . . . . . . . . . . . . . . . 17 2.3.3 Digitization Results . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.4 Digitization Accuracy . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Kinematics Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.1 Time Synchronous Averaging . . . . . . . . . . . . . . . . . . 24 2.4.2 Curve Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3 Longitudinal Hover Dynamics 27 3.1 Introduction and Background . . . . . . . . . . . . . . . . . . . . . . 27 3.1.1 Wing Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1.2 Control Parameters . . . . . . . . . . . . . . . . . . . . . . . . 29 3.1.3 Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1.3.1 Translational Lift and Drag . . . . . . . . . . . . . . 31 3.1.4 Wingstroke-Averaged Forces and Moments . . . . . . . . . . . 33 3.2 Derivation of Linearized Flight Dynamics about Hover . . . . . . . . 34 3.2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.2 Hover Equilibrium Trim Solutions . . . . . . . . . . . . . . . . 35 3.2.3 Inclusion of Perturbation States in the Quasi-Steady Form . . 38 3.2.4 Velocity Components . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.4.1 Local Velocity Component due to Flapping . . . . . 39 3.2.4.2 Local Velocity Component due to Body Velocity . . 40 3.2.4.3 Local Velocity Component due to Rotation Rates . . 40 3.2.5 E ects of Modi ed Velocity . . . . . . . . . . . . . . . . . . . 41 iv 3.2.5.1 Angle of Attack Modi cation . . . . . . . . . . . . . 41 3.2.5.2 E ect on Dynamic Pressure . . . . . . . . . . . . . . 42 3.2.6 Perturbation Equations . . . . . . . . . . . . . . . . . . . . . . 43 3.2.7 Homogeneous System Dynamics (the A matrix) . . . . . . . . 44 3.2.8 Control Derivatives (the B matrix) . . . . . . . . . . . . . . . 44 3.2.8.1 Heave and Pitch Dynamics . . . . . . . . . . . . . . 45 3.2.8.2 Signi cance of the Control Matrix B . . . . . . . . . 46 3.3 Simulation and System Identi cation . . . . . . . . . . . . . . . . . . 48 3.3.1 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.3.1.1 Haltere Modeling . . . . . . . . . . . . . . . . . . . . 49 3.3.2 System identi cation method . . . . . . . . . . . . . . . . . . 50 3.3.3 Model structure . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4.1 Non-Parametric Handling Qualities Identi cation . . . . . . . 53 3.4.2 Heave Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.2.1 E ect of Flap Amplitude . . . . . . . . . . . . . . . . 58 3.4.2.2 E ect of Flap Frequency . . . . . . . . . . . . . . . . 58 3.4.3 Longitudinal Dynamics . . . . . . . . . . . . . . . . . . . . . . 61 3.4.3.1 Pitch Dynamics . . . . . . . . . . . . . . . . . . . . . 64 3.4.3.2 Forward and Heave Velocity Dynamics . . . . . . . . 66 3.4.4 Identi ed Longitudinal Dynamics Veri cation . . . . . . . . . 66 3.4.5 Accuracy of Identi ed Model . . . . . . . . . . . . . . . . . . . 67 3.4.6 Bare Airframe System Identi cation . . . . . . . . . . . . . . . 69 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4 Lateral-Directional Hover Dynamics 73 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2.1 Wing Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2.2 Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.3 System Identi cation . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4 Wingstroke Averaged Forces and Moments . . . . . . . . . . . . . . . 77 4.4.1 Control Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4.1.1 Roll Moment due to Di erential Stroke Amplitude d 78 4.4.1.2 Roll Moment due to Di erential Stroke Bias o ,d . . 79 4.4.1.3 Roll Moment due to Di erential Stroke Plane Incli- nation d . . . . . . . . . . . . . . . . . . . . . . . . 79 4.4.2 Physical Basis for Passive Aerodynamic Mechanisms . . . . . 80 4.4.2.1 Sideslip Damping . . . . . . . . . . . . . . . . . . . . 80 4.4.3 Active Feedback Mechanisms . . . . . . . . . . . . . . . . . . 83 4.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.5.1 Yaw Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.5.2 Roll Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.5.2.1 Linear System . . . . . . . . . . . . . . . . . . . . . 86 4.5.2.2 Physical Mechanism Predicting Roll Damping . . . . 86 v 4.5.3 Full Lateral Dynamics . . . . . . . . . . . . . . . . . . . . . . 87 4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5 Model Validation 93 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.2 Unmodeled Aerodynamics and Velocity Perturbations . . . . . . . . . 95 5.2.0.1 Rotational Lift and Other Aerodynamic Models . . . 95 5.2.0.2 Reverse Flow . . . . . . . . . . . . . . . . . . . . . . 95 5.2.0.3 Spanwise Flow . . . . . . . . . . . . . . . . . . . . . 96 5.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.3.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.3.2 Reduced-Order Aerodynamic Model . . . . . . . . . . . . . . . 97 5.3.3 Computational Methodology . . . . . . . . . . . . . . . . . . . 98 5.4 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6 Feedback E ects and Parameter Variation 106 6.1 Mechanosensory Feedback Variation . . . . . . . . . . . . . . . . . . . 106 6.1.1 Feedback Variation Study . . . . . . . . . . . . . . . . . . . . 106 6.1.2 Gain Variation Results . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Mechanosensory Feedback on CFD model . . . . . . . . . . . . . . . . 110 6.3 Parameter (Phenotypic) Variation . . . . . . . . . . . . . . . . . . . . 112 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7 Reachability Analysis 115 7.1 Longitudinal Flight Dynamics Model . . . . . . . . . . . . . . . . . . 117 7.2 Reachability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.2.0.1 Controllability operator . . . . . . . . . . . . . . . . 120 7.2.1 Reachable space under unit norm input . . . . . . . . . . . . . 122 7.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 8 Forward Flight Dynamics 128 8.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 8.1.1 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . 130 8.2 Longitudinal Flight Dynamics Modeling . . . . . . . . . . . . . . . . 130 8.2.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 8.2.2 Aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 8.2.3 Rigid Body Dynamics and System Identi cation . . . . . . . . 132 8.3 Lateral-Directional Forward Flight Dynamics . . . . . . . . . . . . . . 133 8.3.1 Kinematics and System Identi cation . . . . . . . . . . . . . . 137 8.3.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 138 8.3.2.1 Roll dynamics . . . . . . . . . . . . . . . . . . . . . . 138 8.3.2.2 Yaw dynamics . . . . . . . . . . . . . . . . . . . . . 140 8.3.2.3 Roll/Yaw Coupling (Proverse Yaw) . . . . . . . . . . 141 8.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 vi 9 Conclusions and Future Work 144 9.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 9.2 Dissertation Contributions . . . . . . . . . . . . . . . . . . . . . . . . 145 9.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 vii List of Tables 3.1 Stroke plane inclination in a turning Calliphorid shows di erential stroke plane inclination perturbation and constant collective stroke plane inclination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Range of parameter variation observed in Calliphorid ies throughout maneuvering ight. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3 Control and handling qualities properties used to evaluate the system in reference to rotorcraft handling qualities requirements. GM and PM represent gain and phase margins, resctively. . . . . . . . . . . . 55 3.4 Nominal input parameters used in evaluating B for the hover lin- earization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.5 Uncertainty estimates for identi ed dynamics model. . . . . . . . . . 68 4.1 Nominal input parameters used in evaluating B for the hover lin- earization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.2 Uncertainty estimates for parameters in the identi ed dynamics model. 90 5.1 Stability derivatives as calculated by the curve- tted model and IBINS.103 6.1 Control and stability derivative results for female Drosophila. . . . . . 113 7.1 Control input performance ranking criteria for u1 though u4. . . . . . 126 8.1 Uncertainty values in the parameters are acceptable, with the excep- tion of the w derivatives that poorly identi ed. . . . . . . . . . . . . . 136 8.2 Uncertainty estimates for parameters in the identi ed dynamics model show that roll damping and roll authority is only weakly a ected by increasing u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 8.3 Yaw control authority increases slightly with forward speed, while yaw damping is relatively una ected. . . . . . . . . . . . . . . . . . . 140 viii List of Figures 1.1 Directions in which insect compound eyes and ocelli (rudimentary eyes) are most sensitive (Parsons et al., 2010) may correspond to dynamic properties of the airframe. . . . . . . . . . . . . . . . . . . . 2 2.1 Experimental test rig and free ight chamber. . . . . . . . . . . . . . 16 2.2 Axes and angle de nitions. . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Example images versus voxel reconstructions. . . . . . . . . . . . . . 19 2.4 Visual hull reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Digitized ight sequences include a (a) straight and level sequence, (b) level left turn, (c) climb initiation (d) yaw motion (sideslip). All images are top views shown at 70Hz except for (c), which is a side view shown at 140Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Stroke angle is well detected by principal component analysis, and multiple wingstrokes in straight and level ight provide ample mea- surement of the nominal kinematic pattern. . . . . . . . . . . . . . . 21 2.7 A spherical exclusion rule dramatically reduces the number of body voxels incorrectly assigned to the wings. Blue indicates body voxels that were previously detected as wing voxels. . . . . . . . . . . . . . . 22 2.8 Modi cations to the wing angle detection algorithm improved the results dramatically. . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.9 A reconstruction of the 45 insect model reconstruction shows an occlusal defect on the left wing. . . . . . . . . . . . . . . . . . . . . . 23 2.10 Time synchronous averaging is used to identify a nominal stroke pat- tern associated with a reference condition (straight and level ight shown). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.11 Wing pitch angle in forward ight is more complex than a smoothed square wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1 Axes and angle de nitions. . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Insect apping kinematics parameters are used as control inputs to model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3 De nitions of stability frame velocities vG = us^x+vs^y+ws^z, rotation rates !B = ps^x+qs^y+rs^z, forces F = Xs^x+Y s^y+Zs^z, and moments M = Ls^x +Ms^y +Ns^z. . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Trim inputs 0 and o ,0 as functions of body hovering angle . . . . . 38 3.5 Wing angle of attack is increased in a descent perturbation and re- duced in a climb perturbation, so a apping wing system in axial climb/descent exhibits damping in heave velocity w. . . . . . . . . . . 42 3.6 Simulator used for haltere-on and bare-airframe input/output studies, showing (highlighted path) the critical modi cation: consideration of rigid body state in wing aerodynamics functions. . . . . . . . . . . . . 49 ix 3.7 This study considered the haltere-stabilized insect mapping inputs at point A to outputs at point B. A low gain controller was included that modi ed the frequency sweep inputs to keep the system near hover. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.8 Identi ed transfer functions to pitch rate q. . . . . . . . . . . . . . . . 54 3.9 Identi ed transfer functions to pitch angle . . . . . . . . . . . . . . . 55 3.10 Identi ed transfer functions to forward velocity u. . . . . . . . . . . . 56 3.11 Identi ed transfer function from stroke amplitude to heave velocity w. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.12 As noted in experimental observations of insects, stroke plane ampli- tude shows a linear correlation to perturbations. . . . . . . . . . . 60 3.13 The use of ap frequency f as an input is correlated linearly over a smaller region than . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.14 The identi ed system ts the stroke plane amplitude to heave velocity relation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.15 Time domain comparison of the linear and nonlinear systems. . . . . 67 3.16 Map of haltere-on system poles with uncertainty perturbations as listed in Table (3.5) preserves the qualitative behavior. . . . . . . . . 69 3.17 Map of haltere-o system (the bare-airframe) poles including uncer- tainty perturbations indicates the traditional VTOL modal structure (Sun and Xiong, 2005). . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1 Control input de nitions. . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 An imposed sideslip increases and decreases relative air ow over re- gions of the wingstroke. . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3 An imposed sideslip increases and decreases drag forces during the advancing and retreating strokes in a manner leading to damping in sideslip v about hover. . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.4 Comparison of calculated Yv with Boeing 747, UH-1 Huey, BO-105 helicopter, and XC-142 tactical transport. All values normalized by Froude number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.5 Yaw input d is well correlated linearly to yaw rate r and the 3 avail- able models show excellent spectral agreement. . . . . . . . . . . . . . 85 4.6 Roll input d drives the roll/sideslip dynamics, but not yaw. . . . . . 87 4.7 The applied roll motion in (a) creates a counter-moment in (c) via the di erence in angle of attacks imposed by the velocity distribution in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.8 Comparison of the poles identi ed for the coupled and uncoupled system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.9 The identi ed system pole structure is preserved under uncertainty perturbations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.1 The model drosophila in the computational domain used for IBINS calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 x 5.2 Stroke-averaged forces (a) and moments (b) as given by the curve- tted experimental model (red) and IBINS (blue) for various pertur- bations in u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.3 Stroke-averaged forces (a) and moments (b) as given by the curve- tted experimental model (red) and IBINS (blue) for various pertur- bations in w. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.4 Stroke-averaged forces (a) and moments (b) as given by the curve- tted experimental model (red) and IBINS (blue) for various pertur- bations in q. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.5 Isosurfaces of Q-criterion of the model Drosophila in hover (a), surge u = 0:38 m=s (b), heave w = 0:38 m=s (c), and pitch q = 2:1 rad=s (d) just after pronation. The disturbed ow due to the previous wingstroke remains only a single chord length away from the returning stroke in all cases, and in uences the loads on the wing. 104 5.6 Longitudinal poles of the model Drosophila in hover as calculated us- ing the curve- tted experimental aerodynamics model (red) and CFD (blue). Poles found via system identi cation with the experimentally- derived model in Faruque and Humbert (2010) are plotted alongside for reference (green). . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.1 Kinematic inputs de nitions. . . . . . . . . . . . . . . . . . . . . . . . 108 6.2 The stability properties of the insect dynamics are improved by hal- tere feedback, and a root locus plot over haltere gain shows pole movement to connect the systems. . . . . . . . . . . . . . . . . . . . . 109 6.3 Stabilization of the CFD-derived system (in red) requires a higher pitch rate gain than than originally modeled (in blue). Equivalent convergence rate requires a gain 2.45 times larger (in black). . . . . . 111 6.4 Comparison of the lateral-directional poles identi ed for the nominal (male) and larger (female) insect. . . . . . . . . . . . . . . . . . . . . 113 7.1 Longitudinal control inputs used in reachability analysis. (A) Stroke plane inclination c, (B) Stroke plane o set o , and (C) Di erential wing angle ud. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 7.2 Controllability ellipsoids for the input combinations illustrate the reachable con gurations under the restriction jjuijj 1. Input com- binations u1 through u3 are pairs of control terms, while u4 considers all 3 control terms ( c; o ; ud). . . . . . . . . . . . . . . . . . . . . . 125 7.3 Controllability of input combinations u1 through u4, as ranked by the determinant or Frobenius norm of the square root of the controllabil- ity gramian. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 8.1 Kinematics in forward and aft strokes. . . . . . . . . . . . . . . . . . 131 8.2 The pole structure of is similar to those in hover, with a real heave mode, and a coupled rotationally dominated triple that includes an unstable oscillatory pair. . . . . . . . . . . . . . . . . . . . . . . . . . 133 xi 8.3 A 4 state linear system o transfer function shows good agreement for all but the heave direction, which shows model structure error. . . 134 8.4 A 4 state longitudinal linear system transfer function shows good agreement for all but the heave direction, which shows model struc- ture error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 8.5 In forward ight u = 2 m/s, rst order linear roll damping is an acceptable spectral description of the roll transfer function. . . . . . . 139 8.6 Time-domain agreement of the response to d input doublets in for- ward ight is excellent. . . . . . . . . . . . . . . . . . . . . . . . . . . 139 8.7 Yaw dynamics transfer function t at u = 0m/s. . . . . . . . . . . . . 141 8.8 Roll rate to yaw rate coupling at 2m/s for di erening wing pitch input ud. Roll motions in forward ight induce smaller yet potentially signi cant yaw rates. The magnitude of the induced yaw rate is a function of ud. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 xii Nomenclature Wing angle of attack rad g Geometric angle of attack rad add Angle of attack perturbation rad ud Wing pitch upstroke/downstroke input rad Stroke plane angle rad Coherence nondim. Wing pitch angle rad Body pitch angle rad i Eigenvalue i Advance ratio nondim. Body hovering angle rad Fluid (air) density kg=m3 Stroke angle rad Fore/aft wingstroke angle bias rad Stroke amplitude rad Controllability operator !B Body angular velocity rad=s A System dynamics matrix A Wing aspect ratio nondim. b0; b1; b2 Planar wing t coe cients B Control input matrix c Geometric mean chord m CL Wing lift coe cient nondim. CD Wing drag coe cient nondim. L(t) Instantaneous lift force N dx; dz Components of rG=t in S coordinates m D(t) Instantaneous drag force N f Flap frequency Hz g Gravitational acceleration 9.8 m=s2 Gxx(!) Autospectral density Gxy(!) Cross-spectral density I[ ] Vehicle moments of inertia kg m2 k;K Gain, Gain matrix m Vehicle mass kg p; q; r Body angular rates !B in S coordinates rad=s p^x; p^y; p^z Stroke plane axes P r Spanwise distance outboard along a wing m r^2 Second moment of wing area squared nondim. ra Spanwise moment arm of lift force m rG=p Vector from center of mass G to point P m rG=t Vector from center of mass G to thorax center m fr^x; r^y; r^zg Right wing axes R R Wing length m xiii R2 Coe cient of determination RRS Rotation matrix from S to R fs^x; s^y; s^zg Stability frame axes S S Wing area m2 t Time sec T Flap period T = 1=f sec u(t) Input vector u^(t) Optimal input vector U Input space u; v; w Body translational rates vG in S coordinates m=s ur; vr; wr Body translational rates vG in R coordinates m=s ut Wing tip velocity m=s v Total ow velocity incident on wing m=s vx; vy; vx Total ow velocity incident on wing v in R coordinates m=s v Flow velocity on wing due to apping m=s v! Flow velocity on wing due to body rotation !B rad=s vG Body velocity m=s x(t) State vector xi Points in state space, i = 0; 1::: X; Y; Z Stability frame forces in S coordinates N Xc Controllability Gramian X Con guration space Abbreviations CFD Computational Fluid Dynamics GM Gain Margin IBINS Immersed Boundary Incompressible Navier Stokes MAV Micro Air Vehicle LTI Linear Time Invariant PM Phase Margin xiv Chapter 1 Introduction 1.1 Motivation and Technical Challenges Flying insects represent simultaneously the highest standard of performance in aerial mobility and the animals that achieve robust ight control with the smallest neurological structures. The prospect of agile control performance that requires small, specialized, computational structures is very attractive from the perspective of micro air vehicle (MAV) design. An understanding of the dynamic properties that insect ight leverages to achieve aerial mobility on a limited computational budget would provide insight into the apping wing MAV control problem. Measurements of insect sensory systems has revealed that insects stabilize ight by sensing composite states measured along highly non orthogonal axes, then fusing those inputs into a control law. The directionality of such measurements, as seen in Fig 1.1, is in opposition to traditional engineering measurement of vehicle states. A popular hypothesis has been that the particular combinations insects mea- sure are more important for ight stabilization than traditional isolated, orthogonal measurements like roll or pitch rate (Taylor and Krapp, 2007). If this is the case, then dynamic models should show that insect sensing is most responsive to mo- tions (or directions) that have the highest feedback requirements. Dynamic models are necessary to determine the feedback requirements along di erent axes, which 1 Figure 1.1: Directions in which insect compound eyes and ocelli (rudimentary eyes) are most sensitive (Parsons et al., 2010) may correspond to dynamic properties of the airframe. would allow researchers to determine if the sensor outputs correspond to charac- teristic ight modes that require feedback. While full-scale aircraft actuators often excite along orthogonal axes, insects do not necessarily do so and actuation along non-orthogonal axes is equally possible. The insect sensory systems may be tuned not only to the motions requiring stabilization feedback, but also to the e ects of its non-orthogonal actuation capabilities, such that the insect?s sensors provide it with feedback on the e ects of its actuation. An insect with access to measurements describing the directions resulting from its actuation would be exceptionally able to regulate this actuation. Again, testing this hypothesis requires a dynamic model that describes the motion resulting from each actuator, in order to evaluate whether the sensory systems are tuned to the results of non-orthogonal actuation. A concise dynamic model would allow researchers to evaluate whether insect sensory systems are speci cally tuned to dynamic properties such as insect ight modes or actuator results. Despite the need for a concise model, previous investi- gations have shown that the ight regime is complicated by complex aerodynamic e ects, nonlinear wing kinematics and body dynamics, and small perturbations in- 2 tegrated over many wingstrokes, as well as adaptive sensing and feedback. Detailed modeling and simulation of the aerodynamic e ects alone is an involved process that typically requires several weeks for each simulated case. The need for tractable insect ight dynamics modeling leads directly to the problem statement of this dis- sertation: Problem Statement: Find reduced-order dynamics models for small- scale dipteran apping ight, and determine a path for creating mod- els useful for estimating sensing and feedback requirements. The derived models should be appropriate for control analysis, in partic- ular, answer questions of stability, rates of convergence or instability, modal participation, and feedback requirements, as well as quantify actuation results (eg, maneuverability). 1.2 Background In recent years, researchers have made much progress into the task of under- standing the aerodynamic basis for and the control architecture involved in insect ight, in particular the fruit y Drosophila melanogaster. Advances in the eld of apping wing aerodynamics have largely relied on the ability of researchers to make detailed observations of the insects? apping behavior. Early observations of tethered Drosophila by Vogel (1967) began to observe variations in certain \stroke parameters" de ning wing kinematic patterns and variations in wing contour. Sev- eral years later, Weis-Fogh (1972) used Vogel?s observations in the development of 3 the \quasi-steady" form rst introduced by Osborne (1951). Brie y, this approach uses curve ts to experimental data to allow the instantaneous lift and drag forces of an insect to be represented by drawing analogy to a similar wing translating at the same angle of attack and the same (steady) velocity. The quasi-steady ref- erences found in insect aerodynamics literature are in contrast to the aero-elastic concept, which refers to ignoring unsteady e ects by retaining only terms dependent steady quantities and does not employ curve ts to experimental data. Because of this, we consider the aerodynamic modeling presented in insect ight literature an experimental aerodynamics model presented in quasi-steady form. The quasi-steady form has been extensively used as a foundation by researchers beginning with Ellington and Dickinson (Ellington, 1984a) to develop the aerody- namic theory used in contemporary understanding and prediction of insect ight. A major use of the theory is in the prediction of baseline forces in order to elucidate the contributions of additional aerodynamic mechanisms, predominantly unsteady e ects (Dickinson and Gotz, 1999). The concept has been applied to determine aero- dynamic contributions by mechanisms such as \clap and ing" movements (Sped- ding and Maxworthy, 1986) and dynamic stall (Dickinson et al., 1999). Even so, the chief contribution of the quasi-steady form to the eld of insect ight understanding has been as a means to reduce kinematic and force data taken from both tethered and in- ight recordings of wing kinematics, allowing reduction of in- ight data to nondimensional coe cients that may be interpreted from the perspective of more traditional aerodynamic mechanisms (Fry et al., 2003). The computational workload an insect must perform to support apping wing 4 ight is still an active area of research. Experimental evidence has indicated phys- iological components involved in ight stabilization, such as tangential cells and descending neurons that translate optic ow estimates to ight motor commands. Despite the measurements made on the sensing and control components, the ight stabilization demands that these structures must support are unknown, but must be addressed via an analysis of the plant they are responsible for controlling (ie, the insect dynamics). Recently, Taylor and Thomas (2003a,b) used measurements of tethered locusts to estimate the stability and control properties associated with dipteran insect apping ight, while Hesselberg and Lehmann (2007) and Hedrick et al. (2009) addressed the yaw dynamics of apping ight. Aerodynamics inves- tigations using computational uid dynamics (CFD) (Sun and Xiong, 2005) have also addressed the longitudinal (motion in the insect?s plane of symmetry) apping ight. A concise understanding of how the control inputs available to a dipteran insect (in the form of kinematic variables) a ect the long term motion of the insect ( ight control) is expected to provide insight into why insects tend to use particular control inputs more readily than others and indicate how coupled control inputs can lead to particular motions. A large body of knowledge exists to provide control analysis on state-space dynamics models (Ogata, 2002), meaning that a process that simpli es the relatively complex nonlinear dynamics of dipteran apping wing insect ight is desirable. State-space models and the process to generate them are goals of this study. The simpli ed linear time invariant dynamics models may then be used to evaluate controllability and design control laws for hovering micro-air-vehicles near hover. 5 Despite the widespread usage of the experimental aerodynamics model and the information that it can provide regarding insect aerodynamics, the theory is normally applied as a model operating at a single point. Placing the model in the context of perturbations from that operating point provides insight into the fundamental dynamic behavior, which can then be used to understand the sensing and feedback requirements for stable ight. These are precisely the goals of this research. This research examines the implication of passive aerodynamic stability mech- anisms associated with apping ight. Euler rigid body dynamics are paired with quasi steady aerodynamcs modeling that includes e ects of perturbations from the hover equilibrium. Results are based on analysis of the analytical equations as well as frequency-based system identi cation of the non-linear simulation. The goal is a set of linearized state-space models valid for small motions about hover that may be used to understand sensing and feedback requirements (and directly provide modal insight as in Taylor and Thomas (2003c); Sun and Xiong (2005)). Such models are derived under the fundamental assumption that it is the averaged forces and moments over the wingstroke that are important up to timescales of the rigid body dynamics, an assumption which is examined by direct comparison to a simulator includes the e ects of force and moment changes throughout a winstroke. Other key assumptions are the use of rigid (planar) wings that do not deform, rigid in- sect bodies ( xed abdomen and leg angles). Additionally, aerodynamic e ects such as reverse ow and spanwise ow are ignored as a side e ect of the experimental aerodynamics model used (the e ect of this assumption is examined in Chapter 5). 6 Results from fruit ies (Drosophila) and blow ies (Calliphoridae) are presented in this dissertation. Wing kinematics and aerodynamics measurements are available for fruit ies, while sensory system measurements are available for blow ies. From the perspective of understanding sensing and feedback requirements of insect ight, these measurements are essential. Though the example insects used in the simulation use fruit y and blow y parameters, the theoretical approach is derived for a general insect exhibiting a timescale separation in hover and the qualitative results are applicable to the translational e ects of dipteran apping wing ight. 1.3 Previous Modeling Work Previous research has attempted to generate models appropriate for estimating the dynamic properties and feedback requirements of dipteran ight. Experimen- tal work includes Taylor and Thomas (2003c,a,b), who used measurements taken of tethered locusts in forced air ow to quantify the locusts? response to body pitch angle and velocity perturbations, and more recent work by Hesselberg and Lehmann (2007) identifying body drag as a dominant e ect, as well as the experimental iden- ti cation of yaw damping by Dickson et al. (2010). High delity aerodynamic mod- eling has been applied for perturbations about trim to quantify the system dynamics of several diptera (crane ies, drone ies, and hawkmoths in hovering ight) through linear modeling (Sun et al., 2010; Sun and Xiong, 2005; Gao et al., 2009). Mod- els that may be developed without intensive computational solvers or experimental modeling are now of interest, as exempli ed by recent analytic modeling e orts 7 (Cheng et al., 2010). 1.4 Approach This research approaches the ight dynamics problem by quantifying the fol- lowing: 1. Wing Kinematics: Automated reconstruction of an insect lmed under high speed videography is used to measure wing motions of freely- ying insects, and stereotyped wing kinematics patterns determined. 2. Reduced Order Aerodynamics: Aerodynamics at the insect scale are a subject of active research. Rather than explicitly solve the Navier Stokes equa- tions, an experimentally derived reduced order aerodynamic model is posed in the context of insect egomotion (body translation and rotation). 3. Rigid Body Dynamics: The frequency separation between the insect body modes and high frequency excitation is exploited to allow classical Euler rigid body dynamics to be applied. 4. System Identi cation: The inputs and outputs of a nonlinear insect ight simulation built using the above components is analyzed in the frequency domain to assess linearity and measure experimental transfer functions. Be- ginning with analytical guesses, state space models are then tted to derive reduced order models. 5. High Fidelity Aerodynamics: The reduced order models are directly com- 8 pared to the result of a high delity numerical solution to the Navier Stokes aerodynamics equations. 1.5 Dissertation Contributions This dissertation presents the following contributions to insect ight: ? The rst detailed measurement of trim and maneuvering wing kinematics for Calliphorid ies. ? Implementation of planar wing tting and body exclusion rules to automated wing motion digitization. ? The rst quanti cation of full path measurement error in automated wing motion digitization routines. ? A de nition of biologically inspired control inputs that can be used for ight control. ? A method of placing the contemporary insect aerodynamics models in the context of insect body translation and rotation ? Transfer functions for longitudinal and lateral-directional hovering ight dy- namics of Drosophila. ? The indication that linear time invariant modeling is su cient to describe Drosophila ight dynamics about hover. 9 ? State space dynamics models representing Drosophila ight dynamics about hover. ? The use of a biologically inspired feedback path to determine that sensing of insect body angular rates is su cient for ight stabilization. ? An investigation of how Drosophila hovering ight dynamics are a ected by rate feedback and morphological changes. ? Longitudinal and lateral-directional transfer functions for Calliphorid ies ma- neuvering about forward (cruise) ight. ? Application of control theory to determine all reachable states for bounded inputs to wing kinematics programs. ? A method of ranking wing kinematic input programs in terms of maneuver- ability. 1.6 Dissertation Organization The organization of this dissertation is as follows. Chapter 2 presents an experimental method of determining the wing motions for dipteran (two-winged) insects. The method is an automated method that uses high speed videography of insects in free ight to generate time histories of the wing orientations. The digitization method is validated against at scale manufactured reference models, and the results are used to identify stereotypical kinematics for insect cruise and maneuvering ight. Digitized results are presented for Calliphorid ies. 10 Chapter 3 de nes biologically inspired control inputs for the measured wing kinematics, and applies these to published Drosophila kinematics. These inputs are paired with an experimentally derived hovering aerodynamics model placed in the context of insect egomotion. The aerodynamics are paired with rigid body dynamics and system identi cation performed to determine a longitudinal ight dynamics model for an example Drosophila-like insect maneuvering about hover. Uncertainty estimates and dynamic properties for the model are determined, as well as the e ect of mechanosensory feedback measured in previous work. Chapter 4 presents a similar study, but this time for the lateral-directional dynamics of the insect about hover. In both cases, the system is found to contain both stable subsidence modes and an unstable oscillatory mode that may be stabilized under the addition of mechanosensory feedback. Chapter 5 examines the e ect of additional un-modeled aerodynamic mech- anisms via comparison to a numerical solution of the Navier Stokes equations of motion. The dynamic properties show qualitative agreement (with translational agreement much better than rotational), but the unmodeled aerodynamics increase the speed of response, in both convergence and divergence. Chapter 6 further examines the e ects of feedback on both the reduced order model and the full aerodynamic solution, and studies how the sexual dimorphism prevalent in many species a ect the dynamic properties of the reduced order ight dynamics model. Feedback in the reduced order model shows a smooth progres- sion from bare airframe unstable to feedback stabilized, while the full aerodynamic solution shows that the unmodeled dynamics increase the feedback required for an 11 equivalent performance speci cation. Finally, the ight dynamics of the female and male contain the same dynamic modes, with small changes in the actual magnitudes, suggesting that ight stabilization requirements for the two sexes is not signi cantly di erent. Chapter 7 is an application of the reduced order model to compute and max- imize the maneuverability of a micro aerial vehicle. Control theory is applied to determine the reachable states under a given set of inputs, and a method of rank- ing kinematic programs in terms of maneuverability is presented. The ranking is applied to the hovering fruit y model developed in Chapter 3, showing dramatic di erences in the size of reachable space. The method is applicable to micro air vehi- cle design studies, where the reduction of actuation count and e ort is an important consideration Finally, Chapter 8 extends the reduced order modeling method developed for diptera in hover to forward ight, using the Calliphorid kinematics presented in Chapter 2. 12 Chapter 2 Wing Kinematics Measurement An essential component of a quantitative insect ight mechanics analysis is a detailed description of the wing motions an insect uses to e ect motion throughout its environment. This chapter is concerned with measuring and parameterizing motion throughout an insect?s wingstroke. The goal of the parameterization is to de ne intuitive control terms that may later be used as inputs for the insect ight control problem. 2.1 Previous Work Previous research into recording kinematics began with Vogel (1967), who began to record gross stroke parameter variations in tethered Drosophila, such as stroke plane angle or amplitude. Both Ellington (1984b) and Zanker and Gotz (1990) continued this work for several insects, using high speed photography with a single camera and a side view to record the wingtip trajectory as projected in the longitudinal frame. Later work by Gotz (1987) introduced a wingbeat analyzer, the working principle of which was a photodiode measuring light levels underneath the insect?s wings to infer both amplitude and frequency. Single camera analysis predominated kinematics recording for some time, such as Ennos (1989) work with dipterans (two-winged insects). Such an approach allows estimation of wingbeat 13 frequency, the amplitude of the wingstroke (as measured peak-to-peak), and the orientation of the plane1 in which the insect aps (not necessarily simultaneously). Digitization was primarily a manual and tedious process, though attempts were eventually made to use software to digitize wing markers and record wing steering muscle inputs as well (Balint and Dickinson, 2004; Sherman and Dickinson, 2003). Much of the work focused on tethered insects for ease of lming, though in some cases the motion of the tether was prescribed to determine the kinematic response to such motion (Sherman and Dickinson, 2003, 2004). Eventually, multiple cameras were used to allow determination of the wing orientation (assumed rigid) in three dimensional space (Fry et al., 2003). Detailed kinematic measurements enabled replay through robotic apparatus (Balint and Dick- inson, 2004). Once the robotic apparatus was equipped to record forces and moments generated by the motions, the rst basis of an aerodynamics model could be gen- erated. The manual three-camera setup became the standard method of recording kinematics, though modeling was limited to either the gross stroke parameters as above or directly replaying the measured kinematics through a robotic apparatus or simulator (Altshuler et al., 2005; Vance et al., 2005; Vance and Humbert, 2010; Vance et al., 2010). The next major upgrade was to apply an automated method to digitizing Drosophila kinematics (Ristroph et al., 2009; Fontaine et al., 2009; Ristroph et al., 2010). 1The actual ap pattern is rarely a plane, but a plane is typically de ned, using either the peak-to-peak wing positions or a linear regression to the wingstroke points. 14 2.2 Insect Subjects and Experimental Apparatus Fruit y (Drosophila) kinematics are available from previous literature (due to their ubiquitous presence in genetic research labs) (Dickinson et al., 1999), but de- tailed free ight kinematics for Calliphorid ies had not yet been recorded. Though fruit y wing kinematics and aerodynamics have previously been characterized, their visual and mechanical sensing capabilities are not yet well known, which would be necessary to understand the ight dynamics model developed in the context of insect sensing. In contrast, sensory feedback measurements for blow ies (Calliphoridae) are readily available, but kinematic and aerodynamic data are not. An experimental test rig shown in Fig 2.1 was constructed to allow detailed measurements of free ight maneuvering in Calliphorid species. The experimental apparatus consists of three Vision Research Phantom V710 high speed video cameras and lighting array, orthogonally-mounted about a 10in x 10in x 8in Plexi-Glass test section into which insects are introduced. A calibrated low speed wind tunnel has been constructed around the test section in order to allow steady wind conditions to be applied, and a regulated supply of compressed air allows repeatable gusts to be introduced through a port in the chamber as well. Calliphorid insects captured at the University of Maryland horse farm were allowed to freely y within the chamber, and a trigger mechanism was built to automatically record sequences when an insect ew within the focal volume. Images were recorded at a resolution of 1280 x 800 pixels with a 30 s exposure and 7002 Hz framerate. A library of Calliphorid ight sequences incorporating maneuvering free ight 15 Figure 2.1: Experimental test rig and free ight chamber. behavior in both quiescent air and gusting conditions has been collected. Digitiza- tion is traditionally accomplished by manually aligning a wireframe model of the insect to each of the three views, a labor-intensive task that can take weeks or months for a single sequence. Research on Drosophila ight at Cornell University has demonstrated success using automated tracking techniques such as Hull Recon- struction Motion Tracking (HRMT) (Ristroph et al., 2009), and this study used a modi ed version of the Drosophila tracking code. 2.3 Kinematics Digitization 2.3.1 Coordinate De nitions The description of the insect apping motion requires a family of axes centered at the insect wing hinge. Approximating the wings as rigid bodies, measured insect kinematics exhibit a roughly planar ap motion which is represented using 2-3-2 Euler angles. De ne by reference to Fig. 2.2a a set of stability axes S = fs^x; s^y; s^zg passing through the insect center of mass G, the stroke plane angle as the angle 16 about the pitch axis to an idealized planar stroke motion, and a coordinate axes set aligned with this plane the stroke plane axes P = p^x; p^y; p^z . De ne R = fr^x; r^y; r^zg a set of axes that move along with the right wing, with r^z = p^z and r^y to extend toward the wing tip as in Fig. 2.2b. Similarly, de ne L = n l^x; l^y; l^z o for the left wing, with l^y extending inboard along the left wing spanline. The additional de nition of the geometric angle with respect to the stroke plane as wing pitch angle provides the notation necessary to describe the orientation of two rigid wings at an instant in time. G ? ? s? x s? z p? x p? z (a) Stroke plane axes/angle , body hovering angle Gs? y p? y p? x r? y ? r r? x (b) Stroke angle r and R axes Figure 2.2: Axes and angle de nitions. 2.3.2 Digitization Algorithm Once a ight sequence has been captured, kinematics digitization consists of three major components: image processing, volume pixel (voxel) extraction, and feature identi cation/model tting. In the image processing stage, the y is isolated from the background and 17 noise in a cropped version of the image (400 by 400). The techniques used in this are background subtraction and masking by a thresholded version of the image. Each of the images is then corrected for alignment and magni cation errors to create orthogonal pro les of the animal. In the voxel identi cation stage, a search space is identi ed for the insect, and the pixels in each camera are compared to determine the three dimensional volume obscured by the y in the capture volume. Mathematically, this approach is analogous to a computer-aided-drawing program taking the intersection of each of the animal?s orthogonal pro les to reconstruct its volume. In the model tting stage, the voxels are grouped into head, body, and wing groups by an iterative distance-minimizing approach. The head and body groups are combined and the centroid of the group used to estimate the body center of mass location. The body orientation is determined by the identi ed semi-minor and semi- major axes of a principal component analysis conducted on the body voxels. The wing orientations are estimated by rst isolating the wings from their hinges, then conducting principal component analysis and a planar t to the outboard section of the wing voxels. The above three steps begin with three-view images like those shown in Fig 2.3(a-c) and generate as output volume reconstructions like those shown in Fig 2.3(d-f) and Fig 2.4. 18 (d) top (e) front (f) side Figure 2.3: Example images versus voxel reconstructions. 40 45 50 55 60 65 70 75 5 10 15 20 25 120 125 130 135 140 y x z Figure 2.4: Visual hull reconstruction. 19 (a) Straight/Level (b) Left turn (c) Climb (d) Yaw Figure 2.5: Digitized ight sequences include a (a) straight and level sequence, (b) level left turn, (c) climb initiation (d) yaw motion (sideslip). All images are top views shown at 70Hz except for (c), which is a side view shown at 140Hz. 2.3.3 Digitization Results Digitization was conducted for ight sequences including straight and level ight, a level left turn, initiation of climbing ight, and a yaw motion (saccade) (Bender and Dickinson, 2006b), as seen in Fig 2.5. All of the ight sequences start in forward ight and are in quiescent air. Principal component analysis works quite well for stroke and deviation angle identi cation, shown in Fig 2.6. A number of modi cations were necessary to the automated HRMT approach in order to yield accurate results for Calliphorid wing pitch angles. Previously, wing pitch angle was unreliable, due to body voxels that 20 0 1 2 3 4 5 6 7-90 -60 -30 0 30 60 90 Wingstroke St ro ke an gl e, de g ? r Fit Figure 2.6: Stroke angle is well detected by principal component analysis, and mul- tiple wingstrokes in straight and level ight provide ample measurement of the nominal kinematic pattern. were assigned to the wing grouping and errors in the method of determining a chord vector. A spherical body exclusion rule dramatically reduced the body pixels incor- rectly assigned to wings, as seen in Fig 2.7. The wing angle detection method was also revised. Wing angle estimation previously estimated the wing chord vector from points having the most distance to each other. Instead, each wing is now modeled as a plane and the algorithm computes a least-squares best t to the function x = b0 + b1y + b2z: (2.1) Both modi cations markedly improved the wing angle detection, as seen in Fig 2.8. 21 50 60 70 80 30 35 40 45 50 105 110 115 120 125 130 135 140 x y z Figure 2.7: A spherical exclusion rule dramatically reduces the number of body vox- els incorrectly assigned to the wings. Blue indicates body voxels that were previously detected as wing voxels. 0 0.005 0.01 0.015 0.02 0 W in g 0 0.005 0.01 0.015 0.020 45 In iti al 0 0.005 0.01 0.015 0.020 45 Time, sec Pl an ar Right Left 90 ?90 90 90 ? ? ? Figure 2.8: Modi cations to the wing angle detection algorithm improved the results dramatically. 22 10 20 30 40 50 60 70 30 40 50 60 70 10 20 30 40 50 60 70 80 y x z Figure 2.9: A reconstruction of the 45 insect model reconstruction shows an occlusal defect on the left wing. 2.3.4 Digitization Accuracy Previous work has estimated accuracy of a digitizing method using synthetic views of a computer graphics insect (Ristroph et al., 2009). A series of rapid- prototype insect models with wings in known positions were manufactured and ac- curacy of the technique estimated by placing the model in the capture volume and digitizing the resulting images. The advantage of such a technique over a typical calibration with synthetic images is that it provides a very realistic estimate of the error in the entire measurement system, including imaging errors, lens distortion, camera misalignment. These e ects are not quanti ed with a synthetic image cali- bration, which provides an error estimate that only includes the digitization method. A calibration sequence showing an occlusion (a portion of space obscured in all three cameras) is shown in Fig 2.8. The occlusal defect on the left wing raises the error in the wing angles from a mean of 7 to a mean of 8:3 . 23 0 0.2 0.4 0.6 0.8 1?70 ?30 0 30 60 90 Time, t/T Angle, de g ? r ? r TSA ? r TSA ?g Figure 2.10: Time synchronous averaging is used to identify a nominal stroke pattern associated with a reference condition (straight and level ight shown). 2.4 Kinematics Reduction 2.4.1 Time Synchronous Averaging For insects in an established ight condition, repeated wingstrokes (see Fig 2.6) o er a method of improving the wing kinematics angle estimation by o ering multiple measurements of the nominal wingstroke. This improvement was realized through time synchronous averaging, a method (often used for geartrain diagnosis) which aligns each waveform of a periodic signal and averages each wingstroke at each point during the wingstroke (Dalpiaz et al., 2000). Time synchronous averaging allowed curve ts to be applied to a single stereotypical curve t with a standard error of less than 2 , as seen in Fig 2.10. 24 2.4.2 Curve Fitting Previous Drosophila kinematics work had assumed a planar stroke, sinusoidal wing motion in , and a smoothed square wave for g. The assumption of planar stroke was retained for Calliphora, but both stroke angle and wing angle showed deviations from this pattern, as seen in Fig 2.10, and sinusoidal curve ts often resulted in a t with an unacceptably low coe cient of determination R2. Adding a \stretch parameter" d to the sinusoid as (t) = cos(2 ft+ 1)(d cos( ft+ 2) 2 + 1) + o (2.2) and modeling as a combination of a sinusoid and smoothed square wave (t) = A1[:5 tanh(b1 sin[2 ft+ e]) + :48] sin(2 ft+ + e) 2 +A2[:5 tanh(b2 sin[2 ft+ c]) + :48] tanh(3 sin(2 ft+ c)) + d (2.3) improved t quality to yield R2 > 0:9. Wing stroke angle and wing pitch angle ts are shown in Fig 2.6 and 2.11. 2.5 Summary In this chapter, an automated approach was used to extract kinematics from high speed videography. A body exclusion rule and planar wing tting was added to improve the quality of the wing orientation estimates. The method was validated us- ing manufactured reference models, representing the rst published uncertainty data 25 0 0.2 0.4 0.6 0.8 10 30 60 90 Wingstroke Wing pitch angle, de g ? r Fit Figure 2.11: Wing pitch angle in forward ight is more complex than a smoothed square wave. for three-camera insect wing orientation measurement to include the full extraction path. Trim (stabilized forward ight) and maneuvering kinematics were identi ed for Calliphorid ies, a genus for which detailed kinematic data was previously un- available. Time synchronous averaging was used to determine stereotyped reference kinematics, and curve ts were populated with biologically inspired kinematic pa- rameters which are used as inputs to a reduced order apping ight dynamics model developed in the remainder of this dissertation. 26 Chapter 3 Longitudinal Hover Dynamics 3.1 Introduction and Background We rst analyze longitudinal apping ight mechanics about hover, using as an example Drosophila melanogaster, or the common fruit y. Fruit ies are abun- dant in research settings, owing to their rapid lifecycle that allows tractable genetic studies. Correspondingly, they have one of the most well-understood genetic struc- ture amongst insects. For the goal of ight dynamics research, they are easy to keep and induce to y, so that both hovering kinematics and experimental aerody- namic measurements may be found in published literature (Dickinson et al., 1999). Additionally, there is a large phenotypic variation amongst di ering communities and for di ering sexes. This variation will enable the study of the e ects of large morphological variations in Chapter 6. This section presents the hovering Drosophila wing kinematics used in sim- ulation, a review of the experimental aerodynamics model, and other governing equations used for analysis and simulation. 3.1.1 Wing Kinematics Section 2.3.1 de ned wing stroke plane inclination , stroke angle , and a wing pitch angle g for each wing, and introduced the coordinate frames used in 27 this study, and Fig 3.1 is repeated here for clarity. For a apping Drosophila in or G ? ? s? x s? z p? x p? z (a) Stroke plane axes/angle , body hovering angle Gs? y p? y p? x r? y ? r r? x (b) Stroke angle r and R axes Figure 3.1: Axes and angle de nitions. near hover, the ap angle r undergoes a harmonic motion su ciently represented as a sinusoid (Dickinson et al., 1999) r(t) = r cos(2 frt) + o ,r ; (3.1) where r gives the amplitude of each wingstroke, o ,r the deviation of the point about which the wing oscillates, and fr the right wing ap frequency. The geometric angle of attack g exhibits a harmonic motion roughly resembling a modi ed square wave, which allows the advancing and retreating strokes to both share a positive angle of attack. 28 3.1.2 Control Parameters As postulated by Vogel (1967) and later quanti ed experimentally by Fry et al. (2003), insects modulate the time forces and moments applied to wingstrokes by modi cation of several wingstroke parameters. The control inputs considered in this chapter are the biologically-motivated choice of ap frequency f in Hz, the ap amplitude as de ned in section 3.1.1, stroke plane angle , and the mean position (center) of wing oscillation o . In addition to the mathematical de nitions in section 3.1.1, several control parameters may be seen graphically in Fig. 3.2. Right and left wing variables may be interpreted as longitudinal and lateral control inputs using the transformation 2 6 6 6 6 6 6 6 6 6 6 4 c d c d 3 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 4 0:5 0:5 0 0 0:5 0:5 0 0 0 0 0:5 0:5 0 0 0:5 0:5 3 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4 r l r l 3 7 7 7 7 7 7 7 7 7 7 5 ; (3.2) where [ ]c represents a collective (symmetric) input and [ ]d represents a di erential (lateral-directional) input. Parameter variations are generally small, and the varia- tions remain remarkably small even for aggressive maneuvers such as fast 90 colli- sion avoidance maneuvers known as saccades (Fry et al., 2003), but are nonetheless fundamental for the control of the insect. As an example, Table 3.1 and 3.2 shows nominal and maneuvering parameters for Calliphorid ies maneuvering throughout the sequences in Fig 2.5. Since this chapter addresses longitudinal motion, we de- 29 Sequence ( r; l) ( c; d) S+L (19:5; 21:4) (20:5; 0:9) Turn (15:3; 24:0) (19:7; 4:4) Table 3.1: Stroke plane inclination in a turning Calliphorid shows di erential stroke plane inclination perturbation and constant collective stroke plane inclination . Parameter Collective Di erential Frequency 211-215 Hz 0 Hz Amplitude 107-112 0.8-21 O set 0.3-6 16-17 Table 3.2: Range of parameter variation observed in Calliphorid ies throughout maneuvering ight. scribe the behavior of the right wing and assume symmetry in the left wing frame and the subscripts [ ] may be suppressed until Chapter 4. (a) Stroke plane angle input (b) Stroke o set o and amplitude Figure 3.2: Insect apping kinematics parameters are used as control inputs to model. 3.1.3 Aerodynamics A variety of e ects, predominantly unsteady, are known to be active during an insect?s ight. A thorough treatment of these e ects is outside the scope of this section and may be found in Ansari et al. (2006); Sane (2003). Instead, this 30 treatment reviews the largest contribution to in- ight insect forces: \translational" lift and drag. 3.1.3.1 Translational Lift and Drag Wing \translational lift" is the largest component (approximately 65-85%) of an insect?s lift production in hover and the most straightforward of the lift mecha- nisms known to be active, but includes a number of unsteady e ects via experimen- tal coe cients. The translational component of insect lift can be represented using (Ellington, 1984a) L(t) = 1 2 Sjut(t)j 2r^22CL[ (t)]; (3.3) where the instantaneous lift force L is written as a function of the air density , the wing area S, tip velocity ut, nondimensional second moment of area r^22, and an experimentally-determined lift curve slope (Sane and Dickinson, 2002) CL[ (t)] = 0:225 + 1:58 sin(2:13 g 7:2): (3.4) The second moment of area r^22 may be de ned in terms of the normalized chord c^ = 1=2A c=R and normalized radius r^ = r=R as r^22 = R 1 0 c^ r^ 2 dr^ (Ellington, 1984c). This model is referred to in the insect ight community as \quasi-steady" because the only time dependent variables in this equation form are the kinematic variables ut(t) and (t), yet we know that insect aerodynamics are intrinsically unsteady (Dickinson and Gotz, 1999). Consequently, the nondimensional CL[ (t)] 31 is a curve t to experimental data that hides a number of unsteady e ects, such as the e ect of starting and stopping vortices (Wagner e ect) (Dickinson, 1996) and the mechanism of delayed stall. Though these are unsteady e ects, the CL[ ] curve t is considered to be a function purely of the instantaneous angle of attack, hence the equations appear in the same form as an aeroelastic analysis that neglects such e ects. Because of this, an aerodynamic model presented in quasi-steady form may either neglect unsteady e ects or treat them via modi cation of the CL[ (t)] as done in this dissertation. The other nondimensional term, r^22, addresses the nonuniform spanwise lift distribution and is a function of the wing?s shape (Ellington, 1984c). A key feature of the aerodynamic model is the reduction of forces by the dynamic pressure 12 u 2 t , where the tip velocity ut = _ R is a function of the wing ap angle derivative _ and the wing length R. Translational drag is the component of force production acting opposite the direction of wing motion (in the stroke plane) and has a parallel representation to the translational lift described above. Introducing the experimentally-determined nondimensional coe cient CD[ ] to describe this force, one may model translational drag as (Dickinson et al., 1999) D(t) = 1 2 Sjut(t)j 2r^22CD[ (t)]; (3.5) CD[ (t)] = 1:92 1:55 cos (2:04 g 9:82) : (3.6) 32 3.1.4 Wingstroke-Averaged Forces and Moments Flight dynamics often involves a system (airframe) forced by high frequency forces, creating small high frequency motions (vibration) which are typically ignored and gross translation/rotation (the ight path). The solution of interest is the ight path, or asymptotic solution, the speed of which is limited by the airframe inertia. By reference to results from \two-timing" averaging theory (Cole, 1968), one can represent the asymptotic solution of a system whose response is on slow timescales with periodic excitation on fast timescales as the e ect of the excitation averaged over the fast timescale, subject to speci c constraints on the system (Kevorkian, 1966; Cole, 1968), an approach that works well for systems with body dynamics slower than 1/10th the dynamics of the forcing function (Deng et al., 2006a). Deng et al. applied this concept to the design of ight control for a mechanical apper (Deng et al., 2006b; Schenato, 2004) by viewing the averaging process as a nonlinear map from wrench e orts to state outputs. In this study, we investigate the details of creating a map from kinematic input parameters (controls) to state outputs and the dynamic properties of the resulting system, determining the modal (asymptotic) airframe response time scales via computed frequency responses. Under the assumption that the aerodynamic timescales are an order of magni- tude faster than the body dynamics, the wing aerodynamic forces are averaged over one wingstroke to determine the motion of the insect body. Moreover, nding the ef- fective force over each wingstroke provides a mathematically straightforward means to impose the constraint of the control variables being xed over each wingstroke. 33 Such a constraint is useful because insects, unlike larger animals such as birds or y- ing mammals, do not posses the computational bandwidth to control wing motions more than once per wingstroke (Autrum, 1958). The e ective, or mean, lift force over one wingstroke is Lavg(t) = 1 T Z T 0 L(t)dt; (3.7) where T is the period of the wingstroke given by T = 1=f . An analytic solution to this lift integral would provide an algebraic expression for the e ective lift force acting over a wingstroke, dependent on: the wing geometry (size and shape), the environment (density), the wing aerodynamic performance (via lift curve slope), the longitudinal state variables, and a subset of the control inputs. The general inte- gration is a numerical task due to the numerous \nested" trigonometric functions. Restricted solutions to the integral may be found by recognition that the lineariza- tion performed later does not require a global solution to this integral, but rather the exact solution at a point and its functional dependence on the state and wing ap coordinates (partial derivatives). The drag equation must similarly be averaged over one wingstroke. 3.2 Derivation of Linearized Flight Dynamics about Hover In this section, the experimental aerodynamics model is extended to include perturbation velocities and the local and instantaneous forces are approximated by wingstroke-averaged forces to provide estimates of the vehicle stability and control 34 derivatives. 3.2.1 Governing Equations Assuming a rigid insect body B, the motion of the body- xed stability frame S is a function of the applied forces and moments and can be written as a system of ordinary di erential equations (Nelson, 1989). The longitudinal portion of these equations are X = m ( _u+ qw rv) +mg sin( ) (3.8a) Z = m ( _w + pv qu) mg cos( ) cos( ) (3.8b) M = Iyy _q Ixz( _r + pq) (Iyy Izz) qr; (3.8c) where X, Z, and M are the aerodynamics forces and moment, !B = ps^x+ qs^y + rs^z the rotation rate of frame S, and vG = us^x + vs^y + ws^z the inertial velocity of the center of mass G as seen in Fig. 3.3. 3.2.2 Hover Equilibrium Trim Solutions Using the wingstroke averaged forces and moments, under the assumption of longitudinal motion (p0 = 0, r0 = 0,v0 = 0, 0 = 0), and writing each variable as a nominal condition (notated as [ ]0 and a small perturbation [ ], the general equation for motion along the sx axis may be written for hover (where the nominal 35 G s? z s? x s? z Y, v Z , w X , u L, p N , r M , q Figure 3.3: De nitions of stability frame velocities vG = us^x + vs^y + ws^z, rotation rates !B = ps^x + qs^y + rs^z, forces F = Xs^x + Y s^y + Zs^z, and moments M = Ls^x +Ms^y +Ns^z. w0 = 0 and q0 = 0) as X0 + X mg [sin( 0) + cos( 0) sin( )] = m [ _u] : (3.9) At the trim condition, the small perturbation terms [ ] = 0, so eqn (3.9) becomes X0 = mg sin( 0): (3.10) Calculation of a trim condition requires speci cation of a nominal hovering trim angle 0. Without restriction, we can de ne the stability axes such that 0 = 0. To do so, de ne the pitch angle through which the body axes (aligned with the insects longitudinal morphological axis) must be pitched through in order to become coincident with the stability axes as the \body hovering angle" (shown in Fig. 3.1), then project the moment of inertia tensor through the single axis rotation, and de ne the distance from G to thorax as rG=t = dxs^x+dzs^z, where a symmetrically-located 36 mass center is chosen to preserve longitudinal motion. The hover condition along s^x reduces to X0 = 0. A similar derivation of the hover condition along s^z and about s^y gives Z0 = mg and M0 = 0. Hover equilibrium now reveals multiple solutions for trim inputs. In unac- celerated hover with a level stroke plane, vertical force equilibrium requires that Z0 = 2Lavg(t) = mg, which may be trimmed via either ap frequency or stroke am- plitude inputs. Moment equilibrium reveals that an insect similarly has redundant control inputs via stroke plane angle, stroke angle o set, or a coupled input of both. For a given stroke plane angle input 0, the insect may enforce pitch axis trim via o ,0 = sin 1 sin 0 dz cos 0 dx ra cos(2 0) ; (3.11) where ra is the spanwise distance outboard to the point at which the lift and drag forces may be considered to act. Eqn (3.11) demonstrates that a solution can only exist for 0 6= (2n 1) 4 , n = 1; 2; 3:::. The dependence on body hovering angle may be seen in Fig. 3.4. A detailed estimation of the distance ra, which is involved only in the pitch dynamics, may be addressed via computational uid dynamics studies and is outside the scope of this investigation. The point has been shown experimentally to be roughly constant with respect to changes in angle of attack (Dickson et al., 2008). Initially, the centroid may seem a wise choice, but since lift forces are proportional to the square of the local velocity, the actual location of this point is more subtle. Application of blade element theory with a xed spanwise blade pitch yields ra = 37 ?30 ?20 ?10 0 10 20 30 ?30 ?25 ?20 ?15 ?10 ?5 0 5 10 15 20 Trim stroke offset , de g Trim stroke plane angle , deg =60? =75? =90? ?0 ? o ? ,0 ? ? ? Figure 3.4: Trim inputs 0 and o ,0 as functions of body hovering angle . 3=4R, a value commonly used in helicopter results (and the reason why the twist of a linearly twisted helicopter blade is speci ed at this point), which has been used as a rst approximation for this point. 3.2.3 Inclusion of Perturbation States in the Quasi-Steady Form The current practice of posing the aerodynamics model without any state dependency is useful from a data reduction standpoint but does not provide insight into how perturbations of each of the state variables from the operating point a ect the aerodynamic forces and moments produced on the vehicle. The recent discovery of \ apping counter-torque" (Hedrick et al., 2009) as a means to arrest a turning maneuver is a result of the consideration of perturbation velocities for yaw rotation, which is now examined for the general longitudinal case. We consider inclusion of the perturbation velocities vG = us^x + vs^y + ws^z, 38 and rotation !B = ps^x + qs^y + rs^z. Intuitively, one would expect no aerodynamic dependency on inertial position or orientation.1 Both vG and !B have two major e ects on the wing lift and drag production: (a) they change the local velocity (and hence dynamic pressure) at an airfoil section, and (b) they also change the local angle of attack at each airfoil section. 3.2.4 Velocity Components In order to determine the incremental lift and drag acting on an elemental wing section, one must consider the total velocity incident that wing, de ned as v = v +vG+v!, where v , vG, and v! are the velocity components due to apping, body translational rate, and body rotational rate, respectively. To facilitate numerical simulation, the components are resolved in the wing frame as v = vxr^x+vyr^y+vz r^z. 3.2.4.1 Local Velocity Component due to Flapping The component of velocity due to wing apping at a distance r along the wing spanline is _ r acting along the r^x axis or v = _ r r^x. The original model evaluates this function only at the wing tip for ut = _ R. 1Note that aircraft are occasionally modeled (in reduced order models used for control design) with a dependency on roll angle that arises due to the sideslip velocity v ensuing from such an angle, as a estimate is often more readily available from an instrumentation system and is then used to e ectively estimate v. Theoretically, a formulation including v is well-posed. 39 3.2.4.2 Local Velocity Component due to Body Velocity In the rotating right wing frame R, the body translation velocity vG has the components vG R = urr^x+vrr^y +wrr^z, related by the rotation matrix RRS de ned in eqn (3.13). vG R = RRS vG S (3.12) RRS = 2 6 6 6 6 6 6 4 c rc s r c rs s rc c s rs s 0 c 3 7 7 7 7 7 7 5 (3.13) Similarly, the velocities seen on the left wing are given by RLS, the compliment of RRS. Though a function of the ap angle r and stroke plane angle , vG R is constant with respect to wing spanwise distance. 3.2.4.3 Local Velocity Component due to Rotation Rates The body rotation rate !B enters in a slightly more complicated fashion. The velocity of a point P due to rotation rate !B about the center of mass G can be expressed as v! = !B rG=p, where rG=p is the vector from G to P . Unlike the velocity components due to body translation, the velocity components due to body rotation are functions of the distance r outboard. The choice of stability axes as a basis for ! requires a transformation RRS to express it in the wing frame for numerical computation. 40 3.2.5 E ects of Modi ed Velocity The original aerodynamics model (Sane and Dickinson, 2002) made the as- sumption that vG and v! were both negligible in comparison to the apping veloc- ity v , The inclusion of the additional components formulated as a modi cation of the wing?s 2D angle of attack and local wing frame velocity forms the basis for the dipteran apping wing dynamics analysis. 3.2.5.1 Angle of Attack Modi cation The local velocity components have a marked e ect on the local angle of attack on each wing?s airfoil section. The change in the 2D angle of attack may be written as add = arctan vz vx : (3.14) The e ect of the angle of attack modi cation is shown graphically in Fig. 3.5a and 3.5b, which demonstrates that an imposed (positive) heave velocity (descent) in 3.5b corresponds to an increase in angle of attack of the wings and corresponding increase in lift acting to oppose the motion, while the converse is also true. Then, heave velocity damping is found by consideration of the perturbation velocities, in much the same manner as apping countertorque provides yaw damping via di erential velocities. 41 ? g v x (a) Hover ? g v x ? a dd v v z (b) Descent Figure 3.5: Wing angle of attack is increased in a descent perturbation and reduced in a climb perturbation, so a apping wing system in axial climb/descent exhibits damping in heave velocity w. 3.2.5.2 E ect on Dynamic Pressure Modi cation of velocity components also a ect the dynamic pressure 12 v 2 at each location. For example, any forward speed u results in an increase in dynamic pressure during the advancing stroke and a movement outboard of the location of zero dynamic pressure during the retreating stroke. Moreover, the wing component velocities are functions of both r and and consequently this spanwise location shifts throughout the wingstroke. The total dynamic pressure under addition of velocity components is qtotal = 1 2 v2x + v 2 y + v 2 z : (3.15) For the ight conditions of many dipteran insects (Vogel, 1967), the dynamic pres- sure in eqn (3.15) is dramatically a ected, up to several times the nominal dynamic pressure calculated using the wing kinematics, while the angle of attack as calculated by eqn (3.14) is a ected by several degrees. 42 3.2.6 Perturbation Equations Having found trim conditions, one may subtract eqn (3.10) from eqn (3.9) and solve for _u to obtain _u = X m g cos( 0) ; (3.16) where X has the physical interpretation of being the perturbation force due to state or control perturbations from equilibrium (trim) values. For traditional linearized analysis, a control model is developed by separating out the linear e ects of each of the longitudinal variables u, w, , q, as X m = Xu u+Xw w +X +Xq q + Xc m ; (3.17) where X[ ] is de ned to be X[ ] = 1m @X @[ ] . Similarly, the control term Xc may be expressed in terms of each of the input controls as Xc m = Xf f +X +X o o : (3.18) By expressing the perturbation forces in the form, a (time-invariant) linear sys- tem can now be written in standard state space form _x = Ax + Bu where x = 43 [ u w q ] T , u = [ f o ] T , and A = 2 6 6 6 6 6 6 6 6 6 6 6 6 4 Xu 0 0 g Zu Zw 0 0 Mu Mw Mq 0 0 0 1 0 3 7 7 7 7 7 7 7 7 7 7 7 7 5 , B = 2 6 6 6 6 6 6 6 6 6 6 4 0 0 X X ,o Zf Z 0 0 0 0 M M ,o 0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 5 . (3.19) 3.2.7 Homogeneous System Dynamics (the A matrix) Basic dynamics analysis begins with the unforced di erential equation _x = Ax. Without consideration of the perturbation velocities, di erentiation with respect to u, w, and q of the e ective X term yields that the stability terms Xu, Xw, and Xq are identically zero for any ight con guration because the wing?s lift and drag components are functions of the kinematics, wing shape, and environment, not the state variables. All higher derivatives also do not contain any dependence on the state variables so the conclusion is not limited to a linear analysis. The experimental aerodynamics model in absence of perturbation velocities suggests that some form of feedback is necessary for an insect to possess stability, which would signi cantly increase the computational workload the insect must support. 3.2.8 Control Derivatives (the B matrix) The insect aerodynamic forces and moments at hover equilibrium depend ex- plicitly on the wing kinematic functions, so the control derivative estimates may 44 be found analytically. Intuitively, for perturbations about the nominal hover case 0 = 0, the magnitude of the lift vector (via f or ) does not a ect the X force while a change in direction of the lift vector via does. Mathematically: as the X aerodynamic force contains a sin 0 term and 0 = 0 for the hovering insect, the _u equation has a nonzero X term shown in eqn (3.20), which evaluates to g because Z force equilibrium for the trim condition has already speci ed that the nominal value of the lift vector must be Z = W = mgc 0 , and X = 1 m @X @ , where c 0 = cos and s 0 = sin . X = 1 m @ @ 2 2 (r^22SR 2)(f 2 2 sin )CL;avg = 2 2 (r^22SR 2)(f 20 2 0)CL;avg m = gc 0 (3.20) 3.2.8.1 Heave and Pitch Dynamics The derivation has focused on the s^x axis, but the heave s^z and pitch dy- namics parallel the s^x results. The control input matrix B found by this method is shown in eqn (3.21), where the nominal integral lift L0 has been de ned in eqn (3.22). B = 2 6 6 6 6 6 6 6 6 6 6 4 2g c 0 f0 tan( 0) 2g c 0 0 tan( 0) gc 0 0 2g c 0 f0 2g c 0 0 gc 0 tan( 0) 0 0 0 2Iyy L0 dz 2 Iyy L0 cos h sin 1( dxra ) i 0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 5 (3.21) L0 = 2(r^22SR 2)(f 20 2 0)Cl (3.22) 45 For the biologically-motivated choice of a level stroke plane angle 0 = 0 and hov- ering at 0 = 0, eqn (3.21) can be simplifed to B = 2 6 6 6 6 6 6 6 6 6 6 4 0 0 g 0 2g f0 2g 0 0 0 0 0 2IyyL0dz 2 Iyy L0 cos h sin 1( dxra ) i 0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 5 : (3.23) 3.2.8.2 Signi cance of the Control Matrix B The control matrix B provides direct insight into the physics of how control inputs a ect insect ight motion. In particular, the control matrix exhibits a decou- pling which both suggests a means of control and is consistent with the biological intuition formed by dipteran insect observations. The rst row implies that an insect in hover uses its stroke plane angle to modulate its fore/aft movements, necessitating a coupled adjustment in the input o to maintain equilibrium. This coupled behavior is exactly as observed in free ight insects (Taylor, 2001). Moreover, by appropriate choice of wing length, the insect may lose the use of the input o . Initially it appears that extending/trimming an insect?s wings (to a ect ra) could force an insect to hover at a di ering pitch angle and that the insect would lose the ability to generate pitch moments at a particular wing length. However, since ra is normally several times larger than dx, the insect?s wings would have to be drastically shortened to achieve the dx = ra condition and the lift generated by the trimmed wings would be severely impaired before the 46 capability to generate pitch moments is destroyed. Forward Velocity u: The control a ecting the forward velocity u is the inclination of the stroke plane angle, suggesting a \helicopter-like" mode of regulat- ing forward velocity is physically reasonable. The X term re ects inclination of the lift vector{the g term is related to the trim input solution for equilibrium in hover. Note there is some cross-coupling into the pitch rate via nonzero M . Heave velocity w: A hovering insect has at least two means of controlling its vertical velocity w, both of which are related to changing the magnitude of the lift vector: wingbeat frequency f and ap amplitude . Both of these terms enter in a complementary fashion, and show the expected e ect: an increase in frequency or amplitude causes the insect to rise. Pitch angle : No direct control over the pitch axis is provided; instead control of the pitch is through integration of the pitch rate, as a direct consequence of Eulerian mechanics. Pitch Rate q: Pitch rate is primarily controlled by the o set (mean) term in the harmonic ap function. As expected, de ecting the mean position of the wings forward results in a positive (nose-up) pitch rate. Moreover, the strength of the term is related to the ratio of the distance along the s^x from thorax to the center of mass as compared to the wing length R. Among other things, this provides a means to immediately estimate via a photograph the pitch control authority a particular insect has, much like a xed-wing aircraft?s wing loading gives an immediate means to estimate its velocity. The simple insect pitch control estimate is particularly useful for small insects where characterizing the insect?s pitch inertia is di cult. Pitch 47 rate is also secondarily controlled by the stroke plane angle term which primarily generates a response in the forward velocity, so the o set term is an appropriate choice for controlling pitch rate and is a trend observed in nature (Taylor, 2001). 3.3 Simulation and System Identi cation In this section, a simulator is created using the perturbation velocities aero- dynamics model developed in Section 3.2 and a frequency-based model parameter identi cation approach for a simulated insect parameters described. 3.3.1 Simulation A simulation environment seen in Fig. 3.6 has been created including the aero- dynamics model developed in Section 3.2 and the 6 degree-of-freedom rigid body equations of motion (the longitudinal portions of which are seen in eqn (3.8a) to eqn (3.8c)) in order to capture the most essential rigid body motion of the insect forced by translational lift dynamics. The aerodynamic model was derived using Drosophila data, so the simulator is populated with parameters from this species. Wing mass is an order of magnitude smaller than body mass, and research by Sun and Xiong (2005) and Taylor and Thomas (2003c) has indicated that the e ects of high-frequency wing mass oscillations are small. Accordingly, the simulator repre- sents the inertial properties of the overall insect by the insects body mass and inertia tensor matrix, and body forces with a constant gravitational eld and traditional aerodynamics modeling (Full and Koehl, 1992). 48 K in em a tic s En gi n e w in g m o tio n W in g Ae ro dy n am ic s R ig id B o dy Eq u at io n s o f M ot io n fo rc es , m o m en ts H al te re M od el in g fe ed ba ck ? co n tr o l in pu ts A B B o dy Fo rc es? st at e K in em a tic s En gi n e W in g Ae ro dy n am ic s R ig id B o dy Eq u at io n s o f M ot io n H al te re M od el in g K in em a tic s En gi n e W in g Ae ro dy n am ic s R ig id B o dy Eq u at io n s o f M ot io n Bo dy Fo rc es H a lte re M o de lin g K in em a tic s En gi n e W in g A er o dy n a m ic s R ig id Bo dy Eq u a tio n s o f M o tio n Figure 3.6: Simulator used for haltere-on and bare-airframe input/output studies, showing (highlighted path) the critical modi cation: consideration of rigid body state in wing aerodynamics functions. The simulation includes synthesis of physically realistic wing kinematics tra- jectories, using the bio-inspired analytic expressions presented earlier to solve for the position of a wing at simulation time. The kinematics engine?s bio-inspired pa- rameters are then viewed as the control inputs and used to drive the insect motion as desired. The wing used in the simulation is a symmetric planar wing with a planform characterized nondimensionally via a second-order shape model. The use of such a wing is justi ed analytically and by reference to experimental data and high-speed imagery (Fry et al., 2003). (Flexibility e ects in the species studied are an active area of research and not included in the simulation at this time but will be included when signi cant results appear.) 3.3.1.1 Haltere Modeling Experimental evidence has repeatedly demonstrated that components of the insect?s control system, particularly the portions concerning vision, do not have 49 bandwidths that would allow them to actively stabilize high-frequency behavior aris- ing from the wing. Instead, the higher frequency angular rate loops have been closed using haltere-based rate feedback. Halteres are a pair of oscillating aerodynamically- ine ective \hind wings" that have been shown to be theoretically (Thompson et al., 2008) and experimentally (Derham, 1714; G Nalbach, 1993) capable of acting as a biological rate gyro (Nalbach, 1994). Previous work has characterized the haltere frequency response patterns in both equilibrium (Sherman and Dickinson, 2003) and aggressive ight motions known as saccades (Bender and Dickinson, 2006a) and shown that they can be modeled as bandpass lters with their outputs summed with visual feedback (Sherman and Dickinson, 2004; Bender and Dickinson, 2006b). In the simulation, simple bandpass lters on the angular rates are used in a feedback loop to model the haltere angular rate feedback 3.3.2 System identi cation method We would like to characterize the rigid body motions about hover in a model amenable to control analysis and design tools. In this section, we present a frequency- based linear system identi cation that was conducted on the simulation using Com- prehensive Identi cation from FrEquency Responses (CIFER) (developed by NASA- Ames (Tischler, 1992)). Frequency-based system identi cation identi es a system model only over a range of frequencies so care must be taken to select the frequencies that are relevant for the desired application. The desired application is understand- ing sensing and feedback requirements for ight control design, so capturing the 50 B ar e A irf ra m e In se ct Si m u la tio n H al te re Fe ed ba ck Sw ee p G en er at o r ? b ? R ec o rd Ti m e H ist o ry A B H ig h fre qu en cy dy n am ic s lo o p ? ? ? ?x i R s i v b ? b ? ? ? ? Figure 3.7: This study considered the haltere-stabilized insect mapping inputs at point A to outputs at point B. A low gain controller was included that modi ed the frequency sweep inputs to keep the system near hover. basic rigid body motions (i.e., up to 20Hz) is desired. The haltere-on airframe motion was the primary system to be identi ed, or the high frequency system in Fig. 3.7 (shown in a dashed box) with inputs taken at point A and outputs taken at point B. The identi cation procedure for the haltere-on system via simulations conducted with the controller and mechanosensory feedback running was designed by considering the haltere feedback as a stability augmentation system (SAS). By considering the halteres as a SAS system, procedures for bare-airframe identi ca- tion from SAS-on ight tests (Tischler and Cau man, 1999) could be used. In most simulated ight test cases, open loop design of the frequency sweeps was su cient to keep the insect near hover, but a low-gain external feedback loop was used pri- marily for haltere-o (bare airframe) identi cations, where the controller-modi ed frequency sweeps were considered as the input. With a ap frequency of approximately 200 Hz (1256 rad/s), Drosophila have dynamic motions at such frequencies. Obviously, wing motion predominates at this frequency, but in a simulation with little or no viscous damping e ects on the 51 body, small body motions can be observed, relating to the force dependence on this ap frequency. In particular, the ap frequency and potentially higher order harmonics of it (2/ ap, 3/ ap, and so on) appear in the state history output. The movement of the center of pressure fore and aft generates a pitching moment that leads to oscillatory behavior. The addition of rotational damping could partially mitigate this e ect, but the smaller oscillatory motions observed at these frequencies elicit haltere responses since the vision-based control system response trails o at high angular rates. Signi cant statistical correlations from input to output was not observed above approximately 30 Hz (188rad/s) and frequency responses above 90 Hz were discarded. Limiting the frequencies identi ed has the advantage of removing the high frequency noise and numerical artifacts during the identi cation process. Low frequency identi cation is limited by the length of the time history pre- sented. In theory, the lowest frequency that could be identi ed in a time history of length Tmax is fmin = 1=(2Tmax). In practice, accurate identi cation requires several cycles of the longest mode. All expected or observed frequencies were faster than the lowest possible frequency in the 14 second (2800 wingbeat) trials used in this study: fmin = 0:04Hz (0.22 rad/s). Instead, based on observed modes and numerical capability, the search for a model was restricted to above 0.1 Hz (0.6 rad/s). 3.3.3 Model structure In consideration of linearized longitudinal dynamics about hover, a traditional ight dynamics model takes the form _x = Ax +Bu (Franklin, 2002), where A and 52 B are most generally de ned as in eqn (3.19). Applied to a helicopter or xed-wing airplane, several of the terms in the A matrix are zero. For example, the term Mw may usually be discarded in hover. In forward ight, linearization of the disc angle of attack shows that Mw is related to the body angle of attack via the forward velocity u0 and may not be discarded. 3.4 Discussion 3.4.1 Non-Parametric Handling Qualities Identi cation One way to characterize the control requirements of the haltere-on airframe dynamics is from the perspective of a human tasked with piloting the vehicle, since determing the computational capabilities of the insect not yet established. Several metrics pertinent to human-piloted rotorcraft small-perturbation handling qualities analysis have been quanti ed and evaluated from calculated frequency responses shown in Fig 3.8-3.11 without the need for a state-space model. Table 3.3 shows the properties important for control and handling qualities evaluation, from which it is clear that the pitch and heave dynamics have higher bandwidths than the fore/aft dynamics. As speci ed in ADS-33 (US Handling Qualities Requirements for Rotorcraft), a gain margin of 6dB and phase margin of 45 was used to calculate the bandwidths (Aviation Engineering Directorate, 2000). Care must be exercised in interpreting the gain and phase margins as the calculation of several of these items involves frequencies beyond the range of good coherence. ADS 33 speci es a number of maneuvers known as mission task elements 53 10 30 50 Magnitude, dB Measured Identified ?270 ?180 ?90 phase, de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (a) Stroke plane angle to q 20 35 50 Magnitude, d B Measured Identified ?450 ?360 ?270 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (b) Stroke o set o to q Figure 3.8: Identi ed transfer functions to pitch rate q. (MTEs) for which the helicopter is evaluated. With respect to pitch control, the relatively high bandwidth (compared with full scale rotorcraft) of the system means ADS-33 meets Level 1 handling qualities requirements for the hover and low speed small amplitude pitch-related MTEs, including target acquisition and tracking. However, were the insect airframe to be humanly piloted, the indeterminate gain bandwidth values are a cause for concern because they indicate a susceptibility to pilot-induced-oscillation. With respect to forward motion, recall that was recom- mended as fore/aft control. In general, forward speed regulation is shown to be the slower bandwidth output, which is likely related to the fact that fore/aft motion is primarily accomplished by tilting the insect?s thrust vector, followed by its body. Such a behavior is analogous to a helicopter tilting its rotor disc, followed by its fuselage, and is slower than other dynamics modes. Experimental evidence indicates that an insect?s neural structure responds with a small time delay, about 3-4ms (Autrum, 1958), and a detailed model of wing 54 ?20 0 15 30 Magnitude, dB Measured Identified ?180 ?90 0 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (a) Stroke plane angle to ?5 15 30 Magnitude, dB Measured Identified ?180 ?90 0 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (b) Stroke o set o to Figure 3.9: Identi ed transfer functions to pitch angle . structural dynamics could add an additional delay. While the low phase margin in the primary to u transfer function suggests that introducing a time delay on the order of 5-10ms could a ect the stability of the system, the coherence in this frequency range is degraded and further investigation is necessary. Function GM PM Gain BW Phase BW Delay to 1 37 1 6.2Hz 0 o to 1 23 1 4.5Hz 0 to w 1 106 1 1 0 to u 5dB 4 3.5Hz 2.2Hz n/a o to u 1 36 2.2Hz 2.2Hz 18 ms Table 3.3: Control and handling qualities properties used to evaluate the system in reference to rotorcraft handling qualities requirements. GM and PM represent gain and phase margins, resctively. 3.4.2 Heave Dynamics By observation that axial ight testing is decoupled from the other inputs, the system was modeled as a 1D system representing the heave dynamics and a system 55 ?15 0 15 Magnitude, dB Measured Identified ?180 ?90 0 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (a) Stroke plane angle to u ?15 0 15 30 Magnitude, dB Measured Identified ?180 ?90 0 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (b) Stroke o set o to u Figure 3.10: Identi ed transfer functions to forward velocity u. in R3 and the heave dynamics identi ed rst in standalone axial simulations. For an insect in hover at some pitch angle 0, with an equilibrium stroke plane angle of 0, vertical force equilibrium provided a means to solve for vertical trim and a corresponding heave model of the form _w = Zww + Zf (f f0) + Z ( 0) + Z ( 0) + Z o ( o o ,0): (3.24) Analytic linearization of the quasi-steady form yielded an expression for the terms in the B matrix, B = 2g cos 0 f0 2g cos 0 0 g tan 0 cos 0 0 : (3.25) 56 ?25 ?5 10 Magnitude, dB Measured Identified ?90 ?45 0 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s Figure 3.11: Identi ed transfer function from stroke amplitude to heave velocity w. For the example Drosophila-like insect, the numerical evaluation of the B matrix using the trim parameters shown in Table (3.4) yields B = 0:098 14:584 0 0 : (3.26) control input value f0 200 Hz 0 77:08 0 0 90 ra 1:6mm Table 3.4: Nominal input parameters used in evaluating B for the hover lineariza- tion. Frequency-based system identi cation of the heave dynamics under excitation via modulating the apping amplitude by 1 degree about the nominal (trim) value yields _w = [ 4:862]w + [ 17:1941] : (3.27) 57 Visible in eqn (3.27) is an aerodynamic damping term Zw which is of aerodynamic origin and was predicted due to the modi ed angle of attack due to the perturbation velocities. (See Fig (3.5) for a physical explanation of this e ect.) A simulation with- out consideration of perturbation velocities does not exhibit positive Zw damping, indicating the e ect is due to the modi cation. 3.4.2.1 E ect of Flap Amplitude The identi ed system in eqn (3.27) includes a control term re ecting the fact that heave control in hover is achieved through direct modulation of the aerody- namic thrust. Experimental evidence indicates that insects do have the capability to measure angular rates and they respond to a mechanical disturbance in pitch rate with a change in stroke plane amplitude and stroke o set o that is nearly linear (Thompson et al., 2008). It is theorized that at high angular rates the am- plitude response falls o , but physical limitations mean measurements have not yet been conducted involving a disturbance faster than approximately 800 /s to verify this conclusion. For this reason, the nonlinear simulation models a high frequency cuto on hatere rate feedback. 3.4.2.2 E ect of Flap Frequency In contrast to ap amplitude, an insect?s ap frequency response to a pitch rate is not a linear relationship and there is dispute over how insects use ap frequency as an input. In the linear analysis presented earlier, ap frequency enters in the 58 same functional manner as ap amplitude, but insects often choose to modulate wingbeat frequency independently of the wingbeat amplitude, commonly varying it in conjunction with other control inputs (Dudley, 1995) (Altshuler et al., 2005) (Vance et al., 2005) . Experimental work has indicated that the dependence of lift on is linear and nearly quadratic on f , but insects decrease ap frequency when approaching peak force outputs. The response exhibits a maximum near 200 /s, and decreases as the angular rate increases from that point (Sherman and Dickinson, 2003), which is presumed to be a physiological limitation (Taylor, 2001). The frequency ranges appropriate for a linear system model are determined by reference to the coherence [0; 1] of the input/output pair, which quanti es the degree to which an input-output pair are related by linear system dynamics. Coherence is de ned via (!) = jGxy(!)j2 Gxx(!)Gyy(!) ; (3.28) where 1 indicates entirely linear system dynamics and 0 indicates an input/output pair with no linear correlation. For the stroke plane amplitude to heave relation, the frequency band of linearity covers much of the range of frequencies expected in an insect control problem. Conversely, the frequency range of acceptable coherence in the ap frequency to heave velocity relation is limited to the higher frequency regions shown in Fig. 3.13. The reduction in linear e ectiveness for stroke plane amplitude modulation (indicated by the negative slope in the Bode magnitude plot) and the complementary increase in ap frequency modulation linear e ectiveness (indicated by the coherence increase from 0 to near 1 and the positive slope in the 59 Bode magnitude plot) may explain why insects tend to use frequency and amplitude in a di erent manner. A complementary set of control inputs mirrors a paradigm observed in insect vision sensing, where visual and mechanical (haltere) feedback mechanisms cover di ering portions of the frequency spectrum. Complementary actuation dynamics are often observed in nature, but this means that a transfer function t that includes both control inputs is a compromise between accuracy of the two. ?25 ?5 10 Magnitude, dB ?90 ?45 0 phase,de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s Figure 3.12: As noted in experimental observations of insects, stroke plane amplitude shows a linear correlation to perturbations. The input controls were chosen based on biological observations, but the lim- ited ap frequency band of acceptable coherence restricts its linear application to relatively fast inputs. Any vehicle, including biological structures, must also be de- signed for optimal operation at or near a particular frequency, limiting the ability to use ap frequency as a control term. For these reasons, ap frequency is not recommended as a primary control term but is included in this study to illustrate its complementary control e ect. Under the restriction of high frequency excitation, 60 ?30 ?15 0 Magnitude, dB ?180 ?145 phase,de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s Figure 3.13: The use of ap frequency f as an input is correlated linearly over a smaller region than . the system can be re-identi ed including the e ect of the control input f . During this identi cation, the Z term was xed to the analytic result and the system was determined to be _w = [ 4:7888]w + 114:4228 14:584 2 6 6 4 f 3 7 7 5 : (3.29) Such a system still shows acceptable t from stroke amplitude to heave velocity, as seen in Fig. 3.14. A time domain comparison of the nonlinear and identi ed system appears later in Figure 3.15(b). 3.4.3 Longitudinal Dynamics While the axial climb (heave) dynamics may be adequately described in iso- lation from the other state variables, an accurate description of the insect?s motion along its s^x axis and pitch motion requires consideration of the coupled system. In 61 ?25 ?5 10 Magnitude, dB Measured Identified ?90 ?45 0 phase,de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s Figure 3.14: The identi ed system ts the stroke plane amplitude to heave velocity relation. this section, system identi cation tools are used to reduce the nonlinear Drosophila- like dynamics in , u, and q to a linearized system (about hover) that may be directly applied to control design. For the system identi cation, simulated ight tests were conducted including frequency sweeps applied to each of the channels in turn with o -axis controller portions turned o (or the entire controller o ) to prevent cross-correlated inputs that can interfere with identi cation. Frequency as an input was discarded for the reasons discussed earlier. The identi ed longitudinal system dynamics matrix for the haltere-on system 62 is identi ed to be Aid = 2 6 6 6 6 6 6 6 6 6 6 6 6 4 12:33 0 0 9:810 0:0 4:651 0 0 547:0 0:0 33:26 0 0 0 1 0 3 7 7 7 7 7 7 7 7 7 7 7 7 5 : (3.30) The system dynamics matrix Aid has the eigenvalues 1 = 38:56, 2;3 = 3:51 11:26i, 4 = 4:65, and the eigenvectors seen in eqn (3.31), con rming the observa- tion that the heave dynamics are uncoupled from the pitch/translational dynamics. 1 = 2 6 6 6 6 6 6 6 6 6 6 4 :0097 0 :9996 :0259 3 7 7 7 7 7 7 7 7 7 7 5 , 2;3 = 2 6 6 6 6 6 6 6 6 6 6 4 0:0578 \ 0:3619 0 0:9948 0:0843 \ 1:8730 3 7 7 7 7 7 7 7 7 7 7 5 , 4 = 2 6 6 6 6 6 6 6 6 6 6 4 0 1 0 0 3 7 7 7 7 7 7 7 7 7 7 5 (3.31) The control input matrix is identi ed to be Bid = 2 6 6 6 6 6 6 6 6 6 6 6 6 4 0 20:67 12:69 17:36 0 0 0 2826:0 6020:0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 7 7 5 : (3.32) Each of Aid?s diagonal elements (Xu, Zw, and Mq) are negative, indicating 63 viscous damping along/about that axis. Thus, the experimentally-derived aerody- namics model applied to rigid body dynamics predicts viscous damping in the ab- sence of body drag forces, i.e., as a fundamental part of apping wing ight. Such a passive stabilization mechanism is highly advantageous for reducing the computa- tional and energetic workload, and parallels the recent work characterizing \ apping counter-torque" (Hedrick et al., 2009), or damping about the s^z axis. 3.4.3.1 Pitch Dynamics The most relevant stability derivatives for pitch dynamics are M , M o , and Mq. A single axis pitch dynamics model involving these terms can help illustrate the role of the pitch sti ness and control terms. For the linear system, apply the principle of superposition to write the pitch rate output as the sum of two SISO systems with step inputs applied to and o such that q(t) = (1 eMq t) M Mq + M o Mq o : (3.33) Much of helicopter handling qualities analysis and control design is focused on increasing the magnitude of Mq, because Mq describes the rate at which a pitch rate reaches a steady-state value for a given step input, or the \quickness" of the vehicle?s pitch response to a pilot input. A major advantage of a modern \hingeless" rotor design is the resulting larger bare airframe jMqj, from about 0.5 for a UH-60 to 5-6 for an aerobatic Bo-105 (He ey et al., 1979b). The addition of a ight control system increases jMqj and is desired for a UH-60 to maintain a reasonable 64 pilot workload. The estimation of Mq = 33:3 indicates that the insect pitch response is largely kinematic rather than dynamic. Such a nding implies that mechanosensory feedback mechanisms can markedly reduce the neural demands of insect ight control, which may help explain how insects are able to achieve high levels of maneuverability on a very small computational budget. The control terms M and M o a ect the magnitude of the steady-state pitch rate reached for a given input. A relative comparison of the two derivatives shows that a perturbation in the o input is more than twice as e ective at creating a pitch rate. This fact, coupled with the fact that the input also leads to a forward speed response twice the size of a o input, veri es the theoretical conclusion that o is the appropriate primary control term for pitch dynamics. Direct comparison of analytic versus identi ed M o is complicated high sensitivity to parameter vari- ations in the estimated term ra, but M has the analytical estimate (using the ra approximation) 4660, as compared to the identi ed value 2826: The pitch response to forward speed de ned by Mu is an important result, as it leads to a nose-up motion in the presence of a forward velocity, and hence a restoring force to resist the increase in speed. Mu > 0 is a criteria for static stability. In a helicopter, a positive Mu is created by rotor blade excursions from the disc plane, but the positive Mu result indicated here was reached without providing the wings the ability to move out of the stroke plane de ned by the input. 65 3.4.3.2 Forward and Heave Velocity Dynamics The forward velocity dynamics are represented by the identi ed Xu, X o , X and an additional inertial term formed by inclination of the gravity vector that remains g throughout this analysis. Xu indicates that apping wing kinematics provide static speed stability even in the absence of body drag forces. The X control e ectiveness term was estimated analytically as g but identi ed to be 20.7. The X and X o terms were the derivatives identi ed with the most uncertainty, with insensitivities of 3 and 5.5%, respectively, indicating that additional dynamic e ects not considered in this analysis may be relevant for the description of this degree of freedom. 3.4.4 Identi ed Longitudinal Dynamics Veri cation To verify the model structure used, additional candidate dependencies not included in the analytic model were included, such as Zu and Mw and delays in the input channels to allow for unmodeled states. In all cases, the added terms were identi ed to zero, small, or had only minor e ects on the transfer function and simulated time histories. The responses of the identi ed and nonlinear systems to a variety of inputs were compared. The identi ed system shows tracking during parallel simulations, and a test involving simultaneous actuation in all input controls is shown in Fig. 3.15. 66 0 1 2 3 4 5 6 7 8 ?0.03 ?0.02 ?0.01 0 0.01 0.02 0.03 B od y f w d v el oc ity u , m /s time, seconds Nonlinear Identified (a) Fore/aft velocity 0 1 2 3 4 5 6 7 8 ?0.06 ?0.04 ?0.02 0 0.02 0.04 0.06 B od y h ea ve v el oc ity w , m /s time, seconds Nonlinear Identified (b) Heave velocity 0 1 2 3 4 5 6 7 8 ?40 ?30 ?20 ?10 0 10 20 30 40 B od y p itc h r at e q , d eg /s time, seconds Nonlinear Identified (c) Pitch rate 0 1 2 3 4 5 6 7 8 ?3 ?2 ?1 0 1 2 3 B od y p itc h a ng le , d eg time, seconds Nonlinear Identified ? (d) Pitch angle Figure 3.15: Time domain comparison of the linear and nonlinear systems. 3.4.5 Accuracy of Identi ed Model Besides evaluating the accuracy of the identi ed model via how well the iden- ti ed and experimental transfer functions agree (see Figs 3.8-3.11) and through comparison of simulation results in section 3.4.4, the accuracy of the each of the parameters in the identi ed model may be estimated as part of the system identi- cation procedure. One method for comparing the quality of identi cation of each parameter is via investigation of how variation in that parameter a ects the global 67 Parameter Value Standard Dev. Insensitivity Xu -12.3 4.0% 1.1% Zw -4.7 9.1% 3.8% Mu 547 3.6% 0.9% Mq -33.3 4.6% 1.1% X o 12.7 12.2% 5.9% X 20.7 7.6% 3.0% Z -17.4 4.8% 2.0% M o 6028 3.6% 1.2% M -2826 3.8% 1.1% Table 3.5: Uncertainty estimates for identi ed dynamics model. cost function describing the di erence of the transfer function t and calculated frequency response, a relative estimate known as the \insensitivity." Values of in- sensitivity less that 10% are considered reasonable, and the identi ed values ranged from less than 1% for Mu to 5.9% for X o . A good estimate of the standard de- viation in each parameter to be expected in repeated ight tests and identi cation sequences may be determined through Cramer-Rao (CR) bounds. Table 3.5 shows uncertainty estimates for the estimated parameters, where the standard deviation upper bounds are derived from the lower bounds provided by CR bounds (Tischler and Cau man, 1999). The uncertainty does not destroy the qualitatitive behavior of the eigenvalue map, shown in Fig. 3.16, which shows the presence of a fast and subsidence mode and an oscillatory pair. The slow subsidence mode associated with the smaller real root is the heave damping pole, while the other three re ect coupled motion of the pitch and fore/aft dynamics. 68 ?45 ?40 ?35 ?30 ?15 ?10 ?5 0 5 10 15 Real Imaginar y ?6 ?4 ?2 0 ?15 ?10 ?5 0 5 10 15 Real Nominal poles Uncertainty spread Figure 3.16: Map of haltere-on system poles with uncertainty perturbations as listed in Table (3.5) preserves the qualitative behavior. 3.4.6 Bare Airframe System Identi cation An identi cation of the haltere-on system revealed a decoupled, stable system that indicated a control system would not need to provide stabilization feedback. Identi cation of the bare-airframe system may also be conducted, which reveals a decoupled heave pole, a fast subsidence mode, and an unstable oscillatory pair, as seen in Fig. 3.17. The nding that haltere-based rate damping is su cient to sta- bilize the unstable oscillatory pair associated with pitch dynamics has experimental support (Fraenkel, 1939; Casas and Simpson, 2007). Since the oscillations are at approximately 2 Hz, the stability of the pair is important from the perspective of quantifying the control requirements{for example, an oscillatory instability at 2 Hz would be unacceptable for a human-piloted aircraft. The bare airframe pole struc- ture seen in Fig. 3.17 is the most common pole structure associated with hovering 69 ?32 ?30 ?28 ?26 ?24 ?22 ?15 ?10 ?5 0 5 10 15 Real Imaginar y ?8 ?6 ?4 ?2 0 2 4 ?15 ?10 ?5 0 5 10 15 Real Nominal poles Uncertainty spread Figure 3.17: Map of haltere-o system (the bare-airframe) poles including un- certainty perturbations indicates the traditional VTOL modal structure (Sun and Xiong, 2005). vertical takeo and landing (VTOL) aircraft, including a hovering Harrier VTOL aircraft (Franklin, 2002), many helicopters which exhibit a lightly unstable Phugoid mode in hover (He ey et al., 1979a,b), and the CFD-based analysis (Sun and Xiong, 2005) of a bumblebee. That the same pole structure results from complex computa- tional uid dynamics analysis and from a relatively simple application of curve- tted aerodynamics to a rigid body simulation is signi cant, but the observation of the same pole structure across di ering scales and vehicle con gurations (rotorcraft, xed wing VTOL, and dipteran apping wing) is even more remarkable, because dynamic similarity of the two systems is not suggested by traditional scaling argu- ments such as Froude number comparison (Wolowicz et al., 1979). Instead, linear system analysis has been used in this study to reveal similar passive aerodynamic mechanisms on di erent scales. 70 3.5 Summary In this chapter, the curve- tted insect aerodynamics model (in quasi-steady form) was posed in the context of body translation and rotation by adding the con- cept of perturbation velocities to formulate the wing aerodynamics as functions of the insect (or vehicle) state. The resulting aerodynamic forces and moments were applied to 6 DOF rigid body equations of motion. The equations were rst analyzed via a wingstroke-averaged force and linearization method in analytic form to give controllability term estimates. A numerical simulation of a dipteran insect was then constructed. Geometric, aerodynamic, and inertial properties of a Drosophila-like insect were used in the simulator to create time histories that were then analyzed in the frequency domain, and two linear systems identi ed from the frequency re- sponses, showing that heave dynamics and fore/aft/pitch models are decoupled in the neighborhood of hover. The perturbation velocity analysis indicates that, for the purpose of controlling velocity along the s^z axis, ap frequency and ap amplitude are complementary inputs, but that frequency has only a limited band of linear controllability and is more di cult to use in a linear controls context, while amplitude is much more readily used and is well captured by a linear rst order model. It also indicates that an insect modulating the bio-inspired input parameters suggested may e ect motion in each longitudinal state (excepting pitch angle, controlled via integration of pitch rate). The identi ed system indicates a dipteran insect of these dimensions and wing motions has ample control authority over its longitudinal motion, including 71 pitch dynamics control that signi cantly exceeds modern full scale helicopter results. Though the system is highly maneuverable, haltere-on pitch control is predicted to qualify for full scale Level 1 small perturbation analysis, indicating the insect is presented with a relatively low computational workload. The implication is that a apping wing micro air vehicle design strategy may leverage inherent dynamic properties of the apping wing aerodynamics to reduce size, weight, and power requirements. Frequency based system identi cation of the haltere-on system revealed a sta- bilized system with a slow and fast subsidence mode, and a stable oscillatory mode. Identi cation of the haltere-o (bare airframe) system revealed a slow and fast sub- sidence mode, and an unstable oscillatory mode. In both cases, the slow subsidence mode was associated with the decoupled heave dynamics. While this is the rst look at the haltere-on modal structure, the bare-airframe modal structure is consistent with the previous CFD estimate of a dipteran insect bare-airframe modal struc- ture. The implication that rate damping due to halteres is su cient to stabilize the unstable oscillatory pair signi cantly reduces the computational workload expected for apping wing insects and provides insight into how insects achieve unparalleled maneuverability with limited neural processing. 72 Chapter 4 Lateral-Directional Hover Dynamics 4.1 Introduction Chapter 3 analyzed the longitudinal dynamics of a apping wing insect about hover using curve- tted aerodynamics models (Dickinson, 1996) and Euler rigid body dynamics applied to fruit ies. The results indicated slow and fast subsidence modes, as well as an unstable oscillatory mode that could be stabilized via halteres providing pitch rate feedback. This chapter continues to examine the implication of passive aerodynamic stability mechanisms associated with apping ight, again using hovering fruit ies, but characterizing instead lateral-directional motion. The goal is a simpli ed dynamics model of apping ight lateral motion about hover, which may be readily interpreted via traditional full scale aircraft dynamics and linear control analysis techniques. The underlying mechanics are again posed as nonlinear Euler rigid body dynamics paired with quasi steady aerodynamics model- ing that includes the e ects of perturbations from the equilibrium. The results show that the insect aerodynamics model placed in context of state perturbations is able to predict dynamic behavior without the use of a more detailed and computationally intensive 3D CFD study such as (Ramamurti and Sandberg, 2002, 2007). The organization of the chapter is as follows. Section 4.2 brie y reviews wing kinematics and aerodynamics with asymmetric inputs. Section 4.3 describes the sys- 73 tem identi cation procedure using a nonlinear simulation environment encoding the perturbation velocities concept. Section 4.4 shows how the perturbation velocities concept can be used to develop estimates of the stability and control derivatives, and Section 4.5 presents identi ed models and shows how these results suggest passive aerodynamic mechanisms that may reduce an insect?s ight control requirements about hover. 4.2 Background The insect aerodynamic theory and the governing equations for the analysis and simulation used for lateral dynamics analysis are analogous to Chapter 3 and reviewed only brie y here. 4.2.1 Wing Kinematics As in Chapter 3, the right and left wing motion is quanti ed through a set of 2-3-2 Euler angles. Each wing stroke is planar, with r and l being the inclination of this plane on the right and left sides. The wing position within this plane is notated as r and l, and the inclination (twist) of each wing relative to this plane is the geometric angle of attack g. The right and left wing variables seen in Fig. 4.1 are interpreted as longitudinal and lateral control inputs [ ]c and [ ]d using the transformation in Section 3.1.1. Asymmetric inputs were zero for the symmetric motion analysis in Chapter 3 and the subscripts suppressed, but both collective and di erential inputs are relevant for lateral-directional motion. 74 ?r ?l (a) Right and left ap amplitude r and l. St ro ke Pl an e Ti lt (b) Right and left stroke plane amplitude r and l. Figure 4.1: Control input de nitions. 4.2.2 Aerodynamics The aerodynamics used in this analysis is the insect aerodynamics model de- veloped using experimental ts to Drosophila and Robo y data (Dickinson et al., 1999). While this model proposes lift and drag purely as functions of the instan- taneous wing motion without in ow e ects, the aerodynamic forces are also func- tions of the body motion throughout the surrounding uid. Chapter 3 extended the quasi-steady formulation via a revised de nition of the wing tip speed to be v = v + vG + v!, where v , vG, and v! are the velocity components due to ap- ping, body translation, and body rotation, respectively. The wing angle of attack must also be modi ed as = g + arctan vz vx , where v = vxr^x + vyr^y + vz r^z is the total velocity expressed in the right wing frame. 4.3 System Identi cation To facilitate the identi cation of an equivalent linear system, the simulator developed in Chapter 3 was used in conjunction with Comprehensive Identi cation 75 from FrEquency Responses (CIFER). The goal of this process was to identify a linear system that would allow predictions of the apping wing ight control properties. Frequency sweeps (\chirp" signals) were applied to the inputs d and d to excite the lateral/directional dynamics of the insect. Based on the spectral response of the input Gxx(!) and output Gyy(!) (calculated by the chirp z transform), the transfer function G may be found via G(!) = Gyy(!) Gxx(!) : (4.1) As a linearity measurement, coherence 2 [0; 1] was introduced in Section 3.4.2(de- ned in eqn (3.28)) to quantify the degree to which the magnitude and phase of an input/output pair may be described by linear system dynamics. Our desired result is a model of the form _x = Ax+Bu (Franklin, 2002), where A is most generally taken to be A = 2 6 6 6 6 6 6 6 6 6 6 4 Yv 0 g 0 L0v L 0 p 0 L 0 r 0 1 0 0 N 0v N 0 p 0 N 0 r 3 7 7 7 7 7 7 7 7 7 7 5 (4.2) 76 and B = 2 6 6 6 6 6 6 6 6 6 6 4 Y d 0 L d 0 0 0 0 N d 3 7 7 7 7 7 7 7 7 7 7 5 : (4.3) 4.4 Wingstroke Averaged Forces and Moments In this section, the process of wingstroke averaging developed in Chapter 3 is applied to demonstrate how linearization of rigid body equations of motion forced by experimentally-derived insect aerodynamics can provide estimates of the control input terms and also indicate physical mechanisms providing aerodynamic damp- ing. The linear systems representing rigid body motion were derived with the un- derstanding that we are concerned only with the low frequency rigid body modes describing the ight path of the vehicle, i.e., that the averaged forces and moments (over a wingstroke) were the essential part of the external forcing function. Ap- plication of this principle provides insight into what physical mechanisms provide passive aerodynamic damping. 4.4.1 Control Sensitivity A system identi cation requires an initial estimate of the relevant parameters. The control sensitivity terms, or the entries of the B matrix (eqn (4.3)), were es- timated via the wingstroke averaging procedure. The roll sensitivity to the inputs d, o ,d, and d is derived here; other terms may be found in like manner. 77 4.4.1.1 Roll Moment due to Di erential Stroke Amplitude d For the case of zero stroke plane angle = 0, the roll moment can be written in terms of the left and right wing lifts as Lroll = (Ll Lr)ra, where Ll and Lr are de ned as Lx = 2(r^22 S R 2)(f 2x 2 x)Cl with x = l; r. Though the angle of attack on each wing is di erent in general, they share the same lift coe cient during a perturbation of d while hovering. Using d as de ned in eqn (3.2), the roll moment is Lroll = 2 (r^22SR 2)ra [(fc fd)2( c d)2 cos( o ,c o ,d) (fc + fd)2( c + d)2 cos( o ,c + o ,d)]CL (4.4) Linearizing with respect to d and substituting in the nominal conditions d = fd = o ,d = 0 to obtain @Lroll @ d = 4 2 (r^22SR 2)CL ra cos( off;c)f 2 c c: (4.5) Using the nominal conditions for Drosophila in Table 4.1, @Lroll@ d = 8:84 10 9CL N m/rad. For g = 45, CL = 1:8, and a roll inertia of 5:9 10 13Nms 2, L d = 26; 969 rad/s 2. control input value f0 200 Hz 0 77:08 0 0 90 ra 1:6mm Table 4.1: Nominal input parameters used in evaluating B for the hover lineariza- tion. 78 4.4.1.2 Roll Moment due to Di erential Stroke Bias o ,d Similarly, computation of the derivative with respect to o ,d and substitution of the nominal conditions yields @Lroll @ o ,d = 2 2 (r^22SR 2)raf 2 c 2 c sin( o ,c); (4.6) indicating that for a nominal collective stroke o set of 0, a di erential stroke o set does not lead to a roll moment via changing lift. 4.4.1.3 Roll Moment due to Di erential Stroke Plane Inclination d The body frame Z force due to right wing lift is Zr = Lr cos( ) (and similarly for the left side). Thus, the roll moment may be found after replacing Lr with Lr cos( ) to get Lroll = (Ll cos( l) Lr cos( r))ra. De ning d and c as in eqn (3.2) such that r = c + d and l = c d, one may rewrite the roll moment as Lroll = [Ll cos( c d) Lr cos( c + d)]ra: (4.7) By recognition that Ll and Lr are invariant under perturbations, eqn (4.7) may be di erentiated with respect to d directly to get @Lroll @ d = raLr sin( d + c) + raLl sin( d c); (4.8) 79 and evaluated at c = 0, @Lroll @ d c=0 = ra sin( d)(Lr Ll): (4.9) Equation (4.9) is zero at the reference d = 0, revealing that di erential stroke plane inclination about c = d = 0 does not create a roll moment. 4.4.2 Physical Basis for Passive Aerodynamic Mechanisms Wingstroke averaging applied to rigid body equations of motion and the insect aerodynamics model including perturbations indicates several physical mechanisms that provide aerodynamic damping. Sideslip damping is used as an example, and the other terms can be found via a similar process. 4.4.2.1 Sideslip Damping Consideration of imposed velocities has led to the prediction of a ap damping term (Hedrick et al., 2009). In this section, an imposed sideslip velocity is expressed via the velocity and angle of attack perturbations to predict sideslip damping, or static stability in this degree of freedom. Consider the case of positive sideslip along s^y (a relative wind from the right). The total motion of the wing is now the sum of apping and vehicle translation motions. Figure 4.2 shows that sideslip velocity increases the air ow over some portions of the wingstroke (darker) and decreases it over other portions (lighter), where the inboard region of reverse ow has been neglected. For the nominal case 80 of a level stroke plane with no heave velocity, a sideslip component of velocity does not change the angle of attack. ? g s? y A dv an ci ng s tro ke s? x R el at iv e w in d + ++++ + + + + ++ + +++ ++ + + + + + + + - - - - - - - - - - - - - - - -- - - -- - - - -- -- - - - - - - -- - - - - - -- - -- - - -- -- - - - -- In cr ea se d ve lo ci ty D ec re as ed v el oc ity + + + + - - - - - (a) Advancing s? y R et re at in g st ro ke s? x R el at iv e w in d ? g + +++ + + + + + ++ + + + + ++ + + + + + + + - - - - - - -- - - -- - -- - - - -- - - - -- --- - -- - - -- - - - - --- - - - - - - - -- - -- - - In cr ea se d ve lo ci ty D ec re as ed v el oc ity + + + + - - - - - (b) Retreating Figure 4.2: An imposed sideslip increases and decreases relative air ow over regions of the wingstroke. The translational model of drag accounts for velocity via a quadratic depen- dence on the local velocity. The drag is increased in those regions of increased ow and reduced in those regions of reduced ow as seen in Fig. (4.3), which shows that in both the advancing and retreating strokes, the perturbation in drag leads to a net force in the s^y direction, acting to oppose the motion. Analytically, instantaneous right wing drag force under a sideslip motion ex- 81 ? g s? y A dv an ci ng s tro ke s? x R el at iv e w in d N om in al d ra g Perturbed dra g (a) Advancing ? g s? y R et re at in g st ro ke s? x N om in al d ra g Perturbed dra g R el at iv e w in d (b) Retreating Figure 4.3: An imposed sideslip increases and decreases drag forces during the advancing and retreating strokes in a manner leading to damping in sideslip v about hover. pressed in stability frame, with r = 0 and o ,r = 0, is [FD;r]S = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1=2 S r^22 CD (2 r sin (2 fr t) fr R + sin ( r cos (2 fr t)) v) 2 cos ( r cos (2 fr t)) 1=2 S r^22 CD (2 r sin (2 fr t) fr R + sin ( r cos (2 fr t)) v) 2 sin ( r cos (2 fr t)) 0 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 S (4.10) for the advancing stroke and the negative of eqn (4.10) for the retreating stroke. The left wing has a similar expression. The wingstroke averaged side force is d dv [Yavg]S = d dv Z T 0 [FD;r]S + [FD;l]S dt: (4.11) Interchanging the order of integration and di erentiation, the right and left wing terms contribute equally. De ne sideslip damping Yv = 1m @Yavg @v to write 82 Yv = 4 S r^22 RCD m [cos ( c) sin ( c) c] (4.12) The values given by eqn (4.12) directly compare against other air vehicles as in Fig. (4.4). 0 20 40 60 80 100?1.2 ?1 ?0.8 ?0.6 ?0.4 ?0.2 0x 10 ?3 Amplitude , degrees Y v as f( ) Insect B 747 UH?1H BO?105 XC?142 Si de sli p da m pi ng Y v =F r; s ? 1 ? c ? c Figure 4.4: Comparison of calculated Yv with Boeing 747, UH-1 Huey, BO-105 helicopter, and XC-142 tactical transport. All values normalized by Froude number. 4.4.3 Active Feedback Mechanisms The ight stabilization systems of many dipteran insects include halteres, aerodynamically-ine ective hind wings that are used as angular feedback sensors (Nalbach, 1994). The kinematic response d due to mechanosensory feedback has been characterized for tethered Drosophila (Sherman and Dickinson, 2003). The measurements were used to model the mechanosensory feedback as a bandpass lter on linear roll rate feedback to d, in a similar manner to the pitch rate feedback used in Chapter 3, The yaw rate feedback to d output has not been characterized exper- 83 imentally and was not included in the simulated insect. In general, rate feedback acting on a system may be written as u = Kx, where the gain matrix K has only terms dependent on angular rates. While rate feedback alone will not stabilize an unstable system, rate feedback has the e ect of additional damping in the system. Application of yaw rate feedback would further reduce the bandwidth demands for the ight stabilization and control task via the increase in the yaw damping term N 0r due to the additional damping. 4.5 Results and Discussion 4.5.1 Yaw Dynamics In this section, the yaw dynamics results are compared to previous modeling work. Consider the yaw dynamics case previously discussed in Hedrick et al. (2009), where a linear yaw dynamics model of the form _r = a r was calculated, where a is formed from the wing?s geometric and kinematic param- eters as a = R4 c r^33 2 f CF sin( )(d ^=dt^) I 0zz : Here, the term formed by the force coe cient CF , geometric wing angle g, and nondimensional ap speed d ^=dt^ has the mean value CF sin( g)(d ^=dt^) 6:0. Using parameters from the simulated insect, the FCT model predicts an aerodynamic 84 damping coe cient of N 0r = 167 /s, somewhat higher than both the analytic and the system identi cation results. When the simulated insect was tethered with only a yaw degree of freedom, an equivalent linear system was identi ed as _r = 73:0 r 28; 100 d: (4.13) Despite the di erences in the stability derivative, the associated transfer functions from d to r as calculated from the simulation, FCT (Hedrick et al., 2009), and identi ed models may be seen to have very good spectral agreement in Fig. 4.5. 40 45 55 Magnitude, d B Measured Identified FCT Model ?270 ?235 ?180 phase, de g Measured Identified FCT Model 100 101 102 0 0.5 1 Coherenc e Freq, rad/s Figure 4.5: Yaw input d is well correlated linearly to yaw rate r and the 3 available models show excellent spectral agreement. 85 4.5.2 Roll Dynamics 4.5.2.1 Linear System Application of d sweeps con rms that d is the roll input control as seen in Fig. 4.6,1 with no linear component to the yaw output (as indicated by the low coherence). The degradation in coherence seen in the other plots at approximately 4-5 Hz is related to a resonant phenomenon in the roll axis that yielded motion outside a linear region. This motion appears in the identi ed system as a lightly damped oscillatory pair. A linear roll dynamics system may be identi ed as 2 6 6 6 6 6 6 4 _v _p _ 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 4 9:86 0 9:81 9830 177 0 0 1 0 3 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 4 v p 3 7 7 7 7 7 7 5 + 2 6 6 6 6 6 6 6 6 4 16:5 0 12200 0 0 0 3 7 7 7 7 7 7 7 7 5 2 6 6 4 d d 3 7 7 5 : (4.14) 4.5.2.2 Physical Mechanism Predicting Roll Damping The quasi steady model with perturbations predicts roll damping via di eren- tial angle of attack consideration. The angle of attack modi cation is similar to the heave damping case discussed in the longitudinal analysis, and is seen in Fig. 4.7. 1Identi ed system shown in transfer function plots is the nal full lateral/directional coupled system in section 4.5.3. 86 0 15 30 45 Magnitude, d B Measured Identified ?270 ?180 ?90 phase, de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (a) d to p ?15 0 15 30 Magnitude, d B Measured Identified ?360 ?270 ?180 phase, de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (b) d to ?15 0 15 Magnitude, d B Measured Identified ?450 ?360 ?270 ?180 phase, de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (c) d to v 0 15 30 60 Magnitude, d B Measured Identified ?270 ?180 ?90 phase, de g Measured Identified 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (d) d to r Figure 4.6: Roll input d drives the roll/sideslip dynamics, but not yaw. 4.5.3 Full Lateral Dynamics The two identi ed systems in eqn (4.13) and (4.14) describe the behavior of the isolated systems, but the general system included in eqn (4.2) and (4.3) includes numerous cross coupling terms that may be important. The inclusion of the cross 87 N et m om en t a ct s to op po se th e m ot io n R ol l m ot io n Li ft p er tu rb at io n du e to lo ca l ch an ge in a ng le o f a tt ac k A B C Figure 4.7: The applied roll motion in (a) creates a counter-moment in (c) via the di erence in angle of attacks imposed by the velocity distribution in (b). coupled terms does change the state space representation of the system to be 2 6 6 6 6 6 6 6 6 6 6 4 _v _p _ _r 3 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 6 6 4 9:69 0 9:81 0 9720:0 177:0 0 2:07 0 1 0 0 167:0 462:0 0 71:0 3 7 7 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4 v p r 3 7 7 7 7 7 7 7 7 7 7 5 + 2 6 6 6 6 6 6 6 6 6 6 6 6 4 8:24 0:0 12300:0 0 0 0 0 28100:0 3 7 7 7 7 7 7 7 7 7 7 7 7 5 2 6 6 4 d d 3 7 7 5 : (4.15) During this identi cation, the earlier terms were allowed to vary freely as well, which results in slight changes in the numerical values. As an example, N 0r shifts from -73 in the yaw-only model to -71 in the full lateral model. One may com- pare these changes to the error bounds calculated later in Table 4.2 to determine if the changes are signi cant. Despite the change in state space appearance, model structure determination applied to the system indicates that the additional cross terms are poorly identi ed and should not be included in this model. The system pole structure (which determines the qualitative nature of the body motion) is also remarkably consistent between the decoupled (eqn (4.13) and (4.14)) and coupled (eqn (4.15)) system, as indicated by Fig. 4.8, which represents the rst formal inves- 88 ?150 ?100 ?50 0 ?20 ?10 0 10 20 Imaginar y Real Decoupled Poles Coupled Poles Figure 4.8: Comparison of the poles identi ed for the coupled and uncoupled system. tigation of the lateral dynamics poles of a Drosophila-like insect. Both cases show a fast and slow subsidence mode, and an oscillatory mode close to instability, which suggests signi cantly reduced insect neural dependence for ight control. The pole structure observed in Fig. 4.8 is also consistent with other manned ight vehicles, such as the UH1 Huey helicopter (Franklin, 2002). The fact that qualitatively sim- ilar stability behavior (pole structure) is observed between two vehicles where both the scale and the lift production mechanisms are quite di erent is remarkable: there is no dynamic reason to expect such similarity, and the term that would normally be matched for dynamic similitude (Froude number) is not matched in this case (Wolowicz et al., 1979). While the eigenvalues (system poles) of a linear system determine the qual- itative nature of the body response, the eigenvectors indicate the direction of this response. For the union of the roll/sideslip and yaw system (eqn (4.13) and (4.14)), 89 Parameter Value Standard Dev. Insensitivity Yv -9.69 5.1% 1.8% L0v -9720 6.3% 1.2% L0p -177 7.1% 1.3% N 0r -73 9.0% 2.1% L d -12,300 5.9% 1.2% N d -28,100 7.8% 1.8% Table 4.2: Uncertainty estimates for parameters in the identi ed dynamics model. the poles ( 1; 2,3; 4) = ( 73; 3:4 22:9i; 180:5) have the eigenvectors 1 = 2 6 6 6 6 6 6 6 6 6 6 4 0 0 0 1 3 7 7 7 7 7 7 7 7 7 7 5 2,3 = 2 6 6 6 6 6 6 6 6 6 6 4 0:018; 173 0:999 0:043; 98 0 3 7 7 7 7 7 7 7 7 7 7 5 4 = 2 6 6 6 6 6 6 4 0:000 1:000 0:006 3 7 7 7 7 7 7 5 ; (4.16) indicating roll and yaw subsidence modes 1 and 4 along with an oscillatory roll- dominated roll/sideslip mode 2; 3. An estimate of the uncertainty in the identi ed stability derivatives may be derived from the Cramer-Rao bounds using the methods in Tischler and Remple (2006). The resulting uncertainties, shown in Table 4.2, do not quantitatively change the identi ed system pole structure, as indicated in Fig. 4.9, which shows the nominal poles and the poles perturbed by normally-distributed noise corresponding to the uncertainty predictions. 90 ?200 ?150 ?100 ?50 0?25 ?20 ?15 ?10 ?5 0 5 10 15 20 25 Real Imaginar y Nominal poles Uncertainty spread Figure 4.9: The identi ed system pole structure is preserved under uncertainty perturbations. 4.6 Summary In this chapter, the insect aerodynamics model including perturbations was analyzed to determine reduced order control models that may be used to estimate the control properties associated with lateral-directional ight. The models were derived with the goal of describing the low order rigid body motion of an insect, (not high frequency structural dynamics), so frequency-based system identi cation tools and wingstroke-based averaging were both used in deriving them. In each case, a curve tted aerodynamics was used in the context of untethered dynamic motion of the insect body, rather than a detailed 3D ow solution found during a prescribed body motion, such as three-dimensional ow elds computed in (Ramamurti and Sandberg, 2007) or in the Model Validation of Chapter 5. For the example insect considered (Drosophila-like insect), the analysis indi- cates two stable real poles, and two very lightly-damped and nearly unstable complex 91 poles. One of the real poles is found to correspond to the yaw dynamics of the insect, showing a decoupling of motion about this axis near hover, allowing comparison to Hedrick et al. (2009). The other three poles describe a coupled roll/sideslip motion that may be excited by a roll input, de ned in this study as a di erential modulation of the wingstroke amplitude. The yaw dynamics input was found to be a di erential stroke plane angle (inclination). The models were also used to estimate the e ect of pole additional o -axis dynamic cross-coupling and uncertainties, which show that the dynamic structure is preserved in these cases. Haltere-based roll rate feedback was integral to restraining the modeled insect near equilibrium and was included in these models, which allows the models to be driven directly via the control compo- nents that have signi cantly higher latencies (and thus lower bandwidth), such as the visual feedback system Humbert and Hyslop (2010). The dynamic analysis suggests that mechanosensory rate feedback is an inte- gral part of the animal?s control strategy. The study?s analysis also suggests that inherent passive aerodynamic mechanisms due to di erences in angle of attack and dynamic pressure, can act to stabilize the vehicle and thus reduce the ight stabi- lization/control requirements. From the perspective of developing robotic apping wing micro air vehicles, the ability to leverage passive aerodynamic mechanisms to assist in ight stabilization is attractive to reduce the size, weight, and power requirements of the control system that must be carried onboard. 92 Chapter 5 Model Validation 5.1 Introduction An important problem in apping micro air vehicle (MAV) design is a method of rejecting wind gusts and maneuvering through unknown wind gradients in a rapidly changing environment. Fundamental to solving this problem is the level of sensing and feedback requirements inherent in apping wing ight. This dis- sertation has derived reduced order models of apping wing ight using simpli ed aerodynamics, and this chapter compares reduced order modeling results against more traditional modeling, using computational uid dynamics (CFD) modeling in collaboration with Mac MacFarlane and Brandon Bush (Bush, B. et al., 2010). The high delity CFD solutions incorporate more complex aerodynamics such as unsteady ow, wing/wake interaction, and wake capture e ects. The reduced order longitudinal dynamics model for a hovering Drosophila (fruit y) presented in Chapter 3 had bare airframe instabilities but could be stabi- lized with the addition of a biological pitch rate feedback. The model was developed using a comparatively simple aerodynamics model. Euler rigid body dynamics were then used to simulate a hovering insect and a system identi cation performed to model the inputs and outputs. A CFD-derived solution could signi cantly improve the delity of the aerodynamics model, at the expense of dramatically increased 93 complexity. The goal of this chapter is to directly compare stability derivatives as well as pole locations derived using each method. The longitudinal dynamics of the insect are described in the standard state space form _x = Ax +Bu, where x = [ u w q ]T is the vector of longitudinal perturbation states and u is the vector of inputs to the system. The stability matrix, A, is comprised of stability derivatives: linear changes in a stroke-averaged force or moment due to perturbations from hover. A = 2 6 6 6 6 6 6 6 6 6 6 6 6 4 Xu 0 0 g Zu Zw 0 0 Mu Mw Mq 0 0 0 1 0 3 7 7 7 7 7 7 7 7 7 7 7 7 5 (5.1) To calculate the stability derivatives, the stroke-averaged body forces, X and Z, and pitching moment, M , were measured for various perturbations from hover. As before, only the forces and moments averaged over the course of a wingstroke are considered important to the rigid body dynamics, due to the high frequency of the wingstroke with respect to the frequency of the body dynamics.(Deng et al., 2006b; Kevorkian, 1966) The stroke-averaged forces and moments can be calculated using both the curve- tted aerodynamics model of aerodynamics and CFD, providing a stability matrix derived from each method. This allows a direct comparison of the linearized system generated by curve- tted aerodynamics to the linear system produced by fully unsteady aerodynamics, as given by CFD. 94 5.2 Unmodeled Aerodynamics and Velocity Perturbations 5.2.0.1 Rotational Lift and Other Aerodynamic Models There are aerodynamic e ects not adequately represented in the translational model in Section 3.1.3. In particular, force peaks appear at the start and end of the wing strokes where the translational representation shows only small forces due to the small velocities while the wing is reversing direction. Wing incidence rapidly reverses at these instants, so a convenient representation is to write this additional lift peak as a function of the wing incidence derivative _ g. While rotational e ects can account for up to 35% of the lift generated(Dickinson et al., 1999), modeling lift in this manner requires a theoretical or experimental determination of the wing axis of rotation that has not been demonstrated with the precision necessary for accurate predictive capability in the general case of dipteran insect ight (Bennett, 1970). Other force production models involve a variety of unsteady e ects such as added mass e ects (approximately 5% of lift forces). This study addressed the e ects of translational lift, which applies to dipteran insect ight in general. 5.2.0.2 Reverse Flow In contrast to the curve- tted model?s ut, the total ow incident on the wing v is no longer necessarily opposite to the ap direction. An example is in forward ight during the retreating stroke where inboard sections of the wing generate negative lift. A region of reverse ow is a common problem encountered on helicopter rotors where the magnitude of its e ects are quanti ed using a nondimensional speed known as the 95 advance ratio = u=ut = u_ R . For a Drosophila-like insect apping at f = 200Hz, the peak angular rate is max( _ ) = (2 f) 1755rad/s, leading to a peak wing velocity at the tip of max(ut) = 3:7m/s, and hence a minimal advance ratio at the insect?s preferred forward speed of 2m/s (Vogel, 1967), the corresponding advance ratio of 0.54 for a level stroke plane angle is the lowest advance ratio in \cruise." For comparison, a typical helicopter can y no faster than approximately max = 0:35 to 0:40 (Johnson, 1994) and reverse ow e ects are likely a signi cant factor in forward ight. In this analysis, the insect is restricted to motions less than ~10 cm/s or < 0:03 at peak angular velocity. Helicopter analysis routinely discards reverse ow below = 0:1, so reverse ow e ects are expected to be secondary and this e ect was assumed negligible in the reduced order aerodynamic model presented in Chapter 3. 5.2.0.3 Spanwise Flow Where the traditional quasi-steady aerodynamics formulation assumed vy = 0, vy is now a periodic function of . Spanwise ow is a term that is often neglected in the wing frame airfoil lift and drag calculations, even in extensive detailed numerical simulations of traditional helicopters via blade element theory, because published airfoil data is commonly in 2D and one may invoke the \independence principle" (Jones and Cohen, 1957). In performance calculations, an estimate of the rotor?s pro le drag along this axis is usually included (Leishman, 2006). As a consequence of the 2D nature of the airfoil lift and drag representation, spanwise ow was neglected 96 in the reduced order aerodynamic model presented in Chapter 3.. 5.3 Methodology 5.3.1 Kinematics Wing kinematics are de ned about each hinge as in Section 3.1.1, apping in a plane inclined from body axes by . Motion of the wing within the plane is termed stroke angle , and rotation relative to the plane is the geometric angle of attack g. (t) = cos(2 ft) + off (5.2) g(t) = max tanh[2:7 sin(2 ft+ 180 )] (5.3) where and max are stroke and wing pitch amplitudes, f the apping frequency, off a stroke bias term, and is a relative phase di erence between and g. In the nominal stereotyped Drosophila hovering kinematics, and off are both zero, = 74:9 , f = 200 Hz, max = 45 , and = 66 (Dickinson et al., 1999). 5.3.2 Reduced-Order Aerodynamic Model While apping wing aerodynamics are both complex and unsteady, the reduced order aerodynamic model presented in Section 3.1.3 and used to develop the models presented approximates wing forces as functions of the wing tip speed and angle of 97 attack. Lift L is de ned by L = 1 2 U2tipSr^ 2 2CL( ) (5.4) where is air density, S is the surface area of the wing, and r^2 is the non-dimensional second moment of wing area as de ned in (Ellington, 1984c). Drag is de ned simi- larly in Section 3.1.3, and the lift and drag coe cients are measured experimentally by Dickinson and Gotz (1999) for a Drosophila wing. Other terms suggested by Sane and Dickinson (2002), such as Kramer e ect, added mass, or wake capture are not included in the present calculation. Perturbations from hover a ect the tip speed and angle of attack. 5.3.3 Computational Methodology Detailed ow computations around the hovering Drosophila are provided by an immersed boundary incompressible Navier-Stokes solver (IBINS) previously val- idated for apping Drosophila wings via experimental comparison (Bush, B. and Baeder, J., 2008; Bush, B. et al., 2010). The insect is modeled as three bodies moving through a Cartesian mesh with hyperbolic grid spacing, as seen in Fig. 5.1. During the computational solution, similtude requirements require both the Reynolds? numberRe = 120 and reduced frequency f^ = 0:19 to be maintained.(Shyy, W. et al., 2008) Reynolds number consistency is necessary for aerodynamic simil- tude and is determined by the mean wing chord c, average wing tip speed Uref, and 98 Figure 5.1: The model drosophila in the computational domain used for IBINS calculations. viscosity . Re = cUref = 120 (5.5) Reduced frequency (related to Strouhal number) is introduced if the Navier-Stokes equations are nondimensionalized via apping frequency. Reduced frequency de- scribes vortex shedding (unsteadiness) and can be shown to be a geometric ratio involving stroke amplitude and aspect ratio, as in eqn (5.6). k = f c Uref = c 4R = 0:19 (5.6) 5.4 Computational Results State perturbations were applied to the model for the longitudinal states, and stroke-averaged forces and moments were calculated. The trim force along the body frame s^x (fore-aft) axis is similar, while discrepancies exist along the body-frame s^z axis, most likely a function of the previously unmodeled unsteady aerodynamics. 99 For example, wake capture e ects shown in Fig 5.5 have a signi cant e ect on the ow structure surrounding the wing while maneuvering. For each force or moment, the derivatives with respect to states u;w; and q are estimated via a linear regression (shown as a slope) and scaled by mass or inertia to calculate the stability derivatives populating the A matrix. ?1 ?0.5 0 0.5 1?2 ?1.5 ?1 ?0.5 0 0.5 1x 10 ?5 Surge Velocity, m/s Stroke?Averaged Forces, N X u QS = ?5.49 Z u QS = 0.0 X u CFD = ?5.55 Z u CFD = ?0.285 XQS ZQS XCFD ZCFD (a) ?1 ?0.5 0 0.5 1?5 0 5x 10 ?9 Surge Velocity, m/s Stroke?Averaged Moments, N m M u QS =730.1 M u CFD =4554.4 MQS MCFD (b) Figure 5.2: Stroke-averaged forces (a) and moments (b) as given by the curve- tted experimental model (red) and IBINS (blue) for various perturbations in u. Table 5.1 summarizes the stability derivatives for the curve- tted insect aero- dynamics model and the IBINS aerodynamic calculation. Comparison of the deriva- tives all shows the same sign. The magnitudes of the translational derivatives agree 100 ?1 ?0.5 0 0.5 1?2 ?1.5 ?1 ?0.5 0 0.5 1x 10 ?5 Heave Velocity, m/s Stroke?Averaged Forces, N X w QS = 0.00122 Z w QS = ?2.49 X w CFD = ?0.736 Z w CFD = ?3.70 XQS ZQS XCFD ZCFD (a) ?1 ?0.5 0 0.5 1?5 0 5x 10 ?9 Heave Velocity, m/s Stroke?Averaged Moments, N m M w QS =931.4 M w CFD =1730.6 MQS MCFD (b) Figure 5.3: Stroke-averaged forces (a) and moments (b) as given by the curve- tted experimental model (red) and IBINS (blue) for various perturbations in w. well, as seen in Figs 5.2{5.4. Xu in particular shows only a 6% di erence. However, the moment derivatives show signi cant deviations, suggesting that the unmodeled aerodynamic e ects in uence the rotational dynamics more than the translational dynamics. Rotational lift, in particular, is known to be a signi - cant component of ight forces and causes force peaks that occur at stroke reversal where the wing has the largest moment arm and is most able to e ect pitch mo- tions (Dickinson et al., 1999). Mu is most poorly estimated (a factor of 6) by the experimentally- tted model. 101 ?2 ?1 0 1 2?2 ?1.5 ?1 ?0.5 0 0.5 1x 10 ?5 Pitch/pitch Rate, rad/s Stroke?Averaged Forces, N Xq QS = 0.00310 Zq QS = 0.000540 Xq CFD = 0.0148 Zq CFD = ?0.00105 XQS ZQS XCFD ZCFD (a) ?2 ?1 0 1 20 0.5 1 1.5 2 2.5 x 10?9 Pitch Rate, rad/s Stroke?Averaged Moments, N m Mq QS =?6.37 Mq CFD = ?15.74 MQS MCFD (b) Figure 5.4: Stroke-averaged forces (a) and moments (b) as given by the curve- tted experimental model (red) and IBINS (blue) for various perturbations in q. A comparison of each of system pole locations may be found by reference to Fig 5.6, showing again fast and slow subsidence modes with an unstable oscillatory pair. The decoupled pole Zw shows excellent agreement. However, the deviations in the moment derivatives have doubled both the rate of convergence in the fast subsidence mode and the rate of instability in the unstable oscillatory pair. The change in rate is due to the deviation in Mu, and is particularly signi cant because the speed of the unstable pair dictates the level of feedback required to stabilize the system. 102 Stability Derivative Curve-Fitted CFD Xu -5.49 -5.55 Zu 0.0 -0.285 Mu 730.1 4554.4 Xw 0.00122 -0.736 Zw -2.49 -3.70 Mw 931.4 1730.6 Xq 0.00310 0.0148 Zq 5.40e-4 -0.00105 Mq -6.37 -15.74 Table 5.1: Stability derivatives as calculated by the curve- tted model and IBINS. 5.5 Summary In this chapter, the curve- tted insect aerodynamics model used to derive reduced order ight dynamics models was compared to a numerical solution of the Navier-Stokes equations for incompressible ow, which includes unsteady and other unmodeled aerodynamic e ects. Finite di erencing through 5 points was used to estimate the derivatives with respect to the states u, w, and q and thus obtain stability derivatives for the derivatives. Translational stability derivatives agree well, while there are signi cant deviations in the moment derivatives. The system pole locations again show both a fast and slow subsidence mode and an unstable pair. However, the un-modeled aerodynamics increase the rate of both convergence and divergence in the system. In particular, the increased rate of divergence increases the feedback requirements for ight stabilization, the implications of which will be addressed in the next chapter. 103 Stroke Plane (a) Stroke Plane (b) Stroke Plane (c) Stroke Plane (d) Figure 5.5: Isosurfaces of Q-criterion of the model Drosophila in hover (a), surge u = 0:38 m=s (b), heave w = 0:38 m=s (c), and pitch q = 2:1 rad=s (d) just after pronation. The disturbed ow due to the previous wingstroke remains only a single chord length away from the returning stroke in all cases, and in uences the loads on the wing. 104 ?60 ?40 ?20 0 20?40 ?20 0 20 40 Real Imaginar y QS Poles CFD Poles Simulated Poles (Faruque & Humbert (2010a) Figure 5.6: Longitudinal poles of the model Drosophila in hover as calculated using the curve- tted experimental aerodynamics model (red) and CFD (blue). Poles found via system identi cation with the experimentally-derived model in Faruque and Humbert (2010) are plotted alongside for reference (green). 105 Chapter 6 Feedback E ects and Parameter Variation The longitudinal and lateral models developed in Chapters 3 and 4 are derived for a nominal Drosophila and for both bare airframe and haltere mediated dynam- ics. Section 6.1 explores the dynamic properties of a model with varying levels of feedback between these two extreme cases. Chapter 5 showed that di erences in the reduced-order model are primarily in rate of divergence, which increases the feedback requirements. Section 6.2 applies mechanosensory feedback to the haltere model to determine the magnitude of additional feedback required. While the model?s dynamic properties were shown to be preserved under the estimated uncertainty within the stability and control derivatives, morphological di erences between individuals of a species also a ect the dynamic properties of the model. Section 6.3 quanti es this variation by exploiting Drosophila?s sexual dimorphism to allow comparison of models derived for the nominal male insect and a dramatically larger female. 6.1 Mechanosensory Feedback Variation 6.1.1 Feedback Variation Study The desired outcome of this chapter is an understanding of how mechanosen- sory feedback can a ect reduced order insect ight dynamics models. First, analytic 106 functions were written for the stereotyped wing kinematics of a dipteran insect. The kinematic functions were used along with an insect aerodynamics model in quasi- steady form in order to predict the instantaneous forces and moments applied to an insect body (Dickinson et al., 1999). The quasi-steady aerodynamics form, nor- mally taken to be operating at a point, was extended to apply in the presence of egomotion. The equations of motion were posed as classical rigid body equations of motion, and linearized solutions determined. To investigate the e ect of mechanosensory feedback and the resulting control requirements, the insect models were derived both with and without control feedback provided by aerodynamically ine ective \hind wings" known as halteres, which are thought to encode rate information (Thompson et al., 2008). The haltere model used was the angular rate feedback with bandpass ltering used in Chapters 3 and 4. 6.1.2 Gain Variation Results The results of Chapters 3 and 4 include reduced order models of the insect ight dynamics. The longitudinal motion of the example Drosophila-like insect was written as _x = Ax + Bu, where x is the vector of longitudinal states composed of forward ight velocity u, heave motion w, pitch rate q, and pitch angle ; and u is the vector of kinematic inputs composed of stroke amplitude , stroke plane inclination , and stroke bias o as de ned in Figure (6.1) The system dynamics 107 (a) Stroke plane angle C o n tr o l I n pu t Se n s it iv it y D ro so ph ila lo n gi tu di n a l c on tro l i n pu ts ab o u t h o ve r: Co n tr o lla bi lit y: Pi tc h ra te dr ive n by In de pe n de n t h e a ve co n tro l in pu t Fu lly Co n tr o lla bl e Fo re /a ft m ot io n dr ive n by ? o ff ? ? (b) Stroke angle amplitude, o set. Figure 6.1: Kinematic inputs de nitions. and control input matrices A and B may be written for the haltere-on case as A = 2 6 6 6 6 6 6 6 6 6 6 4 12:33 0 0 9:81 0:0 4:651 0 0 547:0 0:0 33:26 0 0 0 1 0 3 7 7 7 7 7 7 7 7 7 7 5 (6.1) B = 2 6 6 6 6 6 6 6 6 6 6 4 0 20:67 12:69 17:36 0 0 0 2826:0 6020:0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 5 : (6.2) The same procedure applied to the haltere-o dynamics also reduces the longitudinal motion to a linear time invariant system. The ease with which a dynamic system (in this case an insect) naturally returns or diverges to a reference ight condition is described mathematically via the rigid body modes (characteristic motions) of the system and the associated modal poles. 108 ?40 ?30 ?20 ?10 0 10 ?15 ?10 ?5 0 5 10 15 Real Imaginar y Haltere?off Poles Haltere?on Poles Figure 6.2: The stability properties of the insect dynamics are improved by haltere feedback, and a root locus plot over haltere gain shows pole movement to connect the systems. In a pole-zero diagram, pole locations along the x-axis (real axis) represent the speed with which the particular motion either returns to the reference trajectory or diverges from it, and the y-axis (imaginary axis) locations quantify the amount of oscillation in the motion. The pole locations of the longitudinal system, seen in Fig 6.2, show that the insect with no haltere feedback (seen in red) has an unstable mode, but the pole locations of the haltere-on system (seen in blue) are stable. The nonlinear haltere model used in simulation was a band pass lter tted to data by Sherman and Dickinson (2004), but a simpli ed model of the haltere pitch rate feedback is u = kq, where k is the gain of the haltere. Given the haltere o system dynamics Ao , the system dynamics matrix under increasing gain on the 109 halteres may be written as Aon = Ao BK, where K = 2 6 6 6 6 6 6 4 0 0 0 0 0 0 0 0 0 0 k 0 3 7 7 7 7 7 7 5 (6.3) and the haltere-on matrix A is recovered at full haltere gain. In Fig 6.2, the simpli- ed model for pole locations for under varying haltere gains indicates the progression from the unstable haltere-o system to the stable haltere-on system, and that the os- cillation in the pitch and fore/aft degrees of freedom is reduced, while the decoupled pole (corresponding to heave damping) is una ected by rate feedback. Given that unstable systems are inherently more di cult to control than sta- ble systems, linear systems analysis indicates that the addition of mechanosensory feedback in the form of haltere rate feedback can signi cantly reduce the insect?s sensing and feedback requirements by providing a stabilizing e ect on the motion in the pitch and fore/aft degrees of freedom. 6.2 Mechanosensory Feedback on CFD model Section 6.1 investigated the e ect of gain reduction on a linear system, showing a progression from stabilized (haltere-on) to the bare airframe instability. Dipteran insects possess mechanosensory feedback beyond the haltere feedback that has been measured experimentally and was modeled in the nonlinear simulation of Section 3.3. Since the computational uid dynamics investigation in Chapter 5 suggests 110 ?80 ?60 ?40 ?20 0 ?30 ?20 ?10 0 10 20 30 Real Imaginar y Bare airframe Modeled feedback High gain feedback Figure 6.3: Stabilization of the CFD-derived system (in red) requires a higher pitch rate gain than than originally modeled (in blue). Equivalent convergence rate re- quires a gain 2.45 times larger (in black). that reduced-order modeling underestimates the instability, the feedback required to stabilize the CFD-derived system could signi cantly exceed the feedback investi- gated for reduced order modeling. More recent work has suggested that other forms of sensing are thought to provide rate feedback as well, in particular ocelli (Rowell and Pearson, 1983). Ocelli are rudimentary eyes (illumination level detectors) that are well-suited for high angular rate measurement. In order to quantify the additional mechansensory feedback required to sta- bilize the CFD system, a similar approach to section 6.1 was taken, using the lon- gitidunal dynamics model with the haltere gain modeled as a linear feedback law u = Kq such that Acfd,on = Acfd BK. Eigenvalue computation then shows that the modeled haltere pitch rate gain is insu cient to stabilize the CFD-derived sys- tem, and is only su cient to move the unstable pole locations from Re( ) = +12:1 to +3.5, as seen in Figure 6.3. 111 In particular, achieving marginal stability (Re( ) = 0) requires a gain 1.57 times greater, but marginal stability is rarely implemented in practice. To achieve an equivalent convergence rate as the reduced order model (Re( ) = 3:5), the CFD- derived system requires a pitch rate gain 2.45 times the originally modeled gain. This additional feedback suggests that the experimental measurements of haltere- mediated feedback have only quanti ed 41% of the mechanosensory feedback that is active on the insect, and that further experimentation to isolate other forms of feedback is necessary. 6.3 Parameter (Phenotypic) Variation The genotypic limits on species such as Drosophila allow remarkable pheno- typic plasticity. Environmental factors (e.g. ambient temperature) are relevant both during development and while mature, and morphology is highly dependent on both genetic limits and environmental factors. In animals kept at temperatures of 14 to 21 C, body mass varied from .9 mg to 1.8mg (Karan et al., 1998). The insect used for hover modeling, with its mass of m = 1:02 mg and wing length of R = 2:12 mm, uses parameters based on the males presented in the results of Karan et al. (1998); Crill et al. (1996). To compare how the results are a ected by the larger (approximately 150%) morphology presented by a female insect, lateral-directional hover dynamics for a female Drosophila-like insect having a mass mf = 1:70 mg and Rf = 2:355 mm was also system identi ed, and the results, shown in Table 6.1, indicate that the parameter variation generated by this relatively large morphological change is 112 Parameter Value Standard Dev Insensitivity Yv -13.44 6.63% 2.95% L0v -9658 .07% .03% L0p -178 0.2053% .09% N 0r -70.55 0.1648% 0.08% L d -340.73 2.9% 1.3% N d -28,080 3.9% 1.9% Table 6.1: Control and stability derivative results for female Drosophila. ?200 ?150 ?100 ?50 0?25 ?20 ?15 ?10 ?5 0 5 10 15 20 25 Real Imaginar y Nominal (male) poles Larger (female) poles Figure 6.4: Comparison of the lateral-directional poles identi ed for the nominal (male) and larger (female) insect. contained primarily in the control derivatives. Accordingly, the modal behavior of the system shows a relatively minor variation, as exhibited by the pole shift from ( 1; 2,3; 4) = ( 73; 3:4 22:9i; 180:5) to ( 70:6; 5:2 22:2i; 181:1), seen graphically in Fig. 6.4. 6.4 Summary In this chapter, the e ect of a gradual reduction in mechanosensory feedback was considered, showing a progression from the stable haltere on system to the 113 unstable bare aiframe system. The feedback requirements of a system including detailed aerodynamic models were shown to be greater, requiring an increase in feedback gain by 2.45 times to achieve an equivalent convergence rate speci cation. Finally, the degree to which individual morphological di erences within a species a ect the dynamic properties of apping ight were addressed by using the species? sexual dimorphism to consider insects at either end of the commonly observed size range for Drosophila. Di erences were found primarily in the size of the control derivatives, but the modal structure was retained. 114 Chapter 7 Reachability Analysis This dissertation has developed a set of kinematic programs suitable for ma- neuvering an insect or apping wing micro air vehicle (MAV), and ight dynamics models to relate those inputs to the rigid body motion of the vehicle. In several cases, multiple inputs exist to e ect the same desired trajectory. MAV design presents size, weight, and power limitations that reduce actuation capability. Given more inputs than necessary to control a dynamic model, one may consider the question of which inputs are \optimal" in some sense, either from the perspective of minimizing actu- ation e ort for a given motion, or maximizing the set of all reachable states under a restricted input. This chapter introduces a control theoretic framework to quantify the reach- able states for a given set of inputs, and applies it to the hovering fruit y model developed in Chapter 3. Controllability as an application of operator theory is the basis for determining the reachable con gurations under a class of inputs. The ex- pressions for reachable states may then be used to solve a least-squares optimization problem over all possible function inputs. Previous analysis of insect-inspired apping wing locomotion has examined wing kinematic trajectories from the perspective of maximizing lift(Avadhanula et al., 2003; Ansari et al., 2008) or minimizing required power (Berman and Wang, 115 2007). With the introduction of new tools to extract the detailed wingstroke to wingstroke kinematics of insects from high speed videography (Fontaine et al., 2009; Ristroph et al., 2009), a number of species-speci c control strategies for maneuvering have been identi ed. In addition, with the development of micro-scale vehicles that can potentially generate lift forces greater than their weight (Wood, 2008), stability and control aspects of the problem have become an important research need. Despite the critical need, the inherent complexity of small-scale apping ight aerodynamics has obscured a control-theoretic analysis of both biologically relevant and engineered wing kinematic perturbation strategies. While the detailed aerody- namic mechanisms involved in small-scale ight are still an area of active research (Ramamurti and Sandberg, 2007), recent e orts have in fact yielded several ap- proaches for extraction of reduced-order linear time-invariant ight dynamics, either for single degree of freedom experimental cases (Hesselberg and Lehmann, 2007), direct analytic methods (Doman et al., 2010a), or more general computationally (Sun et al., 2010) and spectrally derived models as presented in Chapters 3 and 4. Such formulations are amenable to application of linear control analysis tools, and should provide the next level of insight. Reachability (or more traditionally, controllability) characterizes the amount of control one has over the state of a system through the choice of the input. This is an important topic for small-scale apping wing MAV designers for several rea- sons. Size, Weight, and Power (SWaP) constraints are very stringent at this scale, and reductions in complexity that promote weight reduction or robustness are en- couraged. In addition, these vehicles are intended to operate in gusty and possibly 116 cluttered environments, and a high level of platform maneuverability and actuation authority will be crucial to achieving robust ight path control in the face of these uncertainties. This chapter explores the reachable state space associated with biologically- inspired kinematic control strategies seen in fruit y longitudinal motion about hover (see Chapter 3 for development of the model). In Section 7.1 a frequency-based system identi cation methodology for identifying the stability derivatives of a small- scale apping microsystem about hover is outlined, along with the control derivatives for biologically relevant wing kinematic perturbations for maneuvering. Section 7.2 applies controllability analysis tools to interpret these biologically-relevant control strategies for micro-air vehicle design, using the example of an MAV with Drosophila- like parameters. 7.1 Longitudinal Flight Dynamics Model Previous chapters have outlined a method to formulate linear time invariant (LTI) ight dynamics models of the form _x = Ax+Bu: (7.1) In this case A 2 Rn n and B 2 Rn p represent the state and input matrices, u(t) 2 U Lp2[0;1) the input time history and x(t) 2 X L n 2 [0;1) the state history of the model. The nominal (trim) kinematics are assumed to be a periodic oscillation (t) 117 contained in a stroke plane inclined an angle from horizontal, with the wing in- tersecting the stroke plane at an angle g(t). Biologically relevant control inputs considered in this study are shown in Fig 7.1 and de ned as (a) Stroke plane incli- nation (Fry et al., 2003) c: a tilting of the stroke plane generating pitch moment and forward force, (b) Stroke plane o set (Fry et al., 2003; Oppenheimer et al., 2010; Doman et al., 2010b) o : a fore/aft shift of the wing sweep used primarily to generate pitch moment, and (c) Asymmetric wing angle (Ristroph et al., 2009) ud: an upstroke/downstroke asymmetry in the angle of the wing relative to the stroke plane used primarily to generate forward force. ?c = 12(?r + ?l ) ?ud = 1 2 (?u ? ?d)?o? = 12(?o?r + ?o?l ) ?o?l ?o?r ?d ?u ??r +?l CBA Figure 7.1: Longitudinal control inputs used in reachability analysis. (A) Stroke plane inclination c, (B) Stroke plane o set o , and (C) Di erential wing angle ud. The stability and control derivatives in the ight dynamics model from Chapter 3 were selected to maximize coherence in the low frequency regions (up to 20 Hz), while discarding the small periodic high frequency motion. Note that with the haltere feedback about the pitch axis, the longitudinal system is stable and the matrix A is Hurwitz. The general state space model structure for the longitudinal dynamics includ- 118 ing the control inputs listed above is given by 2 6 6 6 6 6 6 6 6 6 6 4 _u _w _q _ 3 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 4 Xu 0 0 g 0 Zw 0 0 Mu 0 Mq 0 0 0 1 0 3 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4 u w q 3 7 7 7 7 7 7 7 7 7 7 5 + 2 6 6 6 6 6 6 6 6 6 6 4 0 X X o X ud Z 0 0 0 0 M M o M ud 0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4 c c o ud 3 7 7 7 7 7 7 7 7 7 7 5 :(7.2) The general state space model includes an additional input, stroke amplitude c. Stroke amplitude a ects only heave ( w) dynamics and is decoupled from the other states and inputs. 7.2 Reachability Analysis For the longitudinal ight dynamics (7.2), the goal is to have a rigorous frame- work in which to quantify the e ectiveness of the biologically relevant control strate- gies (Fry et al., 2003; Ristroph et al., 2009) that have been described in the previous section. As the heave ( w) dynamics and the collective stroke amplitude control input ( c) are decoupled from the remaining states in the linearized model, the pitch/fore/aft dynamics ( u; q; ) and control inputs ( c; o ; ud) are consid- ered without loss of generality. Application of the controllability rank test for all possible combinations of control inputs reveals that the pitch/fore/aft system is fully controllable with any pair of the remaining control inputs, which motivates the examination of three input 119 pairs: u1 = 2 6 6 4 c o 3 7 7 5 ; u2 = 2 6 6 4 c ud 3 7 7 5 ; u3 = 2 6 6 4 o ud 3 7 7 5 : (7.3) In order to minimize actuator e ort and maintain a small number of controls, we desire the pair of inputs that maximize the span of reachable states x0 resulting from any arbitrary input u(t) 2 Lp2( 1; 0] of unit norm. 7.2.0.1 Controllability operator Consider the linear system _x = Ax+Bu; (7.4) with A 2 R3x3 and B 2 R3x3. De nition 7.2.1 (Reachability) The system fA;Bg is de ned as reachable over the interval [t0; t1]; t1 > t0 if for every pair of states x0; x1 2 X , there exists a control u(t) 2 Lp2[0;1) such that the solution of _x = Ax+Bu; x(t0) = x0 satis es x(t1) = x1. Following Corless and Frazho (2003); Dullerud and Paganini (2000) for a time invari- ant system, t0 may be chosen arbitrarily. Choosing t0 = 1, t1 = 0, the solution 120 to eqn (7.4) is x(0) = e A1x0 + Z 0 1 e A Bu( ) d : (7.5) De ne the operator from U = Lp2[0;1) to X as u = Z t1 t0 eA(t1 )Bu( ) d and for our choice of integration limits, operates from Lp2[0;1) to X as u = Z 0 1 e A Bu( ) d (7.6) Operator takes as an input a time history u(t) and outputs a state x0. More speci cally, it returns the nal state x(0) corresponding to eqn (7.4) with initial condition x( 1) = 0 and forced by u(t). This view motivates a consideration of what input u(t) returns a particular desired state x(0) = x1. If the system is controllable, then u(t) exists by direct ap- plication of the de nition. However, maps an in nite dimensional space Lp2[0;1) onto a nite dimensional space X , thus its kernel has in nite dimension and u(t) is not unique. Given the non-uniqueness, one may consider the case of what u(t) is \optimal." Consider the following least squares optimization problem for a (not necessarily controllable) system: Given the initial condition x( 1) = 0 and a desired nal state x1, nd an input u^(t) that satis es 121 jju^(t)jj2 = inf Z 0 1 jju(t)jj2dt : x^(0) = Z t1 1 e A Bud where x^(0) is the unique vector in X satisfying (7.7) jjx1 x^(0)jj = inf fjjx1 x(t1)jj : _x = Ax+Bu and x( 1) = 0g : In the general case, the minimization problem is to nd the input of least norm that drives the system as close as possible to the nal state. The terminal condition in eqn (7.7) is always satis ed for a controllable system, where a u(t) is known to exist and thus x^(0) = x1. For the controllable pair fA;Bg, then (i) the matrix =: Xc is nonsingular, and (ii) for any x1 2 X , the input u^(t) = cX 1 c x1 is the element of minimum norm in the set fu 2 Lp2[0;1); u = x1g : For a detailed proof of (i) and (ii), see Dullerud and Paganini (2000). 7.2.1 Reachable space under unit norm input In the process of MAV control design, a measure of how \far" inputs may drive the system in the con guration space can assist the choice of a control strategy. Mathematically, the con guration space that is reachable under unit norm input, 122 expressed as f u : u 2 Lp2[0;1) and jju(t)jj 1g is equivalent to Ec = n X 1 2 c xc : xc 2 X and jjxcjj 1 o : To verify the equivalence, one may show that each set is contained in the other. Ec de nes an ellipse in Cn whose geometric properties are determined by the in nite- time controllability gramian Xc for the system (7.2), Xc = c c = Z 1 0 eA BB eA d 0; (7.8) which can be computed given the matrix pair (A;B) via the Lyapunov equation, AXc +XcA +BB = 0: (7.9) The principle axes and lengths of Ec are given by the eigenvectors and eigenvalues of X 1 2 c , respectively, which motivates two control input ranking criteria. The rst is the Frobenius norm of X 1 2 c , jjX 1 2 c jjF = r trace h (X 1 2 c ) X 1 2 c i ; (7.10) whose geometric interpretation is the square root of the summed squares of the axes lengths of Ec. Second, since Xc 0, one can also consider the volume det(X 1 2 c ) of Ec as a non-negative measure of its size. Choosing control degrees of freedom 123 that maximize either of these measures then corresponds to maximizing the set of reachable states over the choice of control degrees of freedom. The controllability ellipsoids of the input pairs are shown in Fig. 7.2A-C, indicating that the reachable space for a unit norm input increases signi cantly from input pair u1 through u3. For comparison, the reachable con guration space for a unit norm input on all three control terms ( c; o ; ud) is shown as u4 in Fig. 7.2D. The results of applying the two ranking criteria to the control input selections u1 through u4 are shown in Table 7.1 and plotted in Fig. 7.3. Out of the three pairs, clearly u3 = ( o ; ud) provides the most authority over the longitudinal dynamics. In terms of the reachable volume measure, u3 provides a 672% increase over u1 and a 38% increase over u2. In addition, utilizing all three control inputs u4 only provides for a modest 6% increase in reachable volume over u3. Similar conclusions follow from the Frobenius norm ranking criteria; u3 provides 94% and 11% increases over u1 and u2 respectively, whereas u4 adds a 2% improvement over u3. The ellipsoidal interpretation also yields important information regarding the resulting system?s controllability along particular directions in state space. In the case of Fig. 7.2, the rotational dominance of the control inputs (and modes) is evident in the larger reachable con guration space along the pitch rate/angle axes, while the range of forward speed is more limited. 124 DC BA 2 3 1 4 ?1 0 1 ?20 0 20 ?2 ?1 0 1 2 B od y P itc h A ng le (rad ) ? B od y P itc h A ng le (rad ) ? ?1 0 1 ?20 0 20 ?2 ?1 0 1 2 ?1 0 1 ?20 0 20 ?2 ?1 0 1 2 ?1 0 1 ?20 0 20 ?2 ?1 0 1 2 Pitch Rate (rad/s) Fwd Speed (m/s)q u u u u u Figure 7.2: Controllability ellipsoids for the input combinations illustrate the reach- able con gurations under the restriction jjuijj 1. Input combinations u1 through u3 are pairs of control terms, while u4 considers all 3 control terms ( c; o ; ud). 0 10 20 30 40 50 u2 u3u1 u4 kX 12c kF det(X 12c ) Input Combination Si ze Figure 7.3: Controllability of input combinations u1 through u4, as ranked by the determinant or Frobenius norm of the square root of the controllability gramian. 125 Input u1 u2 u3 u4 det(X 1 2 c ) 5.87 32.65 45.31 48.06 jjX 1 2 c jjF 17.04 29.80 33.03 33.69 Table 7.1: Control input performance ranking criteria for u1 though u4. 7.3 Summary This chapter introduces a control-theoretic framework and two performance measures to quantify the state reachability for a given choice of biologically inspired wing kinematic control parameters, allowing for a set of candidate control strategies to be ranked. The work presented here leverages the Drosophila melanogaster linear time-invariant (LTI) ight dynamics models previously developed in this disserta- tion. For the example insect-size micro air vehicle (MAV) considered, all four of the input combinations provided full controllability, but the reachable space is dramat- ically improved through proper selection of the input combinations. The kinematic set that provided the most controllability involved a stroke bias term ( o ) and an angle of attack di erence in the upstroke and downstroke ( ud), which improved the reachable volume of state space 672% over the least controllable set. Moreover, while the reachable space is dramatically improved over each of the input pairs, only a slight advantage is found by combining all three inputs; adding a stroke plane angle degree of freedom ( c) to the most e ective pair resulted in only a modest increase (6%) in the volume of reachable state space. The framework and performance measures introduced in this note provide a means to appropriately choose kinematic inputs that minimize the required control 126 energy and maximize the achievable state space of the system. These tools can be applied to reduce actuator complexity to promote robustness and weight reduction, resulting in improved size, weight, and power (SWaP) requirements for micro-air vehicle ight stability and control. For MAV design, factors other than control energy (such as actuator geometry) may be the limitation on kinematic actuation, and a system-level approach must be used to determine the limiting factor. The framework introduced in this chapter is directly applicable for a systems level model to be used in MAV design studies. 127 Chapter 8 Forward Flight Dynamics This chapter uses the Calliphorid kinematic measurements of Chapter 2 and the approach developed in Chapter 3 and 4 to determine forward ight dynamics models for dipteran ight, as opposed to the hover-oriented fruit y models presented previously. 8.1 Background This dissertation has investigated hovering ight dynamics for dipteran ap- ping ight. An example Drosophila melanogaster was built up using nominal kine- matics for hover, an experimentally-derived aerodynamic model, and frequency- based system identi cation on a rigid body dynamics simulation. Analysis of the models led to conclusions about sensing and feedback requirements in hover. How- ever, considerable research has been conducted on insect visuomotor responses in an e ort to understand insect navigation and guidance algorithms. A concise forward ight dynamics model would allow researchers to place the feedback measurements into the context of navigation and guidance, and to understand why insects mea- sure non orthogonal quantities. With experimental models under development that include rigid body rotation and translation (Dickson et al., 2010), researchers are approaching the capability to compare forward ight data with modeling results. 128 Chapters 3 and 4 presented a method of deriving ight dynamics models about hover, and a natural extension of this work could be to look at the ight dynamics about another reference condition. This chapter uses the same approach to in- vestigate the longitudinal and lateral-directional dynamics of blow ies in forward ight, and present an expectation for the results of high delity experimental and computational measurements. Biological kinematics may also play a signi cant role in ight stabilization. Ex- perimental evidence has indicated that dipteran insects in forward ight may have turn rate limitations that are theorized to be correlated to aerodynamic mechanisms (Bueltho et al., 1980). While xed-wing ight is governed by aerodynamic mecha- nisms such as wing, fuselage, and tail interaction (Nelson, 1989), a dipteran apping vehicle?s primary aerodynamic component is its wings and small changes in wing kinematics have a dramatic e ect on the maneuvering of the vehicle. In forward ight, many insects employ di ering wing pitch angles in the fore and aft strokes, typically reducing the morphological angle of attack in the forward stroke and in- creasing it in the retreating stroke (Zanker and Gotz, 1990; Taylor and Thomas, 2003a,b). This upstroke/downstroke asymmetry may be an attempt to preserve a constant angle of attack or a physiological limitation (Zanker, 1990), but it may also have rami cations on the control requirements of forward ight. This chapter will also examine the e ect of the wingstroke perturbation in lateral-directional ight. 129 8.1.1 Previous Work Previous work on forward ight dynamics modeling has been limited. Data were collected for tethered Orthoptera over which air ow was forced (Taylor and Thomas, 2003a,b), to simulate forward ight. Linear time invariant system theory was used to analyze the resulting behavior but the e ect of the insect?s control structures (which were active) were not addressed. Experimental apparatus capable of both apping and egomotion are under current development, and this is expected to soon become a research focus. There is a need for reduced order models with which to interpret the data collected in this emerging research area. 8.2 Longitudinal Flight Dynamics Modeling This section presents the longitudinal ight dynamics modeling, including the e ect of the more complex g function de ned in Chapter 2. 8.2.1 Kinematics For this study, untethered Calliphorid ies in straight and level forward ight at a mean cruise of 2.17 m/s were used to determine stereotypical forward ight wing kinematics, as measured in Chapter 2. Note that the wingstroke has di ering wing pitch on the fore and aft strokes. Mean g values in each stroke are 58:7 and 54:0 , shown in Fig 8.1. These kinematics are de ned as trim inputs, and pertur- bations about those kinematics are used to maneuver a simulated insect in forward ight. Again, symmetry allows right and left wing motions to be parameterized by 130 Fwd Fwd Fwd mean ?60 ?40 ?20 0 20 40 60 80 100 Angle, de g 0 0.1 0.2 0.3 0.4 0.5 Time ?g ?g ? (a) Fwd Aft Aft Aft mean ?g ?g ? ?60 ?40 ?20 0 20 40 60 80 100 Angle, de g 0 0.1 0.2 0.3 0.4 0.5 Time (b) Aft Figure 8.1: Kinematics in forward and aft strokes. longitudinal (collective) motions. 8.2.2 Aerodynamics There are indications that forward ight will require further modi cations to the experimental aerodynamics model (Dickson and Dickinson, 2004). In this dis- sertation, the aerodynamics model has been placed in the context of insect body translation and rotation. Here we derive control-oriented models under the assump- tion that the aerodynamics model including egomotion is valid for cruise ight. When high delity data is available, either via computational solutions like those in Chapter 5 or experimental equipment, the validity of this assumption will be investigated and a similar aerodynamics model validated in forward ight may be applied. 131 8.2.3 Rigid Body Dynamics and System Identi cation Rigid body dynamics and system identi cation is not a ected by the lineariza- tion about a di erent reference ight condition, thus simulation and system identi- cations may be conducted as before. The Calliphorid insect was rst trimmed at o ,c = 14:375 , Vnorth = 2:18m=s; Vup = :06m=s; = 1:58 ; u = 2:18m=s;. Then, 60 second frequency sweeps were applied including frequency content up to 50 Hz. Longitudinal input sweeps were conducted for c, o ,c, and c. However, system identi cation of the unstable system dynam- ics was complicated by the need to apply high gain feedback to stabilize the pitch dynamics, which reduces the signal to noise ratio of the identi ed transfer function by introducing correlated signals into the input. The identi ed forward ight longitudinal dynamics model is A = 2 6 6 6 6 6 6 6 6 6 6 4 1:32 0 0 9:81 25:2 24:7 0 0 1011:0 341:0 12:6 0 0 0 1:0 0 3 7 7 7 7 7 7 7 7 7 7 5 (8.1) B = 2 6 6 6 6 6 6 6 6 6 6 4 71:9 13:7 2:86 320:0 60:9 0 0:0357 2211:0 5955:0 0 0 0 3 7 7 7 7 7 7 7 7 7 7 5 : (8.2) This system has poles at 4:8072 17:0857i, -32.6466, and -15.5768, mirroring 132 ?25 ?20 ?15 ?10 ?5 0 5 10 ?15 ?10 ?5 0 5 10 15 Real Imaginar y Numerical Identified Figure 8.2: The pole structure of is similar to those in hover, with a real heave mode, and a coupled rotationally dominated triple that includes an unstable oscillatory pair. the previous structure, with two real subsidence modes and unstable oscillatory pair, as seen in Fig. 8.2. The transfer function plots (see Fig 8.3 and 8.4) and t convergence criteria indicate that the traditional model structure leads to error in the heave w dynamics. Table 8.1 presents state-space derivatives for the plots, uncertainties lower than 10% for all but the heave and axes, indicating again that model structure re nements will be necessary. 8.3 Lateral-Directional Forward Flight Dynamics In this section, lateral models are derived for forward ight. The models include di ering wing pitch angles on the advancing and retreating wingstrokes, an input strategy observed in experimental studies (Zanker, 1990). This section introduces reduced order models for roll and yaw dynamics in forward ight and 133 ?30 ?15 0 15 30 Mag, d B ?180 ?90 0 Phase, de g 101 102 0 0.5 1 Coherenc e Frequency, rad/s (a) o to u ?15 0 15 30 Mag, d B ?90 ?45 0 45 90 135 180 Phase, de g 101 102 0 0.5 1 Frequency, rad/s Coherenc e (b) o to w 20 30 40 50 Mag, d B ?180 ?135 ?90 0 Phase, de g 101 102 0 0.5 1 Frequency, rad/s Coherenc e (c) o to q ?15 0 15 30 Mag, d B 0 90 180 Phase, de g Simulation Identification 101 102 0 0.5 1 Frequency, rad/s Coherenc e (d) o to Figure 8.3: A 4 state linear system o transfer function shows good agreement for all but the heave direction, which shows model structure error. 134 ?30 ?15 0 15 30 Mag, d B 0 90 180 270 Phase, de g 101 102 0 0.5 1 Coherenc e Frequency, rad/s (a) to u ?15 0 15 30 Mag, d B ?90 ?45 0 45 90 135 180 Phase, de g 101 102 0 0.5 1 Frequency, rad/s Coherenc e (b) to w 0 15 30 45 Mag, d B 0 90 190 Phase, de g 101 102 0 0.5 1 Frequency, rad/s Coherenc e (c) to q ?15 0 15 30 Mag, d B ?90 0 90 Phase, de g Simulation Identification 101 102 0 0.5 1 Frequency, rad/s Coherenc e (d) to Figure 8.4: A 4 state longitudinal linear system transfer function shows good agreement for all but the heave direction, which shows model structure error. 135 Parameter Value Cramer Rao % Xu -1.605 5.638 Zu -8.590 5.756 Zw -7.433 7.300 Mu 569.1 1.749 Mw 17.29 87.98 Mq -12.39 4.110 X 13.20 5.267 X o -2.945 4.469 X 33.04 11.46 Z -44.06 5.315 Z o -55.00 2.763 Z -4.661E+03 177.0 M -2.175E+03 5.492 M o 6.129E+03 1.462 M 2.816E+03 17.58 Table 8.1: Uncertainty values in the parameters are acceptable, with the exception of the w derivatives that poorly identi ed. investigates how the upstroke/downstroke kinematic asymmetry can act to reduce the sensing and control requirements in forward ight by appropriately coupling lift and drag components in a turn. To derive the equations of motion, the insect body was modeled as rigid, allowing Euler rigid body dynamics to be applied. After a trim solution was determined that would provide a forward ight condition, the aerodynamic forces and moments as calculated by the experimental aerodynamics model were applied to the vehicle, which was constrained to motion at varying forward speeds from u = 0m/s to 2m/s but allowed to move about the axis of interest (roll or yaw) as dictated by aerodynamics to yield linear models. In contrast to earlier studies, no visual or mechanosensory feedback terms were included and the modeling represents bare airframe response. 136 8.3.1 Kinematics and System Identi cation The nominal angle of attack pattern used to derive a lateral-directional ight modely is a modi ed square wave. Insects in forward ight commonly apply a decrease in wing pitch in the forward stroke and an increase in the aft stroke, which is modeled as a bias (or o set) term ud applied to the angle of attack pattern in forward ight. Experimental studies (Zanker, 1990; Zanker and Gotz, 1990) have indicated that diptera typically use negative values of ud in forward ight. The wing pronation/supination for this study is given by g;r(t) = 45 tanh [2:7 sin (2 ft+ )] + ud,r (8.3a) g;l(t) = 45 tanh [2:7 sin (2 ft+ )] + ud,l; (8.3b) where r and l subscripts correspond to right and left wings. With suitable parameterizations for stroke and wing pitch angles, the control inputs may be considered the stroke amplitude , o set o , and pitch asymmetry ud. As before, right and left inputs are decomposed into collective and di erential inputs using eqn (3.2). To facilitate reduced order modeling, linear models were generated from time histories created by the simulation described in Chapter 3. The time histories used for analysis were generated by application of frequency sweeps from 0.1 Hz to 32 Hz to the d and d inputs to excite roll and yaw dynamics, respectively. The spectral content of the input Gxx(!) and output Gyy(!) (calculated by the chirp z 137 transform may be used to nd the transfer function G describing the input/output relationship. The linear models generated represent the relationship between the signals in Fig 3.7 and are directly applicable to estimating the sensing and control feedback requirements of dipteran forward ight. 8.3.2 Results and Discussion 8.3.2.1 Roll dynamics Passive aerodynamic mechanisms act to damp out roll rate via changes in angle of attack, as discussed in Section 4.5.2.2. Application of equation (4.1) to the chirp signals simulated for the example Calliphora insect allows computation of the transfer function seen in Fig 8.5, computed for the trimmed forward ight speed of 2m/s. Our desired result for roll dynamics is a model of the form _x = Ax+Bu: d dt 2 6 6 4 p 3 7 7 5 = 2 6 6 4 Lp 0 1 0 3 7 7 5 2 6 6 4 p 3 7 7 5+ 2 6 6 4 L d 0 3 7 7 5 d : (8.4) The excellent coherence over a large portion of the frequency range of interest (0.1 to 32 Hz) indicates largely linear behavior, and a linear system has been t to the transfer function in Fig 8.5. Time domain veri cation conducted using repeated doublets, typically a demanding input signal, shows excellent agreement, seen in Fig 8.6. The roll dynamics terms are only slightly changed when the identi cation was repeated at di ering forward velocities. Table 8.2 summarizes the small variations 138 30 40 50 60 70 Magnitude, dB Nonlinear Linear ?270 ?225 ?180 Phase, de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (a) Roll rate p 0 20 40 60 Magnitude, dB 0 30 60 90 Phase, de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (b) Roll angle Figure 8.5: In forward ight u = 2 m/s, rst order linear roll damping is an accept- able spectral description of the roll transfer function. ?0.2 ?0.1 0 0.1 0.2 ? d, de g 0 0.2 0.4 0.6 0.8 1?1.5 ?1 ?0.5 0 0.5 Roll rate p, rad/ s Time, sec Nonlinear Linear (a) Roll rate p ?0.2 ?0.1 0 0.1 0.2 ? d, de g 0 0.2 0.4 0.6 0.8 1?3 ?2 ?1 0 Roll angle ?, de g Time, sec Nonlinear Linear (b) Roll angle Figure 8.6: Time-domain agreement of the response to d input doublets in forward ight is excellent. 139 Speed u Parameter Value Standard Dev. Insensitivity 0m/s Lp -7.07 8.0% 3.0% L0 d -13400 5.0% 1.9% 1m/s Lp -7.37 8.0% 3.0% L0 d -13692 5.0% 1.9% 2m/s Lp -7.72 8.0% 3.0% L0 d -13822 5.0% 1.9% Table 8.2: Uncertainty estimates for parameters in the identi ed dynamics model show that roll damping and roll authority is only weakly a ected by increasing u. Speed u Parameter Value Standard Dev. Insensitivity 0m/s Lr -5.71 8.6% 3.5% L0 d -2995 4.8% 1.9% 1m/s Lr -5.26 8.7% 3.6% L0 d -3160 4.7% 1.9% 2m/s Lr -6.16 8.2% 3.2% L0 d -4294 4.8% 1.9% Table 8.3: Yaw control authority increases slightly with forward speed, while yaw damping is relatively una ected. in forward ight and uncertainty estimates. 8.3.2.2 Yaw dynamics Experimental evidence has indicated that yaw damping is well modeled by a rst order linear relationship (Dickson et al., 2010; Hedrick et al., 2009). A similar procedure to part 8.3.2.1 was applied to nd transfer functions for yaw motions in forward ight, seen in Fig 8.7. Equivalent linear damping and control e ectiveness terms were also found for the three reference ight speeds (Table 8.3). The results indicate that while yaw damping is relatively una ected by forward ight speed, yaw control e ectiveness increases with speed. 140 20 40 60 Magnitude, dB Nonlinear Linear ?270 ?225 ?180 Phase, de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (a) Yaw rate r 0 20 40 60 Magnitude, dB 0 30 60 90 Phase, de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s (b) Yaw angle Figure 8.7: Yaw dynamics transfer function t at u = 0m/s. 8.3.2.3 Roll/Yaw Coupling (Proverse Yaw) While single input single output models can yield signi cant insight into how roll and yaw damping change with forward speed, aerodynamic roll/yaw coupling may exist in forward ight. As an example, the adverse yaw tendencies of xed-wing aircraft increase pilot workload, which has motivated the development of di erential ailerons and mechanical or electrical interconnections to reduce adverse yaw tenden- cies. In particular, an insect that must support additional neurological structure to counter adverse yaw tendencies may be less able to survive robustly in the wild. If passive aerodynamics can assist in turn coordination, the sensing and control feedback necessary for dipteran apping wing ight is reduced. Roll/yaw coupling in forward ight shows that a roll motion in forward ight induces a yaw motion in the same direction, acting to coordinate the turn, as seen in Fig 8.8. Furthermore, the magnitude of the induced yaw motion varies from 1/500 of the roll rate at ud = 6 to 1/3 the roll rate as a function of the input ud = 6 . 141 The physical mechanism of turn coordination can be understood by reference to the lift and drag polars of the quasi-steady aerodynamics form (Sane and Dickinson, 2002), where a di erential change in angle of attack on each wing also a ects the drag components of each wing to create a yaw moment in the same direction. The nding that an insect with negative ud values (wing pitch reduced on the forward stroke) in forward ight receives turn coordination from passive aero- dynamic behavior suggests that the kinematic perturbations seen in forward ight also act to reduce the control requirements during turning behavior. ?90 ?45 0 45 90 Phase, de g 100 101 102 0 0.5 1 Coherenc e Freq, rad/s ?30 ?15 0 Magnitude, d B ?6 deg ?3 deg 0 deg 3 deg 6 deg Figure 8.8: Roll rate to yaw rate coupling at 2m/s for di erening wing pitch input ud. Roll motions in forward ight induce smaller yet potentially signi cant yaw rates. The magnitude of the induced yaw rate is a function of ud. 8.4 Summary This chapter investigated the longitudinal and lateral-directional ight dynam- ics of dipteran forward ight. Due to instabilities, the identi cation was conducted 142 while a stabilizing controller was running. The longitudinal model showed good agreement for the rotational inputs, but the model structure that was appropriate for hover did not meet convergence criteria in forward ight, particularly in the heave dynamics. Nonetheless, the pole structure associated with the longitudinal dynamics was qualitatively similar to that in hovering ight, with two stable real poles and an unstable oscillatory pair. Lateral-directional studies examined roll and yaw dynamics during constrained forward ight. In agreement with previous work, linear models were found to be su cient to describe the roll and yaw dynamics for the example Calliphora insect (Dickson et al., 2010; Hedrick et al., 2009). Roll and yaw damping were similar for the three ight cases considered, u = 0; 1; 2m/s, while yaw control authority increased with forward speed. The study also examined the upstroke/downstroke asymmetry seen in the wing pitch angle during forward ight. The insect aero- dynamics model predicts that the kinematic input provides a roll-to-yaw coupling through passive aerodynamic mechanisms. The roll-to-yaw coupling can assist in providing turn coordination, reducing the active control that an insect?s neurological structure must support. The kinematic input ud may be similarly used in design to reduce the ight control demands of a dipteran apping wing micro air vehi- cle (MAV) operating in forward ight, and thus reduce its size, weight, and power requirements. 143 Chapter 9 Conclusions and Future Work 9.1 Conclusions In this dissertation, a path for developing reduced order ight dynamics mod- els for dipteran apping wing insects was presented. The analysis began with ex- perimental measurement and characterization of wing kinematic motions through automated processing of high speed video. The wing kinematic motions were used with an empirically derived insect-scale aerodynamics model including state pertur- bations. Rigid body dynamics and system identi cation were applied to generate reduced order models for insects in hover and in forward ight, using as examples Drosophila melanogaster and Calliphorid ies. The e ects of mechanosensory feed- back on the ight dynamics were also determined, showing that biologically observed rate feedback paths are appropriate for ight stabilization. Results were compared to a numerical Navier-Stokes aerodynamic solution to determine the e ect of the un- modeled aerodynamics, nding that the unmodeled aerodynamics increase the rate of response, and thus the required control. Within-species morphological di erences were studied, nding very similar ight dynamics between large females and small males. A framework for determining the reachable con guration space associated with kinematic input programs, and methods for ranking the input programs in terms of highest maneuverability was presented. 144 The work presented is both relevant to understanding the sensing and feed- back paths that are active in ying insects and to the design of micro air vehicles. Insects solve the ight stabilization and control problem in a manner that requires only small, specialized control structures that perform the essential feedback require- ments by measuring composite quantities along non-orthogonal axes. An engineering means to replicate this on micro air vehicles will provide a dramatic increase in the mobility of small scale aerial robotics. 9.2 Dissertation Contributions This dissertation examines insect ight dynamics modeling and contributes the following to the eld: ? Kinematics: The rst detailed measurement of both trim and maneuvering wing kinematics for Calliphorid ies were made. Curve ts to the nominal kinematics provide stereotyped kinematics, while perturbations to the kine- matic parameters serve as control inputs for ight dynamics modeling. While previous research quanti ed only errors inherent in the digitization measure- ment, reference models were used to quantify the accuracy of the full measure- ment path. A kinematic perturbation ud was quanti ed in forward ight. ? Aerodynamics: This dissertation introduced a method of placing the con- temporary insect aerodynamics models in the context of insect body transla- tion and rotation. Previously, the aerodynamics models were posed in quasi- steady form and did not include any state dependence, which prevents dynamic 145 analysis. When placed in the context of insect egomotion, the aerodynamic model now allows estimation of dynamic behavior. ? Flight Dynamics Models: Transfer functions for longitudinal and lateral- directional hovering ight dynamics of Drosophila were generated by spectral transformation of input and output time histories of a nonlinear simulation. Linear time-invariant (LTI) state space models were tted to the transfer func- tions as well, allowing convenient analysis and simulation, and directly identify the form and rate of feedback required for stabilization. The use of LTI mod- eling was supported by linearity measures (coherence). Models were tted up to 500 rad/s (89 Hz) for both fruit and blow ies, though coherence was degraded in many transfer functions over 200 rad/s (30Hz). ? Model Validation: Collaborative work with the authors of Bush, B. et al. (2010) compares the results of the dynamics modeling with the numerical so- lutions to the Navier-Stokes equations of motion, nding excellent agreement in translational derivatives, while improvement is suggested in rotational mod- eling. ? Feedback: The dynamic properties of the ight dynamics models under biologically-modeled rate feedback was investigated, nding that haltere feed- back is su cient to stabilize the insect, but that unsteady aerodynamics intro- duce instabilities that require 2.45 times the gain originally modeled to achieve an equivalent level of stability. 146 ? Morphology: Hovering ight dynamics are found to be similar for large and small Drosophila, suggesting that individual morphological changes, which are commonly observed, do not present signi cantly di erent control demands. ? Forward Flight: This dissertation presents transfer functions for longitudinal and lateral-directional Calliphorid ies maneuvering about forward (cruise) ight. The pole structure for Calliphorid longitudinal forward fright dynamics is qualitatively similar to hovering Drosophila, with two stable real poles (one of which is a decoupled heave pole) and an unstable oscillatory pair. The kinematic perturbation ud gives rise to proverse roll/yaw coupling in forward ight, which can reduce feedback requirements. ? Reachability: Results from control theory are applied to calculate the reach- able states with a given choice of inputs. Quantifying the reachable space leads to a method of ranking inputs and allows for intelligent choice of inputs for an MAV to maximize the reachable space. 9.3 Future Work This research has led to several additional paths. Several are obvious continu- ations of the current research, such as the continued work on the numerical Navier Stokes solver to compare the remainder of the longitudinal and lateral models, and more general forward ight dynamics models. Many dipterans are rarely found hovering (Ellington, 1984b), and previous work has already demonstrated forward ight does a ect the aerodynamics (Dickson 147 and Dickinson, 2004). Changes in the dynamic properties of apping wing ight as an insect transitions to forward ight may be analogous to the changes observed as rotorcraft transition to forward ight. Kinematics have recently been collected for Drosophila in untethered forward ight and Asilidae and Tabanidae are also under investigation. Several ongoing MAV designs incorporate the controllability measures devel- oped to reduce the actuation requirements, including an at-scale robotic apping wing platform incorporating all the required degrees of freedom for controlled ight. Perhaps the most signi cant impact on MAV design is the nding that a biologi- cal form of rate feedback is appropriate for ight stabilization, which has not only guided the development of a haltere-based sensor for ight stabilization, but has suggested control readily manufactured using micro-electro-mechanical machines (MEMS) technology in tandem with at-scale robotic ight platforms. Work is underway to implement a biologically inspired controller for apping wing ight stabilization on an integrated circuit, an accomplishment that will unlock a level of aerial mobility that was previously not yet possible on such minimal processing hardware. Analytic work on the ocelli and compound eye sensitivity functions are under- way and will allow researchers to perform an analysis to uncover the con gurations which best encode ight modes and actuation results. Biologically relevant time delays and phase characteristics for each of the sensor models will be included in this formulation, allowing the frequency range fractionation observed in sensing to be represented. Once a comprehensive model for the combined compound eye/ocelli 148 system has been built, an accurate control-theoretic observability analysis may be conducted, including the e ects of the di ering frequency properties of the sensing modalities. The hypothesis that insects measure quantities that take advantage of the airframe dynamic properties can then be evaluated via determining the pointing directions that minimize measures of the stat estimate quality, such as the norm of the estimate covariance. 149 Bibliography D. L. Altshuler, W. B. Dickson, J. T. Vance, S. P. Roberts, and M. H. Dickin- son. Short-amplitude high-frequency wing strokes determine the aerodynamics of honeybee ight. PNAS, 102(50):18213{18218, 2005. S. A. Ansari, R. Zbikowski, and K. Knowles. Aerodynamic modelling of insect-like apping ight for micro air vehicles. Progress In Aerospace Sciences, 42:129{172, 2006. S. A. Ansari, K. Knowles, and R. Zbikowski. Insectlike Flapping Wings in the Hover Part 1: E ect of Wing Kinematics. AIAA Journal of Aircraft, 45(6), 2008. H. Autrum. Electrophysiological Analysis of the Visual Systems in Insects. Experi- mental Cell Research, 14 Suppl. 5:426{439, 1958. S. Avadhanula, R. J. Wood, E. Steltz, J. Yan, and R. S. Fearing. Lift Force Improve- ments for the Micromechanical Flying Insect. In IEEE/RSJ Intl. Conference on Intelligent Robots and Systems, pages 1350{1356, 2003. Aviation Engineering Directorate. Aeronautical Design Standard Performance Spec- i cation Handling Qualities Requirements for Rotorcraft. Technical Report ADS- 33E, United States Army Aviation and Missile Command, March 2000. C. N. Balint and M. H. Dickinson. Neuromuscular control of aerodynamic forces and moments in the blow y, Calliphora vicinae . Journal of Experimental Biology, 207:3813{3838, 2004. J. A. Bender and M. H. Dickinson. A comparison of visual and haltere-mediated feedback the control of body saccades in Drosophila melanogaster . Journal of Experimental Biology, 209:4597{4606, 2006a. J. A. Bender and M. H. Dickinson. Visual stimulation of saccades in magnetically tethered Drosophila. Journal of Experimental Biology, 209:3170{3182, 2006b. L. Bennett. Insect Flight: Lift and Rate of Change of Incidence. Science, 167(3915): 177{179, 1970. G. J. Berman and Z. J. Wang. Energy-minimizing kinematics in hovering insect ight. Journal of Fluid Mechanics, 582:153{168, 2007. H. Bueltho , T. Poggio, and C. Wehrhahn. 3-D Analysis of the Flight Trajectories of Flies (Drosophila melanogaster). Zeitschrift for Naturforschung, C 35(1122): 811{815, 1980. Bush, B. and Baeder, J. Force production mechanisms of a apping mav wing. In AHS International Conference on Aeromechanics, San Fransisco, CA, January 23{25 2008. 150 Bush, B., MacFarlane, K., Baeder, J., and Humbert, J.S. Development of immersed boundary code with application to mav stability analysis. In 27th Army Science Conference, Orlando, FL, November 29 { December 2 2010. J. Casas and S. J. Simpson. Insect Mechanics and Control. Academic Press, 2007. B. Cheng, S. N. Fry, Q. Huang, and X. Deng. Aerodynamic damping during rapid ight maneuvers in the fruit y drosophila. The Journal of Experimental Biology, 213:602{612, 2010. J. D. Cole. Perturbation Methods in Applied Mathematics. Blaisdell Publishing Company, Waltham MA, 1968. M. J. Corless and A. E. Frazho. Linear Systems and Control: An Operator Perspec- tive. Marcel Dekker, Inc. New York, NY, 2003. W. D. Crill, R. B. Huey, and G. W. Gilchrist. Within- and between-generation e ects of temperature on the morphology and physiology of Drosophila melanogaster. Evolution, 50(3):1205{1218, 1996. G. Dalpiaz, A. Rivola, and R. Rubini. E ectiveness and sensitivity of vibration processing techniques for local fault detection in gears. Mechanical Systems and Signal Processing, 14(3):387 { 412, 2000. ISSN 0888-3270. doi: DOI:10.1006/ mssp.1999.1294. X. Deng, L. Schenato, and S. S. Sastry. Flapping Flight for Biomimetic Robotic Insects: Part I-System Modeling. Transactions on Robotics, 22(4):789{803, 2006a. X. Deng, L. Schenato, W. C. Wu, and S. S. Sastry. Flapping Flight for Biomimetic Robotic Insects: Part I-System Modeling. Transactions on Robotics, 22(4):776{ 788, 2006b. W. Derham. Physico-theology. Journal of Experimental Biology, 1714. M. H. Dickinson. Unsteady mechanisms of force generation in aquatic and aerial locomotion. American Zoologist, 36:537{554, 1996. M. H. Dickinson and K. Gotz. Unsteady Aerodynamic Performance of Model Wings at Low Reynolds Numbers. Journal of Experimental Biology, 174:45{64, 1999. M. H. Dickinson, F. O. Lehmann, and S. P. Sane. Wing Rotation and the Aerody- namic Basis of Insect Flight. Science, 284(1954), 1999. W. Dickson, A. Straw, and M. H. Dickinson. Integrative Model of Drosophila Flight. AIAA Journal, 46(9), 2008. W. B. Dickson and M. H. Dickinson. The e ect of advance ratio on the aerodynamics of revolving wings. Journal of Experimental Biology, 207:4269{4281, 2004. 151 W. B. Dickson, P. Polidoro, M. M. Tanner, and M. H. Dickinson. A linear systems analysis of the yaw dynamics of a dynamically scaled insect model. Journal of Experimental Biology, 213(17):3047{3061, 2010. doi: 10.1242/jeb.042978. D. B. Doman, M. W. Oppenheimer, and D. O. Sigthorsson. Wingbeat Shape Mod- ulation for Flapping-Wing Micro-Air-Vehicle Control During Hover. Journal of Guidance, Control, and Dynamics, 33(3):724{739, 2010a. D. B. Doman, M. W. Oppenheimer, and D. O. Sigthorsson. Dynamics and Control of a Biomimetic Vehicle Using Biased Wingbeat Forcing Functions: Part II - Controller. In Proceedings of the 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Jan 2010b. AIAA. R. Dudley. Extraordinary Flight Performance of Orchid Bees (Apidae euglossoni) Hovering in Heliox . Journal of Experimental Biology, 198:1065{1070, 1995. G. E. Dullerud and F. Paganini. A Course in Robust Control Theory. Springer, New York, NY, 2000. C. P. Ellington. The Aerodynamics of Hovering Insect Flight IV: Aerodynamic Mechanisms. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences (1934-1990), 305(1122):79{113, 1984a. C. P. Ellington. The Aerodynamics of Hovering Insect Flight III: Kinematics. Philo- sophical Transactions of the Royal Society of London. Series B, Biological Sciences (1934-1990), 305(1122), 1984b. C. P. Ellington. The Aerodynamics of Hovering Insect Flight II: Morphological Parameters. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences (1934-1990), 305(1122), 1984c. R. A. Ennos. The Kinematics and Aerodynamics of the Free Flight of some Diptera. Journal of Experimental Biology, 142(1):49{85, 1989. I. A. Faruque and J. S. Humbert. Dipteran Insect Flight Dynamics Part 1: Lon- gitudinal Motion about Hover. Journal of Theoretical Biology, 264(2):538{552, 2010. E. I. Fontaine, F. Zabala, M. H. Dickinson, and J. W. Burdick. Wing and body mo- tion during ight initiation in Drosophila revealed by automated visual tracking. Journal of Experimental Biology, 212:1307{1323, 2009. doi: 10.1242/jeb.025379. G. Fraenkel. The function of the halteres of ies (Diptera). In Proceedings of the Zoological Society of London, volume 109, pages 69{78, 1939. J. A. Franklin. Dynamics, control, and ying qualities of VSTOL aircraft. American Institute of Aeronautics and Astronautics, 2002. Reston, VA. 152 S. Fry, R. Sayaman, and M. H. Dickinson. The Aerodynamics of Free-Flight Ma- neuvers in Drosophila. Science, 300(495), 2003. R. J. Full and M. A. R. Koehl. Drag and Lift on Running Insects. Journal of Expermental Biology, 176:89{101, 1992. R. H. G Nalbach. The halteres of the blow y Calliphora I: Kinematics and dynamics. Journal of Comparative Physiology, 173:293{300, 1993. N. Gao, H. Aono, and H. Liu. A Numerical Analysis of Dynamic Flight Stability of Hawkmoth Hovering. Journal of Biomechanical Science and Engineering, 4(1): 105{116, 2009. K. G. Gotz. Course-Control, Metabolism and Wing Interference During Ultralong Tethered Flight in Drosophila melanogaster . Journal of Experimental Biology, 128:35{46, 1987. T. L. Hedrick, B. Cheng, and X. Deng. Wingbeat Time and the Scaling of Passive Rotational Damping in Flapping Flight. Science, 324:252{255, 2009. R. K. He ey, W. F. Jewell, J. M. Lehman, and R. A. Van Winkle. A Compilation and Analysis of Helicopter Handling Qualities Data Volume One: Data Compi- lation. Technical Report CR-3144, volume 1, National Aeronautics and Space Administration, August 1979a. R. K. He ey, W. F. Jewell, J. M. Lehman, and R. A. Van Winkle. A Compi- lation and Analysis of Helicopter Handling Qualities Data Volume Two: Data Analysis. Technical Report CR-3144, volume 2, National Aeronautics and Space Administration, August 1979b. T. Hesselberg and F. O. Lehmann. Corrigendum: Animal Flight Dynamics II. Longitudinal Stability in Flapping ight. Journal of Experimental Biology, 210: 4139{4334, 2007. J. S. Humbert and A. M. Hyslop. Bio-Inspired Visuomotor Convergence. Transac- tions on Robotics, 26:121{130, 2010. W. Johnson. Helicopter Theory. Dover Publications, 1994. R. Jones and D. Cohen. Aerodynamics of Wings at High Speeds. In A. Donovan and H. Lawrence, editors, High Speed Aerodynamics and Jet Propulsion, volume 7, pages 36{48, 1957. D. Karan, J. P. Morin, B. Moreteau, and J. R. David. Body size and developmental temperature in Drosophila melanogaster : analysis of body weight reaction norm. Journal of Thermal Biology, 23(5):301{309, 1998. J. Kevorkian. The Two Variable Expansion Procedure for the Approximate Solution of Certain Nonlinear Di erential Equations. In J. B. Rosser, editor, Lectures in Applied Mathematics, Volume VII, volume 7, pages 206{275, 1966. 153 J. G. Leishman. Principles of Helicopter Aerodynamics. Cambridge University Press, 2006. G. Nalbach. The halteres of the blow y Calliphora II. Three-dimensional orga- nization of compensatory reactions to real and simulated rotations. Journal of Comparitive Physiology, 175:695{708, 1994. R. C. Nelson. Flight Stability and Automatic Control. McGraw Hill, Inc, 1989. K. Ogata. Modern control engineering. Prentice Hall, 2002. Upper Saddle River, NJ. M. W. Oppenheimer, D. B. Doman, and D. O. Sigthorsson. Dynamics and Con- trol of a Biomimetic Vehicle Using Biased Wingbeat Forcing Functions: Part I - Aerodynamic Model. In Proceedings of the 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Jan 2010. AIAA. M. F. M. Osborne. Aerodynamics of Flapping Flight with Application to Insects. Journal of Experimental Biology, 28:221{245, 1951. M. M. Parsons, H. G. Krapp, and S. B. Laughlin. Sensor fusion in identi ed visual neurons. Current Biology, (20):1{5, 2010. R. Ramamurti and W. C. Sandberg. A three-dimensional computational study of the aerodynamic mechanisms of insect ight. Journal of Experimental Biology, 205:1507{1518, 2002. R. Ramamurti and W. C. Sandberg. Computational Investigation of the three- dimensional unsteady aerodynamics of Drosophila hovering and maneuvering. Journal of Experimental Biology, 201:881{896, 2007. L. Ristroph, G. J. Berman, A. J. Bergou, Z. J. Wang, and I. Cohen. Automated hull reconstruction motion tracking (HRMT) applied to sideways maneuvers of free- ying insects. Journal of Experimental Biology, 212(9):1324{1335, 2009. doi: 10.1242/jeb.025502. L. Ristroph, A. J. Bergou, G. Ristroph, K. Coumes, G. J. Berman, J. Guckenheimer, Z. J. Wang, , and I. Cohen. Discovering the ight autostabilizer of fruit ies by inducing aerial stumbles. PNAS, 107:4820{4824, 2010. C. H. Rowell and K. G. Pearson. Ocellar Input to the Flight Motor System of the Locust: Structure and Function. Journal of Experimental Biology, 103:265{288, 1983. S. P. Sane. The aerodynamics of insect ight. Journal of Experimental Biology, 206: 4191{4208, 2003. S. P. Sane and M. H. Dickinson. The aerodynamic e ects of wing rotation and a revised quasi-steady model of apping ight. Journal of Experimental Biology, 205:1087{1096, 2002. 154 L. Schenato. Analysis and Control of Flapping Flight: From Biological to Robotic Insects. PhD thesis, EECS Department, University of California, Berkeley, 2004. A. Sherman and M. H. Dickinson. A comparison of visual and haltere-mediated equi- librium re exes in the fruit y Drosophila melanogaster . Journal of Experimental Biology, 206:295{302, 2003. A. Sherman and M. H. Dickinson. Summation of visual and mechanosensory feed- back in Drosophila ight control. Journal of Experimental Biology, 207:133{142, 2004. Shyy, W., Lian, Y., Tang, J., Viieru, D., and Liu, H. Aerodynamics of Low Reynolds Number Flyers. Cambridge University Press, New York, 2008. G. R. Spedding and T. Maxworthy. The generation of circulation and lift in a rigid two-dimensional ing. Journal of Fluid Mechanics, 165:247{272, 1986. M. Sun and Y. Xiong. Dynamic Flight Stability of a hovering bumblebee. Journal of Experimental Biology, 208:447{459, 2005. M. Sun, J. Wang, and Y. Xiong. Dynamic ight stability of hovering insects. Acta Mech Sin, 2010(26):175{190, 2010. G. K. Taylor. Mechanics and aerodynamics of insect ight control. Biological Review, 76:449{471, 2001. G. K. Taylor and H. G. Krapp. Sensory systems and ight stability: what do insects measure and why? Advanced Insect Physiology, (34):231{316, 2007. G. K. Taylor and A. L. R. Thomas. Animal Flight Dynamics II. Longitudinal Stability in Flapping ight. Journal of Experimental Biology, 214:351{370, 2003a. G. K. Taylor and A. L. R. Thomas. Corrigendum: Animal Flight Dynamics II. Longitudinal Stability in Flapping ight. Journal of Experimental Biology, 221: 671, 2003b. G. K. Taylor and A. L. R. Thomas. Dynamic ight stability in the desert locust Schistocerca gregaria. Journal of Experimental Biology, 206:2803{2829, 2003c. R. A. Thompson, M. F. Wehling, and J. Evers. Evaluation of the haltere as a biologically inspired inertial rate measurement sensor. In Guidance, Navigation, and Control Conference, 2008. M. B. Tischler. Frequency-Response Method for Rotorcraft System Identi cation: Flight Applications to BO-105 Coupled Rotor/Fuselage Dynamics. Journal of the American Helicopter Society, 37:19{44, 1992. M. B. Tischler and M. G. Cau man. Comprehensive Identi cation from FrEquency Responses, Vol 3. Technical Report TR-94-A-018, US Army Aviation and Troop Command, September 1999. 155 M. B. Tischler and R. K. Remple. Aircraft And Rotorcraft System Identi cation: Engineering Methods With Flight-test Examples. American Institute of Aeronau- tics and Astronautics, 2006. Reston, VA. J. T. Vance and J. S. Humbert. Mechanisms of gust rejection in the honey bee, Apis mellifera. Technical Report 55.6, Society for Integrative and Comparative Biology, Jan 5 2010. J. T. Vance, J. B. Williams, M. M. Elekonich, and S. P. Roberts. The e ects of age and behavioral development on honey bee (Apis mellifera) ight performance. The Journal of Experimental Biology, 212:2604{2611, 2005. J. T. Vance, I. A. Faruque, and J. S. Humbert. The e ects of di erential wing stroke amplitude and stroke o set on insect body movements during perturbed ight conditions. In 34th Annual Meeting of the American Society of Biomechanics, Providence, RI, August 18-August 21 2010. S. Vogel. Flight in Drosophila: II. Variations in Stroke Parameters and Wing Con- tour. Journal of Experimental Biology, 46:383{392, 1967. T. Weis-Fogh. Energetics of Hovering Flight in Hummingbirds and in Drosophila. Journal of Experimental Biology, 56:79{104, 1972. C. H. Wolowicz, J. S. Bowman, Jr, and W. P. Gilbert. Similitude Requirements and Scaling Relationships as Applied to Model Testing. Technical Report TR-1435, National Aeronautics and Space Administration, August 1979. R. J. Wood. The First Takeo of a Biologically Inspired At-Scale Robotic Insect. Transactions on Robotics, 24(2), 2008. J. M. Zanker. The wing beat of drosophila melanogaster. i. kinematics. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 327 (1238):pp. 1{18, 1990. ISSN 00804622. J. M. Zanker and K. G. Gotz. The wing beat of Drosophila melanogaster II: dy- namics. Phil Trans R Soc Lond, 327:247{272, 1990. 156