Aeromechanics of Highly Compliant Structures: Bat wings, compliant membranes and flexibly mounted flat plates by Arnold J. Song B.A., English, Dartmouth College, 1999 B.S., Aerospace Engineering, Case Western Reserve University, 2005 Sc.M., Engineering, Brown University, 2011 Thesis Submitted in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in the School of Engineering at Brown University May 2013 c Copyright 2013 by Arnold J. Song Abstract of “Aeromechanics of Highly Compliant Structures: Bat wings, compliant mem- branes and flexibly mounted flat plates,” by Arnold J. Song, Ph.D., Brown University, May 2013 We present a study of the aeromechanics of highly compliant structures in the context of animal flight, with special attention paid to bat flight. Bats are unique among animal flyers because of their highly articulated wings that are composed of a thin skin membrane and a skeletal structure that contains many long and slender bones. Therefore, significant shape change occurs throughout the wingbeat cycle as a consequence of a strong coupling between the aerodynamic forces and the wing structure. Using high speed videography and photogrammetry techniques, the in-flight shape, motion and articulation of the wing were measured for wind tunnel flights of several individuals of Cynopterus brachyotis for wind speeds ranging from U = 2.6 − 6.6 m/s. The inboard portion of the wing membrane exhibited large, anisotropic strains with the membrane area increasing to nearly three times the minimum area over a wingbeat cycle, but with little to no flight speed dependence. In addition, the leading edge was found to nearly align with the oncoming flow for flight conditions tested perhaps limiting or controlling flow separation that would lead to the formation of a leading edge vortex (LEV). The aeromechanics of the bat wing is distilled into a system consisting of two compo- nents: 1) a compliant membrane and 2) an elastic structure with a resonant frequency near the natural vortex shedding frequency. The steady aerodynamic and unsteady aerome- chanical behavior of a latex membrane wing were measured in a series of wind tunnels tests, which showed the enhanced lift behavior due the adaptive cambering of these elastic membrane wings. We modulate the strength and stability of the flat plate leading edge vortex using only one degree of freedom (pitching) in both passive and forced plate motion. We measure the resultant aerodynamic moment as an indicator of the vortex strength and near-plate residence time. We conclude with a description and implementation of a cyberphysical flat plate that is mounted to a virtual spring-damper system that enables software control of the torsion spring stiffness and damping. This dissertation by Arnold J. Song is accepted in its present form by the School of Engineering as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Kenneth S. Breuer, Advisor Recommended to the Graduate Council Date Sharon M. Swartz, Reader Date Shreyas Mandre, Reader Approved by the Graduate Council Date Peter Weber Dean of Graduate School iii The Vita of Arnold J. Song Arnold J. Song was born on April 25, 1977 in Columbus, Ohio. He resided in San Jose, California, until graduating from Bellarmine College Preparatory in 1995. Upon graduation from high school, he attended Dartmouth College in Hanover, New Hampshire, and received a Bachelor of Arts degree in English in 1999. He attended Case Western Reserve University in Cleveland, Ohio, and received a Bachelor of Science degree in Aerospace Engineering in 2005. Following his undergraduate work at Case Western, Mr. Song began pursuing graduate studies at Brown University in Providence, Rhode Island, where he worked toward a Doctor of Philosophy in Engineering investigating the aeromechanics of highly compliant structures. During his studies at Brown, he was awarded the Simon Ostrach Fellowship in 2008. In 2010, he obtained a position as Research Mechanical Engineer at the Cold Regions Research and Engineering Laboratory (CRREL) in Hanover, NH, pursuing research in flow-driven transport of granular media, discrete element modeling of soil dynamics and metamorphism of snow. Upon the completion of his graduate studies, he will continue with his work at CRREL. iv Preface and acknowledgments This thesis is the culmination of many years of effort that would not have been possible without the help and support of a great many people to which I am indebted. First and foremost, I would like to express my deep gratitude for the guidance and support of my advisor, Professor Kenny Breuer. He allowed me to be adventurous in my scientific pursuits and provided me with opportunities and experiences that will always be a part of my professional endeavors. But, I will be most grateful for the life lessons that he taught me and his encouraging words that helped me stay the course during a trying period. Without a doubt, this work would not have been possible without his patience and compassion. I would like to thank my reviewers, Professors Sharon Swartz and Shreyas Mandre, for the insight and suggestions that have made this work better. I would like to especially thank Professor Swartz for opening my eyes to the wonderful and intriguing world of bats and biological flyers. Rarely is science pursued in a vacuum; I must thank the many undergraduates, graduate students and postdoctoral fellows that kept my research machinery humming along with discussions, late night food runs, a helping hand or the offering of terrible puns: Professor Jeff Guasto, Rye Waldman, Dr. Tatjana Hubel, Dr. Charles Peguero, Emily Israeli, Ricky Galvao, Max Tuttman, Benjamin Strom, Kyohei Onoue, Dr. Jen Franck, Qian Bian, Shawn Kitchner, Matt Novick, Professor Peter Huang, Dr. Dan Riskin, Dr. Pepe Iriate-Diaz and Professor Dave Willis. As this chapter in my life comes to a close, I am awestruck by the level of sacrifice borne by my mother and father in unwavering support of my education and pursuit of happiness. v Words cannot pay tribute to their love and faith as I wander and bump my way through life. Some of my fondest childhood memories are of sitting in the car with my father just outside the fence of the local municipal airport where we would listen to the chatter on the ATC scanner as we watched the Cessnas, Beechcraft and Pipers take off and land as the lazy summer sun slowly set. I cannot help but think that this work is a testament to my father’s passion for flight and the hours that we spent dreaming of riding through the clouds. My partner, my joy and my rock through this journey is Elizabeth Carpenter-Song. I am forever thankful for your warm and gentle support that you give to me each and every day. And last, but certainly not least, to my daughter Madeline, thank you for sleeping through the night to give your dad a little time to write. vi Dedication To my father, who never gave up on me. “If you put your mind to it, you can accomplish anything.” - George McFly, Back to the Future vii Contents Preface and acknowledgments v Dedication vii 1 Introduction 1 2 In-vivo measurement of the extreme shape morphing of bat wings 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.1 Compliant membrane wings . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.2 Leading edge alignment . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Steady forward flight measurements . . . . . . . . . . . . . . . . . . 11 2.2.2 Wing kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.3 Skin stretching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.4 Camber and angle calculations . . . . . . . . . . . . . . . . . . . . . 15 2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Skin stretching behavior . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Directionally preferential strain . . . . . . . . . . . . . . . . . . . . . 20 2.3.3 Leading edge alignment . . . . . . . . . . . . . . . . . . . . . . . . . 21 viii 2.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3 Aeromechanics of membrane wings with implications for animal flight 25 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Analytical Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Experimental Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.1 Compliant Membrane Models . . . . . . . . . . . . . . . . . . . . . . 33 3.3.2 Stereo Photogrammetry . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.3 Lift and Drag Measurements . . . . . . . . . . . . . . . . . . . . . . 35 3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4.1 Static membrane deflection . . . . . . . . . . . . . . . . . . . . . . . 38 3.4.2 Lift and Drag performance . . . . . . . . . . . . . . . . . . . . . . . 43 3.4.3 Dynamic membrane motions . . . . . . . . . . . . . . . . . . . . . . 49 3.5 Concluding remarks and implications for mammalian flight . . . . . . . . . 55 4 Experimental study of vortex-induced forces on a pitching flat plate 59 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1.1 Leading and trailing edge vortex formation and shedding . . . . . . 61 4.2 Materials and methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2.1 Elastically mounted flat plate . . . . . . . . . . . . . . . . . . . . . . 65 4.2.2 Servo-controlled flat plate . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3.1 Passive pitching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3.2 Wake velocity measurements . . . . . . . . . . . . . . . . . . . . . . 72 4.3.3 Forced oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 ix 4.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5 Vortex-induced vibrations of a flat plate mounted to a virtual spring and damper 90 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.2 Cyberphysical system design and setup . . . . . . . . . . . . . . . . . . . . . 94 5.2.1 Control scheme description . . . . . . . . . . . . . . . . . . . . . . . 96 5.2.2 Damping model identification . . . . . . . . . . . . . . . . . . . . . . 98 5.3 Cyberphysical system implementation . . . . . . . . . . . . . . . . . . . . . 100 5.3.1 Free oscillation test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.3.2 Vortex-induced vibrations – Preliminary results . . . . . . . . . . . . 103 5.4 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.4.1 Virtual moment of inertia . . . . . . . . . . . . . . . . . . . . . . . . 106 5.5 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 x List of Tables 5.1 Linear damping model fit results . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Quadratic damping model fit results . . . . . . . . . . . . . . . . . . . . . . 99 5.3 Free oscillation parameters – Varying stiffness . . . . . . . . . . . . . . . . . 103 5.4 Free oscillation parameters – Varying damping . . . . . . . . . . . . . . . . 104 xi List of Figures 2.1 Analysis mesh construction from markers . . . . . . . . . . . . . . . . . . . 12 2.2 Wing surface area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Variation of wing joint angles with respect to time . . . . . . . . . . . . . . 13 2.4 Camber extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.5 Plagiopatagium area variation over normalized wingbeat for one individual 18 2.6 Anisotropic strain in the plagiopatagium for Cynopterus brachyotis . . . . . 20 2.7 Camber extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.8 Angles with respect to the relative velocity at the wing . . . . . . . . . . . . 22 3.1 Examples of mammalian flight . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Model for compliant membrane wing subject to aerodynamic loading . . . . 29 3.3 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.4 Static deformation due to aerodynamic loading . . . . . . . . . . . . . . . . 39 3.5 Maximum camber as function of aerodynamic loading . . . . . . . . . . . . 39 3.6 Comparison of Tip and Root Deflection . . . . . . . . . . . . . . . . . . . . 42 3.7 Aerodynamic Performance of Rigid versus Compliant Wings . . . . . . . . . 43 3.8 Hysteresis in lift and drag due to compliance . . . . . . . . . . . . . . . . . 47 3.9 Effect of Aspect Ratio on Aerodynamic Performance . . . . . . . . . . . . . 48 xii 3.10 Spatial pattern of the unsteady vibration of the compliant wing (α = 28◦ , Re = 1.2 × 105 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Phase map of membrane drumhead modes (eigenmodes) . . . . . . . . . . . 50 3.12 Unsteady membrane vibrations . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.1 Resonant frequency dependence on rod length . . . . . . . . . . . . . . . . . 67 4.2 Passive flapper model diagram . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.3 Realtime video based motion tracking at 200 Hz . . . . . . . . . . . . . . . 68 4.4 Servo controlled pitching model . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.5 Wake velocity spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.6 Effect of aspect ratio on pitching intensity . . . . . . . . . . . . . . . . . . . 75 4.7 Net aerodynamic work per cycle . . . . . . . . . . . . . . . . . . . . . . . . 78 4.8 Torque versus angular position hysteresis . . . . . . . . . . . . . . . . . . . 82 4.9 Re = 6.7 × 104 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.10 Re = 8.0 × 104 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.11 Re = 9.3 × 104 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.12 Aerodynamic work, Ur = 10.0 . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.13 Transition to supercritical behavior for higher Reynolds number flows, α0 = 14◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.14 Explanatory model for positive net aerodynamic work in separated regime. 88 5.1 Description of a flat plate mounted to a virtual spring-damper system . . . 94 5.2 Block diagram of virtual spring-damper feedback loop . . . . . . . . . . . . 97 5.3 Determination of viscous damping . . . . . . . . . . . . . . . . . . . . . . . 98 5.4 Free oscillation tests – Varying stiffness . . . . . . . . . . . . . . . . . . . . 102 xiii 5.5 Free oscillation tests – Varying damping . . . . . . . . . . . . . . . . . . . . 104 5.6 Vortex-induced vibrations of flat plate for Re = 1 × 105 for several values of the virtual spring stiffness, k = 0.42, 0.83, 1.66 (k0 , 2k0 , 4k0 ) . . . . . . . . . 105 xiv Chapter 1 Introduction Micro-air vehicles (MAVs) are quickly becoming a widely used platform for surveillance, inspection and remote sensing applications. Their small size, with a wingspan dimensions of small bird down to large insects and relatively low cost compared to conventional air- craft allow for the deployment in areas that may either be hazardous, with undesirable conditions (e.g., extreme temperatures, lack of oxygen, etc.) or inaccessible to personnel or conventional aircraft. For example, a fleet of MAVs could be used to inspect sections of the 800 km long trans-Alaska pipeline (which transfers crude oil from the Prudhoe Bay to the southern port of Valdez, Alaska) with more detail and efficiency than with conventional air- or land-borne vehicles. These vehicles can navigate roadless terrain and potentially fly or hover in a confined airspace underneath the pipeline while replenishing power periodically at pump stations to increase deployment intervals. This would reduce the need to send personnel out to inspect the pipeline in severe, life-threatening conditions (average January temperature at Prudhoe Bay is −30◦ C. Being air-borne and having a small form factor, the potential applications for MAVs are numerous, however, maneuverability is a key factor in what will distinguish these vehicles from other unmanned systems. 1 2 MAV design falls in three primary categories: 1) fixed wing vehicles, 2) rotorcraft and 3) flappers. Fixed wing flyers are essentially scaled down versions of conventional aircraft with special attention to airfoil and wing design for the small scale (. 15 cm) and low speed (. 20 m/s) of these vehicles; the Reynolds number for these flyers is . 104 . Like conventional aircraft, these vehicles are well-suited for fast, forward flight along with a capacity to glide. The lift, however, is dependent on the forward velocity of these vehicles making a hover flight mode difficult to achieve because the lift would depend solely on a redirection of the propulsive force with the wings contributing nothing to weight support. Lift generated through thrust vectoring is the flight approach taken with rotorcraft where a set of typically high aspect ratio wings (blades) are rotated and oriented for directional flight. This arrangement is very effective for hovering in quiescent conditions, however, in a directional flight mode, the rotorcraft require complex control schemes or multi-rotor configurations to account for asymmetry in lift production as a result of the retreating and advancing phases of a blade’s rotation cycle. In some cases, the retreating blade can experience stall due to the reduced relative blade velocity in forward flight [1], which is a phenomenon that will be more problematic at the low Reynolds number of MAVs. Flapping flight combines the forward flight performance of fixed aircraft with the hovering capability of rotorcraft using the same wing architecture for both modes of flight. In addition, a flapper could be configured to perform in a glide mode periodically for increased operational efficiency. We turn to the flapping flight of biological flyers for inspiration and insight into flight mechanisms that might optimize low-Reynolds number flight. Bats are remarkable flyers; they capture insects on the wing, perform 180 degree turns in the span of just a couple of wingbeats, and sustain hovering flight while feeding on nectar. There are many species of bat that are small and highly maneuverable in confined spaces 3 and can easily transition from fast forward to hovering modes of flight many of the desirable attributes for MAV design. We attribute the flight capabilities of bats to their unique wing structure that enables dramatic wing shape morphing; the wings of bats being composed of an elastic skin membrane [2] that is manipulated through articulation of a compliant skeletal structure [3]. Due the complexity in the flight of bats, it is difficult to parse out the contributions of the membrane and skeletal compliance on their aerodynamic performance. In addition, the challenges in studying bat flight is the inherently unsteady nature of flapping flight and the intermediate Reynolds number regime (Re ≈ O(104 )) where many complex aerodynamic phenomena, such as transition to turbulence and laminar separation, are both present and hard to predict [4, 5, 6]. To address these issues, we have thus embarked on a systematic study with an attempt to isolate different morphological features present in mammalian flight to examine the the role that each plays. In this dissertation study, we have decomposed the salient components of the bat wing into two distinct types of wings: 1) a compliant (extensible) membrane wing and 2) a thin, flat plate that is subject to passive (flow-induced) or forced (flow-independent) motion. The key anatomical feature of bats that distinguishes these flyers from birds and in- sects is a set of wings that is primarily composed of a compliant skin membrane with low extensional stiffness and negligible bending stiffness [2]. These membrane wings composed of thin, yet tough skin are observed not only in bats but also mammalian gliders such as flying squirrels and sugar gliders [7, 8]. The extraordinary flight agility and maneuverabil- ity of these mammalian flyers and gliders can be attributed in part to the high degree of articulation in the wing (from the shoulder and elbow, but particularly the wrist and finger which shape the outboard portion of the wing)[9]. Also, the relatively low modulus of bat 4 wings results in strong coupling between the wing shape and aerodynamic forces where small changes in wing geometry can have significant effects on aerodynamic performance. For example, the aerodynamic forces acting on a compliant wing will cause the wing to push out and increase in camber, which typically lead to increases in the lift and drag on the wing and can also increase the robustness of the wing’s lift characteristics at angles of attack near stall [10]. When flow separation occurs on the membrane wing, the decrease in aerodynamic forces also decreases the airfoil camber potentially avoiding a catastrophic drop in lift due to stall and giving the animal time to adjust wing mechanics for stall conditions. We approach the issue of skeletal compliance as a mechanism for the formation and con- trol of vortex structures that can enhance lift production in flapping flight. These vortical structures that form in the near wake of the leading edge are referred to as leading edge vortices (LEVs) and can increase lift values over those observed in conventional aerodynam- ics due to the increase in suction peak magnitude [11, 12, 13, 14, 15]. Several studies have examined the time scale of LEV formation and stabilization due to translation, pitching or rotational motion or a combination of these degrees of freedom [16, 17, 18, 19]. Previous studies on LEV stabilization mechanisms have focused primarily on the advection of circu- lation out of the LEV to balance the circulation that is being added from the leading edge separation, which is a three-dimensional mechanism that relies on centripetal acceleration of the fluid. We propose and examine a LEV control mechanism that is two-dimensional and controls the unsteady generation of circulation at the leading edge with only a pitching motion that varies the angle of attack and the relative flow velocity at the leading edge of the plate. The unique structure of the mammalian flyers wings suggests that further study of the aerodynamics and biomechanical mechanisms in these animals may provide useful infor- 5 mation and inspiration for the design of MAVs that utilize active and passive morphing of structurally compliant wings. Dissertation outline This dissertation contains experimental studies of the aeromechanics and vortex-induced vibrations of highly compliant structures with an emphasis on the implications to animal flight. The content and relevant contributions are as follows: Chapter 1. Introduction. A description of aeromechanics and vortex-induced vibra- tions in relation to animal flight and micro-air vehicle design. Chapter 2. “In-vivo measurement of the extreme shape morphing of bat wings” by Arnold J. Song, Sharon M. Swartz, Daniel K. Riskin, Jose Iriate-Diaz and Kenneth S. Breuer. The experiments were conceived by Jose Iriate-Diaz, Swartz, Daniel K. Riskin and Breuer. The wind tunnel measurements were conducted and processed by Iriate-Diaz and Riskin with assistance from postdoctoral, graduate and undergraduate members of the Swartz and Breuer laboratories at Brown University. This study was part of an effort that was the first of its kind to capture a time-resolved three-dimensional reconstruction of a bat’s wing shape in vivo and in flight. Song developed the wing morphing analysis that examined the possible mechanisms for the large area changes observed during the wingbeat cycle of the Cynopterus brachyotis. This analyisis also demonstrates leading edge alignment of the wing that may prevent the formation of a LEV. Chapter 3. “Aeromechanics of membrane wings with implications for animal flight” by Arnold J. Song, Xiaodong Tian, Emily Israeli, Ricardo Galvao, Kristin L. Bishop, Sharon M. Swartz and Kenneth S. Breuer. AIAA Journal, Vol. 46, No. 8, pp. 2096-2106, 2008. Experiments were conceived by Bishop, Tian, Israeli, Galvao, Swartz and Breuer. The lift 6 experiments were conducted and analyzed by Israeli and Tian. The membrane shape mea- surements and analysis were conducted by Galvao and Song. The analytical model for the membrane deformation in response to aerodynamic forcing was developed by Song. The results were compiled and reported by Song and Breuer. This study was first to charac- terize the aerodynamic force generation and adaptive cambering behavior of a compliant, extensible membrane wing. As of April 2013, this paper has been cited 59 times. Chapter 4. “Experimental study of vortex-induced forces on a pitching flat plate” by Arnold J. Song and Kenneth S. Breuer. The experiments were conceived by Song and Breuer. The experimental design, measurements and analysis of the wake velocity and plate motion were conducted by Song. The results were interpreted by Song and Breuer. This study examines how the growth and stability of a leading edge vortex might be controlled through modulation of the circulation that feeds into the near wake of the leading edge through pitching of flat plate. Chapter 5. “Vortex-induced vibrations of a flat plate mounted to a virtual spring and damper” by Arnold J. Song, Benjamin W. Strom, Kyohei Onoue and Kenneth S. Breuer. Song conceived the experiments, cyberphysical apparatus and control design with assistance from Strom in the apparatus fabrication. The plant characterization and free oscillation measurements were conducted by Song. The vortex-induced vibration demonstration was conducted by Onoue. The cyberphysical system was implemented and characterized by Song. The design of this virtual spring and damper system builds upon work that was previously done by Hover and coworkers in their development of their virtual cable testing apparatus [20]. That effort along with the recent work of Mackowski and Williamson [21] and Lee and Bernitsas [22] are for the VIV of flexibly mounted cylinders in relatively slow moving water. Our cyberphysical system was designed for the higher frequencies associated 7 with VIV dynamics of a flexibly mounted flat plate in air. Chapter 2 In-vivo measurement of the extreme shape morphing of bat wings Arnold J. Song, Sharon M. Swartz, Daniel K. Riskin, Jose Iriate-Diaz and Kenneth S. Breuer 2.1 Introduction 2.1.1 Compliant membrane wings The key anatomical feature of bats that distinguishes these flyers from birds and insects is a set of wings that is primarily composed of a compliant skin membrane with low exten- sional stiffness and negligible bending stiffness [2]. The mechanical behavior of a bat’s wing membrane is not only complex, it is quite different from other examples of mammalian skin. Previous studies have looked at the wing membrane’s anatomical structure and constitu- tive properties. For example, Holbrook and Odland [23] closely examined the anatomical structure of wing membrane skin using light and transmission electron microscopy showing that the membrane is composed of a mesh-like network of collagen and elastin fibers. The elastin gives the membrane its “stretchiness” and the collagen limits the strain at high lev- els of stress; the strain-limiting property of the collagen is consistent with the findings of 8 9 Swartz et al. [2], who performed ex-vivo tensile tests of skin specimens from the different regions of the wing. They observed a nonlinear, “J-shaped” stress-strain relationship for the skin characterized by relatively large strains at low stress values, then as the stress level is increased, the skin transitions to a state where the strain limiting properties of collagen would dominate and result in small differential strain for large values of stress. In addition, unlike other mammalian skins, which are typically isotropic, they found that the membrane stiffness is highly anisotropic; the stiffness in the chordwise direction can be many times higher than in the spanwise direction. These studies give us insight into the unique struc- ture and mechanical properties of the membrane skin, but are performed ex-vivo and out of the context of flight. The mechanical compliance of these membrane wings can be more complex due to the strong coupling between the wing shape and aerodynamic forces where small changes in wing geometry can have significant effects on aerodynamic performance. Some of these changes in the geometry are due to active morphing of the wing, e.g. the changing surface area as the wing extends out during the downstroke, while other changes in geometry are due to the aeroelastic coupling between the wing membrane and aerodynamic forces. For example, the aerodynamic forces acting on a compliant wing will cause the wing to push out and increase in camber, which typically lead to increases in the lift and drag on the wing [10]. This feedback between the membrane’s elastic response and aerodynamic forces continues until the two are balanced. This study examines the time-resolved wing shape and motion of three individuals of Cynopterus brachyotis in steady, forward flight (U = 2.6 − 6.6 m/s). Using multi-camera video recordings of the flight kinematics of these animals, we quantify how much the wing membrane expands and contracts during the wingbeat cycle, correlate this behavior with 10 the extension and retraction of the wing bones, and demonstrate the anisotropic strain behavior of the plagiopatagium. 2.1.2 Leading edge alignment Being composed of skin, the wing membrane of a bat is very thin relative to the chord dimension and the aerodynamically sharp profile of the wing’s leading edge makes it sus- ceptible to flow separation. This separated flow, when associated with rapid changes in the angle of attack that are typical of flapping flight, can lead to the formation of strong leading edge vortices (LEVs). The flow acceleration due to these vortices enhances lift from static levels [24] . There is evidence that LEVs are formed and can contribute to large values of lift in insects, birds, and most recently Glossophaga soricina during very low speed and hovering maneuvers [25, 24, 26]. However, the generation of LEVs would be counterpro- ductive for fast, forward flying animals because of the high drag associated with the leading edge separation. The preferred flight mode of C. brachyotis is fast forward flight, therefore we would expect that LEVs would be suppressed to minimize drag, possibly by active control of the leading edge orientation through movement of the thumb. To determine whether C. brachyotis utilizes active leading edge control for the suppression of LEVs, we extract airfoil profiles from the 3-dimensional kinematic data and examine it orientation with respect to the oncoming flow. 11 2.2 Materials and Methods 2.2.1 Steady forward flight measurements We investigated the in-flight wing morphing behavior during steady forward flight of three adult female lesser dog-faced fruit bat s (Cynopterus brachyotis). These experiments were performed in the wind tunnel facility at the Harvard University Concord Field Station. This open-return wind tunnel facility has a rectangular test section, with dimensions: 1.4×1.2×1.2 m (L×W×H) and had wire screens located at the front and rear of the test section prevented the animal from flying into the nozzle in the inlet and the fan in the outlet. For 97% of the wind tunnel test section, the measured flow velocity is within 2.5% of the mean wind speed and the turbulence intensity of the tunnel is less than 1.28% [27]. We capture the in-flight motion and three-dimensional wing shape of the wing by video tracking each individual bat. Each animal was marked with an array of 54 markers dis- tributed over the left wing, as shown in Figure 2.1. Self-adhesive retroreflective markers (3M Scotchlite) were placed on the skin membrane while the skeletal members and body reference points were marked with a non-toxic acrylic white paint. The motion of these markers were recorded by three synchronized, high-speed CMOS cameras (Photron USA, San Diego, CA, USA) recording at a frame rate of 1000 Hz. The time-varying locations of the wing markers in the multiple 2-dimensional video views were digitized and mapped to 3-dimensional space coordinates via Direct Linear Transformation (DLT) [28] using custom MATLAB software. For each experiment run, the marked bat was held by a handler and released from the rear of the test section through a side port. The typical flight path of the bat was against the flow from the rear to the front of the test section with the animal coming to rest at the 12 Figure 2.1: Analysis mesh construction from markers. The bat was marked with an array of 54 markers distributed over the left wing. The surface of the wing was discretized into triangular elements formed by the connective edges between the markers, as shown above. The colors denote the regions considered for analysis: proximal plagiopatagium (red), distal plagiopatagium (blue), propatagium (white), dactylopatagium III-IV (green), and dactylopatagium IV-V (pink). The numbered names of the dactylopatagia refer to the bordering digits. inlet screen. 2.2.2 Wing kinematics The wingbeat cycle is defined by the motion of the wrist relative to the shoulder joint where the maximum and minimum wrist heights mark the beginning and end of the downstroke, respectively. Direct comparison of individual wingbeat cycles can be tricky since the cycle duration will vary slightly from wingbeat to wingbeat. To account for variations in cycle duration and upstroke-to-downstroke ratio, we use the lower and upper reversal points as key points and shrink or stretch each upstroke (lower to upper reversal point) to span from 0 to 0.5 and each downstroke (upper to lower reversal point) to span from 0.5 to 1, which also artificially fixes the upstroke-to-downstroke ratio to unity. This is only for convenience’s sake for the analyses and does not suggest that the two phases of the wingbeat have equal length, in fact the upstroke is typically shorter than the downstroke (see Fig. 2.2 ). Several joint angles were measured to distinguish the biomechanical contribution to the 13 Figure 2.2: Wing surface area. The variation of membrane surface area over several wingbeats at 4.1 m/s. The shaded regions indicate the downstroke portion of the wingbeat. The change in the surface area for a particular region is normalized by the the mid-downstroke surface area, denoted by A0 . The surface areas for all of the wing membrane regions are at or near their minima at the beginning of the downstroke and expand until the middle of the downstroke. The area of the plagiopatagium continues to increase while the areas of the dactylopagatia and propatagium stay nearly constant through the transition from downstroke to upstroke until the middle of the upstroke. Figure 2.3: Three joint angles (elbow and wrist angles between digits IV-V and III-IV) vary over the wingbeat cycle and drive the expansion and contraction of the wing membrane. Over several wingbeats, we observe changes of only several degrees in the angles of the handwing (θ2 and θ3 ). In contrast, the elbow angle (θ1 ) varies over a range of 60 degrees during the wingbeat cycle. The large range of plagiopatagium surface area can be partially attributed to this pronounced elbow motion. However the opening and closing of the elbow joint does not explain why the plagiopatagium area continues to increase when the elbow joint begins to close. This expansion of the wing membrane can only be attributed to aeroelastic effects near the lower reversal point. 14 wing morphing from stretching and deformation due to aerodynamic forcing. The elbow angle (θ1 ), the angle between digits IV and V (θ2 ) and the angle between digits III and IV (θ3 ) (see Fig. 2.3) were calculated from the dot product of a pair of vectors with coincidental a·b origins at the elbow (θ1 ) and the wrist (θ2 , θ3 ) using the definition θ = arccos( |a||b| ). 2.2.3 Skin stretching The surface of the wing was discretized into triangular elements formed by edges connecting markers, as shown in Figure 2.1. For the surface area analysis, the wing was divided into 5 regions according to the natural subdivision of the wing membrane between bony members. Due to the size of the plagiopatagium, the arm wing region is further subdivided into two parts. However, more detailed analyses are possible by interrogating smaller groupings of elements or even individual elements.Each triangular element can be represented by two vectors with the same origin. If these two vectors are denoted by u and v, then the area of the triangle formed by these two vectors is 1 A = |u × v|. 2 The variation of the surface area of each of the membrane regions is expressed as a “surface strain”, i.e. the relative change in area with respect to a reference area as A0 . Ideally, this reference area would be the resting area of the skin. However, it is extremely difficult, if not impossible, to define the resting state and geometry for the animal’s wing. Instead, we choose to use the mid-downstroke membrane area as the reference area, A0 . To investigate the directional preference of the membrane strain, we measure the time evolution of the length of two linear elements oriented spanwise and chordwise in the pla- giopatagium region of the wing (see Fig. 2.6). The linear elements are located in the interior 15 of the membrane such that they do not intersect bones. The true strain,  = log(l/l0 ), is calculated with reference to the minimum lengths, l0 , of the elements. 2.2.4 Camber and angle calculations For conventional airfoils with nonzero thicknesses, the camber line is the curve which lies halfway between the upper and lower surfaces. However, for the bat wing the skin thickness is on the order of 0.1% of the mean chord length. Therefore, the airfoil is approximated as having zero thickness, i.e. the wing’s airfoil profile and the camber line are one and the same. The airfoil shape, or camber line, is determined by the intersection of the three-dimensional reconstruction of the wing surface and a plane. For this analysis, we define two construction vectors: 1) a vector that originates at the shoulder and runs to the wrist, which we define as the leading edge vector, a, and 2) the relative flow velocity at the midpoint of vector a, which we define as vrel , as seen in Fig. 2.4. The intersection plane is then defined by two vrel a vrel × a αle camber line α LE chord TE vrel Bottom view Side view Figure 2.4: Camber extraction. The airfoil shape, or camber line, is determined by the intersection of the three-dimensional reconstruction of the wing surface and a plane. The intersection plane is then defined by two vectors: 1) the cross product of the leading edge vector with the relative velocity, a × vrel , and 2) the relative velocity of the leading edge point, vrel . The angle of attack, α, is defined in the traditional sense as the angle of the chord line relative to the oncoming flow direction and we define a leading edge angle of attack, αle , as the angle between the line drawn from the leading edge to the next fitting point and the oncoming flow direction. 16 vectors: 1) the cross product of the leading edge vector with the relative velocity and 2) the relative velocity of the leading edge point. Generally, this intersection plane lies between between the tracking markers and intersects the lines that connect the tracking markers; the camber line is a third-order polynomial fit to these intersection points. The angle of attack, α, is defined in the traditional sense as the angle of the chord line relative to the oncoming flow direction and we define a leading edge angle of attack, αle , as the angle between the line drawn from the leading edge to the next fitting point and the oncoming flow direction. Large positive values of this angle, αle , would suggest that the flow over the wing is likely separated. 2.3 Results and discussion 2.3.1 Skin stretching behavior We begin with a description of how the area of each wing membrane region varies over the wingbeat for C. brachyotis flying at a moderate speed (U = 4.1 m/s). The membrane “area strain” is normalized with respect to the mid-downstroke area, A0 , due to the consistency of this value over multiple wingbeats and the decreased reliance on interpolation during this portion of the wingbeat. In addition, previous work which utilized a quasi-steady model for the generation of aerodynamic forces assume that the mid-downstroke wing geometry is representative of the wing shape for analysis [8, 29] and the changing area of the wing mem- brane is a result of a combination of aeroelastic effects and the arms and fingers spreading out. Here, we quantify the effect of digit, forearm, and arm motion by measuring the joint angles at the elbow and between the metacarpals and use this information in conjunction with the membrane expansion and contraction to infer the significance of aeroelastic effects 17 during different phases of the wingbeat cycle. We begin in the middle of the downstroke (see Fig. 2.2), where we have set the nor- malized membrane areas are equal to unity. As the wing moves through the middle of the downstroke to the lower reversal point, the plagiopatagium expands to 122-130% of its mid-downstroke area and continues to expand through the transition from downstroke to upstroke and into the beginning of the upstroke (10% of the upstroke portion of wingbeat cycle). After this point, the armwing membrane begins to decrease in area and continues to contract until the beginning of the downstroke (10% of the downstroke portion of the wingbeat cycle). However, this expansion and contraction of the membrane does not occur in phase with the elbow flexion and extension. If the elbow extension and flexion were the only contributors to the size of the membrane then we would expect the maxima and minima of the membrane area and the elbow joint angle, θ1 , to coincide, but we actually observe the membrane reaching its maximum area after the elbow has started to flex (75% of the downstroke) (see Fig. 2.3). This indicates that the plagiopatagium continues to expand even though the elbow is beginning to flex (decreasing θ1 ). There are two reasons that the membrane would continue to expand despite the flexion of the elbow: 1) the aerodynamic forces have increased such that the membrane’s aeroelastic response is greater than its elastic recoil response due to elbow flexion, and/or 2) the membrane stiffness has decreased. It is unlikely that the lift and drag would increase at this stage in the wingbeat, with or without a decrease in the membrane stiffness, because the wing begins to lose forward velocity as it nears the end of the downstroke and actually begins to move backwards during the upstroke. This forward deceleration reduces the wing’s airspeed therefore the lift and drag coefficients would need to increase just to maintain the same area if the elbow was not flexing and the membrane stiffness remained constant. This 18 increase in lift and drag would result from an increase in angle of attack, an increase in the camber, or a combination of the two. Because the elbow is flexing, the more likely reasons that the membrane area continues to increase for a bit longer is a combination of an increase in angle of attack and an increase in camber due to a decrease in membrane stiffness, which could result from a relaxing of the plagiopatagiales. The plagiopatagium clearly expands and contracts to a large degree throughout the wingbeat cycle, but does this behavior change as a function of flight speed? One might suspect that the wing surface area expansion would increase with flight speed due to an increase in dynamic pressure resulting in an increase of aerodynamic forces. But, we when compared the surface area of the plagiopatagium over a range of forward flight speeds, U=2.6, 4.1, and 6.6 m/s (see Fig. 2.5), we found that the deformation behavior of the plagiopatagium was strikingly similar over the speeds tested, despite a six-fold increase in the dynamic pressure (q = 21 ρu2 ) from the lowest to highest flight speeds. The bat’s lift Figure 2.5: Plagiopatagium area variation over normalized wingbeat for one individual. Since the duration of each wingbeat cycle generally varies from wingbeat to wingbeat, we normalize the time scale of the upstroke by the total duration of that particular upstroke (time of upstroke varies from 0 to 1), and likewise with the downstroke (time scale of downstroke varies from 1 to 2). The surface area as a function of the normalized wingbeat is shown for the plagiopatagium. These data are exhibit strikingly similar trends in area over a wide range of speeds (U=2.6, 4.1, and 6.6 m/s) despite a six-fold increase in the dynamic pressure from the lowest to highest flight speeds. 19 requirement is only the animal’s weight, so this does not change with flight speed, but the thrust required to overcome drag (D = CD qA, where CD and A are the drag coefficient and the wing’s planform area, respectively) is speed dependent. Therefore, the lack of significant variation in membrane stretching behavior with an increase in flight speed implies that the wing camber changes very little as a function of speed. This is significant because any increase in wing camber, while increasing lift, would counterproductively increase the drag. Therefore, it is more likely other flight parameters are varied, such as angle of attack, flapping amplitude, and wing twist, are the main sources of the additional thrust needed at higher speeds. But how does the animal keep the membrane area nearly constant with relatively large increases in the dynamic pressure? The stiffness of the membrane either varies or is great enough that large changes in the pressure differential across the membrane result in only small strains. This would be consistent with actuation of the plagiopatagiales, which are chordwise aligned muscles in the wing membrane, that increases the membrane stiffness or the membrane is stretched to the strain limited region of the “J-curve”. In contrast to the plagiopatagium, the area of the handwing does not vary much during the wingbeat cycle (see Fig. 2.5). At its largest, the dactylopatagium has only increased by 8-10% from the midstroke area and reaches an area plateau well before the completion of downstroke (67% of downstrok). From this point until middle of the upstroke, the dactylopatagium maintains a near constant area. The angles between the digits, θ2 and θ3 , also show little variation over the wingbeat cycle. There is a small dip in the value of θ2 to 7 degrees during a short period near the upper reversal point, otherwise the angle is nearly constant, ranging from 11-14 degrees. The angle between digits III and IV, θ3 , also remains at a nearly constant value of 5-7 degrees. The close correlation between the spreading of θ2 and θ3 and the expansion of the dactylopatagia suggests that main mechanism stretching 20 these two regions of wing membrane is the spreading of the fingers under direct muscular control. 2.3.2 Directionally preferential strain CL 1 0.8 Strain, log(l/l0) 0.6 0.4 0.2 0 0 50 100 150 200 250 300 350 Time [ms] Figure 2.6: Anisotropic strain in the plagiopatagium for Cynopterus brachyotis. The strain is calculated for line elements in the plagiopatagium that are not attached to any skeletal members. The length of the line element is normalized by the minimum length. The maximum spanwise strain is nearly 2 times than the maximum strain in the chordwise direction, which implies a great membrane stiffness in the chordwise direction. We see that the wing membrane of the Cynopterus brachyotis undergoes dramatic sur- face area changes during the wingbeat cycle, which is mostly due to the large motions of the skeleton extending the wing in the spanwise direction (as seen in Fig. 2.3). As a result, one would expect that the stiffness in the spanwise direction to be significantly lower in 21 the spanwise direction to allow for the wing extension. We interrogated the strain behavior of two linear element that are roughly aligned with the wing’s chord and span (as seen in Fig. 2.6). The strain is initiated at the beginning of the downstroke and the strain be- havior of the two elements increase at about the same rate. However, halfway through the downstroke, the chordwise strain plateaus while the spanwise strain continues to increase until reaching a maximum value of 0.75 significantly higher than the chordwise maximum strain of 0.48. The similarity in strain behavior in the first half of the downstroke and last half of the upstroke suggests that the membrane behaves as if the stiffness is isotropic, but between middownstroke and midupstroke the membrane stiffness transitions into an anistropic regime. 2.3.3 Leading edge alignment Figure 2.7: Camber extraction. The camber line is defined by the line formed by the intersection of the wing and a cutting plane located at 50% of the distance from the shoulder. In figure (b), the filled circles indicate where the cutting plane crosses through the edges of the triangular elements. The camber line is rotated so that the relative wind velocity vector is aligned with the abscissa. Each colored line indicates the camber line at a particular instant during the downstroke, e.g., the blue camber line is at 20% of the downstroke and the dark green indicates the camber line at 50% of the downstroke. The nonzero camber of the wing results in a difference between the angle of attack (the angle between the velocity vector and the chord line), α, and the leading edge angle (the 22 angle between the velocity vector and a line drawn from the leading edge to a point on the airfoil contour at 20% of the chord length), αle (see Fig. 2.7). Measurement of α and αle were conducted for three individuals of C. brachyotis: individual 1 (3 flights), individual 2 (1 flight), and individual 3 (1 flight). For each flight, mean values were determined from a cropped wingbeat cycle (25% into the downstroke to 25% into the upstroke) to focus on the portion of the wingbeat cycle when the wing is extended. Figure 2.8: Angles with respect to the relative velocity at the wing. The average angle of the leading edge and chord line over several wingbeats for several individual bats are plotted as a function of forward flight speed with error bars indicating the standard deviation from the mean angle. The leading edge angle and the angle of attack with respect to the relative velocity are plotted for a specified wing location (50% of distance from the shoulder to the wrist). Despite the variation in forward flight speed, the mean leading edge angle is small and close to zero. This near zero leading edge angle is maintained even when the angle of attack of the wing is relatively high, which are 11 and 7 degrees at 2.6 and 4.1 m/s respectively. To test whether this species of bat might be actively preventing stall, we have measured the average angle of attack and the angle of the leading edge, which we have defined as the line between the leading edge and the next fitting point. These angles are averaged over at least three wingbeats for three forward flight speeds, U = 2.6, 4.1 and 6.6 m/s, with error bars indicating the standard deviation from the mean angle (see Fig. 2.8). The mean angle of attack for the wingbeat cycle is a decreasing function of forward flight speed, ranging 23 from 17 degrees (U=2.6 m/s) to 6 degrees (U=6.6 m/s). This is not surprising since the lift requirement of the animal, i.e. weight support, is independent of flight speed. From aerodynamic theory, we know that the lift force, L, is a function of the flight velocity squared and the coefficient of lift, CL , is a function of angle of attack: 1 L = ρU 2 ACL (α), 2 where the air density and wing’s projected area are denoted as ρ and A, respectively. Therefore, to maintain level flight at a higher speed, the angle of attack necessarily decreases. In contrast, the mean leading edge angle remains small and close to zero. This near- zero leading edge angle is maintained even when the angle of attack of the wing is relatively high, with values of 17 and 10 degrees at 2.6 and 4.1 m/s, respectively. This suggests that the animal’s wing motion during the downstroke is tailored such that the leading edge is nearly always aligned with the oncoming flow. This alignment of the leading edge to the relative velocity may be important to delay or prevent separation over the armwing, therefore preventing a decrease in lift and increase in drag which accompanies the onset of separation. 2.4 Concluding remarks In this chapter, we have presented results of in vivo measurement of the time varying wing shape of three C. brachyotis flying at speeds ranging from U = 2.6 − 6.6 m/s. We show how the plagiopatagium continues to expand even though the wing begins to fold suggesting a relaxation of the plagiopatagiales, which are the muscles that run chordwise in the armwing. In addition, the nearly identical membrane stretching behavior over the range 24 of speeds indicates that the membrane stiffness is large enough that the six-fold increase in dynamic pressure does not result in increased membrane area during the downstroke. The work presented in the following chapter shows that membrane prestrain may be responsible for the speed independent deformation behavior. This pre strain could be achieved through actuation of the plagiopatagiales. In addition to modulating the membrane stiffness by the plagiopatagiales, the wing’s leading edge is nearly always aligned with the oncoming flow that suggests active alignment of the leading edge through movement of the thumb or that the wing motion of this animal varies such that the leading edge is always aligned with the oncoming flow. This alignment of the leading edge with the oncoming flow indicates that flow separation is being actively prevented for the preferred level, cruising flight mode of C. brachyotis. For low Reynolds number flight (Re < 1 × 105 ), the active control of the leading edge portion of the wing may be particularly important not only for flapping flight but steady flight for small scale vehicles, i.e. MAVs, to prevent the formation of separation. Chapter 3 Aeromechanics of membrane wings with implications for animal flight Arnold J. Song, Xiaodong Tian, Emily Israeli, Ricardo Galvao, Kristin L. Bishop, Sharon M. Swartz and Kenneth S. Breuer. AIAA Journal. Vol. 46, No. 8, pp. 2096-2106, 2008. 3.1 Introduction Wings composed of thin, compliant skin membranes are used in the natural world by several vertebrate species, perhaps most notably by flying and gliding mammals, such as bats, flying squirrels, and marsupial gliders. These animals exhibit extraordinary flight capabilities with respect to maneuvering and agility that are not observed in other species of comparable size. Birds, which have been studied extensively, have relatively rigid wings with limited degrees of freedom, while insect flight, which occurs at much lower Reynolds numbers, can be characterized by the relatively simple articulated flapping motion of effectively rigid wings[30, 31, 32]. In contrast, bats have a high degree of articulation in the wing (the elbow, wrist, and finger joints) [9] and, more relevant to the current investigation, the wings of bats are composed of a highly anisotropic, compliant skin membrane (Figure 3.1). These morphological features may be key in enabling bats to fly in such a remarkable fashion, yet little is known of the aerodynamic performance of vehicles with either highly articulated or 25 26 Figure 3.1: Examples of mammalian flight. Photographs of Cynopterus brachyotis (left) and Pteropus vampyrus (right, photo courtesy of Richard Wainwright) during flight. The unique wing morphology of bats include their skeletal structure and thin wing membranes. Note the large degree of camber in the wing. The unique aerodynamic properties of a compliant membrane wing are believed to play significant role in the high degree of maneuverability exhibited by bats in flight. compliant wings. The study of the full complexity of mammalian flight is challenging but is nevertheless considerable progress has been made recently. Most studies have focussed the animal’s kinematic motion [33, 34], their overall lift and drag characteristics [9, 33, 7] or on the wake characteristics behind the animal[33, 35]. However, a significant challenge in working with live animals (and something that is somewhat unfamiliar to the traditional aerodynamics community), is the requirement that the testing be accomplished in a safe and humane manner. Furthermore, it is difficult, if not impossible, to isolate the individual contribution of the various, interdependent aspects of an animal’s morphology to its overall aerodynamic performance. One solution to these difficulties is to test, using more traditional engineering models, key features of the biological system with the goal of understanding how each feature influences the aerodynamic performance of the complete animal. Features that might be isolated include the use of low aspect ratio (LAR ) wings, the thin compliant membrane lifting surface and the low Reynolds number regime (Re < 2 × 105 ), where many complex aerodynamic phenomena, such as transition to turbulence and laminar separation, 27 are present, yet extremely hard to predict [36, 4]. Focussing on the effects of low Reynolds number, Shyy et al.[5] utilized XFOIL, a two- dimensional analysis code, in a computational comparison of the low Reynolds number performance of several airfoils, with a particular interest in the effects of airfoil thickness and camber. For Reynolds numbers based on chord length, Re = 7.5 × 104 − 2.0 × 106 , the well known NACA 0012 and the CLARK-Y airfoil profiles served as verification of the XFOIL predictions and as the baseline for comparison with two airfoils optimized for low Reynolds number flight. The thinner and more cambered airfoils yielded higher lift-to-drag 3/2 ratios and power efficiency, CL /CD , in the range of Reynolds numbers tested. This was particularly striking at the lowest Reynolds number, Re = 7.5 × 104 , where the power efficiency of the designed airfoil was nearly 3 times that of the NACA 0012. Additionally, the authors examined the effects of the deformability of an airfoil. The top of a CLARK-Y airfoil was modeled in XFOIL as a flexible membrane which adopted an equilibrium shape for a particular angle of attack and freestream velocity. The authors concluded that the aerodynamic performance of the flexible and rigid airfoils were comparable. However, the airfoil with the flexible membrane top was less sensitive to oscillations in the freestream velocity. Experimental studies by Mueller and coworkers [37, 38] reported on experiments with rigid, thin LAR airfoils of varying camber and planform shapes for Re ≈ 105 . These series of experiments characterized the effects of aspect ratio and camber on lift, drag, and aerodynamic efficiencies and suggested that the wingtip vortices are responsible for the non-linear regions of the lift curve for wings with aspect ratios less than 1.25. These flow structures play a role similar to the tip vortices on a delta wing, providing vortex lift, particularly at high angles of attack, contributing to an increase in maximum lift, CLmax , 28 and to a delay of stall to higher angles of attack, α. There have been several previous studies of the aerodynamic characteristics of the chang- ing shape of flexible yet inextensible airfoils, i.e. sails, using thin-airfoil theory [39, 40, 41]. Thwaites and Nielsen find that the equilibrium shape is only a function of a parameter, 2 λ = 2ρU∞ c/T, (3.1) using Thwaites’ definition. However, the results of these studies were limited to low angles of attack and small values of camber which comprise only a small subset of conditions for compliant wings. It is experimentally shown in this paper that compliant wings can reach camber values of over 10% of the chord length. Additionally, in vivo camber measurements of bats flying in a wind tunnel, not presented in this paper, have shown that the camber of the skin membrane reaches similar values of near 10% during the downstroke, the portion of the wingbeat during which a majority of the lift is produced. Subsequent studies have coupled numerical studies with experiments focussing on the aerodynamics of inextensible, flexible airfoils[42, 43]. Lorillu et al. show good agreement between computations and experiments, however the computed results require that separation point to be empirically determined. Several studies on the effects of compliance have reported on numerical techniques suit- able for this highly coupled aeroelastic problem, and on improved performance of flexible membrane wings for micro air vehicles (MAVs) [36, 5, 44, 45, 46, 47]. In this paper, we continue this exploration, presenting both analytical and experimental results concerning both the static deflection and the unsteady behavior of LAR membrane wings. In particu- lar, experimental observations have been correlated with some relatively simply analytical 29 predictions that help to explain how the aeroelastic coupling is mediated. 3.2 Analytical Considerations z ds x Leading T T Trailing Edge β β Edge U Fn ε0=0 0.1 ε0=0.1 Camber, z/c ε0=0.2 0.05 ε0=0.3 0 0 0.01 0.02 0.03 0.04 0.05 Weber number, We Figure 3.2: Model for compliant membrane wing subject to aerodynamic loading. Left: the membrane is modeled as a one-dimensional projection onto the span-normal plane, subject to a uniform load. Right: Camber as a function of the nondimensional aerodynamic loading, characterized by the Weber number, We = F c/Et. The camber increases rapidly at low aerodynamic loads, flattening as the load increases. Increased values of pre-strain,0 , result in a more gradual deflection. Before describing the experiments, it is instructive to explore the predictions of a very simple theoretical model for a compliant wing subject to aerodynamic loading. The static aeroelastic deformation of a compliant membrane represents a balance between the aero- dynamic forces generated by the airfoil shape and the tension in the membrane. A simple 30 geometric model capturing this is presented in figure 3.2a. The model considers the two- dimensional projection of a membrane of thickness, t, and undeflected length, c, pinned between stationary leading and trailing edges. When subjected to aerodynamic loading, denoted by Fn in figure 3.2a, the membrane will deflect in the chord-normal direction, z, with contact angles, βLE and βT E . Any aerodynamic forces generated by the leading and trailing edge structures are assumed to negligible in this model. For the membrane to be in static equilibrium, the membrane tension, T , per unit span must balance the aerodynamic forces, Fn , generated by the airfoil: Fn = 2T sin(β) = c∆p , (3.2) where ∆p represents an idealized uniform pressure difference between the lower and up- per surfaces of the wing. Increases in the membrane tension result from the elongation or stretching of the membrane when subjected to aerodynamic loading. Additionally, the membrane has a contact angle, β, with the leading and trailing edge. Therefore, only the vertical component of the membrane tension at each end, T sin(β), supports the aerody- namic loading. Taking ds to be a differential element of membrane length, the elongation of the membrane, ∆c, is given by Z TE ∆c = ds − c, (3.3) LE where the integral on the right hand side of the equation is the arc length of the deflected membrane. For the simple model presented in this paper, the membrane is assumed to behave like a linearly elastic material with a modulus, E, which is independent of strain. In reality, both biological and man-made membranes are known to exhibit nonlinear J-shaped 31 stress-strain behaviors for which E is strain dependent [48]. This behavior is ignored for the current analysis. Therefore, we express the membrane tension per unit width as only a function of the elongation length, ∆c: Z TE  ∆c Et T = T0 + Et = T0 + ds − c , (3.4) c c LE where the initial tension present in the undeflected membrane is denoted as T0 . Substituting this result into equation (3.2) and normalizing the aerodynamic loading by the product of the modulus and membrane thickness we obtain an expression for the non-dimensional aerodynamic loading, expressed as a Weber number, We:   Fn ∆pc T0 ∆c We = = =2 + sin(β) = 2 (0 + ) sin(β). (3.5) Et Et Et c For simplicity, we assume that the membrane adopts a parabolic shape, z(x) = 4zmax x(1 − x). This is the exact solution for a membrane subjected to uniform loading, and is an assumption that will be justified by our experimental results, in section 3.4.1. With this as- sumption, the membrane length is thus only a function of camber (defined as the maximum extension of the membrane, zmax , which occurs in this model at mid-chord), and hence, using equation (3.4) we see that the membrane tension is also solely a function of camber. The contact angles at the leading and trailing edges, βT E and βLE , are generally different. However, the assumed parabolic airfoil shape is symmetric, therefore βLE = βT E . Using equation (3.5), we can determine the aerodynamic loading that can be supported by the membrane geometry with a prescribed camber value. For each value of the camber, zmax , the elongation length, ∆c is calculated by integrating the contour of the camber line (equa- tion (3.3)). This elongation length determines the tension which, with the contact angle, 32 β, determines the aerodynamic loading that can be supported by membrane. Figure 3.2b presents the membrane camber as a function of the aerodynamic load, for varying values of initial pre-strain, 0 = T0 /Et. The camber-enhancing effect is most pronounced when there is no initial pre-strain, i.e. 0 = 0, in which case we see that initially the membrane quickly adopts camber for very small aerodynamic loads. As the membrane deflects further, the camber dependence on aerodynamic loading weakens and reaches an asymptotically linear dependance. The rapid initial rise in camber can be attributed to the fact that the undeflected membrane offers no resistance to any aerodynamic load (since a membrane has zero bending stiffness) and thus starts to camber immediately. At low loads, the small curvature of the membrane results in small values of both sin(β) and membrane elongation, ∆c. Therefore, the chord- normal contribution of the tension is small and the membrane quickly “balloons” out. As the membrane deflects further from the chord line, this effect diminishes. This behavior is strongly dependent on the membrane pre-strain, 0 , which allows the membrane to support larger aerodynamic loads at lower values of camber. Therefore, the membrane camber shows a weaker dependence on the Weber number for increasing pre-strain values, 0 . 3.3 Experimental Procedure All experiments were conducted in the Brown University low-speed wind tunnel. This is a closed-return facility with a constant-speed, variable-angle axial fan. The test cross section measures 61 by 58 cm (width by height) and has excellent flow quality, with a freestream tur- bulence measured to be less than 0.1% (0.1 - 10,000 Hz). The experiments were conducted at freestream velocities ranging from 8-21 m/s (measured by a pitot tube connected to a Baratron pressure transducer (MKS model 398HD), corresponding to Reynolds numbers 33 based on the chord length of Re = 7 × 104 − 2 × 105 . 3.3.1 Compliant Membrane Models Compliant membrane half-wings were designed and manufactured for these experiments. The rectangular wings were composed of a compliant latex membrane held between two stainless steel posts mounted on an aluminum base plate. The steel posts defined the leading and trailing edges (figure 3.3). The membrane material was inserted through a slit aligned with the centerline of the post and then secured using spring steel clamps. With the clamps in place, the leading and trailing edges were approximately parabolic in shape with a maximum thickness of 3.7 mm. Using two sets of posts, with lengths 6 and 9 cm, Latex Membrane with chord separations of 10 and 13 cm, we were able to fabricate wings with aspect ratios of AR = 0.9, 1.4 and 1.8. Side view Test section ceiling 2-axis force balance LE/TE clamps s/2 Latex membrane U∞ c End-on view Deflected Leading membrane edge α Trailing U∞ Chord edge Figure 3.3: Experimental setup. The compliant membrane wings were composed of thin latex sheeting clamped at the leading and trailing edges. The edges at the root and tip were unconstrained and could deflect freely. The leading and trailing edge clamps were mounted to a plate which was attached to the sting arm of a two-axis force balance. 34 Each half-wing was assembled by aligning the inner edges of the posts with reference lines on the membrane. The distance between the reference lines was equal to or shorter than the distance between the inner edges of the posts resulting in an unpre-strained or pre-strained membrane, respectively. Thus, for the model with a pre-strain value of 4%, the membrane was stretched so that it was 4% longer than its original relaxed state, i.e. the distance between the inner edges of the posts. We should note that the assembly technique did not allow for precise control of the tension in the wing membrane. As a result, direct comparisons between wings of different membrane thicknesses are to be taken as qualitative, since the cambering behavior is a strong function of the membrane pre-strain, 0 , therefore slight variation in 0 can have a marked effect on the cambering behavior, as discussed previously. 3.3.2 Stereo Photogrammetry A stereo photogrammetry system was developed to measure the instantaneous deflections of the membrane wings. Two high-speed CMOS cameras (Photron PCI-1024) were placed outside the wind tunnel test section with an approximately 60 − 70◦ angle of separation between the respective optical axes. The cameras have a resolution of 1024 x 1024 pixels, which in conjunction with the imaging optics (Nikon 60 mm Micro Nikkor lenses) provided a spatial resolution of 0.20 mm per pixel. For this series of experiments, the unsteady motion of the membrane was recorded for 1.0 seconds at 1000 Hz (1000 frames) at each speed and angle of attack. Custom software, using Direct Linear Transformation (DLT), was used to recover the coordinates of each marker location in the 3D object space from two 2D images at each time step.[49] A calibration cube consisting of a three-dimensional grid of markers was used to 35 generate the calibration coefficients needed to reconstruct coordinates in the object space via DLT. Using a test object that was moved in a controlled manner, we determined that the measurement uncertainty of photogrammetry system was less than ±35µm for in-plane displacements and ±40µm for out-of-plane displacements, with the reference plane taken to be the chord line. 3.3.3 Lift and Drag Measurements In a separate series of experiments the wing was mounted on a gimballed, two-axis force balance for measurement of the lift and drag. The aerodynamic forces were determined from the differential output voltages from load cells (OMEGA LCKD series subminiature compression load cells), ∆Vx , ∆Vz , placed at the ends of platforms corresponding to each measurement axis. The differential arrangement minimizes the common-mode noise due to vibration of the apparatus. The entire force balance assembly was mounted to a turntable which allowed adjustment of the angle of attack. As the angle of attack was varied, the orientation of the orthogonal platforms remained fixed with respect to the chord line of the wing. Force Balance Calibration At any given speed and angle of attack, the forces in the streamwise and spanwise directions may be determined from a linear combination of the four load cell voltages. Ideally, the output voltages of two load cells at the ends of each of the measurement platforms would correspond in an equal and opposite direction in response to a force. However, each load cell had slightly different gain, and small misalignments between the platforms of the force balance were unavoidable. For this reason, the voltage differentials, ∆Vx and ∆Vz , of both 36 load cell pairs were used to determine both Fx and Fz . The chordwise and chord-normal directions are denoted as x and z, respectively. The relationship between the output voltage of the load cells and forces can be expressed in matrix form:       Cxx Cxz   ∆Vx   Fx    =  (3.6)      Czx Czz ∆Vz Fz The coefficient matrix, C, was determined from a series of calibration measurements performed by applying static forces in both the chordwise and chord-normal directions. For each combination of applied loads, [Fx ; Fz ], the differential load cell voltages, ∆Vx and ∆Vz , were measured. The calibration coefficients, Cij , were then determined using a generalized linear regression of the overdetermined linear system: A·C=F (3.7) where      ∆Vx,1 ∆Vz,1 0 0   Fx,1        .. .. .. ..  ..         . . . .   Cxx   .                   ∆Vx,N ∆Vz,N 0 0    Cxz    Fx,N  A= ,C =  , and F =   , (3.8)           0 0 ∆Vx,1 ∆Vz,1 Czx   Fz,1                  .. .. .. ..     .  ..    . . . .   Czz            0 0 ∆Vx,N ∆Vz,N Fz,N and the number of calibration measurements is denoted by N . We can decompose the 37 matrix of measured voltages, A, using singular value decomposition (SVD): A = USVT (3.9) The dimension of the voltage matrix, A, is 2N × 4. We define S to be a 2N × 4 diagonal matrix of singular values of the orginal matrix A. U and V are then defined to be matrices of dimensions 2N × 2N and 4 × 4 respectively, which have orthogonal columns, so that: UT U = I, and (3.10) VT V = I. (3.11) The calibration coefficients, C, are then solved for using: C = A−1 F = VS−1 UT F. (3.12) Because the inverse of S is used to determine C, artificially large coefficient values may result from the very small entries of S, i.e. singular values of negligible importance. We employed a threshold technique, in which which the spuriously large values of S−1 were set to zero. A value was considered to be spuriously large if it was 50 times higher than the minimum value of S−1 . This calibration method was found to be accurate within 1.5% for loads larger than 1 N, and no worse than 5.0% for calibration loads less than 1 N. The uncertainty in the force measurement is the most significant contributor to the overall uncertainty of the measurement apparatus. Uncertainties in the velocity, temperature, and ambient pressure measurements were all found to be less than 0.5%. Therefore, the uncertainties in the lift and drag coefficients for Re = 1.4 × 105 were found to be no worse than 5% for values of 38 CL and CD less than 0.5 and less than 1% for larger values of CL and CD . 3.4 Results and Discussion These experiments represent our first efforts in an ongoing investigation of compliant mem- brane wing aerodynamics in the context of mammalian flight. As a result, these initial experiments were deliberately broad in their scope in an attempt to identify as many of the critical features of compliant membranes salient to mammalian flight and MAV design. The results are organized as follows: In the following section, the static deflection of the mem- brane is discussed. Following this, lift and drag measurements are presented in section 3.4.2. Finally the dynamic motion of the membrane is presented in section 3.4.3. 3.4.1 Static membrane deflection For a nonzero value of the angle of attack, α, a compliant wing deflects and adopts a cambered shape. Stereo photogrammetry was used to reconstruct the three-dimensional geometry of the wing at several wind speeds and angles of attack. A typical time-averaged shape is shown in figure 3.4 and one sees several characteristic features, including the overall membrane deflection due to aerodynamic loading, as well as enhanced deflection at both the root and wing tip. This last feature is discussed is more detail in the next section. As the load increases (either due to increased air speed or increased angle of attack), the membrane stretches further, adopting higher camber which results in increased lift (according to thin airfoil theory, the lift curve is shifted uniformly upwards as the camber increases). Eventually, the membrane tension balances the aerodynamic load, resulting in the equilibrium wing shape illustrated in figure 3.4. At these moderate angles of attack, the wing deflection was observed to be quite uniform 39 0.15 Camber, z/c 0.1 0.05 0 0 1 0.5 0.8 0.6 Span, y/b 0.4 0.2 Chord, x/c 0 1 Figure 3.4: Static deformation due to aerodynamic loading. The deformed membrane geometry is shown for α = 24◦ , Rec = 1.4 × 105 . The static surface geometry of the membrane is obtained from the time-averaged position of marker array. Note that the deflection of the membrane tip is greater than the root deflection due to the presence of a tip vortex. ε0=0.04 ε0=0.11 −1 10 Camber, z/c 1 1 −2 10 −3 10 −3 −2 −1 10 10 10 Weber number, We Figure 3.5: Maximum camber as function of aerodynamic loading. The maximum camber is plotted as a function of the Weber number, We = 12 ρU 2 cCL /(Et), and compared with the theoretical model. Data presented is a composite assembled from several wings measured at 10 speeds (U = 7 − 16 m/s) and 12 angles of attack (α = 4 − 34◦ ), The wings represened two pre-strains (0 = 0.04, 0.11), and two thickness (t = 0.10, 0.25 mm). We only plot data for low angles of attack (i.e. pre-stall) where the lift slip is approximately linear, and behavior is similar to that predicted by classical thin airfoil theory. For such cases, the lift coefficient can be reasonably estimated to be CL = 2πα. There is excellent agreement between the experimentally observed camber behavior and the analytical model. 40 along the span (with the exception of the root and tip regions), and well-approximated by the parabolic assumption made in the theoretical model presented earlier. That model suggests that the maximum camber of a membrane wing is completely characterized by two parameters, the pretension, 0 , and the Weber number, CL 21 ρU∞ 2 c We = , (3.13) Et where c and t are the chord length and membrane thickness, respectively. Figure 3.5 presents data from wings of two different latex thicknesses (t = 0.1 and 0.25mm), fabricated with two values the pre-strain (0 = 0 and 0.05), and tested over a wide range of freestream velocities (U∞ = 7 − 16 m/s) and angles of attack (α = 4 − 34◦ ). With a few exceptions, we see that there is excellent agreement between the experimentally-measured camber and the analytical predictions. In order to calculate the Weber number from the measurements, we have restricted ourselves to incidence angles below that where stall was observed (see section 3.4.2), and assumed that the lift force is proportional to the angle of attack: CL ∝ α. For a two-dimensional airfoil, the constant of proportionality is 2π. However, for lower- aspect ratio wings, the lift slope will be smaller, whereas the effect of compliancy will be to increase the lift slope. As a suitable compromise and for simplicity, we used a lift slope of 2π, and indeed, we see (figure 3.5) that this choice matches the data extremely well, with the exception of a few outliers at relatively low values of the Weber number. To compare experimental camber measurements with theory, the measurements were fit to the theoretical model by adjusting the value of the pre-strain, 0 . As noted earlier, we were not able to accurately measure the pre-strain during wing fabrication. However, it was found that the best-fit data corresponded to pre-strains of 0 = 0.04 and 0.11, which 41 were larger than the assumed values of zero and 0.05, although the discrepancies are well within the bounds of our fabrication uncertainty. The data outliers at low Weber number correspond to two regimes. The data points falling below the theoretical predictions are at the lowest speed tested, and for these cases, we conclude that the lift force generated was lower than that predicted by our simple linear lift slope assumption, and the camber is consequently lower, resulting in the poor agreement with the theory. The reason for this lower camber might be due to some non-ideal bending stiffness at low aerodynamic loading which is not accounted for by the membrane theory. In contrast, the measurements that lie above the predicted values (both both pre-strains tested) are from measurements taken at high speeds, but at low angles of attack. Here, the lift generated is likely higher than what would be predicted by the linear lift slope assumption - a conjecture that is strongly supported by the lift force measurements at high speeds (discussed in the next section). A key component of the favorable lift characteristics of a low-aspect-ratio wing is the influence of tip vortices over a greater portion of the wing. These tip vortices generate low pressure regions over a significant portion of the wing, a phenomenon distinct from wings with larger aspect ratios where these same tip vortices are confined to a smaller proportion of the lifting surface. Thus, for LAR wings, a greater fraction of the lift generated is due to tip vortex effects. In addition, the downwash from the tip vortices provides the chord- normal momentum needed for the flow to remain attached, thereby allowing a LAR wing to operate at higher angles of attack than would be possible for a wing of higher aspect ratio at these moderate values of Reynolds number (Re ≈ 60 − 130 × 103 ) [6]. Although in the present experiments we cannot measure the pressure or the tip vortex structure, the shape of the membrane provides a window into the pressure distribution over the wing. In this half-span model, the tip and root were unconstrained (figure 3.4) and 42 0.085 Root Tip Normalized Deflection, z/c 0.08 0.075 0.07 0.065 0.06 0.055 20 25 30 35 Angle of Attack, degrees Figure 3.6: Comparison of Tip and Root Deflection. The influence of the tip vortex is particularly strong at low angles of attack as indicated by the increased tip deflection, relative to the root. However, the effect is diminished at stall. we observed a significant difference between their deflections (figure 3.6). These differences qualitatively indicate the strength of the tip vortex, and its influence on the inboard portion of the wing. As the angle of attack rises, it is presumed that a strong vortex is generated at the leading edge. As the angle of attack is increased even further, this leading edge vortex will become unstable and burst. This proposed bursting event leads to an increase in pres- sure over the inboard region of the wing therefore weakening the tip vortices by decreasing the upper and lower surface pressure differential. The combination of the decrease in vortex strength and core size translates into a diminished contribution of vortex lift to the overall lift of the wing. We see this reflected in the relative deflections of the tip and root in fig- ure 3.6. Prior to stall, the wing begins to experience a decrease in lift slope corresponding to a change in the geometry. The inboard camber continues to increase presumably due to the increasing potential lift and drag. However, as the angle of attack is increased even further, the tip deflection begins to decrease or deflate, reflecting a weakening of the tip vortex. This deflation of the tip continues with increasing angle of attack until a sharp drop 43 in the tip deflection is observed. This tip deflation is strongly correlated with the drop in CL and the lift slope at these higher angles of attack (discussed in the next section), and is presumed to be due to the bursting of the tip vortex at the onset of separation. These conjectures need to be confirmed with more detailed flow measurements - an experiment planned for the near future. 3.4.2 Lift and Drag performance 1.5 1.5 Rigid 0.25 mm 0.10 mm 1 Drag coefficient, CD 1 L Lift coefficient, C 0.5 0 0.5 −0.5 Rigid 0.25 mm 0.10 mm −1 0 −20 −10 0 10 20 30 40 50 60 −20 −10 0 10 20 30 40 50 60 Angle of Attack, degrees Angle of Attack, degrees (a) Lift performance (b) Drag performance Figure 3.7: Aerodynamic Performance of Rigid versus Compliant Wings. The lift and drag coefficients of rigid and compliant wings of aspect ratio 1.4 are compared over a range of angles of attack. The compliant wing exhibits significant lift enhancement at low angles of attack and enhanced maximum lift. Both effects are due to the camber induced by the aerodynamic loading. The stall characteristic of the thinner and more compliant wing is much less severe than the less compliant membrane wing and the rigid wing, due to automatic de-cambering at the onset of separation. The drag of the compliant wings is uniformly higher than that of its rigid counterpart. (Re = 1.4 × 105 ) The typical aerodynamic performance of the compliant wing is shown in figure 3.7, which shows the lift and drag coefficients at different angles of attack for three different wings - a rigid wing (a stainless steel sheet), and compliant wings of two thickness (t = 0.10 and 0.25mm). The rigid wing exhibits lift behavior typical of thin, low-aspect-ratio wings (figure 3.7a). The lift slope is nearly constant until a rather abrupt drop in lift, i.e. stall, at 44 higher angles of attack (between 20-30 degrees). In contrast, the compliant wings exhibit several behaviors that are quite different. The less-compliant wing (t = 0.25 mm) has a higher lift slope and reaches a higher value of CLmax than the rigid wing. Also, stall is delayed until much higher range of angles of attack, between 30-40 degrees. The most compliant model (t = 0.10 mm) continues this trend, reaching an even higher CLmax and exhibiting a very soft stall behavior in the range of α = 30 − 40 degrees, before lift is gracefully lost at angles in excess of 40 degrees. The drag performance (figure 3.7b) also highlights the differences between the three models. The rigid wing drag performance can be characterized as a symmetric drag bucket with a slight hiccup at the onset of stall. Both compliant wings generate higher drag than the rigid wing over the entire range of angles tested, with the rise in drag becoming more pronounced as the compliance increases. These behaviors are to be expected, given what we have observed about the wing defor- mation as the angle of attack increases. As the incidence angle rises, the camber increases, hence the composite lift curve will have a lift slope that is higher than its rigid counterpart, because it is “sampling” the lift from a family of wings each with increasing camber. The higher lift is necessarily accompanied by higher (induced) drag, although it is also possible that the more compliant wings also exhibit more drag due to losses associated with unsteady membrane motion. Our theory predicts that a compliant wing will deform dramatically for low levels of aerodynamic loading, i.e. small values of the Weber number, and that the rise in camber versus α slows as the aerodynamic load continues to increase (figure 3.5). This is confirmed by the rapid increase in CL for the thinner latex membrane (t = 0.10 mm), and the lower value of the lift slope at higher angles of attack. The soft stall behavior may also be under- stood by recognizing that as the laminar separation bubble begins to grow, the pressure on 45 the upper surface rises leading to an automatic decambering of the membrane (figure 3.6), backing the wing away from fully separated flow. This “aerodynamic feedback” allows the wing to operate at higher angles of attack, preventing a sharp decrease in lift as the angle of attack is increased. At large angles of attack (> 20◦ ), flow separation for these models likely originates from the leading edge (since the wings have very thin leading edges), and the angle that the leading edge makes to the flow (for a fixed geometric angle of attack) will increase as the speed (and camber) increase. Thus, it is possible that a leading edge separation at low speed might be alleviated at higher speed as the effective angle of attack at the leading edge is decreased. In the case of bat flight, which motivates these studies, the animals are observed to adjust the angle of the leading edge during the wing beat cycle, by directing the propatagium (the flap of skin in front of the shoulder, elbow and wrist) with their thumb [50]. We conjecture that the ability to adjust this leading edge flap, coupled with the adaptive cambering of a compliant membrane both serve to suppress leading edge separation over a wide range of speeds and angles of attack. Too much camber in thin, rigid airfoils can result in a decrease in the lift curve slope. Null and Shkarayev [51] conducted experiments on the effect of camber on the aerodynamics of cambered rigid wings at low Reynolds numbers (Re = 5.0 × 104 − 1.0 × 104 ) and reported that wings with 12% camber showed a dramatically reduced lift curve slope when compared with wings having 3%, 6% and 9% camber. In contrast, we do not observe a reduced lift slope in the present results, despite camber values in excess of 20% with the most compliant model tested (t = 0.10 mm). In addition, the more compliant model, which achieved larger cambers, demonstrated higher values of CLmax - in stark contrast to the rigid, cambered wings of Null and Shkarayev which showed more than a 40% decrease in CLmax for a wing with 12% camber. However, the increased lift of the more compliant wing comes with 46 a price. As one would expect, the more cambered wing had a larger drag penalty and subsequently the lift-to-drag ratio suffered as result. It is notable that the two compliant wings, which reached values of maximum camber of zmax = 10.5% and 23.8% (t = 0.1, 0.25 mm respectively) exhibited relatively similar lift characteristics between α = 10 − 20◦ , with the only difference between the two wings being the drag performance. The ability of the membrane wings to maintain lift at high angles of attack despite very large values in camber, coupled to their more gentle stall characteristics, suggests that these configurations may have significant advantages for vehicles operating in gusty conditions, or with highly unsteady flight patterns - both characteristics of biological flight. Hysteresis effects If one looks carefully at the most compliant membrane lift and drag data, we notice that the symmetry around α = 0 appears to be disturbed. This is not experimental inaccuracy, but is rather due to hysteresis in the membrane deformation. The wing is symmetric, and so has identical aerodynamic performance at both positive and negative angles of attack. At zero angle of attack, the wing, theoretically, has zero camber. However, this is an unstable condition, and one that is never observed in practice, since the slightest perturbation will cause the wing to break symmetry and to adopt camber. Figure 3.8 illustrates this behavior in both the lift and drag coefficients, measured by starting at a large negative angle of attack, increasing through zero and on to large positive angles, and then coming back down again. One sees that the wing maintains camber as the incidence angle crosses through zero. At the angle at which the lift finally reaches zero, the wing “snaps” through and adopts the opposite camber. As one would expect, the rigid wing does not show any hysteresis around the zero angle of attack point. 47 0.4 0.6 0.35 0.4 0.3 Lift coefficient, CD L Lift coefficient, C 0.2 0.25 0 0.2 −0.2 0.15 0.1 −0.4 0.05 −0.6 0 −10 −5 0 5 10 −10 −5 0 5 10 Angle of Attack, degrees Angle of Attack, degrees (a) Lift performance (b) Drag performance Figure 3.8: Hysteresis in lift and drag due to compliance. The lift and drag of rigid and compliant wings are compared near zero angle of attack. The wings were tested from negative angles, through zero to positive angles, and then back again. While the CL and CD of the rigid wing are independent of the direction of testing, the compliant wing exhibits a significant hysteresis in both CL and CD , due to the “snap-through” behavior of the membrane as it passes through α = 0. (Re = 1.4 × 105 ) Although not shown here, both rigid and compliant wings exhibit a similar pattern of hysteresis around the stall point, common for low-Reynolds number wings. CL (α) on ascent was observed to lie below CL (α) on descent. The effect was more pronounced for the compliant wing. We assume that in the post-stall state the flow is characterized by leading edge separation due to the bursting of a leading edge vortex. As α decreases from the post stall state, this leading edge vortex will reform due to turbulent reattachment. For the rigid wing, persists to a slightly lower angle on descent. However, in the case of the compliant wing, the ability of the wing to adjust its camber leads to a more subtle transition over the stall region. However, this has not been confirmed experimentally and needs to be explored with more detailed measurements. 48 1.5 1 1 Drag coefficient, CD 0.5 0.8 L Lift coefficient, C 0 0.6 −0.5 0.4 −1 AR=0.9 0.2 AR=0.9 AR=1.4 AR=1.4 AR=1.8 AR=1.8 −1.5 0 −20 −10 0 10 20 30 40 50 60 −20 −10 0 10 20 30 40 50 60 Angle of Attack, degrees Angle of Attack, degrees (a) Lift performance (b) Drag performance Figure 3.9: Effect of Aspect Ratio on Aerodynamic Performance. For each of the configurations shown, the membrane pre-strain was zero (0 = 0), the membrane thickness (t = 0.15mm) and the dynamic pressure per unit span (ρU 2 c) is the same for all of the wing configurations. The lift slope (a) decreased with the aspect ratio, exhibiting behavior consistent with finite wing theory of rigid airfoils. For small angles of attack (α < 15◦ ), the drag behavior was nearly the same for the three aspect ratios. However, for higher angles of attack, the drag increased with the aspect ratio. This may have implications for gliding mammals (AR ≈ 1) which favor these higher angles of attack in flight. Effects of wing aspect ratio Figure 4.6 presents the lift and drag performance of three compliant wings of varying aspect ratios (AR = 0.92 − 1.78). All three wings had the same membrane compliance and very comparable values of pre-strain. We see that their overall behavior is similar, exhibiting the rapid rise in CL at low values of α, followed by a nearly linear increase in CL until they experience a soft stall at angles of attack between 30 and 50 degrees. Consistent with finite wing theory, the wing with higher AR have higher lift slopes, attributable to the diminished influence of the wing tip vortex for the higher AR wings. However, the lower AR wings exhibits better lift performance at larger, post-stall values of α, again due to the now beneficial effects of the tip flow. This is particularly evident for the AR = 0.92 wing, where CL reaches a plateau near α = 20◦ which remains almost constant until α ≈ 50◦ . The drag coefficients for all three models are nearly identical at low angles of attack 49 (α < 20◦ ) resulting in increased range and power efficiencies for the higher AR wings. The drag for the higher AR wings increases more sharply for high angles of attack than the lower AR wings, once again due to the diminished role of the tip vortex flow, which may help to postpone separation for LAR wings to higher angles of attack. The wing with the lowest aspect ratio (AR = 0.92) exhibits the best high-α efficiency (CL /CD ) of all the configurations, an observation that may have implications for mammalian gliders which have aspect ratios very close to 1.0, and thus may be able to acheive higher lift and more efficient performance at the very high angles of attack they tend to favor during flight[7, 52]. 3.4.3 Dynamic membrane motions The compliant membrane wing is unique because of its ability to rapidly adapt in response to changing flow conditions. This also means that the membrane is sensitive to slight variations in the flow about the wing, such as periodic vortex shedding, perturbations in the pressure distribution, etc. We imagine therefore, that the membrane wing is in constant motion, and in some regimes this motion might be quite substantial. Since the leading edge of the wing is aerodynamically sharp, one assumes that it is sufficient to induce the formation of leading edge vortices at sufficiently high angles of attack, and one might expect that these shed vortices induce membrane motions. The focus of this section of the discussion is on characterizing and understanding these membrane motions. High-speed stereo photogrammetry was used to measure the unsteady motion of the membrane wings. We found that, for every testing condition, some membrane vibration was present. By computing the time-averaged root-mean-squared (rms) amplitude of the motion over the entire membrane, a spatial structure for these vibrations was obtained over the full range of operating angles of attack and wind speeds. A typical vibrational mode shape 50 0.8 0.6 Span, y/b 0.4 0.2 0.2 0.4 0.6 0.8 Chord, x/c Figure 3.10: Spatial pattern of the unsteady vibration of the compliant wing (α = 28◦ , Re = 1.2 × 105 ). The membrane undergoes unsteady vibrations with a dominant frequency (in this case 110 Hz), and with a spatial structure related to a drumhead mode, characteristic of the driven vibrations of a membrane under tension. 130 5 120 3 Reynolds number, 10 110 4 100 3 90 2 80 70 1 60 0 4 8 12 16 20 24 28 32 Angle of Attack, degrees Figure 3.11: Phase map of membrane drumhead modes (eigenmodes). The number of modes observed in the membrane spatial vibration is shown as a function of speed and angle of attack, with a region close to stall exhibiting the most complex dependancies, presumably due to active vortex shedding from the leading edge. 51 is shown in figure 3.10, in which light and dark regions indicate high and low vibrational amplitudes, respectively. The spatial structures are observed to adopt a standing wave with well-defined nodes and anti-nodes. These structures resemble the eigenmodes that derive from a one-dimensional membrane under tension, subject to Dirichlet boundary conditions at the leading and trailing edges and Neumann boundary conditions at the tip and root. There are, of course, deviations from the ideal membrane mode shapes, but the analogies are quite apparent. This general eigenmode structure was visible in all cases tested, although the harmonic order of the eigenmode (i.e. the number of peaks in the spatial structure) was found to depend quite strongly on the operating conditions. We have constructed a phase map of the harmonic order, providing a qualitative picture of the vibrational modes of the membrane (figure 3.11). At lower speeds the membrane motion is always in the lowest mode - a “breathing” mode in which the entire membrane vibrates in unison. However, higher-order modes and larger amplitude vibrations were excited beyond a critical Reynolds number of approximately 80,000, and in a range of angles of attack falling between 16 − 28◦ . The confinement of the higher eigenvalue modes to this range of Reynolds numbers and angles of attack suggests that the membrane is responding to a resonant forcing by the flow at these conditions. In particular, the angles of attack coincide with the onset of stall, when the flow is beginning to separate and one assumes that there is active vortex shedding on the upper wing surface. In contrast to the results of Galvao et al.[53], where the rms magnitudes were nearly uniform in the spanwise direction, the higher magnitudes of the rms values, i.e. the stronger vibrations, tended to be near the tip. This difference is believed to be a result of the more than 50% increase in the aspect ratio compared to the wing tested by Galvao, resulting in a decrease of the proportion of the wing influenced by the tip vortices. At near-stall angles of attack, the 52 tip flutter induced by flow leakage from the lower to upper surface may be an additional source of vibrational forcing, while the leading edge vortex shedding may act to dampen these vibrations, resulting in the lower rms magnitudes in the inboard region. To gain some insight on the interaction between the flow structures and the natural frequencies of the membrane, we can model the wing as a planar membrane governed by the wave equation: ∂2z a∇2 z = , (3.14) ∂t2 where a is the membrane wave speed. Since the membrane tip and root are unconstrained, the chordwise tension is much greater than the spanwise tension, allowing us to consider only one-dimensional solutions, z(x). The most energetic vibrational response of the system occurs when it is driven at a frequency close to membrane’s natural frequency, f0 , and its harmonics, fi = if0 . Following classical membrane theory [54], the natural frequency, f0 , is dependent on the wave speed, a: a f0 = , (3.15) 2c which in turn is proportional to the square-root of the membrane tension, T : T To + EA a2 = = , (3.16) ρ ρ where ρ is the mass per unit chordwise length for a membrane with cross-sectional area, A. We can estimate these natural frequencies by extending our membrane model described earlier. Assuming once again a parabolic membrane deflection, equations 3.2-3.5 can be 53 used to show that the membrane tension varies linearly with the Weber number: WeEt T = , (3.17) 8ρzmax and that the nondimensional frequency, k0 , is related to the Weber number and camber by r f0 c 1 We k0 = ∼ . (3.18) U U zmax We have already observed (figure 3.5) that the camber, zmax increases at some rate that is slower than linear growth, zmax ∼ Weη , (η < 1), in which case we can predict that the natural frequency of the membrane should decrease with increasing Weber number: r r 1 We 1 We k0 ∼ ∼ ∼ U −η/2 . (3.19) U zmax U Weη This behavior is indeed confirmed by our measurements, shown in figure 3.12. Here we show frequency spectra of the membrane motion as a function of the Weber number (for a fixed angle of attack, α = 24◦ ). The spectra were computed from the chord-normal displacement of the membrane, sampled at the point of maximum vibrational amplitude, and are plotted such that dark shades represent frequencies with high vibrational energy. The solid lines are the predicted membrane frequencies, computed using the measured camber, zmax , and equation (3.19). At low Weber numbers, the membrane tends to be excited at the second eigenmode, 2f0 , corresponding to a Strouhal number of approximately 0.8. This selection mechanism might be due to a resonance between the natural frequencies of the membrane and of vortex shedding from the leading edge. As each vortex is released, it imparts a pressure pulse on 54 3 −35 2.5 Reduced frequency, fc/U 2 −40 5f0 1.5 dB 4f0 1 3f0 −45 2f0 0.5 f0 0 −50 0 0.02 0.04 0.06 0.08 Weber number, We Figure 3.12: Unsteady membrane vibrations. Waterfall chart illustrating the membrane vibration response as a function of wind speed at a fixed angle of attack (α = 24◦ ). For each value of the Weber number, the ensemble-averaged power spectral density is estimated from the Fourier transformation of the unsteady time signal of the individual membrane markers. The power densities are plotted as an intensity map, with the higher power densities denoted by darker regions. The frequency is normalized by freestream velocity and chord: k0 = f c/U . These power spectra are plotted versus the Weber number, defined as We = 12 ρU 2 cCL /(Et). Superposed on the spectra are lines of the natural frequency, f0 , and the higher harmonics (integral p multiples f0 ), of the membrane predicted using membrane theory in which f0 = T /ρ/(2c), where T and ρ are the membrane tension and the mass per unit chord length respectively. The tension is estimated from the observed camber of the wing at each point tested. the upper surface of the membrane, and if the shedding and vibrational frequencies are close, one would expect the membrane to respond as a forced oscillator. However, as the Weber number increases (i.e. the velocity increases), we see an interesting and unexpected result: the dominant response of the membrane jumps (in an almost quantized manner) to higher and higher frequencies, each time closely matching the frequency predicted by our analysis, but associated with higher membrane eigenmodes. The spatial response of the membrane at each speed (figure 3.11) closely matches this trend. The reason for this migration to higher frequencies remains somewhat of a puzzle. We see two possible explanations. The is that the forcing frequency due to vortex shedding (from the leading edge or the laminar separation bubble) remains constant. As the non- dimensional natural frequency of the membrane falls, as predicted by the theory, the mem- 55 brane responds to the forcing through some nonlinear resonance mechanism such that the membrane response is at the closest integer multiple of the forcing frequency present in the flow. While this explanation may be valid, a second, more likely possibility is that the forcing frequency also increases with speed (or Weber number). In many flows, one expects that the Strouhal number for the vortex shedding from a bluff body is nearly constant over a large range of Reynolds (or Weber) numbers. However, this result is only strictly valid for high Reynolds number flows, and recent experiments by Yerusevych et al.[55] of symmetric airfoil wakes have shown that the Strouhal number of the wake vortex shedding can increase significantly at low Reynolds numbers, attributable to the presence of a separation bubble which can exert a significant influence on the subsequent shear layer instabilities. This effect could well be occurring in the present case, but without the benefit of supporting velocity field measurements, it is impossible to fully explain these observations. 3.5 Concluding remarks and implications for mammalian flight The results obtained in this series of experiments include some subtleties that are quite diverse and there are several features of the flow over compliant wings that are quite sur- prising and unexpected. As compared to a rigid wing, the compliant membrane wings generate more lift with higher lift slopes because the camber of the wing increases dynam- ically with increasing aerodynamic loading. In addition, the most compliant membrane model also exhibits a less-severe decrease in lift at the onset of stall than the rigid and less compliant wing due to the de-cambering of the wing as the flow starts to separate. Lastly, stall is observed to be delayed in the compliant wings by as much as 10 degrees. As has been pointed out before [45], to understand the effects of compliance, a key nondimensional parameter must be added to the description of the flow. The Weber number effectively de- 56 scribes the interaction between aerodynamic loading (either by dynamic pressure or angle of incidence) and the membrane elasticity, and this must be considered when performing scale tests or making comparisons with animal flight. Although the compliant membrane wings have the advantage of increased lift, they also exhibit increased drag, although the overall effect on lift to drag is still generally improved. Indeed, the membranes which had been pre-strained (0 > 0) exhibited better aerodynamic efficiencies than the model with no pre-strain. Despite these complexities, two simple analytical models predicts the static and dynamic deformations with surprising accuracy. These models can certainly be improved. In particular the static model would benefit from the use of a higher-order model (for ex- ample thin airfoil and lifting line theory) to predict the aerodynamic loading as a function of angle of incidence. The dynamic model does a good job of predicting the membrane frequencies, but still needs a means to predict the frequency of fluid forcing due to leading edge vortex shedding. Bat flight, like most animal flight, is characterized by low-Reynolds number flight regimes, and we expect that the effects of Reynolds number will be critical to understanding the aerodynamic behavior of these membrane wings. Clearly, the details of the transition to turbulence, the factors that influence leading edge separation and turbulent re-attachment, are strongly influenced by the low-Reynolds number flight environment although a detailed study of these issues was not the focus of this study. However, for these experiments, which were conducted at a variety of wind speeds using a model of fixed dimension, we varied both the Reynolds number and the Weber number by changing the freestream velocity, and thus it can be difficult to tease these effects apart. However, over the narrow range of speeds tested, the primary effects observed are clearly due to Weber number effects, rather than Reynolds number effects. This assertion is supported by the generally excellent 57 collapse of the wing camber data (figure 3.5) over the range of We tested. Having said this, low Reynolds number effects are clearly critical and deserve close attention in future experiments. The motivation for this study originated in a desire to better understand the funda- mentals of the aerodynamics of gliding mammals, such as flying squirrels or sugar gliders. Recent studies of the mechanics of living mammalian gliders [52] provide an opportunity to evaluate the extent to which these physical models capture important features observed in gliding mammals. In those experiments, flying squirrels were observed to glide at average angles of attack of α = 42.5◦ , while sugar gliders averaged α = 44.2◦ . At these angles of attack, the flying squirrels and sugar gliders produced lift coefficients averaging 2.12 and 1.48, respectively, with drag coefficients averaging 0.98 and 1.08. This corresponds to range 3/2 efficiencies (CL /CD ) of 2.26 and 1.39, and power efficiencies (CL /CD ) of 3.33 and 1.70 for the squirrels and sugar gliders, respectively. The lowest aspect ratio model in the present study (AR = 0.92) is the most similar geometry to these species, and at similar angles of attack, the physical models demonstrated a lift coefficient of 0.95, drag coefficient of 0.81, range efficiency of 1.17 and a power efficiency of 1.13. In Bishop’s study, the sugar gliders exhibited aerodynamic performance similar to this AR = 0.92 compliant membrane wing. In contrast, flying squirrels appear to generate substantially more lift with only a slight increase in drag (which translates into higher range and power efficiencies). Flying squirrels have morphological specializations which allows the use of the thumb to control the angle of the leading edge of the wing. Sugar gliders do not have this morphological feature, and their wings more closely resemble the simplified physical models tested in this study. We propose that this anatomical feature may be operated in a manner similar to a leading edge flap on human-engineered vehicles. This articulated leading edge flap allows the flying 58 squirrels to maintain attached flow and consequently to achieve high values of CL at high incidence angles. It should be noted that bats also possess this articulated leading edge flap feature. Mammalian gliders typically glide at relatively high values of the glide angle, i.e. the angle of descent with respect to the horizontal. In addition, these animals glide at high angles of attack which in combination with large values of glide angle put the animals’ body angles near horizontal, i.e. perpendicular to gravity. In this flight configuration, the lift will contribute to the animal’s thrust in addition to counteracting gravity. Similarly, the drag contributes to weight support. Therefore, increased drag due to the higher camber of the compliant wings is not as undesirable as one might first think. In the experiments presented in this paper, the highest lift-to-drag ratio was observed for models for which the camber remains relatively small (below 14% of the chord). Similarly, Bishop[52] found that wing camber of flying squirrels and sugar gliders averaged 10% of their chord length, with the camber exceeding 14% in only a few of the flights of the flying squirrels. Finally, the compliant wings permit relatively efficient flight for a wide range of extreme angles of attack (α > 30). This kind of aerodynamic performance would be desirable in natural flight, as well as for micro air vehicles, where such aerodynamic performance is nec- essary to ensure robust flight characteristics over a wide range of flight conditions, including gusts and variations due to sudden maneuvers to avoid obstacles and other hazards. These experiments represent a small window into the subtle and complex performance of flexible aerodynamic structures and certainly much more needs to be investigated. Chapter 4 Experimental study of vortex-induced forces on a pitching flat plate Arnold J. Song and Kenneth S. Breuer 4.1 Introduction In 1973, Torkel Weis-Fogh described a clap and fling mechanism for lift generation that helped to explain the discrepancies between the large lift coefficients needed for weight support in the hovering flight of Encarsia formosa (a small parasitic wasp) and the lift coefficients he estimated using conventional aerodynamics [56]. To satisfy the Kelvin cir- culation theorem requirement that circulation is neither created nor destroyed, the bound vortex of an impulsively started airfoil is balanced by an equal strength but oppositely signed vortex formed in the wake. Initially, this starting vortex is located close to the airfoil and until there is sufficient distance between the airfoil and the starting vortex, the lift values will be smaller than the values estimated from steady aerodynamics. This unsteady effect, known as the Wagner effect [57], presented a problem when Weis-Fogh tried to es- timate the aerodynamic forces produced during an insect’s wingbeat using a conventional, two-dimensional analysis; there simply was not enough time for lift to build to a level suffi- 59 60 cient for weight support. Using high-speed imaging of the wing motion of Encarsia formosa, Weis-Fogh observed something in the wings’ motion that could result in an unsteady lift mechanism that is not subject to the time lag associated with the Wagner effect. According to his proposed mechanism, the two wings “clap” together at the beginning of each wing beat cycle and, as the they pull apart, fluid rushes into the space between wings generating circulation around each wing. The bound vortices around each wing have circulations of opposite senses and act as starting vortices for the other, so vorticity need not be shed into the wake to satisfy the Kelvin circulation theorem. A striking aspect of this lift mechanism is that it is valid for inviscid flow. Lighthill then extended Weis-Fogh’s theory to account for the viscous effects that are significant for low Reynolds number flyers [58]. He described how these viscous effects would induce flow separation and the formation of a vortex at the wing’s leading edge, which if remained attached to the wing would enhance lift production by increasing the circulation of the wing’s bound vortex. The discovery of this vortex-based mechanism capable of producing high values of lift and thrust in flapping flight has spawned a number of animal flight studies that have confirmed the presence of these leading edge vortices (LEV) in the flight of insect, birds and bats [30, 15, 59, 26, 25, 14]. These vortices can enhance lift but are transitory by nature and may only stay near the wing’s surface for short periods of time before they pinch off and are advected away [11, 60]. For LEVs to be a lift mechanism of any import for flapping wing flight, these flow structures need to remain attached to the wing’s surface to provide lift enhancement for time scales relevant to the wingbeat cycle. Several studies have examined different types of wing motion to further understand the mechanisms that influence the formation and evolution of LEVs: translation (linear motion in a horizontal plane), flapping (rotational motion in a vertical plane) and sweeping (rotational motion 61 in a horizontal plane) [16, 18, 59, 19, 13, 61]. Notably, Lentink and Dickinson [19] show that wings with sufficiently low radii of gyration (in a horizontal plate), as compared to the mean chord length, can maintain stable LEVs and hypothesize that the main stabilizing mechanism centripetal acceleration of the fluid transporting circulation out of the inboard portion of the LEV. Rather than advecting vorticity from the LEV, we propose a growth and stabilization mechanism where the pitching motion of the flat plate modulates the rate of circulation accumulating in the wake. The aim of this study is to show how we can influence the strength and stability of a LEV on a flat plate with an oscillatory pitching motion. 4.1.1 Leading and trailing edge vortex formation and shedding The vortices that form at the leading and trailing edges of a flat plate have a specific time scale for formation and subsequent shedding that depends on the interaction between the shear layers of opposite sign originating from the leading and trailing edge separation points. Circulation from one shear layer will accumulate into a vortex that forms in the wake region close to the body. This vortex will continue to grow until reaching a critical size and strength that is large and strong enough to entrain the opposite signed circulation from the other shear layer. Enough of this oppositely signed circulation will cause the vortex to “pinch-off” and shed from the body [62, 63]. This process will repeat with a vortex forming and shedding from the trailing edge and will continue in alternating fashion between the two shear layers. The Strouhal number, St, characterizes this vortex shedding time by scaling the vortex shedding frequency, fs , by the the free stream velocity divided by the shear layer separation 62 distance, U/(c sin α): fs c sin α fs c˜ ∆tc St = = = . (4.1) U U ∆ts The Strouhal number can also be thought of as the ratio between the time that it takes the fluid to travel the shear layer separation distance, ∆tc , and the vortex shedding period, ∆ts . The shear layer distance, c˜ = c sin α, is the flow-normal distance between the trailing and leading edges. For example, c˜ = 0 when the plate is oriented parallel to the flow and the trailing and leading edges are aligned with each other (α = 0◦ ) and c˜ = c when the plate is perpendicular to the flow (α = 90◦ ). Fage and Johansen [64] observed that the characteristic vortex shedding frequency for a rigid, two-dimensional flat plate was virtually constant with a value of St = 0.15 for angles of attack ranging from α = 30 − 90◦ . For smaller angles of attack (α = 12 − 30◦ ), as α decreased, the Strouhal increased to a maximum of St = 0.22, similar to the Strouhal number observed for vortex shedding behind a cylinder. The strength of these vortices and their shedding rate depend on the circulation flux into the shear layers from the leading and trailing edges. Fage and Johansen [64] estimated that the rate of circulation shed into the wake, dΓ/dt, from the leading and trailing edge separation points was nearly constant and was proportional to the dynamic pressure at the edge of the boundary layer near the separation point: 1 dΓ/dt = Ub2 (4.2) 2 where Ub is the velocity at the edge of the boundary layer at the separation point. To maximize the vortex lift that results from the formation of an LEV, we would like to capture as much circulation in the leading edge vortex before it is shed downstream. Bearman [] combines equations 4.1 and 4.2 to estimate the maximum strength of a leading edge vortex, 63 Γv ,       dΓ 1 2 c˜ 1 2 1 Γv = β ∆ts = β U =β U (4.3) dt 2 b U St 2 b fs where β is the fraction of the circulation that survives in a wake vortex, c˜ is shear layer separation distance, St is the Strouhal number and fs is the vortex shedding frequency. The circulation survival coefficient, β, for a flat plate is approximately equal to 0.6 and is typically less than unity due to the mixing of oppositely signed circulation. We assume that one way to increase vortex lift is to increase the amount of circulation captured in a leading edge vortex. Using Eq. 4.3, we see that we can increase the leading edge vortex strength, Γv , by either increasing the separation point velocity, Ub , or reducing the shedding frequency, fs . But, we know that the shedding frequency for a static flat plate for relatively high angles of attack, α = 30 − 90◦ , is nearly constant and for lower angles of attack the Strouhal number increases, which is the opposite effect desired. For a static flat plate, our only option to increase the vortex strength would be to increase Ub , e.g., using blowing control at the leading edge. We propose that oscillatory pitching of the flat plate may induce vortex shedding be- havior that deviates from the Strouhal scaling (i.e., a shedding frequency determined by the nearly constant Strouhal number observed for a static flat plate). We define the reduced velocity, Ur , a parameter that describes the unsteady pitching of the plate by normalizing the oscillation period, 1/f0 , by the time that it takes a parcel of fluid to travel one chord length, U/c, as follows: U Ur = (4.4) f0 c where f0 is the plate’s oscillation frequency. The pitching motion modifies the circulation shedding rate, dΓ/dt, by increasing or decreasing the flow velocity at the separation points, 64 Ub , depending on the pitch direction. The flow velocity at the leading edge will increase for a flat plate pitching in a decreasing α direction, which increases dΓ/dt, and vice versa when the flat plate pitching in an increasing α direction. If the pitching amplitude is large enough and the flow conditions are right, the vortex shedding and plate motion coordinate at a single oscillation frequency, which is referred to as locked-in resonance between the fluid and structure [65]. Vortex-induced motion of bluff bodies and airfoils has been examined theoretically and experimentally by numerous investigators [66, 67, 68, 69, 63, 70, 71, 65, 72, 73]. We focus on the studies of Theordorsen and Garrick [68] and von Karman and Sears [69], which both use linear thin airfoil theory to characterize the aerodynamic force production resulting from sinsoidal oscillations of a thin flat plate under the assumption that the pitching motion was composed of small deviations from zero angle of attack with a wake shed parallel to the freestream. These linear analyses showed that for a plate with only a pitching degree of freedom, any free oscillations of the airfoil would be damped out by the fluid for all oscillation frequencies as a result of the phase lag between the pitching moment and the plate’s motion. It was then concluded that the initiation of self-sustaining airfoil motion (the locked-in state) requires nonlinear mechanisms, such as flow separation, to be in play [67]. In this study, we experimentally examined the time scale of vortex shedding and the momentum exchange from the fluid to the plate for passive and active pitching of a flat plate due to the nonlinear mechanism of flow separation and wake vortex dynamics. For the passive pitching experiments, the plate was mounted to an elastic support and the pitching dynamics resulted from a matching between the fluid and elastic forces. We refer to this type of fluid-structure interaction as the passive pitching of the flat plate because the motion 65 depends on the fluid forcing. For the elastically mounted flat plate, the locked-in oscillation frequency and amplitude are prescribed by the resonant frequency of the torsion spring and the plate’s moment of inertia. We varied the flow conditions while measuring the amplitudes of the plate vibrations and wake unsteadiness to infer the strength of the momentum transfer between the fluid and the plate since sustained oscillatory motion requires energy input from the fluid to plate to overcome mechanical damping. But in the experiments with the elastically mounted plate, the oscillation frequency and amplitude were determined by the mechanical properties (torsion spring stiffness and moment of inertia) of our particular model and will not generally be the frequency and am- plitude for which maximum momentum transfer occurs. Therefore, we conducted another set of experiments where we prescribed the plate motion (active pitching that is indepen- dent of elastic properties) using a rotary servo motor to examine the relationship between the pitching motion and vortex shedding for an expanded parameter space of pitching am- plitudes and frequencies. We follow this introductory section with a description of the elastically mounted and servo-controlled experiments, then a report of the experimental results and conclude with a discussion of the relationship between flat plate pitching motion and the strength of shed vortices. 4.2 Materials and methods 4.2.1 Elastically mounted flat plate The vortex induced vibrations experiments were conducted in the wind tunnel facility at Brown University. This wind tunnel, manufactured by Engineering Laboratory Design, 66 Inc. (Lake City, MN) is a closed return configuration with a rectangular cross-section with dimensions, 61×61×122 cm (width×height×length), and can generate wind speeds up to 55 m/s. The fan is powered by a constant speed motor, which is housed within the wind tunnel, and the wind speed is controlled by varying the pitch angle of fan blades. A water cooled heat exchanger upstream of the settling chamber maintains an air temperature near 20◦ C by removing the heat produced by the fan motor and from viscous heating. The model tested was an aluminum flat plate with dimensions: 10 × 10 × 0.16 cm (chord × span × thickness) with sharp leading and trailing edges, which help facilitate comparison between the wind tunnel experiments and numerical simulations by pinning the separation points to the leading and trailing edges. The mass and moment of inertia (MOI) of the flat plate were 14.7 grams and 4.7 kg-m, respectively. For all the experiments presented, the moment of inertia (MOI) and mass of the flat plate were held constant while varying length of the torsion rod and plate aspect ratio to keep the resonant frequencies constant for this study. The MOI and mass of the flat plate were held constant by removing material from the interior region of the plate and covering the hole with a thin Mylar film that added negligible mass to the model. The film was inextensible and impermeable to maintain the plate’s aerodynamic integrity. The plate was attached to a slender 3.2 mm diameter aluminum rod with its axis aligned with the plate centerline with the other end of the rod was clamped to the test section sidewall, as seen in Fig. 4.2. In this configuration, the flat plate model was subject to heaving, or linear translation of the plate, in the direction perpendicular to the chord, and pitching, with the resonant frequencies of these vibrational modes being controlled by the rod length between the clamp and the flat plate, as shown in Fig. 4.1. The slender rod has a circular cross section and has two modes of deformation: twisting 67 Figure 4.1: Resonant frequency dependence on rod length. The slender rod has a circular cross section and has two modes of deformation: twisting and bending. In this figure, we show the resonant pitching and heaving frequencies of the flat plate measured from free vibration tests. and bending. For the twisting mode, the value of the torsional spring constant, kt is: GJ kt = , (4.5) l where G, J, and l, are the shear modulus (for aluminum, G = 26 GPa), polar area moment π 4 of area (for a cylinder, J = 2 r ), and the rod length, respectively. The resonant torsional frequency, f0 , is: 1 p f0 = kt /I (4.6) 2π where I is the moment of inertia of the flat plate. The flat plate’s vortex induced motion of the model was captured using a high frame rate camera (Allied Vision Technologies, Stadtroda, Germany) placed with an end-on view of the plate wingtip. Two retroreflective tape markers were placed at the leading and trailing edges of the wingtip, as seen in Fig. 4.2, and custom software utilizing the OpenCV computer vision library tracked the marker locations in realtime at a rate of 200 Hz. The retroreflective markers contrast sharply with the rest of the model and darkened background 68 Flat plate Tracking marker Elastic rod High-speed video camera Tracking marker (a) Model schematic (b) Tracking image Figure 4.2: Passive flapper model diagram. The flat plate model consists of 1.6 mm thick aluminum flat plate with square edges attached at the center to a slender 3.2 mm diameter rod. The other end of the rod is clamped and the natural pitching frequency of the model can be adjusted by either lengthening or shortening the rod length. The motion of the model is captured by a high-speed camera placed to have an end-on view of the plate edge and the two markers place at the leading and trailing edges. The marker positions were determined in real-time at a rate of 200 Hz. 42 α0=30° 0.4 40 α0=38° 0.3 Stream−normal, y/c 0.2 38 Pitch angle [deg] 0.1 36 0 34 −0.1 −0.2 32 −0.3 30 −0.4 28 −0.4 −0.2 0 0.2 0.4 0.6 0 0.05 0.1 0.15 0.2 Streamwise, x/c Time [s] (a) Tracking results (b) Pitching motion Figure 4.3: Realtime video based motion tracking at 200 Hz. The model was a rigid flat plate, therefore a simple line fit to the two end marker locations provides information about the location and orientation of the plate in response to aerodynamic forcing. In (a), we show the end-plate profiles of two mean angles of attack, α0 = 30 and 38 degrees, at a chord-based Reynolds number of 115 × 103 . As seen in (b), the pitching amplitude increases by more than an order of magnitude. 69 making a simple threshold filter an effective method to create a binary image that extracts the markers from the original video frame. The OpenCV function blobs then searches for contours, which are the borderlines between the true and false binary pixels, within the image and returns a list of blobs, or closed contours, whose properties can be queried. For this series of experiments, we queried the centroid location for each blob and recorded these coordinates for each frame. The plate’s angular position was calculated from the slope of   −yT E the line connecting the leading and trailing edge markers, where α = tan−1 xyLELE −xT E , and can be plotted in time as in Fig. 4.3. The wake velocity was measured by a calibrated hot-wire anemometer placed at a loca- tion 50 cm (5 chord lengths) downstream. The hotwire was positioned laterally at midspan at a height equal to the leading edge height. 4.2.2 Servo-controlled flat plate The experiments with an elastically mounted flat plate demonstrated the conditions where strong coupling between the elastic vibrations and vortex shedding occurred. We conducted a companion set of experiments with a servo-controlled flat plate for which the pitching frequency and amplitude were prescribed to determine the conditions of optimum coupling between the flat plate and vortex shedding independent of the mechanical parameters that govern an elastic plate’s dynamics. The servo-controlled model was also a flat plate with chord length, chord = 10 cm, but with a higher aspect ratio of 4.3. The plate’s axis of rotation ran spanwise at the half-chord location and the pitching motion was driven by a brushless servo motor (SM233AE, Parker Compumotor, Rohnert Park, CA) with a torque constant value of Kt = 0.415N-m/A and a peak output torque of 3.4 N-m. The motor and model were coupled with a reaction torque sensor (TFF400, 2.8 N-m maximum torque, 70 Tunnel wall Flat plate model Tunnel wall Angular Servo motor Torque sensor position sensor Shaft coupling (a) Model setup Function generator Torque sensor, τ αfunction αerror PID Servo +- controller motor Angular position sensor, α (b) Position following control Figure 4.4: Servo controlled pitching model. We conducted experiments where the motion of a flat plate was prescribed and actuated by a servo motor. The model was attached to a shaft that ran from an angular position sensor on one end, through the model, to a torque sensor on the other end. The driving servo motor was coupled to the model shaft with a torque sensor, which measured the torque applied by the servo motor to the model. The servo was operated in a position following mode where the angular position of the plate was prescribed by a sinusoidal control signal of specified frequency and amplitude. Futek, Irvine, CA), which measured the torque applied by the motor to the flat plate. The angular position was measured at the end opposite of the motor with an angular position sensor (P500, Positek, Gloucestershire, UK). The servo motor control loop was closed with a position-following proportional-integral-derivative controller using the voltage output of a function generator as the control signal, as shown in Fig. 4.4. 4.3 Results and discussion 4.3.1 Passive pitching In this series of experiments, we measured the the vortex-induced motion of two elastically mounted flat plates with low-aspect-ratios, AR = 1 and 2. The flat plates were inclined for 71 a range of angles attack (AR = 1: α = 20 − 70◦ and AR = 2: α = 0 − 60◦ ) for wind speeds ranging from U = 10.3 to 20 m/s (Re = 6.7 × 104 − 1.3 × 105 ). The motion of the plate was recorded using high frame rate video tracking and the vortex shedding behavior was characterized by hot-wire anemometer measurements in the wake. The sharp leading edge of our flat plate models encourages flow separation leading to the formation and shedding of the vortices that transfer energy from the fluid to the model through variations in the pressure field. When circulation captured in each vortex is small, the pressures felt by the flat plate will be small compared to the elastic restoring forces (e.g., for low speed flows, low angles of attack or very large torsional stiffness). For these subcritical cases, where locked-in coupling of the plate and vortex shedding has not yet occurred, the pitching amplitudes are very small (∆α < 1◦ ). Therefore, we can treat the coupling between the vortex shedding and plate as being applied in only one way from the fluid to the structure where the fluid motion is independent of the body motion. We will refer to this one way coupling as the subcritical case. As the circulation captured by each shed vortex increases (increasing the amplitude of the pressure fluctuations felt by the flat plate), the plate motion will increase until transitioning to a state where sustained large amplitude motion of the flat plate results from the coordination between the plate’s pitching motion and the pressure fluctuations on the plate. In this supercritical state, the mutual feedback between the fluid and the structure need not be a state where the resonant structural vibrational frequency and the “natural” frequency of the fluid (or the static Strouhal frequency) are matched. In fact, our experiments show when the plate motion is sufficiently large (we observed pitching amplitudes approaching ∆α = 7◦ ), the vortex shedding frequency is driven by the resonant frequency of the torsion spring, which can be quite different than the equivalent static Strouhal shedding frequency. 72 4.3.2 Wake velocity measurements We utilize waterfall plots in Figure 4.5 to visualize the power spectra of the velocity fluctu- ations as a function of reduced velocity, Ur = U/(f0 c), where f0 and c denote the resonant frequency of the torsion spring and the plate chord length, respectively. Each power spec- trum, which is calculated from unsteady velocity measurements from a hotwire anemometer placed 5 chord lengths downstream, is plotted as intensity values with the white and black coloring indicating high and low values of velocity unsteadiness and the dominant frequency for each reduced velocity value is distinguished by a red marker (location of maximum power) and a red curve provides reference for the location of the resonant frequency of the flat plate, f0 . For the AR = 1 plate, the coupling was subcritical (i.e., one way coupling from the fluid to the plate) up to relatively high angles of attack (up to α = 35◦ ), where the dominant vortex shedding frequency followed a scaling consistent with a constant Strouhal number value, St = 0.08 (Fig. 4.5(a)). These constant Strouhal number shedding frequencies were lower than the resonant frequency of the flat plate for the entire range of velocities tested for the subcritical values of angle of attack (α = 30◦ and 35◦ ). Not only were the dominant frequencies of the vortex shedding mismatched to the structural resonant frequency, the power contained at these peak frequencies was relatively weak and was not large enough to excite any appreciable plate motion. Even when the vortex shedding frequencies approached the resonant structural frequency near Ur = 5, the strength of the velocity fluctuations actually decreased slightly suggesting that a transition to supercritical coupling would not occur even if the dominant shedding frequency and resonant structural frequencies were matched at these angles of attack (α < 35◦ ). As we increased the angle of attack to α = 38◦ (Fig. 4.5(b)), the vortex shedding 73 (a) Subcritical behavior, α = 30 and 35◦ (b) Transition to supercritical behavior, α = 38◦ Figure 4.5: Wake velocity spectra. We examine the nature of the wake unsteadiness using this waterfall plot that visualizes the power spectra of the wake velocity fluctuations as measured by a hotwire anemometer placed 5 chord lengths downstream. Each power spectrum is plotted vertically at each value of reduced velocity as intensity values with the white and black coloring indicating high and low values of velocity unsteadiness and the dominant frequency for each reduced velocity value is distinguished by a red marker (location of maximum power) and a red curve provides reference for the location of the resonant frequency of the flat plate, f0 . The plot underneath the waterfall plot indicates the power magnitude of the dominant wake velocity frequency. 74 behavior transitioned from having a well-defined characteristic shedding frequency that followed a constant Strouhal number scaling to shedding that contained a broader band of higher frequencies without a clearly defined dominant frequency. The broad band frequency content of the wake velocity fluctuations suggests that the bluff body vortex shedding has transitioned from laminar or transitional vortex shedding consisting of large scale flow structures in the wake to smaller scale structures characteristic of turbulent shedding. We observed a threshold at Ur = 3 where the coupling between the vortex shedding and the plate motion exhibited a transition from a subcritical to a supercritical state made clear with a dramatic change in the vortex shedding behavior. When the coupling shifted from subcritical to supercritical behavior, nearly all of the fluctuation energy was contained in a shedding frequency matched to the resonant structural frequency indicating that the vortices were being shed in coordination with the plate motion at a frequency that deviated from the vortex shedding frequency for the body in stasis (the constant Strouhal number line at St = 0.08). After lock-in occurred for our elastically mounted flat plate, the strength of the velocity fluctuations quickly reaches a maximum at Ur = 3.2 and subsequently decreases with increasing velocity. Aspect ratio We present results that demonstrate how aspect ratio affects the transition between subcrit- ical and supercritical coupling. We measure the flat plate pitching motion using high frame rate video and track markers at the leading and trailing edges to determine the plate’s time- varying angle of attack. In Figure 4.6, the root-mean-squared value of the unsteady flat plate motion is plotted as a function of the initial angle of attack and the reduced velocity. We refer to the angle for the flat plate at a zero stress state when there are no aerodynamic 75 (a) Aspect ratio = 1 (b) Aspect ratio = 2 Figure 4.6: Effect of aspect ratio on pitching intensity. As the aspect ratio of wing increases, the flow become increasing two-dimensional because of the diminished influence of the tip flow. The vibrational power is plotted for two flat plates with aspect ratios, AR=1 and 2. The moment of inertia is identical for the two flat plates to maintain the same natural pitching frequency. For the higher aspect ratio plate (AR=2), high amplitude pitching initiates at much lower angles of attack and is maintained for a wider range of angles of attack. 76 or elastic stresses felt by the plate as the initial angle of attack, α0 . For the AR = 1 case (Fig. 4.6(a)), the threshold reduced velocity was approximately Ur = 3 below which the plate and vortex shedding was subcritical. As the velocity was increased, supercritical coupling was initiated for a small range of angles of attack (∆α ≈ 2◦ for Ur = 3.4). The range of angles of attack where supercritical coupling occurs increased from a ∆α ≈ 10◦ for Ur = 3.6 to a ∆α of more than 20◦ at the highest velocity tested (Ur = 4.6). We then doubled the aspect ratio while keeping the moment of inertia constant and observed an expansion in the range of angles of attack where supercritical coupling occurred (see Fig. 4.6(b)). The threshold reduced velocity remained the same at approximately Ur = 3, but supercritical behavior was initiated for significantly lower angles of attack to values as low as α0 = 6◦ at Ur = 5.7. The upper boundaries of the supercritical regime, however, for the two aspect ratios were nearly the same. We believe that the earlier transition (α0,crit < 35◦ ) to supercritical behavior for the AR = 2 plate was a result of the reduction tip vortex influence for the inboard regions of the higher aspect ratio plate. The low Reynolds number (Re = O(104 − 105 )) aerodynamics of these very low-aspect-ratio wings, with AR ≈ O(1) can be very different from conventional aerodynamics. When the wing’s aspect ratio is small, the tip vortices have influence over a greater proportion of the total span. The downward directed flow induced by these tip vortices has the effect of reducing the LAR wing lift slope from the infinite wing value of 2πα. But, the tip vortex downwash also has the beneficial effect of increasing the LAR wing’s resistance to flow separation by literally pushing the boundary layer flow against the wing’s upper surface. Torres and Mueller [6] showed that for LAR wings the onset of stall is delayed to higher values of α and the severity of the loss of lift at stall is attenuated. For larger angles of attack after stall, the lift values will level off and then gradually de- 77 crease rather than experiencing a sharp drop in lift observed for 2-dimensional wings. As the aspect ratio of a flat plate increases, these tip vortices exert influence over a smaller proportion of the planform area and the tip vortex downwash has a diminished effect on separation suppression. The transition back to a subcritical state as the angle of attack is increased further forms the upper boundary of the supercritical regime. This return to subcritical coupling suggests that the plate’s motion is no longer able to overcome the mismatch between the “natural” vortex shedding and the resonant structural frequencies. Aspect ratio did not have much effect on the shedding frequency once the flow was fully separated, therefore this upper boundary of supercritical behavior was nearly the same for both plate geometries. 4.3.3 Forced oscillations To examine how the vortex shedding behavior varies as a function of pitching frequency and amplitude, freestream velocity and initial angle of attack, we conducted experiments where the pitching motion of a flat plate was prescribed using servomotor actuation and measured the torque applied to the plate that was needed to follow the desired motion. As described previously, we calculated the net aerodynamic work per cycle and used this value as an indicator of effective transfer of momentum from the fluid to the plate. When the net aerodynamic work value is positive (the plate is gaining energy from the fluid), this is equivalent to the limit cycle oscillations of the elastically mounted flat plate in a supercritical state. 78 The servo-controlled motion can be described by: α(t) = ∆α sin(2πf0 t) + α0 (4.7) α(t) ˙ = 2πf0 ∆α cos(2πf0 t) where ∆α is the pitching amplitude, α0 is the mean angle of attack (equivalent to the zero stress condition for a torsion spring) and f0 is the pitching frequency. W = ∮ τ dα α Figure 4.7: Net aerodynamic work per cycle. The amount of work per cycle needed to pitch the plate is the path integration of torque versus angular position, which includes the torque due to the plate inertia and aerodynamic forces. However, we are only interested in the net aerodynamic work applied to the plate by the fluid, which is indicated by hysteresis in the torque/angular position phase space. If the there is no hysteresis, the net work per cycle is equal to zero because the work done on the plate by the fluid during the pitching up phase is the same as the work done on the plate for the pitching down phase. The plate is extracting energy from the fluid when the loop travels clockwise, as opposed the needing to add energy to the system to overcome the fluid damping when the loop travels counterclockwise. The amount of angular momentum transferred from the fluid to the flat plate per cycle, W , can be determined by integrating the torque applied to the plate, τ , versus angular position, α, path for each cycle, as seen in Fig. 4.7, I W = τ dα (4.8) where the net work transferred from the fluid to the plate is positive if the loop is travelled clockwise or negative work for the counterclockwise direction. When the aerodynamic work 79 is positive, the direction of momentum transfer is from the fluid to the plate and this excess energy the plate extracts from the flow must be dissipated in the motor to maintain the simple harmonic motion. Conversely, when the aerodynamic work is negative for the plate’s pitching cycle, the motor must add additional energy to counterbalance the fluid damping. We use this work per cycle metric as a proxy for the vortex strength using the rationale that vortices with large circulation values will induce a larger moment on the plate, τ , therefore transferring a larger value of work, W , to the plate. In Figure 4.8, we overlay the cycle trace of the measured moment coefficient, Cm , versus the angle of attack over the static moment values for several initial angles of attack, α0 = 0◦ , 30◦ , 60◦ and 90◦ . For a very low pitching frequency, or large reduced velocity (Ur = U/(f0 c), e.g. Ur = 200 shown in Fig. 4.8(a), the dynamic moment loops closely follow the static values as long as the pitching motion remained below (α0 = 0◦ ) or above (α0 = 60◦ and 90◦ ) the separation point at approximately α = 37◦ . The α0 = 30◦ case exhibited hysteresis in the moment versus angle of attack behavior where the moment was higher for the pitch down (α decreasing) half of the cycle compared with the moment for the pitch up (α increasing) portion of the cycle. The pitching motion travels the hysteresis loop in the counter-clockwise position indicating that the additional energy was needed to maintain the prescribed pitching motion. When we increased the pitching frequency, or decreasing the reduced velocity, by an order of magnitude to Ur = 20 (Fig. 4.8(b)), the torque versus angular position loop slopes steepened due to the increased angular momentum of the plate and began to open up for a wider range of α0 indicating that the plate motion itself was inducing separation behavior that deviated from that of a static plate. We then calculate the work per pitching cycle as a function of mean angle of attack for open loop control of the pitching motion, as seen in figures/ch4 4.9-4.11. The work produced or consumed by 80 the fluid was measured for three wind speeds, Re = 6.7 × 104 , 8.0 × 104 and 9.3 × 104 , and four pitching frequencies, f0 = 8, 10, 12, 14 Hz corresponding to reduced velocities, Ur = 16.7 to 7.1. As the pitching amplitude increased for each value of the reduced frequency, the peak value of extracted work shifted towards higher values of mean angle of attack from approximately α0 = 10◦ to 18◦ ; this behavior was consistent for the range of reduced frequency and Reynolds number considered. For Re = 6.7 and 8.0 × 104 , the net amount of work transferred from the fluid to the plate was negative (i.e., additional work provided by the servo motor was needed to sustain the plate’s motion). But, for Re = 9.3 × 104 , the plate extracted a net positive amount of work from the fluid for reduced velocities, Ur = 16.7, 14.3 and 11.1 with the peak mean angle of attack ranging from α0 = 14-16◦ . When we hold the the reduced frequency constant at Ur = 10.0, as in Fig. 4.12, we observed for a small pitching amplitude (∆α = 2◦ in Fig. 4.12(a)) that the net work per cycle was negative and relatively flat over the range of mean angles of attack tested with the exception of a small peak near α0 = 8◦ . As the pitching amplitude was increased to ∆α = 4◦ , the net work remained negative with the peak energy shifting slightly to α0 = 10◦ . This peak in energy became more prominent relative to the cycle average energy at angles of attack above and below the peak where there was a difference of approximately 5 mJ between the peak and the lower minimum (α0 < 10◦ ) and a difference of nearly 4 mJ between the peak and the upper minimum (α0 > 10◦ ). This trend continued as the pitching amplitude was increased up to ∆α = 10◦ at which point the peak energy had shifted to α ≈ 15 − 16◦ with a differential of 15 mJ between the peak energy and the lower minimum and a 10 mJ difference between the peak energy and upper minimum. These trends were independent of the range of Reynolds number tested (Re = 6.7 × 104 − 9.3 × 104 ) for the smaller pitching amplitudes up to ∆α = 6◦ . For the larger pitching amplitudes, ∆α = 8◦ and 10◦ , 81 the peak energy values for Re = 9.3 × 104 were nearly 5 mJ larger than the lower wind speeds suggesting that the change in the energy transferred from the fluid to the plate is Reynolds number dependent. We suspect, however, that the apparent Reynolds number dependence of the aerodynamic work is due to a flaw in the apparatus design rather than a true change in the vortex shedding behavior and subsequent strength over the pitching cycle. For this series of experiments, the torque sensor coupled the shaft section connected to the motor with the section containing the model, as illustrated in Fig. 4.4. Misalignment between the shaft sections would result in spurious torque readings. Although great care was taken to minimize misalignment between the two shaft sections, the speed dependence of the aerodynamic work indicates that shaft misalignment occurred and may have been due increased axial loads from higher lift and drag values for increased wind speeds. Subsequent studies have revised the experimental apparatus to eliminate the effect of shaft misalignment on the torque measurements. 82 (a) Reduced velocity, Ur = 200 (b) Reduced velocity Ur = 20 Figure 4.8: Torque versus angular position hysteresis. For large values of reduced frequency , or low pitching frequencies, (Ur = 200, as in (a)), the pitching motion hardly influences the aerodynamic moment. The dynamic torque closely tracks the static torque behavior for α0 = 0, 60, and 90◦ , therefor one would expect the flow around the pitching plate to look similar to the static plate. However, hysteresis is observed for α0 = 30◦ near the separation point of α = 37◦ . As pitching frequency is increased (e.g. as in (b), Ur = 20), the torque versus angular position loops begin to open up indicating aerodynamic work is being done. 83 −3 −3 x 10 x 10 5 5 2° 2° 4° 4° 6° 6° 0 0 8° 8° 10° 10° Work per cycle [J] Work per cycle [J] −5 −5 −10 −10 −15 −15 −20 −20 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Mean angle of attack [degrees] Mean angle of attack [degrees] (a) Ur = 7.1 (b) Ur = 8.3 −3 −3 x 10 x 10 5 5 2° 2° 4° 4° 6° 6° 0 ° 0 8 8° ° 10 10° Work per cycle [J] Work per cycle [J] −5 −5 −10 −10 −15 −15 −20 −20 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Mean angle of attack [degrees] Mean angle of attack [degrees] (c) Ur = 10.0 (d) Ur = 12.5 Figure 4.9: Re = 6.7 × 104 . The work produced or consumed by the fluid was measured for three wind speeds, U = 8, 10, 12 m/s, and four pitching frequencies, f0 = 8, 10, 12, 14 Hz. For each wind speed and pitching frequency, the mean angle of attack of the maximum work extracted from the fluid shifts with increasing pitching amplitude, from 2◦ to 8◦ from approximately α = 10◦ to α = 18◦ . These peaks shift very slightly to lower mean angles of attack as the wind speed goes up. In addition, the work extracted from the fluid is only positive for U = 14 m/s at the highest pitching amplitude 10◦ . The magnitude of the peak work extracted stays nearly constant across pitching frequencies until f0 = 14 Hz where the peak extracted work falls below zero. 84 −3 −3 −3 x 10 x 10 x 10 5 5 5 2° 2° 2° ° ° 4 4 4° 6° 6° 6° 0 ° 0 ° 0 8 8 8° 10° 10° 10° Work per cycle [J] Work per cycle [J] Work per cycle [J] −5 −5 −5 −10 −10 −10 −15 −15 −15 −20 −20 −20 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Mean angle of attack [degrees] Mean angle of attack [degrees] Mean angle of attack [degrees] (a) Ur = 8.3 (b) Ur = 10.0 (c) Ur = 12.5 −3 x 10 5 2° 4° 6° 0 8° 10° Work per cycle [J] −5 −10 −15 −20 0 5 10 15 20 25 30 Mean angle of attack [degrees] (d) Ur = 14.3 Figure 4.10: Re = 8.0 × 104 . See caption in Fig. 4.9 for details. −3 −3 −3 x 10 x 10 x 10 5 5 5 2° 2° 2° ° ° 4 4 4° 6° 6° 6° 0 ° 0 ° 0 8 8 8° 10° 10° 10° Work per cycle [J] Work per cycle [J] Work per cycle [J] −5 −5 −5 −10 −10 −10 −15 −15 −15 −20 −20 −20 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Mean angle of attack [degrees] Mean angle of attack [degrees] Mean angle of attack [degrees] (a) Ur = 10.0 (b) Ur = 11.1 (c) Ur = 14.3 −3 x 10 5 2° 4° 6° 0 8° 10° Work per cycle [J] −5 −10 −15 −20 0 5 10 15 20 25 30 Mean angle of attack [degrees] (d) Ur = 16.7 Figure 4.11: Re = 9.3 × 104 . See caption in Fig. 4.9 for details. 85 −3 −3 x 10 x 10 5 5 0 0 Work per cycle [J] Work per cycle [J] −5 −5 −10 −10 −15 −15 Re=6.7×104 Re=6.7×104 Re=8.0×104 Re=8.0×104 4 Re=9.3×10 Re=9.3×104 −20 −20 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Mean angle of attack [degrees] Mean angle of attack [degrees] (a) Ur = 10.0, ∆α = 2◦ (b) Ur = 10.0, ∆α = 4◦ −3 −3 x 10 x 10 5 5 0 0 Work per cycle [J] Work per cycle [J] −5 −5 −10 −10 −15 −15 Re=6.7×104 Re=6.7×104 Re=8.0×104 Re=8.0×104 Re=9.3×104 Re=9.3×104 −20 −20 0 5 10 15 20 25 30 0 5 10 15 20 25 30 Mean angle of attack [degrees] Mean angle of attack [degrees] (c) Ur = 10.0, ∆α = 6◦ (d) Ur = 10.0, ∆α = 8◦ −3 x 10 5 0 Work per cycle [J] −5 −10 −15 Re=6.7×104 Re=8.0×104 Re=9.3×104 −20 0 5 10 15 20 25 30 Mean angle of attack [degrees] (e) Ur = 10.0, ∆α = 10◦ Figure 4.12: Aerodynamic work, Ur = 10.0. 86 −3 −3 x 10 x 10 4 4 Re=6.7×104 Re=6.7×104 4 Re=8.0×10 Re=8.0×104 2 4 2 Re=9.3×10 Re=9.3×104 0 0 Work per cycle [J] Work per cycle [J] −2 −2 −4 −4 −6 −6 −8 −8 −10 −10 6 8 10 12 14 16 18 20 6 8 10 12 14 16 18 20 Reduced velocity, U/(f0c) Reduced velocity, U/(f0c) (a) α0 = 14◦ , ∆α = 2◦ (b) α0 = 14◦ , ∆α = 4◦ −3 −3 x 10 x 10 4 4 Re=6.7×104 Re=6.7×104 Re=8.0×104 Re=8.0×104 2 2 Re=9.3×104 Re=9.3×104 0 0 Work per cycle [J] Work per cycle [J] −2 −2 −4 −4 −6 −6 −8 −8 −10 −10 6 8 10 12 14 16 18 20 6 8 10 12 14 16 18 20 Reduced velocity, U/(f0c) Reduced velocity, U/(f0c) (c) α0 = 14◦ , ∆α = 6◦ (d) α0 = 14◦ , ∆α = 8◦ −3 x 10 4 Re=6.7×104 Re=8.0×104 2 Re=9.3×104 0 Work per cycle [J] −2 −4 −6 −8 −10 6 8 10 12 14 16 18 20 Reduced velocity, U/(f0c) (e) α0 = 14◦ , ∆α = 10◦ Figure 4.13: Transition to supercritical behavior for higher Reynolds number flows, α0 = 14◦ . For low pitching amplitudes (∆α = 2◦ and 4◦ ), the net work values versus were independent of Reynolds number and showed the same trend versus reduced velocity. For larger amplitude pitching (∆α > 6◦ ), the highest Reynolds number case (Re = 9.3 × 104 ) consistently transferred more work from the fluid to the plate and actually transitions to positive work for Ur = 11.6, 13.7 and 16.7. 87 In Figure 4.13, we show how the net work transferred to the plate changed as a function of reduced velocity for several values of Reynolds number (Re = 6.7×104 , 8.0×104 and 9.3× 104 ). We choose to drill down at α0 = 14◦ because we observed the largest variation in work values as a function of Reynolds number at this mean angle of attack. For low pitching amplitudes (∆α = 2◦ and 4◦ ), the net work values versus were independent of Reynolds number and showed the same trend versus reduced velocity. For larger amplitude pitching (∆α > 6◦ ), the highest Reynolds number case (Re = 9.3 × 104 ) consistently transferred more work from the fluid to the plate and actually transitions to positive work for Ur = 11.7, 13.7 and 16.7. This transition from negative to positive work appears to be similar to the transition at the threshold reduced velocity, Ur = 3 for the elastically mounted plate, but we observed this transition at pitching frequencies nearly 4 times lower with the forced pitching experiments near Ur = 10. This result indicates that supercritical coupling can occur for large mismatches between the pitching amplitude and the static vortex shedding frequencies if the flow velocity and pitching amplitude are sufficiently large. 4.4 Concluding remarks We conducted experiments with a flat plate that pitched in response to the aerodynamic moment generated by the plate geometry and angle of attack or forced to oscillate at a predetermined amplitude and frequency. Our aim with this study was to understand how the timing between the vortex shedding and the plate motion might optimized to maximize the momentum transferred from the fluid to the plate. From the passively and actively pitched plate experiments, we show that supercritical coupling requires that the flow suspected that the vortex shedding does not contain enough energy to excite the resonant pitching mode of the plate until it was inclined to α > 35◦ and the flow velocity reached a threshold value 88 Figure 4.14: Explanatory model for positive net aerodynamic work in separated regime. As the plate pitches down (t = π/2 to π), a leading edge vortex (LEV) forms and grows until the plate reaches the pitching cycle minimum angular position. The upward pitching motion after the reversal point (t = π to 3π/2) helps stabilize the LEV by inducing a downward flow, which reduces the effective angle of attack. When the plate reaches the end of the cycle (t = 2π), the pitching down motion induces an upward flow at the leading edge increasing the circulation; the LEV grows unstable and is shed from the plate. The LEV is stronger during the pitching up portion of the cycle than during the pitching down phase, i.e. the LEV is “helping” the plate pitch up more than it is resisting the pitching down, therefore the plate extracts a net positive amount of work from the fluid. of Ur = 3. We also observed an expansion in the range of angles of attack where the plate exhibited supercritical behavior due to an increase of aspect ratio from AR = 1 to 2, which we attribute to a reduction in three-dimensional effects on flow separation. In the forced pitching experiments, supercritical behavior was not initiated until the flow velocity was equivalent to Ur ≈ 10. This delay in the transition from sub- to super-critical shedding behavior runs counter to the trend observed for flexibly mounted plate where the transition occurred earlier as the aspect ratio was increased. This discrepancy is likely due to an apparent increase in the damping as a result of shaft misalignment. This flaw in the experimental apparatus has been eliminated in subsequent setup versions. Future work would include quantitative flow visualization, e.g., particle image velocimetry (PIV), of the pitching process to verify this hypothesis. We present an idealized model demonstrating the timing between the plate motion and the vortex shedding to explain what we view as the main difference between the subcritical (net negative work per cycle) and supercritical (net positive work per cycle) states. In Fig. 89 4.14, we illustrate the pitching and vortex shedding cycle for supercritical conditions. As the plate pitches down (t = π/2 to π), a leading edge vortex (LEV) forms and grows until the plate reaches the minimum angle of attack. Then, the plate reverses direction and begins to pitch up. The pitching up motion of the plate from t = π to 3π/2 acts to stabilize the LEV by reducing the rate of circulation being added to the vortex. When the plate reaches the end of the cycle (t = 2π), the pitching down motion increases the circulation shedding rate until the LEV grows unstable and is shed from the plate. The LEV is stronger during the pitching up portion of the cycle than during the pitching down phase resulting in a positive amount of work transferred from the fluid to the plate. This work was preliminary study of the flow conditions and body parameters (aspect ratio, resonant vibrational frequency, angle of attack, etc.) that contribute to the effective momentum transfer from the fluid to the flat plate. Future work should include mea- surements that capture the wake flow field, e.g., PIV, coupled with measurements of the resulting aerodynamic forces to better understand the mechanisms of vortex formation and shedding induced by body motion and would improve our understanding of flapping flight in the animal world and aid in the development of future engineered flapping flyers. Chapter 5 Vortex-induced vibrations of a flat plate mounted to a virtual spring and damper Arnold J. Song, Benjamin Strom, Kyohei Onoue, Kenneth S. Breuer 5.1 Introduction The work presented in the previous chapters examined how the aerodynamics of a compliant structure depended on the nature of the coupling between the structure and the fluid. In the last chapter, we focused on the unsteady aerodynamic forcing of a flat plate resulting from vortex shedding and how passive oscillations (where the plate was mounted to a torsion spring) or forced oscillations (where the plate was mounted to a servomotor and forced to follow a prescribed pitching path) affected the coupling between the plate and the fluid. The dynamics of an elastically mounted flat plate that pitches in response to fluid forces, τf , can be described by Iα ¨ (t) + bα(t) ˙ + k(α(t) − α0 ) = τf (t) (5.1) where the moment of inertia, viscous damping coefficient and torsional stiffness are denoted 90 91 by I, b and k, respectively; α˙ and α ¨ are the first and second derivatives with respect to time of the plate’s angular position, α, and the initial zero torque angle of attack is denoted by α0 . The amount of momentum transferred from the fluid to the flat plate is dependent upon the the matching between the vortex shedding and the structure’s dynamics; changes in this coupling are indicated by variation in the pitching amplitude. For the passively pitching plate, the structure’s resonant frequency, f0 , is determined by the torsional stiffness and 1 p moment of inertia of the structure, according to the relationship: f0 = 2π k/I. The coupling between the fluid and the plate depends the structure’s resonant frequency, which we can modify by varying the torsional stiffness, k, or the plate’s moment of inertia, I. If we normalize Eq. 5.1 by the plate’s moment of inertia, b k τf (t) ¨ (t) + α(t) α ˙ + (α(t) − α0 ) = (5.2) I I I we see that increasing the moment of inertia to modify f0 decreases the acceleration that the plate experiences due to fluid forces, τf , while also decreasing the effect of the viscous damping. Therefore, we would prefer changing f0 through increasing or decreasing the torsion spring stiffness to keep inertial effects consistent while varying f0 . To avoid inertial complications, we can perform forced oscillation experiments that allow us to examine changes in the momentum transfer from the fluid to the plate as we vary the pitching amplitude and frequency. In these types of experiments, however, the plate motion is independent of the fluid forcing, therefore lacks interplay between the structure dynamics and fluid forces that uniquely characterize fluid-structure interaction phenomena [20]. To preserve the essential character of the fluid-structure interactions, we developed an experiment that will allow us to vary the structural parameters that govern the fluid- 92 structure interaction easily and efficiently for broad and detailed sweeps of the structural parameter space. In this cyberphysical system, the structural dynamics are simulated in a computer simulation (the cyber part of the system) using the measured plate position, motion and fluid forces (in real, not virtual, space) as inputs for the virtual system. The results from this real-time simulation are then used to control a physical plant, e.g., in our case, the simulation determines the amount of torque a servo motor applies to our flat plate (the physical part of the system). An appealing aspect of the experimental workspace of a cyberphysical system is that we retain the essence of fluid-structure interactions and the structure’s dynamic response to the fluid forces is determined virtually allowing us to select arbitrary structural properties for the flat plate motion. With this system, we can study the autorotation (k = 0) of a flat plate with an arbitrary damping function (e.g., b = b(α, α, ˙ t)) or we can study an elastically mounted plate that is also being pitched actively, such that α0 (t) is a function of time, by only changing the functions that define k and b in Eq. 5.1, with the change in structural parameters being done in software with a push of the button. We are aware of at least three different cyberphysical systems that have been developed to study vortex-induced vibrations (VIV). All of these previous efforts examined VIV in water for a cylinder that oscillates perpendicular to the flow. The virtual cable testing apparatus (VCTA) developed at MIT by Hover and coworkers was the first cyberphysical system that utilized a servomotor coupled to virtually realize the vortex-induced dynamics of an elastically mounted cylinder [20]. For this cyberphysical system, the cylinder’s mass, damping and spring stiffness were virtually determined and the intended dynamics of the cylinder in response to the external fluid forcing were determined by a computer simulation using the cylinder’s position and force measurements as input. These experiments were conducted for a relatively low speeds at Reynolds numbers ranging from 7.2×103 −11.5×103 , 93 which are appropriate for slow moving ocean current flows over mooring lines. Lee and coworkers [22, 74] developed a virtual spring-damper system for VIV for cylinders that could be used in power harvesting applications. This series of experiments were conducted for Reynolds numbers ranging from 4 × 104 − 12 × 104 , with a maximum Reynolds an order of magnitude higher than the MIT experiments. Mackowski and Williamson also developed a cyberphysical elastically mounted cylinder [21] in water for which they demonstrated excellent agreement between the cyberphysical and purely physical experiments for Re = 4000. They demonstrate the potential of cyberphysical systems to accurately replicate the fluid-structure interaction of a purely physical experiments, but their results were still for a limited range of elastic parameters and relatively small structural resonant frequencies, f0 = O(≤ 1). Our cyberphysical system simulates the VIV behavior of a flexibly mounted flat plate with resonant structural frequencies that are an order of magnitude higher than observed in the previous efforts described above. In addition, our experiments are conducted for air flows for Reynolds number of order 104 to 105 , a range relevant to animal flight and micro air vehicles. In the next section, we describe the design and implementation of a cyberphys- ical system composed of a flat plate and servomotor that simulates the VIV behavior of an elastically mounted plate, with resonant frequencies up to 35 Hz. This description includes the physical apparatus design, the control scheme design and the identification of the fric- tional and viscous losses in the experimental apparatus. In section 3, we present preliminary results from free oscillation experiments that demonstrate the ability to virtually control the stiffness and damping; we compare the cyberphysical behavior to a numerically modeled response. We also present preliminary measurements of the VIV of a flat plate for several values of torsional stiffness at Re = 1 × 105 . The cyberphysical system implemented in this 94 Figure 5.1: Description of a flat plate mounted to a virtual spring-damper system. This figure illustrates the cyberphysical system used in our experiments that examine the vortex-induced pitching motion of a flat plate mounted to a virtual torsion spring and damper. The elastic-viscous behavior of the flat plate results not from a physical torsion spring and damper, but from a servomotor applying torque to the plate according to a control scheme that emulates the spring stiffness and damping specified in the virtual system. The titanium flat plate has a chord length of 15 cm with a span of 42 cm (plate AR = 2.8) and is 0.1 mm thick. chapter had virtual stiffness and damping. We conclude this chapter with a discussion of a cyberphysical system where the structure’s moment of inertia is also virtually specified. 5.2 Cyberphysical system design and setup We designed and implemented a cyberphysical system composed of a flat plate that pitches with one degree of freedom about an axis that runs along the plate’s half chord line. The pitching dynamics are governed by a torsional spring and damper that are virtually realized by a servomotor that actuates the plate’s rotational motion according to a control scheme 95 that applies a torque level dependent upon the plate’s position and velocity. We model the dynamics of this cyberphysical system as a forced harmonic oscillator, with the forcing torque, τf , being supplied by the fluid Ip α ˙ + kv (α[t] − α0 ) = τf [t] ¨ [t] + (bp + bv )α[t] (5.3) where the subscript, v, indicates virtual values of the viscous damping, b, and torsion spring stiffness, k, and the subscript, p, indicates the physical values of the plate’s moment of inertia, I, and damping. We define a reference angle of attack, α0 , that is equivalent to a relaxed, zero torque state for an elastic torsion spring and a variable equal to the angular excursion away from the reference angle of attack, or the difference between the angle of attack, α, and the reference angle of attack, α0 , which we denote as ∆α: ∆α = α(t) − α0 . (5.4) The torque produced by the virtual components is applied by the servo motor to the shaft, τa , τa = bv α˙ + kv ∆α = τf − (Ip α ¨ + bp α) ˙ (5.5) where Ip α ¨ is the rotational inertia of the flat plate and bp α˙ is the damping due to frictional and viscous loss mechanisms present in the physical system. The coefficients of damping, bv , and torsion spring stiffness, kv , governing the plate dynamics are mathematical abstractions allowing us to specify arbitrary values or functions for these parameters, therefore allowing us to test a floppy plate with low spring stiffness and low resonant frequency and a plate with a very large spring stiffness and high resonant frequency with the same setup by merely 96 changing these virtual parameters. Figure 5.1 illustrates the cyberphysical system used in experiments where we measured the vortex-induced pitching response of an elastically mounted flat plate. Our flat plate model is composed of a titanium plate with dimensions, 10 × 42 × 0.01 cm (chord × span × thickness) giving the plate an aspect ratio equal to 4.2. The plate is affixed to a 0.64 cm diameter steel shaft along the plate’s half chord, which is coupled to the servomotor on one end and a rotary encoder at the opposite end. The shaft is supported by ball bearings mounted to a frame that is rigidly attached to the wind tunnel. The plate is actuated by a brushless servo motor (Model: SM233AE, Parker Compumotor, Rohnert Park, CA) with a torque constant value of Kt = 0.415 N-m/A and a peak output torque of 3.4 N-m. A reaction torque sensor (Model: TFF400, Futek, Irvine, CA), also functioning as the mount for the servomotor to the frame, measures the torque applied by the motor, τa and a rotary encoder (Model: HB6M Optical Encoder, US Digital, Vancouver, WA) that outputs two digital quadrature pulse trains. The encoder signals are then converted to analog position (Model: EDAC2 Encoder Digital to Analog Converter, US Digital) and velocity signals (Model: ETACH2 High Speed Encoder to Analog Tachometer, US Digital) corresponding to the angular position, α and velocity, α, ˙ respectively. 5.2.1 Control scheme description Figure 5.2 illustrates the data acquisition and control scheme used for our cyberphysical system. The control scheme is composed of two feedback loops: 1) a current following loop with a servo drive (Model: DPRALTE-020B080, Advanced Motion Controls, Camarillo, CA) that supplies and monitors the current and voltage to the servo motor and 2) a loop that outputs a target torque level based on angular position and velocity measurements. 97 Dynamics simulation loop (1 kHz) α0 DAQ Current following loop (1 kHz input/output) (20 kHz) α α ~ α -+ k Vα τ Vτ Vτ Servo I + Vτ D/A controller/ α• . α + τ amplifier b Vα• Vα servo Vα Rotary α encoder motor A/D Vα• Vα• Tachometer Vτ Torque τ sensor Figure 5.2: Block diagram of virtual spring-damper feedback loop The flat plate accelerates in response to the torque applied by the motor, the rotational inertia, the system losses (e.g., bearing friction) and the fluid forces. The angular position, angular velocity and torque are sampled and recorded at 1 kHz by a data acquisition system (model: PCI-6259, National Instruments, Austin, TX) that is integrated with the cyber portion of our system; this includes the data logging and the dynamics simulation, which are controlled by a Matlab Simulink program. The angular position and velocity signals are received by the data acquisition system as analog voltages and are first converted into radians and radians per second to perform the simulation with units corresponding to physical values. The stiffness and damping gains are then applied to the position and velocity according to Eqn. 5.5 to determine the applied motor torque, τa . This value is converted to a voltage level dependent on the motor torque constant, Kt [N-m/A], and the I conversion constant for the input to the servo drive, Vτ [A/V]. This voltage signal that sets the target current (torque) level for the servo motor is the link between the cyber and physical sides of this dynamical system. 98 The servo drive and motor make up the physical portion of our cyberphysical system. The servo drive receives an analog voltage level from cyber component that specifies the amount of current to supply to a servo motor. The torque output for our motor is linearly proportional to the current supplied, i.e., the motor applies 0.415 N-m per ampere supplied. An internal proportional-integral-derivative (PID) controller in the servo drive attempts to supply the motor with the specified amount of current by measuring the current supplied to the servo motor and adjusting this current according to the PID gain settings for the loop. The update rate of the current following loop is 20 kHz (50 µs sample interval time). The current (torque) loop PID gains are software adjustable so that the user may tune the response characteristics of the motor, e.g., current output rise time, percent overshoot, etc. 5.2.2 Damping model identification 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 Torque [N−m] Torque [N−m] 0 0 −0.05 −0.05 −0.1 −0.1 −0.15 −0.15 −0.2 −0.2 −150 −100 −50 0 50 100 150 −150 −100 −50 0 50 100 150 Angular velocity [rad/s] Angular velocity [rad/s] (a) Linear damping model (b) Quadratic damping model Figure 5.3: Determination of viscous damping. We determined the damping losses inherent to the physical components of the cyberphysical system, τd , by measuring the torque needed to rotate the shaft (without the plate installed) at a constant angular velocity, α. ˙ We compare the linear (a) and quadratic (b) formulations of the damping present in the physical system. For implementation of the cyberphysical system, it is important to characterize the losses inherent to the system, which result from friction and viscous damping from the bearings and servo motor. We determined the damping losses inherent to the physical components of 99 Table 5.1: Linear damping model fit results bp α˙ = β1 α˙ + sgn(α)τ ˙ c Direction β1 [N-m/(rad/s)] τc [N-m] α˙ > 0 (+) 1.63 × 10−4 1.1 × 10−2 α˙ < 0 (−) 1.42 × 10−4 4.2 × 10−3 Table 5.2: Quadratic damping model fit results bp α˙ = β2 α˙ 2 + β1 α˙ + sgn(α)τ ˙ c Direction 2 β2 [N-m/(rad/s) ] β1 [N-m/(rad/s)] τC [N-m] α˙ > 0 (+) −4.10 × 10−7 2.10 × 10−4 9.9 × 10−3 α˙ < 0 (−) 1.34 × 10 −6 2.87 × 10−4 1.8 × 10−3 the cyberphysical system, τd , by measuring the torque needed to rotate the shaft (without the plate installed) with a constant angular velocity, α. ˙ The damping torque is composed of a Columb component, τc , which is opposes motion with a constant torque value that is independent of angular velocity, and a viscous component, β(α), ˙ which opposes motion with a torque that is a function of the angular velocity. We model the damping due to the internal mechanisms of the servo motor, bearings, and encoder as a sum of the Coulomb and viscous damping: τd = bα˙ = β(α) ˙ + sgn(α)τ ˙ c where the viscous damping, bα, ˙ can be described by a linear or quadratic function bp α˙ = β1 α˙ + sgn(α)τ ˙ c (linear) (5.6) bp α˙ = β2 α˙ 2 + β1 α˙ + sgn(α)τ ˙ c (quadratic) (5.7) where β1 and β2 are the coefficients of the linear and quadratic components. In Fig. 5.3, we compare the linear (Fig. 5.3(a)) and quadratic(Fig. 5.3(b)) formulations of the damp- ing present in the physical system. The damping versus angular velocity data were fit to the linear and quadratic viscous damping function with the results given in Table 5.1 and 100 5.2. The damping behavior was not symmetric. For the linear viscous damping model, the damping coefficient was slightly higher for velocities in the positive (α increasing) direc- tion and the Coulomb damping magnitude was more than 2.6 times larger in the positive direction. In contrast, for the quadratic damping model, the viscous damping components were higher for velocities in the negative (α decreasing), however, the Coulomb damping was more than 5.5 times larger in the positive direction. The asymmetry observed in the damping is likely due to an asymmetry in the electromotive behavior of the motor, since it is extremely unlikely to have a directionally dependent Coulomb damping in a rotary system. We used the linear damping model in our control scheme for this implementation of the cyberphysical system. 5.3 Cyberphysical system implementation 5.3.1 Free oscillation test We verified our ability to vary the torsion spring stiffness and achieve predictable behavior through a series of free oscillation experiments for several values of stiffness (0.25k0 , 0.5k0 , k0 , 2k0 and 3k0 , where k0 = 0.415N − m which is the amount of torque generated per am- pere supplied to the motor). These free oscillation tests were conducted with the shaft only, without the plate installed, to exclude the effects of the fluid in our verification of the cyber- physical dynamics. The shaft was released from ∆α = 40◦ = 0.70 rad for the low stiffness cases (k < k0 ) and ∆α = 40◦ = 0.70 rad for the higher stiffness cases (k ≥ k0 ), then allowed to oscillate and ring down until the motion damped out completely. For the cyberphysical system, no virtual damping was added or subtracted from the inherent damping and was assumed to be the values determined for the linear viscous damping model in the previous 101 section. In Fig. 5.3, we compare the measured motion of the cyberphysical system with numerical simulations of the flat plate dynamics using the parameters in Table 5.3. There is excellent agreement between the simulated and measured results for all of the stiffness values except for the lowest value, 0.25k0 , where we observed the cyberphysical plate oscil- lations at the correct frequency, but damping out faster than the simulation. For the other stiffness values, we observe the cyberphysical plate motion mirror the simulation behavior, but there is a slight difference in behavior where the cyberphysical plate oscillations take a little bit longer to completely damp out. 102 Time Series Plot: Time Series Plot: 0.6 0.8 Simulation Simulation Experiment Experiment 0.4 0.6 0.4 0.2 0.2 Position [rad] Position [rad] 0 0 −0.2 −0.2 −0.4 −0.4 −0.6 −0.6 −0.8 −0.8 1 1.5 2 2.5 1 1.5 2 2.5 Time[s] Time[s] (a) Reduced stiffness (b) Reduced stiffness k = 0.25k0 = 0.10 N-m/rad k = 0.5k0 = 0.21 N-m/rad Time Series Plot: Time Series Plot: 0.4 0.4 Simulation Simulation Experiment Experiment 0.3 0.3 0.2 0.2 0.1 0.1 Position [rad] Position [rad] 0 0 −0.1 −0.1 −0.2 −0.2 −0.3 −0.3 −0.4 −0.4 3 3.5 4 4.5 3 3.5 4 4.5 Time[s] Time[s] (c) Baseline stiffness (d) Increased stiffness k = k0 = 0.42 N-m/rad k = 2k0 = 0.83 N-m/rad Time Series Plot: 0.4 Simulation Experiment 0.3 0.2 0.1 Position [rad] 0 −0.1 −0.2 −0.3 −0.4 3 3.5 4 4.5 Time[s] (e) Increased stiffness k = 3k0 = 1.25 N-m/rad Figure 5.4: Free oscillation tests – Varying stiffness We compared the measured and simulated oscillation behaviors of several stiffness values, k = 0.25k0 , 0.5k0 , k0 , 2k0 and 3k0 , where k0 = 0.42 N-m was taken to be the baseline stiffness value for this series of tests. 103 Table 5.3: Free oscillation parameters – Varying stiffness Torsional stiffness, k [N-m] Resonant frequency, f0 [Hz] 0.10 (0.25k0 ) 4.2 0.21 (0.5k0 ) 5.9 0.42 (k0 ) 8.4 0.83 (2k0 ) 11.8 1.25 (3k0 ) 14.5 We conducted similar free oscillation tests while varying the amount of viscous damping that was either added or subtracted from the baseline damping (the damping that is inherent to the experimental apparatus), as shown in Fig. 5.5. We observe excellent agreement between the simulated and cyberphysical response even when the oscillation amplitude becomes very small when we increased the damping and relatively good agreement when the damping was reduced. But, as seen in Table 5.4, the damping values of the measured and simulated cyberphysical system were not always matched to achieve the agreement in the oscillation behavior. For the increased damping cases, we only added Coulomb damping to the simulated results, where we added τC = 2.0 × 10−3 and 1.5 × 10−2 when applying an additional 1.5 × 10−4 and 1 × 10−2 N-m/(rad/s), respectively, to the viscous damping coefficient. For the decreased damping cases, we also had to add a small amount of viscous damping in addition to the Coulomb damping to achieve agreement between the measured and simulation behavior. 5.3.2 Vortex-induced vibrations – Preliminary results This series of measurements was conducted at Re = 1 × 105 for several values of virtual spring stiffness, k = k0 , 2k0 and 4k0 , with no virtual damping added or subtracted from the system, so only the inherent damping due to frictional and viscous losses in the physical system contributed to the plate dynamics. For each run (corresponding to an individual 104 Time Series Plot: Time Series Plot: 0.4 0.3 Simulation Simulation Experiment Experiment 0.3 0.2 0.2 0.1 0.1 Position [rad] Position [rad] 0 0 −0.1 −0.1 −0.2 −0.2 −0.3 −0.3 −0.4 −0.4 3 3.5 4 4.5 3 3.5 4 4.5 Time[s] Time[s] (a) Add damping (b) Add damping β1,+ = 1.5 × 10−4 N-m/(rad/s) β1,+ = 1.0 × 10−3 N-m/(rad/s) Time Series Plot: Time Series Plot: 0.4 0.4 Simulation Simulation Experiment Experiment 0.3 0.3 0.2 0.2 0.1 0.1 Position [rad] Position [rad] 0 0 −0.1 −0.1 −0.2 −0.2 −0.3 −0.3 −0.4 −0.4 3 3.5 4 4.5 3 3.5 4 4.5 Time[s] Time[s] (c) Subtract damping (d) Subtract damping β1,+ = −2 × 10−4 N-m/(rad/s) β1,+ = −5 × 10−4 N-m/(rad/s) Figure 5.5: Free oscillation tests – Varying damping. Table 5.4: Free oscillation parameters – Varying damping Damping model: bp α˙ = (β1 + β1,+ ) α˙ + (τc + τc,+ ) Cyberphysical Simulation β1,+ τc,+ β1,+ τc,+ [N-m/(rad/s)] [N-m] [N-m/(rad/s)] [N-m] 1.5 × 10−4 0 1.5 × 10−4 2.0 × 10−3 1.0 × 10−3 0 1.0 × 10−3 1.5 × 10−2 −2.0 × 10 −4 0 −2.3 × 10−4 −2.5 × 10−3 −5.0 × 10−4 0 −5.3 × 10 −4 −5.5 × 10−3 105 2 2 10 10 k k 2k 2k 4k 4k Pitching amplitude rms [degrees] Pitching amplitude rms [degrees] 1 1 10 10 0 0 10 10 −1 −1 10 10 −2 −2 10 10 −10 0 10 20 30 40 50 −20 −10 0 10 20 30 40 50 Angle of attack [degrees] Angle of attack [degrees] (a) Vibrational amplitude, αrms (b) Vibrational amplitude, αrms Increasing angle of attack Decreasing angle of attack 0 0 10 10 k k 2k 2k 4k 4k −1 10 −1 10 Torque rms [N−m] Torque rms [N−m] −2 10 −2 10 −3 10 −3 −4 10 10 −10 0 10 20 30 40 50 −20 −10 0 10 20 30 40 50 Angle of attack [degrees] Angle of attack [degrees] (c) Torque amplitude, τrms (d) Torque amplitude, τrms Increasing angle of attack Decreasing angle of attack Figure 5.6: Vortex-induced vibrations of flat plate for Re = 1 × 105 for several values of the virtual spring stiffness, k = 0.42, 0.83, 1.66 (k0 , 2k0 , 4k0 ). value of spring stiffness), the plate is initially positioned so that α0 = 0◦ , then α0 was first increased in 2 degree increments, then decreased in 2 degree increments. We plot the vibrational amplitude, ∆αrms (Figs. 5.6(a) and 5.6(b)), and the unsteady torque amplitude, (τm )rms (Figs. 5.6(c) and 5.6(d)), as a function of the time-averaged angle of attack, α ¯. For the lowest stiffness (k = k0 ), the plate immediately jumps to a relatively high angle ¯ = 23◦ , and the pitching amplitude remains below 0.5 degrees as the mean of attack, α angle of attack increases to nearly 50◦ . As we increased the stiffness by a factor of two, k = 2k0 , we observe similar behavior except the initial jump is not quite as large and the 106 ¯ = 12◦ , but the vibrational amplitude remained relatively plate immediately snaps to α ˜ rms < 0.5◦ . For α small, α ¯ ≥ 23◦ and greater, the unsteady behavior of the k0 and 2k0 cases followed the same trend with increasing mean angle of attack, with 2k0 spring vibrating with a slightly smaller amplitude, and remained in a subcritical state. The highest stiffness tested, k = 4k0 , exhibited quite different behavior and showed a transition from subcritical to supercritical coupling. This spring stiffness provided enough restoring torque to maintain low values of α ¯ for low values of α0 . As the offset angle was increased, the plate transitions to ¯ = 7◦ with a pitching amplitude of 2 degrees and reaches a maximum a supercritical state at α ¯ = 31◦ . Further increasing α amplitude of 20 degrees at α ¯ results in a transition back to a subcritical state. The plate remains in an subcritical state as α ¯ decreases with the pitching amplitude mirroring the behavior of the lower values of stiffness. As α ¯ is decreased further, ¯ = 7◦ , but immediately transition the plate transitions briefly to a supercritical state at α ¯ = 5◦ . This series of measurements demonstrates back down to a subcritical amplitude at α that we can initiate VIV behavior for our cyberphysical system that is consistent with a purely physical elastically mounted flat plate. 5.4 Future work 5.4.1 Virtual moment of inertia The cyberphysical system discussed and implemented above specified a virtual spring stiff- ness and viscous damping to determine the elastic-viscous response of the plate to fluid forcing. We now consider the case where the plate’s moment of inertia is also virtual where the expression for the forced harmonic oscillator would take the form Iv α ˙ + kv (α[t] − α0 ) = τf [t] ¨ [t] + (bp + bv )α[t] (5.8) 107 where, as before, the subscript, v, denotes virtual values of moment of inertia, I, viscous damping, b, and torsion spring stiffness, k, and the physical damping inherent to the ex- perimental apparatus due to frictional and viscous losses are indicated by bp . We cannot measure the torque exerted by the fluid to the plate, τf , directly. Using the same arrange- ment as in the previous series of experiments, we use a reaction torque sensor to measure the torque the motor applies to the plate, which is an amalgam of the physical plate’s inertia, the physical damping and the fluid forces, τm = τf − Ip α ¨ − bp α˙ (5.9) ⇒ τf = τm + Ip α ¨ + bp α˙ (5.10) which we can then substitute into Eq. 5.8 to obtain (Iv − Ip )α ¨ + (bv − bp )α˙ + kv α = τm . (5.11) If we are able to either estimate or measure the angular acceleration, α ¨ , then we can modify the control scheme that we used for the virtual spring-damper system by using Eq. 5.8 to calculate the target motor torque output. But, this method requires an additional sensor that either measures the angular acceleration directly (using an angular accelerometer) or indirectly (using a multi-axis linear accelerometer). Or the acceleration can be estimated by passing the angular position or velocity signals through a differentiation filter, how- ever, this approach will tend to amplify the noise and introduce phase lag issues, so direct measurement would be preferable. The approach of Hover et al. [20] used in the implementation of their VCTA system only needed direct measurement of the position and velocity to simulate the virtual dynamics. 108 For this method, we express Eq. 5.11 as a system of equations:         bv −bp d  α˙  − Iv −Ip − Ivk−I v p  α˙  1  Iv −Ip  τm   =   +    (5.12) dt         α 1 0 α 0 0 which we integrate to estimate the target plate velocity. So, instead of sending a target torque set point for the motor’s current following control loop, this approach estimates a target velocity for a velocity following loop. 5.5 Concluding remarks The aim of this chapter was to describe the design and implementation of a cyberphysical system and to demonstrate the feasibility of using such a system for studying VIV of a flat plate in air. We showed that we could predictably vary the torsional stiffness and damping to control the dynamics of the flat plate and that we could achieve VIV behavior that is consistent with flat plate mounted to a torsional spring-damper system. This type of system shows great promise as a versatile and efficient environment to study many different types of fluid-structure interaction phenomena, especially with the integration of the ability to specify a virtual value of the plate’s moment of inertia through software. Bibliography [1] McCroskey, W. & Fisher, R. K. (1972), Detailed aerodynamic measurements on a model rotor in the blade stall regime, Journal of the American Helicopter Society 17, 20-30. [2] Swartz, S. M., Groves, M. S., Kim, H. D., & Walsh, W. R. (1996), Mechanical prop- erties of bat wing membrane skin, Journal of Zoology 239, 357-378. [3] Swartz, S. M., Bennett, M. B., & Carrier, D. R. (1992), Wing bone stresses in free flying bats and the evolution of skeletal design for flight, Nature 359, 726-729. [4] Lissaman, P. (1983), Low-Reynolds-number airfoils, Annual Review of Fluid Mechanics 15, 223-239. [5] Shyy, W., Klevebring, F., Nilsson, M., Sloan, J., Carroll, B., & Fuentes, C. (1999), Rigid and flexible low Reynolds number airfoils, Journal of Aircraft 36, 523-529. [6] Torres, G. & Mueller, T. (2001), Aerodynamic characteristics of low aspect ratio wings at low Reynolds numbers, in Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, edited by Mueller, T. AIAA. Reston, VA. chapter 7, pp. 115–141. [7] Bishop, K. (2006), The relationship between 3-D kinematics and gliding performance in the southern flying squirrel, Glaucomys volans, Journal of Experimental Biology 209, 689-701. [8] Norberg, U. M. (1990), Vertebrate flight: mechanics, physiology, morphology, ecology and evolution, Springer-Verlag. Berlin. [9] Swartz, S. M., Bishop, K. L., & Ismael-Aguirre, M.-F. (2005), Dynamic complexity of wing form in bats: Implications for flight performance, in Functional and evolutionary ecology of bats, edited by Akbar, Z., McCracken, G., & Kunz, T. H. Oxford University Press. Oxford. [10] Song, A., Tian, X., Israeli, E., Galvao, R., Bishop, K., Swartz, S., & Breuer, K. (2008), Aeromechanics of membrane wings with implications for animal flight, AIAA Journal 46, 2096-2106. [11] Saffman, P. & Sheffield, J. (1977), Flow over a wing with an attached free vortex, Studies in Applied Mathematics 57, 107–117. [12] Maxworthy, T. (1979), Experiments on the weis-fogh mechanism of lift generation by insects in hovering flight. part 1. dynamics of the ’fling’, J. Fluid Mech. 93, 47-63. 109 110 [13] Ellington, C., Berg, C., Willmott, A., & Thomas, A. (1996), Leading-edge vortices in insect flight, Nature 384. [14] Hubel, T. & Tropea, C. (2010), The importance of leading edge vortices under sim- plified flapping flight conditions at the size scale of birds, Journal of Experimental Biology 213, 1930-1939. [15] Dickinson, M. H. & Gotz, K. G. (1996), The wake dynamics and flight forces of the fruit fly Drosophila melanogaster, Journal of Experimental Biology 199, 2085-2104. [16] Ringuette, M., Milano, M., & Gharib, M. (2007), Role of the tip vortex in the force generation of low-aspect-ratio normal flat plates, Journal of Fluid Mechanics 581, 453. [17] Milano, M. & Gharib, M. (2005), Uncovering the physics of flapping flat plates with artificial evolution, Journal of Fluid Mechanics 534, 403–409. [18] Pitt Ford, C. & Babinsky, H. (2013), Lift and the leading-edge vortex, Journal of Fluid Mechanics 720, 280–313. [19] Lentink, D. & Dickinson, M. H. (2009), Rotational accelerations stabilize leading edge vortices on revolving fly wings, Journal of Experimental Biology 212, 2705–2719. [20] Hover, F., Miller, S., & Triantafyllou, M. (1997), Vortex-induced vibration of marine cables: Experiments using force feedback, Journal of Fluids and Structures 11, 307- 326. [21] Mackowski, A. W. & Williamson, C. H. (2011), Developing a cyber-physical fluid dy- namics facility for fluid–structure interaction studies, Journal of Fluids and Structures 27, 748–757. [22] Lee, J., Xiros, N., & Bernitsas, M. (2011), Virtual damper-spring system for VIV experiments and hydrokinetic energy conversion, Ocean Engineering 38, 732-747. [23] Holbrook, K. A. & Odland, G. (1978), A collagen and elastic network in the wing of the bat, Journal of Anatomy 126, 21-36. [24] Ellington, C. P., vandenBerg, C., Willmott, A. P., & Thomas, A. L. R. (1996), Leading- edge vortices in insect flight, Nature 384, 626-630. [25] Muijres, F. T., Johansson, L. C., Barfield, R., Wolf, M., Spedding, G. R., & Heden- strom, A. (2008), Leading-edge vortex improves lift in slow-flying bats, Science 319, 1250-1253. [26] Warrick, D. R., Tobalske, B. W., & Powers, D. R. (2005), Aerodynamics of the hovering hummingbird, Nature 435, 1094-1097. [27] Hedrick, T. L., Tobalske, B. W., & Biewener, A. A. (2002), Estimates of circulation and gait change based on a three-dimensional kinematic analysis of flight in cockatiels (nymphicus hollandicus) and ringed turtle-doves (streptopelia risotia), Journal of Ex- perimental Biology 205, 1389-1409. [28] Shapiro, R. (1978), Direct linear transformation method for three-dimensional cine- matography, Research Quarterly 49, 197-205. 111 [29] Watts, P., Mitchell, E. J., & Swartz, S. M. (2001), A computational model for es- timating the mechanics of horizontal flapping flight in bats: Model description and validation, Journal of Experimental Biology 204, 2873-2898. [30] Ellington, C. P. (1984), The aerodynamics of hovering insect flight. Parts I-VI, Philo- sophical Transactions of the Royal Society of London Series B-Biological Sciences B 305, 1-181. [31] Sane, S. P. & Dickinson, M. H. (2002), The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight, Journal of Experimental Biology 205, 1087-1096. [32] Lehmann, F., Sane, S., & Dickinson, M. (2005), The aerodynamic effects of wing-wing interaction in flapping insect wings, Journal of Experimental Biology 208, 3075-3092, (http://dx.doi.org/10.1242/jeb.01744). [33] Tian, X., Iriarte-Diaz, J., Middleton, K., Galvao, R., Israeli, E., Roemer, A., Sullivan, A., Song, A., Swartz, S., & Breuer, K. (2006), Direct measurements of the kinematics and dynamics of bat flight, Bioinspiration & Biomimetics 1, S10-S18. [34] Riskin, D. K., Willis, D. J., Iriarte-Diaz, J., Hedrick, T. L., Kostandov, M., Chen, J., Laidlaw, D. H., Breuer, K. S., & Swartz, S. M. (2008), Quantifying the complexity of bat wing kinematics, Journal of Theoretical Biology 254, 604-615. [35] Hedenstrom, A., Johansson, L., Wolf, M., von Busse, R., Winter, Y., & Spedding, G. (2007), Bat flight generates complex aerodynamic tracks, Science 316, 894-897, (http://dx.doi.org/10.1126/science.1142281). [36] Lian, Y. & Shyy, W. (2007), Laminar-turbulent transition of a low Reynolds number rigid or flexible airfoil, AIAA Journal 45, 7, (http://dx.doi.org/10.2514/1.25812). [37] Pelletier, A. & Mueller, T. (2000), Low Reynolds number aerodynamics of low-aspect- ratio thin/flat/cambered-plate wings, AIAA Journal 37, 825-832. [38] Torres, G. & Mueller, T. (2004), Low-aspect-ratio wing aerodynamics at low Reynolds numbers, AIAA Journal 42, 523-529. [39] Nielsen, J. (1963), Theory of flexible aerodynamic surfaces, Journal of Applied Me- chanics 30, 435-442. [40] Thwaites, B. (1961), The aerodynamic theory of sails. I. Two-dimensional sails, Pro- ceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 261, 402-422, (http://dx.doi.org/10.1098/rspa.1961.0086). [41] Jackson, P. (1983), A simple model for elastic two-dimensional sails, AIAA Journal 21, 153-155. [42] Lorillu, O., Weber, R., & Hureau, J. (2002), Numerical and experimental analysis of two-dimensional separated flows over a flexible sail, Journal of Fluid Mechanics 466, 319-341, (http://dx.doi.org/10.1017/s0022112002001283). [43] Newman, B. & Low, H. (1984), Two-dimensional impervious sails: Experimen- tal results compared with theory, Journal of Fluid Mechanics 144, 445-462, (http://dx.doi.org/10.1017/s0022112084001683). 112 [44] Shyy, W., Berg, M., & Ljungqvist, D. (1999), Flapping and flexible wings for biological and micro air vehicles, Progress in Aerospace Sciences 35, 455-505. [45] Shyy, W., Ifju, P., & Viieru, D. (2005), Membrane wing-based micro air vehicles, Applied Mechanics Reviews 58, 283-301. [46] Lian, Y. S., Shyy, W., Viieru, D., & Zhang, B. N. (2003), Membrane wing aerodynamics for micro air vehicles, Progress in Aerospace Sciences 39, 425-465, Review. [47] Stanford, B., Sytsma, M., Albertani, R., Viieru, D., Shyy, W., & Ifju, P. (2007), Static aeroelastic model validation of membrane micro air vehicle wings, AIAA Journal 45, 2828-2837, (http://dx.doi.org/10.2514/1.30003). [48] Swartz, S., Groves, M., Kim, H., & Walsh, W. (1996), Mechanical properties of bat wing membrane skin, Journal of Zoology 239, 357-378. [49] Mikhail, E., Bethel, J., & McGlone, J. (2001), Introduction to Modern Photogrammetry, John Wiley & Sons, Inc. [50] Iriate-Diaz, J., Riskin, D. K., Willis, D. J., Breuer, K. S., & Swartz, S. M. (2011), Whole-body kinematics of a fruit bat reveal the influence of wing inertia on body accelerations, Journal of Experimental Biology 214, 1546-1553. [51] Null, W. & Shkarayev, S. (2005), Effect of camber on the aerodynamics of adaptive- wing micro air vehicles, Journal of Aircraft 42, 1537-1542. [52] Bishop, K. (2007), Aerodynamic force generation, performance and control of body ori- entation during gliding in sugar gliders (Petaurus breviceps), Journal of Experimental Biology 210, 2593-2606. [53] Galvao, R., Israeli, E., Song, A., Tian, X., Bishop, K., Swartz, S., & Breuer, K. (2006), The aerodynamics of compliant membrane wings modelled on mammalian flight mechanics, in 24th AIAA Applied Aerodynamics Conference. [54] Rayleigh, J. W. S. (1894), The Theory of Sound, MacMillan. London. 2nd edition. [55] Yerusevych, S., Sullivan, P., & Kawall, J. (2006), Coherent structures in an airfoil boundary layer and wake at low Reynolds numbers, Physics of Fluids 18. [56] Weis-Fogh, T. (1973), Quick estimates of flight fitness in hovering animals including novel mechanisms for lift production, Journal of Experimental Biology 59, 169-230. [57] Wagner, H. (1925), Z. angew, Math. Mech 5, 17. [58] Lighthill, M. (1973), On the Weis-Fogh mechanism of lift generation, Journal of Fluid Mechanics 60, 1-17. [59] Dickinson, M. H., Lehmann, F., & Sane, S. P. (1999), Wing rotation and the aerody- namic basis of insect flight, Science 284, 1954-1960. [60] Huang, M.-K. & Chow, C.-Y. (1982), Trapping of a free vortex by Joukowski airfoils, AIAA Journal 20, 292-298. 113 [61] Bomphrey, R. J., Lawson, N. J., Harding, N. J., Taylor, G. K., & Thomas, A. L. R. (2005), The aerodynamics of Manduca sexta: Digital particle image velocimetry anal- ysis of the leading-edge vortex, Journal of Experimental Biology 208, 1079-1094. [62] Gerrard, J. (1966), The mechanics of the formation region of vortices behind bluff bodies, Journal of Fluid Mechanics 25, 401-413. [63] Roshko, A. (1954), On the development of turbulent wakes from vortex streets, Tech- nical Report 1191. National Advisory Committee for Aeronautics. [64] Fage, A. & Johansen, F. (1927), On the flow of air behind an inclined flat plate of infinite span, Proceedings of the Royal Society of London. Series A 116, 170-197. [65] Blevins, R. D. (1990), Flow-Induced Vibrations, Van Nostrand Reinhold. [66] Bearman, P. (1984), Vortex shedding from oscillating bluff bodies, Annual Review of Fluid Mechanics 16, 195–222. [67] McCroskey, W. (1982), Unsteady airfoils, Annual Review of Fluid Mechanics 14, 285-311. [68] Theodorsen, T. & Garrick, I. (1940), Mechanism of flutter: A theoretical and experi- mental investigation of the flutter problem, Technical Report 496. National Advisory Committee for Aeronautics. [69] von Karman, T. & Sears, W. (1938), Airfoil theory for non-uniform motion, Journal of Aerospace Sciences 5, 379-390. [70] Dowell, E. H. & Hall, K. C. (2001), Modeling of fluid-structure interaction, Annual Review of Fluid Mechanics 33, 445–490. [71] Dowell, E. H. (2012), Nonlinear oscillations of a fluttering plate., AIAA Journal 4. [72] Jones, M. & Shelley, M. (2005), Falling cards, Journal of Fluid Mechanics 540, 393– 426. [73] Williamson, C. & Govardhan, R. (2004), Vortex-induced vibrations, Annu. Rev. Fluid Mech. 36, 413–455. [74] Lee, J. & Bernitsas, M. (2011), High-damping, high-Reynolds VIV tests for energy harnessing using the VIVACE converter, Ocean Engineering 38, 1697-1712.