ABSTRACT Title of dissertation: AN INVESTIGATION OF A MAV-SCALE FLEXIBLE FLAPPING WING IN FORWARD FLIGHT: FLOW FIELD AND AIRLOADS EXPERIMENTS WITH COUPLED CFD-CSD David Benjamin Mayo Doctor of Philosophy, 2014 Dissertation directed by: Dr. Inderjit Chopra Department of Aerospace Engineering Experiments were systematically executed in conjunction with a coupled com- putational fluid dynamics (CFD) and computational structural dynamics (CSD) aeroelastic analysis for a micro air vehicle (MAV)-scale flexible flapping-wing in forward flight. 2-D time-resolved particle image velocimetry (PIV) and force mea- surements were performed in a wind tunnel at a flow speed of 3 m/s. Rigid and flexible wings were tested. The wings underwent pure flapping kinematics, held at a fixed wing-pitch angle at the root. Chordwise velocity fields were obtained at equally spaced spanwise sections along the wing (30% to 90% span) at two instants during the flap cycle (mid-downstroke and mid-upstroke) for the reference Reynolds number of 15,000. The flowfield measurements and averaged force measurements were used for the validation of the 3-D aeroelastic model. The objectives of the combined efforts were to understand the unsteady aerodynamic mechanisms and their relation to force production on a flexible wing undergoing an avian-type flapping motion. The coupled CFD-CSD analysis combined a compressible Reynolds Averaged Navier Stokes (RANS) solver (OVERTURNS) with a multi-body structural solver (MBDyn) to resolve the complex, highly vortical, three-dimensional flow on a flexible wing. The coupled CFD-CSD predicted and experimental flowfields showed comparable results. A vorticity summation approach was used to calculate the circulation of the leading edge vortex (LEV) from both the PIV measurements and the numerical simulations to further validate the CFD-CSD code. The time-averaged vertical and horizontal aerodynamic forces measured from the experiment using a miniature force balance were also used to validate the force predictions from the CFD-CSD model. Pertaining to the flow physics, the flapping motion induces large angles of attack along the wing span. This motion causes the outboard sections to stall during downstroke. However, during the upstroke, the outboard sections operate at very low or even negative angles of attack. It was seen that the dynamic twisting produced by the flexible wing helped in decreasing the effective angle of attack (nose-down) during the downstroke and increased the angle of attack (nose-up) for the upstroke. This temporal and spanwise variation of wing pitch angle affected both lift and drag, and primarily helped the wing produce positive thrust during both upstroke and downstroke, which is not possible with a rigid wing undergoing pure flap at a constant pitch angle. Both PIV and CFD-CSD studies showed that the LEV stayed attached for a longer duration of the flap cycle during downstroke for the flexible wing compared to the rigid wing, especially towards the outboard sections. Also, during the upstroke, the LEV strength for the flexible wing was significantly higher than that for the rigid wing. This is due to the effective change in angle of attack. Concerning performance for the flexible wing, the optimal CZ/CX occurs at the flap frequency 10 Hz and initial set pitch angle 12◦. INVESTIGATION OF A MAV-SCALE FLEXIBLE FLAPPING WING IN FORWARD FLIGHT: FLOW FIELD AND AIRLOADS EXPERIMENTS WITH COUPLED CFD-CSD by David B. Mayo Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2014 Advisory Committee: Dr. Inderjit Chopra, Chair/Advisor Dr. James E. Hubbard Jr. Dr. James Baeder Dr. Darryll J. Pines Dr. Amr Baz © Copyright by David Benjamin Mayo 2014 To my wife and my parents ii Acknowledgments First, I must express my heartfelt appreciation to Dr. Inderjit Chopra for affording me the opportunity to work on a very exciting and dynamic research topic and the use of the Micro Air Vehicle Laboratory and Alfred Gessow Rotorcraft Center facilities. His support and leadership helped me to grow within the group and prepared me for the many challenges that lay ahead. I would also like to thank my committee Dr. James E. Hubbard Jr., Dr. Darryll J. Pines, and Dr. James Baeder for their mentorship during the doctoral process. There is no finer group of academic mentors and experts. I would like to thank the aerospace engineering staff for their help assistance and making me feel like I was part of a family. Special thanks to student and friends; Dr. Moble Benedict and Dr. V.T Nagaraj for assisting with experiments, giving guidance, and providing expertise. James Lankford has been a great friend and colleague academic. I thank him especially for the motivating conversations, friendship, constructive criticism, and always maintaining our “Team Flapping-Wing” spirit even when things were not going so well. Additional thanks goes to Dr. Joseph Milluzzo, Dr. Anish Sidney, Dr. Bharath Govindarajan, and Joseph Ramsey for their belief in me when others did not. I would like to acknowledge my colleagues in the lab for their meticulous attention to detail, insight, and all the help they provided throughout the experimentation and analysis process. Finally, I must thank my wife/best friend, Talitha Hampton-Mayo for bringing out the best in me. Her love, encouragement, and unrelenting support served as a iii source of encouragement through this challenging process. I would like to thank my parents (Jessie and Conchita Mayo) and my 6 siblings (especially my oldest brother Dr. Jessie B. Mayo Jr. who passed away on April 16, 2013). My hope is that they can be inspired by this milestone in my life knowing how much they have inspired me to rise above our humble beginnings. My family’s unconditional love anchored my self-confidence and enabled me to accomplish all that I have. Thank you God for giving me an amazing family, friends, professors, and opening doors that I never knew existed. You saw fit to bring me through the violence and destruction of war and combat as a US Marine (mentally and physically in tact). I pray that this milestone in my life can be an example to others. Once a Marine, Always a Marine... Semper Fidelis. iv Table of Contents List of Tables viii List of Figures ix List of Abbreviations xiv 1 Introduction 1 1.1 Motivation for Micro Air Vehicle (MAV) Development . . . . . . . . . 1 1.1.1 MAV Design Specifications . . . . . . . . . . . . . . . . . . . . 2 1.2 Fixed-Wing MAVs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Rotary-wing MAVs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Flapping Wings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4.1 Bio-inspired Flapping MAVs . . . . . . . . . . . . . . . . . . . 19 1.4.1.1 Entomopters . . . . . . . . . . . . . . . . . . . . . . 20 1.4.1.2 Ornithopters . . . . . . . . . . . . . . . . . . . . . . 25 1.5 Challenges to the Development of MAVs . . . . . . . . . . . . . . . . 28 1.5.1 Vehicle Control/Actuation . . . . . . . . . . . . . . . . . . . . 29 1.5.2 Energy Sources . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.5.3 Propulsion and Platform Integration . . . . . . . . . . . . . . 32 1.5.4 Low Reynolds Number Effects . . . . . . . . . . . . . . . . . . 34 1.6 Unsteady Lift Enhancement Mechanisms . . . . . . . . . . . . . . . . 42 1.6.1 The Leading Edge Vortex . . . . . . . . . . . . . . . . . . . . 43 1.6.2 Other Lift Augmenting Mechanisms . . . . . . . . . . . . . . . 48 1.7 Review of Flapping Wing Experimental Studies . . . . . . . . . . . . 52 1.8 Review of Flapping Wing Numerical Studies . . . . . . . . . . . . . . 59 1.9 Need for New Experimental Data . . . . . . . . . . . . . . . . . . . . 64 1.9.1 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . 67 1.10 Research Objective and Approach . . . . . . . . . . . . . . . . . . . . 68 1.11 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . 69 v 2 Methodology: Numerical Approach and Experimental Setup 71 2.1 Numerical Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 2.1.1 Description of Flow Solver . . . . . . . . . . . . . . . . . . . . 72 2.1.2 Description of Grid System and Simulation . . . . . . . . . . . 74 2.1.3 Grid Sensitivity Study . . . . . . . . . . . . . . . . . . . . . . 75 2.1.4 Structural Model and Coupling . . . . . . . . . . . . . . . . . 76 2.2 Experimentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 2.2.1 Flapping Mechanism and Wing Design . . . . . . . . . . . . . 79 2.2.2 Particle Image Velocimetry Setup . . . . . . . . . . . . . . . . 84 2.2.3 Wind Tunnel Specifications . . . . . . . . . . . . . . . . . . . 89 2.2.4 PIV Experiment Challenges . . . . . . . . . . . . . . . . . . . 92 2.2.5 Experimental Uncertainty . . . . . . . . . . . . . . . . . . . . 93 2.3 Structural Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 2.3.1 Wing Covering Material Tensile Tests . . . . . . . . . . . . . . 101 2.3.2 Wing Frequency Verification Tests . . . . . . . . . . . . . . . 106 2.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3 Comparison of Predicted and Experimental Results 109 3.1 Comparison of Structural Dynamic Solver to Experiments . . . . . . 109 3.1.1 Flexible Wing Static Loading Comparison . . . . . . . . . . . 109 3.1.2 Flexible Wing Dynamic Tests . . . . . . . . . . . . . . . . . . 110 3.1.2.1 Comparison of Flexible Wing Frame Fundamental Frequencies . . . . . . . . . . . . . . . . . . . . . . . 110 3.1.2.2 Comparison of Flexible Wing Frame Flapping Motion115 3.1.3 Comparison of Aerodynamic Solver (CFD-CSD) and PIV Results118 3.1.3.1 Rigid and Flexible Wing Flowfield Comparison . . . 118 3.2 Comparison of Sectional LEV Strength . . . . . . . . . . . . . . . . . 138 3.2.1 Rigid Wing LEV Circulation (PIV and CFD) . . . . . . . . . 138 3.2.2 Flexible Wing LEV Circulation (PIV and CFD-CSD ) . . . . 143 3.3 Airloads Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 3.3.1 Rigid Wing Comparison of Aerodynamic Forces from Flowfield 146 3.3.2 Comparison of Predicted and Averaged Direct Force Measure- ments for the Rigid and Flexible Wing . . . . . . . . . . . . . 149 3.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 4 Rigid and Flexible Wing Comparison (Aerodynamic Force Production) 153 4.1 Understanding Lift and Thrust Production for a Rigid Flapping Wing 154 4.1.1 Rigid Wing Force History . . . . . . . . . . . . . . . . . . . . 154 4.1.2 Analysis of the 3-D Rigid Wing FlowField over a Flap Cycle (Q-Criterion) . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 4.1.3 Rigid Wing Analysis of the 3-D Flowfield Over a Flap Cycle (Spanwise Vorticity Contours) . . . . . . . . . . . . . . . . . . 160 4.1.4 Rigid Wing 2-D Analysis of Wing (Sectional Vertical and Horizontal Force Coefficients along the Span) . . . . . . . . . 167 vi 4.1.5 2-D Analysis of Rigid Wing (Sectional Lift and Drag Coeffi- cients along the Span) . . . . . . . . . . . . . . . . . . . . . . 173 4.1.6 Section Summary . . . . . . . . . . . . . . . . . . . . . . . . . 177 4.2 Rigid and Flexible Wing Comparison . . . . . . . . . . . . . . . . . . 178 4.2.1 Aerodynamic Force Production . . . . . . . . . . . . . . . . . 178 4.2.2 3-D Analysis of Wing (Aerodynamic Force Production over a Flap Cycle) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 4.2.3 3-D Rigid and Flexible Wing Flowfield Comparison (Q-Criterion)185 4.2.4 3-D Rigid and Flexible Wing Flowfield Comparison (Spanwise Vorticity Contours) . . . . . . . . . . . . . . . . . . . . . . . . 187 4.2.4.1 2-D Analysis of Wing (Sectional CZ and CX along the Wing Span)) . . . . . . . . . . . . . . . . . . . . 193 4.2.4.2 2-D Analysis of Wing (Sectional Cl and Cd along the Wing Span) . . . . . . . . . . . . . . . . . . . . . . . 197 4.2.5 Section Summary . . . . . . . . . . . . . . . . . . . . . . . . . 202 5 Concluding Remarks 204 5.1 Conclusions from Rigid Wing Analysis . . . . . . . . . . . . . . . . . 205 5.2 Conclusions from the Comparison of Rigid and Flexible Wings . . . . 208 5.3 Major Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 5.4 Suggestions for Future Work . . . . . . . . . . . . . . . . . . . . . . . 213 5.5 Acknowledgment of Sponsors . . . . . . . . . . . . . . . . . . . . . . 215 Bibliography 216 vii List of Tables 1.1 MAVs described by DARPA [1]. . . . . . . . . . . . . . . . . . . . . . 4 1.2 Typical parameters for wing kinematic for large hovering insects [2]. . 16 1.3 Characteristics of existing entomopter designs [3]. . . . . . . . . . . . 25 2.1 Time-resolved PIV test matrix. . . . . . . . . . . . . . . . . . . . . . 93 2.2 Tensile test results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 2.3 Flexible wing properties. . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.1 Comparison of flexible chordwise member natural frequencies. . . . . 115 4.1 Chart of averaged CZ and CX for flap frequency and pitch angle (rigid wing) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 4.2 Chart of UP, UT, V and φ along the span (rigid wing). . . . . . . . . 168 4.3 Chart of CZ and CX for flap frequency and pitch angle (flexible and rigid (in parenthesis) wing respectively). . . . . . . . . . . . . . . . . 182 4.4 Chart of UP, UT, V and φ along the span (flexible wing). . . . . . . . 194 viii List of Figures 1.1 Defined mission scenarios for MAVs [4,5]. . . . . . . . . . . . . . . . . 3 1.2 Fixed-wing MAVs (AeroVironment) [6]. . . . . . . . . . . . . . . . . . 7 1.3 Ducted single rotor, single rotor (with vanes), and coaxial MAVs. . . 9 1.4 MAV quad rotors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5 Unconventional rotating-wing MAVs (University of Maryland). . . . . 13 1.6 a) Typical hovering flight wing motion. b) Wing path for “insect” type flyers: 1) Pronation; 2) Downstroke; 3) Supination; 4) Upstroke. Modified from original CTA-MAST Program Report, used with per- mission. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.7 Schematic diagram of coordinate systems and wing kinematics for angles corresponding to Table 1.2 [7]. . . . . . . . . . . . . . . . . . . 16 1.8 Wing kinematics for various “avian” type flyers (a) albatross, fast gate; (b) pigeon, slow gate; (c) horseshoe bat, fast gate (d) horseshoe bat, slow gate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.9 Perspective views of avian wing movements during normal forward flight [8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.10 a) Upstroke and b) downstroke forces on a flapping wing in forward flight [9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.11 University of Toronto Mentor Entomopter. [10] . . . . . . . . . . . . . 20 1.12 Thrust Augmented Entomopter [11]. . . . . . . . . . . . . . . . . . . 21 1.13 Light-weight MAVs: Concept drawing of the Micromechanical Flying Insect (MFI), the prototype Microbot, and 3-D printed hover prototype. 23 1.14 Flapping-wing MAVs developed at the University of Maryland, a) Daedal and b) “64 Gram Flapper”. . . . . . . . . . . . . . . . . . . . 25 1.15 AeroVironment’s Nano Hummingbird [12]. . . . . . . . . . . . . . . . 26 1.16 Commercial Kinkade ParkHawk onithopter on UMD Morpheus Lab test mount [13]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.17 MAV-scale ornithopters . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.18 Energy storage systems energy densities [14]. . . . . . . . . . . . . . . 32 1.19 Comparison of efficiency parameters for hover-capable MAV-scale fliers [15]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 ix 1.20 Airfoil performance as a function of Reynolds number [16]. . . . . . . 40 1.21 Laminar separation bubble [17]. . . . . . . . . . . . . . . . . . . . . . 41 1.22 3-D conical LEV structure: a) LEV cross-section; b) three-dimensional structure and spanwise flow (note that le is the leading edge, te is the trailing edge, dss is the dividing stream surface) [18]. . . . . . . . . . 44 1.23 Airloads and flow structures on a 2-D airfoil during the dynamic stall process [19]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 1.24 Illustration of wing with flow restricted by baffles, and vorticity with overlaid velocity field (Note: Blue denotes a vortex) [20]. . . . . . . . 47 1.25 Illustration of the clap-and-fling/Weis-Fogh mechanism [21]: A) at end of half-stroke wings moving toward each other; B) wings continue approach to each other; C) fluid is moved out by wings; D) wings rotation about trailing edge; E) fluid pulled into the void left by the wings; F) LEV formation begins as wings move apart. . . . . . . . . . 49 1.26 Illustration of the wake capture mechanism [21]: (A) After translation, wing generates vorticity at both the leading and trailing edges, (B) Vortices induce a strong velocity field in the intervening region, (C/D) wing comes to a halt and reverses stroke (D/E), wing translates in the opposite directions [21]. . . . . . . . . . . . . . . . . . . . . . . . 51 1.27 Robofly oil tank test setup [22]. . . . . . . . . . . . . . . . . . . . . . 55 1.28 Passive pitch, bi-stable single-degree-of-freedom flapping mechanism [23]. 56 2.1 Wing mesh embedded in background Cartesian mesh. . . . . . . . . . 75 2.2 Variation of predicted aerodynamic force coefficients for different mesh resolutions (k = 0.32, θ = 0◦). . . . . . . . . . . . . . . . . . . . . . . 77 2.3 Four-bar linkages and flapping mechanism. . . . . . . . . . . . . . . . 82 2.4 a) Rigid and b) flexible wing construction. . . . . . . . . . . . . . . . 83 2.5 Schematic depicting the general setup of a 2-D PIV system and image processing methodology. . . . . . . . . . . . . . . . . . . . . . . . . . 86 2.6 Schematic PIV experiment setup in wing tunnel test section. . . . . . 87 2.7 Actual PIV experiment setup in wing tunnel test section. . . . . . . . 88 2.8 Schematic of low speed wind tunnel. . . . . . . . . . . . . . . . . . . 90 2.9 Variation in wind tunnel turbulence with wind tunnel speed. . . . . . 91 2.10 Typical test setup for motion tracking system. . . . . . . . . . . . . . 96 2.11 Schematic of VICON test for isolated bending and torsion testing. . . 98 2.12 Large deformation of a beam under point load at free tip [24]. . . . . 99 2.13 Torsion test schematic. . . . . . . . . . . . . . . . . . . . . . . . . . . 102 2.14 Image of wing membrane material in during tensile test. . . . . . . . 104 2.15 Tensile tests for wing covering. . . . . . . . . . . . . . . . . . . . . . . 104 2.16 VICON motion tracking “shaker” test for modal analysis. . . . . . . . 107 3.1 Left side, images of static wing loading tests with VICON markers and Right side, load deflection curves for chordwise ribs. . . . . . . . 111 3.2 Raw signal of tip (TE) displacements and Fast Fourier Transform from motion tracking experiment. . . . . . . . . . . . . . . . . . . . . 113 x 3.3 Mode shapes for wing frame structural members. . . . . . . . . . . . 114 3.4 Comparison experiment and predicted (MBDyn) wing frame tip dis- placements for 6 Hz flapping. . . . . . . . . . . . . . . . . . . . . . . 117 3.5 Perspective view schematic of the rotating wing setup. . . . . . . . . 119 3.6 Comparison of representative chordwise PIV and CFD vorticity con- tours at t = T/4, (40◦ flapping amplitude, θ = 24◦, 4 Hz flapping frequency, V∞ = 3 m/s, k = 0.32). . . . . . . . . . . . . . . . . . . . . 120 3.7 Comparison of representative chordwise PIV and CFD vorticity con- tours at t = T/4, (40◦ flapping amplitude, θ = 12◦, 6 Hz flapping frequency, V∞ = 3 m/s, k = 0.48). . . . . . . . . . . . . . . . . . . . . 122 3.8 Comparison of representative chordwise PIV and CFD vorticity con- tours at t = T/4, (40◦ flapping amplitude, θ0 = 12◦, 4 Hz flapping frequency, V∞ = 3 m/s, k = 0.32). . . . . . . . . . . . . . . . . . . . . 123 3.9 Comparison of representative chordwise PIV and CFD vorticity con- tours at t = T/4, (40◦ flapping amplitude, θ0 = 24◦, 6 Hz flapping frequency, V∞ = 3 m/s, k = 0.64). . . . . . . . . . . . . . . . . . . . . 124 3.10 Comparison of representative chordwise PIV and CFD vorticity con- tours at t = T/4, (40◦ flapping amplitude, θ0 = 0◦, 8 Hz flapping frequency, V∞ = 3 m/s, k = 0.64). . . . . . . . . . . . . . . . . . . . . 125 3.11 Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, rigid wing ,θ = 0◦ at 4 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000).129 3.12 Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, rigid wing, θ = 12◦ at 6 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000).131 3.13 Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, rigid wing, θ = 24◦ at 4 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000).133 3.14 Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, flexible wing, θ0 = 12◦ at 4 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 3.15 Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, θ0 = 24◦ at 6 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). . . . . . 135 3.16 Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, flexible wing, θ = 0◦ at 8 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000).137 3.17 Distribution of averaged LEV circulation, Γv along the span, rigid wing, Ω =4 Hz, (40◦ flap amplitude, V∞ = 3 m/s) with error bars showing standard deviation. . . . . . . . . . . . . . . . . . . . . . . . 139 3.18 Distribution of averaged LEV circulation, Γv for a variation in flap frequency, rigid wing, y/b = 0.5 (40◦ flap amplitude, V∞ = 3 m/s) at 50% span. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 xi 3.19 Distribution of averaged LEV circulation, Γv along the span, flexible wing, Ω =4 Hz, (40◦ flap amplitude, V∞ = 3 m/s). . . . . . . . . . . . 144 3.20 Distribution of averaged LEV circulation, Γv for a variation in flap frequency, flexible wing, y/b = 0.5 (40◦ flap amplitude, V∞ = 3 m/s) at 50% span. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 3.21 Schematic explaining CZ calculation. . . . . . . . . . . . . . . . . . . 148 3.22 Variation of vertical force coefficient (Cz) from rigid wing averaged PIV and instantaneous CFD flowfield data (50% span location, 4 Hz flapping). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 3.23 Rigid and flexible wing comparison of averaged experiment and pre- dicted vertical (FZ) and horizontal forces (FX) for various flap frequencies.151 4.1 Variation of total vertical force coefficient (CZ) and total horizontal force coefficient (CX) from CFD data (rigid wing). . . . . . . . . . . . 156 4.2 Iso-surfaces of Q-Criterion colored with vorticity magnitude contour at various instances of the flap cycle (rigid wing). . . . . . . . . . . . 159 4.3 Vorticity magnitude along the span at various instances of flap cycle (rigid wing). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 4.4 Vorticity contours with overlaid velocity along the span: left side, mid-downstroke and right side, mid-upstroke (rigid wing), 4 Hz, θ = 24◦.166 4.5 X-Vorticity magnitude along the span from CFD data (rigid wing), 4 Hz, θ = 24◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 4.6 Rigid wing sectional Cz and Cx, at mid-downstroke and mid-upstroke, normalized by V∞. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 4.7 Rigid wing Cz and Cx, at mid-downstroke and mid-upstroke, normal- ized by the resultant V. . . . . . . . . . . . . . . . . . . . . . . . . . 172 4.8 Schematics of force and velocity vectors for the airfoil cross section. . 174 4.9 Lift coefficient, Cl, and drag coefficient, Cd, along the span, at mid- downstroke and mid-upstroke. . . . . . . . . . . . . . . . . . . . . . . 175 4.10 Rotation of the lift and drag vectors toward the wing tip. . . . . . . . 177 4.11 Rigid and flexible wing comparison of predicted vertical (FZ) and horizontal force (FX) for single flap cycle with flap frequency of 6 Hz and θ0 = 12◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 4.12 Variation of total vertical force coefficient (CZ) and total horizontal force coefficient (CX) from coupled CFD-CSD data (flexible wing). . . 183 4.13 Rigid and flexible wing comparison of iso-surfaces of Q-Criterion colored with vorticity magnitude at various instances of the flap cycle. 188 4.14 Superimposed rigid and flexible showing deformations, θ0 = 12◦, Ω =6 Hz, (40◦ flap amplitude, V∞ = 3 m/s). . . . . . . . . . . . . . . . . . 189 4.15 Chordwise vorticity contour with overlaid velocity field at 30%, 50%, 70%, and 90% span locations at mid-downstroke (flexible and rigid wing). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 4.16 Chordwise vorticity contour with overlaid velocity field at 30%, 50%, 70%, and 90% span locations at mid-upstroke (flexible and rigid wing).192 xii 4.17 Sectional Cz and Cx, at mid-downstroke and mid-upstroke, 6 Hz, θ0 = 12◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 4.18 Variation of geometric pitch angle along the span during the mid- upstroke and mid-downstroke at an initial pitch angle of 12◦ for 4, 6, 8, 10 Hz flap frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . 199 4.19 Variation in geometric pitch angle at the mid-downstroke and mid- upstroke (6 Hz). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 4.20 Lift coefficient, Cl, and drag coefficient, Cd, along the span at mid- downstroke and mid-upstroke. . . . . . . . . . . . . . . . . . . . . . . 201 xiii Nomenclature A Flap amplitude Ae Effective area AR Aspect ratio b Wing span c Wing chord Cd Sectional drag coefficient Cl Sectional lift coefficient CX Horizontal force coefficient CZ Vertical force coefficient EI Bending stiffness FM Figure of merit FX Propulsive Force FZ Vertical Force GJ Torsional rigidity J Advance ratio k Reduced frequency, k = Ωc2V∞ M Image magnification P Pressure PL Power loading r Radial dimension rc Vortex core radius Re Chord Reynolds number = ρV c/µ St Struhaul number t Time, s T Period of flap cycle u, v, w Flow velocity in the x, y, z directions UP Parallel (to flap plane) component of resultant velocity UT Normal (to flap plane) component of resultant velocity V Resultant velocity experienced by the wing Vr Radial velocity Vθ Swirl velocity V∞ Free stream velocity x, y, z Cartesian coordinate system xiv Greek Symbols α Angle of attack δ change in specimen length  Normal strain Γv Vortex circulation µ Dynamic viscosity ω Vorticity, rad/s Ω Frequency of flapping, rad/s φ Flap induced angle Ψ Phase angle in the flapping cycle σ Normal stress ρ Density θ Pitch angle xv List of Abbreviations CCD Charged Coupled Device CFD Computational Fluid Dynamics CMOS Complementary Metal Oxide Semiconductor CSD Computational Structural Dynamics c.g. center of gravity DAQ Data Acquisition DARPA Defense Advanced Research Projects Agency DPIV Digital Particle Image Velocimetry FFT Fast Fourier Transform fps Frames-Per-Second FSI Fluid Structure Interactions FV Flow Visualization GTOW Gross Takeoff Weight IVP Initial Value Problem LEV Leading Edge Vortex MAV Micro Air Vehicle MBDyn Multibody Dynamics MLS Moving Least Squares MSD Multibody System Dynamics N-S Navier Stokes RC Radio Controlled Re Reynolds Number URANS Unsteady Reynolds-Averaged Navier-Stokes xvi Chapter 1: Introduction 1.1 Motivation for Micro Air Vehicle (MAV) Development Because of the continued increased use and efficacy of unmanned aircraft, research and development of micro air vehicles (MAVs) continues for a variety of military and civilian applications. In general, MAVs are autonomous flying vehicles with a size constraint (typically palm-size). They are small light-weight systems, which offer the attributes of portability, stealth, and enhanced flight agility. Their relative size also provides the benefit of a vehicle with comparatively low production and operational costs compared to larger scale unmanned aerial vehicles [25–27]. These characteristics make MAVs an ideal platform for reconnaissance and surveillance missions in confined and hazardous areas. MAVs have both a military and civilian role. This is because hover-capable/low- speed flight vehicles with gust tolerances can perch and observe, offer greater situa- tional awareness for military covert scouting/reconnaissance, and search and rescue missions inside buildings, caves or tunnels. As it pertains to civil applications, MAVs can be useful for missions such as investigating a chemically dangerous area or a radioactive environment, border security, law enforcement, disaster response, 1 homeland defense, pipeline patrol, etc. Figure 1.1 illustrates the mission scenarios defined by the Defense Advanced Research Projects Agency (DARPA) and the Army Research Lab [5]. The need for such vehicles that can handle these scenarios has been around for decades, and today there is even more of a need for this capability that MAVs can provide. 1.1.1 MAV Design Specifications Most major MAV research started in 1997 when DARPA requested the devel- opment of what was then a new type of unmanned vehicle with a size constraint [4]. The program’s intention was to spark innovation in creating a new type of small vehicle. Initially, a 15 cm size constraint, an endurance of about an hour, and a gross takeoff weight (GTOW) of less than 100 grams were defined for such a vehicle. However, DARPA initiatives expanded to require a vehicle which could fit into a cylinder 7.5 cm in diameter and up to 7.5 cm long [28]. Table 1.1 lists some other MAV design specifications. Furthermore, DARPA defined a second, lighter class of vehicles called “Nano-Air” (NAVs), which should have a GTOW of less than 20 grams. While the classification of “MAVs” has grown today to encompass vehicles both smaller and larger than 15 cm, enormous scientific interest continues to drive the development of small, autonomous MAVs. The size constraint presents several challenges in the areas of low Reynolds number aerodynamics, sensors and actuators, propulsive systems design, energy stor- age, and system integration. Although these challenges are inherent to MAV design 2 (a) 1997 DARPA illustration of MAV mo- bile immersion sensor mission scenario. (b) 1997 DARPA illustration of Urban Op- erations Missions scenario for MAVs. (c) 1997 DARPA illustration of ”Over- The-Hill Reconnaissance” mission sce- nario for MAVs. (d) 2014 CTA-MAST illustration of perimeter defense and building search mission scenarios. Figure 1.1: Defined mission scenarios for MAVs [4,5]. 3 Table 1.1: MAVs described by DARPA [1]. Criteria Specification Requirements Vehicle dimension Weight, g (GTOW) 100 Size, m (Maximum dimension) <15.24 Performance parameters Speed, m/s (Maximum flight speed) 15 Endurance, min (Loiter time on station) 60 Range, km (Operational range) 1–10 Altitude, m (Operational ceiling) <150 Weapons/Sensor package Payload, g, (Mission dependent) 20 Budget Cost, $ (Maximum cost) 1,500 and practical solutions will be the foundation of building successful MAV concepts, research results reported in this dissertation are focused mainly on aerodynamics issues. Currently, the two main types of MAVs being developed are rotating-wing and fixed-wing. Though these two types of types have been the most commonly examined concepts in the quest for successful MAVs, small bio-inspired platforms such as flapping-wings are also an option for meeting the demands of the aforementioned mission scenarios (e.g., highly agile/maneuverable, operation in confined, obstacle- rich urban environments, compact, and hover-capable). These types of MAVs are described in the next section. 1.2 Fixed-Wing MAVs MAVs of the fixed-wing variety are more technologically developed than other MAVs. These vehicles have the advantage of being less complex than other MAVs as there are often fewer moving parts associated with the design. These platforms also have high cruise speeds, more range, and very good stability characteristics. Also, like 4 their large-scale high endurance counterparts (e.g., modern bombers or commercial aircraft), they have the advantage of being relatively efficient. These types of vehicles are often deployed on long endurance missions as they have characteristics that are best utilized for a perimeter defense or reconnaissance type mission. However, research continues to further improve these vehicles. Numerous experimental and numerical studies of low-Reynolds number airfoil and wing designs have been performed [25,29–31]. Specifically, Mueller has deter- mined the aerodynamic characteristics of many MAV-scale low-aspect ratio wings varying wing planform, aspect ratio, and Reynolds number in a low-speed wind tunnel [32]. These wing model tests have proven very useful for MAV development in that they have shown the effects of a range of wing parameters at the MAV-scale. In the area of MAV platform testing, the University of Arizona built four MAV wind-tunnel models to test the effects of varying wing camber [33]. It was observed that at the MAV-scale, a low camber wing yields an optimal lift-to-drag ratio for high-speed, and efficient flight, and highly cambered wings gave the best low-speed performance because near stall the wings experienced only mild pitching moments while maintaining high lift-to-drag ratios. In the area of flexible wings, a very popular field of research, the University of Florida used wind tunnel testing to compare flexible wings to a rigid baseline wing [34]. Airloads, 3-D visual image correlation for deformation measurements, and flow visualization were examined. It was observed that pitch stability could be improved by using wing perimeter reinforced flexible wings. This resulted in an increase in maximum lift coefficient. Both the University of Arizona and and 5 the University of Florida conducted tests on MAVs they developed. Current re- search innovation pertaining to fixed-wing MAVs has been in the area of active morphing. This is a technique where the wing lift distribution is controlled using spatially distributed actuators, conformal sensors, and aeroelastic deformations of light-weight/flexible/wing membranes to improve performance [35–38]. There have been several successful fixed-wing MAV prototypes developed over the past two decades. Figure 1.2 shows three images from AeroVironment’s fixed-wing MAV program. In 1998, they produced the 80 g Black Widow (see Fig. 1.2(a)). This vehicle has a 6 inch wing span with an endurance of 30 minutes. Later, in 2002, an endurance record of 107 minutes was established by the WASP [39]. WASP has a wing span of 33 cm and weighed 170 grams (see Fig. 1.2(b)). The latest prototype (2006) is called WASP Block III, selected by the U.S. Air Force’s Battlefield Air Targeting Micro Air Vehicle (BATMAV) program. Each successive aircraft experiences an increase in aspect ratio to reduce drag and an increase in battery size. WASP III (see Fig. 1.2(c)), has a wingspan of 72 cm and a weight of 430 grams [6]. Other major fixed-wing MAV prototypes include the University of Florida’s 4.5 inch “bendable-wing” MAV [34, 40, 41], the Black Kite (National Aerospace Laboratories, India) [42], and NRL’s MITE [43]. There are some design challenges and disadvantages to the use of fixed-wing MAVs. These platforms use low aspect ratio (AR) wings designed to maximize wing area for a constrained wingspan. With this, there is a considerable increase in the amount of induced drag on the vehicle. This is also the case with large-scale airplane wings but the effects are more severe for fixed-wing MAVs. Typically, MAV 6 (a) Black Widow (AeroVironment 1998) [44]. (b) WASP Block I (AeroVironment 2002) [45]. (c) WASP Block III (AeroVironment 2006). Figure 1.2: Fixed-wing MAVs (AeroVironment) [6]. 7 fixed-wing aspect ratios are in the range of 1 < AR < 2, while large-scale endurance UAVs have an aspect ratio of above 20. Another major disadvantage is in the inability of these aircraft to fly at very low speeds or to hover. This is because they must maintain a high forward velocity so that their wings can generate lift. Also, there is a large turning radius for these vehicles. Therefore, even with the many innovations directed at improving fixed-wing aircraft performance and endurance, these vehicles are not well suited for missions in confined areas [46]. 1.3 Rotary-wing MAVs Currently, rotating-wing MAVs represent the majority of the vehicles being developed at the MAV scale. This is because rotating-wing MAVs have the ability to hover. Most micro-rotorcraft use their rotors to provide both lift and thrust, though a ring-wing shroud can augment lift in hover. To reduce the complexity of the rotor system, many rotary-wing MAV designs use fixed-pitch rotors and use control tabs in the rotor downwash or incorporate multiple rotors to control the vehicle. The main types of rotary-wing MAVs are the single main rotor, coaxial, and quadrotor systems. At the MAV scale, single main rotor designs with a mechanism for anti-torque or a coaxial configuration are often used for low-speed flight (see Figs. 1.3). A single main rotor configuration would require a tail rotor as a separate device for anti-torque. To eliminate the need for a tail rotor and to reduce the overall size of the single rotor design, vanes in the downwash can be used to counteract the torque. Many single 8 (a) University of Maryland, “TiShrov” [47]. (b) University of Kansas, shrouded single rotor “XQ-138” [48]. (c) University of Maryland, Coaxial “Mi- cor” [48,49]. Figure 1.3: Ducted single rotor, single rotor (with vanes), and coaxial MAVs. 9 main rotor concepts are basically scaled-down conventional helicopters and suffer from comparatively lower efficiency ratings. Ducts and shrouds have been utilized to improve aerodynamic efficiency for the coaxial design. The duct acts to accelerate the inlet airflow and reduces tip loss effects to provide increased thrust. At the MAV scale, the thrust gains need to exceed the duct weight. Ducts also add profile drag to the vehicle and large pitching moment in forward flight. In general, because of their inherent compactness, coaxial designs appear more attractive than vehicles with a single rotor, but neither are particularly able to handle gusty environments. The advantage of the single and coaxial designs is in the general familiarity in construction and control of the vehicle. Another widely researched rotary-wing concept is the quad rotor (see Fig. 1.4). These vehicles have four rotors and offer a wide range of advantages such as good flight stability and agility, large operational envelopes, and improved thrust-to-weight ratios. However, disadvantages are in the area of size. Vehicle weight can be an issue because it can be difficult to scale down the motors and the four rotors consume more power than a vehicle with fewer rotors. There are several unconventional rotating-wing concepts being researched (e.g. flapping rotors, powered Samaras, and cyclocopters). Referred to as Flotors (flapping rotors), ornicopters, or rotopters are bio-inspired designs that combine the motion of vertically flapped wings with the rotary-wing motion of helicopters. Figure 1.5(a) shows the micro rotor experimental test stand used to demonstrate the effects of high frequency blade flapping on average rotor forces at the University of Maryland. This test rig is capable of independent flapping and rotating motion [50]. Blade 10 (a) University of Maryland/Daedalus Flight Sys- tems Micro Quad. (b) Stanford University Mesicopter. Figure 1.4: MAV quad rotors. rotation and superimposed blade flapping showed improvements in both figure of merit (FM), ratio of ideal power required to hover to actual power required, and maximum thrust. However, one disadvantage is in increased design complexity. Also, there are periodic inertial forces that are directly proportional to the mass of the blades and improvements in efficiency often do not justify the increase in weight and complexity of the system. Another bio-inspired design, the Samara, is based on the winged seed of the maple tree. When these seeds fall to the ground, they generate lift by auto-rotating. An MAV concept was designed based on the same principle. To hover, a motor driven propeller is added to maintain the rotation of the vehicle. Figure 1.5(b) shows the University of Maryland’s Samara MAV design [51]. This vehicle is controlled by manipulating the propeller RPM/throttle and the collective pitch. The disadvantage of this type of MAV is that it needs a launch space, the payload (cameras, weapons, or sensors) should not be sensitive to the rotation of the vehicle, and stabilizing the 11 vehicle can be challenging. In a cycloidal rotor system (Ref. [52]), the span of the blades run parallel to the axis of rotation of the rotor, i.e., a form of paddle-wheel type of concept (see Fig. 1.5(c)). The effective disk area in this case can be defined as the product of rotor span and rotor diameter. The pitch angle of each blade on a cyclorotor is varied periodically such that the angle of attack is such that the blades create lift at both the top and bottom parts of the rotational cycle. Varying in configuration and size/weight, the cyclorotor concepts developed at the University of Maryland have achieved efficiency values within the range of 0.42 < FM < 0.65 for range of effective disk loadings, whereas conventional rotor systems operated at MAV scales (when operated with the same or similar values of disk loading) achieve a somewhat lower figure of merit of between 0.4 and 0.5 [15]. The optimized MAV-scale cycloidal rotor could achieve a higher power loading than conventional rotor in hover. For comparison, it cannot be overlooked that rotating-wing MAVs have lower aerodynamic efficiency ratings and in some cases be more difficult to control compared to their fixed-wing counterparts, especially in gusty environments. The trade-off is that fixed-wing vehicles cannot hover. With these limitations, it is important to consider flight vehicles, which can hover, handle gusty environments, and fly at very low-speeds. Inspired by natural flyers, flapping-wing MAVs may present one solution. The next section discusses some of the recent advancements in flapping wing designs. 12 (a) Flapping rotor test rig, “Flotor” [50]. (b) Robotic Samara [51]. (c) Cyclocopter MAVs [52]. Figure 1.5: Unconventional rotating-wing MAVs (University of Maryland). 13 1.4 Flapping Wings Bio-mimetic flapping wings have received much interest because of perceived high aerodynamic efficiencies and gust tolerance at the MAV scale [53]. In nature, the four biological groups that have evolved for flapping flight are pterasaurs (vertebrates), and insects (invertebrates). Vertebrate (bird and bat) wings are evolved through the modification of an existing limb so they have muscles for control. Insect wings consist of thin non-living membranes and rods, and they are controlled by pulling on those rods. These wings are passive structures controlled from the wing base. These two types of flapping motion are nearly vertical (“avian” type) and horizontal (“insect” type) with respect to the body frame. This is not a fixed physical rule because there are a few examples of birds and insects using both vertical and horizontal wing motions. Nevertheless, bird-type and insect-type MAVs are often called ornithopters and entemopters respectively. The “insect” type wing kinematics are most often observed in hovering flight as the wings move at a very high wing beat frequency at low Reynolds numbers. For hover-type wing kinematics, the stages of the wing flap are the upstroke translation, wing rotational stages, and downstroke translation (see Fig. 1.6). Here, the leading edge of the wing is denoted by a dot on the wing section, and it is the same during both half-strokes. These stages occur in the horizontal plane relative to the body frame. The wing sweeps through the air at relatively large pitch angles. The wing rotational stages occur at the end of each half-stroke and are called pronation and supination. This is where a reversal in wing orientation occurs and the upper 14 Figure 1.6: a) Typical hovering flight wing motion. b) Wing path for “insect” type flyers: 1) Pronation; 2) Downstroke; 3) Supination; 4) Upstroke. Modified from original CTA-MAST Program Report, used with permission. surface becomes the lower surface. This ensures that the wings reach positive pitch angles on both parts of the translational motion. There can be a small amount of stroke deviation, where the wing moves perpendicular to the stroke plane; however, its purpose is not well understood [7]. As seen in Fig. 1.6, both wings move symmetrically, and the motion is described using kinematic parameters like stroke amplitude, wing beat frequency, wing angle of attack, stroke plane angle, downstroke/upstroke ratio, wing tip trajectory, and timing for wing rotation. Typical wing parameters for large insects are shown in Table 1.2. Typically, for most large animals, hovering flapping stroke amplitude is 120◦ to 180◦ or greater. The geometric angle of attack of hovering insects at mid-stroke is 30◦ to 45◦ [2]. In forward flight, the wing motion is more birdlike as the stroke plane tilts forward. Figure 1.7 describes the insect wing and body movement. The body kinematics can be represented by the body angle χ and the stroke refers to angle 15 Table 1.2: Typical parameters for wing kinematic for large hovering insects [2]. Kinematic Parameter Typical Value Wing-beat Frequency, n 20− 40 Hz Stroke Amplitude, Φ 120◦ Angle of Attack, α 30◦ Body Angle, χ 50◦ → 10◦ (flight speed: 0→ max) Stroke Plane Angle, β 10◦ → 50◦ (flight speed: 0→ max) Stroke Timing, d/u 1− 1.1 φ: position angle θ: elevation angle α: angle of attack of the wing β: stroke-place angle χ: body angle Figure 1.7: Schematic diagram of coordinate systems and wing kinematics for angles corresponding to Table 1.2 [7]. formed by the maximum and minimum sweep positions. The stroke plane angle and body angle vary with flight speed. The wing-beat kinematics can be described by the position angle φ, elevation angle θ, and the feathering angle α (the angle of attack). Unlike the “insect” wing kinematics, “avian” kinematics are associated with flapping mostly in a vertical plane while sweeping wings forward and back. Wing twist along the span and wing folding allow for the adjustment of the wing’s wetted area [7]. This flapping motion is important because it generates both vertical and propulsive forces that enable flight. Figure 1.8 shows the wing tip path for various 16 Figure 1.8: Wing kinematics for various “avian” type flyers (a) albatross, fast gate; (b) pigeon, slow gate; (c) horseshoe bat, fast gate (d) horseshoe bat, slow gate. birds and bats in forward flight. Typically, figure eights and ovals are the two tip path shapes. The wing tip path can vary based on flight mode (e.g., take-off or landing) and patterns are highly dependent on flight speed. The sequence of motions is illustrated in Fig. 1.9. For “avian” type flyers, variable wing pitching plays a very important role during flight. This pitching is implemented using a very complex musculoskeletal system [54, 55]. Figure 1.10 shows the forces generated during forward flight by a flapping wing with positive wing camber during downstroke and upstroke. During downstroke, the wing is always at a positive angle of attack with significant upward inflow near the wing tip. The relative inflow tilts the lift vector forward so that it produces thrust. This is experienced more heavily near the wing tip. The pitch remains nose down (into the flapping direction) so the flow stays attached. This is because nose down pitching is the primary mode for creating a positive angle of attack during the upstroke. This positive angle of attack is achieved by creating significant upward twist in the wing. This is necessary because the downward inflow angle, due to the flapping motion, increases with increasing span position. If a 17 Figure 1.9: Perspective views of avian wing movements during normal forward flight [8]. 18 Figure 1.10: a) Upstroke and b) downstroke forces on a flapping wing in forward flight [9]. positive angle of attack is achieved, the lifting force will be upwards and backward, increasing the drag on the wing. However, if the angle of attack during upstroke is negative, the lift vector will be negative and forward, creating thrust but also causing negative lift. The motion of the wing tip is very important for increasing thrust and decreasing drag, while the interior of the wing primarily produces lift. There is a significant amount of flow near the wing tip where the flap motion is large [9, 56]. 1.4.1 Bio-inspired Flapping MAVs Observations of biological flyers in nature and systematic gust studies [57,58] on real insects have shown high maneuverability and gust tolerance capability for flapping wings. This is because of their ability to generate large, almost instantaneous, control forces and moments over just a couple of wing strokes. The following section highlights two current flapping-wing MAVs (entomopters and ornithopters). 19 Figure 1.11: University of Toronto Mentor Entomopter. [10] 1.4.1.1 Entomopters The entomopter is a versatile insect-like robot, often with the capability of flying and/or crawling. Though many entomopter type MAVs lack the ability to crawl, insect-like flying has been integrated into many designs. The first entomopter to demonstrate controlled free flight for both translation and hover was the SRI Institute/University of Toronto Mentor (see Fig. 1.11). The four wings are arranged in an X-configuration. This vehicle is controlled by onboard electronics that move fins and gyros to direct the downwash. The original Mentor design was powered by an internal combustion engine and had a 36 cm wingspan. The weight was 580 g. This design flew for one minute and was able to achieve hover for 20 seconds. A later version of the vehicle only weighed 440 g and used a Nickel Cadmium battery to power a brushless motor [10]. At the University of Maryland, a hybrid, hover-capable insect-based MAV 20 concept was proposed [11]. This vehicle would be called the Thrust Augmented Entomopter (TAE). The estimated weight would be 160 grams and the wingspan would be 218 mm. The proposed configuration would utilize flapping wings and a pusher propeller for thrust in forward flight. The wings would be configured for minimum power requirements in hover, and would not flap during forward flight. This MAV was at a concept stage and was never built. However, this concept would be the motivation for future entomopter designs. Figure 1.12 shows the concept drawing of this MAV design. (a) Thrust Augmented Entomopter in Airplane Mode. (b) Thrust Augmented Entomopter in Insect Mode. Figure 1.12: Thrust Augmented Entomopter [11]. 21 The Micromechanical Flying Insect (MFI) project was developed at the Uni- versity of California Berkeley [59]. This vehicle was designed to mimic a blow fly, Calliphora (see Fig. 1.13(a)). The design specifications were that the mass would be 100 mg, the wing length would be 11 mm, and the wing beat frequency would be 150 Hz. Actuation would be done by using piezo bimorphs to connect the wings to a four-bar mechanism. This would only utilize about 10 mW of actuator power. The outcome of this study was that this vehicle did generation the adequate lift equal to weight (500µN) [59–63]. An offshoot of the MFI program’s research was the Harvard Microrobotic Fly (see Fig. 1.13(b)). This was expected to produce the smallest entomopter to demonstrate tethered flight [64]. This vehicle used only one piezoelectric actuator and weighs 56 milligrams. The wings were rigid, but a root joint allowed feathering in response to aerodynamic and inertial forces. The flapping and feathering strokes were limited to 100◦ using mechanical stops. The vehicle was able to exert a vertical force that was 3.6 times its own weight. This vehicle has demonstrated the ability to achieve thrust equal to weight, however, only a small number of controlled flights have been reported and the vehicle was powered from the ground. Ongoing modifications are being made to this vehicle to create an autonomous flight-capable vehicle weighing about 120 mg. The use of 3-D printing technology has greatly expanded possibilities for wing design, allowing wing shapes that mimic those of real insects or virtually any other shape. This reduces the time of a wing design cycle to a matter of minutes. Figure 1.13(c) shows a flapping-wing prototype produced at Cornell University. This 22 (a) UC Berkley MFI [59]. (b) Harvard University Microbot [64]. (c) 3-D printed flapping-hovering insect prototype [3]. Figure 1.13: Light-weight MAVs: Concept drawing of the Micromechanical Flying Insect (MFI), the prototype Microbot, and 3-D printed hover prototype. 23 research focused on developing an extremely light-weight hovering MAV using 3-D printed wings and mechanical parts. This prototype has a mass of 3.89 grams and has demonstrated an 85-second passively stable untethered hovering flight [3]. There have been many hovering MAVs produced and Table 1.3 shows the characteristics for many of these vehicles. Figure 1.14 shows two flapping-wing MAVs developed at the University of Maryland that have been successfully flight tested. Specifically, Fig. 1.14(a) shows the Daedal prototype, which was developed as a collaboration between Daedalus Flight Systems, LLC and the University of Maryland’s Autonomous Vehicle Laboratory. This 12 gram MAV is capable of free-flight at a wing-beat frequency of 25 Hz. Each wing’s length is 91 mm and the wing skin is made of 5µm thick Mylar sheet [65, 66]. The Daedal prototype did demonstrate both hover and free-flight; however, vehicle flight control was limited by actuator latency. The vehicle’s endurance was only 3–4 minutes. Currently, there are plans to reduce the weight to reach a GTOW of 10 grams. Figure 1.14(b) shows the University of Maryland’s “64 Gram Flapper” flapping MAV. The actual weight is 63.8 grams and control for this MAV was enabled by using a wing-based control system. This vehicle has successfully demonstrated both hover and controlled flight. However, the “64 Gram Flapper” experienced some issues with control and over-heating. DARPA sponsored a Nano Air Vehicle program centered around developing an ornithopter technology demonstrator. Through this program, AeroVironment developed a vehicle called the Hummingbird to satisfy the requirements (see Fig. 1.15). This vehicle is arguably the most notable flapping wing MAV built to date. The 24 (a) University of Maryland Daedal [65]. (b) University of Maryland “64 Gram Flapper” [67]. Figure 1.14: Flapping-wing MAVs developed at the University of Maryland, a) Daedal and b) “64 Gram Flapper”. vehicle is comparable in size to a natural hummingbird in that the mass is 19 grams and the wingspan is 16.5 cm. The AeroVironment Hummingbird has the ability to hover for several minutes and has a maximum forward flight speed of about 6.7 m/s [68]. Table 1.3: Characteristics of existing entomopter designs [3]. Design Year Mass (g) Span (cm) Wings Endurance (s) Features Mentor 2002 580 36 4 >60 Gas Powered DelFly II 2006 16.07 28 4 480 Camera, R/C van Breugel 2007 24.2 45 8 33 Passively Stable Microbot 2007 0.06 3 2 N/A Piezoelectric DelFly Micro 2008 3.07 10 4 N/A Camera, R/C NAV 2009 10 7.5 2 20 Active Wing Pitch 3-D Print 2010 3.89 14.3 4 85 3-D Print Wings 1.4.1.2 Ornithopters Avian-based flight utilizes a nearly vertical, up-down, flapping motion. Often, the wings are comprised of flexible membrane skins covering a musculoskeletal system controlled by small synchronized variations of angle of incidence. Similar to 25 Figure 1.15: AeroVironment’s Nano Hummingbird [12]. fixed-wing vehicles, ornithopters cannot hover, and they require an initial airspeed for taking off. Also, commercial ornithopters can have wing spans on the order of one to two meters and are often much heavier than the MAV design requirements. However, ornithopter-type MAVs have not been able to replicate the performance of birds because of their inability to copy the morphing structures of avian wings. Nevertheless, ornithopters do provide insight into low Reynolds number flight, and they are used by both hobbyists and researchers. One successful commercial ornithopter is the Kinkade model (see Fig. 1.16). The Kinkade Parkhawk series of ornithopters has a wing span of 107 cm and weighs about 454 grams. Flight speeds range from 2.5 m/s to 9 m/s. The total stoke amplitude is 55◦ and experiment flap frequencies were 2-6 Hz. The wings on this vehicle employ passive morphing behavior during flight which allows for variations in the local bending and twisting angles. Several flight tests (using motion tracking) 26 Figure 1.16: Commercial Kinkade ParkHawk onithopter on UMD Morpheus Lab test mount [13]. (a) Delft Technical University’s DelFly I [70]. (b) AeroVironment’s Microbat [71]. Figure 1.17: MAV-scale ornithopters . and predictive analyses have been performed on this vehicle at the University of Maryland’s Morpheus Laboratory [13,69]. There have been several avian-type MAVs produced specifically for research (see Fig. 1.17). Figure 1.17(a) show’s the DelFly ornithopter created by Delft University [70]. The vehicle’s wing beat frequency is 6 Hz; wingspan is 350 mm. The weight is 17 grams, and it can cruise at speeds of 1.8 m/s. However, the endurance is only 12 minutes. Figure 1.17(b) shows Aerovironment’s Microbat weighs about 12 grams, which uses microelectromechanical systems (MEMS) technology. The wing span is 14 cm but like many avian-based ornithopters, this vehicle can not hover [71]. 27 The flapping-wing concept generates much interest within the aerodynamics community because of the relatively strong, unsteady, three-dimensional vortical flow structures that are produced on small span, low aspect ratio wings. Although some successful flapping MAVs have been produced, flapping-wing research has lagged fixed and rotating wing research, and current flappers have lower hovering efficiency ratings compared to rotating wing vehicles. The flapping-wing must be thoroughly investigated from the perspective of fundamental aerodynamics. Specifically, there is a need for greater understanding of the fundamental lift and thrust mechanisms, as well as it requires systematic performance measurements that can be correlated with predictions in order to develop a viable mechanism for an MAV. An understanding of the design challenges may then allow more efficient and reliable MAVs to be designed. 1.5 Challenges to the Development of MAVs There are several technical barriers to producing mission capable MAVs. Dif- ficulties in the theoretical analysis, experimentation, replication of natural wing kinematics, and numerical simulation of the capabilities of the various competing MAV concepts have led to a wide range of assessments and opinions of their actual capabilities. There are also issues such as small-scale power generation and storage, navigation and communications, and efficient propulsion and aerodynamics [14]. Conventional rotors have been scaled down for MAVs, but have suffered from consid- erably lower aerodynamic efficiencies than those achieved at larger scales [1, 72,73]. 28 The same lower efficiencies are also seen for flapping wings. Research efforts have produced several hover-capable flapping-wing MAVs with varying performance like the Microbot [61,64], Daedal [65], Delfly [70], and Nano Hummingbird [68]. However, these vehicles have largely been developed with a focus on an operational vehicle design rather than an in-depth understanding of the fundamental aerodynamics associated with flapping MAVs. This approach has led to the development of several MAV prototypes, but these vehicles are only capable of flight with limited payload and in controlled environments. Challenges facing MAV development are discussed below. 1.5.1 Vehicle Control/Actuation Manipulation of the control surfaces is a crucial task for MAV design, especially at the small scale. This is because actuation can define the maximum possible performance. The two main criteria for an actuation concept are the power to weight ratio and the technical feasibility. A further constraint in the selection of an actuator concept for flapping MAV applications is wing beat frequency (speed of actuation). There are two types of actuator categories: rotary and linear. Rotary type actuators include brushed DC motors, brushless DC motors, micro internal combustion engines (ICEs), and piezoelectric motors. Linear actuators include piezoelectric ceramic, shape memory alloy (SMA), magnetostrictor, solenoid, and electroactive polymers. Some of these actuation concepts are under development and have not yet matured for MAV applications, but the most promising are electric motors and solenoids, 29 SMA wires, and piezo elements. Electric motors and solenoids are based on the same physical principles, with the difference being that solenoids deliver a linear motion, while motors rotate. Electric motors and solenoids can be scaled down to 1-2 grams; however, solenoids have a lower energy density and reaction time in comparison to electric motors. Miniaturized electric motors do suffer from a 30% decrease in efficiency when compared to large motors [39]. In general, the efficiency of a rotating motor with gearing has been preferable to a linear solenoid. SMA wires imitate natural muscles such as the ones of birds and performance is very good within the expected range. This material can contract when heated and expand when cooled. Some SMA actuators are commercially available for micro-scale applications, and a theoretical analysis was conducted for the BAT-MAV program [74]. However, the disadvantage is in the speed and fatigue life of such actuators. This is because of the slow time associated with heat dissipation. With this, SMA actuators are not fast enough for application in flapping MAVs [75]. Another interesting technology is piezo based actuators. These are actuators that use piezoelectric material that responds mechanically when an electric charge is applied. Depending on the size, they can be miniaturized to a weight of less than 1 g [76]. They can also achieve very high frequencies. The disadvantage of such elements is the high voltage (typically greater than 100 V) needed to induce the contraction of the material [77,78]. In summary, solenoids must have a very high strain rate but a low energy density. SMA wires are interesting because they mimic the behavior of natural muscle, 30 but their performance is not yet suitable [14]. They do not have the bandwidth required for high frequency flapping. They may become an option in future if the performance is improved. Piezo elements are inadequate due to the high voltage and small deflections. Finally, micro-scale electric motors are currently the most suitable (light-weight and efficient) devices for the construction of a flapping wing mechanisms and therefore often used for MAVs. 1.5.2 Energy Sources All MAVs must draw power from a fuel source to function (e.g., to power motors and sensors). Figure 1.18 shows that hydrogen and gasoline currently have the highest energy density of the energy sources reviewed. Even though fossil fuels have a higher energy density than current battery technology, the disadvantages are a increased noise signature (internal combustion engines) and increased emissions. The main issue is the lower efficiency (¡ 5%) of the small scale engines. The energy sources for robotic flyers includes conventional batteries, fuel cells, ultra capacitors, and solar cells. Most small-scale flight vehicles use batteries. Various consumer electronics industries have driven up the power density of batteries. Mainly driven by the cell phone industry, lithium polymer batteries are the only commercially available technology that can satisfy the requirements of insect-sized MAVs [79, 80]. However, the capacity of lithium polymer batteries decreases at high discharge rates. The properties of the capacity derating depend heavily on the battery design and manufacturing parameters. As a result, it is 31 Figure 1.18: Energy storage systems energy densities [14]. difficult to estimate battery performance without a specific battery in mind; at the same time, it is important to consider derating because the high power requirements of hovering MAVs inevitably translate to high battery discharge rates. The derating relationship may be nonlinear. Promising new technologies include micro solid oxide fuel cells [81] (where chemical fuel is transformed into electric current through an electrochemical process), lithium batteries with silicon nanowire anodes [82], ultra capacitors, solar cells, and lithium air batteries [83]. 1.5.3 Propulsion and Platform Integration The devices by which most MAVs are propelled are electric motors, as they are light-weight, controllable, reliable, and efficient. Motor overheating is one disad- vantage, though this can be mitigated by proper sizing. Small internal combustion engines have been employed, but efficiency, noise, and system integration can be 32 problematic. Some promising alternatives include hybrid gas/electric systems [84,85] and miniature gas turbine generators [86], but these concepts have not yet been applied to MAVs because of the size limitations and difficulty in integrating the propulsion system into the rest of the vehicle components. In general, platform integration of MAV components is a very challenging endeavor. The main consideration is to keep the vehicle weight low while maintaining the functionality to complete a given mission. Because smaller structures and hence stiffer components are required, MAVs tend to be very rigid unless intentionally made flexible. One example of the incorporation of vehicle flexibility is at the University of Florida where interstitial (overlapped) wing materials were used [40]. The multi-functionality of components becomes very important when vehicle space is limited. Also, innovation in design becomes crucial when antennae or sensors need to be integrated into the wings, or the power source needs to be integrated with the fuselage structure (as examples) [87]. The sustained hover endurance goal of one hour has eluded MAV researchers and designers, especially for flapping wings. Advancements in MEMS, improvements in battery technology, and clever ideas for platform integration have advanced the state of the art. However, motor efficiency and the battery specific energy and power are not at the level to meet the endurance goals. To improve endurance, there is an opportunity to optimize the aerodynamic design of MAVs for operation in a low Reynolds number regime. Considering the typical missions MAVs will have to perform, it is not difficult to conclude that these key performance parameters are the maximum lift producing capability (Cl,max), aerodynamic efficiency (lift per unit 33 power), maneuverability (control forces per unit body inertia), and robustness to external disturbances such as wind gusts, both in an open-loop (bare airframe) and closed-loop sense. Both rotating- and fixed-wing MAVs are highly vulnerable to the nonsteady aerodynamic environment, but flappers appear to stand out because the low Reynolds number effects are often overcome by leveraging unsteady aerodynamic effects. The next subsection contains a more in-depth discussion of low Reynolds number issues. 1.5.4 Low Reynolds Number Effects Flapping MAVs operate in the Reynolds number range of 104–105. In this range, there are aerodynamic sensitivities to wing profile and overall geometry, low lift- to-drag ratios, increased susceptibility to flow separation, increased boundary layer thicknesses, and higher relative viscous losses [7,19]. These unfavorable aerodynamic characteristics are seen when scaled down conventional rotors show considerably lower aerodynamic efficiency than those achieved at larger scales [1, 72]. Reference 15 discusses the aerodynamic efficiency for various hover-capable MAVs as well as avian and entomological flyers examined on the basis of their absolute maximum attainable efficiency. Though various other types of hover-capable MAV designs have also been developed, including biomimetic vehicles, most of the flapping wing MAVs still do not exceed the hovering efficiencies that have been obtained with their conventional counterparts in both hover and forward flight. A comparison of the aerodynamic efficiency for many different hover-capable 34 MAVs, as well as natural fliers, was performed based on the maximum attainable efficiency in the form of power loading [15]. The power loading, PL, of a flight vehicle is the ratio of the thrust, T , produced by the vertically lifting aerodynamic surfaces per unit of power, P , expended to produce thrust. To predict the absolute minimum power required to hover, ideal flow conditions (i.e., incompressible, inviscid, one-dimensional, and quasi-steady) are assumed. The power loading for a hovering vehicle can be related to wing loading or rotor disk loading, T/Ae, and the figure of merit, FM , by the following equation: PL = √ 2ρ FM √ T/Ae (1.1) In this equation, ρ is the fluid density; Ae is the effective area of the reciprocating wing stroke plane or rotor disk; and FM is the figure of merit, a measure of efficiency given by the ratio of actual to ideal power required. To achieve maximum hover endurance, a special effort must be made to obtain the best possible performance in hover, i.e., maximize power loading. For hovering vehicles, as effective disk loading decreases, power loading increases rapidly (see Eq. 1.1). The result is that a hover-capable vehicle must be very light and/or have a large wingspan or rotor disk area. MAV size constraints ensure particular attention is given to a reduction in vehicle weight through careful engineering design. In hover, thrust equals weight, so reducing vehicle weight results in a reduced disk loading, T/Ae. A flapping wing MAV concept that can achieve lower effective disk loadings through a reduction in vehicle weight can, for a fixed figure of merit, also achieve 35 Figure 1.19: Comparison of efficiency parameters for hover-capable MAV-scale fliers [15]. improved power loadings. Figure 1.19 shows measured and computed disk/wing and power loading data for various full-scale and small-scale hovering vehicle concepts with lines of constant FM given by momentum theory. For small-scale and large-scale rotary-wing vehicles, the theoretical maximum attainable levels of hovering efficiency are greater than or comparable to those achieved by natural fliers (when compared on the basis of equivalent wing stroke area/disk loading). However, current MAVs have generally achieved relatively low hovering efficiencies. Currently, large-scale rotors typically exhibit maximum figure of merit (FM) values between 0.7 and 0.8, whereas MAV- scale rotors have FM values near 0.4–0.5 [19]. Flapping MAVs have even lower efficiency values, FM < 0.2. As such, there is room for improvement in MAV 36 design to increase overall efficiency. To this end, an improved understanding of the aerodynamic issues associated with flight operations at the small-scale is needed, with special attention to performance and efficiency parameters. A key concern today, is the fact that, even though many of these studies, especially by Dickinson [20, 22, 88–90], shed light on the complex aerodynamics of insect flight operating at very low Reynolds numbers (Re < 1000), it is not very clear whether these conclusions would still hold true at MAV-scale Reynolds numbers, which are at least an order of magnitude higher than 1,000. In that case, the usefulness of these conclusions in formulating the design principles for an actual flapping-wing MAV are questionable. Therefore, systematic studies need to be performed at MAV-scale Reynolds numbers (Re > 10, 000). Also, these studies should not be limited to understanding the fundamental physics, but they should also be able to derive standard quantitative aerodynamic efficiency metrics (such as thrust/power and figure of merit) for a flapping wing in hover. This would clearly show how efficient or inefficient the flapping wing concept really is, in comparison to other hover-capable concepts such as a conventional helicopter rotor. One key challenge in developing efficient and robust FMAVs is overcoming aerodynamic issues attributed to flight at Reynolds numbers on the order of 104, such as low lift-to-drag ratios and increased susceptibility to flow separation [1, 91]. Though birds and insects have successfully adapted to such aerodynamic conditions and fly well in the lower Reynolds number regime, these aerodynamic challenges contribute to the low efficiency observed in MAVs [7]. Reynolds number (Re) is defined as the ratio of vehicle inertial force to viscous 37 force. Reynolds number is a non-dimensionalized parameter that can be seen in Navier-Stokes equations. Navier-Stokes equations are the fluid flow equations and have the following form for an incompressible fluid [17]: ∂(~u) ∂(t) + ( ~u · ~5 ) ~u = ~5p+ 1 Re ~5 2 ~u (1.2) where, 5 is the spatial gradient operator normalized by some length scale (L), u is the velocity vector normalized by some velocity scale (U), t is the time normalized by the convective time scale (L/U), p is the pressure (normalized by the dynamic pressure ρU2), and ρ is the fluid density. The Re is defined as: Re = ( ρUL µ ) (1.3) Here, ρ is again the fluid density, and U is the free-stream velocity. The reference length of the object is represented by L, which is usually chord length of a wing, and µ is the fluid viscosity. Flow separation and transition are highly sensitive to Reynolds number, pressure gradient, and the disturbance environment. The aerodynamic characteristics of a vehicle, (e.g., wing design) in turn affect the static, dynamic, and aeroelastic stability of the entire vehicle. Therefore the management of the boundary layer for a particular low Reynolds number vehicle design is critical. Figures 1.20 shows the variation in airfoil section Clmax, Cdmin, and (Cl/Cd)max as a function of Reynolds number. The general trend is that performance is greatly reduced as the chord Reynolds number decreases below 100,000. Again, MAVs operate 38 at Reynolds numbers less than 105. This means that as vehicle control (thrusting and lifting) surfaces are made smaller and flow velocities are reduced, these vehicles become inherently less efficient. In the range between 1, 000 < Re < 10, 000 the boundary layer flow is laminar and it becomes difficult to cause transition to turbulent flow. For chord Reynolds numbers between 10,000 and 30,000, the boundary layer is completely laminar. Experience with hand-launched glider models indicates that when the boundary layer separates, it does not reattach. Figure 1.20(a) shows the variation in lift coefficient Cl with Reynolds number. This parameter is dependent on the local airfoil camber which affects the relative angle of attack. Airfoils with large camber have high lift coefficients at small and even negative angles of attack, but they also achieve stall at angles of attack lower than small cambered or symmetric airfoils. Figure 1.20(b) shows the drag coefficient (Cd) as a function of Reynolds number. The two components of drag are the profile (form drag, skin friction and interference drag) and induced drag (drag required to maintain lift). Induced drag increases when the aspect ratio is decreased (specific to MAVs). Figure 1.20(c) shows the lift-to-drag ratio for an airfoil. This is a measure of the aerodynamic efficiency in that aerodynamic efficiency is optimal at the maximum lift-to-drag ratio. The two factors that affect lift-to-drag ratio include wing planform and local airfoil shape. For Reynolds numbers < 105, wings must be designed to reduce the profile drag by avoiding the formation of laminar separation bubbles. A laminar separation bubble is shown in Fig. 1.21 and is caused by an adverse pressure gradient. This 39 (a) Maximum lift coefficient. (b) Maximum drag coefficient. (c) Maximum lift-to-drag ratio. Figure 1.20: Airfoil performance as a function of Reynolds number [16]. 40 Figure 1.21: Laminar separation bubble [17]. adverse pressure gradient is the reason the boundary layer starts to separate. The laminar separated shear flow is unstable and transitions to a turbulent separated shear flow. The turbulence then transports momentum from the free-stream, across the shear layer, and down towards the surface. When the momentum transport is sufficient, the turbulent boundary layer is considered to be reattached to the surface, thus closing the separation bubble [19]. These bubbles become larger as the Reynolds number decreases, usually resulting in a rapid deterioration in performance, i.e., substantial decrease in L/D. Low Reynolds number airfoils will experience flow separation and stall at low angles of attack. However, a sharpened leading edge of a thin flat (or curved arc) airfoil acts as a boundary layer trip. This causes the flow to become turbulent once it passes the leading edge of the airfoil. Even with these types of modifications, the lifting performance of optimized wings are not enough to handle the abrupt change in performance brought on at low Reynolds numbers. This is the reason that many 41 current MAVs utilize curved thin plates or membrane wings. MAV designers are interested in the 30,000 to 70,000 Reynolds number range because this is where many model-scale aircraft operate. In this range, it has been seen that thick airfoils (i.e., 6% and above) can suffer from degraded performance as a result of laminar separation [32]. Near the 50,000 Reynolds number, the free shear layer after laminar separation does not transition to turbulent flow in time to reattach. Thin airfoil sections (i.e., less than 6% thick) operate reasonably well at the upper end of this regime [92]. For the Reynolds number range 70, 000 < Re < 200, 000, where some small radio controlled airplanes also operate, laminar flow can be observed on the wings and airfoil performance is improved. 1.6 Unsteady Lift Enhancement Mechanisms Like helicopters, flapping wings generate a very complex aerodynamic flowfield with large vortical structures on the wings. In general, dynamically pitching airfoils allow for the flow to remain attached to angles of attack greater than what would cause stall during static conditions. This is the case for both rotary- and flapping- wing MAVs, but especially for FMAVs. This mean that the unsteady nature of the flow can be used to enhance the lift characteristics for flapping wings. With careful engineering, designers can bring forth the next-generation of MAVs to bridge the gap between rotating- and fixed- wing vehicles. A key is that flapping-wing micro air vehicles (FMAVs) must produce high maximum lift coefficients by generating strong three-dimensional vortices on small, low aspect ratio wings. 42 Many studies have been conducted in relation to the shedding of concentrated vorticity from the leading edge of a flapping wing as a significant source of lift enhancement at the MAV scale. However, steady and quasi-steady aerodynamic theory, experiments, and modeling (e.g., computational fluid dynamics or CFD) have not yet fully explained the ability of flapping-wing flyers to generate the values of lift needed to achieve hover and to fly forward at the observed speeds [18, 20, 22]. It has been established that MAVs suffer from degraded performance, but biological flyers perform well in this low Reynolds number aerodynamic environment. This is done by leveraging several unsteady aerodynamic mechanisms to improve performance. These mechanisms are discussed in the next two subsections. 1.6.1 The Leading Edge Vortex The leading edge vortex (LEV) is the most significant source of lift enhancement for a flapping wing. A LEV forms on the upper surface of the wing when flow separates over the sharp leading edge and rolls up into a coherent vortex. The lift increase occurs because the presence of a LEV lowers the local pressure on the upper surface of the wing, thus increasing the overall lift. Figure 1.22 shows an illustration of a two-dimensional cross section of the LEV, and the three-dimensional segment of the LEV, showing the axial flow component. Here the streamlines roll up into the vortex as it spirals towards the wing tip (spanwise direction) [18]. The spiraling flow structure is caused by spanwise fluid transport. Ellington was the first to examine the LEV using flow visualization on a live 43 Figure 1.22: 3-D conical LEV structure: a) LEV cross-section; b) three-dimensional structure and spanwise flow (note that le is the leading edge, te is the trailing edge, dss is the dividing stream surface) [18]. tethered hawkmoth at chord Reynolds numbers near 2,000 [18]. His pioneering research also explored the flow on flapping wings at Reynolds numbers from 4,000 to 7,000 and sparked much follow-up research [18,93,94]. It was subsequently concluded that because insects fly with wings at a high angle of attack and at low Reynolds numbers, they must utilize the lift-enhancing effects of a LEV to produce sufficient lift for flight [95, 96]. Experiments were continued on mechanical models to study the LEV and it was shown that the prominent flow structure generated on these wings was a three-dimensional LEV. The main observations were that 1) the LEV was smaller towards the wing root and larger outboard, 2) near the wing tip, the LEV separated and joined with the tip vortex, and 3) LEV formation resembled the process of dynamic stall. Dynamic stall occurs when the angle of attack of an airfoil increases rapidly, causing a delay in the onset of flow separation and the formation of a vortex on the upper surface of the wing [19, 97]. The vortex grows in strength before detaching from the leading edge and convects over the wing and into the wake, resulting 44 Figure 1.23: Airloads and flow structures on a 2-D airfoil during the dynamic stall process [19]. in overall wing stall. Figure 1.23 shows a schematic of the flow morphology and unsteady airloads during the dynamic stall process for a two-dimensional airfoil. The five stages indicate the various levels of flow separation, LEV growth, LEV shedding/convection, and flow reattachment as the wing changes angle of attack [98]. Subsequent experiments have confirmed the lift-generating characteristics of LEVs on flapping wings [18, 20, 20, 22, 89, 90, 99]. Based on these studies, it is 45 concluded that an attached LEV is the key high-lift source on a flapping wing. While experiments on live birds and insects can be extremely challenging, the use of mechanical models that mimic the flapping kinematics of the wings seen in nature has allowed for further and more detailed studies of the LEV [73, 100, 101]. Such mechanical models make it easier to systematically study the effects of variations in Reynolds number, angle of attack, and wing aspect ratio on the behavior of the LEV. Although most of these experiments have been conducted on rigid flapping wings, they have provided valuable insight into the overall aerodynamic effects produced by LEVs. Dickinson and colleagues independently observed a stable LEV on a dynamically scaled robotic fruitfly wing model operating in an oil tank at much lower Reynolds numbers (Re = 100 - 1000) [20, 22, 89, 90]. Large lift and drag coefficients were measured by the force balance at the wing root, and the magnitude of these forces were directly related to the strength of the LEV, clearly showing that the LEV is the leading contributor to the aerodynamic forces on a flapping wing. The important consideration was in how flapping wings could produce such a coherent stable vortical structure, where the flow would not separate on a steady translating 2-D wing at such high angles of attack. Ellington postulated that this was because there is a strong spiraling spanwise flow, which would transport vorticity from the LEV along the wing length to the tip vortex that prevents the LEV from shedding [18,94]. Dickinson showed that at Reynolds numbers matching the flow relevant to most insects (Re < 1, 000), flapping wings do not produce a spiraling vortex, and the stability of the LEV was attributed to the decrease in effective angle of attack 46 Figure 1.24: Illustration of wing with flow restricted by baffles, and vorticity with overlaid velocity field (Note: Blue denotes a vortex) [20]. because of the downward flow induced by the strong tip vortices [20]. Birch and Dickinson further investigated the stability of the LEV on a dynamically scaled fruit fly wing (Re = 160). To assess the effect of spanwise flow on the LEV, they applied fences and baffles on the wing’s upper surface. This would restrict the spanwise flow. Even after placing these baffles, the flow still did not detach from the wing. It was stated that there must be some Reynolds number dependence pertaining to LEV stability. Figure 1.24 show some flowfield results of the LEV and a schematic of the wing. Recent studies by Lentink and Dickinson indicate that the LEV is stabilized by the ‘quasi-steady’ centripetal and Coriolis accelerations of the fluid particles close to the wing’s surface because of the propeller-like sweep of the wing, especially when the radius to chord ratio is low (lower Rossby number) [101]. This study also showed that LEV stability is not dependent on Reynolds number, at least over the range 47 most relevant for insect flight (100 < Re < 14, 000). Although early analyses of flapping wings has identified high maximum vertical force coefficients, the LEV is not the only lift augmenting mechanism for a flapping wing. 1.6.2 Other Lift Augmenting Mechanisms The other unsteady mechanisms for a flapping wing are 1) clap and fling, 2) rotational lift, 3) added mass, and 4) wake capture. Studies have exposed a general understanding of the unsteady aerodynamic forces produced by flapping wings [99,102,103]. The Weis-Fogh or clap and fling mechanism is a phenomena initiated when two wings move close together approaching pronation or supination. This fluid dynamic process was discovered to be prevalent for small insects [104,105]. In general, fluid is pulled into the void left by the wings as the move away from each other. An illustration of this phenomena is shown in Fig. 1.25. Figure 1.25(A-C) shows the clap phase of the mechanism. Here, the leading edges of the wings come together during wing pronation/supination until they are in close proximity and almost parallel to one another. The fluid is pushed out by the wings. The wings then rotate about the trailing edge, and as the wings separate, a vortex pair forms on the leading edges. During the fling phase, shown in Fig. 1.25(F-D), the wings separate, and the fluid rushes into the growing gap increasing the circulation. Lehmann et al. [106] examined the lift augmentation produced by the clap-and-fling mechanism on a scaled fruit fly wing at a chord. They compared force measurements from a single flapping-wing 48 Figure 1.25: Illustration of the clap-and-fling/Weis-Fogh mechanism [21]: A) at end of half-stroke wings moving toward each other; B) wings continue approach to each other; C) fluid is moved out by wings; D) wings rotation about trailing edge; E) fluid pulled into the void left by the wings; F) LEV formation begins as wings move apart. configuration to that of a dual flapping-wing configuration. The clap and fling mechanism increased the lift by as much as 17% (and the LEV circulation in the fling phase by 32%) for Reynolds number between 50–200. However, it was found that as Reynolds number increased, lift enhancement resulting from the clap-and-fling motion began to decrease. A second mechanism of lift augmentation is circulation from wing rotation (flip). This is caused by the rotation of the wings and is especially prevalent at the end of the hovering stroke (pronation and supination). This mechanism is caused by Kramer’s effect, which is closely related to delay/dynamic stall [19, 21]. The wing goes through a large rotation at a very high angular rate towards the end of 49 the stroke. The mechanism of circulation generated by rotation is similar to the Magnus effect. Though this effect cannot be applied directly to the flat surface of a flapping wing, where the establishment of the Kutta condition is necessary [17], this mechanism is enhanced when the wings are flexible (favorable deformation) [107]. The third mechanism is wake capture (also called wing-wake-interaction), caused by the reciprocating nature of the wing flapping. Figure 1.26 shows a schematic of the wake capture phenomena. This happens immediately after stroke- reversal (see Fig. 1.26(AF)). At the ends of the stroke (pronation and supination), the wing reverses direction, the incoming free-stream velocity may also be increased by the induced effects of the incoming shed wake or wakes (see Fig. 1.26(C-D)). It is even possible for a flapping wing to come into contact with a shed LEV at the beginning of its half-stroke (see Fig. 1.26(E-F)). Therefore, a better understanding of this interaction may lead to a MAV with more highly optimized wing kinematics that could take advantage of this mechanism of lift production. Work done by Dickinson et al. [22] shows that lift enhancement could be achieved for their flapping-wing mechanism. Confirmation that the enhanced lift force due to wake capture was shown by observing that the lift peak continued even after the wing motion was terminated. This showed that the shed wake was the energy source for lift generation. It was also shown that lift increases from the wake capture mechanism are dependent on the wing kinematics at the end of each half-stroke. Finally, there are added or “virtual” mass effects. As the wing accelerates (because of the rotational and translational motion), it encounters a reaction force due to the inertia of the accelerated fluid [104,108,109]. This reactive, non-circulatory 50 Figure 1.26: Illustration of the wake capture mechanism [21]: (A) After translation, wing generates vorticity at both the leading and trailing edges, (B) Vortices induce a strong velocity field in the intervening region, (C/D) wing comes to a halt and reverses stroke (D/E), wing translates in the opposite directions [21]. force falls outside standard circulation based steady-state analysis, but can be modeled by additional mass. Concerning MAVs, Ansari et al. [108] noted that flapping wings cause high accelerations and rapid stroke reversals making forces significant even though the mass of the wing may be small. In general the inertia of the flapping wing is increased by the mass of the accelerated fluid [104,106,110]. Because added mass forces typically occur at the same time as the circulatory forces, it is usually difficult to measure them in isolation [110]. The added mass force is often modeled in quasi-steady terms, using a time-invariant added-mass coefficient, and any time dependence is implicit due to the time course of wing acceleration [88]. In general, all unsteady phenomena, except the LEV, persist for only a short duration during the wing stroke and also depend on the timing of rotation. Therefore, they do not contribute much to the time-averaged wing forces. Also, very high rotational lift because of wing rotation was only observed at very low Reynolds 51 numbers (Re < 1000) and is not of much consequence for flapping-wing MAVs that operate above that range. Therefore, for flapping wing studies, which are geared towards MAV development, the strength and stability of the LEV should be the focus because the LEV does contribute the most to the time-averaged forces independent of the small variations in wing kinematics. 1.7 Review of Flapping Wing Experimental Studies While many studies of the unsteady aerodynamics of pitching and plunging air- foils predate interest in flapping wing aerodynamics [111–113], the flows produced by insects and birds (and biomimetic mechanisms) have presented many new challenges in the areas of flow measurement and simulation. In the early 1900s, Knoller and Betz were the first to identify that positive propulsive thrust is produced due to pure plunging motion of a wing in freestream. This is because the wing motion results in a change in the effective angle of attack during the up-stroke and downstroke [114,115]. As mentioned previously, later work (e.g., References 20,89) has shown that a strong LEV forms on such a wing and likely contributes to lift production, but there is still a need for additional research on flapping wings to understand how they produce lift. There have been several studies conducted on live insects and birds [116,117]. Mainly, PIV and flow visualization have been used to analyze the flowfield on live insect and bird flapping-wings. Bomphery et al. [118] used two-dimensional PIV to examine the velocity fields produced on a tethered hawkmoth. They performed PIV at several locations in the wing stroke and compiled these results into a time- 52 resolved history of the flowfield. PIV measurements were taken of the flowfield of a free-flying hummingbird by Warrick et al. [119]. They examined the complex wake structure as well as how the shed vortices evolved and diffused into the far-field wake. Concerning wing flexibility, Wang et al. showed that dragonflies improved their aerodynamic efficiency by producing chordwise camber, thereby significantly enhancing the thrust [120]. Other research on biological flyers has come from the zoology and biology community and has focused on the wing behavior. For example, it was stated that hover-capable birds and insects (e.g., hummingbirds and bumble bees) flap at wing beat frequencies on the order of 30 to 90 Hz [121]. It was shown that the wing motion throughout the stroke is generally in the horizontal plane. Larger flyers such as song birds, eagles, and condors (wing beat frequencies 3 – 20 Hz) [121] often flap their wings mostly in the vertical plane while changing their wing planform to accommodate different flight modes. Over the flap cycle, dynamic twist (temporal change in twist) is utilized to vary the angle of attack along with the flapping motion. Later, it was discovered that large changes in wing pitch angle are necessary to produce lift for hover. Many of the lift enhancement techniques discussed earlier (e.g., leading edge vortices, Weis-Fogh mechanism, and rotational circulation) used to increase thrust production over the stroke [99,102,103,122] were discovered by observing natural flyers. However, measurements on live animals are extremely challenging because specimens can vary from test to test, fatigue of the specimen can be an issue, and it is not possible to perform direct power measurements on live animals. These issues contribute doubt to the reported efficiency levels of biological 53 flappers. Therefore, researchers have turned to (often simplified) mechanical models that mimic the flapping behavior. The inherently three-dimensional nature of the flow structure, and the overall dependence of the flow on the wing kinematics, geometry, and wing-wake interactions, make a simplified experiment useful for understanding the flow physics. The reason for many of the simplified flapping mechanism is because it is difficult to reproduce small scale, complex reciprocating motions at high flapping frequencies. There have been many simplified rigid revolving wing models used to mimic various aspects of the wing stroke, especially the translational part of the stroke [101,123–128]. These efforts have typically focused on a simplified wing model (e.g., a rigid flat plate) and kinematics (e.g., pitching, plunging, translating, rotating) in an attempt to identify the flow structures that produce lift, and shed light on the wing geometry and kinematics that most efficiently produce these structures. One such model is Dickinson’s robofly (see Fig. 1.27) [22]. This flapper was designed to model the kinematics of a fruit fly by flapping in mineral oil at Reynolds numbers O(100). Both wings were capable of rotational motion about three axes and were controlled by six stepper motors for force sensors and digital particle image velocimetry measurements. The main focus was on rigid wings. Revolving-wing experiments in oil and water tanks do aid in gaining a better understanding of the vortical structures produced on flapping wings. However, the major technical barrier for such experiments is the contamination of the inertial force with the measured aerodynamic forces. The inertial force could be much larger than the aerodynamic forces. Therefore, experiments conducted for the same kinematics in 54 Figure 1.27: Robofly oil tank test setup [22]. air as well as in inertial can be performed so that the vacuum forces can be separated from the total forces. This allows for the isolation of the aerodynamic forces [129]. Even though this may sound conceptually straightforward, implementing such a technique and obtaining accurate aerodynamic data is quite challenging because the inertial forces could be much larger than aerodynamic forces, especially in the plane of flapping. Therefore, developing a good experimental procedure for inertial force subtraction is useful. This is the primary reason why experiments should be carried out in air, even though experiments may be executed more easily in a higher density medium such as water or oil. Many notable flapping wing studies have been conducted in air at the University of Maryland on a biomimetic flapping-wing device [73,100]. This device was operated at a mean chord Reynolds number of about 15,500 and is shown in Fig. 1.28. This device was capable of single-degree-of-freedom flapping. Averaged air loads were 55 Figure 1.28: Passive pitch, bi-stable single-degree-of-freedom flapping mechanism [23]. measured using a load cell with piezo-resistive strain gauges. There were two key observations. The inertial loads constituted the major portion of the total loads acting on the wings tested and the propulsive thrust decreased at higher frequencies. Flow visualization and PIV experiments showed that there was a continual formation, convection and shedding process of the non-stable LEV, with several LEVs formed on the upper surface of the wing during a single wing stroke [73]. Force and power measurements performed on dynamically scaled insect-wing models have revealed some of the key physical phenomena contributing to the large aerodynamic forces on insect wings [20, 22, 89, 90, 101]. Birch and Dickinson used PIV to quantify the flow around the wing throughout the entire wing stroke of a dynamically scaled fruit fly wing in an oil tank [89]. More recently, Hart et al. [130] experimentally measured the instantaneous lift and drag of a pitching and plunging small aspect ratio rigid wing. This study included PIV measurements of the flowfield at one spanwise location (a very sparse data set). Here, it was shown that MAV-scale 56 flapping wings, which generally have a very low aspect ratio, experience significant 3-D effects that play an important role in maintaining flight. Though this study does not explain the flow physics in much detail, the data is still useful for CFD validation. Yuan et al. [131] conducted root based flapping experiments on a rigid wing that mimicked avian wing flight kinematics. The focus was on obtaining instantaneous force measurements over a flap cycle. Rigid-wing experiments such as those performed by Yuan et al. are useful as valuable insight can be gained regarding the fundamental flow physics [128,132]. However, in nature, the wings of birds and insects are very light-weight and highly flexible. Rigid flapping wings do not mimic the aeroelastic effects inherent to the wing motions seen in nature or those produced by flight-capable MAVs [133]. For experiments, mechanical models make it easier to systematically study the effects of variations in Reynolds number, angle of attack, and wing aspect ratio on the behavior of the LEV. A detailed understanding of the unsteady, and three-dimensional aerodynamics and their complex interaction on the flexible wing and the flow is critical to understanding flapping MAV performance. Some flowfield and airloads experiments on flexible wings have been conducted in hover [134–136], but very few studies have been carried out on flexible flapping wings in forward flight. There have been even fewer experiments performed in the lower Reynolds number regime on flexible flapping wings for both hover and forward flight. Tanaka et al. showed using experiments, hoverfly based wings (at the insect scale) could exploit intrinsic compliance, but a rigid wing could be more suitable for producing large lift [137]. For this study, there was no flowfield analysis, no forward flight cases, no analysis on the effect of wing flexibility on drag, and 57 no modeling component. Ho et al. and Heathcote conducted rigid and flexible flapping wing wind tunnel experiments [138–140]. The main parameters varied were flap frequency and freestream velocity. While flexibility along leading edge leads to deteriorating lift performance, flexibility near the trailing edge actually improves propulsive thrust [138]. Heathcote performed one of the few aeroelastic experimental studies for forward flight conditions. Force and flowfield measurements were made on a heaving rectangular wing (AR = 6) in a water tunnel at Reynolds numbers between 10,000–30,000 [139]. Only the spanwise flexibility (highly flexible, moderately flexible, and rigid) was varied and the wings were stiff in the chordwise direction. The conclusion of this study was that a limited degree of flexibility was very advantageous, for a wing with intermediate flexibility, propulsive thrust was increased by 50%, and for the highly flexible wing, propulsive thrust was reduced. This experiment has been used for validation in many studies because it is one of the very few studies involving structurally characterized flexible wings. Malhan et al. [141] also carried out experiments on flexible flapping wings in air and showed that wing flexibility plays an important role in aerodynamic performance. More recently, Malhan et al. validated computational fluid dynamics (CFD) with average force measurements for a MAV-scale flapping wing in forward flight [142]. However, detailed flowfield studies were not performed to analyze the flow state along the wing during the stroke. It was highlighted that it would be useful to conduct flowfield experiments where the unsteady phenomena occur, especially during the dynamic process. Vortex shedding was caused by the flapping wing motion, but the cause of the dynamic stall conditions is not always clear. 58 Though valuable insight can be gained into the fundamental flow physics of flapping wings [128], most of the flow experiments have utilized rigid wings, which do not mimic the aeroelastic effects inherent to the wing motions seen in nature or those produced by flight-capable MAVs [133]. Even with the limited number of MAV-scale experiments on flexible wings in forward flight, it is too time consuming to test every parameter variation. Therefore this is an area that needs more exploration. Detailed flowfield experiments are extremely challenging at the small scale and some attempts have been made at overcoming those challenges. However, the approach should be to analyze the flowfield with airloads for various wing types so that forces can be placed in the context of what is happening in the flowfield. The low number of previous studies underscores the need for more experiments (flowfield and force measurements) in conjunction with computational fluid/structural dynamics simulations that account for the complex unsteady aerodynamic effects of the flapping wing problem. Validation of modeling efforts using experimental data will serve to provide confidence in determining the effects of combinations of Reynolds number and flapping frequency resulting in high effective angles of attack incurred during a flexible wing’s motion. The next section contains a discussion of the modeling efforts in the area of flapping wings. 1.8 Review of Flapping Wing Numerical Studies A detailed understanding of the complex, unsteady, and three-dimensional aerodynamics and the complex interaction is critical to an understanding of flapping 59 MAV performance; however, the governing equations that describe these processes are non-linear and are challenging to resolve analytically. Attempts to model such unsteady loads produced on flapping wings in forward flight using classical aerodynamic analysis have been marginally successful. A few attempts have been made at estimating the forces produced on a flapping wing theoretically. Minotti worked to estimate these forces generate by a flapping wing by using empirical coefficients [143]. The assumptions were that the flow was two-dimensional and inviscid. The model correlated relatively well with experiment, but real flowfield conditions are far more difficult for numerical solvers to handle. This can be attributed to difficulties in modeling the highly nonlinear unsteady aerodynamic effects produced by flow separation, shedding of concentrated vorticity, and wing- kinematics associated with forward flight [56,144,145]. Because of this, few analytical and numerical studies have been carried out on flapping wings in forward flight using complex non-linear calculations. There have been various attempts at building semi-empirical extensions based on classical theories; however, issues with modeling turbulence and the implementation of moving, overlapped grids present challenges that require more empirical data [120,146–149]. Using modified strip theory, DeLaurier developed an aerodynamic model for large aspect ratio flapping wings in forward flight [150]. This model included the effects of the mean camber, skin friction drag, and sectional mean angle of attack. Partial leading edge suction, stall behavior, and vortex-wake effects were also addressed. Many of the assumptions were based on the fact that the wing’s aspect ratio is assumed large enough that the flow over each wing section is essentially 60 chordwise. This model was used to calculate average propulsive thrust and vertical forces, power required, and propulsive efficiency of a flapping wing in steady level flight. However, this model does not hold for low aspect ratio wings (the MAV scale). Later, Singh produced an unsteady rigid-wing model to analyze insect wing aerodynamics in hover [151]. This model was based on the indicial function and was specifically modified to handle translation and rotational circulation based on thin airfoil theory. The induced flow velocity was determined using only momentum considerations (disk theory). Thin airfoil theory addresses the non-circulatory forces and the Wagner function is used to account for the starting vortex produced during both the rotational and translational circulation. The leading edge vortex was determined by using Polhamaus leading edge suction analogy for delta wings at high angles of attack [152]. The Kussner function is used to address the shed wake and tip vortex effects. It was assumed that the inflow was uniform. This model was validated against the Robofly oil tank data at very low frequencies [22]. From this model, there was a significant difference in the experiment and predictions at pronation and supination. Computational analyses of flapping wings conducted by Lui and Kawachi, examining the aerodynamic characteristics of the hawkmoth and fruit fly wing shapes, focusing on the spanwise flow inside the core of the LEV [153]. Evaluation of the characteristics of the LEV showed that the fruit fly wing shape maintained a stable LEV throughout its translational motion, whereas for the hawkmoth wing, the LEV was shed as the wing decelerated during the downstroke of the wing motion [41]. These studies were only conducted in hovering flight. Sun et al. computationally 61 examined the aerodynamics of various insects, including the fruit fly in hover and in forward flight, reporting calculations of power requirements and effects of wing kinematics [148, 149]. While these studies assessed flight dynamics, a potential simplified model was used and the analysis was limited by the 2-D linear form of the equations of motion. Using CFD models to understand the rigid flapping-wing flowfield [148,153,154] has thus far provided only limited additional physical insight. This is mainly because of the difficulties in modeling turbulence and in preserving the three-dimensional wake vorticity to relatively old ages in terms of wing flapping cycles. These studies underscore the need for computational fluid dynamics and structural dynamics simulations that account for the complex unsteady aeroelastic effects of flexible flapping wings. Well-validated simulations using experimental data, both flowfield and force measurements, are also necessary. In the literature, much of the research has focused on experimentation or computational analysis of rigid flapping wings rather than more realistic flexible flapping wings. This is because adding aeroelastic wing deformation produced by aerodynamic and inertial loading throughout the flapping cycle adds a level of complexity. Real flapping MAV wings are highly anisotropic with varying levels of chordwise and spanwise flexibility. In order to generate sufficient forces, the wings operate at high flap frequencies and significant flap and pitch amplitudes. It has been shown that wing flexibility does improve its flight performance and modeling this aspect is important. Only a few aeroelastic solvers for flapping wings have been produced. The extreme wing flexibility and low wing weight cause large non-linear deformations. 62 Therefore, it is very important to model these flexible wings accurately, with non- linearities in both the structures and fluids. Newman presents research in the area of linear membrane aerodynamics and solution methods [155]. However, Thwaites was the first to analyze the nonlinear equation associated with the aerodynamics of the wing membrane [156]. This two-dimensional sail theory enabled a variety of experimental, theoretical, and numerical approaches to membrane aerodynamics. Later, researchers expanded the analysis to include viscosity effects, nonlinearities, and three dimensionality. A tethered moth was simulated by Smith et al. [157] using an unsteady aero- dynamic panel and finite element methods. This study highlighted the importance of including the wake in the unsteady analysis of flexible wings. The method used a combination of beam and membrane elements. However, the model did not address the separation at the leading edge. As a result, there was not good correlation of forces when compared to experiments. Singh also developed an aeroelastic model [151]. At the University of Maryland an in-house structural solver was coupled to an unsteady aerodynamic solver. The structural solver was based on linear finite element plate analysis. It was validated with prior research on rotating plates and in-house experiments. This analysis was used to simulate experiments of flexible isotropic plates undergoing insect kinematics in hover but only small flap amplitudes and low wing angles of attack were examined. The main conclusion from this study was that aerodynamic loads cannot be neglected while computing the wing response. Kim et al. [158] also used modified strip theory [150] in their model. This model addresses the issues associated 63 with dynamic stall effects at high angles of attack. A reduced multi-body dynamics model for a rectangular wing was used to account for structural flexibility. This simulation was validated with experiments at low flapping frequencies. Malhan et al. validated CFD with average force measurements for a flapping wing in forward flight [142, 145, 159]. However, detailed flowfield studies were not performed to analyze the flow state along the wing, especially for regions where unsteady flow phenomena such as vortex shedding or stall occur. During the dynamic process of the flapping wing motion where these conditions are more prevalent, the cause is not always clear. Therefore, flowfield validation is necessary to determine the effects of combinations of Reynolds number and flap frequency resulting in high effective angles of attack incurred during the wing motion. Though much computational work has given insight into the flapping wing phe- nomena [140,160], much of this work is limited by a 2-D flow environment, simplistic aerodynamic modeling (based on unsteady lifting line airfoil theory, unsteady panel methods), simplistic structural models (based on linear/non-linear beam models, linear plate/membrane models), or all of the above. Also, many of the simulation results have not been validated against experiments. 1.9 Need for New Experimental Data Some technological advances in micro manufacturing have illuminated the intricacies of natural wing kinematics and have shown that the wing motion can be reproduced on bio-mimetic vehicles. However, there is a dearth of experimental 64 data on flexible wings in forward flight. This is because the goal of most of the previous studies was to understand only the aerodynamics of flapping wings and, therefore, utilized rigid wings with prescribed kinematics. However, all the biological wings and even the expected MAV wings are flexible and undergo large deflections while flapping, due to aerodynamic and inertial forces. Even though the rigid wing experiments could help understand the aerodynamics, a realistic flapping wing should be analyzed as a coupled aeroelastic problem. Thus, these studies need to be extended to flexible wings. Specifically, most of the previous experiments were conducted using dynamically scaled mechanical wing models flapping in water or oil at extremely low frequencies (< 0.5 Hz) [20,22,88–90,101]. Such an approach was followed because the Reynolds number of a natural flyer or even an MAV could be matched at very low frequencies in a high density medium such as water or oil. This could greatly reduce the inertial force contamination of the measured fluid dynamic forces (eliminating the need for inertial force subtraction) and would also reduce the structural loads on the flapping mechanism. Also, in most of these experiments, each degree of freedom in the wing kinematics was actuated by independent servo motors and not using a continuously rotating motor with a linkage mechanism. Therefore, another reason for performing these low frequency experiments in water or oil is because of the fact that these servo motors do not have the bandwidth to generate the required wing kinematics in air at high frequencies. Even though this approach is more convenient and even more accurate, than high frequency experiments in air, this would limit such studies to rigid wings. These low frequency motions can simulate 65 high levels of fluid dynamic forces (and match the Reynolds number) owing to the higher density of oil and water (almost 800 times higher than air), the inertial forces would be several orders of magnitude smaller than the fluid dynamic forces. This is because the inertial forces are proportional to the square of frequency. For a realistic flexible wing flapping in air, the inertial forces have a stronger contribution than aerodynamic forces. Therefore, in order to produce similar deflections on a flexible wing, along with matching the Reynolds number and the reduced frequency, it is also important to match the ratio of aerodynamic to inertial forces in these dynamically scaled experiments, which is not feasible if the experiments are carried out at such ultra-low frequencies in oil or water tanks. A flexible wing must be tested in air at high frequencies to understand flapping wing aerodynamics. This has yet to be done in much of the literature and only a small number of flight viable flight capable flapping MAVs have been produced. Although a few flight-capable flapping vehicles of varying performance levels have been produced, there are still minimal experimental measurements, especially in the area of flowfield measurements. For most of these studies, there has not been any structural characterization of the wings. Currently, there is a need for experimental data on structurally characterized flexible, light-weight wings. Of particular interest, for both avian and insect flapping flight, is the creation and manipulation of flow unsteadiness when these wings are simultaneously pitched (feathering) and flapped at a particular wing beat frequency during the reciprocating motion. These parameters become very important for understanding the flow physics on the realistic flapping wings. 66 1.9.1 Chapter Summary This chapter contains a discussion of the motivation and background for MAVs. Advancements in MAV design, experiment techniques, and computational efforts are highlighted along with many challenges associated with developing MAVs. Specifically, low Reynolds number effects and the flowfield inherent to MAVs were also described. In this flow regime, there are lower maximum lift coefficients, increased boundary layer thickness, increased shear drag, and lower max L/D compared to conventional aircraft. Details of the complexity of the aerodynamic environment are explored in the areas of flowfield unsteadiness, wing-wake interaction, shed vortices, flow separation, and dynamic stall. Three dimensional effects are enhanced by the increased influence of root/tip vortices. These characteristics underscore the challenges associated with analyzing and producing efficient MAVs. In the area of experiment, the concerns associated with complex wing kinematics and the difficulties in making measurements on small low aspect-ratio wings with complex geometry/structural designs origins was presented. With modeling, the 3-D unsteady low Reynolds number flow invalidates the 2-D assumptions. In general, many reduced order models and methods become difficult to justify, especially when large non-linear deflections are the norm. In general, there is a dearth of research in the area of structurally well characterized flexible, light-weight flapping wings. The literature lacks detailed experimentation and analysis of how the flowfield affects force production mechanisms, especially in the context of biologically-inspired flexible flapping wing MAVs in forward flight. 67 This chapter is concluded with the research objectives and an outline of the dissertation. 1.10 Research Objective and Approach The goal of this study was to assess the lift and drag of a flexible MAV-scale flapping wing and associate them with the flowfield response. In-house motion tracking, flowfield, and airloads experiments were executed in an open-circuit wind tunnel on rigid and flexible flapping wings. The flexible wings were structurally characterized by performing experimental testing and these structural properties and test data were used to compare with an aeroelastic solver. The primary focus was to examine the dynamics of a flapping flexible wing and to understand the contribution of flapping kinematics along with wing deflections (primarily dynamic twist) to the vertical and horizontal aerodynamic forces during a flap cycle. The objectives of this research are listed below: 1. Design and characterize the structural parameters of a flexible wing 2. Develop effective flow diagnostic techniques to permit high-fidelity measure- ments in the flowfield of a bio-mimetic flapping wing 3. Analyze flowfield in forward flight conditions (flexible and rigid wings) 4. Evaluate flight parameters in forward flight conditions 5. Apply aeroelastic model to study the flowfield response to wing flexibility 68 1.11 Organization of the Dissertation In the present work, the temporal evolution of the flow features associated with rigid and flexible flapping wings in forward flight were quantified and the predictive capability of a computational flow solver coupled with a structural analysis tool was compared using flowfield measurements (using PIV) and direct force measurements. Rigid wing experiments are used as a baseline. Force and PIV experiments are performed using pure flapping kinematics. Due to the torsional compliance of the wing, a spanwise variation in pitch angle (twist) was produced, in addition to the temporal variation (dynamic twisting). The parameters defining the wing kinematics are wing beat frequency and root pitch angle variations, which were examined systematically. Pertinent flowfield results are presented with experimental and computational aerodynamic forces. Specifically, analysis of the LEV, downstream wake, and airloads associated with the wing are examined using the CFD-CSD model and experiment results. Utilizing non-intrusive flow and force measurement techniques along with an aeroelastic solver to understand the influence of wing flexibility on force production and the flow physics of a flapping wing MAV in forward flight, • Chapter 1 contains information about the motivation for MAV development and highlights the pertinent past research contributions. • Chapter 2 describes the research methodology and approach. Specifically, the areas of discussion are the experimental and analytical methods. 69 • Chapter 3 shows a validation/comparison of the aerodynamic and structural solver with experiments. • Chapter 4 contains the findings related to the flapping wing fundamental flow physics. • Chapter 5 summarizes the key conclusions from this work. • Chapter 6 presents the major contributions of this dissertation research con- ducted at the University of Maryland and suggests potential follow-up research. 70 Chapter 2: Methodology: Numerical Approach and Experi- mental Setup This chapter contains a discussion of the approach and methodology associ- ated with the numerical aeroelastic solver used in conjunction with particle image velocimetry, force measurements (airloads analysis), and motion tracking (for struc- tural properties). Wind tunnel specifications, experiment challenges, and experiment uncertainty are also presented. 2.1 Numerical Approach An in-house CFD solver, OVERTURNS, is used to model the aerodynamics of a flapping wing. This solver is a Reynolds Averaged Navier Stokes Solver and has been used extensively for rotary wing applications, but has been extended to low Reynolds number flapping wing flows. The focus is to simulate the flowfield of flapping wing MAVs. However, only one wing is modeled in this study and its surface can be treated as a solid wall. The aerodynamic solver is coupled to a multibody dynamics solver to model the structural dynamics. The numerical solution algorithms are described in this section. The flow domain that is being studied is initially identified 71 and the details of the mesh system, aerodynamic solver, structural solver, and the connectivity approach are discussed. 2.1.1 Description of Flow Solver The 3-D fluid dynamics of the flapping wing were computed using a compressible, structured, overset Reynolds Averaged Navier-Stokes (RANS) solver, OVERTURNS [49]. This overset structured mesh solver utilizes the diagonal form of the implicit approximate factorization method developed by Pulliam and Chaussee [161] with second-order accuracy in time. The advantage of using an overset structured grid is that different grids can be generated independent of each other and can be placed in the region of interest without any distortion [49]. Due to these advantages, the current work employs overset meshes. To calculate the inviscid terms, a third-order Monotone Upstream-Centered Scheme for Conservation Laws [162] with Roe flux difference splitting [163] and Korens limiter [164] is used. A second-order central differencing method is implemented in order to calculate the viscous terms of the compressible RANS equations. Given the relatively low Mach numbers of the flapping wing cases, a time accurate, low-Mach pre-conditioning dual time-stepping scheme described by Buelow et al. [165] and Pandya et al. [166] is applied. The low-Mach preconditioning not only improves the convergence of the simulations, but also the solution accuracy [49]. The Spalart-Allmaras [167] turbulence model is employed for closure of the RANS equations. The equation is given by: ∂ν¯ ∂t + V · (∇ν¯) = 1 σ [∇ · ((ν¯ + ν)∇ν¯ + cb2(∇ν) 2] + cb1S¯ν¯ − cw1fw [ ν¯ d ]2 (2.1) 72 These equations relate the Reynolds stresses to the mean strain. The constants in this equation are cb1 = 0.1355, cb2 = 0.622, cw1 = cb1 κ + 1+cb2 σ , κ = 0.41, σ = 2 3 , and cw2 = 0.3, cw3 = 2.0. The turbulent eddy viscosity, νt, is obtained by solving this partial differential equation for a related variable, where the two quantities are related by νt = ν¯fv1. Here, fv1 is a function of ν¯ and the molecular viscosity, ν. The equation for fv1 is fv1 = χ3 χ3 + c3v1 (2.2) with χ = ν¯ν , cv1 = 7.1, fv2 = 1 − χ 1+χfv1 , fw = g [ 1+c6w3 g6+c6w3 ]1/6 . In the equation, g = r + cw2(r6 − r) and r = min( ν¯S¯κ2d2 , 10.0). S¯ is equal to S + ν¯ κ2d2fv2 and S is |ω|. The distance from the wall is d, and V is the mean flow velocity [167]. Closure for all variables is provided by loosely coupling the Navier-Stokes equations and obtaining the turbulent eddy viscosity. From this, the shear stress in the moment and energy equations can be evaluated. The turbulence model has the advantages of ease of implementation, computational efficiency, and numerical stability. An implicit hole cutting method developed by Lee [168], and further refined by Lakshminarayan [49], is employed to find the connectivity information between the overset meshes used in the simulations. Previous studies [169–171] suggest that a Large-Eddy Simulation (LES) may more accurately predict the highly vortical and separated flow present in low-Reynolds number MAV-scale CFD problems. However, in comparison to RANS, the LES approach would yield a penalty of higher computational costs and increased grid 73 refinement. Validation of the current solver against the LES results shown by Gordnier et al. showed good agreement [171]. Therefore, RANS was used for all the simulation results presented in this work. 2.1.2 Description of Grid System and Simulation An avian wing has very complex wing motions that are a combination of flapping, pitching, and wing folding/bending. The two main degrees of freedom, flapping and pitching, are considered for this study. The wing planform used during simulation was based on the rectangular, low-aspect ratio wing used during experimentation. In the computational model, the wing has an aspect ratio of 1.67 and a flat-plate cross section consisting of a thickness-to-chord ratio of 1.0%. In the simulation, the wing is offset from the flapping axis by 0.33 chords to account for the 20% root cutout in the experimental setup. An overset system of meshes, consisting of a C-O type airfoil mesh for wing and Cartesian background mesh is used for the computation. The wing mesh has 307×161×75 grid points in the wrap-around, spanwise and normal directions, respectively. To better capture the shed wake and vorticity at these locations, points were more tightly grouped toward the leading edge, trailing edge, wing root, and wing tip. All simulations were run using at least 6 sub-iterations per time-step to minimize linearization errors with farfield boundary conditions being implemented 20 chord lengths away from the wing surface. When the grid has been generated and grid motions are prescribed, the flowfield properties at each grid point within the overset mesh system can be obtained by solving the 74 Figure 2.1: Wing mesh embedded in background Cartesian mesh. conservation laws of physics for fluid flow. An example of the CFD grid can be seen in Fig. 2.1. 2.1.3 Grid Sensitivity Study A grid sensitivity study was executed for the rigid wing simulation. All computational grid sensitivity simulations were carried out on the flat plate airfoil described in the previous section. In order to determine the appropriate simulation grid size, a single flapping wing case with second-order spatial discretization was performed on three different C-O meshes. The nominal mesh consisted of 267×141×85 nodes in the wrap-around, spanwise and normal directions. A coarse mesh and fine mesh consisted, respectively, of 133×113×65 and 397×169×105 in the wrap-around, spanwise and normal directions. Each mesh was used to simulate the 4 Hz flapping 75 case at a pitch angle of 0◦. A total of 1440 time-steps per cycle were carried out for a duration of 4 flap cycles. Figure 2.2 shows a comparison of the predicted aerodynamic force coefficients for the three grids tested during the 4th flap cycle. It can be seen that the aerodynamic coefficient values obtained from the coarse, nominal, and fine meshes show similar agreement. Note that the results from the coarse mesh show greater oscillations in the vertical force time-history near instances of peak aerodynamic force production. This is because the coarse mesh does not have enough nodes to completely resolve the leading edge vortex formed. Consequently the nominal mesh was used to obtain all the results discussed in the subsequent sections to resolve the flow features without increasing computation time. 2.1.4 Structural Model and Coupling The structural analysis portion of the numerical simulations was conducted using MBDyn, which is an open-source general purpose multibody software developed at Politecnico di Milano [172, 173]. MBDyn integrates the equations of motion temporally for a constrained mechanical systems (in this case for flapping wing). The equations associated with the linear and angular motion of a set of independent rigid bodies (Newton-Euler equations) are explicitly constrained by a set of algebraic equations. These equations express the kinematic constraints between the bodies. Modifications were made to enable the use of non-linear 2-D structural elements (shells) and allow for generic Fluid-Structure Interaction functionality. This solver’s implementation of non-linear beam elements, in conjunction with the newly added 76 (a) Vertical lift coefficient, CZ. (b) Horizontal lift coefficient, CX. Figure 2.2: Variation of predicted aerodynamic force coefficients for different mesh resolutions (k = 0.32, θ = 0◦). 77 non-linear shell elements [174], gives it the capability to accurately model the large, non-linear deflections and complex structural dynamics of highly flexible wings [141]. The structural analysis is based on the direct integration in time of an Initial Value Problem that describes the dynamics of a set of arbitrarily interconnected bodies. Compliant connectivity is realized in the form of non-linear finite elements, using the aforementioned beam and shell elements, while the kinematic constraints are modeled by algebraic relationships enforced using Lagrange multipliers. In order to facilitate the data exchange between the CFD solver (OVERTURNS) and the CSD solver (MBDyn), a Python computational framework was implemented. This coupling algorithm allows for the exchange of structural deformations (kinematic mapping from the CSD to the CFD solver) and aerodynamic loads (from the CFD to the CSD solver) at the fluid-structure interface (i.e., the wing surface). Tight- coupling is used for the data exchange, meaning that information is passed during each sub-iteration of the non-linear problem solution procedure for a given time step. However, in an effort to reduce computational cost, it is relaxed as soon as it is allowed by the relative compliance of the coupled problem. While coupling of the two solvers was done using both a first order and staggered second order accurate approach, which satisfies the discrete geometric conservation law, it was observed that both approaches gave similar results. It should be noted that the domains and meshes of the two solvers are non-conformal. These incompatible domains were mapped via an original scheme that preserves the work exchanged between the structural and the aerodynamic domains. This mapping is based on a Moving Least Squares fitting of the discretization of the interface between the two domains using a 78 compact support that consists of Radial Basis Functions. The full simulation results were compared to flowfield and force experiment data. 2.2 Experimentation This section contains a discussion of flapping wing mechanism design, rigid and flexible wing designs, and the high-speed flow measurement techniques and force measurements. The main experiment technique for measuring the flow was time-resolved digital particle image velocimetry (DPIV). PIV was used to acquire velocity fields, and the components of vertical (lift) and propulsive or drag force were measured using a force transducer. These flow interrogation techniques allowed for the quantification of the unsteady mechanisms. 2.2.1 Flapping Mechanism and Wing Design A flapping mechanism was used to provide the required flapping kinematics to the wing (see Fig. 2.3(a)). A schematic of the four-bar linkages is shown in Fig. 2.3(b). The link AD is the driving link and BC is connected to the wing. If the link AD rotates by θ, the driven link rotates by φ (wing flap angle) which is calculated as follows: x = L2 cos θ − L1 (2.3) y = L2 sin θ (2.4) 79 φ = tan−1 y x − cos−1 L24 + x 2 + y2 − L23 2L4 √ x2 + y2 (2.5) For the mechanism, the lengths of the links were chosen to closely replicate sinusoidal motion. The baseline linkage lengths were L1 = 0.42 cm, L2 = 5.52 cm, L3 = 6.41 cm, and L4 = 1.6 cm. The wing would be attached to L4. Flapping kinematics were designed to be close to harmonic for the flapping amplitude of 40◦. The four-bar was driven by a flywheel attached to the shaft of a motor. A 260-W outrunner motor manufactured by MODEL MOTORS (Czech Repub- lic; AXI 2217) was used to rotate the flywheel. A 3:1 planetary gear box was used to provide the required torque. A potentiometer was used to measure the instantaneous angular position of the wing as it flaps. The flapping frequency was measured from the 1/rev signal obtained using a Hall switch and was verified using a digital laser tachometer, stroboscope, and the PIV system timing hub. The flapping-wing mechanism is fixed to the bottom of the wind tunnel such that the flow is parallel to the wing plane when the angle of attack is zero. A wing attachment was designed such that the wing can be set at any desired geometric angle of attack. Only pure flap tests were executed (no active wing pitching). This means that the wing is held fixed at the root and is flapped at a set initial pitch angle. As the wing flaps, there is a spanwise variation in pitch at one instant of time for the flexible wing, so the wing sections rotate by different amounts along the span due to torsion flexibility. Also, the spanwise twist varies as a function of time. This is referred to as dynamic twisting and is caused by the inertial and aerodynamic forces acting on the wing 80 during flapping. Rigid wings were fabricated using unidirectional carbon rods, which are 1.5 mm in diameter as shown in Fig. 2.4(a). High-speed videos of these wings while flapping showed that there was no discernible deformation. These wings had a rectangular planform with a span of 12.7 cm and a chord of 7.62 cm with a 20% (3.05 cm) root cutout. For PIV experiments, a thin layer of matte black paint was used to reduce reflections. The wing mass was 8.2 grams. Similar to the rigid wings, torsionally flexible wings were also fabricated using unidirectional carbon rods, which are 1.5 mm in diameter as shown in Fig. 2.4(b). The chordwise ribs (Root, Mid, Tip) were constructed from unidirectional single layered carbon fiber strips. These wings also had a rectangular planform with a span of 12.7 cm and a chord of 7.62 cm with a 20% (3.05 cm) root cutout. The wings were rigidly attached to the four-bar mechanism at the root tab. For PIV experiments, the wings were covered with heat resistant film to prevent laser burn, and a thin layer of matte black paint was used to reduce reflections. The wing mass was 6.47 grams. A six-component balance was used to measure the flapping loads. This load cell, the ATI Nano 17, weighs less than 10 grams and was attached to the root of the wing to measure the loads acting on the wing directly without any interference from the inertial loads generated by the rest of the mechanism. The first natural frequency of this balance is 7200 Hz. The force transducer was connected to a National Instruments (Austin, TX) USB DAQ device (NI USB-6251). Raw data from the sensors were transmitted to a computer using a software called LabVIEW. These 81 (a) Flapping mechanism. (b) Four-bar linkages. Figure 2.3: Four-bar linkages and flapping mechanism. 82 (a) Rigid wing. Leading edge spar Ro ot r ib Mid rib Tip rib 5 in 3 in (b) Flexible wing. Figure 2.4: a) Rigid and b) flexible wing construction. 83 voltages were later post-processed with the ATI-supplied calibration matrices to generate the observed forces and moments. Since instantaneous force measurements would include significant inertial forces, only the time-averaged force measurements were used for validating the analysis. 2.2.2 Particle Image Velocimetry Setup PIV is a measurement technique that is non-intrusive and gives a planar velocity field within a desired region of interest. This technique is especially useful for MAV applications where the scale of the flow is such that hot wire anemometers or wing pressure taps are not feasible. In general, uniformly dispersed tracer particles are introduced to the fluid medium to track the flow. The three main components of the digital PIV system are the high repetition pulsed laser (with optics), PIV camera, and data acquisition system (including a sychronizer). The laser produces a planar laser light sheet. Particles within the region of interest are illuminated by the laser sheet. The digital camera is orthogonal to the light sheet. The camera is uniformly focused over the region of interest. The seeded flowfield is then illuminated twice within a small pulse separation time, ∆t = t2 − t1, which is on the order of microseconds. The camera is then digitally synchronized with the laser to capture two successive images (frame A and frame B). After a time series of image pairs is obtained, the raw data is hyper-streamed to a processing computer. For the selected pulse separation time, each image pair contains the pixel displacement of the particles in that region of interest. During processing, the images 84 are divided into smaller interrogation windows. A cross-correlation technique is used to match the intensities of the particles between the two frames within each small window, and then computes the magnitude and direction of the pixel displacement (∆x,∆y) [175]. This process was performed utilizing a Fast Fourier Transform (FFT) algorithm in the frequency domain [176]. A deformation grid algorithm is used to improve particle correlation [177]. Once the pixel displacement was found, the local velocity could be calculated using the equation { U V } = 1 M∆t { ∆x ∆y } (2.6) where M is the image magnification. The image magnification is calculated within the PIV software from a master image with a calibration plate placed in the laser plane to calibrate the ratio of pixels to millimeters. The cross-correlation process was then repeated for all interrogation windows within the image, the outcome being a planar velocity field. The interrogation windows were sized so there were at least 5-10 seed particles within each and the particles moved only about one quarter of the window length. One of the assumption associated with the PIV cross-correlation technique is that the particles within each interrogation window translate linearly between image pairs. With this in mind, the interrogation window size, as well as the value of ∆t, was chosen to best resolve the relatively steep velocity gradients induced by the flapping wing. The PIV system was acquired from LaVision Inc and Fig. 2.5 outlines the PIV process methodology. 85 Figure 2.5: Schematic depicting the general setup of a 2-D PIV system and image processing methodology. A schematic of the PIV test setup is shown in Fig. 2.6. A low-speed, open- circuit/closed-section wind tunnel was used to conduct the experiments. The test section cross sectional dimensions are 20 in × 28 in, and the tunnel’s maximum speed was 100 mph. At the tunnel inlet, the flow was seeded using a six-jet atomizer allowing for sub-micron sized tracer particles (<0.2 µm) to be introduced into the wind tunnel inlet. Bis (2-ethylhexyl) sebacate was used as the seeding fluid. PIV images were acquired in chordwise planes for the span locations y/b = [0.3 0.4 0.5 0.6 0.7 0.8 0.9] using the time-resolved measurement system. A double-pulsed high-speed laser illuminated the seeding particles (Litron LDY304 Nd:YLF laser, wavelength of 527 nm with 30 mJ/pulse at 1 kHz rep-rate). Two mirrors were used to redirect portions of the laser sheet back towards the flapping wing to reduce the shadow cast 86 Figure 2.6: Schematic PIV experiment setup in wing tunnel test section. by the wing. A high-speed camera (Phantom V311, 4MPx, 3,250 fps) was positioned above the test section and oriented downward to view the wing through a circular acrylic window. The actual tunnel test section is shown in Fig. 2.7 and the camera view is from the top of the wind tunnel. The flapping mechanism was mounted vertically from the bottom of the test section with the wing cantilevered. To avoid any obstruction to the flow, the flapping mechanism was installed outside the test section with only the wing exposed to the flow. The wing pitching mechanism was located at the root, permitting an unobstructed field of view for the camera. PIV experiments were conducted with the wing flapping at the pitch angles θ0 = 0◦− 24◦ in increments of 3◦, for 0, 4, 6, 8, and 10 Hz for a freestream velocity of V∞ = 3 m/s. The PIV images contained the pixel displacement of the seed particles for the selected pulse separation time. The seeding density was adjusted to provide at least 10 particles in each interrogation window [178]. The raw PIV images were 87 Figure 2.7: Actual PIV experiment setup in wing tunnel test section. then divided into smaller 48-by-48 pixel interrogation windows (with a 50% overlap) and a deformation grid correlation method was used for processing [177]. A cross- correlation FFT based procedure was used to match the intensities of the particles between the two images within each small window, providing the magnitude and direction of the pixel displacements [179]. This was carried out with interrogation windows of size 24-by-24 for better resolution of the higher velocity gradients within the flow. Spurious vectors in the results were determined using a Gaussian peak method with a signal-to-noise ratio of 1.3. The data was sent through a local vector validation process by examining a neighboring set of 3-by-3 vectors using a universal median test. Vectors identified as spurious were then removed. PIV data (e.g., 50 image pairs corresponding to the downstroke and upstroke respectively) were phase-averaged from each experimental case. Concerning alignment, great care was taken such that the PIV camera was 88 oriented at the center of the wind tunnel test section and the camera viewing plane was perpendicular to the wing span at the mid-stroke location. A PIV mask of the static wing was created within the frame. This mask was used as an indicator of the wing chord at mid-downstroke and mid-upstroke during flapping. Time-resolved PIV data was acquired during the wing flapping motion at 1300 frames per second and only the image pairs that fit the custom mask (produced when the wing was static) were used to represent the mid-upstroke and mid-downstroke positions. 2.2.3 Wind Tunnel Specifications The wind tunnel experiments were performed using an open circuit, low- turbulence wind tunnel located at the University of Maryland (schematic shown in Fig. 2.8). The main components of the wind tunnel are the screen, nozzle (contraction section), test section, diffuser, and the fan. Downstream of the test section is a diffuser that slows the flow as it approaches the impeller. The impeller is driven by a variable speed electric motor. The tunnel inlet is where seed particles were passed through a honeycomb structure, which aided in removing much of the larger turbulence and eddies that would enter the test section. The wind’s exit was released into the open air. DPIV was used to measure the background turbulence within the wind tunnel at the center of the tunnel at a height of 10 in (tunnel center), 5 in, and 2 in from the tunnel floor of the test section. An interrogation window size of 2 in×2 in was utilized to quantify the freestream fluctuations. Ten seconds of data was obtained at a ∆t of 50µs. Though the target freestream velocity was 3 m/s, 89 Figure 2.8: Schematic of low speed wind tunnel. the tunnel velocity was set at 1, 3, 6, 9, and 23 m/s. The turbulence intensity in the test section has been determined by dividing the standard deviation of the wind tunnel speed by the average freestream velocity [180]. Figure 2.9 shows the percent turbulence for various locations in the wind tunnel with varying tunnel speed. The higher levels of turbulence were seen at a height of 2 in, which is close to the tunnel wall. Close to the wall, as the tunnel speed is increased, the turbulence increases as well (due to wall effects). The wind tunnel fluctuations remained below 2% in all cases. A Fast Fourier Transform was applied to the freestream velocity, V∞, data and there was no discernible frequency peak. Additionally, the velocity fluctuation signal magnitude was extremely low. Also, the low intensity of the fluctuations was such that it did not have any effects on the flow structures because wing flap speeds and flow feature velocities were several orders of magnitude greater than the tunnel fluctuations. Therefore, tunnel fluctuations are negligible. The two main parameters associated with flapping flight are the reduced frequency and Strouhaul number. The reduced frequency is a measure of flow 90 Figure 2.9: Variation in wind tunnel turbulence with wind tunnel speed. unsteadiness. It is the ratio between forward velocity and flapping velocity. This is expressed in an equation as follows: k = Ωc 2V∞ (2.7) where Ω is the flapping frequency, c is the wing chord, and V∞ is the freestream velocity. Flow is considered unsteady if k > 0, but can generally be considered quasi-steady for a reduced frequency 0 < k < 0.03, where unsteady effects are not significant [7]. For 0.03 < k < 0.1, flow can be considered moderately unsteady, and beyond k = 0.1, flow is considered fully unsteady. Reduced frequency is typically between 1 and 10 for small insects [7]. The reduced frequencies, k, for this study corresponding to flap frequencies, Ω, are 0, 0.32, 0.48, 0.64 and 0.80 respectively. 91 Strouhal number describes the oscillating flow mechanism and is also utilized to characterize flapping in forward flight. Strouhal is related to reduced frequency by St = kA pic (2.8) where A is the flap amplitude. This data set is applicable to low Reynolds number flyers where the Strouhal number is within the range of 0.2 < St < 0.4 [181]. Conversely, the advance ratio is given by J = V∞ 2piΩA (2.9) At first glance, this dimensionless parameter may appear to be the inverse of the reduced frequency. Ho et al. claimed that this parameter accounts for the three- dimensional nature of the flow [138], where the breakpoint between quasi-steady and unsteady flow occurs when J = 1. For J > 1, the flow can be considered quasi-steady while J < 1 corresponds to unsteady flow regimes. Most insects operate in this unsteady regime. This experiment set explores the set of advance ratios J = [2.4317, 1.6207, 1.2154, 0.9722]. The main experiment parameters for this dissertation are listed in Table 2.1. 2.2.4 PIV Experiment Challenges Many of the challenges in performing PIV measurements on flapping wings resides in the ability to view the flow associated with 3-D unsteady low Reynolds number aerodynamic phenomena occurring on a small low aspect ratio wing under- 92 Table 2.1: Time-resolved PIV test matrix. Freestream, V∞ 3 m/s Reference Reynolds number, Re 15,000 Flap amplitude ±40◦ Span location 30% – 90% Initial Pitch Angle, θ0 0◦ – 24◦, ∆θ = 3◦ Flap frequency, Hz 0, 4, 6, 8, 10 Reduced frequency, k 0, 0.32, 0.48, 0.64, 0.80 Strouhal number, St 0, 0.13, 0.2, 0.26, 0.33 going large non-linear deformations with complex wing kinematics. Specifically, it is difficult to acquire the optimum laser placement because every wing position and phase angle gives at least some surface reflections from the wing. These spurious reflections can, in some cases, obscure the flow features or otherwise limit the ability to track the flow features in time because they exist for a very short duration of time. Wing masking and background image subtraction were found to be useful techniques for removing unwanted reflections before processing the PIV images. In particular, intense reflections from the wing planform made the cross-correlations near the wing surface extremely difficult. These reflections were slightly reduced by the appropriate selection of camera aperture setting and laser intensity level. Masking was used to remove areas in the region of interest that were not desired for processing and background subtraction involved generating a composite minimum intensity image, which was then subtracted from each individual PIV image. 2.2.5 Experimental Uncertainty The preset wing angle of attack and spanwise position are accurate to within 0.25◦ and 1 mm respectively. Again, the flapping frequency was verified using a digital 93 laser tachometer, a stroboscope, and the PIV system timing hub. The two main sources of error associated with DPIV are the mean bias error (bias) and root-mean- square (RMS) error (RMS) where the total error (tot) is equal to the sum of bias and RMS. The fundamental source of bias and RMS error is the implementation of cross correlation. Other contributors include the sub-pixel accuracy of the correlation peak-finding scheme and noise within the particle images themselves. Both the mean-bias and RMS errors have been found to be on the order of 0.1 pixel for cross correlation and the error associated with particle image pattern matching are an order of magnitude less than those of cross correlation [182]. For the current PIV setup, bias = −0.01 due to image pixel loss. The RMS error, produced from the calibration, was kept below -0.3 pixels. Therefore, the total error associated with the PIV measurements is tot = −0.04 pixels. Based on an average displacement of 6 pixels, the estimated uncertainty of the measured velocity is 0.667%. There is also a confidence level associated with the motion tracking measurements. Given the size of the retro-reflective markers (3.5 mm in diameter) and the relative proximity of the motion capture camera system to the flapping wing or wing frame, it was estimated that the error in the spatial tracking of the markers was roughly 2.8% of the measured displacements. With the PIV experiments, the LEV and trailing edge vortices were imaged at strategic instances during the wing stroke. This was compared to the predicted flowfield, which produced the LEV formation, convection, and shedding. Force measurements quantified the effect of the LEV on the lift and drag. The results obtained contributes to the understanding of flapping-wing flow physics. From these 94 experiments, needed data is added for the low Reynolds number realm. The next section contains a discussion of how the structural properties for the flexible wing were characterized. 2.3 Structural Properties Characterization of the flexible wing frame structural properties was primarily conducted using the VICON motion tracking system. Tensile tests were used to obtain the properties for the wing covering. In-house VICON motion tracking experiments (static and dynamic) were done to generate more data for the structural model validation. The VICON system allows for the tracking of an object in three-dimensional space. This is done by applying retro-reflective reference markers to a surface and placing high-speed infrared cameras in various positions and orientations around the object. The camera system then tracks the markers in space and time by utilizing triangulation methods. The camera system must be placed such that at least two of the cameras can view each marker throughout the duration of the experiment. A VICON camera system has an operable capture frame rate in the range of 20 to 2000 frames-per-second (fps) depending on the hardware platform, camera choice, and number of cameras being used. For this study, four cameras were strategically placed around the test area at different heights and angles to track the markers. Half sphere retro-reflective markers were placed on the wing. Each marker’s diameter was millimeters. Figure 2.10 shows the VICON test setup and the camera placements. 95 Figure 2.10: Typical test setup for motion tracking system. During motion tracking experiments, it is very important to ensure that the infrared cameras have an unobstructed views of the markers. The total mass of each marker was 0.02 g. After the raw data were captured, marker reconstruction and processing provided the complete displacement of the wing components during static tests and the time-history of the three-dimensional location of each marker during dynamic tests. The structural properties of each of the wing’s members were obtained experi- mentally, and then they were used to build a structural model in the Multi-Body Dynamics solver. Bending and torsion tests were carried out for the chordwise ribs and the leading edge spar. Loads were added to simulate bending and torsion and the wing member deflections were measured using four infrared cameras. The bending and torsion stiffness properties (EI and GJ) were then obtained from these results. Fig. 2.11 shows the marker placement and setup. The determination of the bending 96 and torsion stiffness parameters in this work are based on Euler Bernoulli Beam theory. Figure 2.11(a) shows a schematic of the marker and representative load place- ment to measure the bending of the tip member. The retro-reflective markers are shown on the schematic and were placed at the leading edge, midpoint, and tip of the member. To calculate the bending stiffness from the measured, nonlinear deflections, an equation for bending is given using a homotopy analysis method with explicit formulas for calculating large deflections of a beam [24]. A schematic of the large deformation of a cantilever beam under point load at the free tip is shown in Fig. 2.12. The bending equation of a uniform cross-section beam with large deformation is expressed as follows: ds dθ = P EI (l1 − x), θ(0), θ ′(l) = 0 (2.10) where the arc-coordinate of the neutral axis is s, the horizontal coordinate from the fixed end is x, the beam length is l, the point load at the free tip is P . The bending stiffness is EI, the rotation of cross-section of the beam is θ, and the horizontal distance between the two ends is l1. For the cantilevered member, the axial elongation is smaller than the deflection at the tip. Differentiating the equation with respect to s and then using the dimensionless variables ξ = s/l, the equation 2.10 becomes θ′′ + α cos θ θ(0), θ′(1) = 0 (2.11) 97 (a) Isolated bending. (b) Isolated torsion. Figure 2.11: Schematic of VICON test for isolated bending and torsion testing. 98 Figure 2.12: Large deformation of a beam under point load at free tip [24]. where α = Pl 2 EI . The rotation angle of the cross-section at the tip is θB = θ(1), and the exact vertical displacement of the free tip is given by the equation fB = l − 2[E(µ)− E(φ, µ)] √ EI P (2.12) where E(µ) and E(φ, µ) are the elliptic integrals of the second kind shown below: µ = √ 1 + sin θB 2 φ = arcsin ( 1 √ 2µ ) (2.13) With this, the dimensionless vertical displacement at the free tip is shown below fB l = 1− 2 √ α [E(µ)− E(φ, µ)] (2.14) Additionally, l1 l = √ 2EI P l2 sin θB = √ 2 sin θB α (2.15) The dimensionless horizontal displacement of the free tip is expressed by the following 99 equation: δB l = l − l1 l = 1− √ 2 sin θB α (2.16) If there are large deflections, √ α = K(µ)−F (φ, µ) where K(µ) and F (φ, µ) are the elliptic integrals of the first kind. The solution is obtained using the homotopy analysis method where for large deflections, θB = α 2 f(α) g(α) (2.17) and f(α) =1 + 3.98575× 10−2α2 − 5.41174× 10−2α4 + 5.72575× 10−3α6 + 3.79533× 10−4α8 − 8.87896× 10−6α10 + 2.63041× 10−8α12 − 1.51429× 10−11α14 − 2.29142× 10−15α16 − 3.45006× 10−21α18 − 7.00678× 10−28α20 (2.18) g(α) =1 + 0.131524α2 − 5.99231× 10−2α4 + 2.34466× 10−3α6 + 9.90299× 10−4α8 − 1.37001× 10−5α10 + 3.44172× 10−8α12 − 1.45098× 10−11α14 − 2.26721× 10−15α16 − 3.50731× 10−21α18 − 7.00678× 10−28α20 (2.19) 100 This analysis was used to predict the bending stiffness, EI, for each flexible wing member from the measured deflection data. Figure 2.11(b) shows a schematic of the test setup of a single rib with the load applied on the wing frame for a purely torsional loading. To increase the moment arm, a carbon beam is put on the tip. One load is applied from the top with a pulley while the other hangs freely. To be sure that the load applied is not creating bending, in addition to torsion, the vertical displacement of the beam tip was measured to ensure minimal deflection. Specifically, the vertical displacement of Marker 3 was verified to ensure that Markers 4 and 5 received the same load. Using this method, the shear center is the beam center (Marker 3). The torsion angle is obtained by using the displacement of the ends of the carbon rod. Figure 2.13 shows a schematic of a beam/rib in torsion where θ is the angle of twist, T is the applied torque, L is the length of the beam, J is the moment of inertia, and G is the modulus of rigidity. The assumption is that GJ is constant. The torsion equation is given below. θ′ = T (x) GJ (2.20) 2.3.1 Wing Covering Material Tensile Tests An important area of materials mechanics is the behavior of deformable bodies, especially the details of how a material will behave under load. Tensile and compres- sion experiments allow for the determination of these properties. With tensile tests, 101 Figure 2.13: Torsion test schematic. a load is applied along the longitudinal axis of a test specimen. The specimen will stretch along this axis, and the amount of elongation can be measured. Loads and deformations may be converted to stresses and strains by σ = P A  = δ L (2.21) where σ, is normal stress on a plane perpendicular to the longitudinal axis of the specimen, P is the applied load, A is the cross sectional area,  is the normal strain in the longitudinal direction, δ is the change in the specimen length, and L is the original length. The stress-strain curve gives a direct indication of the material properties. The initial portion of the stress-strain diagram for most materials is a straight line. For the initial portion of the diagram, the stress (σ) is directly proportional to the strain . Therefore, the slope of the straight-line portion is the (Elastic Modulus or Young’s Modulus) and is given by the following: E = σ  (2.22) 102 Tensile tests were performed to obtain the modulus of the wing covering. A hydraulic materials testing system (MTS) was used. Three layered rectangular specimens were clamped at each end and placed in the tensile test machine. The sample was then stretched so that the force and displacement could be measured using the MTS DAQ and a 100 lbf load cell. A 4th order Butterworth low pass filter with a cutoff frequency of 50 Hz was applied to the raw data to eliminate all high frequency noise sources. The test samples included a 10-layer Mylar membrane, a 10-layer heat resistant material, and a 5-layer heat resistant material. Table 2.2 shows data from the tests. Of particular interest is the stress and strain in the elastic region as it was not expected that the aerodynamic forces would exceed the yield stress of the material. The two wing cover materials were Mylar and a heat resistant material. Fig- ure 2.15 shows representative data from the tensile tests. Figure 2.15(a) shows raw and filtered force and displacement data for the heat resistant material in the elastic region. Here, the stress applied to the 10-layer heat resistant material sample is proportional to the strain and the material will return to its original shape when unloaded. Figure 2.15(b) shows the results of 5 separate tests for the 5-layered heat resistant material. The five specimens were identical and were each stretched to different maximum (0.05–0.30 in) to observe the elastic limit, which is the stress that can be applied without resulting in permanent deformation when unlaoded. It is seen that the 5-layer heat resistant sample does return to the original shape up to a yield stress of 3,250 lbfs2 ( = 0.02). The structural properties of each of the flexible wing members was obtained 103 Figure 2.14: Image of wing membrane material in during tensile test. (a) Unfiltered and filtered load and displacement for 10 layer heat resistant wing covering . (b) Stress and strain 5 layer heat resistant wing covering. Figure 2.15: Tensile tests for wing covering. 104 Table 2.2: Tensile test results. Test specimen Thickness (in) Length (in) Width (in) Weight (g/in) E (lbf/in2) Mylar-10 0.0105 4.863 1.532 0.2504 1.1889e+4 Heat Resistant-5 0.012 7.113 1.047 0.2531 1.6429e+5 Heat Resistant-10 0.025 7.41 1.072 0.5647 1.6168e+5 Table 2.3: Flexible wing properties. Member Properties Wing chord (m) 0.0762 Wing span (m) 0.127 Leading edge spar, GJ (N/m2) 6.1 ×10−3 Leading edge spar, EI (N/m2) 2.4 ×10−2 Root rib, EI (N/m2) 4.92 ×10−4 Mid rib, EI (N/m2) 2.4 ×10−4 Tip rib, EI (N/m2) 2.1 ×10−4 Root Rib, GJ (N/m2) 2.05 ×10−4 Mid Rib, GJ (N/m2) 4.86 ×10−5 Tip Rib, GJ (N/m2) 2.02 ×10−4 Membrane, E (N/m2) 5.02 ×109 Membrane, ν 0.45 experimentally. Bending and torsion tests were carried out for the chorwise ribs and the leading edge spar. Loads were added to simulate bending and torsion and a motion tracking system was used to measure the deflections. The bending and torsion stiffness properties (EI and GJ) were then obtained from these results. The modulus of elasticity was determined for the wing covering using tensile tests. The geometric and structural properties of the flexible wing are shown in Table 2.3. These properties were used as inputs to the MBDyn structural model. With these structural component details, a complete MBDyn structural model was used. The spar and ribs were modeled as beam elements and the wing covering was modeled as 2-D shell element. The VICON system was used for shaker frequency experiments. 105 2.3.2 Wing Frequency Verification Tests The predicted natural frequency of each wing member was verified through modal tests using a “shaker”. A shaker is a device that excites an object or structure. The shaker device was a LDS V101 Permanent Magnet Shaker produced by Bruel and Kjaer Inc. This shaker model is capable of achieving an 8.9 N (2 lbf) peak sine thrust. The continuous peak-to-peak displacement is 2.5 mm and the frequency range is 5–12,000 Hz. The shaker was controlled using a Hewlet-Packard Agilent 8656B RF Signal Generator. The shaker was fixed to the leading edge at the root rib, mid rib, and tip rib of the flapping wing frame. A swept sine signal was used as the input for normal mode testing. The result was compared to the prediction values. Shaker tests also allowed for the determination of the mass for each structural member (see section 3.1.1). Figure 2.16 shows the attachment of the wing to the shaker for the modal experiments. 2.4 Chapter Summary This chapter contains a discussion of the numerical and experimental methods utilized in this study. The aeroelastic solver is a coupled OVERTURNS/MBDyn solver modified for low Reynolds number flows, and the details were discussed. Details of the flapping-wing test rig, the PIV system, force measurement technique, wind tunnel, and motion tracking setups are provided. The PIV was used to acquire velocity field measurements in several planes along the wing span. Several challenges associated with performing flow measurements on a flapping-wing MAV were also 106 (a) Root rib mounted. (b) Mid rib mounted. (c) Tip rib mounted. Figure 2.16: VICON motion tracking “shaker” test for modal analysis. 107 described. The VICON motion tracking system allowed for the measurement of the wing structural properties and wing kinematics. Finally, the methods for determining the structural properties of the flexible wing were presented. 108 Chapter 3: Comparison of Predicted and Experimental Re- sults This chapter provides a comparison of the simulation results with the experiment results. First, the structural solver was validated and then the flowfield. Flow features produced from the experiment and the simulation of the flapping wing for various flap frequencies and pitch angles are compared using PIV and CFD-CSD respectively for rigid and flexible cases. Force measurements are also used for simulation comparison. Of particular interest in the present study is a comparison of the formation and convection of the LEV and its effects on the unsteady aerodynamic forces produced on the wing. 3.1 Comparison of Structural Dynamic Solver to Experiments 3.1.1 Flexible Wing Static Loading Comparison Static structural tests were performed on each wing frame member and the whole wing (covering included). VICON markers were placed along the root, mid, and tip ribs (3 span locations). Loads were placed at tips of the members. Figure 3.1 shows representative images of the covered wing, wing loading, and VICON markers 109 (Left column: Figs. 3.1(a), 3.1(c), and 3.1(e)) along with deflection results for the chordwise ribs (Right column: Figs. 3.1(b), 3.1(d), and 3.1(f)). The left column contains the predicted (MBDyn) results for comparison to the experiment deflections. When comparing the results, it is seen that there is good correlation (predicted within < 5% difference compared to measured values). The analysis discussed in section 2.3 for large bending deflections was used to assess the EI for the various flexible wing members. 3.1.2 Flexible Wing Dynamic Tests 3.1.2.1 Comparison of Flexible Wing Frame Fundamental Frequencies For further comparison of the flexible wing simulation to the experiment, the natural (fundamental) frequency of the wing frame was obtained by exciting the wing frame from the root tab. The shaker was given a 15 second sinusoidal ramp to 50 Hz while VICON data was taken for 20 seconds. The displacements of the markers located at the trailing edges of the root, mid, and tip members were measured using the motion tracking system. Figure 3.2 shows the raw signal of the tip displacements and the corresponding Fast Fourier Transform (FFT) to assess the signal for the shaker tests. The FFT is used to convert the time/space data to the frequency domain. The left column (Figs. 3.2(a), 3.2(c), and 3.2(e)) shows the raw displacements from the trailing edge root, mid, and tip ribs. Several clear peaks in the displacement are observed from the raw signal. These are the instances in time where the resonant frequencies of the 110 Wing Root rib loading (a) Root spar load. 2 4 6 8 10 -25 -20 -15 -10 -5 0 Root Rib Experiment (VICON) Predicted (MBDyn) Mid Rib Tip Rib Loading, g Dis pla ce me nt, m m (b) Load at trailing edge of root rib. Wing Mid rib loading (c) Mid spar load. 2 4 6 8 10 -30 -25 -20 -15 -10 -5 0 Root Rib Mid Rib Tip Rib Loading, g Di sp lac em en t, m m (d) Load at trailing edge of mid rib. Wing Tip rib loading (e) Tip spar load. 2 4 6 8 10 -50 -40 -30 -20 -10 0 Root Rib Mid Rib Tip Rib Loading, g Dis pla ce me nt, m m (f) Load at trailing edge of tip rib. Figure 3.1: Left side, images of static wing loading tests with VICON markers and Right side, load deflection curves for chordwise ribs. 111 frame members occur. Figure 3.2(b) shows that the frequency at which the peak signal occurs at about 30 Hz for the tip displacements of the root spar. In Fig. 3.2(d) and 3.2(f) show FFT signal peaks at about 16.5, 23.5 and 31 Hz. There are several peaks observed in the tip displacement signal because the frame is a coupled system. The shaker test for the frame was modeled using the MBDyn solver. The frequency content of the frame is compared to the predicted results and displayed in Table 3.1. The predicted natural frequencies were calculated by computing the eigenvectors and eigenvalues within the structural solver. In comparing the natural frequency content, the percent difference in the predicted natural frequencies and the measured remained between 1.05% and 7.93%. The resonant frequencies gives a great deal of information about the charac- teristics of a system (e.g., mode of vibration is characterized by a modal frequency and a mode shape). Figure 3.3 shows the mode shapes for the wing frame structural members. For the first mode of the root spar (Fig. 3.3(a)), it is seen that the deflection of the root rib does not have a significant effect on the tip and the middle ribs. However, the middle rib has an effect on the tip spar (see Fig. 3.3(b)). This oscillatory coupling of the individual degrees of freedom is also seen as the tip rib affects the leading edge spar and the middle rib (see Fig. 3.3(c)). Because each wing member is interconnected by the wing covering, the natural frequency is not approached when the wing is flapping for the current tests (¡10 Hz). 112 0 5 10 15 20−6 −4 −2 0 2 4 6 Time,s Displ acem ent, m m (a) Root spar displacement. 0 10 20 30 40 500 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Frequency (Hz) |Y(f)| (b) Root spar FFT. 0 5 10 15 20−6 −4 −2 0 2 4 6 Time,s Displ acem ent, m m (c) Mid spar displacement. 0 10 20 30 40 500 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Frequency (Hz) |Y(f)| (d) Mid spar FFT. 0 5 10 15 20−6 −4 −2 0 2 4 6 Time,s Displ acem ent, m m (e) Tip spar displacement. 0 10 20 30 40 500 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Frequency (Hz) |Y(f)| (f) Tip spar FFT. Figure 3.2: Raw signal of tip (TE) displacements and Fast Fourier Transform from motion tracking experiment. 113 (a) Root spar mode shape. (b) Mid spar mode shape. (c) Tip spar mode shape. Figure 3.3: Mode shapes for wing frame structural members. 114 Table 3.1: Comparison of flexible chordwise member natural frequencies. Structural Member Experiment ωn Predicted ωn Root Rib 31.18 Hz 32.01 Hz Middle Rib 23.97 Hz 22.07 Hz Tip Rib 16.89 Hz 17.07 Hz 3.1.2.2 Comparison of Flexible Wing Frame Flapping Motion The VICON system was used to compare experiment and predicted results for the complex kinematics of the flexible wing frame throughout the entire stroke. Again, the tip displacements were tracked using the VICON system and simulations were conducted using MBDyn. Figure 3.4 shows a representative result from the dynamic tests. It is seen that the upstroke and downstroke are the same magnitude, resulting in a “symmetric” kinematic profile. It is seen that the experiment and simulation are periodic; however, there is some higher frequency content in the experiment that is not accounted for in the simulation. For the VICON data, the tip motion is not purely sinusoidal because of the inertial effects of the wing components. This difference in the kinematics is caused by inertial effects. To account for aerodynamic damping in the simulation, Rayleigh damping was used within the simulation to account for damping within the system. Rayleigh damping is an empirical means to account for damping, and was applied to the wing frame simulations. This damping is associated with linear theory, and is used for isolating each mode [183]. This proportional damping is a widely used approach to model dissipative forces in complex structures and it has been used in various dynamic problems. The proportional damping model expresses the 115 damping matrix as a linear combination of the mass and stiffness matrices. C = α1M + β1K (3.1) where α1 and β1 are real scalars. With Raleigh’s proportional damping, the modal damping factors have a special form: α1 + β1ω 2 i = 2ωiξi (3.2) where resonant frequencies, ωi, are obtained through modal analysis, and ξi are damping ratios. Rearranging gives this equation results in the following expression: ξi = α1 2ωi + β1ωi 2 (3.3) Assigning values to α1 and β1, allows for the filtering or retention of the higher mode effects. For this study, β1 was set to 1.0e−4 and the higher harmonic term α1 is set to zero because no tests were performed to achieve the coefficients. Figure 3.4 shows a comparison of the motion tracking experiment and predicted (MBDyn) wing frame tip displacements. The flap frequency for this test was 6 Hz. Even with this simplification, the flap characteristics are captured in that the sinusoidal nature of the wing flapping motion, the flapping phase, and the flap amplitude are comparable. This worked well for the comparison, and the inertial effects seen when the frame was flapping were further damped when the wing covering was present. 116 0 0.2 0.4 0.6 0.8 1−200 −150 −100 −50 0 50 100 150 200 Time (sec) Dis pla cem ent (mm ) TipMidRoot (a) Motion capture experiment tip displacements. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -200 -150 -100 -50 0 50 100 150 200 Time (sec) D is pl ac em en t ( m m ) (b) Predicted tip motion. Figure 3.4: Comparison experiment and predicted (MBDyn) wing frame tip displace- ments for 6 Hz flapping. 117 3.1.3 Comparison of Aerodynamic Solver (CFD-CSD) and PIV Re- sults 3.1.3.1 Rigid and Flexible Wing Flowfield Comparison To visualize the flow, the vorticity field, ~ω = ~∇× ~V , was acquired by computing the velocity gradients of the PIV velocity field using a second order accurate Least Squares differencing method. The Least Squares method was chosen to reduce measurement uncertainty from the PIV [179]. For the contours presented, the wing is at the middle of the downstroke, t = T/4 (mid-downstroke), and the free stream velocity is 3 m/s from left to right. The masked wing in the PIV is shown as black squares. Every 5th and 20th vector has been plotted for the PIV and CFD, respectively. The pitch angles and reduced frequencies presented were chosen to highlight the similarities and differences for comparison of the simulation and experiment at low and high pitch angles and flap frequencies. To define the wing kinematics, Fig. 3.5 shows a perspective view of the wing rotation and the wing’s nomenclature. The solid wing represents an initial time within the wing half stroke and the translucent wing represents the wing at a later instant in time within that same half stroke. The stroke angle is given as the angle Ψ and is defined as the angle through which the wing rotates. The pitch angle, θ, is defined as the angle of incidence of the wing, b is the distance from the rotational axis to the wing tip, y/b is the non-dimension distance along the span, and c is the chord. 118 20% Tip Midspan Span locations Root Rotational axis θ y x z Ψ c b Figure 3.5: Perspective view schematic of the rotating wing setup. Figures 3.6(a), 3.6(b), and 3.7(a) and 3.7(b) show representative rigid wing PIV and CFD results (vorticity contours) for span locations of y/b = [0.3 – 0.9], at t = T/4, mid-downstroke. Here, the salient aerodynamic features in the flow on the flapping wing can be seen at the different chordwise locations, including the leading-edge vortex (LEV), the trailing edge vortex, and small eddies from a turbulent shear layer. Here the blue regions denote locations of increased (positive) vorticity. For the θ = 24◦ (at 4 Hz) and θ = 12◦ (at 6 Hz) cases, a leading edge vortex has formed at all span locations. The reduced frequencies, k, are 0.32 and 0.48 respectively. In Figs. 3.6(a) and 3.7(a), this LEV can be seen as tightly structured near the wing root but greatly expands until approximately y/b = 0.6 which is when the LEV starts to diffuse for the PIV. For the CFD, the LEV grows in size along the chord also out to the y/b = 0.6 span location; however, the size of the LEV is decreased past that point out toward the wing tip (see Figs. 3.6(b) and 3.7(b)). 119 (a) PIV, θ = 24◦, Ω = 4 Hz (b) CFD, θ = 24◦, Ω = 4 Hz Figure 3.6: Comparison of representative chordwise PIV and CFD vorticity contours at t = T/4, (40◦ flapping amplitude, θ = 24◦, 4 Hz flapping frequency, V∞ = 3 m/s, k = 0.32). 120 Note that an increase in flap frequency results in an increase in vorticity magnitude. There are some differences seen in the 3-D vorticity flowfields, e.g., the vorticity near the tip appears to be more diffused for PIV (an effect of averaging). Also, the CFD dpes not capture the same shed vorticity near the wing tip. This is because the prediction is not capturing the exact same behavior caused by the interraction of the LEV and the tip vortex. Therefore, velocity profiles will be presented to quantify the differences. Overall, there is a qualitative correlation between CFD and PIV for all the span locations and many of the flow features seen in the PIV measured flowfield are captured in the CFD-predicted results. Note that an increase in flap frequency, results in an increase in vorticity magnitude. For θ0 = 12◦ (4 Hz), 24◦ (6 Hz), and 0◦ (8 Hz) cases, again a LEV is formed at all span locations (Figs. 3.8 to 3.10). In Figs. 3.8(a) and 3.9(a), this LEV can be seen as tightly structured near the wing root, but greatly expands until approximately y/b = 0.6 and then the LEV starts to diffuse and move further away from the wing’s surface for the PIV. For the coupled CFD-CSD predictions, the LEV grows in size along the chord also out to the y/b = 0.6 span location; however, the size of the LEV is decreased past that point out toward the wing tip (see 3.8(b) and 3.9(b)). There is a good correlation between coupled CFD-CSD and PIV for all the span locations and many of the flow features seen in the PIV measured flowfield are captured in the predicted results. Figures 3.9(a) and 3.9(b) show an increase in flap frequency results in an increase in vorticity magnitude. Though there is good agreement, the coupled CFD-CSD is not capturing the same levels of shedding compared to the experiment at y/b > 0.5. Figures 3.10(a) and 3.10(b) show the wing undergoing large amounts of 121 (a) PIV, θ = 12◦, Ω = 6 Hz (b) CFD, θ = 12◦, Ω = 6 Hz Figure 3.7: Comparison of representative chordwise PIV and CFD vorticity contours at t = T/4, (40◦ flapping amplitude, θ = 12◦, 6 Hz flapping frequency, V∞ = 3 m/s, k = 0.48). 122 (a) PIV, θ0 = 12◦, Ω = 4 Hz (b) CFD, θ0 = 12◦, Ω = 4 Hz Figure 3.8: Comparison of representative chordwise PIV and CFD vorticity contours at t = T/4, (40◦ flapping amplitude, θ0 = 12◦, 4 Hz flapping frequency, V∞ = 3 m/s, k = 0.32). 123 (a) PIV, θ0 = 24◦, Ω = 6 Hz (b) CFD, θ0 = 24◦, Ω = 6 Hz Figure 3.9: Comparison of representative chordwise PIV and CFD vorticity contours at t = T/4, (40◦ flapping amplitude, θ0 = 24◦, 6 Hz flapping frequency, V∞ = 3 m/s, k = 0.64). 124 (a) PIV, θ0 = 0◦, Ω = 8 Hz (b) CFD, θ0 = 0◦, Ω = 8 Hz Figure 3.10: Comparison of representative chordwise PIV and CFD vorticity contours at t = T/4, (40◦ flapping amplitude, θ0 = 0◦, 8 Hz flapping frequency, V∞ = 3 m/s, k = 0.64). 125 flexibility during the downstroke. The reduction in angle of attack, caused by wing flexibility, results in a smaller LEV formed on the wing’s upper surface. The next section contains details about the vorticity contours and normalized velocity profiles for both the CFD and PIV (rigid and flexible wing) at the 50% span location. For the rigid and flexible wing, downstream velocity cuts (i.e., vertical sectional cuts) were taken along the Z-direction at the X/c = 1.1, 1.5, and 2 (X/c = 0 is the leading edge) to give a quantitative comparison of the CFD and PIV wakes. The downstream velocity cuts have been denoted by the letters A, B, C and A’, B’, C’ for the PIV and CFD respectively. The vertical velocity cuts in the wake are 2 chords in length. The velocity cuts show good agreement and the shear layer is captured by the simulation; however, there are slight differences when comparing the shear layer angle and size of the shed vortices. Furthermore, a swirl velocity, Vθ, profile was taken through the LEV for a quantitative comparison at the LEV structure. The center of the vortex was located using the λ1-function. The LEV structure required a more precise definition of the vortex core, so the non-local vortex detection scheme of Graftieaux et al. [184] was employed. This function was used to locate and quantify the rotation rate of local swirl regions. The λ1-function is given by λ1(P ) = 1 N ∑ S sin(θM) where N is the number of points inside a two-dimensional area S centered around point P . The angle θM is formed by the position vector from P to M and the velocity 126 vector at point M . The value of λ1 can range from -1 to 1. |λ1| = 1 when the point P is at the center of an ideal vortex. Note that the value of λ1 depends on the direction of the velocity field, not its magnitude. The center of a vortex is taken to be the local maximum of |λ1|, and the vortex core is the region within the |λ1| ≥ 0.6 contour. Swirl velocity profiles were generated along a diagonal cut through the vortex core. The diagonal dashed line (parallel to the wing) shows the velocity cut through the LEV. The vortex center is denoted by a white dot. A specified number, n, of equally spaced points along the cut (n = 200) were assigned to the line. Along this cut, the u and w velocities corresponding to points along x and z in the velocity field were obtained. Because PIV yields a planar grid of velocity vectors, not all points on the diagonal fall on a node in the velocity field grid. If a point n on the line did not fall on a grid node, the u and w velocity components were interpolated based on the neighboring vectors. To calculate the swirl velocity, the u and w velocity values (from the Cartesian grid) on the diagonal were converted to polar coordinates to give the radial (Vr) and tangential (Vθ) velocities of the vortex. Each swirl velocity profile is a plot of Vθ (normalized by V∞) as a function of position along the diagonal cut. The upper left of the diagonal dashed line corresponds to the bottom of the velocity profile.The vortex is analyzed by determining the core radius, rc, defined as the radial distance from the location where Vθ is a maximum to the center of the vortex, determined from the swirl velocity profiles. Figures 3.11(a) and 3.11(b) show the vorticity fields with overlaid velocity for the CFD and PIV results. Here, the wing is in its mid-downstroke position and the focus is on the mid-span plane. The wing is at 0◦ and the flap frequency is 127 4 Hz. The vorticity gives the structure of the vortex. It is seen that the velocity magnitude is very similar; however, vortex structure is different. Both the averaged PIV and instantaneous CFD vorticity fields show the LEV as a tight concentration of vorticity near the leading edge encompassing about 25% of the chord. There is a small amount of opposite sign vorticity, which shows as red within the LEV for the predicted result. The shear layer appears as more diffused for the PIV. This is a result of averaging and the turbulence model used for the simulation, which was developed for higher Reynolds number flows. Even with these considerations, the velocity fields are quite similar for the PIV and CFD in that the LEV and smaller trailing edge vortices form for both. The LEV is bound to the wing’s upper surface. Figures 3.11(c) and (d) shows the downstream horizontal velocity component profile at three different downstream locations. A reduction in velocity behind the wing is representative of a wing that is producing drag via a decrement in momentum. The normalized u velocity just behind the wing at 1.1 chord length is about half the freestream velocity. In Fig. 3.11(e), which compares the instantaneous normalized vortex swirl velocities (Vθ/V∞) as a function of nondimensional radial distance from the vortex center (rc/c), it is seen that the vortex from CFD and PIV are of similar size and magnitude. Here the vortex size is about 0.23 chord lengths as determined from the peak to peak distance in the swirl velocity profiles. In general, the velocity fields show very good agreement as seen by the velocity profiles. In Fig. 3.12, the leading edge vortex centers (determined by λ1-function) are the centers of the velocity cuts through the vortex. Here, the wing pitch angle is 12◦ and the locations of the leading edge vortex centers for the measurement and CFD 128 Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 D A B C (a) PIV vorticity contour with overlaid velocity field. Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 A’ B’ C’ D’ (b) CFD vorticity contour with overlaid velocity field. −3 −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona l Z−c oordi nate, in ch ords Normalized streamwise velocity, U/Voo A − 1.1 chordB − 1.5 chordC − 2.0 chord (c) PIV downstream U/V∞ velocity distribution. −3 −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona lized Z−co ordin ate, i n cho rds Nondimensionalized streamwise velocity, U/Voo A’ − 1.1 chordB’ − 1.5 chordC’ − 2.0 chord (d) CFD downstream U/V∞ velocity distribu- tion. −3 −2 −1 0 1 2−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 Norm alize d Vo rtex R adius , r/ch ord Normalized Vortex Swirl Velocity, (Vθ/U) D − PIVD’ − CFD (e) PIV and CFD normalized LEV swirl velocity profiles. Figure 3.11: Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, rigid wing ,θ = 0◦ at 4 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). 129 are close in proximity (about 5% difference in straight line distance from the leading edge). Here the LEV is bound to the wing’s surface. Though the recirculation region near the leading edge is not well predicted in the CFD, and the trailing edge vortices have a different size and trajectory. Also, the trend line for the velocity profiles and magnitude of the flow are similar. There are some differences in comparison that are highlighted by the vorticity contours (e.g. vortex structure); however, the peak swirl velocity and the core radius of the LEV for the PIV and CFD are comparable, as shown in the 6 Hz, 12◦ case. It is seen that the trajectories of the vortices shed from the trailing edge have a different trajectory. It can be seen in Figs. 3.12(c) and 3.12(d) that profiles C and C’ do not match. This is because the C profiles cut through a significant region of vorticity present in the downstream wing wake. The trailing edge vortex present in the wing wake of the PIV image was not appropriately captured in the CFD simulation. The magnitudes of the velocity deficit are similar up to about 1.5 chord lengths behind the leading edge. Nevertheless, the leading edge vortex and shed trailing edge vortices are predicted with similar values for vorticity, and velocity magnitude. Figure 3.12(e) shows that the LEV encompasses almost half of the wing’s upper surface; however, the swirl velocity profiles do not quite have the same magnitudes. The higher flap frequency does produce stronger vortices with higher swirl velocities. Figure 3.13 shows the lowest tested flap frequency, 4 Hz, but the highest geometric pitch angle, θ = 24◦. When the wing is in the mid-downstroke position, the vortex has diffused, and flow separation is observed on much of the upper wing’s surface. The remnants of the LEV are present; however, the circulation is not bound 130 Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 D A B C (a) PIV vorticity contour with overlaid velocity field. Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 D’ A’ B’ C’ (b) CFD vorticity contour with overlaid velocity field. −3 −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona l Z−c oordi nate, in ch ords Normalized streamwise velocity, U/Voo A − 1.1 chordB − 1.5 chordC − 2.0 chord (c) PIV downstream U/V∞ velocity distribution. −3 −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona lized Z−co ordin ate, i n cho rds Nondimensionalized streamwise velocity, U/Voo A’ − 1.1 chordB’ − 1.5 chordC’ − 2 chord (d) CFD downstream U/V∞ velocity distribu- tion. −3 −2 −1 0 1 2−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 Norm alize d Vo rtex R adius , r/ch ord Normalized Vortex Swirl Velocity, (Vθ/Voo) D − PIVD’ − CFD (e) PIV and CFD normalized LEV swirl velocity profiles. Figure 3.12: Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, rigid wing, θ = 12◦ at 6 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). 131 to the wing’s surface. The trends in the downstream wake are very similar. However, the velocity deficit is under predicted, and the swirl velocity profiles are dissimilar (see Figs 3.13(b) and 3.13(c)). The variation in the predicted and the measured swirl velocity profiles is a result of the unsteady nature of the vortex. Though the swirling flows near the leading edge are similar, the vortex identification method (determined by λ1-function) has assessed the vortex center of the predicted result to be about 0.25 chord lengths further away from that of the PIV result. There is also far less opposite sign vorticity within the experiment results. The predicted and measured flow fields are similar, however; the nature of the flow at this point is such that the flow features are more aperiodic, which accounts for the difference in the profiles. In general, the flow solver predicted the overall flow features. At higher angles of attack, the flow on the wing’s upper surface becomes more aperiodic. This is the reason some of the flow details are not shown in the averaged flowfields. However, the general trends are reflected by the CFD analysis. These results serve as a comparison between CFD predictions, and experimental measurements show very good agreement when the LEV is attached. Next, the flexible wing results will be discussed. Figures 3.15(a) and 3.15(b) show the vorticity contours for the flexible wing at 6 Hz, for θ0 = 24◦, in the mid-downstroke position. The mid-span location is shown. Figure 3.15(e) compares the instantaneous normalized vortex swirl velocities (Vθ/V∞) as a function of nondimensional radial distance from the vortex center (rc/c); it is seen that the vortex from coupled CFD-CSD and PIV are of similar size and magnitude. 132 Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 D A B C (a) PIV vorticity contour with overlaid velocity field. Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 D’ A’ B’ C’ (b) CFD vorticity contour with overlaid velocity field. −3 −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona l Z−c oordi nate, in ch ords Normalized streamwise velocity, U/Voo A − 1.1 chordB − 1.5 chordC − 2.0 chord (c) PIV downstream U/V∞ velocity distribution. −3 −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona lized Z−co ordin ate, i n cho rds Nondimensionalized streamwise velocity, U/Voo A’ − 1.1 chordB’ − 1.5 chordC’ − 2.0 chord (d) CFD downstream U/V∞ velocity distribu- tion. −3 −2 −1 0 1 2−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 Norm alize d Vo rtex R adius , r/ch ord Normalized Vortex Swirl Velocity, (Vθ/U) D − PIVD’ − CFD (e) PIV and CFD normalized LEV swirl velocity profiles. Figure 3.13: Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, rigid wing, θ = 24◦ at 4 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). 133 Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −1 −0.5 0 0.5 1 A B C D (a) PIV vorticity contour with overlaid velocity field. Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −1 −0.5 0 0.5 1 D’ C’B’A’ (b) CFD vorticity contour with overlaid velocity field. −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona l Z−c oordi nate, in ch ords Normalized streamwise velocity, U/Voo A − 1.1 chordB − 2.0 chordC − 3.0 chord (c) PIV downstream U/V∞ velocity distribution. −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona lized Z−co ordin ate, i n cho rds Nondimensionalized streamwise velocity, U/Voo A’ − 1.1 chordB’ − 2.0 chordC’ − 3.0 chord (d) CFD downstream U/V∞ velocity distribu- tion. −4 −2 0 2−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 Norm alize d Vo rtex R adius , r c/ch ord Normalized Vortex Swirl Velocity, (Vθ/Voo) D − PIVD’ − CFD−CSD (e) PIV and CFD normalized LEV swirl velocity profiles. Figure 3.14: Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, flexible wing, θ0 = 12◦ at 4 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). 134 Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −1 −0.5 0 0.5 1 A B C D (a) PIV vorticity contour with overlaid velocity field. Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −1 −0.5 0 0.5 1 D’ A’ C’B’ (b) CFD vorticity contour with overlaid velocity field. −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona l Z−c oordi nate, in ch ords Normalized streamwise velocity, U/Voo A − 1.1 chordB − 2.0 chordC − 3.0 chord (c) PIV downstream U/V∞ velocity distribution. −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona lized Z−co ordin ate, i n cho rds Nondimensionalized streamwise velocity, U/Voo A’ − 1.1 chordB’ − 2.0 chordC’ − 3.0 chord (d) CFD downstream U/V∞ velocity distribu- tion. −4 −2 0 2−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 Norm alize d Vo rtex R adius , r c/ch ord Normalized Vortex Swirl Velocity, (Vθ/Voo) D − PIVD’ − CFD−CSD (e) PIV and CFD normalized LEV swirl velocity profiles. Figure 3.15: Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, θ0 = 24◦ at 6 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). 135 Figure 3.16 shows a comparison of the PIV experiment and CFD-CSD solver flexible wing results for the 8 Hz flapping case when the wing’s initial pitch angle was set to 0◦. Here, the wing is at the mid-downstroke position and the y/b = 0.5 span location is shown. The CFD-CSD solver is slightly under-predicting the wing deflection in that the wing angles are about −18◦ and −15◦ for the experiment and prediction, respectively. This slight difference in wing angle is likely caused by a nonuniformity in the wing material near the trailing edge of the wing as exhibited by the sharp change in wing angle at about the mid-chord location (see Fig. 3.16(a)). Nevertheless, the concentration of vorticity at the wing’s leading edge are comparable in size for both the experiment and the CFD-CSD prediction. In the downstream wake, the shear layer at the trailing edge for the PIV and the CFD-CSD have a similar trend in length and diffusion characteristics, but they are at slightly different angles. This is reflected in the downstream velocity profiles shown Figs 3.16(c) and 3.16(d). For the CFD-CSD, the profile A’ cuts through the beginning of the trailing edge shear layer, while the profile A intersects the shear layer further downstream. Figure 3.16(e) shows the swirl velocity profiles for the PIV and CFD-CSD velocity fields. The location of the vortex center is predicted and the vortex structures are similar, however, there are some discrepancies in the prediction of the flow ahead of the wing’s leading edge. In general, the trends is predicted. Because the LEV is a continuous vortex filament, chordwise sections of the LEV do give insight into its behavior. In general, during wing flapping, the vortex goes through a periodic growth, convection, and dissipation process. If this vortex could remain coherent and attached (bound), perhaps lift could be improved during the 136 Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −1 −0.5 0 0.5 1 A B C D (a) PIV vorticity contour with overlaid velocity field. Normalized X−coordinate, x/c Norm alize d Z−c oordi nate, z/c −0.5 0 0.5 1 1.5 2 2.5 3 3.5 −1 −0.5 0 0.5 1 C’B’A’ D’ (b) CFD vorticity contour with overlaid velocity field. −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona l Z−c oordi nate, in ch ords Normalized streamwise velocity, U/Voo A − 1.1 chordB − 2.0 chordC − 3.0 chord (c) PIV downstream U/V∞ velocity distribution. −2 −1 0 1 2−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nond imen siona lized Z−co ordin ate, i n cho rds Nondimensionalized streamwise velocity, U/Voo A’ − 1.1 chordB’ − 2.0 chordC’ − 3.0 chord (d) CFD downstream U/V∞ velocity distribu- tion. −4 −2 0 2−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 Norm alize d Vo rtex R adius , r c/ch ord Normalized Vortex Swirl Velocity, (Vθ/Voo) D − PIVD’ − CFD−CSD (e) PIV and CFD normalized LEV swirl velocity profiles. Figure 3.16: Velocity vectors and vorticity field with overlaid velocity profile cuts, for corresponding normalized downstream and swirl velocity profiles, flexible wing, θ = 0◦ at 8 Hz. (Note: y/b = 0.5, t = T/4, Reref = 15, 000). 137 flap cycle. One reason for the persistence of an LEV is the stretching of the vortex as it convects over the wing’s surface; however, the change in effective angle of attack may also be a cause. The vortex strength is related to its persistence. In general, the effects of the vortical structures are observed in the experiment measurement (velocity profiles) and the coupled CFD-CSD is predicting the overall flow features and the general trends. There are some limitations, such as slight differences in the trailing edge vortices, but the limitations associated with the predictive model are not so significant that insight could not be gained from this analysis. These comparisons between coupled CFD-CSD predictions and experimental measurements provide confidence in the CFD-CSD results when the LEV is coherent. However, to further improve the understanding of the flow physics and the capability of the modeling tool, the LEV strength based on the CFD-CSD analysis and PIV flowfields are compared. 3.2 Comparison of Sectional LEV Strength LEV circulation, which is related to the lift, is quantified for the predicted and experimental results. This was done for both the rigid and flexible wing. 3.2.1 Rigid Wing LEV Circulation (PIV and CFD) The velocity field of both the CFD and PIV were used to estimate the strength of the LEV for comparison. The process involves choosing a suitable integration contour (i.e., around the LEV in this case) and numerically evaluating the closed- 138 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Nondimensional span location, y/b LEV Circ ulati on, (Γ LEV ), m 2 s− 1 Avg. Instantaneous PIV − θ = 0oInstantaneous CFD − θ = 0o (a) θ0 = 0◦ 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Nondimensional span location, y/b LEV Circ ulati on, (Γ LEV ), m 2 s− 1 Avg. Instantaneous PIV − θ = 12oInstantaneous CFD − θ = 12o (b) θ0 = 12◦ 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Nondimensional span location, y/b LEV Circ ulati on, (Γ LEV ), m 2 s− 1 Avg. Instantaneous PIV − θ = 24oInstantaneous CFD − θ = 24o (c) θ0 = 24◦ Figure 3.17: Distribution of averaged LEV circulation, Γv along the span, rigid wing, Ω =4 Hz, (40◦ flap amplitude, V∞ = 3 m/s) with error bars showing standard deviation. 139 loop velocity integral. This loop must completely enclose the LEV, but without intersecting any other extraneous circulation. The equation for circulation is given below as Γv = − ∫ ~V · ~ds = − ∫ ∫ (∇× ~V ) · ~dS (3.4) where ~V is the velocity vector, ~ds is the directed line segment at a point on a predefined contour, ∇× ~V is the vorticity, and ~dS is the area enclosed by a curve [17]. For this analysis, a vorticity threshold was set (−400 < ω < 400) so that low magnitude extraneous vorticity would not be included in the integration contour that was used to enclose the LEV. Next the vorticity field is scanned in order to locate pockets of vorticity of a certain size (3 × 3 grid square) within that threshold. These regions are then stored as grid locations. Line contours are drawn around these concentrations of vorticity so that an area integral can be assessed to obtain the circulation for the vortices in the PIV experiment (with error bars) and those generated by CFD. Figures 3.17(a), 3.17(b), and 3.17(c) show the phase-averaged spanwise distribu- tion of LEV circulation over the wing for pitch angles of 0◦, 12◦, and 24◦ respectively for the PIV and CFD. The flap frequency here is 4 Hz. LEV size and strength does vary at different span locations along the wing. This is because an increase in vorticity would also mean that there is an increase in net circulation around the wing, thereby causing increased lift when the LEV is present. The measured peak in the circulation was observed to occur between the y/b = 0.6 and 0.7 span locations with a rapid 140 decrease in circulation toward the wing tip and root. The low values of circulation near the wing root and tip arose because of the high flap-induced velocities there and also because the flow at these locations is partially separated. Even though the CFD simulation is able to predict the circulation of the inboard sections with sufficient accuracy there is some over-prediction in the outboard sections. The difference in the predicted and measured values for the outboard section at θ0 = 12◦ is largely due to the simulation’s inability to accurately resolve the flow when the tip vortex interacts with the LEV. In general, circulation is under-predicted further outboard for this case when there is significant mixing. There is an over-prediction along the span when the wing is at 24◦. This is because the flow is fully separated and difficult to resolve, although there is less of an effect from the tip vortex. Figure 3.18 shows the distribution of LEV strength for a variation in flap frequency for the y/b = 0.5 span location. The pitch angles are θ0 = 0◦, 12◦, and 24◦. Here the general trend is an increase in circulation as the wing flap frequency is increased. There is very good agreement in LEV strength value at the 0◦ pitch angles and lower frequencies. However, as the flap frequency is increased along with pitch angle there is less agreement between the measurements and the simulation. The effect of the tip vortex is more prevalent at higher frequencies making it difficult to resolve because of mixing with the LEV, especially at high reduced frequencies. This is an exacerbation of the effect of the interaction of the tip vortex with the LEV caused by the higher frequency flapping. 141 0 2 4 6 8 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Flap frequency, Hz LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV ExpCFD (a) θ0 = 0◦ 0 2 4 6 8 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Flap frequency, Hz LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV ExpCFD (b) θ0 = 12◦ 0 2 4 6 8 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Flap frequency, Hz LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV ExpCFD (c) θ0 = 24◦ Figure 3.18: Distribution of averaged LEV circulation, Γv for a variation in flap frequency, rigid wing, y/b = 0.5 (40◦ flap amplitude, V∞ = 3 m/s) at 50% span. 142 3.2.2 Flexible Wing LEV Circulation (PIV and CFD-CSD ) The velocity field of both the coupled CFD-CSD solver and PIV were used to estimate the strength of the LEV for comparison. The circulation was also calculated in the manner explained previously in that line contours were drawn around the concentrations of vorticity so that an area integral can be assessed to obtain the circulation for the vortices in the PIV experiment and those generated by the coupled CFD-CSD solver. Figures 3.19(a), 3.19(b), and 3.19(c) show the spanwise distribution of LEV circulation over the wing for pitch angles of 0◦, 12◦, and 24◦ respectively for the flexible wing case. Again, the flap frequency here is 4 Hz. The measured peak in the circulation was observed to occur between the y/b = 0.5 and 0.7 span locations with a rapid decrease in circulation toward the wing tip and root. The low values of circulation near the wing root and tip arise because of the high induced velocities there and also because the flow at these locations is partially separated. Even though the coupled CFD-CSD simulation is able to predict the circulation of the inboard sections with sufficient accuracy, there is some over-prediction in the outboard sections. Figure 3.20 shows the distribution of LEV strength for a variation in flap frequency for the y/b = 0.5 span location. The pitch angles are θ0 = 0◦, 12◦, and 24◦. The general trend is an increase in circulation as the wing flap frequency is increased. There is a very good agreement of LEV strength value at the 0◦ pitch angles and lower frequencies. 143 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Nondimensional span location, y/b LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV − θ = 0oCFD − θ = 0o (a) θ0 = 0◦ 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Nondimensional span location, y/b LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV − θ = 12oCFD − θ = 12o (b) θ0 = 12◦ 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Nondimensional span location, y/b LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV − θ = 24oCFD − θ = 24o (c) θ0 = 24◦ Figure 3.19: Distribution of averaged LEV circulation, Γv along the span, flexible wing, Ω =4 Hz, (40◦ flap amplitude, V∞ = 3 m/s). 144 0 2 4 6 8 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Flap frequency, Hz LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV − θ = 0oCFD − θ = 0o (a) θ0 = 0◦ 0 2 4 6 8 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Flap frequency, Hz LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV − θ = 12oCFD − θ = 12o (b) θ0 = 12◦ 0 2 4 6 8 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Flap frequency, Hz LEV Circ ulati on, (Γ LEV ), m 2 s− 1 PIV − θ = 24oCFD − θ = 24o (c) θ0 = 24◦ Figure 3.20: Distribution of averaged LEV circulation, Γv for a variation in flap frequency, flexible wing, y/b = 0.5 (40◦ flap amplitude, V∞ = 3 m/s) at 50% span. 145 LEV size and strength does vary at different span locations along the wing. This is because an increase in vorticity indicates that there is an increase in net circulation around the wing, thereby causing increased lift when the LEV is present. When the LEV is shed, the vorticity becomes less concentrated and diffused rather rapidly when the LEV starts to detach. This effect is most pronounced further outboard on the wing. 3.3 Airloads Comparison 3.3.1 Rigid Wing Comparison of Aerodynamic Forces from Flowfield Along with the flowfield, it is also important to validate the instantaneous forces predicted by CFD with test data. However, discerning the instantaneous aerodynamic force component on a flapping wing is a challenging task. This is because the total force is dominated by the inertial force component especially in the vertical direction (along Z-axis). Therefore, an attempt was made to calculate the vertical force (Fz) from the measured flowfield. A conventional method for obtaining lift from velocity field is by numerically calculating the circulation using a line integral along a closed contour enclosing the airfoil [17]. The size of the integration contour around the airfoil is increased until the circulation converges to a steady value. Even though such a method could work for the steady attached flows, its applicability is limited for highly unsteady, vortex-dominated separated flows, similar to those observed for flapping wings in the present study. In most of the present test cases, the control volume-based method did not converge to a unique value and, hence, could 146 not be used with sufficient confidence. Therefore, a different method was developed, which utilizes a rectangular control volume around the airfoil (ABCD) as shown in Fig. 3.21. This method requires the use of a momentum-based approach to calculate the vertical force acting on the airfoil. As shown in Fig. 3.21, the moving fluid exerts pressure and shear forces on the wing, which creates a resultant aerodynamic force per unit span, given by R. In turn, by Newton’s third law, the body exerts equal and opposite pressure and shear forces (−R) on the flow. This effect would manifest itself in the flowfield around the airfoil. This method is implemented using the integral form of the momentum equation given below ∂ ∂t ∫∫∫ ρ~V dv + ∫∫ (ρ~V · ~dS)~V = − ∫∫ p ~dS − ~R (3.5) The force acting on the airfoil (~R) can be calculated as follows: ~R = − ∂ ∂t ∫∫∫ ρ~V dv − ∫∫ (ρ~V · ~dS)~V − ∫∫ p ~dS (3.6) The first term on the right hand side of equation 3.6 is the acceleration of the fluid in the control volume and the second term is the momentum flux across the boundary of the control volume. These two terms can be calculated from the PIV measured flowfield. However, the third term ( ∫∫ p ~dS) cannot be obtained from the velocity field. One approach would be to consider a control volume much bigger than the airfoil so that the pressure along the boundary is the ambient pressure (P∞) and the surface integral would average to zero. However, a practical window size used in 147 Figure 3.21: Schematic explaining CZ calculation. a PIV measurement would not be large enough to compute this. Note that the PIV window size is 2 chords × 3.5 chords. Therefore, a hybrid approach is used in the present study where the pressure along the box sides was obtained from CFD. Also note that since the pressure always acts normal to the sides, only the pressure along sides AB and CD contributes to the vertical force (Fz). An alternative experimental approach to obtain pressure could be to use an array of pressure transducers on the ceiling and floor of the wind tunnel test section. Another method could be to compute the pressure from the PIV-measured velocity data, which is a very active area of research. The vertical force coefficient (Cz) obtained from the rigid wing PIV measured flowfield using the present method is compared with Cz obtained from CFD in 148 0 2 4 6 8 10 0 1 2 3 4 Flapping frequency (Hz) Ve rtic al for ce co eff icie nt, C Z 0 deg (PIV) 12 deg (PIV) 24 deg (PIV) 0 deg (CFD) 12 deg (CFD) 24 deg (CFD) Wing pitch angle Figure 3.22: Variation of vertical force coefficient (Cz) from rigid wing averaged PIV and instantaneous CFD flowfield data (50% span location, 4 Hz flapping). Fig. 3.22. This is for the 4 Hz case and data was used for the 50% span location (mid-downstroke). As shown in this figure, there is generally good agreement between the PIV and CFD for a wide range of pitch angles and flapping frequencies. This further validates the CFD analysis and also the present method of calculating forces. 3.3.2 Comparison of Predicted and Averaged Direct Force Measure- ments for the Rigid and Flexible Wing Along with the flowfield comparison, it is important to compare the forces predicted by CFD-CSD with test data. These direct force measurement experiments were performed with the flapping mechanism instrumented with a force transducer at the root. The wing frequency was varied for three pitch angles (θ = 0◦, 12◦, and 24◦). The flap amplitude was 40◦ and the wind tunnel speed was 3 m/s. Figure 3.23 149 shows the results for the averaged forces. Figures 3.23(a) and 3.23(b) show the resulting variation of average FZ and FX for the rigid wing. Here, 0◦ pitch flapping generates a propulsive force. This is a very important observation and will be examined in subsequent chapters. The production of this positive aerodynamic force in the chordwise direction (propulsive force) is called the Knoller-Betz effect [114]. For the θ0 = 0◦ pitch angle, the FZ force is symmetric for the up-stroke and down-stroke, and thereby the average vertical force is nearly zero. As the pitch angle is increased, the propulsive force is decreased and the wing begins to produce more drag. As the pitch angle is increased, the FZ force increases and the FX force goes further negative. Although pure flapping of this rigid wing is not a practical solution for real MAV-scale forward flight, these results provide insight into the effects of flapping on the wing. The flow solver results correlate well with the averaged force measurements. Figures 3.23(c) and 3.23(d) show the resulting variation of average FZ and FX with flapping frequency for flexible wings. It can be observed that the wing with root pitch angle of 0◦ produced almost zero FZ and 9.5 grams of FX at 10 Hz. There is some propulsive force generated at the zero pitch angle (similar to the rigid wing). On increasing the root pitch angle, the vertical force increases, but the propulsive force reduces. For the case with θ0 = 24◦, FZ of 12 grams and FX of 5 grams was generated at 10 Hz. In general, the horizontal force is improved by adding wing flexibility. The details for why this is the case will be discussed next in Chapters 4 and 5. As the pitch angle is increased, the FZ force increases and the FX force decreases. There 150 4 5 6 7 8 9 10−2 0 2 4 6 8 10 12 14 Frequency, Hz Verti cal fo rce, F Z (g) Exp − θ = 0oCFD − θ = 0oExp − θ = 12oCFD − θ = 12oExp − θ = 24oCFD − θ = 24o (a) Rigid wing variation of FZ. 4 5 6 7 8 9 10−4 −2 0 2 4 6 8 10 Frequency, Hz Horiz ontal force , F X (g ) Exp − θ = 0oCFD − θ = 0o Exp − θ = 12o CFD − θ = 12o Exp − θ = 24oCFD − θ = 24o (b) Rigid wing variation of FX. 4 5 6 7 8 9 10−2 0 2 4 6 8 10 12 14 Frequency, Hz Verti cal fo rce, F Z (g) Exp − θ = 0oCFD − θ = 0oExp − θ = 12oCFD − θ = 12oExp − θ = 24oCFD − θ = 24o (c) Flexible wing variation of FZ. 4 5 6 7 8 9 10−4 −2 0 2 4 6 8 10 Frequency, Hz Horiz ontal force , F X (g ) Exp − θ = 0oCFD − θ = 0oExp − θ = 12oCFD − θ = 12oExp − θ = 24oCFD − θ = 24o (d) Flexible wing variation of FX. Figure 3.23: Rigid and flexible wing comparison of averaged experiment and predicted vertical (FZ) and horizontal forces (FX) for various flap frequencies. is good agreement between averaged force measurements and CFD-CSD predicted result for a wide range of pitch angles and flapping frequencies. This does somewhat validate the CFD-CSD analysis for large wing deflections with highly unsteady and separated flows. 3.4 Chapter Summary The multibody dynamics solver was used to model the structure of a MAV flexible flapping wing. Cantilevered wing members were tested in both bending and 151 torsion to obtain the structural properties. Validation of the shell model was executed using static load experiments and analysis on the full wing. The capability of MBDyn to model plate structures with large deformations was successfully demonstrated. In- house experiments were also conducted to measure the deflections of the wing using a VICON system during flapping. Shaker tests were also performed to acquire the wing member natural frequencies for comparison to predicted values, and satisfactory correlation was obtained. The coupled aerodynamic and structural model were also compared to experimental results with satisfactory correlation. Once the 3-D CFD-CSD model was systematically compared with experiments, the model was used to understand the physics of lift and thrust production on a flapping wing in forward flight. The model was used for this instead of experiments because PIV measurements were only conducted at two instants during the flap cycle (mid-downstroke and mid-upstroke). Also, the computational analysis can provide more information such as the temporal and spatial variation of wing forces, which are very important for understanding the physics. 152 Chapter 4: Rigid and Flexible Wing Comparison (Aerody- namic Force Production) In the previous chapter, the results from the 3-D CFD-CSD analysis were compared systematically with the results from the PIV and force transducer wind tunnel experiments. In this chapter, the results of the computational analysis are utilized to better understand the flow physics about a flapping wing in forward flight and, more importantly, the influence of flexibility on the production of lift (FZ) and propulsive thrust (FX). The numerical analysis was used because it was more appropriate than using the experimental results since the PIV measurements were taken at only two points within the flap cycle (mid-downstroke and mid-upstroke). Additionally, the computational analysis is capable of providing more detailed information with respect to the spatial and temporal variations in the aerodynamic forces over the course of a flap cycle. It should be stated that throughout this study, the structural characteristics of the wing were kept constant with the goal being to compare the results of a structurally characterized flexible wing to those of a rigid wing. 153 4.1 Understanding Lift and Thrust Production for a Rigid Flapping Wing 4.1.1 Rigid Wing Force History It is important to understand the variation in forces over the flap cycle. Fig- ure 4.1 shows the variation in vertical force coefficient (CZ) and horizontal force or propulsive thrust coefficient (CX) of the flapping wing over the course of one flap cycle for the different combinations of flap frequency and pitch angle for the rigid wing. Note that, for (CX), a positive value indicates drag force while a negative value indicates propulsive thrust force. The flap frequencies tested during the experiments and numerical simulations (0 Hz, 4 Hz, 6 Hz, 8 Hz, and 10 Hz) correspond to reduced frequencies (k) of 0, 0.32, 0.48, 0.64 and 0.798 respectively. Figures 4.1(a) and 4.1(b) show the variation of CZ and (CX), respectively, over one flap cycle for a pitch angle of 0◦ and flapping frequencies of 0, 4, 6, 8 and 10 Hz. From Fig. 4.1(a) it can be seen that for the static case (0 Hz), as expected CZ is almost zero, but as the flapping frequency is increased the instantaneous CZ is also increased because of the increase in angle of attack induced due to flapping and also due to the higher resultant velocity experienced by the wing. For the 0◦ pitch angle, the upward force during the downstroke is almost exactly the same as the downward force during the upstroke, producing almost zero net vertical force over a flap cycle. However, as shown in Fig. 4.1(b), the wing produces positive propulsive thrust (negative (CX)) during 154 both upstroke and downstroke and the magnitude of thrust increases with flapping frequency. As the pitch angle is increased (12◦ and 24◦), as shown in Figs. 4.1(c) and 4.1(e), there is an asymmetry in the variation of CZ over the flap cycle. A larger magnitude positive CZ is produced during the downstroke in comparison to the negative CZ in upstroke resulting in a net positive vertical force. Also, there is an increase in the mean of CZ over the flap cycle with increased pitch angle. However, as shown in Figs. 4.1(d) and 4.1(f), CX follows the opposite trend because as the pitch angle is increased, the mean (CX) becomes more positive indicating that the wing is producing more drag than thrust. Also, for all the frequencies examined, the average propulsive thrust was always negative (positive (CX)) for 12◦ and 24◦ pitch angles. This is because the wing produces positive (CX) during the downstroke, the magnitude of which is much higher than the negative (CX) during the upstroke. In essence, for the cases with a non-zero pitch angle, in a majority of cases the positive vertical force and negative thrust force is produced during the downstroke of the flap cycle. Conversely, the negative vertical force and positive thrust is produced during the upstroke of the flap cycle. Table 4.1 provides a list of the averaged (CZ) and (CX) force coefficients over the duration of one flap cycle for the flap frequencies and pitch angles tested. Table 4.1: Chart of averaged CZ and CX for flap frequency and pitch angle (rigid wing) 4 Hz 6 Hz 8 Hz 10 Hz θ CZ CX CZ CX CZ CX CZ CX 0◦ -0.02 0.00 0.01 -0.01 -0.01 -0.04 -0.02 -0.07 12◦ 0.53 0.12 0.66 0.13 0.80 0.14 0.54 0.04 24◦ 0.89 0.40 1.09 0.47 1.19 0.50 0.94 0.35 155 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−6 −4−2 02 46 Flap Cycle (t/T)V ertical Force Coeffic ient (C Z) 0Hz4Hz6Hz8Hz10Hz Downstroke Upstroke (a) CZ, θ = 0◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−0.2 −0.15−0.1 −0.050 0.050.1 0.150.2 Flap Cycle (t/T)Ho rizonta l Force Coeffic ient (C X) 0Hz4Hz6Hz8Hz10Hz Downstroke Upstroke (b) CX, θ = 0◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−6 −4−2 02 46 Flap Cycle (t/T)V ertical Force Coeffic ient (C Z) 0Hz4Hz6Hz8Hz10Hz Downstroke Upstroke (c) CZ, θ = 12◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−2 −1 0 1 2 3 Flap Cycle (t/T)Ho rizonta l Force Coeffic ient (C X) 0Hz4Hz6Hz8Hz10Hz Downstroke Upstroke (d) CX, θ = 12◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−6 −4−2 02 46 Flap Cycle (t/T)V ertical Force Coeffic ient (C Z) 0Hz4Hz6Hz8Hz10Hz Downstroke Upstroke (e) CZ, θ = 24◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−2 −1 0 1 2 3 Flap Cycle (t/T)Ho rizonta l Force Coeffic ient (C X) 0Hz4Hz6Hz8Hz10Hz Downstroke Upstroke (f) CX, θ = 24◦. Figure 4.1: Variation of total vertical force coefficient (CZ) and total horizontal force coefficient (CX) from CFD data (rigid wing). 156 4.1.2 Analysis of the 3-D Rigid Wing FlowField over a Flap Cycle (Q-Criterion) In order to gain a deeper understanding of the cause behind these trends, it is important to investigate the three-dimensional flow across the wing through the course of one flap cycle. Throughout the rest of this section, the figures and tables will focus on one case with a flap frequency of 4 Hz and a pitch angle of 24◦. To visualize the vortical flow across the wing, the second invariant of the velocity gradient tensor (Q-criterion) of the flowfield is calculated. Physically it represents the influence of pure shear to that of pure rotation on a given fluid element and can be used to determine the presence of coherent vortical structures within the flowfield [185]. Figure 4.2 shows iso-surfaces of the second invariant of the velocity gradient tensor, (Q-criterion), colored with a magnitude of vorticity contour at various intervals of the flap cycle for this case. Figures 4.2(a) and 4.2(b) represent the initial stages of the downstroke portion of the flap cycle where t = 0 and t = T/8 respectively. It is here that the initial formation of the leading edge, root and tip vortices upon the upper surface of the wing can be seen. Note that in Fig. 4.2(b), along the bottom surface, there is spanwise vorticity forming from the root to the tip of the wing. In Fig. 4.2(c) (t = T/4, mid-downstroke) and Fig. 4.2(d) (t = 3T/8) when the wing is undergoing the second half of the downstroke, there is a large accumulation of vortical structures on the wing’s upper surface that has begun to encompass much of the wing chord. Figures 4.2(e) and 4.2(f) show the flowfield at t = T/2 and t = 5T/8, respectively 157 and depict the end of the downstroke and the beginning of the wing’s upstroke. Globally, there is a noticeable decrease in the vortices about the wing’s upper surface. The amount of vorticity present has become smaller in magnitude and less coherent in comparison to that at t = T/4. This is due to the change in flow direction around the root and tip associated with the change in flap direction. Figure 4.2(g) characterizes the mid-upstroke of the flap cycle. Notice that vortical structures are now beginning to form on the underside of the wing and there are significantly fewer vortical structures present in comparison to the mid-downstroke position (Fig. 4.2(c)). The reason for this is the lower angles of attack experienced by the wing during the upstroke. Lastly, Fig. 4.2(h) shows the latter half of the upstroke portion of the flap cycle. There is a clear formation of a vortex core at the wing tip as the majority of the vorticity has developed on the bottom surface of the wing. For the duration of the upstroke, while there was a clear formation of a tip vortex, the root vortex remained less coherent, suggesting that the difference between the high pressure and the low pressure side is small at the wing root compared to the tip. The plots of Q-Criterion over time in Figs. 4.2 provide a global view of how the 3-D flow around the wing was evolving over the course of a single flap cycle. In this case, the Q-criterion allows for a qualitative examination of where vorticity was being generated and showed how this vorticity migrated within the flowfield. The Q-Criterion does show the overall complexity of the flow and the evolution of the flow structures. However, the iso-surface plots are not able to give a clear view of the vortical flow at different locations of the flapping wing so that it could be related to the forces produced by the wing. 158 (a) t = 0/T . (b) t = T/8. (c) t = T/4. (d) t = 3T/8. (e) t = T/2. (f) t = 5T/8. (g) t = 3T/4. (h) t = 7T/8. Figure 4.2: Iso-surfaces of Q-Criterion colored with vorticity magnitude contour at various instances of the flap cycle (rigid wing). 159 4.1.3 Rigid Wing Analysis of the 3-D Flowfield Over a Flap Cycle (Spanwise Vorticity Contours) In order to provide a clear view of the change in the flowfield over the wing through the duration of a flap cycle, Fig. 4.3 illustrates the vorticity magnitude at specific spanwise locations along the wing. Within each figure, the chordwise sections shown, from root to tip, correspond to the 30%, 50%, 70%, and 90% spanwise locations respectively. In each image the root is on the left and the tip of the wing is on the right. In Fig. 4.3(a), as the wing initiates the downstroke of the flap cycle (t = 0), the initial formation of the LEV at the various spanwise locations is observed. In Fig. 4.3(b) (t = T/8), there is a significant growth in the size of the LEV along the span as the wing begins to progress through the downstroke. At the 30% span location, the LEV is still attached to the upper surface of the airfoil. However, on the more outboard sections of the wing, portions of the LEV are beginning to separate from the wing. This is due to an increase in the vertical component of velocity along the span induced by flapping motion of the wing. This increase in vertical velocity leads to an increase in the local angle of attack and resultant velocity from root to tip. Figure 4.3(c) represents the mid-downstroke (t = T/4) and is the point during the flap cycle where the largest magnitudes of resultant velocity and angle of attack occur. At this time, there is a high level of vorticity and separated flow on the upper surface of the wing creating a large region of low pressure. The highest value of vorticity and flow separation takes place at the mid-downstroke and explains the 160 peaks in the magnitude of the CZ and CX in Figs. 4.1(e) and 4.1(f) respectively. In Fig. 4.3(d) (t = 3T/8), there is still a significant region of separated, vortical flow on the topside of the wing. However, from Figs. 4.1(e) and 4.1(f), we notice a decrease in the magnitude of the aerodynamic force generated by the wing. This is because the wing is decelerating as it approaches the end of the downstroke and the magnitudes of the local angles of attack and resultant velocities along the span begin to decrease. At t = T/2, shown in Fig. 4.3(e), the wing has completed the downstroke portion of the flap cycle and is now beginning the upstroke. A moderate level of separated flow is present along the entire wing span, but the extent and magnitude has noticeably decreased between Fig. 4.3(e) (t = T/2) and Fig. 4.3(d) (t = 3T/8). This is because at the beginning of the upstroke, the velocity due to the wing’s motion is essentially zero and the magnitude of the resultant velocity is equal to the free-stream velocity. In Fig. 4.3(f), as the wing begins to progress through the upstroke, there is a shift in the direction of the resultant velocity vector that results in a decrease in the local angle of attack along the span. There is also a significant decrease in the magnitude of vorticity across the wing’s upper surface, and a slight amount of vorticity is beginning to form on the bottom surface of the wing. The differences in vorticity about the upper and lower surfaces result in the near zero values for CZ and CX seen at t = 5T/8 as shown in Figs. 4.1(e) and 4.1(f). At the mid-upstroke (t = 3T/4), the resultant velocity vector is again at its maximum magnitude along the span but is acting slightly downward due to the negative component of the vertical velocity vector. From Fig. 4.3(g), it can be 161 seen that on the inboard portion of the wing, much of the vorticity is located on the topside of the wing because of the positive angle of attack, but going outboard, a majority of the vorticity is located on the underside of the wing owing to the negative angle of attack induced by flapping. Therefore, the flapping motion causes a change in sign of the angle of attack and hence a change in direction of the lift force from root to tip. This is responsible for the small peak in negative vertical force (negative CZ) and positive propulsive thrust (negative CX) at that instance in the flap cycle (refer to Figs. 4.1(e) and 4.1(f)). Lastly, in Fig. 4.3(h), the magnitude of the resultant velocity vector has decreased and there is little to no separation along the wing span as it begins to decelerate. While the level of vorticity present is low in comparison to the other instances of the flap cycle, the majority is located on the wings upper surface and is responsible for the low magnitude of positive vertical force (positive CZ) and negative propulsive thrust (negative CX) at t = 7T/8. It is interesting to note that in Fig. 4.3(h), toward the tip, the vorticity is concentrated on the bottom surface of the airfoil, but toward the root on the upper surface. This again implies a difference in the sign of the local angle attack along the wing span even during this instance of the upstroke. To explain the changes in the flowfield over the rigid wing at particular instances of the flap cycle, Fig. 4.4 shows the vorticity contours at specific spanwise locations along the wing. The chordwise sections shown correspond to the 30%, 50%, 70%, and 90% spanwise locations respectively. Figures 4.4(a), 4.4(c), 4.4(e) and 4.4(g) show the vorticity at the mid-downstroke (t = T/4), and Figs. 4.4(b), 4.4(d), 4.4(f) and 4.4(h) correspond to the mid-upstroke (t = 3T/4). The mid-downstroke and 162 (a) t = 0/T . (b) t = T/8. (c) t = T/4. (d) t = 3T/8. (e) t = T/2. (f) t = 5T/8. (g) t = 3T/4. (h) t = 7T/8. Figure 4.3: Vorticity magnitude along the span at various instances of flap cycle (rigid wing). 163 mid-upstroke are highlighted because at these moments in the flap cycle, the wing is subject to the highest instantaneous local resultant velocities, and it exhibits peaks in aerodynamic force production. During mid-downstroke a majority of the vorticity is concentrated on the upper surface of the wing. At the 30% span location (Fig. 4.4(a)) there is a noticeable formation of vorticity present at the leading edge of the wing. Near the wing root the LEV remains attached to the wing surface, but begins to separate toward the wing tip. In Fig. 4.4(c) and 4.4(e) it can be clearly seen that outboard of the 50% span location, the LEV has detached from the wing’s upper surface and shed. At the 90% span location (Fig. 4.4(g)), the flow becomes significantly chaotic due to tip effects which cause the leading edge vortex to break down and separate. At the mid-upstroke (t = 3T/4), the resultant velocity vector is again at its maximum magnitude along the span but is acting slightly downward due to the negative component of the vertical velocity vector. It can be seen that on the inboard portion of the wing (Fig. 4.4(b)), much of the vorticity is located on the upper surface of the wing because of the positive angle of attack. At the 50% span location (Fig. 4.4(d)), there is almost no vorticity present on the wing surface due to the near zero local induced angle of attack at that section of the wing. Outboard on the wing, (Figs. 4.4(f) and 4.4(h)), a majority of the vorticity is located on the lower surface of the wing because of the negative angle of attack induced by flapping. Therefore, the flapping motion causes a change in sign of the effective angle of attack and hence a change in direction of the lift force from root to tip. This is responsible for the small peak in negative vertical force (negative CZ) and positive propulsive thrust (negative CX) at that instance in 164 the flap cycle. To analyze the degree of three-dimensionality present in the flow, Figs. 4.5(a) and 4.5(b) show chordwise cuts along the length of the wing at the 20% chord location from the LE, during the mid-downstroke (t = T/4) and mid-upstroke (t = 3T/4) of the flap cycle respectively. Here, a qualitative representation of the vorticity and fluid flow along the wing’s span is presented. In Figs. 4.5(a) and 4.5(b), there is a concentration of vorticity at the wing root and tip due to the presence of the root and tip vortices. As expected, the tip vortex is of greater size and magnitude in comparison to the root vortex. Consequently, toward the wing root and wing tip, there was a high magnitude of spanwise velocity on the order of the freestream velocity. However, the magnitude of spanwise velocity at locations further from the root and tip is significantly less. At these spanwise locations, the direction of the velocity vector can be assumed to lie in the 2-D chordwise plane. Also, from the vorticity contours at mid-downstroke, (shown in Fig. 4.5(a)) it can be inferred that both the root and tip of the wing (and therefore every spanwise section of the wing) are at a positive angle of attack because the root and tip vortices have opposite sign vorticity. This is because the flow is leaking from the lower side to the upper side of the wing at both root and tip. However, at the mid-upstroke (Fig. 4.5(b)) the vorticity has the same sense for both root and tip vortices because the angle of attack is positive at the root and negative at the tip. Thus far, the variation of CZ and CX of the entire wing over the course of one flap cycle (shown in Fig. 4.1) was explained using the flowfield around the wing (Figs. 4.2, 4.3, 4.4, and 4.5). It has also been shown that the flow on the 165 Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (a) Mid-downstroke, 30% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (b) Mid-upstroke, 30% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (c) Mid-downstroke, 50% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (d) Mid-upstroke, 50% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (e) Mid-downstroke, 70% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (f) Mid-upstroke, 70% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (g) Mid-downstroke, 90% span. Normalized X−coordinate, x/c Norma lized Z −coord inate, z/c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (h) Mid-upstroke, 90% span. Figure 4.4: Vorticity contours with overlaid velocity along the span: left side, mid-downstroke and right side, mid-upstroke (rigid wing), 4 Hz, θ = 24◦. 166 (a) X Vorticity, Mid-Downstroke. (b) X Vorticity, Mid-Upstroke. Figure 4.5: X-Vorticity magnitude along the span from CFD data (rigid wing), 4 Hz, θ = 24◦. flapping wing is predominantly 2-D except at the wing root and tip. Therefore, in the following section, the vertical and propulsive thrust forces produced by the flapping wing are investigated further using the computed force coefficients at various spanwise sections of the wing. Note that, even though these are 2-D force coefficients, they are obtained from the full 3-D analysis of the flapping wing and hence would include the 3-D effects as well. The following analysis and results will be focused on the mid-downstroke (t = T/4) and mid-upstroke (t = 3T/4) because it is at these instances where the peaks in the magnitude of CZ and CX occur during a flap cycle. 4.1.4 Rigid Wing 2-D Analysis of Wing (Sectional Vertical and Hori- zontal Force Coefficients along the Span) To illustrate the variation in the sectional aerodynamic force along the wing, Figures 4.6 shows the variation in the 2-D sectional values of Cz and Cx along the span for both the mid-downstroke and mid-upstroke respectively for the same case (24◦ pitch angle and 4 Hz flapping frequency). In Figs. 4.6(a) and 4.6(b), it is interesting to point out that the sectional values 167 Table 4.2: Chart of UP, UT, V and φ along the span (rigid wing). Span, y/b 0.30 0.40 0.50 0.60 0.70 0.80 0.90 UP, (m/s) 0.80 1.07 1.34 1.604 1.872 2.139 2.407 UT, (m/s) 3.00 3.00 3.00 3.00 3.00 3.00 3.00 V, (m/s) 3.10 3.185 3.284 3.402 3.536 3.685 3.846 φ, (deg) 14.971 19.622 24.021 28.138 31.961 35.491 38.736 of Cz and Cx are large in magnitude, especially toward the tip of the wing. Note the forces shown in the aforementioned figures were non-dimensionalized by the freestream velocity (V∞). Due to the flapping kinematics of the wing, there is a significant vertical velocity component perpendicular to the freestream direction. Thus the resultant velocity experienced at a spanwise section of the wing is defined as: ~V = ~Vwind + ~Vflapping = √ U2P + U 2 T (4.1) Here, UP is the velocity parallel to flap plane and UT is the velocity component normal to flap plane. The relative inflow angle (or induced angle of attack) is φ = tan−1 ( UP UT ) . Table 4.2 presents the values of UP, UT, V and φ for the 4 Hz flapping case at the mid-downstroke. Here it is shown that the resultant velocity is significantly higher than V∞ and the flap induced velocity φ is as high as 40◦ near the tip (significantly beyond the static stall angle). Figure 4.7 shows the spanwise distribution of Cz and Cx along the span, non- dimensionalized by the calculated resultant velocity, V , at the corresponding spanwise locations, for the 4 Hz flapping case. Note that the magnitude of the force coefficient 168 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Nondimensional Span Location, y/b Vert ical Forc e Co effic ient, C z Mid DownstrokeMid Upstroke (a) Vertical force coefficient, Cz, θ = 24◦. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.5 0 0.5 1 1.5 2 Nondimesional Span Location, y/b Hor izon tal F orce Coe fficie nt, C x Mid DownstrokeMid Upstroke (b) Horizontal force coefficient, Cx, θ = 24◦. Figure 4.6: Rigid wing sectional Cz and Cx, at mid-downstroke and mid-upstroke, normalized by V∞. 169 values has now decreased considerably, especially toward the tip of the wing where the flap induced velocity is higher. It can be seen that the values of Cz and Cx along the span during the downstroke are both positive, resulting in the wing producing positive lift and negative propulsive thrust (such as drag) along the entirety of the wing span. However, that same trend is not seen in the upstroke. During the upstroke, a change in the direction of the Cz and Cx force along the span occurs. Closer to the root, the wing is producing a positive lifting force (positive Cz) and a drag producing force (positive Cx). However, in the outboard section (y/b > 0.45), there is a shift, and the wing is producing a negative lifting force (negative Cz) and a positive propulsive thrust (negative Cx). The main difference between Fig. 4.6 and Fig. 4.7 is the magnitude of the values for Cz and Cx. Normalization by the freestream velocity causes a higher lift coefficient value (compared to the resultant velocity) because the induced velocity caused by the flapping motion is not taken into account. In the literature, the calculation of the aerodynamic forces acting on a flapping wing has been an issue. Previous studies have made assumptions about the flowfield attempting to estimate the angle of attack. Usherwood and Ellington [186,187] assumed a “triangular” downwash distribution to estimate the angle of attack, where the downwash was assumed to vary linearly with the span. Birch and Dickinson [89] calculated an angle of attack using the velocity vectors that were acquired using PIV, where they assumed that the mean orientation of the flow vectors near the lower surface of the wing represented the freestream velocity. In the current study, the forces are normalized by the resultant velocity, the angle of attack is determined by Eqn. 4.2, and for clarity, the lift and 170 drag coefficients are decomposed into vertical and horizontal force coefficients. The flowfields in Figs. 4.3(c), 4.3(g), and 4.4, which represent the vorticity present during the mid-downstroke and mid-upstroke respectively, can be used to provide some information pertaining to the changes seen in the sectional aerodynamic coefficients along the span. In Fig. 4.3(c), mid-downstroke, much of the vorticity is close to the upper surface of the wing creating an increase in Cz and Cx. After the 50% span location (Fig. 4.4(c)), there is a reduction in both Cz and Cx due to the leading edge vortex becoming detached and shedding downsteam. This decreasing trend in the magnitude of Cz and Cx continues until approximately the 70% span location (4.4(e)) where the tip vortex begins to have a more significant influence over the flow. The tip vortex generates a region of low pressure near the tip resulting in the slight increase in lift and drag seen at the 80% and 90% span locations. For the mid-upstroke, there is a more linear, decreasing trend in the vertical and horizontal force coefficients from root to tip. From Fig. 4.3(g), at the 30% span location (see Fig. 4.4(b)), most of the vorticity is concentrated on the wings upper surface resulting in a net positive vertical force (positive Cz) and a net drag force (positive Cx). Near the 50% span location (e.g., Fig. 4.4(d)), the direction of net vertical and horizontal force switches due to changes in the angle of attack across the span. Figure 4.4(h) shows the 90% span location where all of the vortical flow is along the underside of the wing resulting in a peak in the magnitude of negative lifting force (negative Cz) and positive propulsive thrust (negative Cx) along the span. To explain the change in direction of the vertical and horizontal forces, it is 171 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.5 0 0.5 1 1.5 2 2.5 3 Nondimensional Span Location, y/b Ver tica l Fo rce Co effi cien t, C z Mid Downstroke Mid Upstroke (a) Vertical force coefficient, Cz, θ = 24◦. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Nondimensional Span Location, y/b Ho rizo nta l Fo rce Co effi cien t, C x Mid Downstroke Mid Upstroke (b) Horizontal force coefficient, Cx, θ = 24◦. Figure 4.7: Rigid wing Cz and Cx, at mid-downstroke and mid-upstroke, normalized by the resultant V. 172 important to look at the local angle of attack, α, at the various spanwise sections. From equation 4.2, α is defined as: α = θ + φ (4.2) Equation 4.2 shows that α is a function of the pitch angle, θ, and the flap induced angle, φ. During the upstroke, φ is negative and smaller in magnitude than θ toward the root and greater than θ toward the tip. This results in a positive α at the wing root and negative α at the wing tip, thus explaining the difference in the direction of aerodynamic force along the wing span and the sign difference in Cz and Cx about the downstroke and upstroke as a whole. 4.1.5 2-D Analysis of Rigid Wing (Sectional Lift and Drag Coefficients along the Span) The presence of the flap induced velocity, UP, is responsible for a change in the local α along the span and is also thought to be a factor in the increase of the force coefficients. The vertical velocity component, UP, not only changes the magnitude and direction of the resultant velocity vector, but also the lift and drag force vectors. For static wing cases, where α is solely based on the pitch angle of the wing, (i.e., φ = 0), the magnitudes and directions of Cz and Cx correspond directly to those of Cl and Cd respectively. However, for a flapping wing, φ is non-zero and the magnitudes and directions of Cz and Cx are not the same as Cl and Cd. Figure 4.8 illustrates the difference in the magnitude and direction of the vectors for Fz, Fx, lift and drag. 173 (a) Downstroke. (b) Upstroke. Figure 4.8: Schematics of force and velocity vectors for the airfoil cross section. From the schematics in Figs. 4.8(a) and 4.8(b), Cz and Cx can be expressed as a function of Cl and Cd using the following equations: Cz = Cl cos(φ) + Cd sin(φ) (4.3) Cx = −Cl sin(φ) + Cd cos(φ) (4.4) In equations 4.3 and 4.4, the magnitudes of Cz and Cx only depend on the magnitudes of Cl and Cd and φ. The magnitude of Cl and Cd are governed by the local α at a given spanwise section. Figure 4.9 shows a plot of Cl and Cd along the span for the mid-downstroke and mid-upstroke at a pitch angle of 24◦. In Fig. 4.9, a positive or negative angle of attack corresponds to a positive or negative value of Cl. In contrast, the value of Cd is always positive, regardless of the sign of α. For the case where θ is equal to 24◦, near the wing root, where the resultant velocity is approximately equal to the freestream value, the magnitudes of Cl and Cd for the mid-downstroke and mid-upstroke are not comparable to the static wing case at 24◦. The flowfield analysis presented 174 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.5 0.1 0.7 1.3 1.9 2.5 Nondimensional Span Location, y/b Lift an d D rag Co effi cie nts , C l & C d Cl, downstroke Cl, upstroke Cd, downstroke Cd, upstroke −11.5°−7.9° −4.1°−0.1°4.4° 43.6° 48.0° 52.1° 55.9° 59.5° 62.7° −14.7° α = 38.9° α = 9.0° Figure 4.9: Lift coefficient, Cl, and drag coefficient, Cd, along the span, at mid- downstroke and mid-upstroke. before indicated a dynamic stall type of phenomena on the wing due to the unsteady variation of angle of attack at very high reduced frequency. This could be the reason for the delayed stall and unusually high Cl values. Additionally, these large values of Cl and Cd cannot be predicted using a simple quasi-steady analysis. Note that during the mid-downstroke of the flap cycle, the value of Cl never exceeds a value of about 2.4. However, when comparing the maximum values of Cz for the aforementioned case (Fig. 4.6(a)) it can be seen that the maximum sectional value is approximately 3. This is because of the contribution of Cd to Cz resulting from the large induced angle (φ). As shown in Figs. 4.10(a) and 4.10(b), as the magnitude of φ increases toward the wing tip, there is a rotation of the resultant velocity vector (V ) by the same angle φ, upward during downstroke and downward during upstroke. In turn, 175 the lift and drag vectors also rotate by the same angle. Therefore, the magnitudes of Cz and Cx at any span location would depend not only on the magnitudes of lift and drag but also their direction (given by φ) as given by equations 4.3 and 4.4. In Fig. 4.10(a), a schematic of the rotation of the wing lift and drag vectors, the magnitude of φ increases toward the wing tip at the mid-downstroke and results in an increase in the local angle of attack from root to tip. This causes the lift vector to tilt forward and drag vector to tilt upward while moving from root to tip. This would cause the lift to contribute more towards propulsive thrust, while the drag would contribute to positive vertical force and produce less negative thrust in the outboard sections. From Fig. 4.9 the magnitude of Cl decreases toward the wing tip; however, in Fig. 4.7(a) there is an increase in the magnitude of Cz from the 30% to 50% span locations. This increase in Cz is due to the combination of the vertical component of the lift force as well as a vertical component of the drag force causing a net increase in Cz greater than the magnitude of the Cl alone. From Fig. 4.10(b), the variation in φ along the span during the upstroke causes the lift vector to act primarily downward and slightly forward as one goes closer to the wing tip. This rotation of the lift vector results in part of the lift force acting to propel the wing forward (positive propulsive force) as well as pull the wing down (negative vertical force). From Fig. 4.9, it can be seen that the value of Cd during the upstroke remains relatively constant, but the magnitude of the lift vector increases from approximately the 50% span location to the wing tip. This results in a larger component of force acting in the negative x-direction as opposed to the positive x-direction creating a net propulsive thrust (negative Cx). 176 (a) Downstroke. (b) Upstroke. Figure 4.10: Rotation of the lift and drag vectors toward the wing tip. 4.1.6 Section Summary In this section, a combination of 3-D and spanwise 2-D flowfield, along with instantaneous force coefficient values were used to understand the underlying physics behind a wing undergoing pure flapping kinematics at a fixed pitch angle. A key insight gained from this analysis on rigid, low-aspect ratio wings undergoing pure flap kinematics is that much of the positive vertical force and drag is produced during the downstroke portion of the flap cycle while a majority of the propulsive thrust and negative vertical force is generated during the upstroke portion of the flap cycle. The vertical and propulsive thrust distribution on a flapping wing depends on the magnitude and direction of the lift and drag along the wing span. Significantly large lift and drag coefficients and stall delay to very high angles of attack (much higher than the static stall value) were observed along the span of the wing because of the unsteady angle of attack variation at very high reduced frequency leading to a dynamic stall type of phenomena. It was also concluded that pure flap kinematics without active wing pitching 177 would not be able to produce both positive vertical force and propulsive thrust for a flapping wing, thus rendering it impractical for a flapping MAV design; however, to improve the net propulsive thrust of a rigid wing configuration, a negative pitch angle should be utilized during the downstroke and a positive pitch angle during the upstroke. Additionally, to improve the net positive lift, there needs to be an asymmetry in the pitch variation between the upstroke and downstroke. Further improvement in the efficiency of force production would require varying the pitch angle of the wing sections along the span along with the temporal variation (or dynamic twist). This is because different spanwise sections are subjected to different local angles of attack over time, and adjusting the pitch angle along the span would give more control over the angle of attack at which different spanwise sections operate. This could be achieved using a torsionally flexible wing with a suitable spanwise flexibility distribution with the chordwise c.g. behind the elastic axis. If the flexibility and mass distribution are properly tailored, the inertial and aerodynamic forces acting on the wing can be utilized to passively twist the wing without the need for any additional pitching actuator. The next section, Section 4.2, contains a discussion of the comparison of the rigid and flexible wing results. 4.2 Rigid and Flexible Wing Comparison 4.2.1 Aerodynamic Force Production In the previous section, the results from the 3-D CFD model were systematically analyzed. In this section, the results of the coupled CFD-CSD numerical solver will 178 be utilized to better understand the flow physics about a flapping wing in forward flight. Specifically, the influence of flexibility on the production of lift (FZ) and propulsive thrust (FX) will be presented. Throughout this study, the structural characteristics of the wing were kept constant with the goal being to compare the results of a structurally characterized flexible wing to those of a rigid wing. Figure 4.11 shows the variation in the vertical force coefficient (CZ) and the horizontal force coefficient (CX) over the course of one flap cycle for the case with a flap frequency of 6 Hz and a pitch angle of 12◦. Figures 4.11(a) and 4.11(b) also show the time-averaged CZ and CX values during the flap cycle. Note that a positive value for CX represents a drag producing force while a negative value for CX represents a propulsive thrust producing force. For concurrent rigid wing tests [132], it was observed that a majority of the positive vertical force production occurs during the downstroke of the flap cycle while a majority of the propulsive thrust production occurs during the upstroke. This behavior can clearly be seen in Fig. 4.11(a) and 4.11(b). A flexible wing also exhibits similar behavior with regard to vertical force production. Due to a positive geometric pitch angle, there is asymmetry in the force production between the downstroke and upstroke. The asymmetry in vertical force production results in an average CZ of approximately 0.66 for both the rigid and flexible wings. In contrast, in Fig. 4.11(b), there is a significant difference in the net horizontal force produced between the rigid and flexible wings. During the upstroke of the rigid wing (t/T > 0.5), as stated earlier, a majority of the propulsive thrust is produced, and this remains true for the flexible wing also. However, unlike the rigid wing 179 case, the drag force generated is relatively low for a majority of the downstroke of the flexible wing (1 < t/T < T/2). The slight amount of drag producing force generated is of a significantly lower magnitude than that of the rigid wing during the downstroke. When comparing the net horizontal force production over a complete flap cycle, the net CX is 0.13 for the rigid wing and -0.19 for the flexible wing clearly showing a large propulsive thrust for the flexible wing. From the rigid wing analysis, it was determined that the use of rigid wings undergoing pure flap kinematics in MAV designs would not be practical because they are not able to simultaneously produce a significant net lifting and net propulsive force during a flap cycle. However, from Fig. 4.11, introducing flexibility allows for the wing to produce both a net lifting force and net propulsive force over one flap cycle. These findings are motivation to further explore the aerodynamic force production on a flexible wing and, more importantly, how flexibility influences the flow around the wing and the mechanisms responsible for aerodynamic force production. 4.2.2 3-D Analysis of Wing (Aerodynamic Force Production over a Flap Cycle) To gain an understanding of how flexibility influences the performance of a flapping wing in forward flight, it is important to analyze the production of aerodynamic forces over time. Figure 4.12 shows the variation in vertical force coefficient (CZ) and horizontal force coefficient (CX) over the course of one flap cycle for the flap frequencies and pitch angles tested during experimentation. Again, note 180 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−4 −3 −2 −1 0 1 2 3 4 5 6 Flap Cycle (t/T) Ver tica l Fo rce Co effi cien t (C Z ) FlexibleFlexible: Mean CZ RigidRigid: Mean CZ Downstroke Upstroke (a) Variation of CZ. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−1 −0.5 0 0.5 1 1.5 Flap Cycle (t/T) Ho rizo nta l Fo rce Co effi cien t (C X ) FlexibleFlexible: Mean CX RigidRigid: Mean CX Downstroke Upstroke (b) Variation of CX. Figure 4.11: Rigid and flexible wing comparison of predicted vertical (FZ) and horizontal force (FX) for single flap cycle with flap frequency of 6 Hz and θ0 = 12◦. 181 Table 4.3: Chart of CZ and CX for flap frequency and pitch angle (flexible and rigid (in parenthesis) wing respectively). 4 Hz 6 Hz 8 Hz 10 Hz θ0 CZ CX CZ CX CZ CX CZ CX 0◦ 0.01 (-0.02) -0.09 (0.00) 0.00 (0.01) -0.32 (-0.01) 0.01 (-0.01) -0.82 (-0.04) 0.02 (-0.02) -1.69 (-0.07) 12◦ 0.53 (0.53) -0.02 (0.12) 0.70 (0.66) -0.19 (0.13) 1.03 (0.80) -0.63 (0.14) 1.42 (0.54) -1.43 (0.04) 24◦ 0.89 (0.89) 0.22 (0.40) 1.25 (1.09) 0.19 (0.47) 1.74 (1.19) -0.12 (0.50) 2.50 (0.94) -0.91 (0.35) that a positive CX value represents a drag producing force and a negative CX value represents a propulsive thrust producing force. Figures 4.12(a) and 4.12(b) show the change in CZ and CX, respectively, for an initial pitch angle of 0◦ at flap frequencies of 4, 6, 8, and 10 Hz. In Fig. 4.12(a), it can be seen that as the flapping frequency is increased, there is an increase in the magnitude of the instantaneous CZ. This is due to an increase in the resultant velocity and induced angle of attack about the wing caused by the flapping kinematics. The positive lifting force generated during the downstroke of the wing is nearly identical to the downward vertical force generated during the upstroke of the wing, resulting in a net vertical force near zero. However, in Fig. 4.12(b), the wing is generating propulsive thrust, (negative CX), throughout almost the entire flap cycle. The magnitude of the peak instantaneous CX value increases with frequency and the variation of instantaneous CX is nearly symmetrical about the upstroke and downstroke. For the cases where the initial pitch angle is increased (12◦ and 24◦), there is an asymmetry between the downstroke and upstroke resulting in a difference in the aerodynamic force production over a flap cycle. In Figs. 4.12(c) and 4.12(e), the magnitude of the positive CZ generated during the downstroke is larger than the magnitude of the negative CZ generated during the upstroke. This creates a net 182 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−10 −5 0 5 10 Flap Cycle (t/T) Verti cal F orce Coef ficien t (C Z) 4Hz6Hz8Hz10Hz (a) CZ, θ0 = 0◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−4 −3 −2 −1 0 1 Flap Cycle (t/T) Horiz ontal Forc e Co efficie nt (C X ) 4Hz6Hz8Hz10Hz (b) CX, θ0 = 0◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−10 −5 0 5 10 Flap Cycle (t/T) Verti cal F orce Coef ficien t (C Z) 4Hz6Hz8Hz10Hz (c) CZ, θ0 = 12◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−4 −3 −2 −1 0 1 Flap Cycle (t/T) Horiz ontal Forc e Co efficie nt (C X ) 4Hz6Hz8Hz10Hz (d) CX, θ0 = 12◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−10 −5 0 5 10 Flap Cycle (t/T) Verti cal F orce Coef ficien t (C Z) 4Hz6Hz8Hz10Hz (e) CZ, θ0 = 24◦. 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1−4 −3 −2 −1 0 1 Flap Cycle (t/T) Horiz ontal Forc e Co efficie nt (C X ) 4Hz6Hz8Hz10Hz (f) CX, θ0 = 24◦. Figure 4.12: Variation of total vertical force coefficient (CZ) and total horizontal force coefficient (CX) from coupled CFD-CSD data (flexible wing). 183 positive lifting force over the course of a flap cycle, which increases in magnitude with increased pitch angle. In Figs. 4.12(d) and 4.12(f), the increase in pitch angle generates an asymmetry in horizontal force generation between the downstroke and upstroke not seen in the 0◦ pitch angle cases. For the 12◦ pitch angle cases shown in Fig. 4.12(d), the peak magnitude of CX during the downstroke is significantly less than that seen during the upstroke. However, for a majority of the downstroke, the wing is thrust producing throughout all the flap frequencies shown except for the two lowest flapping frequency cases. At the 4 Hz and 6 Hz flapping frequency, the wing is producing a drag force, (positive CX), through the entire downstroke and generates thrust during the upstroke. In Fig. 4.12(f), which presents the 24◦ pitch angle cases, the asymmetry in thrust production between the upstroke and downstroke is even more significant. For all but the 10 Hz flapping frequency, the wing is producing drag during the downstroke and producing thrust during the upstroke. However, these cases, except for the 4 Hz and 6 Hz flapping frequencies, still produce a net negative CX (i.e., propulsive thrust producing) force over the duration of a flap cycle. Table 4.3 provides a list of the calculated net CZ and CX values over a flap cycle for different of flap frequencies and pitch angles presented in Fig. 4.12. Note that for each case within Table 4.3, the time-averaged coefficient value for the flexible wing is listed first, followed by the time-averaged coefficient value for the rigid wing in parentheses. From Table 4.3, it is clear that a majority of the flap frequencies and pitch angles tested for the flexible wing generated both net lifting force and propulsive thrust. One way to assess the performance of the wing is to analyze the lift-to-drag 184 ratio (L/D). Lift-to-drag ratio is a measure of aerodynamic efficiency of the wing. Conventionally, lift-to-drag ratio is the amount of lifting force generated by a wing, or vehicle, divided by the drag force. In this case CZ and CX can be used in replace of lifting force and drag force. However, for this particular flapping wing configuration, L/D can be positive or negative depending on whether the wing is producing drag or propulsive force. The optimal CZ/CX value would need to be negative and close to -1 implying that the configuration is simultaneously producing similar magnitudes of lift and propulsive thrust. From Table 4.3, this occurs at the 10 Hz, 12◦ case which yields a CZ/CX value of -0.993. However, this is an approximation from observing the current data set. Applying an optimization module to the simulation would give a more exact value for the combination of flap frequency and pitch angle for obtaining both high lift and propulsive force during flight. 4.2.3 3-D Rigid and Flexible Wing Flowfield Comparison (Q-Criterion) Unlike the rigid wing analysis, flexible wings operating at particular flapping frequencies and initial pitch angles are able to produce a net positive lifting force and a net propulsive thrust over a flap cycle. It is important to investigate the three dimensional flowfield around the flexible wing in order to acquire a deeper understanding of what is causing the trends seen in Fig. 4.12. Throughout the rest of this section, the case of a wing subjected to a flap frequency of 6 Hz and initial pitch angle of 12◦ will be highlighted. Figure 4.13 shows iso-surfaces of the second invariant of the velocity gradient tensor, (Q-Criterion), for both the rigid and flexible wings 185 at specific instances of the flap cycle. Note that Figs. 4.13(a) and 4.13(c) represent iso-surfaces of Q-criterion for the flexible wing while Figures 4.13(b) and 4.13(d) are representative of those for the rigid wing. The wing root is on the left on each image and the wing tip is toward the right. Figures 4.13(a) and 4.13(b) depict the mid-downstroke of the wing (t = T/4). This is the point in the flap cycle where there is a peak in the magnitude of instantaneous lifting force. There is a significant formation of vortical structures about the top surface of the wing near the leading edge with both root and tip vortices also forming on the top surface of the wing. At this point, the flexible wing has noticeably deformed with the wing pitch angle along the span having decreased in comparison to the rigid wing. While the rigid wing is showing larger leading edge vortices in comparison to the flexible wing, the vortical structures begin to break down towards the wing root and tip as seen in Fig. 4.13(b). At the more outboard sections of the wing, the leading edge and tip vortices are beginning to detach and move downstream. However, in Fig. 4.13(a), the vortical structures at the leading edge, while smaller in magnitude, are more coherent and extend slightly greater across the wing span before breaking down. The root and tip vortices are of a slightly smaller size and appear to be closer in proximity to the wing’s surface. Figures 4.13(c) and 4.13(d) illustrate the vortical flow about the wing during the mid-upstroke. This case was chosen because there is a peak in propulsive thrust production by both the rigid and flexible wings. Due to the change in flap direction, with regard to the downstroke, a majority of the formation of vortical structures is occurring on the underside of the wing. This suggests that during the upstroke of 186 the wing, the angle of attack along the wing span has become negative. Both the tip and root vortices have begun to form on the underside of the wing suggesting that the high pressure side of the wing is now on the wing’s upper surface. The flexible wing has deformed creating an increase in the geometric pitch angle along the span in comparison to that of the rigid wing. While both wings exhibit a relatively large volume of incoherent vortical structures in the immediate wake of the wing, the rigid wing shows a slightly higher amount in comparison to the flexible wing. The Q-Criterion plots of the flexible and rigid wings in Figs. 4.13(a)–4.13(d) provide a means of qualitatively comparing the three dimensional flow about the wing at points in the flap cycle, which display peak aerodynamic force production. It is seen that the introduction of flexibility was globally affecting the fluid flow and evolution of vorticity across the wing over time. However, plots of the iso-surface of Q-Criterion do not allow one to easily distinguish the difference in the formation of vortical structures at various locations along the span because it is difficult to relate the deformations seen along the span of the flexible wing to the flowfield. Therefore, chordwise slices of the flowfield (vorticity contours) are analyzed. 4.2.4 3-D Rigid and Flexible Wing Flowfield Comparison (Spanwise Vorticity Contours) Figures 4.14(a) and 4.14(b) show the superimposed wing shapes at the mid- downstroke and mid-upstroke respectively. The wing shape of the rigid wing is colored blue while the flexible wing is colored gray. In Figure 4.14(a) for the mid-downstroke, 187 (a) t = T/4, downstroke flexible. (b) t = T/4, downstroke rigid. (c) t = 3T/4, upstroke flexible. (d) t = 3T/4, upstroke rigid. Figure 4.13: Rigid and flexible wing comparison of iso-surfaces of Q-Criterion colored with vorticity magnitude at various instances of the flap cycle. 188 (a) Mid-downstroke. (b) Mid-upstroke. Figure 4.14: Superimposed rigid and flexible showing deformations, θ0 = 12◦, Ω =6 Hz, (40◦ flap amplitude, V∞ = 3 m/s). it can be seen that the pitch angle along the span of the flexible wing is significantly different than that of the rigid wing. The difference in pitch angle can be seen to increase along the span with the local pitch angle on the flexible wing decreasing and even becoming negative from root to tip. Conversely, in Figure 4.14(b) for the mid-upstroke, the difference in pitch angle between the flexible and rigid wing is less severe but still shows an increasing trend along the span with the magnitude of the local pitch angle increasing from root to tip. Figures 4.15 and 4.16 illustrate the vorticity contours, with overlaid velocity field, at various spanwise locations along the wing. These images more clearly display the formation of vortical flow over the wing and provide a means to better understand what role flexibility plays in enhancing the aerodynamic performance of a flapping wing. Note that in each contour, the chord-wise sections correspond to the 30%, 50%, 70%, and 90% span locations from wing root to wing tip. Fig- ures 4.15(a), 4.15(c), 4.15(e), and 4.15(g) represent spanwise vorticity contours for the 189 flexible wing while Figures 4.15(b), 4.15(d), 4.15(f), and 4.15(h) are representative of those for the rigid wing in the downstroke. The same can be said about the corresponding contours in Fig. 4.16 for the upstroke. The contours shown in Fig. 4.15 correspond to the mid-downstroke (t = T/4). This is the point in the flap cycle where the wing is subject to the highest instantaneous local resultant velocities and induced angles of attack. For both the rigid and flexible wings, a majority of the vorticity is concentrated on the topside of the wing. For the flexible wing, at the 30% span location (Figs. 4.15(a) and 4.15(b)) there is a noticeably lower magnitude of vorticity present in comparison to the same span location on the rigid wing. This because while both sections are experiencing the same local resultant velocity, the deformation of the flexible wing has caused the pitch angle to decrease and consequently, lower the magnitude of the local angle of attack relative to that of the rigid wing. Comparing the accumulation of vorticity across the span of the rigid and flexible wings: for the rigid wing, the concentration of vorticity present is remarkably larger than that of the flexible wing for the 50% span location. However, outboard of the 50% span, the leading edge vortex along the rigid wing has begun to separate and shed. On the other hand, the leading edge vortex about the flexible wing remains attached further along the wing span until the 90% span location. At this point, the presence of the tip effects causes the leading edge vortex to breakdown and separate. Figure 4.16 shows the spanwise vorticity magnitude during the mid-upstroke of the flap cycle (t = 3T/4). Note that from Fig. 4.11, it is at this instance during the flap cycle that both the flexible and rigid wings generate the peak instantaneous 190 Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (a) 30% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (b) 30% span, rigid. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (c) 50% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (d) 50% span, rigid. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (e) 70% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (f) 70% span, rigid. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (g) 90% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−0.6 −0.4−0.2 00.2 0.40.6 0.81 (h) 90% span, rigid. Figure 4.15: Chordwise vorticity contour with overlaid velocity field at 30%, 50%, 70%, and 90% span locations at mid-downstroke (flexible and rigid wing). 191 Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (a) 30% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 (b) 30% span, rigid. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (c) 50% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 (d) 50% span, rigid. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (e) 70% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 (f) 70% span, rigid. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8−0.6 −0.4−0.2 00.2 0.40.6 (g) 90% span, flexible. Normalized X−coordinate, x/c Norma lized Z −coord inate, z /c −0.5 0 0.5 1 1.5 2−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 (h) 90% span, rigid. Figure 4.16: Chordwise vorticity contour with overlaid velocity field at 30%, 50%, 70%, and 90% span locations at mid-upstroke (flexible and rigid wing). 192 propulsive thrust and negative lifting force. For both wings, toward the wing root, there is nearly zero vorticity present because the magnitude of the local induced angle of attack is relatively small (< 10◦). Further outboard along the span, the concentration of vorticity along the underside of the wing increases due to an increase in the magnitudes of the local resultant velocity and local induced angle of attack. At the 70% and 90% span locations, the flexible wing exhibits a slightly smaller amount of vorticity about the leading edge in comparison to the rigid wing. At the outboard sections of the rigid wing, there is a slight amount of flow separation present, particularly in Fig. 4.16(f). However, the leading edge vortices on the flexible wing remain attached across the span. This is because the deformation of the flexible wing along the wing span causes an increase in the local pitch angle resulting in a decrease of the magnitude of the local induced angle of attack. The decrease in the local induced angle of attack helps to mitigate detachment but still allows for the formation of a substantial leading edge vortex. 4.2.4.1 2-D Analysis of Wing (Sectional CZ and CX along the Wing Span)) Until this point, the focus was on comparing the flowfield of the flexible wing to that of the rigid wing to better understand how introducing flexibility changes the fluid flow about the flapping wing. Qualitatively, it is seen that flexibility and wing deformations cause a change in pitch angle along the span, the formation and location of vorticity along the span including its magnitude is affected. In this section, the variation in the sectional aerodynamic coefficients along the wing span will be investigated at the mid-downstroke (t = T/4) and mid-upstroke (t = 3T/4) as these 193 Table 4.4: Chart of UP, UT, V and φ along the span (flexible wing). Span, y/b 0.30 0.40 0.50 0.60 0.70 0.80 0.90 UP, (m/s) 1.20 1.60 2.01 2.41 2.81 3.21 3.61 UT, (m/s) 3.00 3.00 3.00 3.00 3.00 3.00 3.00 V, (m/s) 3.23 3.40 3.61 3.85 4.11 4.39 4.69 φ, (deg) 21.86 28.14 33.76 38.74 43.10 46.93 50.27 are the instances where the peaks in magnitude of Cz and Cx occur. Figures 4.17(a) and 4.17(b) depict the variation of sectional Cz and Cx, respec- tively, along the span for both the flexible and rigid wings at the mid-downstroke and mid-upstroke. It should be noted that prior to this section, the vertical and horizontal forces produced by the entire wing were non-dimensionalized by the free-stream ve- locity (V∞). However, at the mid-downstroke and mid-upstroke, the velocity induced by the flapping kinematics is at its largest magnitude and contributes significantly to the resultant velocity at a particular spanwise section of the wing. Thus, the sectional aerodynamic coefficients presented are normalized using equation 4.1 For reference, Table 4.4 contains the values corresponding to the magnitude of UT , UP , V and the flap induced angle φ at the various spanwise locations analyzed for a flap frequency of 6 Hz. From Table 4.4, toward the wing root, the influence of the flap induced velocity on the local resultant velocity is relatively small causing only a 7% increase in its magnitude. However, toward the wing tip the flap induced velocity is much higher resulting in an over 50% increase in the magnitude of the resultant velocity in comparison to the free-stream velocity. Given the significant difference between the calculated resultant velocity and free-stream velocity, it is important to consider the 194 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−2 −1 0 1 2 3 4 Nondimensional Span Location, y/b Ver tica l Fo rce Co effi cien t, C z Mid−Downstroke: Flexible Mid−Upstroke: Flexible Mid−Downstroke: Rigid Mid−Upstroke: Rigid (a) Vertical force coefficient, Cz. 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Nondimensional Span Location, y/b Ho rizo nta l Fo rce Co effi cien t, C x Mid−Downstroke: Flexible Mid−Upstroke: Flexible Mid−Downstroke: Rigid Mid−Upstroke: Rigid (b) Horizontal force coefficient, Cx. Figure 4.17: Sectional Cz and Cx, at mid-downstroke and mid-upstroke, 6 Hz, θ0 = 12◦. 195 induced flap velocity when calculating the sectional force coefficients. Shown in Fig. 4.17(a), during the mid-downstroke, the value of CZ along the wing span is consistently positive for both the flexible and rigid wings, corresponding to the wings producing a positive CZ at that point in the flap cycle. Also, the magnitude of CZ along the span during the mid-downstroke is noticeably greater toward the wing root. Note that when normalized by the local resultant velocity, the value of sectional CZ of the flexible wing along the span during the mid-downstroke is relatively constant at approximately 2.5 until the 70% span location. Further outboard of the 70% span location, it is thought that the tip vortex has a more significant effect and influences the aerodynamic force produced in that region. Conversely, during the mid-upstroke, the value of sectional CZ along the span is consistently negative for the rigid and flexible wing, resulting in the wing producing a negative instantaneous CZ at that point in the flap cycle. Also, the rigid wing consistently produces a negative sectional lifting force of greater magnitude than the flexible. For the flexible wing, unlike during the mid-downstroke, the magnitude of negative CZ is increasing along the span. This increase in magnitude continues until the 70% span location where again, the tip vortex is said to begin having a more significant influence over the aerodynamic force produced. In Fig. 4.17(b), during the mid-downstroke of the flexible wing, the sign of the sectional CX value changes from positive to negative after the 50% span location. However, for the rigid wing, the value of CX is consistently positive and is of a much higher magnitude. This shows that the rigid wing is producing a significant amount of drag along the span, but on the flexible wing, the more inboard sections 196 are producing small drag forces while those sections toward the wing tip generate propulsive thrust. These findings suggest that for the flexible wing, there is a change in the direction of the sectional resultant force vector along the wing span. The sum of the propulsive thrust generated by the outboard wing sections is slightly greater than the total drag force generated by the inboard wing sections resulting in an instantaneous, net propulsive thrust force produced at the mid-downstroke. During the mid-upstroke of the flap cycle, both wings are thrust producing (negative CX) along the entirety of the span. Toward the wing root, the sectional CX is of like magnitudes, but moving along the span, the flexible wing produces a notably greater amount of sectional propulsive thrust. For the flexible wing, similar to the trend of sectional CX during the mid-downstroke, the magnitude of CX increases along the span until approximately the 70% span location. After this point, there is a slight deviation from this trend and a decrease in the magnitude of sectional CX due to the presence of the tip vortex. During the upstroke, the flexible wing produces significantly higher propulsive thrust compared to the rigid wing in the outboard sections. 4.2.4.2 2-D Analysis of Wing (Sectional Cl and Cd along the Wing Span) To explain the change in magnitude and direction of the vertical and horizontal forces along the wing span, it is important to discuss the flap induced velocity, UP, at the various span locations. The vertical velocity component influences the magnitude of the resultant velocity vector, V , as well as the flap induced angle, φ. Consequently, 197 UP has an effect on the local angle of attack and the direction of the lift and drag force vectors. Using the classical definitions of lift and drag, the lift force vector is perpendic- ular to the resultant velocity while the drag force vector is parallel to the incoming resultant vector. For a static wing where the wing is not flapping, φ = 0 and the lift and drag force vectors, respectively, coincide with the vertical and horizontal forces. However, if the wing is undergoing flapping kinematics, (i.e., φ 6= 0), the vertical and horizontal forces do not coincide with the lift and drag forces, and can be described again by Equations 4.3, 4.4, and 4.2. Note that for the flexible wing, at a given time instance, the values of θ and φ vary considerably along the wing span. The value of φ is mainly dependent on the prescribed flapping kinematics of the wing and can easily be calculated analytically. However, the value of θ is influenced by the wing flexibility and mass distribution as well as the inertial and aerodynamic loads acting on the wing. Figure 4.18 shows the variation of geometric pitch angle along the span during the mid-upstroke and mid-downstroke at an initial pitch angle of 12◦ for all the flap frequencies tested. We see that as the flapping frequency is increased, the change in magnitude of the pitch angle is significantly increased in comparison to the initial pitch angle of 12◦. During the downstroke of the wing, the direction of inertial loads acting on the wing creates a nose-down wing twist. However, the opposite occurs during the upstroke and the geometric pitch angle increases along the span due to a nose-up wing twist. Note that the observed twist variation along the span was not linear. 198 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−25 −20 −15 −10 −5 0 5 10 15 20 Nondimensional Span Location, y/b Pitc h A ngl e (°) 4Hz, downstroke 4Hz, upstroke 6Hz, downstroke 6Hz, upstroke 8Hz, downstroke 8Hz, upstroke 10Hz, downstroke 10Hz, upstroke Figure 4.18: Variation of geometric pitch angle along the span during the mid- upstroke and mid-downstroke at an initial pitch angle of 12◦ for 4, 6, 8, 10 Hz flap frequencies. Figure 4.19 shows the change in geometric pitch angle at the mid-downstroke and mid-upstroke for the 6 Hz flapping frequency cases. At this frequency, the wing begins to experience significant aerodynamic and inertial loads causing a large degree of deflection along the wing span. For all the cases, there is a significant difference in the magnitude and direction of geometric twist between the mid-downstroke and mid-upstroke. Also, while the geometric twist is nearly symmetrical for the 0◦ initial pitch angle case, as expected, the 12◦ and 24◦ root pitch cases exhibit an asymmetry in the change in wing twist between the mid-downstroke and mid-upstroke. Figure 4.20 shows the plot of Cl and Cd at the mid-downstroke and mid-upstroke along the span of the flexible wing. First, it should be noted that the sign of the angle of attack, (positive or negative), is indicative of a positive or negative lift coefficient. The value of the drag coefficient is always positive, regardless of the sign of the angle 199 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−15 −10 −5 0 5 10 15 Nondimensional Span Location, y/b Pitc h A ngl e (°) 0°, downstroke 0°, upstroke 12°, downstroke 12°, upstroke 24°, downstroke 24°, upstroke Figure 4.19: Variation in geometric pitch angle at the mid-downstroke and mid- upstroke (6 Hz). of attack. Second, the magnitudes of Cd and, more importantly, Cl are significantly different from those one would expect in comparison to a static wing operating at a high angle of attack. This is due, in part, to the unsteady effects generated by the flapping kinematics of the wing, and it would not be possible to predict these unsteady values for the aerodynamic coefficients using a basic quasi-steady analysis. However, the inclusion of non-constant translation, rotation and acceleration could improve the predictive capability of such a quasi-steady model in comparison to our results. In Fig. 4.20, from the curves representing the mid-downstroke, notice that the magnitude of Cl is greater than the magnitude of Cd except at the 90% span location. This becomes an important fact when explaining the flexible wing’s ability to produce thrust during both the upstroke and downstroke. Toward the wing root, 200 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−2 −1.1 −0.2 0.7 1.6 2.5 r/R C l, C d Cl, downstrokeCl, upstrokeCd, downstrokeCd, upstroke Figure 4.20: Lift coefficient, Cl, and drag coefficient, Cd, along the span at mid- downstroke and mid-upstroke. the magnitude of the vertical force coefficient in Fig. 4.17(a) during the downstroke is slightly greater than the value of Cl in Fig. 4.20. Due to the presence of a positive induced flap angle, a component of the drag force is directed in the positive vertical direction enhancing the lifting force. However, the component of Cl directed in the negative x-direction is smaller in magnitude than the component of Cd directed in the positive x-direction, resulting in a drag producing sectional force. However, toward the wing tip where the magnitude of φ is increased, the rotation of the resultant velocity vector is also increased. The change in the rotation of the resultant velocity vector rotates the direction of the lift and drag vectors. This causes the lift vector to contribute more toward propulsive thrust and the drag vector to contribute more toward vertical lifting force. At this point, the component of Cl directed in the negative x-direction is greater than the component of Cd directed in 201 the positive x-direction. This results in the simultaneous production of propulsive thrust and lifting force. 4.2.5 Section Summary A combination of 3-D and 2-D flowfield analysis, in conjunction with spanwise sectional force analysis, was used to further understand the underlying physics behind a flexible flapping wing undergoing pure flap kinematics at a given initial pitch angle. One of the key insights gained is that unlike a rigid flapping wing in pure flap, a flexible wing is capable of producing substantial net positive lifting force and net propulsive thrust over the course of a flap cycle. Similar to a rigid wing, a majority of the positive vertical force is generated during the downstroke of the flap cycle, and downward vertical force is produced during the upstroke. However, for a flexible wing, propulsive thrust is generated during both the downstroke and upstroke with a majority of the propulsive thrust being generated during the upstroke. The enhanced aerodynamic force production of the flexible wing can be at- tributed its ability to dynamically twist while flapping. During flapping, the inertial and aerodynamic loads passively deform the wing creating a nose-down blade twist during the downstroke and a nose-up blade twist during the upstroke. The spatial and temporal variations of the pitch angle produces a more favorable local induced angle of attack along the span. It was seen that for the flexible wing, the decrease in the magnitude of the local induced angle of attack allows the LEV to remain attached to the wing surface further along the wing span. A rigid wing with constant 202 pitch along the span undergoing pure flap kinematics does not allow for control of the local induced angle of attack and exhibits decreased aerodynamic efficiency, especially toward the wing tip. For the flexible wing, the optimal CZ/CX occurs at the 10 Hz, 12◦. The vertical lift and propulsive thrust along the wing span is largely influenced by the magnitude and direction of the sectional lift and drag vectors. The presence of unsteady effects allowed for the development of significantly large lift and drag coefficients. The large magnitudes of the induced angles of attack along the wing span lead to a dynamic stall type of phenomena, contributing to the increase in lift and drag coefficient. 203 Chapter 5: Concluding Remarks Systematic experimental and computational analyses were conducted to un- derstand the flow physics and force production on rigid and flexible flapping wings in forward flight. The experimental study involved detailed flowfield measurements using PIV and direct force measurements. The computational study was performed using a coupled CFD-CSD based aeroelastic analysis. This included a RANS solver coupled to a large deflection geometrically exact structural model (using MBDyn). Avian-based fixed-pitch flap wing kinematics were implemented on a rectangular flex- ible wing in a wind tunnel at a wind speed of 3 m/s. Time-resolved planar PIV and time-averaged force measurements were used to validate the 3-D coupled CFD-CSD results. Time-averaged PIV flowfield results at mid-downstroke, LEV circulation, and total force measurements were compared to the coupled CFD-CSD showing good agreement. With an understanding of the limitations of the computational model, analyses of the wing at various instances in the flap cycle were presented to describe the flow physics that govern the lift and propulsive force production. 204 5.1 Conclusions from Rigid Wing Analysis Pure flap tests were carried out on a rigid wing. From force measurements, it was observed that with a root pitch angle, θ = 0◦, positive propulsive force (FX) was generated. However, the average force in vertical direction (FZ) was zero for V∞ = 3m/s. On increasing the root pitch angle, FZ increased but FX became negative. Upon analyzing the flowfield, the following specific conclusions are drawn from the rigid wing study: 1. The PIV on the rigid wing showed that the LEV was observed to develop mostly in phase along the span of the wing, peaking in strength between the y/b = 0.5 and 0.6 span locations and decreasing quickly in strength at the root and tip of the wing. The root and tip of the wing were shown to trail dominant vortices into the wake. Estimated using PIV, the LEV size and strength does vary at different span locations along the wing. This is because an increase in vorticity would also mean that there is an increase in net circulation around the wing, thereby causing increased lift when the LEV is present. When the LEV is shed, the vorticity becomes less concentrated and, the LEV diffused rapidly, especially further outboard on the wing. 2. CFD simulation results resolved the flow features present for the rigid flapping wing, particularly at the inboard section for low flap frequencies and low pitch angles. However, aperiodicity further outboard at higher frequencies and 205 high pitch angles decreases simulation accuracy. With this, the turbulence, convection, and dissipation parameters for vortex behavior can be modified to more closely reflect vortex strength and trajectory at higher angles of attack and flap frequencies. Overall, CFD results showed good agreement with the experimental data. 3. Flow over the rigid wing was highly susceptible to changes in the flap-induced angle resulting from the flapping motion and variations in reduced frequency, which manifested in the predicted airloads. It is observed that during the wing stroke, the wing does experience stall conditions when the effective angle of attack is very high. However, at the stall angles, Cl, and Cd were significantly higher than the static values, indicating a dynamic stall type of phenomena. It is important to note that for the LEV, the behavior (e.g., formation, residency time, shedding and downstream convection) forms the core problem of dynamic stall. As this vortex convects away from the wing surface, it does alter the chordwise pressure distribution, and it produces transient forces fundamentally different from those in static stall. 4. Based on the computational analysis, the spanwise flow component was not significant except near the wing tip, and therefore, most of the vertical force and propulsive thrust produced could be explained using the magnitude and direction of the sectional lift and drag forces acting on the wing. This was was because of the tip vortex and its interaction with the LEV. 5. Most of the upward vertical force was produced during the downstroke and 206 positive propulsive thrust during the upstroke for the pure flap wing kinematics used. This shows the need for appropriate temporal and spanwise pitch modulation of the wing along with flapping to produce positive vertical force and propulsive thrust during the entire flap cycle. Temporal and spatial pitch modulation could keep the vortex attached while maintaining a propulsive force. It was seen that both positive propulsive and vertical forces are required to optimized forward flight; therefore, pure flap of a rigid wing may not be a viable configuration. With this, the effects of wing flexibility were investigated to understand its effect on those parameters and the flowfleid. 207 5.2 Conclusions from the Comparison of Rigid and Flexible Wings From the rigid wing test, it was seen that the achievement of the best efficiency from a wing, both temporal and spanwise variations of pitch angle are required. Therefore, experiments and analyses were performed for a realistic (highly anisotropic with different flexibilities in the spanwise and chordwise directions) flapping wing. The conclusions from the rigid and flexible wing comparison are summarized: 1. The PIV showed that for the flexible wing, the LEV remains attached to the wing up to the mid-downstroke inboard region; however, it sheds further outboard, especially at higher flap frequency and higher angle of attack. The LEV was observed to peak in strength between the 50% and 70% span locations and develop mostly in-phase along the span of the wing. The LEV decreases in sectional strength moving from mid-span toward the root or toward the tip of the wing. The root and tip of the wing were also observed to trail small vortices into the wake, with the portion of the LEV near the root/tip detaching more quickly in the stroke than further inboard. 2. The downstream wake and LEV swirl velocity profile shows a comparable result for both the coupled CFD-CSD and PIV, which shows that the salient features (LEV, trailing edge vortices, root vortex, tip vortex) of the flow are captured by the CFD-CSD methodology. 3. The effects of wing flexibility are shown in LEV size and strength at the different span locations along the wing, which is reflected comparably in the coupled 208 CFD-CSD and PIV results. In both cases, there is an increase in vorticity while the wing is in its downstroke causing an increase in net circulation around the wing, thereby increasing lift when the LEV is present. When the LEV is shed, the vorticity becomes less concentrated, and it diffuses rapidly as seen further outboard on the wing. 4. Coupled CFD-CSD simulation results resolved the flow features present for the flapping wing. However, flow separation further outboard at higher frequencies and high initial pitch angles decreased simulation accuracy. With this, modifi- cations to the turbulence model can be made to more accurately predict the convection and dissipation of the vortices, particularly further away from the airfoil surface. It would be expected that the predicted solution would more closely reflect the vortex strength and trajectory seen in the experiment at higher angles of attack and flap frequencies. Overall, the computational results showed good agreement with the experimental data. 5. For both rigid and flexible wings, the flow over the wing was highly susceptible to changes in induced angle of attack. This was a direct result of increasing the flapping frequency (increased reduced frequency). Throughout the wing stroke, highly vortical flow was generated by the wing, and stall conditions were observed at high effective angles of attack. However, at the observed stall angles, Cl, and Cd were significantly higher than the expected values for a static wing under the same conditions. This indicates that the flapping motion creates a dynamic stall type phenomena where increased lift is observed. 209 6. Wing flexibility increased the aerodynamic performance in comparison to a rigid wing. Flexibility, more specifically caused dynamic twisting, allowing for temporal and spatial variations in the geometric pitch angle throughout the flap cycle. The variations in pitch angle decreased the magnitude of the induced local angle of attack along the span. In comparison to a rigid wing, the LEV of the flexible wing remained attached over more of the wing span resulting in enhanced aerodynamic force production. 7. Unlike rigid wings undergoing pure flap kinematics at a fixed pitch, flexible wings are capable of simultaneously producing positive lifting force and positive propulsive thrust over a given flap cycle. For the flexible wing, the optimal CZ/CX occurs at the 10 Hz, 12◦. This highlights that the inclusion of dynamic twist provides a means of passive temporal and spanwise pitch modulation in a flapping wing MAV design. In order to further improve aerodynamic efficiency, it is necessary to aeroelastically tailor the wing to create a more optimal twist distribution along the wing span throughout the flap cycle. There are several benefits to using a flexible wing. In general, the flexible wing mechanism is less complex and lighter. Also, by employing dynamic twisting, better pitch angles are obtained for the wing. This leads to greater forces at all spanwise sections as compared rigid wings. Finally, the flexible wing itself is also much lighter than the rigid wing. These results show the beneficial effects of the appropriate amount of wing flexibility. A well tailored flexible wing may be a suitable configuration for use on an MAV. 210 5.3 Major Contributions The main contribution of the current work was to extend the limited body of work on 3-D flapping wings using experiments and an aeroelastic solver to study the performance and flow physics of realistic flexible flapping wing MAVs. The major components of the study are the flowfield experiment, airloads measurements, wing characterization, and the coupled aeroelastic solver. To this end, experiments on realistic flexible wings were executed to determine the physics behind generating vertical and propulsive forces in forward flight. The wings were then structurally characterized to generate a comprehensive data set, currently unavailable in literature, which can be used for the validation of aeroelastic analyses. Also, this experimental data set was used to compare with the predictive results. The major contributions of this work are: 1. An experimental forward flight data was acquired for a set of structurally well-characterized, flexible, lightweight flapping MAV wings. Flow diagnostic techniques were developed to permit high-fidelity measurements of the flowfield. Airloads were also measured for comparison with the aeroelastic solver. 2. This work quantified and identified the vortex strength for the rigid and flexible wing undergoing passive wing deformation. Also, the wing motion and the aeroelastic deformations of a flexible flapping-wing MAV are quantified. It was demonstrated that forward flight performance can be improved, simultaneously producing positive lifting force and positive propulsive thrust over a given flap 211 cycle, by aeroelastically tailoring flexible wings. 3. Flowfield characteristics are correlated to airloads adding to a lacking body of work at the low Re scale (i.e., flow visualization, PIV, characterization of wing properties, analysis of deflections, natural frequencies, airloads). Refinement and validation of comprehensive CFD-CSD analysis capable of simulating realistic anisotropic wings undergoing large non-linear deformations were carried out systematically. The tactical impact is that the insights gained through experimental data and analysis can reduce the amount of time and effort spent on experimental trial and error platform development which translates to aerodynamic improvements of MAVs that will be used to keep military and rescue personnel safe. 212 5.4 Suggestions for Future Work This dissertation is a step towards the development of a fully autonomous, hover-capable flapping wing micro air vehicle. While this research coverers techniques for experimentation and analysis, there are several approaches along which further research can be carried out. 1. More time averaged experimental studies can be conducted, especially in the area of PIV measurements of the spanwise flow (stereoscopic measurements). Targeted variations in Reynolds numbers, aspect ratio, wing planform, and different wing kinematics can be tested. The present study only utilized a limited selection of wing coverings (wing membrane materials) for wing construction; however, more advanced light-weight materials can also be used on the wings. Also assessment of the tip vortex could be used to quantify the lift and drag coefficients to the development of the the tip vortex and LEV. Future work can be expanded with the analysis of the vortex to include a larger range of conditions for LEV evolution and local Reynolds number effects. Also, the amount of wing flexibility is important; however, this analysis only assessed the performance with one set of wing properties. Future studies can explore the point of diminishing returns as it pertains to wing bending and torsional stiffness. 2. Most vortex models have been tested for higher Reynolds numbers; however, a low Reynolds number vortex model would be useful for MAV analysis in the 213 future. LEV analysis can continue to characterize its structure and breakdown toward developing a new applied vortex model. Currently, the value in this work is that it offers a mode for flowfield comparison between the PIV and CFD-CSD analysis. In lieu of pressure and force measurements, which are difficult to execute at MAV scale, the assumptions made in the control volume method for estimating the horizontal forces from the flow field are valid. However, estimating the forces from an unsteady flowfield is an “active area of research” and the control volume method should be explored in more detail with comparison to direct force measurements. 3. To reduce the number of experiments and simulations time, an optimization module can be added to the simulation. Optimization would allow for the maximizing or minimizing of a system by systematically choosing input values from within an allowed set and computing the value of the function. This would aid in better understanding design trade-offs and would aid in finding an optimal aerodynamic design point without more experimentation and simulation cost. The simulation can be expanded to model an entire vehicle. 4. The present work only investigated simplified wing kinematics, avian-type wing kinematics in pure flap. It would be useful to also expand the analysis to more complex avian wing kinematics (e.g., active asymmetric pitching) as well as hover wing kinematics. MAVs are very susceptible to wind gusts so gust tolerance analysis using stability derivatives would be helpful in determining the agility and maneuverability of a vehicle. 214 5.5 Acknowledgment of Sponsors This work was supported, in part, by the Army Collaborative Technology Alliance (CTA) for the Micro Autonomous Systems and Technology (MAST) program Contract No. W911NF-08-2-0004 with Mr. Chris Kroninger (ARL-VTD) as technical monitor and by the L-3 Communications Graduate Research Fellowship. The author would also like to acknowledge with thanks the support of the NACME Alfred P. Sloan Foundation Graduate Research Scholarship, the American Helicopter Society’s Vertical Flight Fellowship, and Aerospace Corporation Graduate Fellowship Awards. 215 Bibliography [1] Pines, D. J., and Bohorquez, F., “Challenges Facing Future Micro-Air-Vehicle Development,” Journal of Aircraft, Vol. 43, No. 2, March-April 2006, pp. 290– 305. [2] Conn, A., Burgess, S., R., H., and Ling, C. S., “From Natural Flyers to the Mechanical Realization of a Flapping Wing Micro Air Vehicle,” IEEE International Conference on Robotics and Biomimetics, Kunming, China, December 17-20 Proceedings of the 2006 2006. [3] Richter, C., and Lipson, H., “Untethered Hovering Flapping Flight of a 3D- Printed Mechanical Insect,” Proceedings of the Alife XII Conference, Odense, Denmark, January 2010. [4] McMichael, J. M., and Francis, M. S., “Micro Air Vehicles - Toward a New Dimension in Flight,” DARPA , July 1996. [5] Piekarski, B., Dologoski, T., and Kroninger, C., “Army Collaborative Tech- nology Alliance (CTA) for the Micro Autonomous Systems and Technology (MAST) Annual Report,” Army Research Laboratory , January 2008-2015. [6] AeroVironment “AeroVironment,” AeroVironment Website, Available from: http://www.aerovironment.com/, Accessed on August 24, 2014. [7] Shyy, W., Lian, Y., Tang, J., Viieru, D., and Liu, H., Aerodynamics of Low Reynolds Number Flyers Cambridge University Press, New York, NY, 2008. [8] Brown, R. H. J., “The Flight of Birds, II. Wing Function in Relation to Flight Speed,” Journal of Experimental Biology, Vol. 30, No. 1, March 1953, pp. 90–103. [9] Azuma, A., The Biokinetics of Flying and Swimming Springer-Verlag, Tokyo, 1992. 216 [10] Zdunich, P., Bilyk, D., McMaster, M., Loewen, D., Delaurier, J., Kornblun, R., Low, T., Stanford, S., and Holeman, D., “Development and Testing of the Mentor Flapping-Wing Micro Air Vehicle,” Journal of Aircraft, Vol. 44, No. 5, September-October 2007, pp. 1701–1711. [11] Tarascio, M., and Chopra, I., “Design and Development of a Thrust Augmented Entomopter: An Advanced Flapping Wing Micro Hovering Air Vehicle,” 59th American Helicopter Society Forum, Pheonix, AZ, May 2003, AIAA 2007-662. [12] AeroVironment “DARPA Awards AeroVironment Phase II Contract Exten- sion for Nano Air Vehicle Development Program,” AeroVironment Press release, Available from: http://www.avinc.com/downloads/ NAVPRLong- DARPAV4.doc.pdf, Accessed on August 8, 2014, June 30 2009. [13] Harmon, R. L., “Aerodynamic Modeling of a Flapping Membrane Wing Using Motion Tracking Experiments,” Masters of science, University of Maryland, August 2008. [14] Petricca, L., Ohlckers, P., and Grinde, C., “Micro and Nano-Air Vehicles: State of the Art,” International Journal of Aerospace Engineering, Article ID 214549. [15] Mayo, D., and Leishman, J., “Comparison of Hovering Efficiency of Rotating Wing and Flapping Wing Micro Air Vehicles,” Journal of American Helicopter Society, Vol. 55, 2010. [16] McMasters, J. H., and Henderson, M. L., “Low Speed Single Element Airfoil Synthesis,” Technical Soaring, Vol. 6, No. 2, 1980, pp. 1–21. [17] Anderson, J. D., Fundamentals of Aerodynamics McGraw-Hill, New York, NY, 2007. [18] Ellington, C. P., van den Berg, C., Willmott, A. P., and Thomas, A. L. R., “Leading-Edge Vortices in Insect Flight,” Nature, Vol. 384, No. 6610, December 1996, pp. 626–630. [19] Leishman, J. G., Principles of Helicopter Aerodynamics, 2nd ed., Cambridge University Press, Cambridge, United Kingdom, 2006. [20] Birch, J. M., and Dickinson, M. H., “Spanwise flow and the attachment of the leading-edge vortex on insect wings,” Nature, Vol. 412, No. 6848, August 2001, pp. 729–733. [21] Sane, S. P., “Review: The aerodynamics of insect flight,” The Journal of Experimental Biology, Vol. 206, No. 663, January 2003, pp. 4191–4208. [22] Dickinson, M. H., Lehmann, F.-O., and Sane, S. P., “Wing Rotation and the Aerodynamic Basis of Insect Flight,” Science, Vol. 284, No. 5422, 18 June 1999, pp. 1954–1960. 217 [23] Singh, B., and Chopra, I., “Wing Design and Optimization for a Flapping Wing Micro Air Vehicle,” American Helicopter Society Annual Forum, Baltimore, MD, 7-10 June 2004. [24] Wang, J., Chen, J.-K., and Liao, S., “An explicit solution of the large defor- mation of a cantilever beam under point load at the free tip,” ScienceDirect: Journal of Computational and Applied Mathematics, Vol. 212, No. (2008), January 2008, pp. 320–330. [25] Azuma, A., Okamoto, M., and Yasuda, K., “Aerodynamic Characteristics of Wings at Low Reynolds Number,” In Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, T. J. Mueller, Ed., Vol. 195 of Progress in Astronautics and Aeronautics. AIAA, Reston, Virginia 20191-4344, 2001, Ch. 17, pp. 341–398. [26] Mueller, T. J., and DeLaurier, J. D., “Aerodynamics of Small Vehicles,” Annual Review of Fluid Mechanics, Vol. 35, 2003, pp. 89–111. [27] Ol, M. V., “Unsteady Aerodynamics for Micro Air Vehicles,” NATO Re- search and Technology Organization Contractor Report TR-AC/323 (AVT-149) TP/332,Final Report of Task Group AVT-149,, 2010. [28] Flynn, A. M., “Flapping Nano Air Vehicles,” MicroPropulsion MicroPropulsion Technical Report Number 100, Ames Research Center and U.S. Army Air Mobility Research and Development Laboratory, Moffett Field, CA, January 2010. [29] Broeren, A. P., and Bragg, M. B., “Unsteady Stalling Characteristics of Thin Airfoils at Low Reynolds Number,” In Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, T. J. Mueller, Ed., Vol. 195 of Progress in Astronautics and Aeronautics. American Institute of Aeronautics and Astro- nautics, Inc., 1801 Alexander Bell Drive, Reston, Virginia 20191-4344, 2001, Ch. 10, pp. 191–213. [30] Pelletier, A., and Mueller, T. J., “Low Reynolds Number Aerodynamics of Low-Aspect-Ratio, Thin/Flat/Cambered-Plate Wings,” Journal of Aircraft, Vol. 37, No. 5, September-October 2000, pp. 825–832. [31] Selig, M. S., Guglielmo, J. J., Broern, A. P., and Giguere, P., “Experiments on Airfoils at Low Reynolds Numbers,” 34th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 15-18 January 1996, no. AIAA in 96-0062. [32] Mueller, T. J., and DeLaurier, J. D., “An Overview of Micro Air Vehicle Aerodynamics,” In Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, T. J. Mueller, Ed., Vol. 195 of Progress in Astronautics and Aeronautics. American Institute of Aeronautics and Astronautics, Inc., 1801 Alenander Bell Drive, Reston, Virginia 20191-4344, 2001, Ch. 1, pp. 191–213. 218 [33] Null, W., and Shkarayev, S., “Effect of Camber on the Aerodynamics of Adaptive Wing Micro Air Vehicles,” 2nd AIAA Flow Control Conference, Portland, OR, July 2004. [34] Ifju, P. G., Stanford, B., M., S., and Albertani, R., “Analysis of a Flexible Wing Micro Air Vehicle,” 25th AIAA Aerodynamic Measurement Technology and Ground Testing Conference, San Francisco, CA, 5–8 June 2006, Vol. AIAA 2006-0000. [35] Lian, Y., and Shyy, W., “Numerical Simulations of Membrane Wing Aerody- namics for Micro Air Vehicle Applications,” Journal of Aircraft, Vol. 42, No. 4, July-August 2005, pp. 865–873. [36] Lian, Y., and Shyy, W., “Laminar-Turbulent Transition of a Low Reynolds Number Rigid or Flexible Airfoil,” 36th AIAA Fluid Dynamics Conference and Exhibit, San Francisco, California, 5-8 June 2006, Vol. AIAA 2006-3051. [37] Stanford, B., Viieru, D., Albertani, R., Shyy, W., and Ifju, P., “A Numerical and Experimental Investigation of Flexible Micro Air Vehicle Wing Deformation,” 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, January 9–12 2006. [38] Guerreiro, N. M., and Hubbard, J. E., “Flight Testing of a Model Aircraft with Trailing-Edge Flaps Optimized for Lift Distribution Control,” AIAA Atmospheric Flight Mechanics Conference and Exhibit, Honolulu, Hawaii, August 2008, AIAA-2008-6709, Ed. [39] Grasmeyer, J., and Keennon, M., “Development of the Black Widow Micro Air Vehicle,” AIAA, 2001. [40] Ifju, P., Flexible-Wing-Based Micro Air Vehicles, Vol. C.H.M.Jenkins, ed., WIT Press, Compliant Structures in Nature and Engineering, chap. 8 C.H.M.Jenkins, ed., Billerica, MA, pp. 171-191 2005. [41] Shyy, W., and Liu, H., “Flapping Wings and Aerodynamic Lift: The Role of Leading-Edge Vortices,” Journal of Aircraft, Vol. 45, No. 12, 2007, pp. 2817– 2819. [42] Sankaranarayanan, S., Antony, R., and Suraj, C. S., “Technology Driven Program for the Development for the Development of a Fixed wing micro Air Vehicle at NALDRIVEN,” MAV Workshop - Challenges and Preferred Outcomes, National Aerospace Laboratories, Bangalore, Vol. 42, No. 4, 13 November 2008, pp. 865–873. [43] Kellogg, J., Bovais, C., Foch, R., McFarlane, H., Sullivan, C., Dahlburg, J., Gardner, J., Ramamurti, R., Gordon-Spears, D., Hartley, R., Kamgar-Parsi, B., Pipitone, F., Spears, W.and Sciambi, A., and Srull, D., “The NRL Micro Tactical Expendable (MITE) Air Vehicle,” Aeronautical Journal, Vol. 106, No. 1062, August 2002, pp. 431–441. 219 [44] Keennon, M., and Grasmeyer, J., “Development of Two MAVs and Vision of the Future of MAV Design,” 39th Aerospace Sciences Meeting and Exhibit, Reno, NV, 8-11 January 2001, AIAA 2001-127. [45] Keennon, M., and Grasmeyer, J., “Development of the Black Widow and Microbat MAVs and a Vision of the Future of MAV Design,” AIAA/ICAS International Air and Space Symposium and Exposition: The Next 100 Years, Dayton, OH, 14-17 July 2003, AIAA 2003-2901. [46] Chopra, I., “Hovering Micro Air Vehicles: Challenges and Opportunities,” AHS International Meeting on Advanced Rotorcraft Technology and Safety Operations (Heli Japan), Sonic City, Ohmiya, Japan, November 2010. [47] Hrishikeshavan, V., “Experimental Investigation of a Shrouded Rotor Mi- cro Air Vehicle in Hover and in Edgewise Gusts,” Doctoral dissertation, University of Maryland, College Park, Maryland, August 2011. [48] Barrett, R. M., “CFD-CSD Coupled Aeroelastic Analysis Of Flexible Flap- ping Wings For MAV Applications,” 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, CA, April 19-22 2004, AIAA 2004-1886. [49] Lakshminarayan, V. K., “Computational Investigation of Micro-Scale Coaxial Rotor Aerodynamics in Hover,” Ph.d. dissertation, University of Maryland at College Park, Department of Aerospace Engineering, College Park, Maryland, January 2009. [50] Fitchett, B. K., Development and Investigation of a Flapping Rotor for Micro Air Vehicles Doctoral Dissertation. University of Maryland, College Park, MD, 2007. [51] Ulrich, E., Pines, D., and Gerardi, S., “Autonomous Flight of a Samara MAV,” 65th Annual Forum of the American Helicopter Society, May 27–29, 2009. [52] Benedict, M., Ramasamy, M., Chopra, I., and Leishman, G., “Performance of a Cycloidal Rotor Concept for Micro Air Vehicle Applications,” Journal of the American Helicopter Society, Vol. 55, No. 022002, January 2010. [53] Sibilski, K., “Dynamics of Micor Air Vehicles with Flapping Wings,” Acta Polytechnica, Vol. 44, No. 2, 2004, pp. 15–21. [54] Videler, J. J., Avian Flight Oxford Ornithology Series. Oxford University Press, Oxford, 2005. [55] Biewener, A. A., “Muscle Function in Avian Flight: Achieving Power and Control,” Philosophical Transactions of the Royal Society, Vol. 366, 2011, pp. 1496–1506. 220 [56] Shyy, W., Klevebring, F., Nilsson, M., Sloan, J., Carroll, B., and Fuentes, C., “Rigid and Flexible Low Reynolds Number Airfoils,” Journal of Aircraft, Vol. 36, June 1999, pp. 523–529. [57] Lian, Y., and Shyy, W., “Aerodynamics of Low Reynolds Number Plunging Airfoil Under Gusty Environment,” 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 8-11 January 2007, Vol. AIAA 2007-71. [58] Vance, J. T., Faruque, I., and Humbert, J. S., “Kinematic Strategies for Mitigating Gust Perturbations in Insects,”. [59] Schenato, L., Deng, X., and Sastry, S., “Flight Control System for a Microme- chanical Flying Insect: Architecture and Implementation,” Proceedings of the IEEE International Conference on Robotics and Automation, Seoul, Korea, May 2001. [60] Sitti, M., “Piezoelectrically Actuated Four-Bar Mechanism with Two Flexible Links for Micromechanical Flying Insect Thorax,” IEEE/ASME Transactions on Mechatronics, Vol. 8, No. 1, March 2003, pp. 26–36. [61] Wood, R. J., “Liftoff of a 60mg Flapping-wing MAV,” Proceedings of the 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Diego, CA, Oct. 29 – Nov. 2 2007. [62] Avadhanula, S., Wood, R. J., Steltz, E., Yan, J., and Fearing, R. S., “Lift Force Improvements for the Micromechanical Flying Insect,” IEEE International Conference on Intelligent Robots and Systems, University of Bristol, Bristol, UK, October 27–31 2003, Vol. 2, pp. 1350–1356. [63] Steltz, E., Avadhanula, A., and Fearing, R. S., “High Lift Force with 275 Hz Wing Beat MFI,” Proceedings of the 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems, San Diego, CA, October 29 - November 2 2007, Vol. 2, pp. 3987–3992. [64] Perez-Arancibia, N., “First Controlled Vertical Flight of a Biologically Inspired Microrobot,” Bioinspiration and Biomimetics, Vol. 6, No. 3, 2011. [65] Samuel, P. D., and Humbert, J. S., “Design and Development of a Bio-Inspired Flapping Wing Micro Air Vehicle,” American Helicopter Society 68th Annual Forum Proceedings, Fort Worth, TX, 1–3 May 2012. [66] Mayo, D., Ramsey, J. P., and Leishman, J., “Particle Image Flow Diagnos- tics of a Flapping-Wing Micro Air Vehcile,” 27th Army Science Conference Proceedings, Orlando, FL, 7–10 December 2010. [67] Coleman, D., Benedict, M., Hrishikeshaven, V., and Chopra, I., “Design and Development of a Hover-Capable Two-Winged Flapping-Wing Micro- Air-Vehicle,” Proceeding of the AHS International Specialists’ Meeting on Unmanned Rotorcraft, Chandler, AZ, January 2015, Vol. 1. 221 [68] Keennon, M., Klingebiel, K., Won, W., and Andriukov, A., “Development of the Nano Hummingbird: A Tailless Flapping Wing Micro Air Vehicle,” 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Nashville, TN, 09 - 12 January 2012. [69] Wissa, A., “Analytical Modeling and Experimental Evaluation of a Passively Morphing Ornithopter Wing,” Doctoral dissertation, University of Maryland, College Park, Maryland, August 2012. [70] Koopmans, J., “Delfly Freeflight Autonomous flight of the Delfly in the wind tunnel using low-cost sensors,” MPhil, Delft University of Technology, June 15 2012. [71] Pornsin-sirirak, T. N., Lee, S. W., Nassef, H., Grasmeyer, J., Tai, Y. C., Ho, C. M., and Keennon, M., “Titanium-Alloy MEMS Wing Technology for a Battery-Powered Ornithopter,”. [72] Hein, B., and Chopra, I., “Hover Performance of Micro Air Vehicles: Rotors at Low Reynolds number,” Journal of the American Helicopter Society, Vol. 52, No. 3, July 2007, pp. 254–262. [73] Ramasamy, M., and Leishman, J. G., “Phase-Locked Particle Image Velocime- try Measurements of a Flapping Wing,” Journal of Aircraft, Vol. 43, No. 6, November-December 2006, pp. 1867–1875. [74] Bunget, G., and Seelecke, S., “BATMAV: A Biologically-Inspired Micro-Air Vehicle for Flapping Flight,” SPIE Conference Proceedings Active and Passive Smart Structures and Integrated Systems, San Diego, CA, 9 March 2008, Vol. 6928. [75] Spinella, I., Mammano, G. S., and Dragoni, E., “Conceptual Design and Simulation of a Compact Shape Memory Actuator for Rotary Motion,” Journal of Materials Engineering and Performance, Vol. 18, No. 5–6, 2009, pp. 638–648. [76] Liu, J., Yang, Z., Zhao, H., and Cheng, G., “Novel Precision Piezoelectric Step Rotary Actuator,” Frontiers of Mechanical Engineering, Vol. 2, No. 3, 2007, pp. 356–360. [77] Mueller, T. J., Kellogg, J., Ifju, P., and Shkarayev, S., “Introduction to the Design of Fixed-Wing Micro Air Vehicles Including Three Case Studies,” Piezo News: Quarterly News Letter of the Piezo Institute, Vol. 1, October 2007, pp. 1–7. [78] Clemons, L., Igarashi, H., and Hu, H., “An Experimental Study of Unsteady Vortex Structures in the Wake of a Piezoelectric Flapping Wing,” 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Reno, NV, 4–7 January 2010, Vol. AIAA 2010-1025. 222 [79] Meuller, D. C., “Energy Storage: Strong Momentum in High-End Batteries,” Credit Suisse, Suisse, 2007. [80] Mandal, S. K., Bhojwani, P. S., Mohanty, S. P., and Mahapatra, R. N., “IntellBatt: Towards Smarter Battery Design,” Proceedings of the 45th Design Automation Conference (DAC08), 2008, pp. 872–877. [81] Evans, A., A., B.-H., Galinski, H., Rupp, J. L. M., Ryll, T., B., S., Tolke, R., and Gauckler, L. J., “Micro-solid oxide fuel cells: status, challenges, and chances,” Monatsh Chemistry Journal, Vol. 140, February 2009, pp. 975–983. [82] Ren, J.-G., Wang, C., Wu, Q.-H., Liu, X., Yang, Y., Hea, L., and Zhang, W., “A silicon nanowirereduced graphene oxide composite as a high-performance lithium ion battery anode material,” Nanoscale, Vol. 6, Jauary 7 2013, pp. 3353– 3360. [83] Kumar, B., Kumara, J., Leesea, R., Fellner, J. P., Rodrigues, S. J., and Abrahamc, K. M., “A Solid-State, Rechargeable, Long Cycle Life LithiumAir Battery,” Journal of the Electrochemical Society, Vol. 157, No. 1, November 13 2009, pp. A50–A54. [84] Mehra, A., Zhang, X., Ayon, A. A., Waitz, I. A., Schmidt, M. A., and Spadaccini, C. M., “Six-wafer Combustion System for a Silicon Micro Gas Turbine Engine,” Journal of Microelectromechanical Systems, Vol. 9, No. 4, 2000, pp. 517–527. [85] Glassock, R. R., “Multimodal Hybrid Powerplant for Unmanned Aerial Systems (UAS) Robotics,” Proceedings of the 24th Bristol International Unmanned Air Vehicle Systems Conference, University of Bristol, Bristol, UK, January 2009. [86] Shan, X. C., Wang, Z. F., Maeda, R., F., S. Y., Wu, M., and Hua, J. S., “A Silicon-based Micro Gas Turbine Engine for Power Generation,” Symposium on Design, Test, Integration and Packaging of MEMS/MOEMS (DTIP 06), Stresa, Italy, April 2006, Vol. 2, pp. 3987–3992. [87] Michelson, R. C., “Novel Approaches to Miniature Flight Platforms,” Pro- ceedings of the Instn Mechanical Engineers: Journal of Aerospace Engineering, Vol. 218, September 2004, pp. 363–373. [88] Sane, S. P., and Dickinson, M. H., “The Aerodynamic Effects of Wing Rotation and a Revised Quasi-Steady Model of Flapping Flight,” The Journal of Experimental Biology, Vol. 205, No. 8, April 2002, pp. 1087–1096. [89] Birch, J. M., Dickson, W. B., and Dickinson, M. H., “Force Production and Flow Structure of the Leading Edge Vortex on Flapping Wings at High and Low Reynolds Numbers,” Journal of Experimental Biology, Vol. 207, No. 7, 2004, pp. 1063–1072. 223 [90] Birch, J. M., and Dickinson, M. H., “The Influence of Wing Wake Interactions on the Production of Aerodynamic Forces in Flapping Fight,” Journal of Experimental Biology, Vol. 206, Jul. 2003, pp. 2257–2272. [91] Chopra, I., “Hovering Micro Air Vehicles: Challenges and Opportunities,” Proceedings of the American Helicopter Society Specialists, Seoul, Korea, 15-17 October 2007, International Forum on Rotorcraft Multidisciplinary Technology. [92] Laitone, E. V., “Wind Tunnel Tests of Wings at Reynolds Numbers Below 70,000,” Experiments in Fluids, Vol. 23, No. 5, November 1997, pp. 45–49. [93] van den Berg, C., and Ellington, C. P., “The Vortex Wake of a “Hovering” Model Hawkmoth,” Philosophical Transactions of the Royal Society of London, Series B, Biologial Sciences, Vol. 352, No. 1351, 1997, pp. 317–328. [94] Willmott, A. P., Ellington, C. P., and Thomas, A. L. R., “Flow Visualization and Unsteady Aerodynamics in the Flight of the Hawkmoth, Manduca sexta,” Philosophical Transactions of the Royal Society of London, Series B, Biologial Sciences, Vol. 352, No. 1351, March 1997, pp. 303–316. [95] van den Berg, C., and Ellington, C. P., “The Three-Dimensional Leading-Edge Vortex of a “Hovering” Model Hawkmoth,” Philosophical Transactions of the Royal Society of London, Series B, Biologial Sciences, Vol. 352, No. 1351, 1997, pp. 329–340. [96] Ellington, C. P., and Usherwood, J. R., “Lift and Drag Characteristics of Rotary and Flapping Wings,” In Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, T. J. Mueller, Ed., Vol. 195 of Progress in Astronautics and Aeronautics. American Institute of Aeronautics and Astronautics, Inc., 1801 Alexander Bell Drive, Reston, Virginia 20191-4344, 2001, Ch. 12, pp. 231–248. [97] McCroskey, W. J., “Unsteady Airfoils,” Annual Review of Fluid Mechanics, Vol. 14, 1982, pp. 285–311. [98] Leishman, J., “Dynamic stall experiments on the NACA 23012 aerofoil,” Experiments in Fluids, Vol. 9, 1990, pp. 49–58. [99] Willmott, A. P., and Ellington, C. P., “The Mechanics of Flight in the Hawkmoth Manduca sexta, II. Aerodynamic Consequences of Kinematic and Morphological Variation,” Journal of Experimental Biology, Vol. 200, 1997, pp. 2723–2745. [100] Tarascio, M. J., Ramasamy, M., Chopra, I., and Leishman, J. G., “Flow Visu- alization of Micro Air Vehicle Scaled Insect-Based Flapping Wings,” Journal of Aircraft, Vol. 42, No. 2, March-April 2005, pp. 385–390. [101] Lentink, D., and Dickinson, M. H., “Rotational Accelerations Stabilize Leading Edge Vortices on Revolving Fly Wings,” Journal of Experimental Biology, Vol. 212, April 2009, pp. 2705–2719. 224 [102] Weis-Fogh, T., “Flapping Flight and Power in Birds and Insects,” Swimming and Flying in Nature, Pasadena, California, 8-12 July 1975, T. Y.-T. Wu, C. J. Brokaw, and C. Brennen, Eds., Vol. 2, Plenum Press, London, pp. 729–762. [103] Maxworthy, T., “The Fluid Dynamics of Insect Flight,” Annual Review of Fluid Mechanics, Vol. 13, 1981, pp. 329–350. [104] Ellington, C. P., “The Aerodynamics of Hovering Insect Flight, I. The Quasi- Steady Analysis,” Philosophical Transactions of the Royal Society of London, Series B, Biologial Sciences, Vol. 305, No. 1122, February 1984, pp. 41–78. [105] Ennos, A. R., “The Kinematics and Aerodynamics of the Free Flight of some Diptera,” Journal of Experimental Biology, Vol. 142, October 1989, pp. 49–85. [106] Lehmann, F.-O., Sane, S. P., and Dickinson, M., “The Aerodynamic Effects of WingWing Interaction in Flapping Insect Wings,” Journal of Experimental Biology, Vol. 208, No. 16, August 15 2005, pp. 3075–3092. [107] Agrawal, A., and Agrawal, S., “Design of bio-inspired flexible wings for flapping wing micro-sized air vehicle applications.,” Advanced Robotics, Vol. 23, 2009, pp. 979–1002. [108] Ansari, S. A., Z˙bikowski, R., and Knowles, K., “Aerodynamic modelling of insect-like flapping flight for micro air vehicles,” Progress in Aerospace Sciences, Vol. 42, No. 2, February 2006, pp. 129–172. [109] Sane, S. P., and Dickinson, M. H., “The Control of Flight Force By A Flapping Wing: Lift and Drag Production,” The Journal of Experimental Biology, Vol. 204, No. 15, August 2001, pp. 2607–2626. [110] Katz, J., and Plotkin, A., Low Speed Aerodynamics Cambridge University Press, 2001. [111] Carr, L. W., McAlister, K. W., and McCroskey, W. J., “Analysis of the Development of Dynamic Stall Based on Oscillating Airfoil Experiments,” NASA NASA-TN D-8382, Ames Research Center and U.S. Army Air Mobility Research and Development Laboratory, Moffett Field, CA, 1977. [112] Theodorsen, T., “General Theory of Aerodynamic Instability and the Mecha- nism of Flutter,” NACA 496, Washington, D.C., 1949. [113] Liiva, J., “Unsteady Aerodynamic and Stall Effects on Helicopter Rotor Blade Airfoil Sections,” Journal of Aircraft, Vol. 6, No. 1, 1969, pp. 46–51. [114] Knoller, R., “Die Gesetze des Luftwiderstandes,” Flugund Motortechnik (Wien), Vol. 3, No. 21, 1909, pp. 1–7. [115] Betz, A., “Ein Beitrag zur Erklaerung des Segelfluges,” Zeitschrift fur Flugtech- nik und Motorluftschiffahrt, Vol. 3, 1912, pp. 269–272. 225 [116] Blick, E. F., Watson, D., Belie, G., and Chu, H., “Bird Aerodynamic Exper- iments,” Swimming and Flying in Nature, Pasadena, California, 8-12 July 1975, T. Y.-T. Wu, C. J. Brokaw, and C. Brennen, Eds., Vol. 2, Plenum Press, London, pp. 939–952. [117] Tobalske, B. W., “Biomechanics of bird flight,” Journal of Experimental Biology, Vol. 210, 2007, pp. 3135–3146. [118] Bomphery, R., Lawson, N., Harding, N., Taylor, G., and Thomas, A., “The Aerodynamics of Manduca Sexta: Digital Particle Image Velocimetry analysis of the Leading Edge Vortex,” The Journal of Experimental Biology, Vol. 208, 2005, pp. 1079–1094. [119] Warrick, D. R., Tobalske, B. W., and Powers, D. R., “Lift Production in the Hovering Hummingbird,” Proceedings of The Royal Society, Vol. 276, 2009, pp. 3747–3752. [120] Wang, H., Zeng, L., Liu, H., and Yin, C., “Measuring Wing Kinematics, Flight Trajectory and Body Attitude During Forward Flight and Turning Maneuvers in Dragonflies,” The Journal of Experimental Biology, Vol. 206, September 2003, pp. 745–757. [121] Alexander, D. E., Natural Flyers: Birds, Insects and the Biomechanics of Flight The Johns Hopkins University Press, Baltimore, Maryland, 2002. [122] Ellington, C. P., “The Aerodynamics of Hovering Insect Flight, IV. Aerody- namic Mechanisms,” Philosophical Transactions of the Royal Society of London, Series B, Biologial Sciences, Vol. 305, No. 1122, February 1984, pp. 79–113. [123] DeVoria, A., Mahajan, P., and Ringuette, M., “Vortex Formation and Satu- ration for Low-Aspect-Ratio Rotating Flat Plates at Low Reynolds Number,” 49th AIAA Aerospace Sciences Meeting and Exhibit, Orlando, FL, 4-7 January 2011, AIAA 2011-396. [124] Kim, D., and Gharib, M., “Experimental study of three-dimensional vortex structures in transalating and rotating plates,” Experimental Fluids, Vol. 49, 2010, pp. 329–339. [125] Ozen, C. A., and Rockwell, D., “Flow Structure on a Rotating Plate,” Experi- ments in Fluids, 2012. [126] Granlund, K., Ol, M., Bernal, L., and Kast, S., “Experiments on Free-to-Pivot Hover Motions of Flat Plates,” 40th AIAA Fluid Dynamics Conference and Exhibit, Chicago, IL, 28 June - 1 July 2010, AIAA 2010-4456. [127] Pitt Ford, C. W., and Babinsky, H., “Lift and the Leading Edge Vortex,” Submitted to the Journal of Fluid Mechanics, January 2013, pp. 0–0. 226 [128] D., M., “Evolution and Breakdown of a Leading Edge Vortex on a Rotating Wing,” 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Grapevine, TX, 7-10 January 2013. [129] Seshadri, P., Biometric Insect Flapping Wing Aerodynamics and Controls for Micro Air Vehicles Master of Science Thesis. University of Maryland, College Park, MD, 2011. [130] Hart, A., and Ukeiley, L., “Low Reynolds Number Unsteady Aerodynamics over a Pitching-Plunging Flat Plate,” 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Orlando, Florida, 4-7 January 2010, AIAA 2010387,. [131] Yuan, W., Lee, R., Hoogkamp, E., and Khalid, M., “Numerical and Experi- mental Simulations of Flapping Wings,” International Journal of Micro Air Vehicles, Vol. 2, No. 3, September 2010, pp. 181–209. [132] Mayo, D. B., Lankford, J. L., Benedict, M., and Chopra, I., “Experimental and Computational Aerodynamic Investigation of Avian-Based Rigid Flapping Wings for MAV Applications,” Fifth Decennial AHS Aeromechanics Specialists’ Conference, San Francisco, CA, January 2014. [133] Benedict, M., Coleman, D., Mayo, D. B., and Chopra, I., “Force and Flowfield Measurements on a Rigid Wing Undergoing Hover-Capable Flapping and Pitching Kinematics in Air at MAV-Scale Reynolds Numbers,” AIAA 54th ASC Structures, Structural Dynamics, and Materials Conference, Boston, MA, 8-11 April 2013. [134] Thiria, B., G.-D. R., “How wing compliance drives the efficiency of self- propelled flapping flyers,” Physics Review, Vol. E 82:015303(R), July 2010. [135] Frampton, K., Goldfarb, M., Monopoli, D., and Cveticanin, D., “Passive Aeroelastic Tailoring for Optimal Flapping Wings,” In Mueller’s Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, Progress in Astro- nautics and Aeronautics, AIAA, Vol. 195, 2001, pp. 473–482. [136] Combes, S. A., and Daniel, T. L., “Into Thin Air: Contributions of Aerody- namic and Inertial-Elastic Forces to Wing Bending in the Hawkmoth Manduca sexta,” Journal of Experimental Biology, Vol. 206, No. 1, 2003, pp. 2999–3006. [137] Tanaka, H., and Wood, R., “Effect of Flexural and Torsional Wing Flexibility on Lift Generation in Hoverfly Flight,” Journal of Integrative and Comparative Biology, Vol. 1, May 2011, pp. 1–9. [138] Ho, S., Nassefa, H., Pornsinsirirakb, N., Taib, Y., and Ho, C., “Unsteady Aerodynamics and Flow Control for Flapping Wing Flyers,” Progress in Aerospace Sciences, Vol. 39, January 2003, pp. 635–681. 227 [139] Heathcote, S., Martin, D., and Gursul, I., “Flexible flapping airfoil propulsion at zero freestream velocity,” AIAA, Vol. 42, No. 11, November 2004, pp. 2192– 2204. [140] Heathcote, S., Wang, Z., and Gursul, I., “Effect of Spanwise Flexibility on Flapping Wing Propulsion,” Journal of Fluids and Structures, Vol. 24, No. 2, February 2011, pp. 183–199. [141] Malhan, R. P., Baeder, J., Chopra, I., and Masarati, P., “CFD-CSD Coupled Aeroelastic Analysis Of Flexible Flapping Wings For MAV Applications,” 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Boston, MA, 8-11 April 2013, AIAA 2013-1644. [142] Malhan, R. P., Benedict, M., and Chopra, I., “Experimental Studies to Understand the Hover and Forward Flight Performance of a MAV-Scale Flap- ping Wing Concept,” Journal of the American Helicopter Society, Vol. 57, No. 022003, February 2012, pp. 1–11. [143] Minotti, F. O., “Unsteady Two-Dimensional Theory of a Flapping Wing,” Physical Review E, Statistical Physics, Plasma, Fluids, and Related Interdisci- plinary Topics, Vol. 66, No. 051907, January 2002, pp. 1–10. [144] Bush, B., and Baeder, J., “Force Production Mechanisms of Flapping-Wing MAV,” American Helicopter Society 64th Annual Forum Proceedings, Montreal, Canada, 29 April-1 May 2008. [145] Malhan, R., Benedict, M., and Chopra, I., “Experimental Investigation of a Flapping Wing Concept in Hover and Forward Flight for Micro Air Vehicle Applications,” Proceedings of the 66th Annual American Helicopter Society Forum, Phoenix, Arizona, May 2010, Vol. 2, pp. 939–952. [146] Bush, B., and Baeder, J., “Computational Investigation of Flapping-Wing Flight,” 37th AIAA Fluid Dynamics Conference and Exhibit, Miami, FL, 25-28 June 2007, no. 4209 in 1. [147] Ramamurti, R., and Sandberg, W., “Simulation of Flow About Flapping Airfoils Using Finite Element Incompressible Flow Solver,” AIAA Journal, Vol. 39, No. 2, February 2001, pp. 253–260. [148] Sun, M., and Tang, J., “Unsteady Aerodynamic Force Generation by a Model Fruit Fly Wing in Flapping Motion,” Journal of Experimental Biology, Vol. 205, 2002, pp. 55–70. [149] Sun, M., and Tang, J., “Unsteady Aerodynamic Force Generation by a Model Fruit Fly Wing in Flapping Motion,” The Journal of Experimental Biology, Vol. 205, 2002, pp. 55–70. [150] DeLaurier, J. D., An Aerodynamic Model for Flapping-Wing Flight, 5th ed., No. (97)964 in pp. 125–130. Aeronaut. J., 1993. 228 [151] Singh, B., “Dyamics and Aeroelastics of Hover Capable Flapping Wings: Experiments and Analysis,” Doctoral dissertation, University of Maryland College Park, May 2006. [152] Polhamus, E. C., “A Concept of the Vortex Lift of Sharp-Edge Delta Wings Based on A Leading-Edge-Suction Analogy,” NASA TN-D-3767, Langley Reseach Center, December 1966. [153] Liu, H., Ellington, C., Kawachi, K., van den Berg, C., and Willmott, A., “A computational fluid dynamic study of hawkmoth hovering,” The Journal of Experimental Biology, Vol. 201, 1998, pp. 461–477. [154] Wang, Z., “Two-Dimensional Mechanism for Insect Hovering,” Physical Review Letters, Vol. 85, No. 10, 2000, pp. 2216–2219. [155] Newman, B. G., “Aerodynamic theory for membranes and sails,” Progress in Aerospace Sciences, Vol. 24, 1987, pp. 1–27. [156] Thwaites, B., “The Aerodynamic Theory of Sails I: Two-Dimensional Sails,” Proceedings of the Royal Society of London. Series A, Mathematical and Phys- ical Sciences, Vol. 261, No. 1306, May 1961, pp. 402–422. [157] Smith, M. J., Komerath, N., Ames, R., Wong, O., and Pearson, J., “Per- formance Analysis of a Wing With Multiple Winglets,” 19th AIAA Applied Aerodynamics Conference, Anaheim, CA, 11-14 June 2001, no. AIAA in 2001- 2407. [158] Kima, D, K., Lee, J. S., Lee, J. Y., and H., H. J., “An Aeroelastic Analysis of a Flexible Flapping Wing Using Modified Strip Theory,” SPIE 15th Annual Symposium Smart Structures and Materials, 2008, Vol. 6928. [159] Malhan, R., Benedict, M., and Chopra, I., “CFD-CSD Coupled Aeroelastic Analysis of Flexible Flapping Wings for MAV Applications: Methodology Validation,” 53rd AIAA Structural Dynamics, and Materials Conference, Honolulu, HI, 2326 April 2007, AIAA 2012-1636. [160] Zhang, X. W., and Zhou, C. Y., “Computational study on the hovering mechanisms of a chordwise flexible wing,” Journal of Aerospace Engineering, Vol. 226, April 2011, pp. 1–14. [161] Pulliam, T., and Chaussee, D., “A Diagonal Form of an Implicit Approximate Factorization Algorithm,” pp. 347–363. [162] Van Leer, B., “Towards the Ultimate Conservative Diff Scheme V. A Second- Order Sequel To Godunovs Method,” pp. 229–248. [163] Roe, P., “Approximate Riemann Solvers, Parameter Vectors and Difference Schemes,” pp. 250–258. 229 [164] Koren, B., “Multigrid and Defect Correction for the Steady Navier-Stokes Equations,” Proceedings of the 11th International Conference on Numerical Methods in Fluid Dynamics, Willamsburg, VA, June 1988, Vol. 13th AIAA Computational Fluid Dynamics Conference. [165] Buelow, P. E. O., Schwer, D. A., Feng, J., and Merkle, C. L., “A Precon- ditioned Dual-Time, Diagonalized ADI scheme for Unsteady Computations,” AIAA paper 19972101, Snowmass Village, CO, July 1997, Vol. 13th AIAA Computational Fluid Dynamics Conference. [166] Pandya, S. A., Venkateswaran, S., and Pulliam, T. H., “Implementation of Pre- conditioned Dual-Time Procedures in OVER-FLOW,” 41st AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, January 2003, Vol. AIAA paper 20030072. [167] Spalart, P. R., and Allmaras, S. R., “A One-equation Turbulence Model for Aerodynamic Flows,” 30th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 6–9 January 1992, Vol. AIAA Paper 1992-0439. [168] Lee, Y., “On Overset Grids Connectivity and Automated Vortex Tracking in Rotorcraft CFD,” Ph.d. dissertation, University of Maryland at College Park, Aerospace Engineering Department, College Park, Maryland, January 2008. [169] Breuer, M., Jovii, N., and Mazaev, K., “Comparison of DES, RANS and LES for the separated flow around a flat plate at high incidence,” International Journal for Numerical Methods in Fluids, Vol. 41, No. 1, February 2003, pp. 357–388. [170] Piomelli, U., “Large-eddy simulation: achievements and challenges,” Progress in Aerospace Sciences, Vol. 35, No. 1, January 1999, pp. 355–362. [171] Gordnier, R., Chimakurthi, S., Cesnik, C., and Attar, P., “Implicit LES Simulations of a Flexible Flapping Wing,” 51st AIAA Structural Dynamics, and Materials Conference, Orlando, FL, 12–15 April 2012, no. 2012-1636 in 1. [172] Masarati, P., Morandini, M., Quaranta, G., and Vescovini, R., Multibody Analysis of a Micro-Aerial Vehicle Flapping Wing Multibody Dynamics 2011, Brussels, Belgium, July 2011. [173] Masarati, P., and Sitaraman, J., “Tightly Coupled CFD/Multibody Analysis of NREL Unsteady Aerodynamic Experiment Phase VI Rotor,” 49th AIAA Aerospace Sciences Meeting, Orlando, Florida, January 2011. [174] Witkowski, W., “4-node combined shell element with semi-EAS-ANS strain interpolations in 6-parameter shell theories with drilling degrees of freedom,” Journal of Computational Mechanics, Vol. 43, No. 2, January 2009, pp. 307–319. [175] Adrian, R. J., “Particle-Imaging Techniques for Experimental Fluid Mechanics,” Annual Reviews of Fluid Mechanics, Vol. 23, January 1991, pp. 261–304. 230 [176] Raffel, M., Willert, C. E., and Kompenhans, J., Particle Image Velocimetry: A Practical Guide Springer, 1998. [177] Westerweel, J., Dabiri, D., and Gharib, M., “The Effect of a Discrete Window Offset on the Accuracy of Cross-Correlation Analysis of Digital PIV Recordings,” Experiments in Fluids, Vol. 23, 1997, pp. 20–28. [178] Tropea, C., and L., Y. A., Springer Handbook of Experimental Fluid Mechanics Springer-Verlag, Berlin, 2007. [179] Scarano, F., “Iterative Image Deformation Methods in PIV,” Measurement Science and Technology, Vol. 13, January 2002, pp. R1–R19. [180] Spedding, G. R., and Hedenstrom, A., “PIV-based investigations of animal flight,” Experiments in Fluids, Vol. 46, December 2009, pp. 749–763. [181] Taylor, G. K., Nudds, R. L., and Thomas, A. L. R., “Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency,” Letters to Nature, Vol. 425, October 16, 2003, pp. 707–711. [182] Huang, H., Dabiri, D., and Gharib, M., “On errors of digital particle image velocimetry,” Journal of Measurement Science Technology, Vol. 8, January 1997, pp. 1427–1440. [183] Amabili, M., Nonlinear Vibrations and Stability of Shells and Plates Cambridge University Press, January 2008. [184] Graftieaux, L., Michard, M., and Grosjean, N., “Combining PIV, POD and Vortex Identification Algorithms for the Study of Unsteady Turbulent Swirling Flows,” Measurement Science and Technology, Vol. 12, No. 9, August 2001, pp. 1422–1429. [185] DaSilva, C. B., and Pereira, J. C., “Invariants of the velocity-gradient, rate-of- strain, and rate-of-rotation tensors across the turbulent/nonturbulent interface in jets,” Physics of Fluids, Vol. 20, 2008. [186] Usherwood, J. R., and Ellington, C. P., “The Aerodynamics of Revolving Wings: I. Model Hawkmoth Wings,” Journal of Experimental Biology, Vol. 205, No. 11, 2002, pp. 1547–1564. [187] Usherwood, J. R., and Ellington, C. P., “The Aerodynamics of Revolving Wings: II. Propeller Force Coefficients from Mayfly to Quail,” Journal of Experimental Biology, Vol. 205, No. 11, 2002, pp. 1565–1576. 231