Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Wing selection and design for flapping flight at small scale
(USC Thesis Other)
Wing selection and design for flapping flight at small scale
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Wing Selection and Design for Flapping Flight at Small Scale by Emma K. Singer A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (MECHANICAL ENGINEERING) May 2022 Copyright 2022 Emma K. Singer I dedicate this dissertation to the family and friends who showed me unwavering support during my graduate studies. As John Donne said, “no man is an island," and I cannot say that I would have made it to the finish line without the wonderful people that I share my life with. In particular, I would like to thank my father Paul, who taught me (from personal experience) that it is okay to take eight years to get a PhD. ii Acknowledgements The author acknowledges the financial support of the Department of Aerospace & Mechanical Engineering in the Andrew and Erna Viterbi School of Engineering at the University of Southern California. I would like to thank my advisor Dr. Geoffrey R. Spedding for his guidance, enthusiasm, and patience while we delved into parts unknown. I could not have asked for a better advisor; I enjoyed every minute of our mutual uncertainty about my research. I would also like to thank Dr. Néstor O. Pérez-Arancibia, who believed in me from the beginning, and gave me the opportunity to work on exciting projects that I never would have imagined being part of. I would like to thank my committee members Dr. Mitul Luhar and Dr. Steve Nutt, for being part of my qualifying and defense committee, as well as providing valuable feedback on this work. Finally, I would like to thank all of my fellow graduate students at the AMSL and the Stratified Flow Laboratory, Chan- ye Ohh in particular, who helped make research engaging and fulfilling during late nights and tight deadlines. iii TableofContents Dedication ii Acknowledgements iii ListofTables vi ListofFigures vii Abstract xi Introduction 1 0.1 Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Chapter1: Design&FabricationoftheFlapping-WingMicroAirVehicle 9 1.1 The Flapping-Wing Microrobotic Insect . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.2.1 Carbon Fiber Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.2.2 Composite Fabrication and the SCM Method . . . . . . . . . . . . . . . . . 22 1.2.3 Advanced Carbon Fiber Composite Fabrication Techniques . . . . . . . . 24 1.2.4 Bimorph Piezoelectric Cantilever Actuators . . . . . . . . . . . . . . . . . 28 Chapter2: Forcemeasurementsatmicro-roboticscale 31 2.1 Clip–Brazing for the Design and Fabrication of millimeter-scale, micronewton- resolution Force Sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.1.1 Fabrication of mm-scale structures . . . . . . . . . . . . . . . . . . . . . . 33 2.1.1.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.1.1.2 Potential Applications . . . . . . . . . . . . . . . . . . . . . . . 38 2.1.2 The use of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.1.2.1 Micro-Force Sensor Design . . . . . . . . . . . . . . . . . . . . . 39 2.1.2.2 Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.1.2.3 Evaluation of the Fabrication Method . . . . . . . . . . . . . . . 51 2.1.3 Calibration And Validation . . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.1.3.1 Experimental Force Validation . . . . . . . . . . . . . . . . . . . 59 2.1.3.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.1.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.2 Measurements & Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 iv Chapter3: TheEffectofChordwiseWingFlexibilityonMicroroboticForceGener- ation 74 3.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.1.1 Experimental Platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.1.2 Force Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.1.3 Particle Image Velocimetry (PIV) . . . . . . . . . . . . . . . . . . . . . . . 94 3.1.4 Practical Considerations for Tomographic-PIV . . . . . . . . . . . . . . . . 96 3.1.5 Vector Field Smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.1.6 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.2.1 Downwash . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 3.2.2 Spanwise Vorticity and Circulation . . . . . . . . . . . . . . . . . . . . . . 113 3.2.3 Force-Torque Time Traces . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 3.2.4 Vortex Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 3.4 Comparison to the FW-MAV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 3.5 Summary and conclusions for flapping model . . . . . . . . . . . . . . . . . . . . 141 OverallConclusions 143 3.6 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 References 146 v ListofTables 1.1 Actuator Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1 Material Properties of Invar-36 and Silver . . . . . . . . . . . . . . . . . . . . . . . 41 2.2 Geometric Parameters of the Structure. . . . . . . . . . . . . . . . . . . . . . . . . 45 2.3 Inertial and Dynamic Parameters of the Structure. . . . . . . . . . . . . . . . . . . 45 2.4 Mechanical Failure Properties - Tensile Test . . . . . . . . . . . . . . . . . . . . . 58 2.5 Measured Average Periodic Lift and Absolute Error by Frequency . . . . . . . . . 60 2.6 Modeled vs Empirical Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.1 Scaling and dimensionless Parameters (Π-Variables) of the system. . . . . . . . . . 81 3.2 Summary of volume self-calibration polynomial fit errors for downwash and span- wise experimental configurations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.3 Force estimatorF z for downwash integrals . . . . . . . . . . . . . . . . . . . . . . 112 3.4 Average values of the coefficients of lift and drag over three periods . . . . . . . . 123 3.5 Expected trailing edge deflection by wing . . . . . . . . . . . . . . . . . . . . . . . 125 vi ListofFigures 1 Types of flapping-wing motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1 The 70-mg Double-Actuator FW-MAV . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Piezoelectric Bimorph Cantilever Actuator . . . . . . . . . . . . . . . . . . . . . . 11 1.3 Transmission Motion Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4 Transmission Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.5 Passive hinge and flapping motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.6 Modeling of the Actuator-Transmission-Wing System . . . . . . . . . . . . . . . . 15 1.7 Dynamic Model of the System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.8 Carbon fiber layups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.9 Out-of-Autoclave Composite Curing . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.10 Equipment and Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.11 Single-Actuator Microrobotic Bee . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.12 SCM-Type Composites used for Microrobotic Designs . . . . . . . . . . . . . . . . 23 1.13 The SCM Fabrication Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.14 File preparation and folding design . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.15 Folding Method vs 5-Layer Composite for Transmission Fabrication . . . . . . . . 25 1.16 Five-Layer Transmission Layup (No Folding) . . . . . . . . . . . . . . . . . . . . . 26 1.17 Example Experimental Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.18 Unibody Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.19 Batch Piezoelectric Cantilever Fabrication . . . . . . . . . . . . . . . . . . . . . . 29 1.20 Actuator Connections and Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 vii 2.1 The Micro-Force Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2 Braze Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3 Force sensor model and mechanism of action . . . . . . . . . . . . . . . . . . . . . 39 2.4 Analytic and numeric Models and their expected frequency response . . . . . . . 43 2.5 Fabrication process of the sensor using folding and clip-brazing . . . . . . . . . . 47 2.6 Inspection with Scanning Electron Microscopy . . . . . . . . . . . . . . . . . . . . 51 2.7 Experimental Setup and Input-Output Verification . . . . . . . . . . . . . . . . . . 52 2.8 Calibration and Validation with Weight Testing . . . . . . . . . . . . . . . . . . . 54 2.9 Tensile Test: Preparation and Results . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.10 Spread of the calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.11 Regression Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.12 Flapping-wing force measurement experiments and comparisons with Computa- tional Fluid Dynamics simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.13 Early Data and Displacement Attenuation . . . . . . . . . . . . . . . . . . . . . . 66 2.14 Instantaneous lift forces for the Single-Actuator Microrobotic Bee . . . . . . . . . 67 2.15 Single-Actuator Microrobotic Bee: instantaneous and average forces . . . . . . . . 68 2.16 Double-Actuator Bee: Attenuation of output flapping angle at high frequencies . . 68 2.17 Single- vs Double-Actuator Microrobotic Bee . . . . . . . . . . . . . . . . . . . . . 69 2.18 System Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 2.19 Double-Actuator Bee: Actuator Displacement vs Flapping Force . . . . . . . . . . 71 2.20 Consistency of Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 2.21 Lift force vs flapping frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2.22 Comparison of Instantaneous Force Measurements at various frequencies . . . . . 73 3.1 Scaled system gearbox and end-effector . . . . . . . . . . . . . . . . . . . . . . . . 75 3.2 Scaled mean chord versus reduced flapping frequency . . . . . . . . . . . . . . . . 77 3.3 Comparison of original wing to scaled model . . . . . . . . . . . . . . . . . . . . . 80 3.4 Tomographic-PIV configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.5 CAD model of the scaled flapping wing system . . . . . . . . . . . . . . . . . . . . 85 viii 3.6 Flapping kinematic profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.7 Schematic of laboratory equipment . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.8 Scaled wing and coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.9 Raw force-torque signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.10 Frequency response ofF c ,F n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.11 Comparison of cutoff frequency on F c ,F n . . . . . . . . . . . . . . . . . . . . . . . 91 3.12 Frequency response of filtered system . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.13 Inertial subtraction of non-wing components . . . . . . . . . . . . . . . . . . . . . 93 3.14 Principles of Tomographic-PIV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.15 Illumination intensity profiles of reconstructed tomographic volumes . . . . . . . 98 3.16 Smoothing synthetic test signals corrupted with Gaussian noise . . . . . . . . . . 102 3.17 Smoothing algorithm sensitivity to sampling rate . . . . . . . . . . . . . . . . . . 103 3.18 The use of the smoothing algorithm on a Tomographic-PIV dataset . . . . . . . . 104 3.19 Choice of cubic spline smoothing parameter . . . . . . . . . . . . . . . . . . . . . 107 3.20 Laser volume illumination relative to wing position . . . . . . . . . . . . . . . . . 109 3.21 Downwash force estimatorF z over two flapping periods . . . . . . . . . . . . . . 111 3.22 Out-of-plane vorticityω z threshold limits . . . . . . . . . . . . . . . . . . . . . . . 114 3.23 Wings 1 & 2: Vorticity plots at quarter-period intervals . . . . . . . . . . . . . . . 117 3.24 Wings 3 & 4: Vorticity plots at quarter-period intervals . . . . . . . . . . . . . . . 118 3.25 Circulation by depth plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 3.26 Instantaneous coefficients of lift and drag for all wings . . . . . . . . . . . . . . . 121 3.27 Cycle asymmetry ofC L ,C D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 3.28 Trailing edge deflection and cantilever model . . . . . . . . . . . . . . . . . . . . . 124 3.29 Coefficients of lift and drag versus flapping kinematics . . . . . . . . . . . . . . . 126 3.30 Vorticity comparison att/T = 0.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 3.31 Vorticity comparison att/T = 0.85 . . . . . . . . . . . . . . . . . . . . . . . . . . 128 3.32 Vorticity isosurfaces att/T = 0, development of LEV and TEV . . . . . . . . . . . 129 ix 3.33 λ 2 -criterion isosurfaces att/T = 0.4 . . . . . . . . . . . . . . . . . . . . . . . . . . 130 3.34 λ 2 -criterion isosurfaces att/T = 0.6 . . . . . . . . . . . . . . . . . . . . . . . . . . 131 3.35 λ 2 -criterion isosurfaces att/T = 1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . 132 3.36 Coefficient of lift time traces: scaled model vs original . . . . . . . . . . . . . . . . 138 3.37 Selected coefficient of lift plots from current literature . . . . . . . . . . . . . . . . 140 x Abstract For small biologically-inspired flying machines such as sub-gram flapping wing micro-air-vehicles (FW-MAVs) which operate in flows with Reynolds number Re≤ 10 3 , unsteady aerodynamic mechanisms are critical in attaining net positive lift, control authority and payload capacity. The aerodynamic performance is sensitive to the wing stiffness, with complex three-dimensional and time-varying flow structures evolving along and behind the wing. This work provides a method for the design and fabrication of a microNewton-resolution dynamic force sensor to quantify flight performance metrics, but the micro-robotic vehicles are known for poor repeatability of the complex fabrication processes and highly coupled mechanical subsystems. It is difficult to study the fluid flow around the wings at-scale due to the physical size and actuation frequency of the system. Therefore this work uses a dynamically-scaled model undergoing prescribed sinu- soidal flapping and pitching motion with zero mean flow, to study the induced flow and forces during a hovering flight condition. Tomographic PIV is used to visualize the three-dimensional vortex formation and force-generating flow structures while a 6-axis force-torque sensor at the wing root captures instantaneous data throughout the flapping period. The coefficients of lift and drag are computed for several cycles of flapping motion for four different wings that vary in chordwise stiffness and the traces are related to the flow structures observed in the induced velocity fields. The results show that chordwise flexible flapping wings in this Reynolds number regime demonstrate higher mean lift and lower RMS drag than their rigid counterparts; therefore chordwise FW-MAV wings have the potential to increase payload capacity and power economy without adding any additional weight to the vehicle. xi Introduction Research on sub-gram microrobotic vehicles is multidisciplinary and covers a wide range of topics from studies on new engineering materials and micron-to-millimeter scale fabrication methods to actuation, power systems, and controller design. The success of this research requires the use of the strongest, lightest materials and efficient power sources, necessitating a broad search across many scientific fields and constant iteration of designs and procedures. The first of its kind, the sub-gram barrier for microrobotic fliers was broken with the research in [1] with the use of an ultra-light carbon fiber frame and bio-inspired gearless flapping-wing mechanisms, but this system is limited by the tethered external power system. Over a decade of research, these directly-powered robots became controllable [2, 3], more reliably designed [4], and more powerful [5]. This perpetual improvement of the flapping-wing micro air vehicle (FW-MAV) system has enabled the newest designs to achieve free (untethered) flight [6, 7], pushing the robotic designs closer to practical use. From mechanical pollination, distributed sensing swarm networks, to search and rescue efforts, the potential impact is significant. However, scaled production, power autonomy, onboard computational power and control structure must be significantly improved to achieve these goals. A major obstacle to achieve free flight is the limitations of the energy density of practical power sources for the robots. Currently, even the most efficient batteries at a physical scale able to be incorporated into the chassis (approximately 2x2x14 mm) are incapable of providing the required energy for sustained flight and would require on-board power electronics which decrease payload capacity given a fixed maximum aerodynamic lift production. It is not unrea- sonable to expect battery performance in general to improve in the near future, to the point at 1 which it is feasibly incorporated into a ground-based-locomotion microrobot, but this weight is prohibitive for the flapping-wing flight vehicles. The lower minimum actuation frequencies for walking/crawling motion compared to flapping enable alternative actuation technology such as shape memory alloys (SMAs) as in [8]. However, locomotion has a similar (or exacerbated) sur- face navigation problem as current, state-of-the-art macro-robots which are known to have only moderate performance on uneven and unstable terrain. Therefore, we consider flight as the pre- ferred motion for microrobotic designs given the previously-stated goals, to avoid ground-based obstacles and the necessity of secondary subsystems for activities such as climbing. The promis- ing new results using photovoltaics [7] and laser-delivered power [6] have been demonstrated as a viable alternative to batteries for the power autonomy problem for sustained flight. Despite the achievements in bio-inspired FW-MAV design and control research, there remain many unanswered questions about the aerodynamics of their flight. Steady-state aerodynamic theory alone is insufficient to explain the high aerodynamic force coefficients during flapping [9], necessitating the existence of unsteady mechanisms in the flow. Flapping flight at small scale in fact relies on a number of these unsteady effects. These include phenomena such as delayed stall [10–13], vortex formation at the stroke reversal [12, 14] and added-mass effects [15, 16], each of which have been proposed as potential explanations of the enhanced lift over that predicted by steady state models. Direct Numerical Simulations (DNS) of the Navier-Stokes (NS) equations require sufficient spatial and temporal resolution so all scales of motion are resolved in space and in time. De- spite increases in the power of modern computers, turbulent flows present an extremely large computational burden that may not even be achievable given the particulars of the fluid state. Augmenting simulations with turbulence models can relax the mesh constraints to reduce the time and computational cost, but without an exact solution for comparison the validity of the results is unknown. Therefore it is essential for simulations to be corroborated with empirical data, such as the velocity fields from particle image velocimetry (PIV), among other methods. 2 Figure 1: Simplified wing planform showing common motion types for flapping-wing experi- ments in three dimension space, which may or may not have a mean or surrounding flow, u. (Left): Plunging; reciprocating linear translation (Middle): Oscillating root-flapping (Right): Os- cillating pitching Improvements in experimental technology as well as computational power means that it is now possible to resolve time-varying, 3D velocity fields of unsteady flow at high spatial resolution. Biological flapping-wing flight involves complex three-dimensional motion and is experimen- tally studied in varying degrees of bio-inspired simplification. From a 2D perspective, the wing motion can be described as a combination of heaving (the translation of the wing within the flow field) and pitching (the oscillating rotation relative to the leading edge). The extension into three dimensions can be split into two branches of research; whether the heaving is pure translation normal to the wing planform or instead rotating about a central sweep axis, or root-flapping, either at the wing root or at an offset. The motion types are shown individually in Figure 1. These topics can be further divided into two sub-categories based on type of flight, depending on whether it is forward flight or hovering, which governs the direction (and existence) of the surrounding flow velocity, u. The governing equations of the fluid are the incompressible 3D Navier-Stokes (NS) equations and the continuity condition, expressed in vector form as ρ ∂u ∂t +(u· ∇)u =− ∇p+µ∇ 2 u ∇· u= 0, (1) 3 where u(x,y,z,t) is the 3D velocity vector, ∇p is the pressure gradient, ρ is the fluid density andµ is the dynamic viscosity coefficient. Substituting non-dimensional variables (denoted with superscript *) into the momentum equation in 1(a), ∂u ∗ ∂t ∗ +(u ∗ · ∇ ∗ )u ∗ =− ∇ ∗ p ∗ + 1 Re ∇ ∗ 2 u ∗ Re= ρU ref L µ (2) the fluid state is shown to be characterized by a single parameter, the Reynolds number Re, given a reference velocityU ref . The relative magnitude of the viscous forces in the system defined by Equation 1 scales with Re, and when Re is the same, then two flow systems which may otherwise differ in velocity and length scales, the systems are dynamically similar. When the momentum balance has an acceleration forcing term such as sinusoidal flapping at a single frequency ω, the non-dimensional form of Equation 1 becomes k π ∂u ∗ ∂t ∗ + u ∗ · ∇ ∗ u ∗ =− ∇ ∗ p ∗ + 1 Re ∇ ∗ 2 u ∗ . (3) This introduces a new dimensionless parameter, the reduced frequency k = ωc/2U [17, 18], where c is the representative wing chord length and U is the reference velocity. A classical measure of the unsteadiness of the flow, the reduced frequency is defined as the ratio between the spatial wavelength of the flow disturbance and the mean half-chord of the wing [19]. The construction of this parameter assumes small amplitude motion, and in fact does not include any measure of amplitude, A at all. For flapping-wing systems that have peak-to-peak stroke angles of up to 120°, it is reasonable to look for another dimensionless parameter that takes this linear travel distance into consideration. For hovering flight with zero mean flow, a Strouhal number St= fA/U can be calculated with the mean wingtip velocity in the denominator. In [20], flapping systems found in nature are found to lie in the range 0.2< St< 0.4, a range in which the FW-MAV in [21] operates. 4 Though high-complexity flows can be computed and measured, it may be unclear how to translate that information into practical flying machines. A fundamental problem is supporting the body weight of the vehicle; the majority of current designs have maximum aerodynamic lift on the order of the body weight without an on-board power supply [22–25]. Since real flapping- wing fliers in nature carry power supply and cargo, this suggests that there is a lot of room for improvement. There are lessons to be learned from the fluid flow in order to increase perfor- mance, but the task of discerning the relevant information from high-resolution experiments and simulations is not straightforward. Existing force data from in-situ and scaled experiments warrant further study of the aero- dynamic forces produced during the flapping-and-pitching motion [21], especially during the initial transient period where the flow begins to develop. A typical micro-robotic wing assem- bly has a passive pitching hinge that results in a variable angle of attack during the flapping motion. This hinge has been extensively studied and optimized to deliver favorable pitching an- gles, while the wing planform itself has been largely assumed to be rigid. Through the use of directional (anisotropic) carbon fiber composites, the wings were designed to be as light as pos- sible while remaining stiff enough in the spanwise direction to maintain the bandwidth to flap at frequencies up to 120 Hz without severe flexing. The degree of deformation along this axis was generally left unquantified, only visually verified through high-speed image capture that the wings achieved the full stroke angle without attenuation. In practice, two major points were observed through this design process; seemingly small human fabrication errors could make the difference between achieving liftoff/flight versus a failed prototype, and high-speed videography of wing experiments showed inconsistent and significant magnitude out-of-plane deformation of the wing surface during motion. Therefore the stiffness of the wing planform itself was directly observed to be a performance parameter of interest, and current literature supports the idea that the bending stiffness of the anisotropic wing design used in [21] likely has a significant effect on the wake structure and lift-increasing vorticity[17, 26]. 5 The effect of wing flexibility on the aerodynamic lift and propulsive efficiency of flapping- wing systems has been a topic of increasing interest in the research community, as improving experimental and computational systems are well-suited to investigate the complex 3D problem as well as the continuing advances of FW-MAV research. The flapping system has been demon- strated to be sensitive to Re [17], aspect ratio AR [27], and frequency ratio f flap /f n [28]. The fundamental mode of vibration f 1 can be modeled using Euler-Bernoulli beam theory for the free vibration of an undamped 1-D cantilever [29] as f n=1 = ω n 2π = β 2 2πL 2 s EI ρA , β n=1 = 1.8751 (4) whereL is the beam length with cross-sectional areaA, moment of inertiaI, densityρ and elastic modulus E. Overall, there is evidence that flexible wings have the potential to improve perfor- mance over their rigid counterparts [30–35]. There is general agreement in the literature that chordwise flexible wings that produce negative camber reduce the wing frontal area for high angles of attack, streamlining the body and decreasing the mean drag coefficient [17, 26, 28, 36, 37], but if and how the lift force is increased for these same systems has not been definitively an- swered. Research such as [30, 38] suggest that critical stiffness ranges exist where flexible wings produce higher lift as well as reducing drag, which would increase the performance of a FW-MAV without adding weight or complexity. Chordwise flexibility was selected as the parameter of interest in this work as the spanwise flexibility is highly coupled with the second-order dynamics of the wing flapping at scale. There- fore we consider the dimensionless rigidity parameter, effective stiffness Π 1 , to quantify the de- gree of aeroelastic bending. This variable is a representation of the ratio of elastic bending to the aerodynamic fluid forces acting on the wing and equal to D s /ρU 2 L 3 , where D s is the bending 6 stiffness coefficient of an isotropic plate. This parameter appears in the dimensionless equation for the forced out-of-plane motion of an isotropic thin plate[17] as Π 1 ∂ 4 w ∗ ∂x ∗ 4 + 2 ∂ 4 w ∗ ∂x ∗ 2 ∂y ∗ 2 + ∂ 4 w ∗ ∂y ∗ 4 − ρ ∗ h ∗ k π 2 ∂ 2 w ∗ ∂t ∗ 2 = f ∗ , (5) where w is the displacement normal to the plate of a point on the mid-surface, given the dis- tributed transverse load f ∗ , chord-normalized plate thickness h ∗ = h /c m , and structure-to-fluid density ratioρ ∗ = ρ s/ρ f . The same reduced frequency term,k, from Equation 2 appears here due to the forcing from sinusoidal flapping. As payload capacity is limited by the capacity for lift generation, a hovering flight mode is studied in this work, with no mean flow around the flapping-wing mechanism. The flexing thin plate wing follows a prescribed flapping Φ(t) and pitchingΨ(t) function for simple sinusoidal motion, Φ(t)=Φ A sin(2πft+δ φ ) (6) Ψ(t)=Ψ A sin(2πft), (7) whereδ φ = 90 ◦ . This implies that the stroke angle starts at the maximum valueΦ A and the pitch angle starts at zero. Measured over the flapping cycle, the force/torque data are used as in [39] to estimate the coefficients of lift and drag ( C L andC D , respectively.) The resultant force and wake data are compared for both the first and second cycle of periodic motion to determine the effect of wing flexibility on the propulsive power of the low-Reynolds-number flapping-and-pitching wing aerodynamic system, and propose aΠ 1 value that engenders a favorable wake and aerodynamic performance for FW-MAVs. 7 0.1 ResearchObjectives Like the scope of micro-robotic research in general, this thesis covers a wide range of topics from the details of microfabrication and force sensing to flow visualization and vortex identification. The reader has thus far been introduced to the background motivations and physical principles that govern the sub-gram FW-MAV system under investigation. The main objective of this work is to improve the design and performance of flapping-wing microrobots; to do this there must be a clear understanding of how the components are made as well as how their performance is modeled and quantified. Critical subsystems such as the flexure transmission and piezoelec- tric actuation are addressed in detail, as well as underlying fabrication methods such as carbon fiber layup design and curing processes. Even with the knowledge and capability to create these robots, there is no commercially-available dynamic force sensor that is able to measure the forces produced over the periodic flapping motion. The cycle-mean aerodynamic lift is on the order of several milliNewtons at a flapping frequency of at least 100 Hz, therefore a high-resolution, high- bandwidth device must be used. Therefore this work also provides a framework for the design and fabrication of a dynamic microNewton-resolution force sensor that can be used to measure the force output of the microrobotic flapping wings as well as individual components such as the actuators. Inherently linked to the main objective of FW-MAV performance improvement, this work also aims to investigate the fluid dynamic phenomena behind the performance that are inaccessible at scale. Fluid velocity field measurements on dynamically-scaled models allow for flow visual- ization and analysis that can provide insight on the fluid-structure interaction that occurs during flapping-wing flight. By studying a hovering flight condition, the focus is on the wings capacity for lift generation, which is currently a critical limitation of the devices. This research aims to guide mechanical design toward practical regions in the wing parameter space for more powerful (higher lift) and/or efficient (lower drag) micro-robotic wings. 8 Chapter1 Design&FabricationoftheFlapping-WingMicroAir Vehicle The FW-MAV is a complex dynamical system involving mechanical actuation, linkages, compos- ite materials, and electrical systems. Most if not all of the parts (depending on the design) are laser cut from specially-made multi-material composites and assembled by hand with the visual aid of a stereo-microscope. For the numerous individual components there are physical phenomena to study as well as fabrication problems to solve and optimize, which requires a stable experimental platform. An established fabrication and characterization process of microrobotic subsystems is as crucial to the success of the research as the quality of the experimental measurements of the FW-MAV in operation. Not only does it create a baseline comparison from which to quantify performance improvement, it allows for engineering problems to be solved without building an entire vehicle, reducing build time and mechanical coupling of components. This chapter describes the development of designs and procedures that produce repeatable experimental results which have been verified to closely match the magnitude of lift produc- tion in [3]. Critical design parameters such as carbon fiber fabrication and piezoelectric actuator design are discussed in detail to guide the reader towards an understanding of the minimum re- quirements for FW-MAV flight. Seemingly smaller details and assumptions that are not normally included in peer-reviewed literature can be defining (and yet hidden) characteristics of successful 9 Figure 1.1: The 70-mg Double-Actuator FW-MAV created at the Autonomous Microrobotic Sys- tems Laboratory (AMSL) at the University of Southern California. designs. This section attempts to thoroughly address the microfabrication process and highlight the common points of failure. 1.1 TheFlapping-WingMicroroboticInsect The microrobotic insect is comprised of actuator, transmission, and wing assemblies held together by a simple carbon fiber frame. Depending on the particular design and number of actuators, these carbon-fiber-based micro-robots weigh between 70 and 90 mg. Sinusoidal motion inputs generated by piezoelectric cantilever actuator deflection drive a flexure transmission which flaps the wings with periodic motion. This flexure transmission replaces “traditional" rotation- and gear-based systems that have bearings and axles as as their weight and frictional forces are pro- hibitively large for the system scale. The output of the transmission drives the leading edge of the wing and a passive pitching hinge allows the wing planform to rotate about the leading edge and experience an induced angle of attack over the flapping cycle. This reciprocating rotation generates aerodynamic lift and allows the robot to overcome the gravitational force of the body weight and achieve flight. 10 Figure 1.2: As shown in the detail view, the central carbon fiber is unidirectional, with the fibers aligned along the longitudinal direction of the cantilever. The appropriate combination of input signals to the three separate layers gives the ability to control the direction and magnitude of the cantilever deflection. A piezoelectric (PZT-5H4E) ceramic bimorph cantilever actuator is used for its fast dynamic response, high resolution, and relatively large deflections for the scale of the system. The trape- zoidal shape is a modified cantilever to maximize deflection while minimizing stress in the beam. The peizoelectric ceramic plates are cured to a central carbon fiber layer, using coplanar SiO 2 segments (shown in white in Figure 1.2) as an electrically inert base and tip material. Extra car- bon fiber is cured on the surface to reinforce the seam between the disparate materials, as well as acting as an adhesive for the custom copper-clad FR4 electrical contact pads for the power tether. Details of this process are explained in the fabrication section in [40]. When the two plates are held at a voltage difference, an input signal sent to the central uni- directional carbon fiber spine splits the ratio of the voltage between the two outer faces of the cantilever. The carbon fiber itself is highly conductive but has a significant interface resistiv- ity, allowing the gap in potential. If a sinusoidal voltage signal is sent to the carbon fiber, when the top plate is at V H the cantilever will deflect +δ max and when the bottom plate is at V L the cantilever will deflect − δ max . 11 Table 1.1: Definitions of Actuator Equations and Variables. Subscripts p and s refer to the piezo- electric layers and the support layer, respectively. V - Applied voltage ρ - Density L - Cantilever length d 31 - Transverse piezo coeff. w - Cantilever width A= E s/E p - Modulus ratio t - thickness B= t s/2t p - Thickness ratio E - Young’s Modulus C= ρ s/ρ p - Density ratio Extensive research on the design and optimization of the PZT actuators can be found in [41]. The piezoelectric bimorph has been thoroughly characterized in terms of expected free-end de- flection, blocked force, and maximum frequency [42] which in theory provides enough infor- mation for a model to predict the necessary size to fully actuate the wing system at the desired frequency. Given the variables in Table 1.1, the free-end deflection is expressed as δ = 3d 31 L 2 V 2(2t p +t s )t p · (B+ 1)(2B+ 1) AB 3 + 3B 2 + 3B+ 1 , (1.1) the blocked force at the tip is computed as F b = 3d 31 E p w(2t p +t s ) 2 V 8Lt p · 2B+ 1 (B+ 1) 2 , (1.2) and the maximum frequency of operation is f max = 3.52(2t p +t s ) 4πL 2 s E p 3ρ p s 1+ 3(2B+ 1) 2 + 4AB 3 4(B+ 1) 2 (BC+ 1) . (1.3) However, as the actuator is over half of the total weight and the best-case lift margin of the FW- MAV is small, the volume minimization problem prevents the use of an over-sized actuator to ensure performance at all frequencies. Further explained in Section 2.1.3.1, under-performance issues such as deflection attenuation due to incorrect cured thickness of the carbon fiber and low tensile modulus of the material are common but and are difficult to diagnose in situ. The robot must flap at frequencies at and above 100 Hz to generate enough lift to take off, so this bandwidth issue is critical in every step of the design. 12 Figure 1.3: Diagram of the kinematic system, showing the relation of input displacement from the actuator, δ, to the output wing angle, θ wing . The detail view shows the composite flexure hinge made from Kapton ® (yellow) which is the micro-scale approximation of a linkage pin. Figure 1.4: Transmission Design. (a) Isolated from the chassis (as it varies from project to project), the CAD render of the transmission shows the carbon fiber links in black, connected by the yellow flexure hinges. (b) The ‘shoulder’ of the wing is made from the same composite as the transmission, to create the passive rotational wing. The wing subassembly is glued onto a tab on the transmission to complete the assembly. 13 Figure 1.5: (a) The flapping motion is a rotation about the ˆ z axis (b) Passive pitching occurs from the periodic oscillation of the flapping wing, causing rotation about the ˆ x axis. The microrobotic flexure transmission is shown in Figure 1.4 which transforms the actuator input into wing motion. The tip of the actuator in Figure 1.2 has a tab to engage the slot in the final link of the transmission (Figure 1.4). The inverse kinematics have been studied in [43, 44] and the output-input relationship isθ wing ≊ 1 /L 3 · δ for smallθ wing andδ. Thistransmissionratio can also be approximated as in [24] as a simple lever mechanism with effort arm L 3 and output arm R cp , which corresponds to the center of pressure on the wing. The other end of the linkage is rigidly fixed in the chassis of the robot, which causes the square top link to rock back and forth with the oscillating input from the cantilever. The design was developed in [43–45] and used extensively in projects such as [1] and [24]; the linkage is made from 100µm flat carbon fiber bonded together with segments of 12.5µm Kapton ® polyimide (DuPont de Nemours Inc.) which function asflexure hinges. These hinges allow the structure to move with approximately equivalent kinematics to a rod-and-pin mechanism with bearings. While the beam deflection input is not linear, the off-axis motion is damped by the polymer hinges and no significant problems have been identified. During flapping (Figure 1.5), an appropriate pitching angle (of attack) (Figure 1.5(b)) of the wings generates a net positive aerodynamic lift force per period. During this passive rotation, the pitching angle is determined by the moment of inertia of the hinge segment and the magnitude and location of the center of pressure on the wing. Hinge dimensional parameters require signifi- cant iteration to tune given that the force on the wing can vary based on robot fabrication quality and age of the material. A single microrobotic bee can take days to fabricate from start to finish, 14 Figure 1.6: Modeling of the Actuator-Transmission-Wing System. (a) Model of the full second- order system with coordinatesδ andθ. (b) Reduced model with equivalent stiffness, mass, and damping, as well as a reduced coordinate system. so this presents an obvious need for a testbed that can have easily-interchangeable parts for rapid iteration. However, loose connections cause unwanted body vibrations and loss of efficiency, and many designs that would maintain the fit while being interchanged do not work at this physical scale or do not hold up to repeated use. These three main mechanical components of the robot, the actuator, transmission and wings are used to formulate the lumped dynamic model of the system. These components exhibit non- linear inertial, geometrical, damping, and elastic behavior, but the system proposed in [25] and shown in Figure 1.6 effectively linearized these components and validated the lumped model through system identification. As in [25], to use a linearized model for the PZT actuator we must assume that the voltage source is ideal (directly translatable/equivalent to some force in the mechanical domain) and that the resonant frequencies of the actuator are not near operating conditions. Figure 1.6(a) shows the complete model of the three components; the forced mass- spring-damper actuator model with position state variable x that is the input to the mechanical transmission T , outputs the wing angle displacementφ into the damped rotational dynamics of the loaded wing. 15 Figure 1.7: Dynamic Model of the System. Given that there is an approximate kinematic relationship between input displacementδ and output angle θ wing , there is only a single degree of freedom in the system. The effect of trans- mission dynamics on the system can be incorporated into the lumped actuator model in Figure 1.6(a), mapped to the actuator reference frame with the transmission ratio. Thus, Figure 1.6(a) can be simplified into the system shown in 1.6(b) where a forced, second-order damped model of the form m eq d 2 δ dt 2 +b eq dδ dt +k eq δ =F a (1.4) with equivalent massm eq , stiffness k eq and damping: b eq m eq =m a + T 2 J φ b eq = T 2 r cp b k eq =k a + T 2 k t . (1.5) The mass and stiffness of the actuator are m a ,k a , respectively,J φ is the rotational inertial moment of the wing,r cp is the distance to the center of pressure on the wing (geometry-dependent),b is the damping from [25], andk t is the stiffness of the transmission. From this second-order system, an equivalent dynamic model can be constructed (Figure 1.7) to relate the input deflection to outputforce. With a linearized system model and feedback (e.g. a sensor to measure the actuator 16 deflection), it is possible to develop controllers for autonomous flight patterns. In the frequency domain, the transfer function is thus T(s)= ˆ δ(s) /f(s)= 1 m eq s 2 +b eq s+k eq . (1.6) 1.2 Fabrication Sub-gram flapping-wing Micro-Air Vehicles (FWMAVs) are created from a range of high-performance engineering materials, selected for characteristics such as strength, weight, elasticity, electrical conductivity, etc. This section details the individual materials, the process of microrobotic com- posite fabrication and assembly, and sub-gram dynamic micro-actuation via piezoelectric ceramic composites. 1.2.1 CarbonFiberFabrication The structural components of the microrobots are made from carbon fiber laminated composites, chosen for the high strength-to-weight ratio. In raw form, bundles of carbon filaments are pressed into a flat “tape," infused with a curing adhesive media, which is called unidirectional carbon fiber prepreg. Prepreg made from High-Modulus (HM) carbon fiber has a tensile modulus in the range [400:600] GPa. Cured under temperature and pressure profiles in an autoclave or vacuum-bag environment, these layers are stacked in sequences of fiber orientations to create laminate sheets. As the tensile strength is along the fiber direction, depending on the arrangement of the layers, the final material can be designed to have particular strengths in certain axial dimensions, as shown in Figure 1.8. The layups are described using degrees of rotation from the first vertically-oriented fiber layer, the 0 ◦ ply, depicted on the plot in Figure 1.8(a). All of our carbon fiber microrobot parts are cut with the laser fromflat sheets, but the stacking sequence varies based on the function of the component. Basic cross-plies, only containing pairs of (-0-90-) ◦ layers, creates an orthotropic material that is well suited for general structural components. Initially, the balanced layup was 17 made 0-90-90-0, but it is quite thick (> 100 µm) and difficult to integrate with the dimensional constraints of the transmission hinges. As experiments progressed and produced more repeatable results, we began to experiment with other layups to improve performance. Although not symmetric about a central plane, the thinner layups 0-90-0 and 0-45-0 are symmetric sequences and will cure without planar curvature but do not have equal modulus in the two fiber directions. This in fact turned out to be useful in diagnosing an early problem with wing span failure, as we were able to eliminate the problem by aligning the 0-direction with the length of the span, providing a strong leading edge and a central − 45 ◦ -ply layer aligned with the wing spars which experience smaller forces. The topic of planar strain and its relationship to curvature (presented in [46]) based on an applied load is fascinating and has many potential possibilities for laminates for microrobotic designs. The designs are restricted to largely rectangular geometry and curved carbon fiber surfaces (not edges) are not achievable with current methods, but from some serendipitous accidents during the learning process, we found that extremely thin asymmetrical layups made from 15 gsm (grams per square meter) UD fiber are flexible and can ‘snap’ back and forth between two stable positions with opposing planar curvature. The use of these layups to engineer bi-stable manifolds is a topic with great promise. The laminates are cured with pressure and temperature to compact the layers and activate the curing of the bonding media. For our application we use epoxy resin-based, aerospace-grade car- bon fiber prepreg for the temperature resistant qualities, high modulus, and low density. Sourcing this material in small quantities is a difficult task; prepregs from Intermediate Modulus (IM) to High Modulus are usually special order in large production runs, and the Ultra-High-Modulus (UHM) prepreg is even more difficult to find (after five years I still have not been able to get any). In addition, the pre-impregnated high-T g epoxy thermoset resin has a shelf life of maximum a few years even when stored in the ideal environment of temperatureT <− 20 ◦ C, therefore even beyond the cost of a full production run of the material, it is not rated for long-term storage. As 18 Figure 1.8: Carbon fiber composite layups: layers of unidirectional carbon fiber prepreg are placed in a “sandwich" orientation and cured to result in a sheet composite with tensile properties in the desired direction. For a composite that lays flat in equilibrium, such as is desired for this work, the layups must have a symmetrical composition in the thickness direction. Figure content adapted from [47]. we learned through experience with various other unidirectional fibers and cure cycles, the mod- ulus of the material is most important for the actuator performance. After switching the carbon fiber in the central layer of the actuators from Toho Tenax Q-A112 prepreg (IM fiber, modulus 275 GPa) to Royal Tencate M40J-based prepreg (HM fiber, modulus 495 GPa) in January 2018, successful jumping/takeoff tests of the FW-MAV went from the outlier to the norm. The physics of this phenomena is discussed with the actuator model. The curing process for these epoxy-based composites is shown in Figure 1.9, for which we adjust the maximum set temperature and ramp times depending on the glass transition curve of the epoxy matrix (between 120 ◦ C and 130 ◦ C for the prepreg used here). This plot is roughly ap- proximated from specification sheets, training instructions, and literature such as [48] and [49]. Designed to control the modulus, compaction, and void content, this process exploits the depen- dency of epoxy viscosity on temperature and pressure. When the viscosity becomes theoreti- cally infinite and the epoxy flow around the fiber bundles is arrested, the formation of branched molecules in the liquid commences, called gelation. Continued exothermic reactions from curing drive the temperature higher and cause vitrification, the increased tension on the surface causes the plane to become smooth and characteristically glassy. 19 Figure 1.9: General curing process recipe for out-of-autoclave carbon fiber prepreg. Temperature ramp initialized att 1 . At timet 2 , the ramp continues until reaching the cure temperature and gel point neart 3 . The vitrification point occurs at t 4 and the temperature dwells until the ramp down att 5 . Crossing the T g -transition at other points results in incomplete gelation and adverse results on the mechanical properties of the final composite material. The actual temperature reaches the set temperature shortly after the gel point of the resin and the flow around the fiber bundles is ar- rested. Seven distinct phases of the process are shown in Figure 1.9, which we initially performed in an Elastomeric Vacuum-Bagging Tool (EVT) that I helped design (with Torr Technologies, Inc.) until we purchased a hydraulic heated platen press in 2016 to scale up the production. From room temperature, pressure is applied to start compacting the layer and remains constant. In Phase I, the pressure is immediately ramped to the maximum value at which point Phase II starts at time t 1 , and a temperature ramp is initialized and dwells (Phase III) at a mid-range temperature to allow for continued compaction while the viscosity and temperature rise. At timet 2 , the ramp continues (Phase IV) at a fixed rate until reaching the cure temperature. At t 3 the resin matrix reaches its gel point. The vitrification point happens at t 4 when the resin passes the glass transition threshold and becomes an amorphous solid. The length of dwell time between t 4 and the initialization of the ramp back down to room pressure,t 5 , depends on the desired amount of curing of the resin for properties such as tensile strength, etc. 20 Figure 1.10: Overview of the SCM process flow. (Left): Individual layers are cut from structural and adhesive material. (Center): Composites are cured with heat and pressure in an elastomeric vacuum-bagging tool (EVT) and oven or in a hydraulic platen press. (Right): Cured SCM com- posites are release-cut with the laser micro-cutter. Figure 1.11: The Single-Actuator Microrobotic Bee. (a) (Clockwise from top) Wing composite sheet, airframe composite sheet, transmission composite sheet. (b) Released components from airframe composite sheet and the trapezoidal bimorph actuator. (c) Released transmission linkage and the folding process. Tabs are left on the initial part to help with the manual folding step and are cut off to have a compact final form. (d) Assembled and wired 70-mg microrobot 21 1.2.2 CompositeFabricationandtheSCMMethod The Smart Composite Microstructures method was pioneered in [50] and continued with the development of the original microrobotic research platform in [45]. A stiff structure is laminated to a flexible layer, which creates a network of rigid members connected by flexible joints. Early work with simple hinges and linkage designs have grown into an entire research subject, with multiple university labs across the country using this method to produce specialized mechanical motion at both the micro- and macro-scale. From research such as [1] and [3], we were able to emulate the designs and methods and replicate the flight experiment platform. Of course this is more simply said than done, and hundreds of microrobotic failures before success highlights the many minute imperfections that can cause a large effect on the system output. For SCM with [60 : 130]µm carbon fiber, we use 7.5µm thin-film sheet adhesive (DuPont FR1500) to bond to various thicknesses of Kapton polyimide sheets (shown in Figure 1.12). Sum- marized in Figure 1.10, the fabrication process starts with cutting individual layers for the com- posite with our laser micro-cutting system with a 10-µm beam spot size. They are then aligned on a jig, and then cured at high temperature temperature (130 ◦ C) and moderate pressure (200 psi) in a machine such as a platen press or out-of-autoclave vacuum-bagging, and finally release-cut from the composite sheet. Individual sheets are aligned on a flat base with a pin alignment system to fix the layers in the correct position during curing. The cured composite is laser micro-machined a second time to release a single piece with contours generated by the individual layers such as flexible connections between subcomponents. These connections, shown as the yellow Kapton segments in Figure 1.12 and Figure 1.13 allow constrained rotation between the rigid carbon fiber links. The SCM method is incredibly precise for components and linked subsystems, but subsequent manual steps can destroy the results of careful planning and design. Even for the 12 simple Single- Actuator Microrobotic Bee shown in Figure 1.11(b-d), the accumulated fabrication error was so difficult that the design was ultimately abandoned. Simple linkages are moderately straightfor- ward to design and fabricate, but with increasing complexity (as in progressing to the 3-DOF, 22 Figure 1.12: SCM-Type Composites used for Microrobotic Designs. Views are cross-sections of the structures. (a) Piezoelectric ceramic cantilever, bonded to the carbon fiber spine during the curing process of the raw carbon fiber. (b) Basic linkage design used for transmission, etc, is constructed from two sheets of cured 0-90-0 carbon fiber bonded to either side of 12.5µm Kapton polymer with thin-film sheet adhesive. (c) Simpler and lighter components like the wings require only a single piece of carbon fiber which can be bonded (with the sheet adhesive) to a polymer like Mylar on only one surface. Figure 1.13: The Smart Composite Manufacturing (SEM) method. (a) Individual layers aligned on flat base with a pin alignment system (b) Cured composite release cut (c) Example of a (now obsolete) transmission linkage created from this process. (d) The linkage is folded and glued into the proper form for integration into a larger frame. 23 dual-actuator design) comes increasing assembly difficulty. The components of the microrobots are manually assembled with tweezers using a stereomicroscope, and as the number of parts rises, so does the possibility of error at each step. Additionally, structures such as the wings and actu- ators are incredibly fragile but they must be precisely positioned or the robot will not function correctly. The overall scale of the microrobotic parts makes visual inspection or quality control in general difficult and the underlying reason for low performance is not easily identifiable. Ideally the components could be studied independently, but the complex dynamic interactions between the components in complete and loaded system can make simplified, isolated experiments not inherently informative about the fully-assembled robot’s dynamics. Therefore the SCM method has been adapted to improve assembly tasks by incorporating limit-stops and simplifying larger substructures into flattened composites that can be later release-cut with the laser in a fully- assembled state. 1.2.3 AdvancedCarbonFiberCompositeFabricationTechniques While the original SCM method included manual folding and glueing of the transmission links into the proper position, the internal 90 ◦ -angles on linkL 3 in Figure 1.3 were extremely difficult to properly set and caused kinematic and dynamic inconsistencies. To make the process both more internally consistent (same person, different iteration) as well as externally consistent (person- to-person), limit-stops were added to the folding connections. Inspired by the folded designs in [4], we introduced “teeth" (as shown in Figure 1.14) in the individual layer files; this pattern creates rectangular cavities along the inside angle of the folding joint and prevents the part from over-rotating. To be clear, when the depth of the tooth pattern is equal to the thickness of the carbon fiber, the walls prevent over-rotation of the links As was done previously, this connection is carefully reinforced with cyanoacrylate adhesive (CA) (shown in Figure 1.15), applied with custom-made micro-tools such as a 52 gage wire clamped in an X-Acto blade handle. Beyond prescribing nominal angles between the members in the linkage, the teeth/folding method is a useful tool to create structures that can be built quickly and easily, with higher rates of success in 24 Figure 1.14: Pin holes and alignment marks for the second cut process are shown in white. Each color represents a layer in the layup. Detail (a) shows the transmission linkage that is cured concurrently with the frame. Detail (b) shows the tooth pattern for the limit-stopped folding process. Figure 1.15: Folding Method vs 5-Layer Composite for Transmission Fabrication. (a) Folding method to create 90-degree corners in the linkage is highly dependent on human ability. (b) For the 5-layer transmission, proper linkage angles and dimensions come from the alignment quality, design and dimensions of the cuts. terms of dimensions and performance. While this three-layer composite is too heavy to make a chassis for the flying robots, it is excellent for making experimental setup components and other complex subassemblies or jigs. The limit-stopped folding technique was proven to create sufficiently constrained designs to achieve liftoff in many jumping tests for both the single-actuator bee as well as the double- actuator bee, but the method is time-consuming and requires extensive skill and experience. I decided to try a proposed design improvement from [24] to simplify the linkage to be made from 25 Figure 1.16: Five-Layer Transmission Layup (No Folding). (a) Initial prototypes of the five-layer transmission that removes the manual folding step from the assembly process, shown partially released. (b) Trapezoidal transmission (left) and rectangular transmission (right). a single composite made from five layers (CF-Kapton-CF-CF-CF-Kapton-CF) instead of three (CF- Kapton-CF). Wary of previous mistakes in terms of trying new designs (that were so different that all of the body dynamics changed) and motivated by the need to create a functional, parts- interchangeable testbed, I did my best to match this 5-layer transmission to the microrobotic insect design that we have proven to work. After testing a purely rectangular linkage (constant cross-sectional area of the links with respect to their length), the shape was changed to the slightly trapezoidal profile introduced in [4] (composite shown in Figure 1.16, distinct designs shown in 1.16(b)). The release-cut laser recipes took many iterations to identify, to avoid any interfering cuts or burns on the Kapton layer. This adjustment nominally aligns the intersection of the two neutral axes of the cantilever actuators over the center of gravity of the chassis and decreases torsion on the flexure hinges in the transmission, and the final dynamic motion is smooth. The success of the 5-layer linkage de- sign provided motivation to create a “unibody" testbed, which is in fact the design that is shown in Figure 1.14. The same base layers as the lower transmission triplet (CF-Kapton-CF) are em- ployed to make a 3-layer foldable frame that is pre-attached to the 5-layer linkage. This structure is shown in various stages from composite to parts to assembly in Figure 1.18. All components of the chassis and transmission are part of this structure except two cross-body braces and the “head." Assembly time of these parts is on the order of minutes with higher dimensional precision 26 Figure 1.17: The CAD model shows an example of a structure designed for testing and character- ization of components. Figure 1.18: Unibody Design. (a) Two unibodies, as released on the right and partially folded on the left. (b) Detail view of the transmission, which is folded into position and secured by the head, the folding teeth, and CA adhesive. (c) Unreleased composite. (d) Partially released composite and partially assembled unibodies. than anything achieved manually with the previous methods. The narrow body creates an access point to the rear face of the actuator which is normally occluded in the Double-Actuator Bee, allowing repairs and rewiring to be done in-place, and the tongue-and-groove connection to the alumina base precisely locates the height of the neutral center axis of the cantilever with respect to the engagement feature on the transmission. The CAD render of the unibody design with the actuator and wing sub-assemblies attached is shown in Figure 1.17. While I have presented the current iteration of this design, as always there is both a history as well as continual room for improvement. More than five significantly different structural designs 27 were prototyped for this setup and failed due to collisions, tolerances, fabrication errors, etc. The design process unfortunately has many steps of 3D modeling, planning the folding mechanisms and sequencing manual steps of the process. The surface profiles are made into the folding pat- terns by connecting the profiles in a 2D drafting software (AutoCAD) which are then exported to the laser marking software. The lack of automation in this process and the use of multiple soft- ware modules naturally results in mistakes. Tolerances on materials such as the flatness of the alumina that are considered normal (say± 0.01µm) can mean the difference between a good fab- rication session and something that ends in the trash. In terms of improvement, initial attempts to make a simplified wing connection feature were unsuccessful and for a useful experimental setup, this issue must be solved. In general, planning to use a liquid adhesive on a location near the transmission carries a large risk of unintentional spreading and irreversible freezing of the transmission. However, mechanical mates are also unfavorable because they add weight and have a very small region to connect to, given that the transmission is only 1.5 mm wide in the direction of the longest dimensions of the link. The current tab geometry is mimicked from the feature on the Double-Actuator Bee and is known to degrade with time and repeated use, but we know that it works. 1.2.4 BimorphPiezoelectricCantileverActuators Although the actuator design has been previously discussed, the fabrication process is extremely complicated and time consuming and merits a short discussion. Based off the work in [51] and [40], we batch-fabricate multiple actuators at a time, curing a specialized composite (Figure 1.12) and cutting out the actuator shapes in the final step, depicted in Figure 1.19. Individual layers such as the PZT are cut with the laser to be assembled in the pin-based alignment stack that is also used for the general SCM process. A “spring-clip" jig fabricated from inert FR-4 holds the PZT in-plane with the alumina, sandwiching the uncured carbon fiber layer. When cured like normal carbon fiber prepreg, this layer bonds the ceramics together and creates a weak electrical connection without the addition of an actual adhesive. The cured composite is aligned to the laser 28 Figure 1.19: Batch Piezoelectric Cantilever Fabrication. (a) Laser-cut individual layers are inserted in the alignment stack. (b) A “spring-clip" jig holds the PZT in-plane with the alumina. (c) Second Alignment and release cuts (d) Manual release of actuators. again and cut halfway through, flipped over and realigned, and finally completely cut through and released from the surrounding material. The actuators are manually removed from the stack and are post-processed with 1200-grit polishing paper on the edges to prevent electrical shorts during operation. This process took almost four years to get to the level where our actuators performed similarly to the ones in [40]. Even starting with the base design of just the three-layer bimorph, the fragile piezoelectric ceramic was difficult to work with and kept free of cracks and defects during the fabrication process. Cutting, curing, and releasing these actuators can easily take two full days, at which point they need to be mechanically post-processed, wired, tested, and inserted into the robot. Between lack of knowledge of the system, overpressure/incorrect curing cycles and learning how to physically handle the material, nearly all of the early actuators that even survived into the testing phase were moderately damaged and often had highly unequal deformations in the two directions of motion. Although significantly more time consuming in the initial stage, we moved to the copper-clad design (again, from [40]) which requires a second cure cycle to bond the copper-clad FR4 electrical connection pads to the alumina base as well as creating mechanical ‘bridges’ to support the co- planar interfaces where the PZT meets the alumina. Seen in Figure 1.20(a), the high-bandwidth (wide base) actuator composite with copper leads connects all of the actuators in parallel for rapid testing and characterization. The Double-Actuator Bee uses a narrower cantilever as shown in Figure 1.20(b) and is shown split from the other half and mounted on a test frame in Figure 1.20(c). 29 Figure 1.20: Actuator Connections and Testing. (a) Addition of cladding layer to the PZT stack layup (shown partially released) The actuator produced by this stack is shown at the top of the image. (b) Smaller actuators for the split-actuator microrobotic bee, shown after the first half of the stack is cut. (c) Released and wired test setup to run 8 actuators synchronously. Conductive epoxy is used to bridge the small gap between the copper-clad FR4 and the surface of the PZT. This testing method is crucial for the success of the flying robots (as opposed to a simplified test system for a single wing), as the 3-DOF microrobot requires a matched pair of actuators. The same signal is sent to all actuators at once for a duration of a few minutes to test the material (for cracks, etc) and drive the actuators to their steady-state behavior for analysis. As I grew experienced with the laser alignment process and my cutting errors decreased, mechanical polishing was no longer necessary to prevent short circuits across the top and bot- tom PZT plates. Instead, we used the same edge treatment scheme used in the aforementioned benchmark research to use high-power laser pulses at the sharp edges of the cantilever to induce melting and remove the stress-concentrating geometry that causes cracking and failure in the material during operation. The refinement of this process in combination with the high modulus carbon fiber now produces significantly higher quality actuators with long lives, with matched motion in both directions and no loaded system attenuation at operating frequencies near 100 Hz. 30 Chapter2 Forcemeasurementsatmicro-roboticscale The following section contains the peer-reviewed journal article, “Clip–Brazing for the Design and Fabrication of millimeter-scale, micronewton-resolution Force Sensors"[21], published in edi- tion XX ofSmartMaterials&Structures from the Institute of Physics (IOP) Publishing. This mate- rial is followed by Section 2.2 which extends the use of the force sensor and examines additional data collected from FW-MAV flapping experiments. 2.1 Clip–BrazingfortheDesignandFabricationofmillimeter- scale,micronewton-resolutionForceSensors To move the frontiers of millimeter-scale robotics, high precision structures of the same scale are required. Current methods to create these types of structures include monolithic fabrication [52–55] robotic assembly [56], and 3D printing [57]. Each method has advantages and draw- backs, but a common theme among them is a high cost for ultra-precise structures of this scale, especially those made from robust materials such as metals. Here we present a new method for fabricating high-precision standalone structures on the mm–cm scale with µm-scale accuracy that maintain tolerances during fabrication and closely match simulated and numeric dynamic predictions. Interlocking metal parts are brazed together with wire clips to fuse the entire struc- ture in a single heating sequence, minimizing tolerance stack-up. Three-dimensional structures 31 Figure 2.1: The brazed dynamic micro-force sensing structure. can be designed through subcomponent moment-of-inertia exploitation to selectively transmit, amplify or attenuate forces in predetermined directions. Spanning many resolutions and ranges, micro-force measurement techniques vary based on application and are chosen for qualities such as time response and dynamic characteristics, load- ing conditions, et cetera. For the given resolution requirement of the time-varying force of in- terest in this work, produced by a sub-gram flapping-wing micro-air vehicle (MAV), there are various measurement methods presented in the literature that would be sufficient. However, to use a common example, amicro-electro-mechanicalsystems (MEMS) accelerometer could be used to measure total force, but the packaging problem is nontrivial and the sensor could not isolate the aerodynamic component. This decoupling is achievied by flapping the microrobot’s wings while it is rigidly fixed to a load transduction device, as the body dynamics are removed from the total force and the aerodynamic component remains as the sole input force. Quantifying the reaction of a load transduction structure to an applied force is performed with a transducer, which generally falls into one of two categories: contact and non-contact measurement. Contact methods such as fiber Bragg sensors [58], strain gages [59], and other piezoresistive transducers such as those in [60, 61] have high dynamic resolution but would add fabrication complexity, as they need to be fabricated and placed precisely at the sub-mm scale. 32 Electro-active transducers such as piezoelectric polymer or ceramic devices [62] can have higher resolution and large operating ranges for forces and frequency [63] but share the same mounting problem as the strain gages. Non-contact measurement techniques include optical sensing by laser interferometry [64] and computer vision [54] which are both difficult to integrate at this time and length scale without expensive equipment. A popular method of contactless measurement is the change in capacitance between a datum reference and a body under an applied force, which has been demonstrated with cantilevers [25, 65] and other more complex structures such as micro- comb drives [66]. This method has high time and spatial resolution and was chosen to be the measurement technique employed here. This work extends the concepts developed in [65, 67, 68] where a dual-cantilever, parallel- beam structure is designed to deform in a predefined manner to map deflection distance to ap- plied force. Clip-brazing increases ease of fabrication, reduces cost, and creates self-contained structures that are not fixed to a large base or substrate. We demonstrated the effectiveness of the method as a high-strength, high-precision micro-joining process by fabricating a dynamic micro-force sensor capable of micronewton resolution for a selectable range of frequencies of the periodic input force, that is validated by mechanical testing and electron microscopy. Section II discusses brazing for the fabrication of flexible mm-scale structures. Section III details the use and validation of the method to fabricate a micro-force sensor that is capable of measuring the dynamic forces produced by an 80 mg microrobot with a resolution of less than 1% of the body mass, or 7.85µN. To demonstrate the utility of this micro-measurement system, this work com- pares computational fluid dynamics (CFD) simulations for a sub-gram flapping-wing microrobot and the empirical force profiles measured in-situ. Sections IV and V contain the discussion and conclusions. 2.1.1 Brazingforthefabricationofflexiblemicro-structures Brazing is the creation of a metallic bond between two surfaces using a lower-melting-point filler metal. External energy is applied to melt the filler which fills the interface to form a connection 33 that is a combination of dissolving, alloying and diffusion reactions [69]. Brazing is performed at temperatures between 540°C and 1620°C [70], depending on the elemental composition of the filler, but always below the melting temperature of the metals to be joined. The wetting and spreading phenomena characterize the interaction of a liquid and a substrate. The wettability, defined by the contact angle of a liquid droplet on a smooth flat surface in a sessile droptest [71], is a measure of the chemically-driven equilibrium state of the three physical phases present (liquid, solid substrate and gaseous atmosphere). By Young’s equation, when wettability is high, the adhesive forces dominate the cohesive forces in the droplet and the molten material spreads over the surface [72], indicating the suitability of a given filler to create an intimate surface bond between two bodies. Spreading describes the amount and rate of transport of the filler material along the surface away from the source, the kinetics and equilibrium of which are driven by a number of complex factors such as the formation of inter-metallic compounds (IMCs) or phases [73, 74] and the chemical diffusion that occurs in advance of the leading edge of the spreading filler [75]. Brazing relies on capillary action, a transport phenomenon related to the wetting and spreading of the filler material and the topology of the interface. At molten temperatures of the filler, the intermolecular adhesion forces dominate the cohesion forces and produce a net force that draws the filler metal deep into microscale features in the connection [69]. This brazing phenomenon occurs predominantly in interfaces that have clearances in the range[0: 0.15] mm [76]. In general, from the subgroup of alloys that are known to wet the chosen component mate- rial, the brazing alloy is selected with consideration to the desired degree of potential chemical interaction and spreading as well as other properties such as strength and coefficient of thermal expansion (CTE). Brazing requires the use of a shieldingflux to chemically prepare the surface and to prevent oxidation from interfering with the bond formation in an open working atmosphere. This chemical cleaning agent is specifically chosen for use with a particular filler and host metal [69], as the compound should be at or near boiling around the melting temperature of the given filler metal. The high working temperatures required for brazing compared to those used in 34 Figure 2.2: Braze types. (a) Open corner joint braze: Starting from a folding-line, the sheet pulls open to make a ledge to rest the filler material on. When heated above the melting temperature, the filler wets the surface to create a clean fillet at the right angle. (b) Slotted tee-joint braze: Starting with two components, the pieces are transition-fit together in a slot (width equal to the sheet thickness), filler metal and flux is added and heated with a butane torch which draws the material into the interface. The result is a filleted structure with orthogonal faces. traditional soldering allow the use of high-melting-point alloy metals and result in strong and ductile bonds that can tolerate complex interface geometries [70, 77]. In some conditions, these connections may be stronger than the parent material [78, 79]itself. As brazing is a thermally- driven process, the required energy can be provided in many ways such as gas torches, furnaces, electrodes, lasers, induction, et cetera [69, 80]. Each method has its advantages; for instance, the precision of a laser or the batch capabilities of a furnace. 2.1.1.1 Methods The creation of precise brazed metallic mm-scale structures relies not only on the fabrication of individual components but also on the ability to make uniform and repeatable connections between them. The solution proposed here is to use interlocking parts, fully dimensionally spec- ified so that when properly assembled, the structure is constrained by design to be in the de- sired form. The high-precision features and connection points are laser micromachined from flat pieces and bent at perforated folding lines into a particular orientation as in [65], while also em- ploying slots like in [25] that hold the subcomponents together (as depicted in the first images of Figures 2.2(a,b), respectively). These connections are brazed together and made permanent. 35 The interlocking features function as a jig during the brazing process, preventing the assembly from shifting without any external support. We present an application of this method to create three-dimensional cuboid shapes constructed entirely with folded parts joined by the clip-brazing method introduced here. The ideal connections are modeled in a computer-aided design (CAD) program and the com- ponents cut with a laser micro-cutting system with a 10µm spot size. Adjusting for the kerf, this results in a graphic resolution of 5µm for features and components. The relatively square edges of the laser cut parts, near-nominal dimensions and precise locating slots allow a folded part to be manually inserted and held with zero degrees of freedom (DOF) via a transition fit. This man- ual assembly step has the most potential to introduce imperfections so careful inspection at high magnification is necessary to confirm that the specimen has not been damaged during handling. Initially the structure is held in place by friction and contact forces, so the folding points and the slots need to be reinforced to preserve the structural accuracy on the order of tens of microns. This objective is achieved with two types of silver brazed joints,opencornerjoint brazes andslottedteejoint brazes. These two braze types are shown in Figure 2.2. The open corner brazes in Figure 2.2(a) are used to reinforce the folding lines at prescribed unsupported angle changes in the structure where the material is weakened after folding. We employ a double-thickness base for the slotted tee-joint braze to anchor beams in a perpendicular orientation and produce a corner fillet seam, as presented in Figure 2.2(b). Thus the µm-scale accuracy is determined by the cutting resolution, near-zero positional tolerance on mating surfaces, and careful heating to allow the natural kinetics of brazing (due to the geometry of the design and the properties of the material) to occur without causing heat damage to the structure. The control of the brazing process is achieved through three principal factors: design for capillary action, choice of brazing alloy and temperature control. In regards to design, the partial cuts in the flat stock that fold to build three-dimensional geometries leave channels that the braze- clip filler alloy should flow through [81] during the heating process to become the corner brazes. 36 At the transition-fit interface for the tee-braze, specified clearance between the plates promotes penetration of the filler and spreading within. Brazing fillers with a high silver content, as used in this work, become highly fluid at molten temperatures and exhibit extensive spreading on surfaces with good wetting [69, 79, 82]. While this spreading is a desirable attribute for the tee-braze with a 36-mm 2 area to fill between the doubled base plates (Figure 2.2(b)), the corner braze would benefit from a minimally-spreading filler. The positioning and volume control of the brazing alloy was achieved through the creation of braze-clips made of folded 30-AWG (AmericanWireGage) (0.254-mm diameter) silver alloy wire and the use of an appropriately-chosen paste flux as defined in Section 2.1.2.2. This tensioned clip slides into the channel created by the folding lines and holds itself in place via friction at the connection point, regardless of how the structure is oriented during the heating and melting process. This allows pre-placing and securing the filler for all connections simultaneously as well as preventing migration when the cleansing flux paste froths at high temperatures, evenly distributing a controlled volume. However, the thin wire gage requirement of the clip-brazing method restricts the choice of available product so we used the same diameter wire for both brazes, the results of which are discussed in Section 2.1.2.3. This brazing process concurrently requires the precise and consistent control of the temper- ature and temperature gradient throughout the structure, to promote the capillary action of the molten filler into the interface [83]. The mechanical properties of the brazed joint such as gap filling ability and process time depend on the geometry of the connection point (area, gap size, et cetera), surface energy interaction between filler-base and type of heat applied to the system[70]. Most brazing research is conducted under a controlled atmosphere and publications such as [75] recognize the lack of research on brazing in open conditions; our work is an example of open-environment brazing with quality joints and reversible oxidation damage. As it was empirically shown to be effective, we selected torch brazing as a method to quickly heat and manually control the local temperatures throughout the structure. While oxy-acetylene is the preferred fuel mixture to produce a feather-like reducing flame for a high-quality braze, butane 37 was chosen for ease of handling in the form of a compact micro-torch. Commercial butane torches are specified to have a flame temperature in air of approximately 1430°C, but recent research on axisymmetric butane flames shows that temperatures vary between [650:1500] ◦ C depending on relative distance to the fuel jet origin [84]. This range serves as an approximation to assist in the placement of the flame during brazing. In addition, small drops of flux are placed at various locations on the body and serve as a visual indicator of the local temperature, as the paste turns glasslike and translucent approximately 200°C below the solidus of the filler alloy. As the molten filler flows toward the heat, once the structure is maintained at the solidus temperature of the filler ( 743°C), the torch is focused on the areas behind the filler until the softer metal reaches its liquidus temperature at 810°C and spreads throughout the joint. 2.1.1.2 PotentialApplications Brazing as a metal-joining method has existed for thousands of years, but the proposed new scaled-down technique can be applied to make a variety of novel three-dimensional rigid or flex- ible millimeter-scale structures, utilizing limited infrastructure and equipment at low cost. Braz- ing is suitable for most metals and some ceramics[80] which can be freely selected based on a given mechanical design constraint. In this case, Invar-36 (FeNi 36 ) was selected as the base mate- rial due to its extremely low CTE,α, of approximately 1.2× 10 − 6 · K − 1 . This manufacturing and design flexibility of material choice could be useful in new research on micro-manipulators and end-effectors, thermally-invariant micro-supports, passive dampers and mechanical sensors, to name a few examples. We implement the proposed fabrication process to make a micro-force sensor, which is a brazed flexible cantilever structure that translates applied force on the distal end to a parallel displacement which is externally measured. The parameter selection of the structure prescribes the sensitivity, resolution and frequency range of the effective sensing region [25, 65] that allows us to validate the proposed approach through a comparison with [85] and measureµN-resolution lift forces generated by a microrobotic insect. 38 Figure 2.3: Force sensor model and mechanism of action. (a) CAD render of the structure. (b) Photograph of the sensor used in this work. (c) Starting at a nominal offset δ nom , forces from the affixed body on the distal face translate into pure vertical motion of the target plate, measured with the capacitive displacement sensor. 2.1.2 TheuseofClip-Brazingtofabricateamicro-forcesensor 2.1.2.1 Micro-ForceSensorDesign Criteria: To demonstrate the fabrication method, we wish to design a structure that transforms periodic forces in the expected range[− 2: 2] mN into a pure translational distance that can be measured by an external sensor. This structure must have a high off-axis rejection ratio, and a loaded input-output lowpass bandwidth of approximately one order of magnitude larger than the maximum periodic frequency (100 Hz) while maximizing the sensitivity inside this range. This specification guarantees that the sensed forces do not have fundamental frequencies that will be near or at structural resonance. Further analysis assesses the allowable frequency measurement range based on acceptable error, as is presented in a later section. The resolution of this sensing system should be smaller than 1% of the gravitational force from the body mass of the affixed microrobot, which for an 80-mg body results in a target resolution of 7.85µN. SensingMethod: Commercially available capacitive displacement sensors have a satisfactory dynamic response to measure time-varying inputs in the desired operating frequency range, and as analog devices they are well-suited for real-time data filtering while retaining high fidelity [86, 87]. These sensors are one of the most widely-used small distance measurement devices due to their contactless nature, sub-nanometer resolution, high speed and low nonlinearity error [88]. 39 Changing the distance between the sensor’s metallic probe head and a flat conductive surface generates a change in capacitance, which is directly translatable into distance from the manu- facturers specifications and readily transformed into force employing a lookup table obtained through a weight calibration procedure. However, for highly sensitive probes, the measurement range is on the order of tens of microns and are affected by any off-parallel position error of the measurement surface. The device used here is the 8-mm-diameter cylindrical Physik Instrumente PISeca D510.021 capacitive probe with a 20-µm range and 0.2-nm resolution, which requires a measurement surface 50% larger than its surface area. A metallic target disc is incorporated into a structure (as identified on the model in Figure 2.3(a) and shown on the finished product in (b)) that is designed to take force input and transform it into a quasi-parallel deflection with respect to and within the measurement range of the sensor, shown in Figure 2.3(c). The capacitive probe is aligned to the target with a 3-DOF bracket assembly. The smaller the probe head, the smaller that additional mass contribution is, but this dimension is restricted by commercial availability. PhysicalDesign The sensing structure must deflect when forces are applied along the desired measurement axis and reject forces along other directions. A natural candidate for this structure is a rigidly fixed cantilever beam, which is significantly more sensitive to applied loads perpendic- ular to its longest longitudinal dimension due to the smallest moment of inertia in that direction. This concept is widely used in designs described in the literature for force sensing[54, 89, 90]. The entire design comprises a double-cantilever structure and a target plate as in [65]. The con- nected pair of beams constrain the natural cantilever rotational deflection at the distal end to restrict the front face of the structure to move approximately parallel to the base for small angles and deflections, to be compatible with the sensing system described in section 2.1.2.1 and Figure 2.3(c). As such, the cantilever pair can be modeled as a single beam with half the load on a roller, see Figure 2.4(a), and the microrobot can be simply physically modeled as in Figure 2.4(b). In order for the structure to function as a force sensor, its natural frequency must be many times higher than the fundamental frequencies of the measured time-varying forces. Otherwise, the measurements could be distorted by the natural response of the cantilever configuration. 40 Table 2.1: Material Properties of Invar-36 and Silver Property Invar-36 Silver Density (kg· m − 3 ) ρ 8100 10500 Modulus (GPa) E 141 76 Poisson’s Ratio ν 0.29 0.27 Loss factor η 0.0007 — CTE (× 10 − 6 · K − 1 ) α 1.2 9.2 Modeled as a mass-spring-damper system with a constant mass and natural frequency in Hz of f n = 1 /2π( k /m) 1 /2 , a stiff beam has a high natural frequency but low sensitivity due to smaller tip deflections as compared to a less stiff beam that deforms more readily but has a lower f n . This phenomenon creates a design trade-off because the measurement range depends on the mass of the robot that generates the forces to be measured, as this constant payload changes the dynam- ics of the dual cantilever-robot system. To select the dimensions of the structure, the fabrication, materials, and static and dynamic factors must be considered. In [65], this problem was investi- gated through an optimization approach; here, we employ simulations and a simplified numeri- cal analysis based on first principles. We chose Invar-36, whose relevant properties are listed in Table 2.1, as the structural material in order to minimize thermal instability during sensing ap- plications. The nominal static payload is 83 mg, and generates a maximum instantaneous force of 2 mN. Mechanics equations provide deflection models for rigidly-clamped cantilevers that de- scribe a structure that deflects linearly within the 20µm range of the capacitive probe (see section 2.1.2.1). The effective total stiffness k tot and total mass m tot (comprising the distributed weight of the body of the sensor, the robot chassis and expected maximum lift) must be determined to obtain a simple analytic model to estimate the peak response and the bandwidth of the sensing system. This is compared against numerical data from finite element analysis (FEA) simulations in the next section. Models: To obtain an analytic model of the sensing structure, we assume the double-cantilever to have a fixed base and a load F 0 acting at its distal end, as depicted in Figure 2.4(a). Because the ends of both beams are rigidly coupled together, they are constrained to move in a purely 41 vertical manner; this boundary condition is modeled as a roller associated with a small beam deflection angle α. From symmetry, this structure can be thought of as two independent beams in parallel with their bases clamped and moving extremes guided; each loaded with F 1 =F 0 /2 as in Figure 2.4(b). Then, from basic solid mechanics, it follows that the vertical deflection of the single beam in Figure 2.4(b) is given by ∂ 2 y 1 ∂x 2 1 = M 1 (x 1 ) EI , (2.1) where M 1 (x 1 )= F 1 (L/2− x 1 ) is the moment distribution along x 1 , and E and I are the elastic modulus and moment of inertia of the beam, respectively. Thus, we have that y 1 (x 1 )= F 1 EI L 4 x 2 1 − 1 6 x 3 1 , (2.2) which allows us to find the boundary condition δ 1 = y 1 (L) = F 1 L 3 /12EI and a simple static relationship betweenF 1 andδ 1 (stiffness constant) with the form k 1 = 12EI L 3 . (2.3) Therefore, we can estimate the total equivalent stiffness of the double-cantilever structure as k tot = 2k 1 = 24EI L 3 . (2.4) Similarly, we can find an equivalent mass of the single beam in Figure 2.4(b), m eq , assuming that it is concentrated atx 1 =L. As done in [65], we match both the boundary conditions and kinetic energies of the cantilever and the equivalent system, i.e., 1 2 m eq ˙ δ 2 1 = ρ 2 Z L 0 ˙ y 2 1 (x 1 )dV = ρA 2 Z L 0 ˙ y 2 1 (x 1 )dx 1 , y 1 (x 1 )=δ 1 3x 2 1 L 2 − 2x 3 1 L 3 , (2.5) 42 Figure 2.4: (a) Diagram for the analytic model showing how the vertical boundary condition is imposed with a roller. (b) Output of the numeric SolidWorks ® simulation for static deflection generated by the point force F a /2= 0.001 N (c) Frequency response: normalized amplitude of deflection (dB) vs sinusoidal input frequency, or a Bode Magnitude plot. The quasi-static region falls below approximately 100 Hz. where ρdV =ρAdx 1 is an element of mass, ρ is the density of the beam’s material, dV is an element of volume andA is the area of the beam’s cross section. Thus, noticing that (2.2) defines a static mapping, regardless of the parameters of the single cantilever (Figure 2.4(b)), we can estimate the associated equivalent mass as m eq = 13 35 ρLA= 13 35 m b , (2.6) wherem b is the total mass of the unloaded single cantilever in Figure 2.4(b). Once again consid- ering the symmetry of the double-cantilever, it follows that the undamped frequency response of the unloaded sensing structure is completely determined by the pair k tot ,2m eq +m front with resonance at q k tot / 2m eq +m front , where m front is the mass of the front plate that joins both cantilever beams. Furthermore, when the aggregated mass of all the constant loads acting on the 43 double-cantilever, m tot , is considered, the undamped frequency response of the system is com- pletely determined by the pair{k tot ,m tot } with natural frequencyω n = p k tot /m tot . In this case, m tot = 26 35 m b +m front +m target +m robot , (2.7) wherem target is the mass of the target plate andm robot is the mass of the robot statically measured with a balance. Note that according to this modeling method, bothk tot andm tot remain constant during the performance of experiments; therefore, the dynamics and properties of the sensing structure, such as bandwidth, also remain constant and they do not depend on the input force (time-varying load) generated by the robot while flapping. To find an input–output mapping from the instantaneous force acting on the double-cantilever, u(t), to the instantaneous deflection at the distal end of the structure, δ(t), one option is to in- troduce a dissipative force term proportional to ˙ δ(t) to obtain m tot ¨ δ(t)+c tot ˙ δ(t)+k tot δ(t)=u(t), (2.8) in which c tot is a parameter that can be tuned experimentally or via simulation, which can be an experimentally complex procedure. Here, for simplicity, we use the hysteretic loss factor η in Table 2.1 to represent the energy dissipation in the structure when experiencing a harmonic excitationu(t)=U 0 e iω 0 t ,i= √ − 1. Namely, U 0 e iω 0 t =m tot ¨ δ(t)+k tot (1+iη)δ(t), (2.9) whereU 0 is the amplitude associated with the fundamental frequency of excitationω 0 . By solving (2.9) for a range of frequencies, the frequency response can be estimated. However, the main 44 Table 2.2: Geometric Parameters of the Structure. Parameter Symbol Value Sheet Thickness (mm) τ 0.153 Beam Length (mm) l 7.0 Beam Width (mm) w 3.0 Inter-Beam Spacing (mm) b 3.7 Diameter of Target Plate (mm) d 10 Table 2.3: Inertial and Dynamic Parameters of the Structure. Parameter Symbol Analytic Numeric Total Mass (mg) m tot 240.0 238.6 Total Stiffness (kN · m − 1 ) k tot 10.19 9.52 Natural Frequency (Hz) f n 1040 1112 Cutoff Frequency (Hz) f c 1615 1724 Sensor Error Limit (Hz) E 1% < 100.9 < 104.7 Sensor Error Limit (Hz) E 2.5% < 167.5 < 171.2 advantage of using (2.9) is that the magnitude of the frequency response can be directly estimated according to |δ(ω 0 )|= U 0 k tot q 1− ω 2 0 /ω 2 n 2 +η 2 . (2.10) Note that sinceη is a property of the material (Invar 36), the estimation method given by (2.10) is entirely analytical. For design and analysis purposes, we employ the set{S ω 0 |f 0 =ω 0 /2π∈[0 : 2000] Hz}. The final dimensions of the sensing structure were selected through an iterative design pro- cess similar to those in [65] and [25]; the corresponding geometric parameters are shown in Table 2.2. These values are also used in the implementation of SolidWorks ® numeric simulations such as that depicted in Figure 2.4(c). The analytically-found inertial and dynamic parameters of the structure corresponding to the geometrical quantities in Table 2.2 are shown in the third column of Table 2.3; the associated magnitude of the frequency response is shown in blue in Figure 2.4(d). In this case, the input–output relationship betweenU 0 and|δ(ω 0 )| is normalized such that the DC-gain is 0 dB. Note that the normalizing gain is simply k tot / p 1+η 2 . In the implementation of the numerical simulations and analyses, the mass distribution of the complete 45 double-cantilever–robot system is modeled with the simplified structure shown in Figure 2.4(c), in which a parallelepipedal solid protrusion with the same weight as that of the studied force- generating microrobot is fixed at the center of the front face of the cantilever. Employing this configuration, in a series of two numerical-based studies, we first perform a simulation of the static structure in the absence of time-varying loads in order to determine the expected nominal static deflection δ nom . Thisδ nom is then employed to unbias the sensor’s output and normalize the frequency response magnitude (the DC-gain is set to 0 dB). Then, we perform a set of dynamical simulations in which the system is excited with sinusoidal loads with frequencies in the range [0 : 2000] Hz. The inputs and outputs of these dynamic simulations are employed to estimate the magnitude frequency response shown in red in Figure 2.4(d). Additionally, the results from both the static and dynamic numerical simulations are employed to calculate the inertial and dynamic parameters of the structure shown in the fourth column of Table 2.3. Note that exactly as in the analytical case, the numerically-estimated magnitude frequency response in Figure 2.4(d) was normalized to set the DC-gain to zero for ease of comparison and analysis. As can be clearly seen in Table 2.3 and Figure 2.4(d), both the analytical and numerically- estimated models predict similar natural frequencies f n , at 1040 Hz and 1112 Hz, respectively. Consistently, both models exhibit similar low-pass bandwidths with cutoff frequencies f c at ap- proximately 1650 Hz (i.e., the normalized amplitude of deflection is − 3 dB at approximately f c = 1650 Hz). For measuring purposes, however, the system is required to operate in the fast quasi- static range, which in Figure 2.4(d) is marked by the shadedred region, in order to avoid contam- ination of the measurements by the resonance responses of the structure. For input frequencies under 105 Hz, the models predict less than 1% error between the true signal and the measurement (E 1% in Table 2.3, row 6) and less than 2.5% below 167.5Hz (E 2.5% in Table 2.3, row 6). In addition, the simulated static deflection δ max for the maximum expected load per beam of F 1 = 0.001N is 0.105µm while the analytic model predicts δ = 0.098µm. Thus, these analyses predict that when the sensing structure is excited by forces with frequencies around 100 Hz (i.e., significantly below the resonance of the system), the factor of safety (FOS) 46 is almost 100 for bi-directional motion and the measurement range is more than 2700 times the nominal resolution of the sensing system. Note that the sensitivity of the sensing system can be easily increased by simply lengthening the cantilever beams to make them more compliant to force; this, however, would also reduce both its resolution and bandwidth. Since the fundamental frequencies of the time-varying forces of interest are less than 100Hz, the sensing mechanism meets the performance criteria set in section 2.1.2.1, provided that the brazed structure behaves as if it were solid Invar-36. The next section explains the procedure for constructing such a body. Figure 2.5: Fabricationprocessofthesensorusingfoldingandclip-brazing: Step 1: Com- ponents are cut from Invar shim stock, showing a detail view of the folding lines. Step2: Folding the structure and inserting the double-cantilever into the doubled base. Step 3: Adding brazing alloy wire clips. Flux paste is used to coat the surface. Step 4: Using the butane torch to melt the filler and fuse the connection points. Step 5: Adhering the target plate with conductive epoxy, with photographs of the final sensor. The plots of the magnitude frequency responses obtained through both the simplified first- principles model and simulations, shown in Figure 2.4(c), and are used to assess the suitability and correctness of the design. In this figure, the ordinate axis indicates the normalized amplitude of deflection in decibels (dB) and the abcissa axis indicates the frequency of the sinusoidal dynamic load. The two models predict similar peak responses for the first mode at a natural frequency f n of approximately 1100 Hz and a low-pass bandwidth of approximately 1650 Hz (i.e. the normalized amplitude of deflection is − 3 dB at f B = 1650 Hz) as shown in Table 2.3. For measuring purposes, however, the system is required to operate in the fast quasi-static range, which in Figure 2.4(c) is marked by the shadedred region in order to avoid contamination of the measurements by inertial resonance forces of the structure. Consistently, for input frequencies under 105 Hz, the models 47 predict less than 1% error between the true signal and the measurement, (E 1% in Table 2.3, row 9), and less than 2.5% below 167.5 Hz. The fundamental frequencies of the time varying forces to be measured are less than 100 Hz, which indicates that the sensor bandwidth meets the design criteria, provided that the brazed structure behaves as if it were solid Invar-36. The next section explains the procedure for constructing such a body. 2.1.2.2 Fabrication Components are cut with a diode-pumped solid state (DPSS) laser beam, (Photonics Industries DCH-355-3) steered by a galvanometer (ScanLab HurryScanII-14). For the correct fabrication of the double-cantilever structure, it is critical that the beams are precisely positioned with interior angles as close to 90 ◦ as possible (angleθ in Figure 2.4(a) near-zero), while remaining fixed in a rigid base. Therefore, the design is divided into two parts: the beams and front face that fold to have proper inter-beam spacing; and the slotted base that justifies the position and angle. This ensures that the internal angles of the assembled structure maintain the nominal design values. Prior research has used methods such as a combination of epoxy and soldering [65] and laser welding [25] to fix the relative positions of the structural elements into a heavy steel base. This work adds to that research by combining the fabrication ease and accessibility of the former and the precision of the latter, with a small footprint that can be easily situated in various locations with only cyanoacrylate (CA) glue (Henkel Loctite 416). The proposed brazing procedure can be simplified into the following steps, as shown in Figure 2.5 and detailed in the following subsec- tions: LinePathsandCutting: From a given set of dimensions of the double-cantilever structure, the design is transformed into a two-dimensional schematic composed of through-cuts and folding- line cuts (shown in the detail view ofStep1 in Figure 2.5); the perimeter is released with through- cuts and the folding lines are scored with only a few repetitions of the pulsed laser path. After cutting, the flat structure is cleaned in a 99% isopropyl alcohol ultrasonic bath. 48 FoldingandPositioning: By applying torque on either side of the folding lines, the material begins to enter a ductile failure mode and tears open cleanly to make connections as shown in Step 2 of Figure 2.5, and in detail in Step 3 of Figure 2.5. The two 90 ◦ folds at the cantilever end face force the shape into the three-sided box that forms the bulk of the structure. The fourth and final side is made by the double-thickness slotted base plate into which the beams are manually press-fit. The beam that supports the target plate is also manually folded into the proper relative angle, with the top face parallel to the main beams. Slight off-parallel fabrication errors can be compensated for with the 3-DOF structure that holds the capacitive probe. The structure is then soaked in a methanol bath to degrease the surface to prepare for brazing. Brazing Filler and Flux We selected 80/16/4 weight (wt.)% Ag-Cu-Zn as the brazing filler alloy, a common metalsmithing material with a melting point of 743°C. This metal falls in a gen- eral category of high-silver-content brazing fillers, known for its favorable wetting on iron-group metals without significantly penetrating the base structure via diffusion [70, 73, 82, 91, 92], to pre- serve the material characteristics of the Invar-36 as much as possible (CTE, strength, stiffness, et cetera). During fabrication, the amount and manual placement of the flux paste (Harris Prod- ucts Stay-Silv Black) must be carefully controlled because when the ratio of the volume of flux to weight of filler is too high, the unpredictable kinetic activity of the boiling flux can transport the filler to an undesirable location on the structure. This can be particularly problematic here because the CTE of silver is nearly 10 times that of Invar-36, potentially making the structure more sensitive to environmental conditions, especially temperature variations. The braze-clip is made from a segment of 30 AWG brazer wire, cut to double the beam width and initially bent in half around a steel shim to make a tight, ‘hairpin’ bending radius. This tensioned clip is main- tained in position by friction, as depicted in Step 3 in Figure 2.5, which allows the flux to boil away without disturbing the filler material and evenly distributes a controlled volume. Note that this symmetrical clip creates an internal fillet which additionally helps to promote the desired boundary condition of pure vertical motion of the front plate. 49 Application of Heat: Butane has a nominal flame temperature in air that is in excess of the melting point of the Invar-36, 1427°C, so the heat must be applied evenly and sparingly to prevent unintentional melting or buckling of the parent material. An adjustable-flame micro-torch was used to precisely control the braze locations (Step 4 of Figure 2.5), performing the operations in seconds to avoid the progressed grain growth that happens in the Invar-36 when heated above 1035°C. When the flux becomes glasslike at approximately 565°C, the torch is focused on the area into which the filler is desired to flow. Once the silver melts, the structure is immediately removed from the heat, and the part is subsequently air-cooled to room temperature. Due to irregular heat transfer and the relatively significant difference in thickness between the base plate assembly and the thin corner brazes, the filler flows at different times when the body is heated uniformly. To account for this phenomenon, the heat is initially focused on the thickest part of the assembly; then, it is steadily refocused on the distal points of the structure. This sequential procedure can result in slight wicking of the filler down the length of the dual- cantilever structure if its temperature exceeds a critical point; however this excess material can be removed once the brazing process is completed. For more precise control, this step can be done in two stages using distinct alloys with different melting temperatures, thus avoiding heat damage. Because the brazing process is not performed under inert gas or in a vacuum due to cost and infrastructure constraints, the surface of the Invar-36 material becomes oxidized and nonconductive. To restore the surface, the structure is soaked in 60°C muriatic acid (37% by mass hydrochloricacid (HCl) solution) for 30 minutes, mechanically polished and oiled to prevent rust. AttachmentofTargetPlate: The Invar-36 target plate is mirror-polished and adhered to the brazed structure with conductive epoxy (MG Chemicals 8100) as shown in Step 5 of Figure 2.5, reinforced with CA adhesive. This ensures that the target plate is electrically grounded as well as rigidly fixed to the dual-cantilever. Although it is possible to join the target plate to the struc- ture during the initial brazing step, it is difficult to maintain in the correct position and angular orientation while heat is applied with the butane torch. In addition, by keeping the target plate as a separate piece, the need for further cleaning and polishing is avoided. 50 Figure 2.6: SEM images of braze cross-sections cold-mounted in epoxy resin; this weakly con- ductive base adds a distracting and irregular background, which has been cropped and replaced with the same shade as the least conductive neighboring regions. Zero image processing has been done besides vector annotation on the actual body of interest, and the cropped border can easily be seen to be distinct from the edges of the Invar-and-Silver cross-section which can be seen at high magnification. (a) Corner braze showing filler metal contact angle and spreading (Detail 1), as well as possible inclusions (Detail 2) and possible intermetallic phase formation (Detail 3). (b) Tee-braze with doubled baseplate, showing similar fillet sizes and large spreading throughout. Detail 4 shows that the spreading even continues beyond the interface and onto the bottom face. The braze on the left is the largest contact angle outlier of the set. (c) Another corner braze with a less-than-90 degree angle but similar fillet/corner braze profile. Spillover from excess filler in the corner is shown in Detail 5., and Detail 6 shows the plastic deformation region from folding. 2.1.2.3 EvaluationoftheFabricationMethod The best approach to evaluate the fabrication method is direct inspection of the completed struc- ture employingscanningelectronmicroscopy (SEM) on cross-sections of the brazes, cold-mounted in epoxy resin and mechanically polished with aluminum oxide and diamond lapping compound to 0.5µm grit, as shown in Figure 2.6. These images show that while the parallel beams appear to be perpendicular with respect to the front face from the far-field, the interior angles are less than 90 deg. This appears to be due to the slight increase in parallel inter-beam distance from the plastic deformation at the folding line location (detail 1 in Figure 2.6(a)), which forces the beams 51 Figure 2.7: (a) Photograph of the experimental setup with capacitive sensor, dual cantilever and flapping-wing microrobot. (b) CAD render of calibration loading setup with a micrometer- mounted support for the weight. (c) System thermal response to the introduction and removal of high-powered lights (d) Calibration loading signal examples for two weights, showing signal-to- noise ratio difference as well as the vibrational disturbances from the loading conditions. to bow inward. We believe that this misalignment can be rectified by introducing a fixed offset applied to the folding line location, employing measurements from these images. These SEM pictures also show the homogeneity, non-porosity and extent of the brazing filler spread, as well as consistent and small contact angles. Estimations of the Young contact angle of the final brazes (an example angle is shown in detail 1 in Figure 2.6(a)) from the SEM images show a range of 15 : 23 deg for six distinct specimens and a single outlier of 37 deg. These angles, all significantly less than 90 deg, show the high wettability [71] of the 80/16/4 wt.% Ag-Cu-Zn filler on Invar-36 [73]. Qualitatively, the alloy spreads deep into the interface in Figure 2.6(b), forming a continuous fillet, consistent with high wetting as defined in [72] and spreading [70, 71, 73]. In all the examined cases, no cracks can be observed in any fillet, and the only morphological irregularities appear in the diffusion zone such as the clumps in detail 3 of Figure 2.6(c). The characteristic asymptotic concavity of the fillets between the perpendicular faces of the tee-brazes are smooth and regular, as the abundant volume of brazer fills the seam between the double- thickness base-plate (detail 4 in Figure 2.6(b)), spreading across the full interfacial area of 36 mm 2 . 52 The corner-brazes (as shown in Figures 2.6(a,c)) are less uniform; small regions of poor adhe- sion or voids such as those at detail 5 in Figure 2.6(a) only appear in these outer brazes, possibly due to the slightly scalloped surface of the interior faces of the folding cuts of the Invar-36 due to the pulsed laser characteristics. Detail 6 in Figure 2.6(a) suggests that some precipitate/IMC formation might occur during the brazing process; however, this photographic evidence is nei- ther conclusive at this scale nor consistent with visual indicators of IMCs such as those in [73, 92] that can result in residual stresses at the joint. The images of the corner-brazes indicate that more control could be achieved by reducing the diameter of the filler wire, as the excess volume spills out of the channel and spreads onto the outer surfaces. This extra material, however, can be mechanically removed as part of the normal cleaning and polishing process after the HCl bath although we note that the abrasion may cause the unintentional loss of the concave surface of the outer corner-brazed connection (see difference between Figure 2.6(a) and (c), (a) has been polished on the outer surface). 2.1.3 CalibrationAndValidation The empirical relationship between output voltage from the capacitive sensor’s signal conditioner and force applied at the end of the cantilever is determined by employing a weight calibration on the experimental setup in Figure 2.7(a,b). The capacitive sensing probe is aligned to the resting position of the unloaded structure using a 3-DOF bracket assembly constructed from a 2-DOF optical flexure mount (Siskiyou Corp. IXF-50.ti) and linear micrometer stage. The robot is ad- hered to the front plate of the dual-cantilever to align the force input to the principle sensing axis of the structure. The calibration is performed with the 80-mg microrobot attached to the struc- ture via a 3-mg carbon fiber truss; the applied force is aligned with the nominal line of action of the lift force. We constructed micro-scale test payloads between 4 mg and 850 mg (measured with a Mettler Toledo XP6 microbalance with 1-µg accuracy) which apply a known, constant gravitational weight when suspended from the “lower" extreme of the robot, as shown in Figure 53 Figure 2.8: (a) Very close matching of linear calibration fits, C1-C3 were taken on three consec- utive days and C4 was taken three weeks later. The sensing system remains stable over several days, stiffening over weeks of heavy use, possibly due to work hardening of the silver or the inherent temporal instability of Invar-36. (b) The use of the micro-force sensor as a static scale: measurement and analysis of the sample voltages and estimation of the true weight, compared with the measured value from a microbalance. Likely due to human error in measurement for small weights, the relative error is less than 4% for masses 28-430 mg, but for largerσ SNR the error is less than 0.3% as the difference in the voltage signal magnitude becomes more significant. 2.7(b). During operation, the analog voltage ouput from the sensor signal conditioner is read us- ing a Mathworks xPC Target system (Version 5.5, MATLAB R2013b) and a National Instruments PCI-6229data-acquisition (DAQ) board employing a 10 kHz sampling rate. Although shown to be sensitive to temperature in Figure 2.7(c), the unfiltered output voltage signal is demonstrated to have a stable equilibrium position for time periods less than approximately 1 second (unloaded zero signal can be seen between t =[1.5 : 1.75] min in Figure 2.7(c)). While this temporal sta- bility in the zero position is sufficient for measurements of multiple flapping periods, it is also susceptible to disturbances which are discussed in more detail in the following sections. To minimize the probability of measurement errors from the temporal instability of capacitive sensors, each calibration load must be applied instantaneously. To achieve this difficult task, each test weight is composed of a small steel washer and a hanger made of thin copper wire, as seen in the close-up of Figure 2.7(b). In this configuration, during suspension, the hanger’s hook is located exactly above the washer’s center of gravity, which also intersects the robot’s yaw axis ˆ z. To suspend each test-weight, a vertical micrometer stage stably holds the washer’s flat bottom exactly under the robot with the hook poised above but not contacting the loading point. When the stage is rapidly lowered, the bottom support for the test-weight is removed, which almost 54 instantly loads the sensing structure; two examples of the resulting Heaviside-like voltage signals are shown in Figure 2.7(d). As desired (and expected), when the sensing structure is unloaded, the voltage measurement returns back to its original pre-loaded value, thus producing a square pulse signal. For each test-weight, the pre-load–post-load voltage difference is measured five times in or- der to verify the repeatability (small standard deviation) of the measurements. The collected data are employed to find a least-squares (LS) fit for the linear mapping from applied force to probe voltage, as shown in Figure 2.8(a) for four different calibrations. Here, it can be seen that the procedure produces consistent results as only small changes are observed in the slopes of cali- brations performed at different times. For example, calibrations C1 to C3, taken over the course of three days, vary less than 3% in system sensitivity between each other; and after three weeks of heavy use and inherent temporal instability, the sensitivity varies less than 15% (C4). The lin- earity of the experimentally-obtained calibration mappings in Figure 2.8(a) is consistent with the fact that both the capacitive sensor and dual-cantilever structure are expected to exhibit linear quasi-static behaviors for small deflections. The inverse of the calibration function is the map- ping that translates the analog time-varying voltage signal v(t), a function of the displacement of the dual-cantilever, back into the force domain. The regression R 2 -value is above 0.997 for all calibrations, andstandarderrors of the mapping are less than 11µN. These statistical parameters provide compelling evidence that the identified first-order mapping adequately describes the re- lationship between the probe voltage and the instantaneous force within the allowable frequency range. As a simple static test, we employ the system shown in Figures 2.7(a,b) to measure unknown payloads which are compared to the weight as measured using the microbalance (Figure 2.8(b)). Ten voltage readings per unknown weight were taken and mapped to milligrams to validate the precision of the sensing system. At the lower limit of the sensing capabilities, the calibra- tion procedure has a small signal-to-noise ratio (SNR dB =( V load,RMS/V noise,RMS ) 2 ) and it is difficult to distinguish the voltage that corresponds to the loaded state from the noise (see Figure 2.7(d)). 55 However, assuming we have sampled a statistically normal population, all mapped weights fell within the95%confidenceinterval (CI) as well as thestandarderrorofthemean (SE) for the regres- sion mapping. As the unknown payloads increase in weight and SNR dB gets larger, the relative error (RE) between true and mapped weight decreases from 4 to 0.3%. This does not imply that the sensor is not accurate for low weights, only that the system is difficult to manually load and unload which is most noticeable at this scale in the form of vibrational disturbances that distort the signal. These issues are assumed to be negligible during the normal working conditions; the robot that generates the forces to be measured is pre-fixed to the sensing structure and these forces are in the design input range. During the performance of force measurement experiments, the analog output signal from the capacitive sensor is susceptible to two types of signal perturbations, drift and noise. These two types of disturbances might significantly impact the determination of the calibration mapping (and hence its sensitivity) as well as its resolution which is proportional to the root-mean-square (RMS) value of the noise signal. Mechanical disturbances and 60 Hz AC line noise can be damped but there is noise from the DAQ process as well. In addition, capacitive drift, a change in the measured capacitance over time between two nominally fixed components, highly dependends on changes in the local environmental conditions such as temperature and humidity [88] [87]. Even slight air currents might cause detectable changes in the sensor voltage readings, an effect that is greatly mitigated by isolating the entire experiment inside an acrylic enclosure. Using the experimental methods discussed here, the drift could be reduced to less than 0.0001 sec (cor- responding to an effective drift in measured force of 0.01 mN· s − 1 ) which results in near-zero distortion of the signals across a small number of flapping periods. Also, the drift is approxi- mately linear in short time scales; therefore it can be accounted for and removed from longer data sample sequences. Although remaining linear, this drift is amplified with the use of the high-powered lights required for high frame rate videography, but remains small and removable while at thermal equilibrium. However, as all materials have a non-zero CTE, the greenhouse effect of the lights on the acrylic enclosure drive the system to a new equilibrium as shown in 56 Figure 2.9: (a) Dimensions for the dogbone specimens for the tensile test and the reinforcement of the folding-line cut. (b) Five brazed specimens, S1-S5, were tested against untreated Invar-36. (c) Tensile test data for the baseline and five brazed samples. Except for the outlier S5, the brazed specimens all exhibited yield strengths and ultimate strengths higher than the normal range for plain silver (shaded red region), but well below the range for Invar-36 (shaded blue region). The specimen with the most silver, S3, had much higher ductility and more than twice the percent elongation at break. Figure 2.7(c). Here, it can be observed that the sensing system responds to intense illumination by decreasing its output from 7 V to 1 V as time approachest=t exp = 55 minutes, at which point the lights are switched off and the system output rises to 5 V by approximately timet=t exp = 85 minutes. This behavior represents a significant combined thermal expansion of the probe and sensing structure. Note that the shape of the signal in Figure 2.7(c) indicates that this thermal response can be modeled by a first order linear time invariant (LTI) model with an approximate settling time of 55 minutes, depending on the temperature of the room. In assessing the quality of fabrication methods, nondestructive inspection such as microscopy and ultrasound can only provide information on the homogeneity and crystal structure of the braze. However, only destructive methods such as tensile or impact testing can be used to directly 57 Table 2.4: Mechanical Failure Properties - Tensile Test Invar S1 S2 S3 S4 S5 0.1%σ y (MPa) 442.2 362.7 360.4 386.8 — — σ TS (MPa) 648.9 415.2 388.7 455.1 238.5 383.9 measure the strength of a brazed joint [70]. Nominal values for material characteristics such as tensile strength come with a strict definition of test specimen dimensions and aspect ratio [69, 93]; however, the workspace size limitations of our laser micro-machining system (as described in section 2.1.2.2) combined with the thickness of the Invar-36 shim stock (0.1524 mm) prevented the use of standard grip dimensions. We prepared typical dogbone specimens with a parallel length of 12.5 mm and width of 3 mm, shown in Figures 2.9(a,b) which were compared to an identically-sized pure Invar-36 specimen. While these tests can give us an order of magnitude of the tensile strength of the connections, they do not accurately replicate the in-situ loading conditions which include both shear and normal components. However, there is no easy method to test the structure as assembled due to its complex geometry and combination of joining types, as well as its designed deformation mode. As the clip brazing method was devised to make the brazer placement straightforward, simply creating a butt- or lap-joint from two separate pieces to assess the strength of the connection is not only difficult but also not necessarily representative of the true, physical case. In fact, it points to why the method is so useful. Without the tensioned clip and self-jigged structure, the filler flows unconstrained and the two halves become misaligned. Therefore, these samples were made in the same way that the clip-brazing folding components are made, with a folding line cut through nearly the entire thickness of the shim stock that is reinforced or backfilled using the silver filler (see Figure 2.9(a)). These tensile tests were performed with a benchtop Instron 5567 at a constant crosshead speed of 2.25 mm/min as per the ISO 6892-1:2016 standard for ambient tensile testing of metal samples. The load and displacement data from the experiment are mapped to engineering stressσ = F test/A o 58 Figure 2.10: Spread of the data for two calibration procedures. (a) Box plot of calibration C1, showing no outliers and decreasing spread of the measurements as the voltage difference in- creases in magnitude, with the exception of the 629-mg payload. (b) More spread in the data and two weights with obvious outlier measurements for the 58- and 343-mg readings. and engineering strain ε = ∆L test/L o as shown in Figure 2.9(c). The reference of pure Invar-36 is marked with the dashed line and has a 0.1% offset yield strength σ y of 422 MPa and tensile strengthσ TS of 649 MPa, which falls within the expected range for nominal Invar-36 properties [93] (region shaded in blue on the stress-strain diagram in Figure 2.9(c)). Three of the five samples exhibited moderate to large ductility and had a σ y between 363 MPa and 387 MPa. Four of the five brazed samples had similar σ TS , between 384 MPa and 455 MPa, falling within the expected failure stress range for pure silver, marked by the red shaded region on the plot. Sample 4 failed at a much lower stress, which might indicate that the braze filler did not adhere well to the parent Invar-36. Sample 3 had an atypically wide ductile region, which could be explained by the slight excess in filler material on that particular sample, sample S3 in Figure 2.9(b). The strength test results are summarized in Table 2.4. 2.1.3.1 ExperimentalForceValidation To demonstrate the suitability of the sensing system for measuring time-varying forces in the micronewton range, we compare the experimental average lift governed by the flapping-wing 59 Figure 2.11: Regression Analysis. Detailed views of the regressions on the calibration datasets and confidence intervals on the mapping. Regression outputs reported with representative statistical parameters Table 2.5: Measured Average Periodic Lift and Absolute Error by Frequency F f 55 F f 60 F f 65 F f 70 F f 75 F f 80 F f 85 F f 90 F f 95 F f 100 Avg F θ 75:85 [mN] 0.14924 0.18758 0.22216 0.25612 0.28735 0.30877 0.34690 0.36492 0.37628 0.39956 Avg F θ 80:98 [mN] 0.14148 0.16887 0.18624 0.21083 0.24165 0.26787 0.29378 0.30893 0.31650 0.33412 θ w [ ◦ ] k h [µN· m] ∆F f 55 ∆F f 60 ∆F f 65 ∆F f 70 ∆F f 75 ∆F f 80 ∆F f 85 ∆F f 90 ∆F f 95 ∆F f 100 θ 75:85 1.75 0.007438 0.001596 0.006312 0.004879 -0.001541 -0.027686 -0.028098 -0.047547 -0.052447 -0.025628 θ 80:98 1.84 0.016594 0.010632 -0.003627 -0.014020 -0.014427 -0.023851 -0.031586 -0.048183 -0.058039 -0.028555 θ 75:85 1.75 0.005538 0.002597 0.006818 0.007722 0.003619 -0.024103 -0.013203 -0.025338 -0.024563 0.018053 θ 80:98 1.84 0.015630 0.007888 -0.007363 -0.012531 -0.011621 -0.021786 -0.021368 -0.027389 -0.037998 0.000926 microrobot shown in Figure 2.12(a) with CFD simulations for flapping frequencies in the range [60:100] Hz. Research in [94] used a quasi-steady model to compare these types of forces and [23, 85, 95] have initial findings on the use of CFD for force simulations. The experimental setup can be seen in Figure 2.7(a-b) and during measurement in Figure 2.12(b), with the experiment viewed from below. In this figure, the experimental forces in millinewtons are computed from the output voltages of the capacitive probe, using the calibration mapping according to the method described in section 2.1.3. The robot’s wings are designed to be nominally rigid and are monolithically fabricated from carbon fiber and 2.5µm Mylar (PET) film, bonded together with curable thin-film adhesive. The microroboticflexure transmission, developed upon the methods in [1] and [4], translates approx- imately linear deflection from a piezoelectric cantilever actuator into an approximately linear wing angle output. The wing is attached to the point of rotation by a rectangular sheet of 7.5-µm 60 Figure 2.12: (a) The 80-mg flapping-wing microrobot with a US penny coin for scale. The si- nusoidal deflection input signals from two piezoelectric cantilever actuators drive the flexure transmission to flap with periodic motion. The passive pitching hinge on the wing generates a net positive lift force per period which is measured by the force sensor. (b) One period of flap- ping as captured from below by a 45 ◦ mirror, note the passive pitching. The ghost image near the leading edge is from the non-first surface reflection. (c) Instantaneous force profile for 75 Hz flapping with both wings engaged (halved to represent a single wing contribution), with raw and filtered data compared against the CFD simulation. (d) Instantaneous force profile for 75 Hz flap- ping with a single wing engaged, compared against the same CFD simulation. (e) Comparison of 40 simulated average lift forces for two overlapping frequency ranges against the empirical measurements. 61 Kapton polyimide film, which acts as a low friction hinge and and enables the passive rotation of the wing about an approximately fixed axis. During periodic flapping, this phenomenon gener- ates lift forces dependent on its pitching profile. Here we applied the method in [85] to numerically solve a simplified model based on the Navier–Stokes equations that describe the fluid-structure interaction problem for this case at each frequency of interest. This numerical solution is plotted against empirically-measured in- stantaneous force profiles in Figure 2.12(c,d), filtered with a zero-phase digital filter employing a cutoff frequency of f c = 1100 Hz. The experiment was conducted twice, for measured empirical wing angle inputs θ 1 =[75: 85] ◦ and θ 2 =[80: 98] ◦ , using a single wing with these restricted angle ranges in order to study the relevant pitching effects in this region. The magnitude of the average aerodynamic lift is highly sensitive to the passive wing hinge stiffness parameter k h , a value nominally 1.84µNm evaluated from the dimensions and material properties. However, we know that this parameter varies (softens) in time so we allow slight tuning of this input to 1.75µNm and see the resulting change in the solution space for simulated lift, as seen in Figure 2.12(e). The model-empirical matching from the distinctk h CFD inputs and two flapping data sets ranged from 0.5% relative error (RE) to 18%. For the high-amplitude test withk h =1.75µNm, the difference was less than 7.8% for all frequencies in the set. 2.1.3.2 Discussion Besides the limitations of the nominal resolution of the capacitive sensor, there are multiple fac- tors that can introduce discrepancies between measurements and simulated output, from imper- fect fabrication of the sensing system and robot to calibration inaccuracies. For example, the inputs and parameters to the CFD simulations are not modeled to capture non-ideal time-varing properties of the materials and geometrical differences between the physical and modeled wings. Similarly, the structure of the dynamical models employed in the simulations does not account for the existence of phenomena such as the flexibility of the carbon fiber chassis structure and wings as well as the out-of-plane deviation of the wing’s leading edge. Figures 2.12(c) and 2.12(d) show 62 Table 2.6: Modeled vs Empirical Data Analytical Numerical Experimental k tot (kN m -1 ) 10.19 9.52 9.55 Sensitivity (V N − 1 ) 98.2 105.0 104.7 Resolution (µN) 4.07 1.91 7.48 Range (mN) 203.7 190.5 190.9 Natural Frequency (Hz) 1040 1112 — the instantaneous force profiles for two cases in which the flapping angle θ is a sinusoidal signal with a frequency of 75 Hz. In Figure 2.12(c), the magnitude of the measured force generated by the robot flapping its two wings is divided by two and compared with the CFD simulation force of a single wing (raw in dashed pink and low-pass filtered in solid mauve). In Figure 2.12(d), the magnitude of the measured force generated by the robot flapping a single wing (raw in dashed cobalt and low-pass filtered in sky blue) is compared with the CFD simulation aerodynamic force. From these plots it can be observed that the halved two-wing-force matches the instantaneous simulation result better than the single wing force. It is not clear why the single-wing force pro- file has a periodic spike that is not predicted by the simulation. We postulate that the loss of the symmetric wing configuration increases off-axis deviation of the transmission and leading edge of the wing [95] which can introduce inertial disturbances in the linkage, a phenomenon that was observed in the videos of the experiments. These forces approximately negate each other and the average over the period is negligibly affected, which justifies the study of a single wing to investigate lift production as it is difficult to perfectly coordinate the flapping of two wings and additionally we do not wish to introduce flow disturbances from a secondary wing. The flexible clip-brazed structure was demonstrated to perform consistently for ten distinct specimens and quantified by stiffness, sensitivity, and range, empirically measured to be within 0.3% and 6.3% of the predicted values. Table 2.6 compares these experimental results with the analytic and numeric simulations for the sensing capabilities of the brazed system, calculated from the data on the specific structure used for the measurements in this work. 63 The stiffness is determined by the amount of deflection of the cantilever divided by the ap- plied force; experimentally this is determined using the manufacturer-specified displacement sen- sitivity of the D510.021 sensor (1µm − 1 ) described in section 2.1.2.1 to transform the voltage to displacement. The sensitivity of the brazed sensing system relates the output voltage of the ca- pacitive sensor given an applied force on the structure in units of N − 1 , which we obtain during calibration with gravitational force. The range is limited by the voltage output of the D510.021 sensor of [0:10] V, calculated as the total voltage divided by the sensitivity. The resolution defines the minimum force that the system is capable of measuring and depends on the resolution of the capacitive sensor, which is specified as 0.2 nm. The minimum forces to create such a displace- ment can be found by employing equation 2.2 and outputs from the numerical CAD simulation. The experimental resolution, calculated by the RMS of the noise signal to be 7.48µN, is almost double the analytically derived value and four times the numeric value. It is large due to noise in the analog signal as previously discussed, yet still satisfies the 1% body-mass criterion of 7.85µN. Note that the experimental value for the natural frequency was not measured due to our inability to physically transmit a proper signal for system identification procedures. As the sensor was modeled as homogeneous Invar-36, the agreement between experiments and simulations indicates that the clip-brazing method creates a robust bond that behaves simi- larly to the parent material. The signal-to-noise ratio ranged from a low of 12.5 dB for the 55 Hz tests to a high of 26.8 dB at 100 Hz, indicating sufficient fidelity of the sensor. 2.1.4 Conclusions We presented a new method for the fabrication of high-precision, millimeter-scale micronewton- resolution force sensors which we have validated with SEM, mechanical testing and agreement with CFD simulations for a range of frequencies in the nominal operating range of the flapping- wing MAV. We laser micro-machined 0.1524-mm Invar-36 sheets into mating shapes that fold into three-dimensional structures that are permanently joined by the novelbraze-clip method in- troduced here. The use of braze-clips produces uniform connections at precise locations, creating 64 useful structures such as the dynamic micro-force sensor shown in this work. We have fabricated clip-brazed structures in a range of dimensions from 2× 1.5× 1.5 mm to 17× 6× 6 mm. These designs are modular and could be fused together in a secondary brazing process with a lower- melting-point braze filler to create new configurations with different capabilities. This process is fast, inexpensive, requires little infrastructure and is scalable in both its physical potential for modularity as well as the ability to be batch-produced. For the dual cantilever in this work, the entire structure can be fabricated in two hours (from cutting to attaching the target plate) with less than $6 USD of consumables including silver, flux, butane, muriatic acid and Invar-36. We provided compelling evidence of the method’s capability to make high-precision struc- tures, demonstrated by the agreement of modeled characteristics (e.g. empirical beam stiffness), destructive tensile testing and use as a static scale, with standard errors of the mapping less than 11µN. Empirical measurements of the average lift forces generated by a flapping-wing micro- robot are shown to have reasonable matching with CFD simulations for two sets of flapping angles employing two distinct wing passive hinge stiffness values for the simulation. For the large-angle set withk h = 1.75µNm, the maximum relative discrepancy was 7.8% with an upper bound of 58µN on the absolute error between simulation and experiment. Although not per- fectly matching, the simulations show that the model can be effectively employed to determine reasonable approximations for wing pitching and lift force generation which can serve to inform the design of control inputs for maneuvering flapping-wing microrobotic flight, which we leave as future work. 2.2 Measurements&Data The journal paper on clip-brazing [21] focused on the novel fabrication method and itsvalidation through comparisons of average and instantaneous forces with CFD simulations. However, it was the culmination of four years of research to arrive at the point at which the empirical force mea- surements resembled the CFD solutions in shape and magnitude. In terms of design, the initial 65 Figure 2.13: Early Data and Displacement Attenuation. Measured capacitive sensor voltage dur- ing flapping in Blue and measured actuator displacement data in Black. The first column of each four-pane image shows the full experiment length for scale, and the second column shows a few flapping periods (a) 75 Hz flapping (b) 85 Hz flapping (c) 100 Hz flapping focus was on the Single-Actuator Bee, the 1-DOF robot that employs a single piezoelectric can- tilever to drive a single transmission and wing pair. Designed to move on a rail system, this robot was created for constrained altitude control experiments. However, in theory a well-constructed specimen should be capable of simple jumping tests (showing takeoff stability and sufficient lift to combat gravitational forces) as well as generating representative aerodynamic lift data on a fixed structure such as the dual-cantilever sensor experimental setup described in Chapter 1, Figure 2.7. Early flapping-wing experiments showed a systematic problem with the actuators. As can be seen in the displacement information in Figure 2.13, the actuator deflection in the positive direction is smaller in magnitude than that in the negative direction. This asymmetry is even more pronounced at high frequencies, in addition to obvious attenuation of the peak-to-peak deflection. These two factors alone severely limit the force production capability, but in addition, 5000 frames-per-second (FPS) video data showed possible collisions in the passive hinges which could explain the noisy data. At first, we just thought we had fabrication issues with the transmission assembly and actu- ator mounting process. I made many bees in an attempt to try to make a ‘perfect’ one but still had pronounced asymmetry and attenuation. Extracting the cause of these problems was unclear 66 Figure 2.14: Instantaneous lift forces for the Single-Actuator Microrobotic Bee because we thought we had appropriate dimensions, and it was particularly strange that the at- tenuation was nearly always exacerbated in the positive direction. We began to realize at this point that the intermediate-modulus carbon fiber we were using was simply not stiff enough to achieve the full flapping angle at the operating frequency of 100 Hz. Even using secondary stiff- ening mechanisms like carbon fiber braces and epoxy coatings, we could not achieve the correct flapping or pitching angles. Presented in the plots in Figure 2.15(a,b), the average forces from the wings was not consistent between robots. In Figure 2.15(a) we see nearly double the forces to the sets in (b). The tests in Figure 2.15(b) were internally consistent while experiencing degradation over time. Figure 2.15(c) shows one of the best measurements of the Single-Actuator Bee, at a flapping frequency of 80 Hz. The average lift per period and the required lift for liftoff are shown by the dashed lines. Although this test does not prove that there is enough lift to take off, this bee actually had several successful jumping tests (conducted at the partially attenuated 100 Hz input f ) before I adhered it to the sensor, just to be sure. As explained in Section 1.2.1, we were able to obtain the correct modulus and FAW carbon fiber prepreg in January 2018. By this time, initial fabrication work on the double-actuator bee (which was started by myself and shared with my labmate Ariel, as was a large amount of this project) had borne fruit. The lack of consistency and visible problems with the wing pitching and off-axis deviation of the leading edge of the wings was the deciding factor to abandon the Single- Actuator Bee. No matter how perfect the robot looked, the measurements never looked remotely similar to baseline references in [94] and preliminary results from in-house CFD simulations. 67 Figure 2.15: Single-Actuator Microrobotic Bee: instantaneous and average forces. Figure 2.16: Double-Actuator Bee: Attenuation of output flapping angle at high frequencies. Left: low-voltage-amplitude stroke produces small wing flapping angles, with the attenuation increasing with frequency (see red marker). Right: for larger stroke amplitude the attenuation at high frequencies is somewhat less and follows a distinct trend. With the introduction of the new material, we decided to use the double-actuator bee for the force measurement experiments. It is important to identify the relationship between the actuator deflection and the output wing flapping angle θ wing . Now considering the Double-Actuator Bee, measuring wing angle with the high-speed camera during flapping at quasi-steady frequencies (1 Hz, 5 Hz...) allows us to determine the kinematic relationship between the input δ and the output θ wing . This linear 68 Figure 2.17: Single- vs Double-Actuator Microrobotic Bee. Raw and filtered force measurements for both microrobots. mapping is used to transform the measured input to the transmission in millimeters to approxi- mate dynamic wing angles. The results of these mappings are shown in Figure 2.16 for a smaller actuator voltage amplitude (200 V) and full-range voltage (240). In Figure 2.16(a), I am not sure how attenuation could become resonance with increasing frequencies for this system with an expected second-order-system response, but currently I can only comment on the outlier nature of the 95 and 100 Hz data points. In Figure 2.16(b), for the higher voltage, tests follow the ex- pected trend, of decreasing amplitude with frequency as we approach the end of the functional frequency range of the loaded system. The instantaneous lift force measurements of the Double-Actuator Bee are significantly dif- ferent than those of the Single-Actuator Bee, which can be clearly seen in Figure 2.17, and look much more like those in [94]. To attempt to understand certain observed phenomena in the data, we look at Figure 2.18, the measured values for displacement, voltage, and driver input voltages for a 95-Hz flapping test are compared on the same time scale, 2.18(a) for the entire duration of the experiment and 2.18(b) during the period when flapping is initialized. In 2.18(a), the long-term stable zero position is clear, showing a few random disturbances (likely from myself pushing but- tons on the driver box) as well as the reactions from the actuator being powered on and off. Note that this was driven with an old Simulink model, before we realized that the near-instantaneous voltage changes like the step in the input signal plot drastically reduce the performance and life of 69 Figure 2.18: (a) Raw voltage output from capacitive sensor, showing a stable zero and repeatable periodic flapping behavior (b) Raw displacement output from the laser distance sensor (c) Ac- tuator driver signals and their relation to the timing of visually evident phenomena in (a) and (b). the actuator, especially for the wide (3.5 mm) ones. The actuator input signals refer to the inputs of our in-house high-voltage piezo driver system, where V in =− 2: V out = 0 and 0V in ≈ 400V out . Looking closer in 2.18(b), we see the clear periodic behavior in the force readings, although the displacement sensor has significant noise. Electrical isolation, cable damping, proper DAQ set- tings, and a low-pass zero-phase filter with cutoff frequency of 1000 Hz all contributed to the smoothing of these readings over time. Looking at the relationship between input displacement and instantaneous force in figure 2.19, the periodicity is clear, with the same displacement corresponding to the same force at each repeat of the cycle. While the zero-phase filter with 1000 Hz cutoff frequency smooths the sawtooth pattern at the beginning of the first downstroke, the raw data suggests that some perturbations are introduced that are not part of the model used in [94] and in [95] (the one used in this work). This might be an effect of the slightly flexible wing as it is modeled as a rigid body and the experimental videos shows clear bending. It also might explain the extra peak per cycle in 70 Figure 2.19: Actuator Displacement vs Flapping Force the Single-Actuator Bee, as an over-rotation of the passive hinge on the second upstroke. Six tests were performed on the same microrobot, the plot color darkening for each subsequent test. It is clear that as the hinge ages, the magnitudes of the instantaneous forces go down, although the profile remains nearly identical. The first two tests involved tuning of the actuator parameters (using the high speed video as reference) and match closely, the third and fourth which were taken in quick succession (within 30 min) are very close, the fifth test was taken 4 hours later without the video lighting, and the final test was eight hours later, at which point the shape was still similar. This suggests that the Kapton hinge was possibly entering steady state. Four full datasets of average periodic lift force measurements are shown in Figure 2.21. The data suggests similar conclusions about a possible steady-state condition of the wing performance in some window of time. We can see that the new wing (shown in red) produces high forces relative to typical values from a ‘mid-life’ robot. Later that day, two sequential tests with several hours in between have very similar measurements. When the test was performed a week later after some intermediate testing, the average forces followed the same trend but at a reduced 71 Figure 2.20: Consistency of Measurements. Six 100-Hz flapping experiments were performed on the same microrobot. magnitude. This suggests that careful attention to the number of cycles that a particular wing has been tested under could help estimate the performance but also save time by having a known lifetime. This reinforces the need to have a more easily replaceable wing feature as has been described earlier in this work. As is discussed in Chapter 3, we ultimately used data from the Single-Wing Test to compute the average lift force per period because a single-flapping-wing case is more ideal in terms of internal reaction forces, out-of-phase behavior of the wings, and potential resonance and flow disturbance from the second wing. However, likely because the actuation system was designed to have both wings flapping at once, there are some unresolved efficiency losses due to forces being transmitted through the body instead of through the transmission. A comparison of the halvedtwo-wing test, theone-wing-powered test, and theone-wing-removed test is shown in Figure 2.22. The test where only one wing waspowered is very similar in shape to the halved two-wing test for the 75, 90, and 100-Hz data, markedly different than the single-wing test. As I could see 72 Figure 2.21: Lift force vs flapping frequency. Four tests of the same microrobot. The first three tests (red, green, blue) were taken in order on the same day, and the black was measured one week later. Figure 2.22: Comparison of Instantaneous Force Measurements at 75, 90 & 100 Hz flapping fre- quencies for the Halved-Two-Wing, One-Wing-Powered and One-Wing-Removed tests. sympathetic motion of the un-powered wing during video testing, I wonder if the larger ‘spikes’ are a disturbance from the reaction motion of off-side of the bee. However, the average forces for the halved two-wing and one-wing case are actually very close (and the one-wing-powered case is not far off either) therefore we decided to report the more ideal case in the paper while still pointing out the differences between the two instantaneous force measurement datasets. 73 Chapter3 TheEffectofChordwiseWingFlexibilityonMicrorobotic ForceGeneration 3.1 Methods To study the wake field and force generation at scale presents a technical challenge. While a custom force sensor was developed in [21] (and reproduced in Section 2.1) for this exact purpose, the physical size and time-scale of the flapping necessitate specialized equipment optimized for high-speed image capture for flow visualisation such as particle image velocimetry (PIV). This challenge is compounded by concomitant micro-robotic issues such as human assembly error and poorly-characterized fabrication repeatability. Therefore the following test case is proposed: a two degree-of-freedom (2-DOF) flapping and pitching mechanism is submerged in water, ob- served by a four-camera tomographic-PIV system and fit with a 6-axis force-torque sensor at the wing root. The wing design from [21] is scaled to have a similar chordwise and spanwise Reynolds number, modeled with a simplified wing planform constructed from an isotropic poly- mer material. The chordwise stiffness of the wing is quantified by the dimensionless parameter, Π 1 , the effective stiffness [17]. Π 1 has been shown to have a significant effect on the wake field [30, 96, 97], and the data in [38] suggests that there is room for improvement of the lift-drag ratio and/or propulsive power if the effective stiffness of the FW-MAV wing is reduced. Therefore a 74 Figure 3.1: End-effector of the mechanical flapping and pitching system, showing core variables and coordinate system. scaled wing with similarΠ 1 to the at-scale design is compared to three other wings ranging from ultra-flexible ( 10 0 ≤ O(Π 1 )≤ 10 1 ) to rigid (10 2 ≤ O(Π 1 )≤ 10 3 ). The fluid flow field of the FW-MAV undergoing sinusoidal flapping and pitching motion in a hovering condition is studied with a dynamically-scaled model in water. This permits a reduction in the characteristic velocity, using commercial stepper motors to drive the system and regular- frame-rate cameras to capture data, while maintaining a mathematically similar flow state to the original state. In this case, we consider the Reynolds number of the wing from two distinct characteristic lengths, the span b and the mean chord c m . Given a flapping frequency f with stroke amplitudeΦ A for the system shown in Figure 3.1, we can formulate the chordwise Reynolds number with the reference velocity as the mean wingtip velocity,U ref = 4Φ A fb/2 as in [17]; Re c = ρLU ref µ = 2ρfΦ A c m b µ , (3.1) and the spanwise Reynolds number withU g as the cycle-mean velocity at the radius of gyration R g as in [39], Re b = ρLU ref µ = ρbU g µ = 4ρfΦ A bR g µ . (3.2) While these two equations use different reference velocities, the two speeds are quite similar 75 as R g is not far from the midspan. This difference is a deliberate attempt to keep the reported values for Re c and Re b defined consistently with the literature from which they are taken. In work such as [22], Re b has been used to identify wing parameter spaces that exhibit unique flow regimes in the wake and characterize the performance of certain designs and prescribed motion patterns. The design and exact scale of the microrobotic bee has evolved over time at the various institutions where related research is being conducted, but the dimensions presented here match the robot from previous chapters in this work. With a wing length b= 14 mm and mean chord c m = 3.9 mm, these FW-MAVs have chordwise Reynolds number equal to 761 and spanwise Reynolds number equal to 2726. As the aim of this work is to study the effect of chordwise wing flexibility on the wake, Re c is more representative of the parameter of interest. It is interesting to note that the FW-MAVs in [4, 25, 43] are most often generally classified as having Reynolds numbers O(10 2 − 10 3 ), but their spanwise Re is several times higher and into the transitional regime, which could mean that the combined effect might produce flow fields that were not predicted for the design. However, this work focuses on mimicking the existing micro-robotic wing designs without delving into the deeper meaning of the difference in Reynolds number classification of the wings by characteristic length. The addition of a force-torque sensor in the scaled model presents a technical challenge to be able to preserve both Re b and Re c , as the additional offset of the sensor from the sweep axis disproportionately increases Re b . However, this is not inherently undesirable as the microrobotic bee design also has a relatively large offset between the sweep axis and the wing center of pres- sure [24]. Choosing the mean wing chord as the characteristic length, in keeping with the main goal of this work, we select the other experimental parameters to mitigate the magnitude of Re b . Conducting the experiment in water as opposed to air changes the density and viscosity of the medium, therefore the similarity conditions are satisfied by scaling the chord length, span, am- plitude and frequency of the flapping. 76 Figure 3.2: Example solution spaces for the dimensionless parameter similarity problem, given a representative span magnification factor of 3. For a simple case where the sweep axis is at the wing root and the scaled span is three times its original length, it is clear from Figure 3.2 that the scaled frequency of flapping needs to be at least two orders of magnitude smaller than the original case for the wing size to be large enough for standard laboratory instrumentation to drive the system and capture data. While commercially- available macro lenses are able to record images at the length scale of the wing from the original design parameters, high-frequency flapping motion is difficult to achieve with common motors (stepper motors, servos, etc) as well as observe/record with traditional flow capture software frame rates. Due to the dichotomy between chordwise and spanwise Re, the scaling equations require some amount of iterative calculations to arrive at a set of system parameters that result in a scaled model that reasonably represents the original system. Attempting to limit the experimental Re b to that of the microrobotic bee (2726) results in extremely short wings that would not produce enough forces to register on the force sensor at the desired low operating frequency. With the gearbox design and the offset of the tool face of the force-torque sensor at 46.5 mm from the sweep axis, Re b increases quickly with increasing span length. Additionally, without a specified planform profile, R g and thus U g cannot be accurately calculated. Therefore, we solve for the scaled wing length using a target value of Re b and placeholder variables f (representing thatR g 77 is some algebraic scaling factor that is proportional to the total wing length), the force sensor offset b 0 and true wing length,b. The reference velocity is thus U g = 4Φ A fR g = 4Φ A f· b 0 +s f (b− b 0 ) (3.3) and the target Re is as in [39, 98], equivalent to Equation 3.2 where the kinematic viscosityν is equal to the dynamic viscosity divided by the fluid density ρ, Re b target = U g b ν . (3.4) Rearranging equation 3.4 forU g and setting it equal to Equation 3.3, we arrive at a second-order polynomial of the formαb 2 +βb+γ = 0 with coefficients α = 2Φ A fs f (3.5) β = 2Φ A f(1− s f ) (3.6) γ =νRe b target . (3.7) The roots of this equation are the possible dimensions for sweep axis to tip length (b 0 +b) that satisfy the target Reynolds number. In this case, the target Reynolds number was set to be 3700 and the scaling factor was set to be 0.72 as in [24]. The scaled experiment flapping frequency was reduced to 0.2 Hz to accommodate the standard stepper motors and the constraints demonstrated in Figure 3.2. The peak-to-peak flap amplitude was halved from Φ A = 120 ◦ toΦ A = 60 ◦ as Re b is sensitive to Φ A and the larger angular travel distance places too strong a restriction on the 78 quantityb 0 +b for a fixed offset b 0 . The resultant estimated wing length,b ′ =b− b 0 = 56.9 mm, was subsequently used in the scaling equations for chordwise Reynolds number similarity: Re c 1 =Re c 2 (3.8) ρ 1 c m 1 U ref 1 µ 1 = ρ 2 c m 2 U ref 2 µ 2 (3.9) 2ρ 1 c m 1 Φ A 1 f 1 b 1 µ 1 = ρ 2 c m 2 Φ A 2 f 2 b 0 + b ′ 2 µ 2 (3.10) c m 1 µ 2 ρ 1 Φ A 1 f 1 b 1 2µ 1 ρ 2 Φ A 2 f 2 b 0 + b ′ 2 =c m 2 . (3.11) This equality resolves to a mean chord length of c m = 24.8 mm; a simple wing planform can be constructed that satisfies this mean value while keeping the span b ′ less than roots of the polynomial equation. Through trial and error in the design, the wing profile was determined to be a half circle of radiusr= 24.5 mm, centered ath= 8 mm below the leading edge of the wing, defined by x=r cos(θ)+b 0 +r (3.12) c(x)=r sin(θ)− h (3.13) wherex is the length along the span andc(x) is the chord profile. This contour is easily integrated from 0≤ θ ≤ 2π to verify that c m is 23.6 mm and the span length b ′ is 48.5 mm. With this simplified geometry the radius of gyration is straightforward to compute from the moment of inertia about the sweep axis,I z and the total areaA TOT , breaking the shape into the partial circle and rectangular segments, with lengthl= 2r; I z = 1 12 (hl 3 )+(lh)· (b 0 +b ′ /2) 2 + 1 8 (πr 4 )+(πr 2 /2)· (b 0 +b ′ /2) 2 (3.14) A TOT =(lh)· πr 2 /2 (3.15) R g = r I z A TOT . (3.16) 79 Figure 3.3: a) At-scale FW-MAV wing, with section view of at the mean chord (shown with the red dotted line). The relative lengths were used to calculate the weighted average effective stiffness, Π 1 . b) Scaled designs with equalc m , (L): transversely anisotropic (R) isotropic This scaling process approximately maintains the chordwise and spanwise Reynolds numbers of the original system. The swept arclength normalized by the mean chord was reasonably pre- served, reducing from 3.75 to 3.20 for the experiment in water. Note that the pitching amplitude Ψ A does not factor into any of the scaling equations therefore the value can be preserved from the original MAV system. Depending on the hinge design, the passive pitching amplitude ranges from 45− 70 ◦ , where the coefficient of lift is maximum for the lower limit but the lift-to-drag ratio C L/C D increases towards the upper limit [99]. The value ofΨ A = 45 ◦ was selected to study the maximum lift case and to be more directly comparable to literature such as [30, 36, 96, 100]. Finally, it was not possible to preserve the wing aspect ratio,A R , which is known to have an effect on the flow field and lift production [38, 101]. This parameter was kept within reasonable bounds for bio-inspired wings, having been reduced from 3.72 to 2. A third dimensionless parameter is used to characterize the response of the transversely- loaded thin plate wing flapping in a fluid. The effective stiffness Π 1 is the ratio of elastic bending forces along the chord to the fluid dynamic forces [102], where bending stiffness coefficient D s 80 Table 3.1: Scaling and dimensionless Parameters (Π-Variables) of the system. f c m R g Φ A A R Re b Re c FW-MAV 100 Hz 3.7 mm 0.0118 mm 120 ◦ 3.7 2726 761 Scaled 0.2 Hz 23.5 mm 0.0716 mm 60 ◦ 2 3030 675 Π 1 h s Material FW-MAV 313 65 µm Mylar and Carbon Fiber Wing 1 5 50.8 µm Mylar Wing 2 81 127 µm Mylar Wing 3 256 254 µm Polycarbonate Wing 4 870 381 µm Polycarbonate is calculated from wing thickness h s , elastic modulus E and Poisson’s ratioν using an isotropic plate model. This dimensionless quantity can be expressed for isotropic materials as Π 1 = D s ρU 2 ref c 3 m = Eh 3 s 12[1− ν 2 ]ρU 2 ref c 3 m isotropic plate , (3.17) whereρ is the fluid density and U ref is the same reference velocity as in Equations 3.1 and 3.2 [17, 38]. There is no clear method to calculate the effective stiffness of the original composite wing due to the asymmetric planform and spar design (Figure 3.3(a)) as well as the anisotropic [0-45- 0] unidirectional carbon fiber layup. The polymer is extremely thin with an unsupported free end and it is not practical to attempt to load a wing like a cantilever to measure deflection. The thickness relative to the wing size is also an obstacle to producing reliable finite-element analysis (FEA) simulations that can approximate bending stiffness. Therefore the value is approximated with a weighted-average approach, considering the cross sectional area of the wing at the mean chord. This assumes that additivity applies and that the materials that make up the composite are operating in the elastic region which is reasonable for the loads expected in this case. As this experiment will compare wings across the flexibility spectrum, the comparison of the results are still valid despite the exact unknown stiffness of the FW-MAV wing. 81 From the designs in [21, 24] there are two wing spars that cross the mean chord, the areas of the mylar and carbon fiber areas are computed as A mylar =b· h mylar = 14(0.0065)= 0.091 mm 2 (3.18) A CF = 2h CF · w spar = 2(0.044)(0.135)= 0.0108 mm 2 (3.19) Therefore the estimatedΠ 1 for the at-scale wing is found using the weighted average W = Σ n i=1 w i X i Σ n i=1 w i (3.20) such that the combined effective stiffness is Π 1 = A m A CF +A m E mylar h 3 mylar 12[1− ν 2 mylar ]ρ fluid U 2 ref c 3 m ! + A CF A CF +A m E CF h 3 CF 12[1− ν 2 CF ]ρ fluid U 2 ref c 3 m ! , (3.21) which evaluates to Π 1 = 313, corresponding to the fully rigid regime defined in [38]. This is consistent with FW-MAV literature that designs the wing planforms to have a high degree of stiffness [22]. A design improvement should target an effective stiffness in the range Π 1 = O(10 0 ): O(10 1 ), when the elastic forces are close in magnitude to the fluid forces, shown by [38] to increase lift and reduce the mean drag. To design the isotropic wings with a specified Π 1 for further comparison, Equation 3.17 is solved forh s such that h s = 3 s 12Π 1 [1− ν 2 ]ρ fluid U 2 ref c 3 m E , (3.22) which allows for a range of stiffnesses to be achieved via appropriate material selection. The parameter space of h s and E is constrained by the requirement that the material be optically transparent for one of the two flow visualization configurations (see Section 3.1.3) as well as com- mercial availability of polymer sheet thicknesses. Generic polycarbonate and polyester (DuPont 82 Mylar ® ) were selected for their optical properties as well as their differences in elastic modulus and Poisson’s ratio, in addition to their suitability to be bonded to carbon fiber with epoxy resin. Working within these constraints, the scaled wing that is most similar to that of the FW-MAV is Wing 3 in table 3.1, with an effective stiffness of 257. While this is 18 % lower than the at- scale wing estimate ofΠ 1 , the weighted-average did not account for the angled wing spar which would reduce the stiffness in the chordwise direction. This wing is compared to three others, a set that spans the flexibility range as defined by [38], from low- Π 1 values that represent extremely flexible wings to high- Π 1 , rigid wings that are nominally analogous to the wing commonly used for flapping wing research. While the effective stiffness of the scaled wing (Wing 3) technically falls under the rigid regime in [38], the planform was observed to deform during experiments and so it was decided to compare it against an even more rigid specimen, the Wing 4 presented in Table 3.1. 3.1.1 ExperimentalPlatform To test these new designs, a milli-Newton-resolution, 6-axis force-torque (F/T) sensor (ATI In- dustrial Automation Nano17 Titanium IP68) attached to a single wing in a symmetric cm-scale flapping-wing mechanism is driven to both periodically flap (here Φ) and pitch (Ψ) (see Figure 3.3) with a predefined sinusoidal input. The 2-DOF mechanical flapper is submerged in the 3 m 3 water tank shown in Figure 3.4 and driven by two stepper motors through a coaxial drive shaft and 3D-printed gearbox. Concurrent with the force measurement, time-resolved 3D wake information is captured with the LaVision DaVis Flowmaster Tomographic-PIV equipment. This volumetric PIV is conducted for two distinct laser volume configurations- a horizontal volume beneath the trailing edge to capture the downwash as well as a vertical volume that bisects the sweep angle at the midstroke to evaluate the development of the LEV and TEV. The mechanical system is designed to have a minimal impact on the fluid inside the tank besides the desired wing motion. As much as possible, the components are kept above the surface of the water on a rigid structure made from extruded aluminum t-slot channel material. The 83 Figure 3.4: Experimental 4-camera tomographic-PIV system. (Left) “Horizontal" laser configura- tion to capture downwash. (Right) “Vertical" laser configuration to capture cross-section of wing chord and developing spanwise vorticity. schematic of the assembly is shown in Figure 3.5. Two NEMA23 stepper motors with integrated drivers are linked to a coaxial drive shaft with 3 : 2 reducing pulley pairs, one for each of the degrees of freedom of the flapping system. This drive shaft is in fact an appropriated tail assembly of an RC helicopter, the Align T-Rex 250, as the carbon fiber torque transmitter (center axle) and its hardware is high quality and conveniently off-the-shelf. The drive shaft extends into the gearbox, where the internal shaft is fixed to a bevel gear which turns the wing pitching axis. This gearset is a 16/64T custom pair designed in Siemens NX software following [103] and fabricated with 3D printing technology. The 16T pinion was made on a SLA-type FormLabs Form 3 from V3 clear resin and the 64T crown gear was made on a PolyJet3D-type Stratasys Connex3 Objet 260 from Vero white. It was observed that making at least one of the gears out of Vero resulted in quieter and smoother operation of the gearbox, despite the lower XY resolution of the Objet printer compared to the Form 3. The crown gear is fastened to the central pitching shaft/axis with set screws, transmitting the motion from the drive shaft and bevel gear to the F/T sensor and the wing root. These components are all housed in an open gearbox made from structural parts printed on the Form 3 printer (except for bearings, standoffs and hardware.) The flapping and pitching mechanism was designed to meet several physical constraints be- yond generating the prescribed motion, which in itself was not straightforward. Aside from the 84 Figure 3.5: Computer-Aided Design (CAD) model of the core components in the flapping system. obvious need for waterproofing, the gearbox needs to be submerged as far as reasonably pos- sible into the water tank to minimize the wake’s interaction with the surface which can cause force disturbances as well as illumination distortions in the PIV frames. However, as the coaxial drive shaft gets longer, the system becomes more flexible and vibrations can become a serious issue. An additional problem is that the force sensor has a delicate cable that exits the mounting plate axially and has a minimum bending radius of 3 inches, contraindicating any sort of design that mounts the sensor on one side of the gearbox only. Therefore the cable must gothrough the bearings and gearbox and out the other side, which completely determines the gearbox size scale, as the end connector of the F/T sensor is 17 mm in diameter. This constraint is the reason that there is a channel in the pitching axle/shaft, to safely route the sensor cable out and behind the gearbox. The kinematics for the system are shown in Figure 3.6, derived from Equations 6 and 7. Slices of the wing chord at mid-span are shown at selected values of Ψ, t/T . A single flapping pe- riod is generally considered to be mid-stroke to mid-stroke, highlighted with the red background between 0≤ t/T ≤ 1. The wing velocity is at its maximum at this time instance; initiating the motion at this speed is not practical with the stepper motors. Therefore flapping is started at 85 Figure 3.6: Flapping kinematic profiles for the prescribed motion. Flapping angle Φ(t) and pitch- ing angleΨ(t) are directly specified by the output of the motor control program, with amplitudes of 30 and 45°, respectively. The normalized time axis represents 5 sec per period given the input frequency of 0.2 Hz. Flapping speed ˙ Φ(t) is included for reference, it is not directly specified in the system. Projections of the pitching wing as viewed from the lab reference frame are included above the plot to familiarize the reader with later plots. time− 0.25t/T , where the sweep angle is maximum (Φ(− 0.25t/T)=Φ A ) and the velocity is zero. This approach has the additional benefit of reducing mechanical noise and startup effects on the data for Cycle 1. Three computers are used to carry out this experiment; one for timed motion control, one for force data acquisition, and one for PIV image capture and laser control. A simplified diagram is shown in Figure 3.7. The motor PC uses a National Instruments NI-DAQ PCIe 6321 to generate digital outputs for the clock, step and direction signals to the two NEMA 23 stepper motors, micro- stepped to output 0.225°/step. This is further reduced by the pulley and gearset to 0.15°/step and 0.0375°/step for flapping and pitching, respectively. The amplitude of the sinusoidal motion functions are converted into functions of motor step count, and the array of digital motor control signals is stored in the DAQ card buffer where it can be deployed in real-time. For the coupled mechanical system in the gearbox, when the stepper motors are engaged and holding position, the act of rotating the flapping axis induces a pitch angle as the outer shaft moves around the torque rod. Therefore, to achieve the desired pitching profile, the flapping function is added to the pitching function and the mechanism moves as desired. A unified code for this procedure was written in MATLAB and has since been shared among graduate students in related laboratories due to its general applicability to prescribed motor motion. 86 Figure 3.7: Simplified schematic of the entire experimental setup. Three Microsoft Windows desktops with two National Instruments NI-DAQ input-output (I/O) breakout boards and four LaVision FrameGrabber camera cards, controlling: a) 24V power supply b) NEMA 23 stepper motors c) Sub-miniature S.P.D.T. limit switches d) Force-torque sensor e) Optical isolator and related circuitry for limit switches f) LaVision Imager sCMOS cameras g) PIVision pulsed Nd:Yag laser. The direction signal for the flapping stepper motor is split outside the motor PC’s SCB-68 DAQ breakout board and routed into the data acquisition computer. This computer has a NI-DAQ PCI 6289 card, which is capable of recording 18-bit analog data at high sampling frequencies (up to 500 kS/s multi-channel) while also leaving channels open for other inputs such as the motor trigger signal. This allows for the proper synchronization of the measured force/torque data and the prescribed motion across the two devices. The final subsystem in this experimental setup is the LaVision PIV computer. This device, in tandem with an external timing unit, synchronizes the image capture of the cameras’ field of view with the illumination of the double-pulsed laser. The PIV system is capable of recording double-frame images at frequencies up to 28 Hz for the double-frame∆t used for these experiments but the recording frequency was selected as 20 Hz to maintain a reasonable expected particle displacement between sequential correlations (explained further in Section 3.1.3). Therefore knowing the kinematic motion frequency, the length of the flapping period in terms of number of frames is known a priori. The image acquisition is initiated before the motors begin motion and the PIV data is synchronized manually with the force/torque data based on the frame in which motion commences. The timing is further verified by repeating 87 Figure 3.8: (Left): Removable wing base that attaches to the tool plate of the ATI F/T sensor, showing the coordinate system of the sensor. (Right): Front view of assembly in the left image, showing the sensor rotated by pitch angle Ψ. Lift L and drag D are shown in the laboratory reference frame. this process for the second cycle, ascertaining that the wing returns to the same position T/f frames later. 3.1.2 ForceBalance The force sensor is located at the wing root, and thus rotates and translates with the wing. The system has no absolute zero position or rotation in the lab reference plane, but centering the forces and torques on the cardinal sensing axes is crucial for subsequent analysis of the signal. In the same vein, the installed limit switches on the stepper motor shafts must be given a reference to be useful. The zero-position is determined from data from kinematic tests where the force sensor is installed without any end effector. As the symmetrical motion should produce symmetrical data for an unloaded system, the step count from an arbitrarily-placed switch trip is adjusted until the F x,y,z signals are symmetric. Data from this moving reference frame must be transformed to be able to calculate coefficients of liftC L and dragC D over the flapping period. Here the same method as in [39, 98] is used, given the coordinate system shown in Figure 3.8. As the wing chord is aligned with the ˆ x direction of the sensor, we denote this asF c , a positive magnitude pointing in the positive lift direction relative to the wing. The wing-normal direction is in the ˆ y direction of the sensor, which is further referred 88 to asF n . Given a wing pitch angle ofψ radians, lift and drag forces in Newtons can be calculated from the instantaneous data as L=F n sin(ψ)+F c cos(ψ) (3.23) D=F n cos(ψ)− F c sin(ψ). (3.24) The lift and drag coefficients can be computed from these values as C L = 2L ρU 2 g S (3.25) C D = 2D ρU 2 g S , (3.26) whereρ is the density of water in kg/m 3 ,U g is the reference velocityU g = 4Φ A fR g (Equation 3.3) in m/s, andS is the surface area of the wing planform in m 2 . Dynamic force sensing is achieved with the ATI Industrial Automation Nano17 Titanium IP68, the smallest commercially-available six-axis force-torque transducer in the world. The SI-50-0.5 calibration has forceF x,y,z resolutions of 1 /80 N and torqueT x,y,z resolutions of 1 /16 Nmm and can be further reduced with dynamic signal filtering. The six-channel raw voltage data is acquired with a National Instruments PCI-6289 DAQ card at a sampling frequency of f s = 10 kHz and converted to SI units with a 6× 6 calibration matrix provided by the manufacturer. The signal is tared at the start of the motion cycle when the flapping mechanism is at rest. The TTL trigger signal from the motor control program on a second computer is routed into the DAQ to time the periodic force traces, shown in Figure 3.9. There is significant noise in the raw sensor data, visible in Figure 3.9. This is believed to be from the mechanical system itself, such as gear meshing and the transmittance of vibration from the stepper motors. As the equations presented here for lift and drag (Equations 3.23 and 89 Figure 3.9: Six cycles of unfiltered F/T data, converted to S.I. units with the calibration matrix and separated into forces F x , F y , F z in Newtons and torques T x , T y , T z in Newton-millimeters. The TTL direction signal for the flapping stepper motor is shown in the third subplot; the length of one square wave plus the time interval until the next pulse is equal to the period (5 seconds.) 3.24) depend solely on the forces in the chordwise and chord-normal direction), we look at the frequency response of these two signals via the fast Fourier transform (FFT) in Figure 3.10. As expected, the fundamental frequency occurs at the flapping frequency f Φ = 0.2Hz for both F c and F n , with harmonics at the algebraic multiples of f Φ . Both plots show the presence of the high-frequency noise that is apparent in Figure 3.9, the vast majority of which is in a frequency range more than two orders of magnitude higher than the fundamental. To filter the data with minimal distortion and phase lag, a Butterworth filter is applied to the data with zero-phase digital filtering using MATLAB’s filtfilt() command. Due to the frequency separation between the fundamental and the noise, a third-order filter is sufficient. Two cutoff frequencies f c are compared in Figure 3.11, each at least ten times the fundamental flap frequency which is considered the minimum value for filtering flapping-wing data [21, 26]. For the force traces, the f c = 10 Hz still shows remnants of noise. This is even more apparent in the torque time traces for the same f c . The cutoff frequency of ten times the fundamental, f c = 2 Hz, retains a similar signal magnitude while showing an exceptional reduction in noise for both force and torque traces. As subsequent mathematical manipulation of this signal is required, the lower cutoff frequency is selected to process all of the F/T data to prevent fluctuations in the result. 90 Figure 3.10: Frequency response of the two F/T sensor signals used in the calculations for C L and C D ; F x = F c and F y = F n , respectively. The system response of the raw periodic data was calculated using the Fast Fourier Transform in MATLAB. Figure 3.11: (Left): Three periods of force and torque data, filtered with a third-order Butterworth filter and cutoff frequency f c = 10Hz. (Right): Three periods of force and torque data, filtered with a third-order Butterworth filter and cutoff frequency f c = 2Hz. All data was filtered in MATLAB. 91 Figure 3.12: System frequency response for filtered F c andF n . Resonance peaks for low-frequency components have remained at their original values and the high-frequency noise has been sup- pressed. Butterworth filter magnitude response is plotted as an overlay in red. The filtered system shows attenuation of nearly all of the high frequency noise, demonstrated in the FFT in Figure 3.12 with the overlay of the filter transfer function. The F/T data as it has been presented up to this point cannot be directly used to computeC L andC D as it includes inertial components from the translating and rotating system [104]. The force sensor in particular has a large mass relative to the rest of the components and hence the inertial effects from the kinematics must be removed. As the wing reference frame experiences acceleration, the simplest way to remove this added force is to deploy the kinematics with no wing attached and subtract this zero-case from the force traces measured with the wing in place. However, there are inertial effects from the mass of the wing itself which also contribute to this overall quantity. Without access to a vacuum chamber, this quantity cannot be measured directly [105]. Here, a similar approach to [12, 26] is taken, where the force traces are recorded with the wing mounting plate and carbon fiber leading edge installed but without the wing planform. The CF rod is responsible for 25-45% of the wing mass, respectively, for Wings 1-4 in this work. 92 Figure 3.13: Inertial subtraction of non-wing components, shown for Wing 1. (Left): Force dif- ference between wing- and non-wing data∆F c and∆F n for the two subtraction methods, force sensor only vs force sensor with C.F. rod leading edge, averaged across 5 tests and shown over 3 periods of motion. (Right): Coefficients C L andC D by inertial subtraction method. Note that the major characteristics of this plot remain constant for all wing types; the difference in C L is predominantly found in the range 0.25≤ t/T ≤ 0.75 and thatC D calculated from the null sub- traction is larger for 0≤ t/T≤ 0.65. This method is compared with the data produced by the zero-case with only the wing mount- ing plate attached. Considering the quantities∆F c and∆F n , the differences between wing force trace data and the zero-case, a systematic difference is visible between the two subtraction meth- ods, shown in Figure 3.13 from data averaged over five tests. In the chordwise direction, it is reasonable that the difference ∆F c is smaller in magnitude for the case with the rod attached, as the presence of the rod produces a higher inertial force than if it were absent. In the wing-normal direction,∆F n differs only during the reverse stroke (integer multiples of 0.25≤ t/T≤ 0.75). The force trace for the CF rod is more symmetric about zero which is a logical consequence of sym- metric kinematics. The effect of this difference can be seen in the C L andC D plots on the right hand side of Figure 3.13. The traces deviate the most in the range 0.25<t/T < 0.75 and the periodic multiples of that time range, where the case with no CF rod produces lower lift and higher drag. Considering the 93 asymmetry in the∆F n plot for the null case; if a negative bias of 7% of the amplitude is applied, the signal becomes symmetric about zero. Subsequently, theC L by subtraction method plot exhibits differences over the whole cycle (instead of just the reverse stroke) except for where the velocity is zero (t/T = 0.25,0.75) which makes physical sense. However, applying an arbitrary bias to only one of the signals when all of the data had been processed and averaged in an identical manner was not determined to be appropriate. Neither case produced a mean drag coefficient of zero, which would be the expected output of a symmetric and reciprocating flapping motion, but this may be due to imperfections in the wing fabrication. The wing itself is not symmetric, as the carbon fiber leading edge is adhered with epoxy to one of the surfaces of the wing material. The wing was centered as much as possible with respect to the pitching axis but inherent asymmetry may be a factor in the measured differences for the different motion directions. Ultimately the C.F. rod subtraction method was selected for its similarities to examples in the literature. 3.1.3 ParticleImageVelocimetry(PIV) Direct observation and measurement provides one of the three lynchpins for understanding com- plex flow physics, together with numerical simulation, and analytical/mathematical models [106]. Particle Image Velocimetry (PIV) is a non-intrusive optical technique for fluid dynamical research that calculates vector fields from images of a fluid system. The availability of these velocity vector fields makes it possible to analyze complicated flow dynamics in terms of familiar mathematical quantities such as force, momentum, and pressure. Despite the several different specialized sub-types of PIV, the underlying principles remain the same: a medium is seeded with neutrally-buoyant tracer particles that are assumed to faithfully follow the fluid flow, illuminated with a light source such as a laser. A pair of images of this fluid is captured, separated by some time interval∆t, such that there is an observable particle motion between the frames. Both images are divided into small interrogation windows. Let us call such a window in Image 1 B, of dimension B x × B y pixels. If the fluid is in motion, we expect to be 94 Figure 3.14: The (simplified) process for tomographic-PIV, shown here with the minimum 3 cam- eras instead of the 4 used in these experiments. Image pairs from each camera are recorded from an illuminated fluid volume with tracer particles, the images are used to create a tomographic reconstruction, and the PIV cross-correlation procedure is conducted on the reconstruction to generate a 3D vector field. find the displaced particle in Image 2 att=∆t within the search boxS of sizeS x = 2u max ∆t+B x by S y = 2v max ∆t+B y . The variance-normalized correlationC(i, j) (Equation 3.27) is calculated for these two areas, using decreasing box sizes [107]. The displacementD for each interrogation window is determined from the location of the maximum correlation peak in the window, from which the velocity is directly calculable with∆t. C(i, j)= Σ B x k=1 Σ B y l=1 I A (k,l)− I A I B (k+i,l+ j)− I B q Σ B x k=1 Σ B y l=1 I A (k,l)− I A 2 Σ B x k=1 Σ B y l=1 I B (k+i,l+ j)− I B 2 (3.27) The concepts explained in the previous paragraph are for a single image-pair, producing a single 2D vector field. This approach can be extended to a case using two cameras and stereo-PIV, which produces a single 3D vector field. Currently, the most advanced iteration of this method is Tomographic-PIV, which employs at least three cameras and a laser illuminationvolume. The data from the multiple cameras is combined in a tomographic reconstruction [108] for each double- pulsed frame at each time instance of the recording. This process, shown in Figure 3.14, computes the cross-correlation between the reconstructed volumes as opposed to images, and the result is a 3D vector field. 95 3.1.4 PracticalConsiderationsforTomographic-PIV Tomographic-PIV requires a calibration error less than 0.1 pixel [109] which is generally not achievable with planar calibration tools. The initial calibration is obtained with images of a dual-plane metal plate with a known grid pattern (LaVision Type 11 or 204-15 depending on the camera configuration) which maps pixels to millimeters. The Volume Self-Calibration (VSC) tool (included with the LaVision DaVis 10 software) is used to correct the initial calibration [110]. A user-guided algorithm uses a sequence of images of the particles in motion to create a disparity map for selected depth planes within the measurement volume and fits a polynomial function to the data points. This procedure generates an estimate of the position of each particle in the measurement field based on the measured intensities from all cameras together and back-projects this position onto each sensor reference frame[110]. The difference between the observed peak intensity loca- tion and the projected one creates a disparity vector for each identified particle. As this process is meant to correct the general mapping function, the field of view is discretized into subvol- umes, so the updated, localized correction function can vary based on local data instead of the entire measurement volume. Iterating the VSC procedure while increasing the maximum allowed number of particles to be detected produces clusters of disparity values for each subvolume that collectively converge to a global average value less than 0.1 pixels in magnitude. This means that for all cameras and frames, the final disparity map will look like a grid with the disparity peak in the center of each subvolume and the disparity vectors will not appear randomly oriented in space. In this case more subvolumes is not always better; there must be a sufficient number of particles in each subvolume to contribute to the Gaussian probability density function disparity maps so the user can determine when convergence has been achieved. When this process is done correctly with data from a well-conducted experiment, VSC can reduce the fit error of the initial calibration over two orders of magnitude. The act of preparing a sample set of particle images for VSC is an excellent quality check to record suitable data for tomographic-PIV itself. If the algorithm fails to converge with an 96 acceptable pixel error, it is an indication that something is wrong with the physical system. A common issue is illumination and level of contrast in the camera frame. It is clear from Equation 3.27 that the intensity of the particles in the image should be high, maximizing the dynamic range of the camera sensor (in this case, 16-bit pixel resolution) while minimizing the background. In a similar vein, the aperture of the camera lenses are often set to the lowest f-stop possible to maximize the light, but this causes the a reduction in depth of field. This issue paired with the required inter-camera angles (up to 30 deg [108]) and improper use of tilt-shift/Scheimpflug lens adapters can produce particle images that are sufficiently out of focus such that the algorithm will not converge. For laser volumes several centimeters thick, special attention must be paid to the image focus outside of the center of the field of view which may be near the boundaries of the physical focus limits of the specific lens hardware in use. To address these physical constraints,a 3D parametric CAD model was constructed in Siemens NX of the tilt-shift adapters as well as the focus region of the lenses to ensure that the camera placement and hardware settings resulted in a full depth of field over the entire laser volume. Once a suitable corrected calibration is achieved, the reconstruction and correlation process shown in Figure 3.14 is implemented. However, this process does not have a time correlation from one cross-correlated vector field to the next, which would be clearly beneficial for a time series of data from a deterministic, physical fluid flow. Therefore the LaVision Motion Tracking Enhance- ment (MTE) tool is used to reduce “ghost particle" intensity which is uncorrelated in time as per the original process [111]. Ghost particles are reconstruction artifacts that do not actually exist in the measurement volume and thus their pseudo-intensity raises the noise floor[108]. MTE itera- tively generates enhanced initial guesses for the reconstruction using time and velocity data from multiple reconstructions which reduces abrupt changes between cross-correlated time instances. As this tool requires several (3− 6 depending on user choice) reconstructions per double-frame, the computational cost scales with the number of MTE frame comparisons. Given that it is time- intensive to process Tomographic-PIV even without MTE (for this 4-camera, double-frame data, 97 Figure 3.15: Z-Intensity profiles of the reconstructed measurement volumes, calculated with the built-in software tool in LaVision DaVis. The red and green plotted values are the intensities from the first and second laser pulses, plotted against the z-depth relative to the zero established during planar calibration. a) Z-intensity profile from the spanwise flow experiments with 23.2 mm of laser depth with SNR>2. b) Z-intensity profile from the downwash experiments with 12.6 mm of laser depth with SNR>2. about 15 min per double-frame on a 16-core computer), as many quality checks as possible should be conducted before committing to the processing time for a complete experimental run. Besides ensuring that the fit error is below 0.1 pixels, a “Z-Intensity Profile" may be calcu- lated as a measure of the signal to noise ratio of matched particles vs. likely ghost artifacts that will occur at the extremes of the volume nearest to-, and furthest from the camera image plane. This process involves computing a tomographic reconstruction (no correlation) across the mea- surement volume depth, deliberately setting the depth limits significantly outside of the laser illumination limits. The illumination intensity is then integrated over the depth for the entire Field of View (FOV); the plot of this integral reveals the level of ghost particles in the system as well as the depth range where the Signal-to-Noise Ratio (SNR) is sufficient. Their magnitude is represented on the plots as the intensity value outside of the laser illumination region. Therefore the SNR is the intensity divided by the ghost particle intensity. Tomographic-PIV vector fields are considered of good quality for SNR> 2 [108]. Figure 3.15 shows the Z-intensity profile for the two experimental runs presented later in this chapter. As the laser does not have a knife-edge filter, there is no sharp cutoff between 98 Table 3.2: Summary of volume self-calibration polynomial fit errors for downwash and spanwise experimental configurations. Spanwise/Vertical Downwash/Horizontal z-Plane [mm] Fit Error [px] F.E. [mm] z-Plane [mm] Fit Error [px] F.E. [mm] -15.5 0.0221 2.5× 10 − 3 -5.3 0.0019 1.9× 10 − 4 -6.5 0.0163 1.8× 10 − 3 0 0.0026 2.6× 10 − 4 2.5 0.0211 2.4× 10 − 3 5.3 0.0035 3.6× 10 − 4 illuminated and dark regions in the fluid. While this is not the preferred state, there is still a significant percentage of the laser volume that has sufficient SNR. Both of these cases meet the minimum requirements for SNR in some volume region, but it is clear that the case in Figure 3.15.b had a superior calibration. The peak intensity is smaller but the noise floor is near zero, resulting in a large region of usable data within the laser volume. The explanation for the difference is likely the laser divergence and power setting during the experiment. With a fixed-power laser source, spreading the laser to a wider volume illumination comes with the obvious cost of loss of coherence and local intensity. Additionally, for the horizontal arrangement (Figure 3.15(b)), the laser volume is entirely below the flapping mechanism and reflections are at a relative minimum. For the vertical arrangement in Figure 3.15(a), the wing is passing through the illuminated fluid, which causes large reflections several times per period, which can cause permanent damage to the CMOS sensors in the cameras. Therefore the pulse power is reduced and the overall experimental illumination is lower, essentially guaranteeing a lower correlation value for the PIV results. The difference is quantified in Table 3.2, showing that the experiment for the horizontal/downwash flow data has a calibration fit error an order of magnitude smaller than the vertical/spanwise flow data. However, fit error of the self-calibration is not a direct measure of the uncertainty in the sys- tem. Due to the correlation procedure, every PIV vector is an estimate including some error, and without a quantified uncertainty analysis, the results remain qualitative. For a well-constructed and executed experiment where illumination, focus, and contrast are expected to be sufficient, the error in the PIV calculations is expected to be predominantly a function of random errors as opposed to bias errors. Factors contributing to random error include the particle number density 99 N (usually given as particles per pixel, PPP), particle diameter in pixels,d, box sizeB and search radius,S (Equation 3.27) [107]. In practice, uncertainty calculations for experimental PIV data are derived from statistics-based methods that rely on either multiple observations of a given flow field or synthetic data for comparison [112]. In this case, neither approach is applicable. From the practical guidelines in [113] a likely uncertainty in PIV displacements is about 0.1 pixels. Since this quantity is fixed in pixels, the uncertainty as a fraction of estimated velocities varies based on the experimental parameters. Given an average velocityu in mms − 1 and calibra- tion mapping of m in pixels per millimeter, the average displacement is(um∆t) pixels, where∆t is the image-pair exposure time. Since the uncertainty in timet is much smaller than in spacex, the relative uncertainty can be approximated as ∆u u ≈ ∆x x = ∆x um∆t . (3.28) For the experiments in this work, the mean downwash velocities are in the neighborhood of 6 mms − 1 , thus employing the calibration mapping of m= 8.84 mm /px and ∆t = 19 ms results in mean displacements of 1 pixel. From the pixel uncertainty ∆x, this corresponds to a likely uncertainty of 10% of the velocity measurement (matching the estimations in [107]). Gradient quantities such as vorticity are harder to estimate but can be expected to be at least double those of the velocity. As the flow fields for the flapping-wing system are complex, almost turbulent, the allowable length of the exposure time∆t is limited. The maximum possible∆t of 19 ms was selected given the inter-frame acquisition rate and peak observed particle velocity, therefore the relative uncertainty could not be reduced further. 3.1.5 VectorFieldSmoothing PIV-derived velocity data are first estimates of the flow state and are contaminated with noise. Some kind of smoothing, or regularization, can be expected to improve the estimates, and this can be very important when measuring gradient quantities such as vorticity. Here, we use the 100 Discrete Cosine Transform - Penalized Least Squares (DCT-PLS) algorithm presented in [114] that detects and replaces outliers in addition to smoothing in the three coordinate dimensions of the data. The smoothed vector field ˆ U is solved for with, ˆ U= IDCT2 Γ◦ DCT2 W m ◦ W bs ◦ U− ˆ U + ˆ U , (3.29) where IDCT2 and DCT2 are the type-II discrete cosine transform and inverse, respectively. The MATLAB functions dct(X,N) and idct(X,N) are adapted in [114] to compute N-dimensional transforms and are used here to solve this equality with fixed-point iteration given an initial value ˆ U=U. Masked data are accounted for with the matrixW m and outliers are suppressed with the bisquare weight matrixW bs of studentized residuals. The filtering matrix Γ Γ kl = 1+S 4− 2cos (k− 1)π) m − 2cos (l− 1)π n 2 ! − 1 (3.30) contains the scalar smoothing parameter S, which can either be selected through the minimiza- tion of the generalized cross-validation score (CGV) [115] or specified directly as an input option to the smoothn() function. S is allowed to vary between 0 (data unchanged) to∞. To show an example in 1-D, two polynomial functions are corrupted with Gaussian noise of 50, 25, and 15% of the signal amplitude. The robust smoothing algorithm selects smoothing values and returns an output very similar to the original signal, as shown in Figure 3.16 for a sampling rate of f s = 200 Hz. The smaller the magnitude of the noise, smaller values ofS are required to drive the signal to- wards the original value. For well-performed experiments with minimal noise, a low smoothing value implies that the underlying data contained only a small percentage of outliers and that the 3D vector field was generally smooth for all components. Additionally, the smoothed output improves with an increase in f s , a natural effect of the increase of data points for the algorithm. This is highlighted in Figure 3.17, a zoomed-in view of the partial Fourier approxima- tion of a square wave from Figure 3.16, where f s is increased to 10 kHz with a Gaussian noise 101 Figure 3.16: Test signals “corrupted" with Gaussian noise to demonstrate the robust smoothing algorithm. Left top: f(t)= sin(t)+ sin(3t)/3+ sin(5t)/5+ sin(7t)/7+ sin(9t)/9. Left bottom: f(t)=e − 0.15t sin(3t). The smoothed output is shown for both functions in the right column. amplitude of 15% of the signal amplitude. However, it can be seen that for this noisy signal to be smoothed to a near-approximation of the original, the value of S is on the order of 10 8 . This smoothing value is indeed large, but noise of 15% of the signal is also significantly higher than expected from a well-conducted PIV experiment. To compare this information to results on actual PIV flow data, various smoothing spline interpolants were tested on a representative three-dimensional dataset from a downwash exper- iment and their performance can be seen in 1D cuts in Figure 3.18. The velocity componentsu,v, andw are normalized by the maximum value over the flapping period and plotted against planar coordinates x,y normalized by the mean wing chord length. The smoothness of the line plots of normalized u versus y and v versus x will determine the smoothness of the calculated out-of- plane vorticity. The raw data are plotted against the fully robust result where the algorithm was allowed to self-select a value of S for each frame, as well as three other cases that had specified values for smoothing. These are compared to the output of the built-in 3x3x3 smoothing vector post-process operation from the LaVision DaVis flow capture software. Even for the case of S= 3 when it could be argued that over-smoothing is occurring, the imposed smoothing parameter is 102 Figure 3.17: Demonstration of the algorithm sensitivity to sampling rate. Detailed view of part of the noisy Fourier square wave approximation signal from Figure 3.16 with Gaussian noise of 15% of the signal amplitude, using data with varying sampling rates. For f s = 10 kHz, the smoothed output is nearly coincident with the original signal. eight orders of magnitude smaller than for the best result from Figure 3.17 for high amplitude noise contamination. All smoothing outputs reasonably follow the raw data, exhibiting the expected increasing attenuation of the peaks as the smoothing parameter S increases. These peaks, such as the one seen in thev/v max plot atx/c m = 4, are believed to be from actual flow structures and not noise as each one has several data points in the local neighborhood that appear to follow a relatively smooth contour. The raw data plot markers shown in black are shown for the maximum grid density of the PIV result. All outputs of the robust algorithm reduce the local magnitude of the velocity components where the gradient is high, more so than the simple spatial smoothing filter from DaVis. For the downwash data in Figure 3.18, the smoothed estimate is rather insensitive to the free choice of smoothing parameter in u and v, though w(x) has some variation that indicates a higher level of smoothing is preferred. However, for the flow experiments that observed the vol- ume that intersects the span (the vertical configuration), the smaller robustly-selected smoothing parameter preserved the complex vorticity field with minimal alteration of the information. To be consistent between both experimental configurations, all data were sent through the same spline 103 Figure 3.18: 1D slice plots of a single plane of a tomographic-PIV dataset, comparing various smoothing parameters to the raw data as well as the built-in smoothing function from LaVision DaVis. Raw data plot markers are plotted for each grid point from the vector field output. The section of zero velocity fory/c m < 3.2 in theu/u max plot arises from the PIV masking function. 104 smoothing scheme which chose the value ofS by minimizing the GCV, with an average value of 0.05. 3.1.6 Analysis Once a 3D vector field has been generated, the velocity components u(x,y,z,t)= are specified at each grid point in the reconstructed volume space. Therefore at each time instance in the series of double-frame captures there are threem× n× p matrices that define the fluid state. The field-of-view dimensions m andn are determined by the pixel resolution of the cameras and the interrogation window overlap in the PIV algorithm; p depends on volume depth from VSC as well as the window overlap. The PIV calibration procedure provides the function that maps the volume grid in pixels into standard SI units. These variables are accessed from the proprietary .vc7 vector field files through the readimx() MATLAB function provided by LaVision which allows for further analysis of the data outside of the DaVis software. These measured quantities can be used directly to visualize directional components of the flow. Examples such as 2D and 3D quiver plots show the distribution of velocities as arrows across a gridded area or volume and provide a useful way to interpret a region of values as opposed to a matrix of numbers. It is common to normalize the velocity components, to contextualize the numbers relative to other system parameters and enhance readability at a glance. For example, normalizing someu i ,v i ,w i by the magnitude of the velocity vector,|u i | 2 = p u i 2 +v i 2 +w i 2 , pro- vides an immediate understanding of how much an individual component contributes to the total velocity at a given point. The velocity components and their spatial and temporal derivatives are then used to estimate forces and fluxes in the system. In this work, the vertical force produced by the wing wake is estimated by measuring the strength of the downwash. As the laser volume is positioned below the wing, the jet produced by the moving wing is expected to be primarily aligned with the depth coordinate axis, ˆ z. Therefore 105 the estimatorF z in Newtons is calculated from the square of the corresponding velocity compo- nentw, multiplied by the fluid density ρ and integrated over the defined wake area (see Section 3.2) on each depth plane p i as F z p i =ρ Z Z w(x,y,p i ) 2 dxdy[N]. (3.31) Given the known grid spacing of the discretized volume, this integral is easily calculated with trapezoidal numerical integration in MATLAB, nesting the function to integrate overx andy with the commandtrapz(y,trapz(x,Fz)) for each time instance over the measurement period. For physical systems, the force time series for any given plane in the measurement volume is expected to be a smooth curve with no discontinuities. In experiments, factors such as image ac- quisition rate, discrete integration method and noise will introduce fluctuations in the calculated value of F z . To extract salient data from the integral time series and present the information in a way that is likely more reflective of the true flow state, a cubic smoothing spline interpolation is used. For a set of observations{x i ,y i = 1,...,n}, the cubic spline estimate ˆ f of function f minimizes p s n ∑ i=1 |y i − f(x i )| 2 + (1− p s ) Z |f ′′ (t)| 2 dt (3.32) where p s is the smoothing parameter. The spline estimate is highly sensitive to the value of p s , particularly in the near neighborhood of the value 1 /(1+ε) [116]. The value of ε is equal to h 3 /16, where h is the average difference for all x i in the set. When p s = 1 /(1+0.01ε), ˆ f follows the underlying data closely, and when p s = 1 /(1+100ε), the spline smoothed function demonstrates moderate smoothing that is expected to still faithfully follow the underlying data. This process is implemented in MATLAB with thecsaps(x,y,p) function. Several p s -values near the critical range are tested on an example dataset generated from Equation 3.31 and plotted for part of the first flapping period in Figure 3.19. Considering the multiplication factor of the second term in Equation 3.32, for(1− p s )= 0.38 (shown as the green line) the resulting ˆ f follows the underlying 106 Figure 3.19: Selection of cubic spline smoothing parameter p s on the time historyF z plots, shown for part of cycle I of Wing 1. Data from individual frame calculations of Equation 3.31 are shown with X-markers. The legend compares values of “1− p s " to remain consistent with the form of Equation 3.32. data well without excessive oscillation. This value for the p s provides satisfactory smoothing for both cycles of all wing types and is used for allF z data in Section 3.2. For the vertical laser arrangement, the fluid flow is resolved around the surface of the wing, providing velocity data where performance-critical vorticial structures are known to appear. By Stoke’s theorem, the vorticity vectors through a surfaceS can be measured by the circulationΓ around a boundary that encloses it, such that Γ= Z Z ∇× u· dS. (3.33) Therefore the strength of a vortex structure that can be contained by a boundary in space can be quantified by Γ, which requires the calculation of the curl of the vector field ∇× u curl(u)=∇× u=ω = ω x ω y ω z = ∂w ∂y − ∂v ∂z ∂u ∂z − ∂w ∂x ∂v ∂x − ∂u ∂y , (3.34) the three components of which represent the strength of the rotational motion in the three co- ordinate directions. Considering the orientation of the laser volume relative to the wing for the 107 spanwise vorticity calculations in Section 3.1.1, the parameter of interest is the third element of Equation 3.34. The axis of rotation forω z is aligned with the span at the mid-stroke which makes it an appropriate metric to measure the strength of the leading edge and trailing edge vortices. Considering the circulation in the span-normal plane, Γ z can be calculated from ω z integrated over an XY area that contains a cross-section of a vorticial structure. The vorticity components of the vector field are calculated with the curl(x,y,z,u,v,w) command in MATLAB and the discretized form ofΓ z from Equation 3.33 is further discussed in Section 3.2. Given a 3D vector field that is large enough to encapsulate entire vorticial structures in the re- solved volume, vortex identification algorithms can be an extremely powerful tool that can clarify complex flow data and reveal phenomena that do not necessarily align with the coordinate sys- tems traditionally used to analyze slices of data. The Lambda-2 (λ 2 ) criterion determines whether a location in a flow field is part of a vortex core and defines the vortex as the connected region in which all interior points are classified as core points. This detection method decomposes the 3D velocity gradient tensorJ J≡ ∇u= ∂ x u x ∂ y u x ∂ z u x ∂ x u y ∂ y u y ∂ z u y ∂ x u z ∂ y u z ∂ z u z (3.35) into symmetric and antisymmetric components: S= J+J T 2 , Ω = J− J T 2 (3.36) and the three eigenvalues of S 2 +Ω 2 are computed for each point in the vector field. At each grid point, if two of the eigenvalues are less than zero, it is considered a vortex core. In standard eigenvalue ordering, asλ 1 ≥ λ 2 ≥ λ 3 , this is the meaning of the method name. Isosurfaces can be plotted from the connected regions with two negative eigenvalues; however, iso level threshold levels must be selected carefully as this method is prone to identifying vorticial structures within the noise floor. 108 Figure 3.20: Laser and camera arrangement relative to the wing for the horizontal/downwash (left) and vertical/spanwise (right) configurations. Example slice planes of the volumetric data are shown for both cases, with the wing location overlaid with the black line. 3.2 Results As previously described, there are two distinct experimental configurations to measure the flow data. Relative to the laboratory reference frame, the horizontal laser is positioned such that the closest illumination plane is 10 mm below the trailing edge of the wing. After volume self- calibration, the 12.6 mm volume was discretized into 20 planes in depth. This experiment was designed to capture the downwash below the hovering wing, to quantify the propulsive lift per- formance of the wing normal to the plane of observation. Therefore, the measurement region was restricted to a partial annulus centered at the projection of the sweep axis, with an additional 15° buffer added to the maximum and minimum wing angle. An example frame of PIV results is shown in Figure 3.20 with the fluid velocity magnitude as the color map. For the vertical laser arrangement in Figure 3.20, the illumination region is between normal- ized span 0.21≤ z/b≤ 0.88 during the mid-stroke. From the VSC, the 23.2 mm volume depth is discretized into 30 planes. In this configuration it is possible to visualize the vortex structure and measure circulation at the mid-span during the power stroke. At maximum and minimum stroke 109 angle, the wingtip is inside the converged laser measurement volume therefore vortex forma- tion during the pitch-overturn motion can be observed. The maximum FOV of the tomographic correlation was retained for vorticity analysis as the flow demonstrated large spatial variations. An example of this phenomena is shown in the normalized velocity magnitude slice plot in Fig- ure 3.20 with a field of view of (8× 8)c mean where the wing is returning over its own wake and creating a flow pattern several times its characteristic length scale. 3.2.1 Downwash Considering the fan-like area below the wing, we estimate the force produced by the wake asF z (Equation 3.31 in Section 3.1.6) in Newtons. This function is integrated over the partial annulus area,dS, for each depth plane p i in the measurement volume: F z p i =ρ ∑ i ∑ j w(i, j,p i ) 2 dS [N]. (3.37) Figure 3.21 shows the integral value from Equation 3.37 for 2.25T , so two full upstrokes and downstrokes are presented in addition to the start-up region. The data are plotted for the center plane of the measurement volume as well as an average of all z-planes. The plot color scheme has the line color becoming darker with increasing wing stiffness. The data was captured at an acquisition rate of 20 Hz for flap frequency of 0.2 Hz, therefore each plot trace is derived from a set of 225 data points. These sets are spline-smoothed across the time interval using cubic spline fitting to enhance plot clarity with a minimal p-value to prevent over-smoothing. As the closest plane of the measurement volume is approximately half the mean chord length below the wing, it is reasonable that F z ≈ 0 at the mid-volume for− 0.25<t/T < 0 while the fluid moves towards the observation window. For both the mid-plane and the mean, the force estimator from the stiffest wing increases earlier in the period than the flexible wings but reaches a lower local maximum. The data are repeatable between comparable cycles; (aside from the 110 Figure 3.21: IntegratedF z for four wing stiffnesses over two flapping periods, plotted for − 0.25≤ t/T≤ 2, cubic spline-smoothed from 225 frames of data. The line plot colors are darker for stiffer wings. (a): F z for the mid-volume plane of the downwash experiment. (b): F z averaged across the entire measurement volume. 111 Table 3.3: Average force estimator values in milliNewtons for four wing types across two periods of downwash data. Mean Force EstimatorF z [mN] Wing Cycle I Cycle II 1 1.8± 0.1 1.3± 0.5 2 1.9± 0.6 2.0± 0.4 3 2.8± 0.6 3.1± 0.7 4 2.4± 0.7 2.9± 0.5 start-up region), 0.25<t/T < 0.75 and 1.25<t/T < 1.75 on the mid-plane plot could be the beginnings of periodicity in the wake integral. The most flexible case, Wing 1, has the most notable difference in instantaneous magnitudes for both plots. Additionally, this wing data exhibits fluctuations on a much shorter time scale than the other wings, which could indicate some non-negligible Fluid-Structure Interaction (FSI). The mean downwash data are summarized in Table 3.3 for the first two cycles of flapping. The uncertainty in the variation of the integral across depth planes is estimated by the precision, the maximum difference between the mean value over the volume and any single value of F z from a given plane. Although the mean value of F z in the volume increases with wing stiffness for both cycles I & II, up to 33% and 120%, respectively, given the variation by plane it is not possible to conclu- sively say that the downwash magnitude is different. For cycle I, the uncertainty ranges mostly overlap; for cycle II, the most flexible wing appears to be an outlier with a significantly reduced downwash. While this may be an artifact of the rapid temporal changes in the velocity field that were previously noted, the overall trends of this information match the conclusions in [28] in that more flexible wings produce a weaker downwash than that of their rigid counterparts. While the usable thickness of the laser volume was 12.6 mm (approximatelyc m /2), the length scale of the flow structures and the time scale of the frame acquisition rate made it difficult to iden- tify any other phenomena beyond the downwashw(x,y,z). Even for the beginning of the flapping period where the wing is moving through still fluid, the TEV passes through the measurement window quickly. Further analysis and/or interpretation should be guided by information deter- mined from the vertical/spanwise laser configuration and F/T data. 112 3.2.2 SpanwiseVorticityandCirculation The vortices formed by the interaction of a moving wing in fluid represent the transfer of energy from the kinematic input into velocity and momentum. The formation and persistence of these coherent flow structures are critical to lift production for flapping-and-pitching flight. As an inclined wing pushes into a fluid medium, a leading-edge vortex (LEV) forms along the wing span and wraps around the rear surface of the wing, contributing to the lift force generated. The trailing-edge vortex (TEV) forms as the surrounding fluid is pushed down and around the distal edge, rotating counter to the LEV, contributing to the drag force. As the wing undergoes complex motion, these structures demonstrate rapid changes in shape and magnitude over time, including separation into separate entities. The total lift and drag force generated by a flapping wing system depends on the size and interaction of these vorticial structures. The strength of such a vortex is measured by its circulation; for the field of view considered here where the XY plane intersects the wing normal to the span at mid-stroke, the LEV and TEV central axes are assumed to be aligned with the Z-axis, when their strength can be measured by the z-component of Equation 3.33; Γ z = Z Z ω z dxdy. (3.38) From the perspective of a slice plane along the span, such as is the case for the vertical laser configuration, cross-sections of these vortices should be clearly evident in the PIV-derived vector fields. Thus, in a 2D slice of the illumination volume, the circulation Γ of some region R can be estimated with the discrete approximation Γ R =Σ i Σ j ω z (i, j) dxdy. (3.39) However, as all PIV data has noise, it is undesirable to include this quantity in the total sum. Therefore the data should be trimmed according to some threshold level or percentage of the maximum observed vorticity. Several such contour levels are shown in Figure 3.22, ranging from 20 to 60% of the global maximum vorticity from the set of four wings over two periods of flapping. 113 Figure 3.22: Threshold levels for spanwise vorticity as a percentage of the global maximumω z , shown for four wings at normalized time oft/T = 0.5. The black line overlaid on the plot is the approximate location of the chord line of the wing at the z-plane of observation, mapped from the 3D kinematics as a point cloud in MATLAB. Note that the chord line is plotted as a linear projection and wing bending is not included. This procedure was performed in the same manner as in [117], including ultimately selecting the same threshold of 0.2ω max . This choice of threshold value does retain some smaller vorticial structures away from the wing which is most visible in Wings 1 and 2 in Figure 3.22. Selected time instances of the thresholded data are shown for the mid-span plane in Figures 3.23 and 3.24 with the chord line overlay as in Figure 3.22. The velocity vector field is shown at 1 /4 the true grid density (3.2 mm spacing.) The colormap of these time series shows the spanwise vorticity, presented over two periods of flapping motion for each of the four wings. Note that the wing and fluid are already in motion at t/T = 0; as explained in Section 3.1.2, the experiment starts att/T =− 0.25 atΦ max due to physical motor constraints. As the fluid is initially still and the motor starts at slow speed, the first 25 frames do not show much of particular interest. 114 From the data in the figures, the first period shows the flow structures that are generated in still water and the second period shows the interaction of the induced flow with the wake. The data were recorded at 20 Hz therefore it is not practical to display the 200 frames per wing type. The quarter-period frames were selected to highlight the vorticial structure that can be critical to flapping-wing performance. The plots for the four wing types share the same thresholding and colormap using a global maximum ω max . For the simple sinusoidal motion prescribed for this experiment, at normalized timet/T = 0.5n,n= 0,1,2... the wing is at maximum pitch angle at the mid-stroke, which is called the power stroke. At normalized timet/T = 0.25m,m= 1,2... the wing begins the overturn motion as the flap angle reaches its maximum value and the pitch angle changes sign. This kinematic transition region has been demonstrated to have a potentially large impact on the wing performance, a phenomenon that remains only partially understood [12, 15, 17, 101]. For each wing, there appears to be good repeatability between the first and second period, despite the fact that the first cycle occurs in quiescent fluid. There are qualitatively similar flow features across all wing types, such as the formation of the LEV and TEV at high pitch angles. The LEV does not shed during the translation phase; all wings exhibit delayed stall [15, 118] for both periods. For the forward power stroke att/T = 0, 1, the LEV of Wing 1 has the highest magnitude across the entire contour, remaining reasonably close to the posterior face of the wing while reaching most of the way to the trailing edge of the wing. The comparable vortex from Wing 2 has distinctly lower maximum vorticity values near the trailing edge and exhibits curvature away from the wing surface at approximately c m /3. As the wings become more rigid (with Wings 3 & 4), the LEV does not extend as far towards the TE. For the first period, both wings show a separation of the LEV into two local maxima (which was also slightly evident for Wing 2). For the second period, both shed small vortices behind the wing. Similar behavior with minor differences is seen in the reverse power stroke at t/T = 0.5, 1.5. Note that the colors of the LEV and TEV are reversed as the reciprocal motion reverses the di- rection of the induced vortices relative to the camera reference plane. For this time instance, the 115 LEV of Wing 1 is even more closely attached to the entirety of the wing for both cycles but shows a curvature not evident in the forward stroke. Wings 2, 3, & 4 show more symmetry across the period and the LEV displays reasonable matching with the beginning of the cycle. For the same time instances, there is an inverted trend with the TEV as chordwise stiffness increases. Att/T = 0, 1, Wings 2 & 3 show a completely detached shear layer behind the trailing edge. It remains attached for Wing 1 att/T = 0 but separates within a few frames. For the stiffer wings, the (mainly drag-associated) vortex becomes larger and more coherent, indicating an in- creased amount of fluid entrained by the flapping wing. All wings have an increased complexity in the vorticial structure below the wing in the second period due to interactions with the initial wake. The initial vortex shedding behavior of the four wings is shown att/T = 0.25. These struc- tures separate from the TEV and persist in the local region. Verified by a frame-by-frame analysis, the vortices shown at this time period are the same ones shown at a nearby location a quarter- period later in the t/T = 0.5 plot. Similarities between the first and second period are easily observed despite the significant difference in the complexity of the flow fields for the two cycles. These regions are columnar vortices that exist on all of the z-planes in the measurement volume and are generally aligned with the camera-normal direction. They exist irrespective of choice of ω z threshold and have been verified with isosurface plotting of |ω| and theλ 2 -criterion (swirling strength) of the raw tomographic-PIV output. Wing 1 produces three of these structures that remain intact even pastt/T = 0.75 when new ones are shed from the TEV. For Wing 2, there are two main vortices, and these interact strongly with the new wake that has higher positive vorticity than that of Wing 1. At the same period- normalized time, the wake from Wing 3 has pushed these structures even farther from the wing, despite the smaller size of the newly-shed vortices. Only one of the four columnar structures is still visible in the frame att/T = 0.75. In contrast, the starting vortex of Wing 4 has the largest diameter and spanwise vorticity. Of the three stable vortices for this wing type, the largest is situated directly under the wing when it is at maximum pitching angle. After the wing passes 116 Figure 3.23: Frames of the thresholded vorticity plots at multiples of 0.25t/T , taken at the mid- volume plane, Wings 1 and 2. 117 Figure 3.24: Frames of the thresholded vorticity plots at multiples of 0.25t/T , taken at the mid- volume plane, Wings 3 and 4. 118 over, this structure is ejected the farthest from the TE and does not lose coherence like those from Wings 2 & 3. These shed vortices have a high degree of interaction with the induced flow over the wings, as evidenced by the disparity between the first and second period vorticity plots. The reasonably repeatable yet distinct flow structures across wing types strongly suggests that the flexibility plays a complex and fundamental role in the detailed vortex dynamics. However, Figures 3.23 and 3.24 only show selected parts of the story, at discrete points in the cycle. A more quantitative analysis follows. To quantify and condense the information in the entire time series, Equation 3.39 was em- ployed to calculate the total circulation in the measurement volume.Γ was separated into regions of positive (Γ + ) and negative (Γ − ) so that the two categories did not negate each other. From the 30 z-plane slices in the volume, three representative spline-smoothed plots are shown for the two circulation types in Figure 3.25 at z/b 0 = 0.21,0.57, and 0.75. Note thatΓ fluctuates due to the reciprocating motion of the flapping wing; the LEV and TEV will switch in sign twice per period. The positive and negative circulation are maximum closest to the wing root and decrease monotonically towards the wingtip. At z/b 0 = 0.21, the integral continues to increase over the observed time window, while planes farther out along the span appear to level off and oscillate at a lower magnitude. Another interesting feature in these line plots is the distinctive behavior of Wing 3. For the downwash data the three flexible wings grouped more closely together than they did to the rigid Wing 4. In this case, Wing 3 has the lowest positive circulation out of all the wings. Considering an average of the circulation integrals across all of the time frames, this characteristic remains. Positive circulation for Wing 3 is up to 50% lower than the others, while Wings 1, 2, and 4 are generally grouped within 5%. This wing also generates smaller negative circulation than the other wings, but less than 8% with an outlier of 14% at the plane farthest from the wing root. From this general analysis, a less obviously apparent result also arises; Wing 2 consistently creates the largest magnitude negative circulation despite havingΓ + similar to the others. This difference is 119 Figure 3.25: Total circulation for the entire measurement volume, separated into positive (Γ + ) in the left column and negative (Γ − ) in the right column. The rows represent the distance from the wing root, at normalized lengthsz/b 0 = 0.21,0.57, and 0.75. maximum at the plane nearest to the wing root, 17%, and decreases to 8% at the other end of the measurement volume. In a global sense, these line plots look similar for all wings despite the observable differences in the time series plots in Figures 3.23 and 3.24. As the circulation plotted above is the summation of allΓ in the viewing window, it is clear that the location and degree of coherence of the vortices are important. 3.2.3 Force-TorqueTimeTraces To compare the lift and drag coefficients for the different wings, the first three periods of F/T data were recorded with five trials per wing type. The inertial component from the rotating gearbox and C.F. leading edge was removed with the method described in Section 3.1.2. Statistical 95 % 120 Figure 3.26: Mean values and 95% confidence intervals of coefficients of lift and drag C L andC D for each of the four wings over three periods of flapping. confidence intervals (C.I.) were generated using Student’s t-Test for small sample sizes. Given these assumptions, the coefficients of lift and drag are presented in Figure 3.26 for three flapping periods. With the uncertainty ranges, the wings appear to have very similarC L andC D for most of the cycle. An obvious difference between wings is that Wing 4, the most rigid, has significantly higher C D during the power stroke. This is particularly interesting as according to the work in [38], Wing 3 and Wing 4 are indistinguishably rigid for the amount of fluid forces expected for these characteristic velocities, and yet we see clearly different drag profiles. Wings 1 and 2 are flexible enough that some bending is expected so the frontal area of the wings is reduced during the power stroke, thus reducing drag before even considering induced flow effects. Despite careful zeroing of the system, the lift data exhibits asymmetry over the cycle. During the power stroke (t/T = 0.5n, n= 1,2,3...), the wing reaches the same maximum pitch angle, alternating in sign/direction. However, there are more significant differences between the wings for the forward stroke as opposed to the return motion. This is shown for Wing 2 in Figure 3.27, where the averageC L andC D curves (from Figure 3.26) are plotted for the forward strokeF(t/T) and the reverse stroke R(1− t/T) for 0≤ t/T ≤ 0.5. The hysteretic behavior is not unique to 121 the first period; cycles II and III show similar behavior for the C D andC L plots. This is consistent for all wings studied here, including the change in the C L plots for the second and third cycle return stroke between 0.35≤ (1− t/T)≤ 0.5. What is initially a double-peak in the reverse power stroke smooths into a single peak near (1− t/T)= 0.45 with 20% lower instantaneous C L for the later cycles. In this case cycles II and III are nearly coincident, while the other wings show a monotonic reduction with increasing cycle number, but all cases have approximately 20% reduction in this peak magnitude by the third period. While there could reasonably be a difference between the first and second cycle due to startup effects, there is no immediate answer to the alternating peak magnitudes on a half-cycle period- icity. A similar phenomenon is evident in the drag plots, where the separation in the plots by wing type are significantly larger for the even integer multiples of the period number. The only apparent differences between sweep directions in regards to the wing itself is the reinforced lead- ing edge, which is adhered to one side of the wing surface. In these experiments, the LE is on the anterior surface of the wing for the forward stroke and the posterior surface during the re- verse stroke. The carbon fiber rod has a diameter of 0.9 mm, which may be causing extra drag or a flow disturbance in only one direction. It should be noted that this is the same way that the microrobotic wings are fabricated; for the wings in [21] the ratio of Mylar planform material thickness h to carbon fiber t LE is approximately 12%. In these experiments, Wings 1-4 have the ratios h/t LE = 5, 14, 28 and 42%, respectively, for the fixed value of t LE . The similar phenomena for varying relative spar thickness suggest this is not a factor contributing to the measured force asymmetry. Table 3.4 shows the mean values of the liftC L and dragC D coefficients as well as C ′ D , the RMS value ofC D . Note that for a symmetrically-flapping wing, we would expect that the mean drag coefficient would be zero for each whole cycle, and yet that is not the case. While the profile of the time trace of C D appears approximately symmetric to the eye, all are slightly shifted from the origin with some positive magnitude. This shift is even more pronounced when the alternate null case inertial subtraction method was used. Beyond the carbon fiber rod addressed in the 122 Figure 3.27: Forward-backward stroke asymmetry of instantaneous lift and drag coefficient plots for Wing 2. Forward motion is plotted with a solid line and reverse is plotted with a dashed line, colors becoming lighter for later periods. Table 3.4: Summary of the average values of the coefficients of lift and drag over three periods. The average precision by performance parameter type for all wings and cycles is C L :± 0.14, C D :± 0.28 andC ′ D :± 0.23. C L C D C ′ D Π 1 I II III I II III I II III Wing 1 5 2.07 2.06 1.87 0.30 0.43 0.43 3.33 3.28 3.40 Wing 2 81 1.97 1.92 1.82 0.63 0.88 0.77 3.46 3.41 3.49 Wing 3 256 1.92 1.87 1.68 0.34 0.49 0.52 3.04 3.01 3.06 Wing 4 870 1.79 1.87 1.89 0.96 1.09 1.09 3.96 3.93 3.93 previous paragraph, it is possible that this effect comes from imperfections in the fabrication of the wings; deviation of the wing leading edge (L.E.) from absolutely normal to the F/T sensor tool face could cause this offset. Despite this potential systematic error, all wings were tested in the same method and used the same inertial subtraction data, therefore the comparison of the results should remain valid. Considering the mean values, the wings have more lift capacity and are more efficient with increasing chordwise flexibility. Despite the apparent differences between the wing performances in Table 3.4, when the preci- sion, max{x i − x}, is used as a measure of uncertainty, the difference in lift performance between the wing types is no longer clear. Considering the precision of the coefficients separated by wing stiffness and by cycle number, C L of Wing 1 and Wing 2 did not have overlapping uncertainties with Wing 4 for the first cycle. All other wing comparison pairs and later cycles did not have a 123 Figure 3.28: Overlay comparison of the most rigid wing with the most flexible, Wing 4 to Wing 1, from raw PIV images. Tip deflection δ is measured directly from the image and compared with a cantilever model that assumes the drag force as a point load acting at the mean chordc m . clear distinction in lift performance. Consistent conclusions were drawn from paired t-tests; a statistically significant difference in the mean lift coefficients at a 95% confidence interval was established between Wing 1 & Wing 4 as well as Wing 2 & Wing 4 for the first flapping period. Again considering the individual precision, Wings 1 & 3 have a smaller C D than that of Wing 4, While the mean drag of Wing 2 appears to be an outlier. For the RMS drag coefficient, C ′ D , all flexible wings had lower drag than that of Wing 4. Aggregating the precision values from all cycles and wings, a general average uncertainty for each performance parameter isC L :± 0.14, C D :± 0.28 andC ′ D :± 0.23. However, using these average values, the limits of mean lift for cycle I of Wing 1 and Wing 4 are exactly coincident. No other wing out-performs the other in terms of C L . The conclusions forC D andC ′ D remain consistent with those previously drawn. Using the force traces from Figure 3.26 it is possible to estimate the expected wingtip deflec- tion for each of the wings, which can be corroborated with image frames from the raw PIV data. The illuminated wing chord for normalized timet/T = 0 where the wing is at maximum velocity 124 ∆F n [mN] δ ∆F n [mm] F C D [mN] δ C D [mm] Wing 1 4.0 10 2.3 5.7 Wing 2 4.5 0.72 2.6 0.41 Wing 3 3.1 0.14 2.2 0.10 Wing 4 4.7 0.06 3.3 0.047 Table 3.5: Expected deflection of wing trailing edge, computed with directly-measured wing- normal force F n as well as drag coefficient. Approximate measured deflection from PIV image frame is 5.3 mm . is shown in Figure 3.28, where the most flexible Wing 1 is overlaid with the image of rigid Wing 4. The cross-sectional view of Wing 1 appears to be wider because out-of-plane deformations are captured by the laser volume. Assuming that the load can be approximated as a point load acting at the mean chordc m , a simple cantilever beam model describes the tip deflection as δ = F d c 3 m 3EI 1+ 3c 2c m (3.40) where E is the modulus of elasticity of the wing material and I is the moment of inertia of the cross-section of the wing at c m . The expected force F d can be taken as the net force in the force/torque sensor ˆ y direction, ∆F n , as it aligns with the wing-normal direction. This quan- tity was previously identified in Figure 3.13. This calculation method is compared to the use of the general equation for lift force, Equation 3.23, using the values forC D att/T = 0 from Figure 3.26, summarized in Table 3.5. Given that Wings 2-4 are considered rigid by [38], it is not surprising to calculate expected deflections less than 3% of the mean chord length. While Wing 3 appears to be somewhat of an anomaly with lower values for F n andC D , the expected tip deflection would remain similar to the reported values even if these variables were comparable with those of the other wings. Depending on the calculation method, Wing 1 is expected to deflect between 17-24% of the mean chord length. This is confirmed by measuring the image in Figure 3.28; using the PIV calibration, the tip displacement normal to the wing surface was computed as 5.3 mm or 0.22c m . The smaller 125 Figure 3.29: Wing kinematics versus normalized time, plotted in tandem with the mean signals from Figure 3.26. Only the C.I. of Wing 4 is included. Notable time instances where significant differences between the wing performances are marked for further reference. deflections could not be visually confirmed from the image frames as they were too small to be calculated from image pixels. Returning to the data from Figure 3.26, the mean values of the coefficients are re-plotted in Figure 3.29 with only the stiffest wing’s confidence interval. This allows a clearer view of the instantaneous differences between the wing types, shown against the wing kinematics plot for reference. Several time instances have been selected for further investigation, where differences in the lift and drag traces indicate that there may be observable differences in the flow structure. At normalized timet/T = 0.4, the force data plot shows similar coefficients of drag but dis- tinct coefficients of lift, of which the rigid wing (Wing 4) is the smallest. The quiver plot with a normalized out-of-plane vorticity colormap is plotted for all wings in Figure 3.30 at the mid- volume/midspan plane. Vorticity is normalized asω z c mean/v RMS , the RMS velocity of the wing at mid-span. From the plots, the separated shear layer at the leading edge of Wing 4 is distinctly less coherent than those on the flexible wings which may be responsible for this reduction in C L . 126 Figure 3.30: Normalized vorticity component ω z c mean/v RMS at normalized time t/T = 0.4 for all wings. At normalized timet/T = 0.85, we have the reverse situation, where all wings have similar C L , and the ultra-flexible Wing 1 has the lowest drag coefficient. As shown in Figure 3.31, the flow field around Wing 1 has the most coherent structures, which may be responsible for this phenomena. It could be that the more coherent structures are more effective at moving the wake away from the wing, but this cannot be determined from this frame alone. 3.2.4 VortexVisualization This unsteady flow is highly complex and the plots in Figures 3.23 and 3.24 are only a single depth plane out of thirty for each moment in time. From the tomographic-PIV data, Figure 3.32 shows vorticity magnitude isosurfaces of the LEV and TEV for the four wings att/T = 0, showing the structures that develop at the beginning of the period. Wing 1 has a LEV that extends almost to the trailing edge at large z/b 0 . The attached TEV is small, showing detached structures of the same rotational direction that have already separated from the main vortex. The LEV of Wing 2 is thicker than that of Wing 1, but also has strong rotational velocity in the chordwise direction 127 Figure 3.31: Normalized vorticity componentω z c mean/v RMS at normalized timet/T = 0.85 for all wings. between the LEV and the wing surface. The TEV has a less contiguous contour overall but has larger vorticial structures. In comparison, the LEV of Wing 3 is centered nearer to the leading edge / pitching axis as opposed to behind the planform as with Wings 1 and 2. This structure is separated from the wing surface by a coherent flow structure with strong ω x andω y components, likely reducing the upwards suction that is normally produced by the LEV [119]. Wing 4 has similar structures that exist between the main LEV and the rear wing surface but they are less coherent and do not seem to align with a cardinal axis. Wing 4 also has a large TEV that is mostly a single, connected isosurface that wraps around from behind the trailing edge and continues up the forward face of the wing. This information corroborates the observations from the center slice plane time series in Figures 3.23 and 3.24, in that for increasing chordwise stiffness, the TEV becomes larger and more coherent during the power stroke, likely the cause of the increasingC D . Att/T = 0.4, the same time instance as Figure 3.30, the isosurface plots in Figure 3.33 verify the existence of the vorticial structures seen in the slice plots. In this case, theλ 2 vortex criterion is used to identify and define the isosurfaces, a measure of the swirling strength of the fluid in 128 Figure 3.32: Vorticity magnitude|ω| isosurfaces for all wings, contour limit of 0.2815 S, ordered in increasing stiffness with wing number in upper left corner of each sub-panel. Colormap cor- responds to out-of-plane vorticityω z , direction into or out of the page. 129 Figure 3.33: 2D projection of swirling strength 3D (λ 2 -criterion) isosurfaces fort/T = 0.4, equidis- tant contours between− 0.03112≤ S≤− 0.00896, subplots numbered according to wing num- ber. Symmetric colormap of out-of-plane vorticity − 0.55≤ ω z ≤ 0.55, representing a global ± max(ω z ) over the measurement period for all wings. Approximate overlay of wing chord po- sition shown with dotted lines. motion. As this process will capture many smaller structures including what we would consider to be within the noise floor, we present the isosurfaces with an equidistant threshold to remove the smaller and possibly erroneous structures. In Figure 3.33, the wing is returning over its wake, visualized with λ 2 -isosurfaces in Figure 3.33. The columnar structures seen in the slice plots are clearly visible underneath the wing and it is immediately apparent from the 3D information how vorticial components besides ω z are influencing the fluid state. For Wings 2 & 3, the central axis of the final shed vortex from the forward stroke is no longer mostly aligned normal to the view pane. This vortex generated by Wing 2 develops a strongω x component that points in front of the wing (to the left in this case). For Wing 3, the final vortex is smaller, with the penultimate close behind. The influence of a ω x component in the former may explain why these two columnar structures seem to wrap around 130 Figure 3.34: 2D projection of swirling strength 3D (λ 2 -criterion) isosurfaces fort/T = 0.6, equidis- tant contours between− 0.03112≤ S 2 ≤− 0.00896, subplots numbered according to wing num- ber. Symmetric colormap of out-of-plane vorticity − 0.55≤ ω z ≤ 0.55, representing a global ± max(ω z ) over the measurement period for all wings. Approximate overlay of wing chord po- sition shown with dotted lines. each other. For Wing 4, behind the separated shear layer of the LEV and closest to the wing root, there exists a vorticial structure strongly aligned with the wing chord which may be responsible for the breakup of the lift-promoting LEV. Att/T = 0.6, just past the mid-stroke, these same shed vortices appear to regain coherence and alignment with an axis normal to the depth plane slices for all wings despite the large amount of swirling behind the wing. The most flexible wing generates the largest TEV which is unusual but is consistent with the force data showing the highest drag at this moment in time. The shed vortices for Wing 1 have the highest degree of coherence out of all the wings but remain closer to the trailing edge- it is possible that this flow interacts strongly with the TEV. The flow field only becomes more complicated with time. By the second period, despite strong similarities in the force trace data, isosurfaces of vorticity magnitude and λ 2 -criterion 131 Figure 3.35: 2D projection of swirling strength 3D (λ 2 -criterion) isosurfaces fort/T = 0.4, equidis- tant contours between− 0.03112≤ S≤− 0.00896, subplots numbered according to wing num- ber. Symmetric colormap of out-of-plane vorticity − 0.55≤ ω z ≤ 0.55, representing a global ± max(ω z ) over the measurement period for all wings. Approximate overlay of wing chord po- sition shown with dotted lines. differ greatly. An example is shown in Figure 3.35. At t/T = 1.1, all wings have similar lift coeffi- cients but up to a 33% reduction in drag coefficient relative to the rigid case. As this time instance is just past the mid-stroke, pitching angle and rotational velocity is near maximum, and the re- duction in wing frontal area is expected to reduce drag [38]. However, the fairly inflexible Wing 3 out-performs its flexible counterparts. This may be related to the highly coherent vortex that forms underneath the wing as shown in Figure 3.35, as well as the small and highly-separated TEV. 132 3.3 Discussion Results from the force-torque sensor and two Tomographic-PIV experiments can be combined to give a more complete picture of the effect of chordwise wing flexibility on force and wake gener- ation. Considering the force estimator over the wake area, stiffer wings produce more downwash than their flexible counterparts. As effective stiffness Π 1 increases, the wake travels into the mea- surement volume at an earlier reference time, which indicates a higher initial bulk velocity. This is particularly interesting compared to the three more flexible wings, as the peaks of the estimator integral F z reach local maxima with a time interval of 0.5t/T - the expected difference between power strokes. For Wing 4, the first peak has almost 3t/4T separation from the second, which is followed by the third at the expected interval. This suggests that there is some sort of physical flow effect in the first cycle that is not present for the other wings. This may be related to the strong starting vortex, causingF z to increase more quickly than for the other wings. The downwash was investigated over a partial annulus region that extended± 15 ◦ beyond the maximum/minimum sweep angle of± 30 ◦ , which by construction did not include the u and v components of the flow. Lateral velocity components would obviously change the direction of the jet, which can in turn effect significant changes in the lift-drag relationship for the flapping wing. In addition, chordwise bending of the wing induces changes in the normal vector for the wing surface near the trailing edge, which is expected to change the flow direction [17]. This was confirmed by the spanwise tomographic-PIV experiments. Vorticial structures that separate from the TEV were shown to have cross-sectional normalized diameters of 0.42c m (Figures 3.23 and 3.24) which is 80% of the converged tomographic-PIV volume for the downwash experiment. With their central rotation axis mostly in-plane to the depth plane slices, the w-component of these vortices can skewF z when only part of the moving structure falls within the measurement volume at a given time point. These data support a general argument that the force estimator parameter alone is not enough to relate the wing stiffness to performance. The monotonic increase in F z with chordwise stiffness for the fluid jet below the wing is consistent with [28] but is in direct contrast to the trends in the force-torque data. The force 133 traces show that each cycle has a distinct C L , C D and C ′ D . The fluid is known to be unsteady and the startup effects are significant, therefore it is not unusual to see this change as the flow develops. This reduction in mean lift coefficient with cycle number is consistent with data from an A R = 2 wing started in quiescent fluid in [101]. From the summary in Table 3.4, the reported values forC L are in a similar range to published values for sinusoidal flapping motion [12, 26, 28, 30, 33, 39, 96, 101]. The flexible wings show a decrease in C L with the period number whereas the most rigid wing has an inverse trend. From the three periods of data for the lift and drag coefficients, the difference between the wing types is most pronounced for Cycle I. Wings 1, 2, and 3 have, respectively, 16, 10, and 7% higher mean coefficients of lift for the first flapping cycle relative to that of Wing 4. This corresponds to a small percentage increase inC L for more flexible wings. The most flexible wing in this study with Π 1 = 5 corresponded to the previously-identified ideal flexibility range in [38], therefore these results are in agreement with the literature that a specific range of Π 1 can be beneficial [30, 38, 96, 97]. However, employing the precision as an estimation of uncertainty, the moderately small in- crease inC L may be eclipsed by the variation of the reported values. Despite this, the uncertainties inC L for Wing 1 and Wing 2 for the first period do not overlap with that of Wing 4 and statistical significance was established at a 95% confidence level that the mean coefficient of lift for wings withΠ 1 = 5 andΠ 1 = 81 is higher than that of a wing withΠ 1 = 870 for the same cycle. The averages are grouped closer together by the third cycle and no longer maintain the monotonic progression, but it is not unreasonable to suggest that microrobotic flight is more often oper- ating in unsteady and dynamic motion as opposed to continuous hovering and therefore more weight should be given to the performance of the earlier cycles. A notable corollary from this data is that the upper limit ofΠ 1 for beneficial flexibility is higher than the O(10 1 ) limit proposed in [38], although an exact number is not known due to the lack of wing test cases in the range 82<Π 1 < 256. The wing pitching motion in [38] is just a step function that switches direction at 134 the maximum and minimum stroke angle, therefore the flow dynamics due to continuous pitch- ing kinematics during the stroke reversal may be responsible for the increase in lift shown here [12, 15, 101]. From the drag force time series, the mean drag coefficient C D is shown to be lower for all flexible wings, though there is no clear pattern associated with reducing chordwise stiffness. This is in agreement with [26, 28, 33, 120, 121] that flexible wings reduce the mean drag. While the drag reduction can be at least partly attributed to the streamlining of the wing area during the translational segment of the flapping period [28], there are likely other fluid dynamic mechanisms that account for drag reductions when the bending is small relative to the wing chord [122] or we would not see such a reduction for the Wing 3 (Π 1 = 256) in this work. As was discussed in Section 3.2, the reason for the non-zeroC D could not be conclusively determined, but the wings were tested and analyzed in a consistent manner with a common inertial subtraction. Using either the signal meanC D or the measure of signal powerC ′ D =C D RMS , the data shows that decreasing chordwise wing stiffness reduces drag. For Cycle I, Wings 1, 2, and 3 have a mean drag reduction of 69, 34 and 65% of that of the rigid Wing 4. These reductions fluctuate for Cycles II and III but maintain the relative ordering ofC D 1 <C D 3 <C D 2 <C D 4 . Comparatively, C ′ D reductions for the flexible wings stays very similar for all three cycles, evaluating to 16, 13, and 23% less than of that of Wing 4. Therefore while the most flexible wing has the lowest mean drag coefficient, the semi-rigid Wing 3 has the lowest RMS drag, and C ′ D3 <C ′ D1 <C ′ D2 <C ′ D4 for all cycles. Considering the reported uncertainties in the measurements, the data support the fact that wings withΠ 1 = 5 andΠ 1 = 256 have lower mean drag compared with an equivalent wing withΠ 1 = 870, but a more firm conclusion cannot be drawn from the overlapping uncertainty betweenΠ 1 = 81 andΠ 1 = 870. For the RMS of the drag coefficient, wings with 5≤ Π 1 ≤ 256 have a reducedC ′ D compared with that ofΠ 1 = 870. This information supports an argument for a connection between chordwise wing flexibility and lift production as well as drag reduction. As the analysis thus far has not revealed any funda- mental changes in the flow structure, it must be that the cumulative sum of small changes in the 135 velocity fields are causing this difference. This may be related to one or more of the phenomena so far identified such as shed vortex size, LEV suction behind the wing [119], and persistence in the proximity of the wing path [120]. Differences in the flow structure for the different test cases are most easily seen in the span- wise vorticity plots from the vertical tomographic-PIV data. The normalized vorticity time series in Figures 3.23 and 3.24 show a repeatable difference in size and coherence of the LEV and TEV of the flapping wings for different chordwise stiffness. As flexibility increases, the LEV is larger and centered more closely to the posterior surface of the wing while stretching further in the chord- wise direction, suggested by [31, 119] to enhance lift via suction. Compared with the large TEV of the rigid wing, the same structure is diminished and/or fully separated for the more flexible wings which may partially explain the drag reduction. The four wings exhibit distinct vortex shedding behavior behind the TEV, with different num- bers of core spanwise vortices that persist into subsequent phases of the flapping cycle. This change in shedding behavior is consistent with [120] for tip deflections that are a significant per- centage of the mean chord length, but are surprising for small differences in expected δ such as between Wing 3 and Wing 4. For the thresholded vorticity data, the stiffest wing produces fewer and larger vorticial structures, the starting vortex of which is the strongest. It is the likely cause of the highest downwash and peak instantaneous lift forces of Wing 4 [28], potentially due to wake capture during the translational motion [39]. It separates at the latest point in the period compared with the other wings, translating slowly in the XY plane, and is positioned directly below the wing during the power stroke. Wing 3 produces similar but smaller structures that detach earlier and form a chain, which also persists under the wing during the power stroke. The more flexible Wing 2 has less coherent and fewer standing vortices over this time period, neither of which are under the wing. This may explain the larger shear layer of the TEV and higher C D andC ′ D . The most flexible wing generates stable vortical structures that persist the longest during the cycle, even after the wing has returned over its own wake. At several critical time 136 frames during the period, the flow from Wing 1 displays a higher degree of coherence than the other wings which may be responsible for the largest mean lift coefficient. As was expected from this complex, unsteady flow, the fluid motion is highly three-dimensional and simplified analysis techniques will always give an incomplete view of the system. Overall, the isosurface visualization of the induced flow from the flapping wings can put the slice-wise analysis of the data in context but does not provide compelling evidence of fluid dynamical phe- nomena that depend on the chordwise wing flexibility. These results have highlighted aspects of the flapping wing dynamics that warrant further study- repeatability of the flow evolution and transients, size and degree of coherence of the stable structures shed from the moving wing, interaction of these structures with the developing flow on the wing, among others. Different approaches to data post-processing (as opposed to smoothing) such as proper orthogonal de- composition (POD) or dynamic mode decomposition (DMD) may be able to distinguish between contributing factors in this complex flow. The results from the force/torque data and the spanwise tomographic-PIV experiment suggest that a flexible wing would be beneficial to wing performance as compared to its rigid counter- part. Information from the downwash experiment shows a larger force estimator for stiffer wings, where the two stiffest wings presented here had the same mean and planar estimations. Given the comparative ease of information extraction from these types of experiments compared to the micro-robotic bee at scale, several suggestions are presented to reinforce the conclusions in this work. Force/torque data should be collected from more than five tests to narrow the statistical confidence interval on the wing types and potentially identify more time instances where critical fluid flow differences may exist. More than one wing specimen for each value of effective stiff- ness should be built and tested to determine the F/T data sensitivity to fabrication imperfections. Multiple experiment iterations should be captured with tomographic-PIV to assess repeatability as well as establish averaged results to make a credible velocity uncertainty estimation. 137 Figure 3.36: Comparison of the coefficient of lift C L over one flapping period from the micro- robotic bee measurements and the dynamically-scaled model, Wing 3,Π 1 = 256. At-scale data is presented in raw form as well as filtered with a zero-phase digital filter, the scaled model data is filtered in the same method. Both cases have a cutoff frequency of ten times the fundamental flapping frequency. 3.4 ComparisontotheFW-MAV The question remains of how the results relate to the original micro-robotic system. The design recommendations are only useful if it can be established that the dynamically-scaled model rea- sonably mimics the at-scale robot. Considering a representative trace of the measured forces of the micro-robot from Chapter 2 in Figure 2.19, the lift force is used in Equation 3.25 to calculate the instantaneous lift coefficient C L over a single flapping period. Both the raw and filtered C L curves are computed for the data acquired on the microNewton-scale force sensor from Chapter 2. These time histories are compared to the mean and 95 % confidence interval of C L from the first period of data from scaled Wing 3 ( Π 1 = 256) which is shown in Figure 3.36. The two lift coefficient traces show qualitatively similar behavior such as the overall shape and periodicity of the data, with high positive instantaneous C L during the power strokes in both flapping directions. Over the cycle period T, the mean C L for the micro-robotic case is 1.68 compared to the value of 1.92 from Wing 3. The oscillations in the raw data have a dominant frequency of approximately 1 kHz, which are not likely to be aerodynamic in origin. As discussed 138 at length in previous chapters, the coupled subsystems in the micro-robot introduce additional dynamics that may cause unintended vibrations. However, several obvious differences are also present. Despite the reasonable magnitude of the mean coefficient of lift, the data collected from the micro-robot itself results in calculated instantaneousC L values that are substantially higher than anything reported in the literature [12, 26, 28, 30, 33, 39, 96, 101]. This may be explained by the coexistence of mechanical modes with the aerodynamic ones, and that the in-situ measurement of the microrobotic forces is dominated by an inertial component despite the extremely small mass of the wing. As the dynamically- scaled experiment in this work used an inertial subtraction procedure, this data is more likely to be representative of the aerodynamic conditions only. In addition, the micro-robot produced negative lift during the stroke reversals compared to the model case. This phenomena occurred often over dozens of the flapping-wing MAV prototypes, such as in [21] and [94]. The work in [101] suggests that it is the relatively high velocity of the wing pitch kinematics during the overturn motion that may cause this downward suction, therefore the controlled pitching of the model wings in this work may have preventedC L from becoming negative. While this does highlight a difference in the current work from the behavior of the FW-MAV, the pitching is controllable and therefore these experiments can be conducted with varying pitch overturn velocities to mimic the true motion of the passive pitching of the micro-robotic system. At the same time, the design and performance of the passive hinge system is still an active field of research, and it may yet be determined that stiffer hinges and lower overturn velocities are more desirable. A selection of several lift coefficient time history plots are shown in Figure 3.37, several of which were used as points of comparison forC L in this work. For high angles of attack (Ψ 0 ≥ 45 ◦ ), the data show cycle-mean lift coefficients 1.5≤ C L ≤ 2.5, with peak values up to 4. The average and maximum values are consistent with the dynamically-scaled model previously presented as opposed to the calculated values from the FW-MAV. With the exception of the constant-sweep- velocity experiments in Figure 3.37(g), theC L curves from sinusoidal flapping experiments exhibit 139 Figure 3.37: Selected coefficient of lift plots from current literature. a) Sinusoidal flapping and pitching [96] b) Lift force - sinusoidal flapping with fixed pitch. Note that in this case they con- sidered T to be a peak-to-peak single stroke [38] c) Sinusoidal flapping and pitching [28] d) Sys- tematic variation of flap and pitch function (harmonic) [98] e) Constant flap velocity and static pitch angle (half cycle) [26] f) Sinusoidal flapping and pitching [30] g) Constant sweep velocity and static pitch angle [12] h) Sinusoidal flapping motion with modified sinusoidal pitching, max angle of attack held constant when sweep velocity is high [33] i) Sinusoidal flapping and pitch- ing, an investigation into the pitch overturn velocity’s effect on lift coefficient for [101] (1) fast overturn (2) moderate overturn (3) slow overturn motion. 140 qualitative similarities to each other and to the data in Figure 3.36. The new data presented here does not display the high magnitude peaks in Figure 3.37(i) or the in the fast-overturn test in Figure 3.37(d) (fuchsia line), but rather two regular periodic global maximum values that repeat for each period. In general, the data from these types of experiments are difficult to interpret and compare to contemporary work as theC L vs t/T plots are sensitive to the large parameter space, especially the wing kinematics. This point is especially clear from Figure 3.37(d), where the effect of varying the shape of the harmonic flap and pitch functions dramatically change the C L time history. The clear difference in the shape of the plotted data highlights an unresolved problem in the literature; what exactly is comparable and how does one use this information to make informed design decisions? Are the differences in the force traces between the scaled experiments and the in-situ measurements too fundamental to compare? While the results presented thus far is promising, this research must be continued to attempt to address these questions. 3.5 Summaryandconclusionsforflappingmodel In this work, I have created a dynamically-scaled experimental platform that matches the chord- wise and spanwise Reynolds numbers of sub-gram flapping-wing micro air vehicles (FW-MAVs). This system enables the study of isolated parameters of the flapping-wing system such as input kinematics and wing design. This greatly reduces the functional overhead for testing new designs and minimizes compounding errors that are inherently embedded in the micro-robotic system. Flapping and pitching motion is prescribed to eliminate the uncertainty that is common with the passive pitching hinge in the system at scale. Given the known issues with the micro-robotic wings described in Chapters 1 and 2, chordwise wing flexibility was selected as the parameter of interest. Using the dimensionless parameter effective stiffness, Π 1 , the original wing design was scaled with a dynamic model and compared against three others in the range 5≤ Π 1 ≤ 870. 141 Considering the expected fluid forces on the wing, this domain ranges from an extremely flexible case to a nominally completely rigid one. Features in the force time traces have been observed that can be related to the flow structure, thus improving our understanding of the fluid dynamic mechanisms that govern this unsteady system. Time-varying, 3D flow field and force measurements show that wings with low effective stiffness can produce higher lift and lower drag than their rigid counterparts for a hovering wing case. The different degrees of chordwise flexibility create distinct flow fields during motion which may have individual and equally practical applications. These results can guide mechanical de- sign toward practical regions in the wing parameter space for more powerful (higherC L ) and/or efficient (lower C D ) micro-robotic wings. In particular, Wing 3 (Π 1 = 256) has the closest effective stiffness to the actual wings in use, which was estimated asΠ FW-MAV = 313. This value is likely even higher as the approximation here was simply a weighted average calculated at the mid-chord, but that would only strengthen the arguments presented thus far. The research completed for this work suggests that reducing theΠ 1 of FW-MAV wings has the potential to improve the design by promoting a large, well- attached LEV while reducing the size and coherence of the TEV, as well as reducing mean aerody- namic drag via wing streamlining. Given the uncertainties in the mean coefficients, wings with Π 1 ≤ 81 displayed favorable increases in the mean coefficient of lift, while wings with Π 1 ≤ 256 demonstrated reductions RMS drag coefficient. 142 OverallConclusions This work has given the supporting information for the design, fabrication and testing of micro- robots at original scale as well as individual components with dynamically-scaled models. In- structions and guidance have been given for the characterization of actuators and wing designs, as well as the method to build a low-cost and reliable micro-Newton-resolution dynamic force sensor to measure and analyze their performance. The uncertainty in the chordwise rigidity of the micro-robotic wing was observed to be a significant factor for the maximum lift generated by the flapping wing system, affected by both design choices as well as fabrication imperfections. This wing parameter was isolated in a dynamically-scaled experiment characterized by the di- mensionless variable effective stiffness, Π 1 , and studied for simple periodic sinusoidal flapping and pitching motion which mimics hovering flight. The prescribed wing pitching motion elimi- nates the passive pitching via flexure hinge that is ubiquitous in the design of the FW-MAVs, thus the results apply to the stiffness of the wing planform only. The cycle-averaged lift and drag coefficients were quantified by instantaneous force-torque data and compared with three-dimensional tomographic-PIV flow visualization and vector field analysis over the same time period. The results show that there is increased lift and reduced drag as chordwise flexibility is increased, with distinct flow fields that evolve from the fluid-structure interactions in each case investigated here. Decreasing the chordwise effective stiffness of current microrobotic wing planforms by two orders of magnitude can increase the mean coefficient of lift by up to 16% while reducing the RMS of the drag coefficient by the same percentage. The modest lift improvements for flexible wings have confirmed recent results in the literature that chordwise effective stiffnesses on the order of O(10 0 − 10 1 ) increase lift performance, although 143 this work has shown that this range can be extended closer to O(10 2 ) for symmetric sinusoidal pitching kinematics. The periodic time traces of the coefficient of lift from the dynamically-scaled model were compared to the in-situ measurements of the micro-robotic insect performance. Qualitative sim- ilarities support the initial assumption that the scaled model reasonably represents the original system, while the discrepancies highlight areas of further investigation. In particular, the induced rotation at the passive hinge of the FW-MAV likely has a large impact on the cycle-mean lift, and that careful design of its nonlinear behavior has the potential to reduce periods of negative lift. Overall, the simultaneous lift increase and drag decrease have implications for payload ca- pacity as well as power economy which are both scarce resources for the FW-MAV systems. The chordwise effective stiffness can be tuned by material selection and incorporated into the lay- ered composite fabrication process with little alteration of current processes. Most importantly, this increase in performance does not come at an additional weight cost. As all flexible wings out-performed the rigid comparison for the parameters studied in this work, there is a potential benefit to reducing the wing chordwise bending stiffness in any amount. 3.6 FutureWork The results of this work can be further confirmed following the several suggestions introduced in Section 3.3. The design recommendations can be taken directly to the micro-robotic platform at scale and implemented to test the at-scale extension of the results determined here. It should not be overly difficult to adapt the SCM fabrication method in Chapter 1 to split the wing design such that ultra-high modulus (UHM) carbon fiber is used for the leading edge while lower modulus CF or other material is used for the wing spars. Lift and drag can be measured separately using the single-axis force sensor from Chapter 2 or from an improved, dual-axis version of the same principle. 144 In addition, determining a method to appropriately quantify the effective stiffness of wings with carbon fiber spars would improve the quality of the interpretation of the results. These directional reinforcements of the planform create an anisotropic stiffness matrix that results in a more complicated deformation pattern than what was studied here. The true wing bending is thus coupled in the chordwise and spanwise direction which likely has a significant effect on the induced flow around the flapping wings. An initial proposition to move in that direction would be to do a similar experiment to what was done here, but with spanwise flexibility only. Initial data from pure chordwise and spanwise bending conditions would doubtless be useful when considering the more complex, combined case. 145 References 1. Wood, R. J. The first takeoff of a biologically inspired at-scale robotic insect. IEEE Trans. Robot.24, 341–347 (2008). 2. Pérez-Arancibia, N. O., Ma, K. Y., Galloway, K. C., Greenberg, J. D. & Wood, R. J. First Con- trolled Vertical Flight of a Biologically Inspired Microrobot. Bioinspir. Biomim. 6, 036009. doi:10.1088/1748-3182/6/3/036009 (2011). 3. Ma, K. Y., Chirarattananon, P., Fuller, S. B. & Wood, R. J. Controlled flight of a biologically inspired, insect-scale robot. Science 340, 603–607 (2013). 4. Ma, K. Y., Felton, S. M. & Wood, R. J. Design, fabrication, and modeling of the split actuator microroboticbee in2012IEEE/RSJInternationalConferenceonIntelligentRobotsandSystems (2012), 1133–1140. doi:10.1109/IROS.2012.6386192. 5. Yang, X., Chen, Y., Chang, L., Calderón, A. A. & Pérez-Arancibia, N. O. Bee+: A 95-mg four-winged insect-scale flying robot driven by twinned unimorph actuators. IEEE Robot. Autom. Lett.4, 4270–4277 (2019). 6. James, J., Iyer, V., Chukewad, Y., Gollakota, S. & Fuller, S. B.Liftoffofa190mgLaser-Powered Aerial Vehicle: The Lightest Wireless Robot to Fly in 2018 IEEE International Conference on Robotics and Automation (ICRA) (2018), 1–8. 7. Jafferis, N. T., Helbling, E. F., Karpelson, M. & Wood, R. J. Untethered flight of an insect- sized flapping-wing microscale aerial vehicle. Nature 570, 491–495 (2019). 8. Yang, X., Chang, L. & Pérez-Arancibia, N. O. An 88-milligram insect-scale autonomous crawling robot driven by a catalytic artificial muscle. Science Robotics 5 (2020). 9. Ellington, C. P. The aerodynamics of hovering insect flight. IV. Aerodynamic mechanisms. Phil. Trans. R. Soc. Lond. B305, 79–113 (1984). 10. Ellington, C. P., Van Den Berg, C., Willmott, A. P. & Thomas, A. L. Leading-edge vortices in insect flight. Nature 384, 626 (1996). 11. Liu, H., Ellington, C. P., Kawachi, K., Van Den Berg, C. & Willmott, A. P. A computational fluid dynamic study of hawkmoth hovering. J. Exp. Biol.201, 461–477 (1998). 146 12. Dickinson, M. H., Lehmann, F.-O. & Sane, S. P. Wing rotation and the aerodynamic basis of insect flight. Science 284, 1954–1960 (1999). 13. Bomphrey, R. J., Srygley, R. B., Taylor, G. K. & Thomas, A. L. Visualizing the flow around insect wings. Phys. Fluids 14, S4–S4 (2002). 14. Sun, M. & Du, G. Lift and power requirements of hovering insect flight. ActaMech.Sin.19, 458–469 (2003). 15. Sun, M. & Tang, J. Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion. J. Exp. Biol.205, 55–70 (2002). 16. Sane, S. P. The Aerodynamics of Insect Flight. J. Exp. Biol.206, 4191–4208 (2003). 17. Shyy, W., Aono, H., Chimakurthi, S. K., Trizila, P., Kang, C.-K., Cesnik, C. E., et al. Recent progress in flapping wing aerodynamics and aeroelasticity. Prog. Aerosp. Sci. 46, 284–327 (2010). 18. Sears, W. R. Some aspects of non-stationary airfoil theory and its practical application.Int. J. Aeronaut. Sci.8, 104–108 (1941). 19. Triantafyllou, M. S., Triantafyllou, G. & Yue, D. Hydrodynamics of fishlike swimming. Annu. Rev. Fluid Mech.32, 33–53 (2000). 20. Taylor, G. K., Nudds, R. L. & Thomas, A. L. Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency. Nature 425, 707–711 (2003). 21. Singer, E. K., Chang, L., Calderon, A. & Perez-Arancibia, N. O. Clip-brazing for the design and fabrication of micro-newton-resolution, millimeter-scale force sensors. Smart Mater. Struct. (2018). 22. Chen, Y., Ma, K. & Wood, R. J. Influence of wing morphological and inertial parameters on flapping flight performance in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (2016), 2329–2336. 23. Chen, Y., Gravish, N., Desbiens, A. L., Malka, R. & Wood, R. J. Experimental and computa- tional studies of the aerodynamic performance of a flapping and passively rotating insect wing. J. Fluid Mech.791, 1–33 (2016). 24. Ma, K. Y. Mechanical design and manufacturing of an insect-scale flapping-wing robot PhD thesis (The School of Engineering and Applied Sciences, Harvard University, 2015). 25. Finio, B. M., Pérez-Arancibia, N. O. & Wood, R. J. System Identification and Linear Time- Invariant Modeling of an Insect-sized Flapping-Wing Micro Air Vehicle in Proceedings of the 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2011) (San Francisco CA, Sep. 2011), 1107–1114. doi:10.1109/IROS.2011.6094421. 147 26. Zhao, L., Huang, Q., Deng, X. & Sane, S. P. Aerodynamic effects of flexibility in flapping wings. J. R. Soc. Interface 7, 485–497 (2010). 27. Shahzad, A., Tian, F.-B., Young, J. & Lai, J. C. Effects of flexibility on the hovering perfor- mance of flapping wings with different shapes and aspect ratios. J.FluidsStruct.81, 69–96 (2018). 28. Kang, C.-k. & Shyy, W. Scaling law and enhancement of lift generation of an insect-size hovering flexible wing. J. R. Soc. Interface 10, 20130361 (2013). 29. Singiresu, S. R. Mechanical Vibrations (Addison Wesley Boston, MA, 1995). 30. Lua, K.-B., Lai, K., Lim, T. & Yeo, K. On the aerodynamic characteristics of hovering rigid and flexible hawkmoth-like wings. Exp. Fluids 49, 1263–1291 (2010). 31. Nguyen, T., Sundar, D. S., Yeo, K. S. & Lim, T. T. Modeling and analysis of insect-like flexible wings at low Reynolds number. J. Fluids Struct.62, 294–317 (2016). 32. Shoele, K. & Zhu, Q. Performance of a wing with nonuniform flexibility in hovering flight. Phys. Fluids 25, 041901 (2013). 33. Du, G. & Sun, M. Effects of unsteady deformation of flapping wing on its aerodynamic forces. Appl. Math. Mech.29, 731–743 (2008). 34. Young, J., Walker, S. M., Bomphrey, R. J., Taylor, G. K. & Thomas, A. L. Details of insect wing design and deformation enhance aerodynamic function and flight efficiency. Science 325, 1549–1552 (2009). 35. Nakata, T. & Liu, H. A fluid–structure interaction model of insect flight with flexible wings. J. Comput. Phys.231, 1822–1847 (2012). 36. Nakata, T., Liu, H., Tanaka, Y., Nishihashi, N., Wang, X. & Sato, A. Aerodynamics of a bio- inspired flexible flapping-wing micro air vehicle. Bioinspir. Biomim.6, 045002 (2011). 37. Addo-Akoto, R., Han, J.-S. & Han, J.-H. Aerodynamic performance of flexible flapping wings deformed by slack angle. Bioinspir. Biomim.15, 066005 (2020). 38. Fu, J., Liu, X., Shyy, W. & Qiu, H. Effects of flexibility and aspect ratio on the aerodynamic performance of flapping wings. Bioinspir. Biomim. 13, 036001. doi:10.1088/1748-3190/ aaaac1 (2018). 39. Bhat, S. S., Zhao, J., Sheridan, J., Hourigan, K. & Thompson, M. C. Effects of flapping-motion profiles on insect-wing aerodynamics. J. Fluid Mech.884 (2020). 40. Jafferis, N. T., Smith, M. J. & Wood, R. J. Design and manufacturing rules for maximizing the performance of polycrystalline piezoelectric bending actuators.SmartMater.Struct.24, 065023 (2015). 148 41. Jafferis, N. T., Lok, M., Winey, N., Wei, G.-Y. & Wood, R. J. Multilayer laminated piezoelectric bending actuators: design and manufacturing for optimum power density and efficiency. Smart Mater. Struct.25, 055033 (2016). 42. Karpelson, M., Wei, G.-Y. & Wood, R. J.Areviewofactuationandpowerelectronicsoptionsfor flapping-wing robotic insects in IEEE International Conference on Robotics and Automation, 2008. ICRA 2008. (2008), 779–786. 43. Wood, R. J.Design,fabrication,andanalysisofa3DOF,3cmflapping-wingMAV inIEEE/RSJ International Conference on Intelligent Robots and Systems, 2007. IROS 2007. (2007), 1576– 1581. 44. Finio, B. M. & Wood, R. J. Distributed power and control actuation in the thoracic mechanics of a robotic insect. Bioinspir. Biomim.5, 045006 (2010). 45. Wood, R. J., Avadhanula, S., Sahai, R., Steltz, E. & Fearing, R. S. Microrobot design using fiber reinforced composites. J. Mech. Des.130, 052304 (2008). 46. Reddy, J., Barbero, E. & Teply, J. A plate bending element based on a generalized laminate plate theory. Int. J. Numer. Methods Eng.28, 2275–2292 (1989). 47. Deshmukh, B. B. & Jaju, S. B. Design and Analysis of Fiber Reinforce Polymer (FRP) Leaf Spring - A Review. Int. J. Eng. Sci. Technol.2, 289–291 (2011). 48. Liu, L., Zhang, B.-M., Wang, D.-F. & Wu, Z.-J. Effects of cure cycles on void content and mechanical properties of composite laminates,Compos.Struct.73, 303–309. doi:10.1063/1. 1851510 (2005). 49. Schmelzer, J. W. P., Zanotto, E. D. & Fokin, V. M. Pressure dependence of viscosity.J.Chem. Phys.122, 074511. doi:10.1063/1.1851510 (2005). 50. Wood, R. J., Avadhanula, S., Menon, M. & Fearing, R. S.Microroboticsusingcompositemate- rials: The micromechanical flying insect thorax in IEEE International Conference on Robotics and Automation, 2003. (ICRA’03).2 (2003), 1842–1849. 51. Wood, R., Steltz, E. & Fearing, R. Optimal energy density piezoelectric bending actuators. Sens. Actuator A Phys.119, 476–488 (2005). 52. Gafford, J. B., Kesner, S. B., Degirmenci, A., Wood, R. J., Howe, R. D. & Walsh, C. J. Amono- lithicapproachtofabricatinglow-cost,millimeter-scalemulti-axisforcesensorsforminimally- invasive surgery in 2014 IEEE International Conference on Robotics and Automation (ICRA) (2014), 1419–1425. 53. Sreetharan, P. S., Whitney, J. P., Strauss, M. D. & Wood, R. J. Monolithic Fabrication of Millimeter-Scale Machines. J. Micromech. Microeng.22, 055027 (2012). 149 54. Liu, X., Sun, Y., Wang, W. & Lansdorp, B. M. Vision-Based Cellular Force Measurement Using an Elastic Microfabricated Device. J. Micromech. Microeng.17, 1281 (2007). 55. Epstein, A. H. Millimeter-scale, MEMS gas turbine engines in ASME Turbo Expo 2003, collo- cated with the 2003 International Joint Power Generation Conference (2003), 669–696. 56. Popa, D. O. & Stephanou, H. E. Micro and Mesoscale Robotic Assembly. J. Manuf. Process. 6, 52–71 (2004). 57. Vaezi, M., Seitz, H. & Yang, S. A review on 3D micro-additive manufacturing technologies. Int. J. Adv. Manuf. Tech.67, 1721–1754 (2013). 58. Melle, S. M., Liu, K.,etal. Practical fiber-optic Bragg grating strain gauge system. Appl.Opt. 32, 3601–3609 (1993). 59. Peiner, E., Tibrewala, A., Bandorf, R., Biehl, S., Lüthje, H. & Doering, L. Micro force sen- sor with piezoresistive amorphous carbon strain gauge. Sens. Actuator A Phys.130, 75–82 (2006). 60. Behrens, I., Doering, L. & Peiner, E. Piezoresistive cantilever as portable micro force cali- bration standard. J. Micromech. Microeng.13, S171 (2003). 61. Dao, D. V., Toriyama, T., Wells, J. & Sugiyama, S. Silicon piezoresistive six-degree of freedom micro force-moment sensor. Sens. Mater 15, 113–135 (2002). 62. Kim, D.-H., Lee, M. G., Kim, B. & Sun, Y. A superelastic alloy microgripper with embedded electromagnetic actuators and piezoelectric force sensors: a numerical and experimental study. Smart Mater. Struct.14, 1265 (2005). 63. Wei, Y. & Xu, Q. An overview of micro-force sensing techniques. Sens. Act. A Phys. 234, 359–374 (2015). 64. Ito, Y., Kim, Y. & Obinata, G. Multi-axis force measurement based on vision-based fluid-type hemispherical tactile sensor in 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (2013), 4729–4734. 65. Wood, R. J., Cho, K. & Hoffman, K. A novel multi-axis force sensor for microrobotics ap- plications. Smart Mater. Struct.18, 125002 (2009). 66. Beyeler, F., Muntwyler, S. & Nelson, B. J. Design and calibration of a microfabricated 6-axis force-torque sensor for microrobotic applications in 2019 IEEE International Conference on Robotics and Automation, 2009. ICRA’09. (2009), 520–525. 67. Shen, Y., Xi, N., Li, W. J. & Tan, J. A high sensitivity force sensor for microassembly: design andexperiments inProceedingsofthe2003IEEE/ASMEInternationalConferenceonAdvanced Intelligent Mechatronics (AIM) 2 (2003), 703–708. doi:10.1109/AIM.2003.1225429. 150 68. Shen, Y., Xi, N., Wejinya, U. C. & Li, W. J. High sensitivity 2-D force sensor for assembly of surface MEMS devices in Proceedings of the 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) 4 (2004), 3363–3368. doi:10.1109/IROS.2004.1389936. 69. Schwartz, M. M. in. Chap. 2 (ASM International, Materials Park OH, 2003). 70. Olson, D. L.ASMHandbook:Welding,Brazing,andSoldering (ASM International, Materials Park OH, 1993). 71. De Gennes, P.-G. Wetting: statics and dynamics. Rev. Mod. Phys.57, 827 (1985). 72. Bonn, D., Eggers, J., Indekeu, J., Meunier, J. & Rolley, E. Wetting and spreading. Rev. Mod. Phys.81, 739 (2009). 73. Yang, Z., Zhang, L., Chen, Y., Qi, J., He, P. & Feng, J. Interlayer design to control interfacial microstructure and improve mechanical properties of active brazed Invar/SiO2–BN joint. Mater. Sci. Eng. A575, 199–205 (2013). 74. Villanueva, W., Boettinger, W. J., McFadden, G. B. & Warren, J. A. A diffuse-interface model of reactive wetting with intermetallic formation. Acta Materialia60, 3799–3814 (2012). 75. Yu, X., Yang, J., Yan, M., Hu, X.-W. & Li, Y.-L. Kinetics of wetting and spreading of AgCu filler metal over Ti–6Al–4V substrates. J. Mater. Sci.51, 10960–10969 (2016). 76. Lugscheider, E. & Partz, K. High temperature brazing of stainless steel with nickel-base filler metals BNi-2, BNi-5, and BNi-7. Weld. J.62, 160 (1983). 77. Singh, M., Shpargel, T., Morscher, G. & Asthana, R. Active metal brazing and characteri- zation of brazed joints in titanium to carbon–carbon composites. Mater. Sci. Eng. A 412, 123–128 (2005). 78. Long, W., Zhang, G. & Zhang, Q. In situ synthesis of high strength Ag brazing filler metals during induction brazing process. Scripta Materialia110, 41–43 (2016). 79. Wang, T., Liu, C., Leinenbach, C. & Zhang, J. Microstructure and strengthening mechanism of Si3N4/Invar joint brazed with TiNp-doped filler. Mater. Sci. Eng. A650, 469–477 (2016). 80. Zhou, Y., Hu, A., Khan, M., Wu, W., Tam, B. & Yavuz, M. Recent Progress in Micro and Nano-Joining. J. Phys.: Conf. Ser. 165, 012012. doi:https://doi.org/10.1088/1742-6596/ 165/1/012012 (2009). 81. Yost, F., Rye, R. & Mann Jr, J. Solder wetting kinetics in narrow V-grooves. Acta Materialia 45, 5337–5345 (1997). 82. Liu, C., Ou, C. & Shiue, R. The microstructural observation and wettability study of brazing Ti-6Al-4V and 304 stainless steel using three braze alloys.J.Mater.Sci.37, 2225–2235 (2002). 151 83. Suganuma, K., Miyamoto, Y. & Koizumi, M. Joining of ceramics and metals. Annu. Rev. Mater. Sci.18, 47–73 (1988). 84. Sharma, S., Sheoran, G. & Shakher, C. Investigation of temperature and temperature profile in axi-symmetric flame of butane torch burner using digital holographic interferometry. Opt. Lasers Eng.50, 1436–1444 (2012). 85. Chang, L. & Pérez-Arancibia, N. O. The Dynamics of Passive Wing-Pitching in Hovering Flight of Flapping Micro Air Vehicles Using Three-Dimensional Aerodynamic Simulations in Proceedings of the AIAA Atmospheric Flight Mechanics Conference, AIAA SciTech 2016 (San Diego, CA, Jan. 2016), 0013. doi:10.2514/6.2016-0013. 86. Fleming, A. J. A review of nanometer resolution position sensors: Operation and perfor- mance. Sens. Actuator A Phys.190, 106–126 (2013). 87. Liu, X., Peng, K., Chen, Z., Pu, H. & Yu, Z. A new capacitive displacement sensor with nanometer accuracy and long range. IEEE Sens. J.16, 2306–2316 (2016). 88. Gao, W., Kim, S., Bosse, H., Haitjema, H., Chen, Y., Lu, X.,etal. Measurement technologies for precision positioning. CIRP Annals 64, 773–796 (2015). 89. Muntwyler, S., Beyeler, F. & Nelson, B. Three-axis micro-force sensor with sub-micro- Newton measurement uncertainty and tunable force range. J. Micromech. Microeng. 20, 025011. doi:https://doi.org/10.1088/0960-1317/20/2/025011 (2009). 90. Viry, L., Levi, A., Totaro, M., Mondini, A., Mattoli, V., Mazzolai, B.,etal. Flexible three-axial force sensor for soft and highly sensitive artificial touch. Adv.Mater.26, 2659–2664 (2014). 91. Shiue, R., Chang, Y. & Wu, S. Infrared brazing Ti 50 Ni 50 and invar using Ag-based filler foils. Metall. Mater. Trans. A44, 4454–4460 (2013). 92. Yoon, J.-W., Noh, B.-I., Kim, B.-K., Shur, C.-C. & Jung, S.-B. Wettability and interfacial re- actions of Sn–Ag–Cu/Cu and Sn–Ag–Ni/Cu solder joints. J. Alloys Compd. 486, 142–147 (2009). 93. Corbacho, J. L., Suárez, J. C. & Molleda, F. Grain Coarsening and Boundary Migration dur- ing Welding of Invar Fe-36Ni Alloy. Mater. Charact. 41, 27–34. doi:https://doi.org/10. 1016/S1044-5803(98)00020-5 (1998). 94. Whitney, J. P. & Wood, R. J. Aeromechanics of passive rotation in flapping flight. J. Fluid Mech.660, 197–220 (2010). 95. Chang, L. & Pérez-Arancibia, N. O. Time-Averaged Dynamic Modeling of a Flapping-Wing MicroAirVehiclewithPassiveRotationMechanisms inProceedingsoftheAIAAAtmospheric Flight Mechanics Conference, AIAA Aviation Forum 2018 (Atlanta, GA, 2018), 2830. doi:10. 2514/6.2018-2830. 152 96. Ryu, Y., Chang, J. W. & Chung, J. Aerodynamic force and vortex structures of flapping flexible hawkmoth-like wings. Aerosp. Sci. Technol.56, 183–196 (2016). 97. Michelin, S. & Llewellyn Smith, S. G. Resonance and propulsion performance of a heaving flexible wing. Phys. Fluids 21, 071902 (2009). 98. Bhat, S. S., Zhao, J., Sheridan, J., Hourigan, K. & Thompson, M. C. Uncoupling the effects of aspect ratio, Reynolds number and Rossby number on a rotating insect-wing planform. J. Fluid Mech.859, 921–948 (2019). 99. Jafferis, N. T., Graule, M. A. & Wood, R. J. Non-linearresonancemodelingandsystemdesign improvementsforunderactuatedflapping-wingvehicles in2016IEEEInternationalConference on Robotics and Automation (ICRA) (2016), 3234–3241. 100. Tian, F.-B., Luo, H., Song, J. & Lu, X.-Y. Force production and asymmetric deformation of a flexible flapping wing in forward flight. J. Fluids Struct.36, 149–161 (2013). 101. Addo-Akoto, R., Han, J.-S. & Han, J.-H. Influence of aspect ratio on wing–wake interaction for flapping wing in hover. Exp. Fluids 60, 1–18 (2019). 102. Kang, C.-K., Aono, H., Cesnik, C. E. & Shyy, W. Effects of flexibility on the aerodynamic performance of flapping wings. J. Fluid Mech.689, 32–74 (2011). 103. Fuentes-Aznar, A., Gonzalez-Perez, I. & Pasapula, H. K. Computerized Design of Straight Bevel Gears with Optimized Profiles for Forging, Molding, or 3D Printing in AGMA 2016 Fall Technical Meeting (2016). 104. Dickinson, M. H. & Gotz, K. G. Unsteady aerodynamic performance of model wings at low Reynolds numbers. J. Exp. Biol.174, 45–64 (1993). 105. Yeo, D., Atkins, E., Bernal, L. & Shyy, W. Experimental investigation of the pressure, force, and torque characteristics of a rigid flapping wing in 50th AIAA Aerospace Sciences Meeting (2012), 849. 106. Scharnowski, S. & Kähler, C. J. Particle Image Velocimetry - Classical operating rules from today’s perspective. Opt. Lasers Eng. 135, 106185. doi:https://doi.org/10.1016/j. optlaseng.2020.106185 (2020). 107. Fincham, A. & Spedding, G. Low cost, high resolution DPIV for measurement of turbulent fluid flow. Exp. Fluids 23, 449–462 (1997). 108. Elsinga, G. E., Scarano, F., Wieneke, B. & van Oudheusden, B. W. Tomographic particle image velocimetry. Exp. Fluids 41, 933–947 (2006). 109. Michaelis, D. & Wieneke, B. Comparison between tomographic PIV and stereo PIV in Pro- ceedings of the 14th International Symposium on Applications of Laser Techniques to Fluid Mechanics. (2008). 153 110. Wieneke, B. Volume self-calibration for 3D particle image velocimetry.Exp.Fluids45, 549– 556 (2008). 111. Novara, M., Batenburg, K. J. & Scarano, F. Motion tracking-enhanced MART for tomo- graphic PIV. Meas. Sci. Technol.21, 035401 (2010). 112. Sciacchitano, A. Uncertainty quantification in particle image velocimetry. Meas. Sci. Tech- nol.30, 092001 (2019). 113. Adrian, R. J. & Westerweel, J. Particle Image Velocimetry 30 (Cambridge University Press, 2011). 114. Garcia, D. A fast all-in-one method for automated post-processing of PIV data. Exp. Fluids 50, 1247–1259 (2011). 115. Craven, P. & Wahba, G. Smoothing Noisy Data with Spline Functions. Estimating the Cor- rect Degree of Smoothing by the method of Generalized Cross-Validaton.NumerischeMath- ematik 31, 79 (1978). 116. De Boor, C. & De Boor, C. A practical guide to splines (Springer-Verlag New York, 1978). 117. Spedding, G., Rosén, M. & Hedenstrom, A. A family of vortex wakes generated by a thrush nightingale in free flight in a wind tunnel over its entire natural range of flight speeds. J. Exp. Biol.206, 2313–2344 (2003). 118. Ellington, C. The aerodynamics of hovering insect flight. III. Kinematics. Philosophical Transactions of the Royal Society of London. B, Biological Sciences 305, 41–78 (1984). 119. Percin, M., Hu, Y., van Oudheusden, B. W., Remes, B. & Scarano, F. Wing flexibility effects in clap-and-fling. Int. J. Micro Air Veh.3, 217–227 (2011). 120. Heathcote, S., Martin, D. & Gursul, I. Flexible flapping airfoil propulsion at zero freestream velocity. AIAA Journal 42, 2196–2204 (2004). 121. Heathcote, S. & Gursul, I. Flexible flapping airfoil propulsion at low Reynolds numbers. AIAA Journal 45, 1066–1079 (2007). 122. Wu, P., Stanford, B., Sällström, E., Ukeiley, L. & Ifju, P. Structural dynamics and aerody- namics measurements of biologically inspired flexible flapping wings. Bioinspir.Biomim.6, 016009 (2011). 154
Abstract (if available)
Abstract
For small biologically-inspired flying machines such as sub-gram flapping wing micro-air-vehicles (FW-MAVs) which operate in flows with Reynolds number Re ≤ 1000, unsteady aerodynamic mechanisms are critical in attaining net positive lift, control authority and payload capacity. The aerodynamic performance is sensitive to the wing stiffness, with complex three-dimensional and time-varying flow structures evolving along and behind the wing. This work provides a method for the design and fabrication of a microNewtonresolution dynamic force sensor to quantify flight performance metrics, but the micro-robotic vehicles are known for poor repeatability of the complex fabrication processes and highly coupled mechanical subsystems. It is difficult to study the fluid flow around the wings at-scale due to the physical size and actuation frequency of the system. Therefore this work uses a dynamically-scaled model undergoing prescribed sinusoidal flapping and pitching motion with zero mean flow, to study the induced flow and forces during a hovering flight condition. Tomographic PIV is used to visualize the three-dimensional vortex formation and force-generating flow structures while a 6-axis force-torque sensor at the wing root captures instantaneous data throughout the flapping period. The coefficients of lift and drag are computed for several cycles of flapping motion for four different wings that vary in chordwise stiffness and the traces are related to the flow structures observed in the induced velocity fields. The results show that chordwise flexible flapping wings in this Reynolds number regime demonstrate higher mean lift and lower RMS drag than their rigid counterparts; therefore chordwise FW-MAV wings have the potential to increase payload capacity and power economy without adding any additional weight to the vehicle.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Dynamic modeling and simulation of flapping-wing micro air vehicles
PDF
Aerodynamics at low Re: separation, reattachment, and control
PDF
Biomimetics and bio-inspiration for moderate Reynolds number airfoils and aircraft
PDF
Passive rolling and flapping dynamics of a heaving Λ flyer
PDF
Boundary layer and separation control on wings at low Reynolds numbers
PDF
Near wake characteristics of towed bodies in a stably stratified fluid
PDF
Synthesis and analysis of high-performance controllers for high-speed autonomous aerobatic flight
PDF
Traveling sea stars: hydrodynamic interactions and radially-symmetric motion strategies for biomimetic robot design
PDF
Pattern generation in stratified wakes
PDF
RANS simulations for flow-control study synthetic jet cavity on NACA0012 and NACA65(1)412 airfoils.
PDF
Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows
PDF
Model based design of porous and patterned surfaces for passive turbulence control
PDF
Novel soft and micro transducers for biologically-inspired robots
PDF
Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles
Asset Metadata
Creator
Singer, Emma Katharine
(author)
Core Title
Wing selection and design for flapping flight at small scale
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Degree Conferral Date
2022-05
Publication Date
04/05/2022
Defense Date
01/27/2022
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
flapping wing,FW-MAV,low Reynolds number,micro air vehicle,micro force sensing,OAI-PMH Harvest,tomographic PIV,wing flexibility
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Spedding, Geoffrey (
committee chair
), Luhar, Mitul (
committee member
), Nutt, Steven (
committee member
)
Creator Email
e.katharine.singer@gmail.com,eksinger@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC110878929
Unique identifier
UC110878929
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Singer, Emma Katharine
Type
texts
Source
20220406-usctheses-batch-919
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
flapping wing
FW-MAV
low Reynolds number
micro air vehicle
micro force sensing
tomographic PIV
wing flexibility